AbstractThis work proposes a novel mathematical model based on the Galerkin time‐domain boundary element method for accurately calculating the lightning current distribution and lightning impulse response of the buried substation grounding grid, in a multi‐layer horizontal layered soil model, by taking into account the soil ionization effect. To improve computational efficiency, the quasi‐static complex image method and its closed form time‐domain Green's function have been introduced into the model that has the ability to analytically calculate the mutual inductance coefficient between the branch currents of any two conductor segments and the mutual resistance coefficient between the leakage currents. The Galerkin time‐domain boundary element method proposed in this work can simulate the transient lightning impulse response of a substation grounding grid buried in the multi‐layer horizontal layered soil.