Numerical Calculation for Lightning Response to Grounding Systems Buried in Horizontal Multilayered Earth Model Based on Quasi-Static Complex Image Method Numerical calculation for lightning response to grounding systems buried in horizontal multilayered earth model based on quasi-static complex image method

For lightning protection design, an exact evaluation of transients electromagnetic field within a complex grounding system has fundamental importance. In fact, the earth electrodes constitute a fundamental part of the electrical apparatus in both industrial and civil structures. Grounding systems should have a suitable configuration in order to avoid serious hazard to humans, and to preserve electrical insulation in electrical and electronic equipment and installations. Moreover, in electrical power installations, the shape and dimensions of the earth termination system, as a part of a lightning protection system, are more important than the specific value of the earth resistance, in order to disperse the lightning current into the earth without causing dangerous overvoltages.


Introduction
For lightning protection design, an exact evaluation of transients electromagnetic field within a complex grounding system has fundamental importance. In fact, the earth electrodes constitute a fundamental part of the electrical apparatus in both industrial and civil structures. Grounding systems should have a suitable configuration in order to avoid serious hazard to humans, and to preserve electrical insulation in electrical and electronic equipment and installations. Moreover, in electrical power installations, the shape and dimensions of the earth termination system, as a part of a lightning protection system, are more important than the specific value of the earth resistance, in order to disperse the lightning current into the earth without causing dangerous overvoltages.
Pioneering but comprehensive work on this subject was conducted in the first half of the twentieth century, which is summarized by Sunde in the well known reference book [1]. Important pioneering work is described also in [2] and [3]. More recent work is summarized in [4]. Recently, computerized analysis methods have been developed based on different approaches, for example, on circuit theory [5]- [8], transmission line theory [9]- [15], electromagnetic field theory [16]- [23], and hybrid method [24]- [31].
Hybrid method has been developed from conventional nodal analysis, which combines the electrical circuit method and the electromagnetic field method. It has been proved to have combined the strong points of both the two methods. Dawalibi earlier discussed how the hybrid method came out of the electromagnetic field method in [24] and further discussed the hybrid method in [25], however, the hybrid method was based on quasi-static electromagnetic field theory, and only discusses steady grounding problem in the frequency domain. Meliopoulos also discusses the hybrid method based on quasi-static electromagnetic field theory in [26]; Huang & Kasten developed a new hybrid model to calculate the current distribution in both the grounding system and the metallic support conductors, while considering the voltage drop along the grounding system conductors [27]. However, the model was also based on quasi-static electromagnetic field theory, meanwhile, leakage currents and network currents within the grounding system are separately considered in the calculation, their mutual coupling influence is neglected, and the capacitive coupling effect of the earth is also neglected. Otero, Cidras & Alamo developed a hybrid method to calculate the current distribution in a grounding system [28] within which the mutual inductive and capacitive coupling influence among these current flowing and leaking along the conductor is considered. However, only a uniform half infinite earth model is considered. This hybrid method was combined with the FFT, and so the transient response from the grounding system was obtained. These confines in the frequency domain promoted the development of a novel mathematical model in [29]- [32], which introduced the quasi-static complex image method (QSCIM) for calculating the current distribution in a grounding system buried in both horizontal and vertical multilayered earth models in the frequency domain. However, the hybrid method can be further developed to numerically calculate the transient response from a grounding system buried in multilayered earth model.
Once the multilayered earth model is adopted, the Green's function of a point source will contain an infinite integral for the Bessel function, a complex image method based on Maclaurin's infinite series expansion have been studied in [33] and [34], which has brought up problem about the convergence of the infinite Maclaurin's series. To avoid this convergence problem, QSCIM is introduced to dealt with the infinite integral, which uses finite exponential terms (usually just 3-4 terms) through the Matrix Pencil approach instead of Maclaurin's series to quickly calculate the Green's function.
In this paper, based on previous works [28]- [31], combined with the FFT, a novel and accurate mathematical model is developed for calculating the harmonic wave currents of lightning currents distribution along the grounding system buried in multilayered earth model in the frequency domain, within which not only the conducting effect of the harmonic wave currents leaking into the soil, but also capacitive and inductive effects between different layers of soil have been considered. Both leakage currents and network currents within the grounding system and their mutual coupling are considered in the calculation. The earth is modeled by a multilayered earth model. To accelerate the calculation, QSCIM and closed form of Green's function were introduced, and the mutual inductive and conductive coefficient have analytical formulae so as to avoid numerical integration.
The maximum frequency of applicability of the method is limited by the quasi-static approximation of the electromagnetic fields. For the usual electrodes, it may be applied up to some hundreds of kHz.

Frequency domain analysis
The transient problem is first solved by a formulation in the frequency domain. The time-domain response is then obtained by application a suitable Fourier inversion technique. The response to a steady state, time harmonic excitation is computed for a wide range of frequencies starting at zero Hz. From this frequency response, a transfer function is constructed for every frequency considered. The transfer function is dependent only on the geometric and electromagnetic properties of the grounding system and its environment.
if i(t) represents the injected current at a point in the grounding system, and x(t) denotes an observed response, then where F and F −1 are the Fourier and inverse Fourier transforms, respectively, W(jω) is the transfer function, and ω is the angular frequency.
The physical model is based on the following assumptions.
• The earth comprises horizontal multilayered media, and the air media are homogenous and occupy half-spaces with a common horizontal plane boundary between the air and earth. • The earth and the grounding electrodes exhibit linear and isotropic, arbitrary characteristics. • The grounding system is assumed to be made of cylindrical metallic conductors with arbitrary orientation. However, they are assumed to be subject to the thin-wire approximation, i.e., the ratio of the length l of the conductor segment to its radius r is l r ≫ 1. In practice, a ratio of about 10 is satisfactory. • Energization occurs by the injection of a current impulse of arbitrary shape produced by an ideal current generator with one terminal connected to the grounding system, and the other to the ground at infinity. The influence of the connecting leads is ignored.

Mathematical model of the equivalent circuit of the grounding system in the frequency domain
A set of interconnected cylindrical thin conductors placed in any position or orientation makes up a network to form the grounding system. The grounding network's conductors are assumed to be completely buried in a conductive N e -layer media (earth) with conductivity σ e and permittivity ε e = ε r e · ε 0 (here e = 1, · · · · ··, N e ). The air is assumed to be a non-conductive medium with permittivity ε 0 = 10 −9 36π F/m. All media have permeability µ = µ 0 = 10 −7 4π H/m. The proposed methodology is based on the study of all the inductive, capacitive and conductive couplings between the different grounding system conductors. First, the electrode is divided into N l pieces of segments that can be studied as elemental units, where the discrete grounding system has N p nodes. A higher segmented rate of the electrode can enhance the model's accuracy but increases its computational time. Therefore, it is necessary to achieve a compromise solution between the two determinants.
The grounding network is energized by injection of single frequency currents at one or more nodes. In general, we consider that a sinusoidal current source of value F j is connected at the jth (j = 1, 2, . . . , N p ) node. A scalar electric potential (SEP) V j of jth node on the grounding network referring to the infinite remote earth as zero SEP is defined. In the same way, we define an average SEP U k on kth (k = 1, 2, . . . , N l ) segment. If the segments are short enough, it is possible to consider U k as approximately equal to the average of the kth segment's two terminal nodes SEP. We define a branch current I k b , branch voltage U k b , and leakage current I k s on the kth (k = 1, 2, . . . , N l ) segment.

Mathematical model of the grounding system in the frequency domain
With the above considerations, according to [28]- [31], the electric circuit may be studied using the conventional nodal analysis method [35], resulting in the following equations: where [F] is an N p × 1 vector of external current sources; [Z b ] is the N l × N l branch mutual induction matrix of the circuit including resistive and inductive effects, which gives a matrix relationship between branch currents [I b ]; [Z s ] is an N l × N l mutual impedance matrix, which gives a matrix relationship between the average SEP [U] and leakage currents [I s ] through the rapid Galerkin moment method [38]. Both [A] and [K] are incidence matrices, which are used to relate branches and nodes. There are rectangular matrices of order N l × N p , for whose elements we refer to [28]- [31]. can also be calculated [28]- [31].
Once the branch currents and leakage currents are known, the SEP at any point can be calculated by where G ϕ (r j ,r i ) is the scalar Green's function of a monopole in the multilayered earth model.
The vector magnetic potential (VMP) A at any point can be calculated bȳ Here,Ḡ A (r j ,r i ) is the dyadic Green's function of a dipole in the multilayered earth model, which will be introduced later.
The electrical field intensity (EFI) at any point can be calculated bȳ Computational and Numerical Simulations 396 Numerical calculation for lightning response to grounding systems buried in horizontal multilayered earth model based on quasi- The magnetic field intensity (MFI) at any point can be calculated bȳ The study of the performance of the grounding system in the frequency domain has been reduced to the computation of [Z s ] and [Z b ] matrices.

Computation of [Z b ] and [Z s ] matrices
From [28]- [31], we know that each segment is modeled as a lumped resistance and self-inductance. Mutual inductances or impedances between branch segments' branch currents or leakage currents are also included in the model: .. X 1,N l : : : : : : : ... : ... : : : : ... : ... : : : The diagonal elements consists of self impedance and self induction, the other elements belong to mutual induction between a pairs of conductor segments. The formula for self impedance and self induction can be found in [28]- [31], here, we give the formula for mutual induction: For an infinite homogeneous conductivity medium, one hasḠ A (r j , The above integral can be analytically calculated, this formula can be found in [27].
Z i,j is the mutual impedance coefficient between a pair of conductor segments in the grounding system. The matrix above includes the conductive and capacitive effects of the earth, and its elements are the mutual impedance coefficients Z i,j .
Note the medium surrounding the point current source here was considered as homogeneous and infinite. However, in any practical case, the earth is represented via a multilayered earth model. The QSCIM can be used to dealt with the multilayered earth model, this will be discussed next.

The closed form of the Green's function of a point source in a horizontal multilayered earth model and the QSCIM
The closed form of the Green's function of a scalar monopole and vector dipole buried in horizontal multilayered earth model will be respectively introduced.

The closed form of the Green's function of a scalar point source in a horizontal multilayered earth model and the QSCIM
The main task of simulation grounding system is calculating the element Z jk of matrix [Z s ] in Eq. (11). When the earth model is considered as horizontal multilayered conductivity media, influences from the interface will be considered, which will lead to infinite integral about Bessel function associated with Green's function of a scalar monopole. However, the element Z jk of matrix [Z s ], which includes the Green's function, can be fast calculated by using the QSCIM.
To perform the calculation, the corresponding Green's function of a scalar monopole must be defined first. The Green's function can be regarded as the SEP of one point located at any place produced by a scalar monopole with unit current δ in a horizontal multilayered conducting medium. For low frequency (50 or 60 Hz and higher harmonic wave) and limited size of the substation, the propagating effect of electromagnetic wave can be neglected, so

Computational and Numerical Simulations 398
Numerical calculation for lightning response to grounding systems buried in horizontal multilayered earth model based on quasi- where i, j = 1, . . . , N s , and δ(r, r ′ ) is the Dirac delta function. δ(ij) is Kronecker's symbol, σ i = σ i + jωε i is the complex conductivity of the ith layer medium. Supposing the monopole is located at origin of the coordinate system, seen from Fig. 1.
If the monopole lies in infinite homogenous medium, its expression of Green's function in the spherical coordinate system is.
where R = |r −r ′ | is the distance between the source point and field point. While its expression in the cylindrical coordinates system (ρ, z) is as follows: (15) where J 0 (λρ) is the Bessel function of the first kind of order zero. For horizontal multilayered earth model, for examples, two-layer earth model, in the cylindrical coordinate system, the general form of Green's functions can be expressed by [1]: (18) where G 10 , G 11 and G 12 express the Green's function for the field point in the air, the top layer and bottom layer, respectively.
The six constants A 10 ∼ B 12 can be determined by employing the interface and infinite conditions of the SEP (ϕ 10 = ϕ 11 , σ 0 Consequently, the expression of G 11 can be given as follows: where , h is the thickness of the top layer, σ 0 ,σ 1 and σ 2 are the complex conductivity of air and two layers, respectively. Now this f (λ) can be developed into an exponential series with finite terms by MP method [37] instead of Maclaurin's series to avoid verbose calculation [38]: where α n and β n are constants to be determined by choosing sample points of function f (λ).
Considering the Laplace transform of Bessel function of the first kind of order v J v (λρ), we have [39] ∞ where Re(v) > −1 and Re(c) > 0.
Set v = 0, so we have Lipschitz integration: By employing Eq. (20) and Lipschitz integration, the closed form of the Green's function (Eq. (19)) can be given as following: where the origin of the coordinate system shown in Fig. 1 has been moved to the surface of the earth, the source point is at (0, z ′ ) and the field point is at (ρ, z), It can be seen that each term except 1 R 0 of Eq. (23) can be regarded as a image scalar monopole source, whose location is indicated by R ni,i=1∽4 and amplitude is α n , which can be seen in Fig. 1. However, in Eq. (23) α n and R ni are usually complex numbers, so that this approach is named as the QSCIM. G 10 and G 12 can be similarly gotten.
The procedure for the closed form of the Green's function of a scalar monopole in arbitrary horizontal multilayered earth model is first to derive the specific solution of the Green's function in each layer based on the general expressions (for example, the general expressions of a scalar monopole in three layers media is given in Eq. (15), Eq. (16) and Eq. (17)), and the interface and infinite conditions of SEP of the monopole is used to decide the unknown constants. Then, by employing the MP method exponential series development and the Lipschitz integration, the final expression of the Green's functions can be obtained.
Once the closed form of Green's function of a scalar monopole in the horizontal multilayered earth model has been gotten, the mutual impedance coefficient Eq. (11) can be fast analytical calculated, which is just same as the Eq. (12).

The closed form of the Green's function of a vector point source in a horizontal multilayered earth model and the QSCIM
The another important task of simulation grounding system is calculating the element M jk of matrix [Z b ] in Eq. (9). When the earth model is considered as horizontal multilayered conductivity media, and influences from the interfaces on mutual induction between two conductors should be considered, which also lead to infinite integral about Bessel function in the Green's function of a vector dipole, the element M jk of matrix [Z b ] includes the Green's function, which can also be fast calculated by using the QSCIM.
Like Green's function of the scalar monopole, to perform the calculation of a vector dipole in horizontal multilayered earth model, its corresponding Green's function should also be defined first. The Green's function can be regarded as the VMP A of a point at any place produced by a vector dipole with unit current δ within a horizontal multilayered medium, whose electromagnetic field can also be regarded as quasi-static field, so the VMP A also satisfies the Poisson equation as: where i, j = 1, . . . , N s ; µ is the permeability of the earth medium. Supposing the dipole is located at origin of the coordinate system, seen from Fig. 1 If the vector dipole is lying in homogenous infinite medium, its expression of Green's function in the spherical coordinate system is.
And its expression in the cylindrical coordinates system (ρ, z) is as follows: Not like scalar monopole source, for horizontal multilayered earth model, vector dipole source must be considered into two cases, which are vertical and horizontal dipoles, respectively. This is because any placed dipole in horizontal multilayered earth model can be decomposed into horizontal and vertical components. For vertical dipole case, its Green's function can be defined as G z A z ; For horizontal dipole case, it own two components, which are horizontal x or y component and vertical z component, so its Green's function can be defined as G x A x and G z A x , respectively. We also take the two-layer earth model as an example, the dipole lies in the earth model. In the cylindrical coordinate system, the general form of Green's functions for the vertical dipole G z A z or horizontal dipole G x A x can be expressed by [1]: where G p A p10 , G p A p11 and G p A p12 express the Green's function of the dipole for the field point in the air, the top layer and bottom layer, respectively; for vertical dipole case, p = z, for horizontal dipole case, p = x.
Just like the deduced procedure for Green's function of a scalar monopole, the six constants A 10 ∼ B 12 can also be decided by employing the interface and infinite conditions of the dipole's VMP. For vertical dipole, there is a f (λ), which can also be developed into an exponential series with finite terms by the MP method. Last by applying the Lipschitz integration, the final expression of G p A p11 can be given as follows.
where the origin of the coordinate system shown in Fig. 1 has been moved to the surface of the earth. R 0 , R ′ 0 , R ni,i=1∽4 are just same as Eq. (23). For horizontal dipole, there is a Green's function for z component G z A x11 left, in the cylindrical coordinate system, the general form of the Green's function for z component in the two-layer earth model can be expressed by [1]: Also like the deduced procedure for Green's function of a scalar monopole, the six constants A 10 ∼ B 12 can be decided by employing the interface and infinite conditions of the horizontal dipole's VMP. The f (λ) in the Green's function can also be developed into an exponential series with finite terms by the MP method. Last by applying the Lipschitz integration's varied form ( ∞ 0 e −cλ J 0 (λρ) dλ λ = ln(c + c 2 + ρ 2 )), the final expression of G z A x11 can be given as follows.
where the origin of the coordinate system shown in Fig. 1 has been moved to the surface of the earth. The R ′ 0 and R ni,i=1∽4 are same as Eq. (23). Other parameters are z ′ 0 = z + z ′ , z n(1−4) = (sign a · z + sign b · z ′ ) + z n , in which sign a = 1 for z n1 and z n3 , sign b = 1 for z n3 and z n4 , sign a = sign b = −1 for others.
It can also be seen that each term except 1 R 0 of Eqs. (30) and (35) can be regarded as a image vector dipole source, whose location is indicated by R ni,i=1∽4 and amplitude is α n , which can be seen in Fig. 1. However, in Eqs. (30) and (35), α n and R ni are usually complex numbers, so that this approach is also named as the QSCIM.
The other Green's function for the dipole G z A z10 , G z A z12 , G x A x10 , G x A x12 , G z A x10 and G z A x12 can be given in the similarly way. Once the closed form of Green's function of a vector dipole in horizontal multilayered earth model has been given, the mutual impedance coefficient Eq. (9) can also be fast analytical calculated just like Eq. (10).

Time-domain solutions
Once the transfer functions W(jω) have been determined for each calculated quantity, for example the electric field or current at specified points, the time-domain solutions can be obtained by direct application of (1). The calculation of the inverse Fourier transform is carried out by an FFT algorithm which is well-suited for the evaluation of time-domain responses.
The transient impedance, an essential parameter in grounding system design, is defined as a ratio of time varying voltage and current at injection point [19]: where i(t) represents the injected current at a grounding grid. This injected current represents the lightning channel current usually expressed by the double exponential function i(t) = I m (e −αt − e −βt ), t ≥ 0, where pulse rise time is determined by constants α and β, while I m denotes the amplitude of the current waveform. The Fourier transform of the excitation function is defined by integral [40]: Integral (37) can be evaluated analytically The frequency components up to few MHz are meaningfully present in the lightning current Fourier spectrum with strong decreasing importance from very low to highest frequencies.
Multiplying the excitation function I( f ) with the input impedance spectrum Z in ( f ) provides the frequency response of the grounding system: Applying the IFFT, a time domain voltage counterpart is obtained. IFFT of the function U( f ) is defined by the integral [40]: As the frequency response U( f ) is represented by a discrete set of values the integral (41) cannot be evaluated analytically and the discrete Fourier transform, in this case the IFFT algorithm, is used, i.e., Implementation of this algorithm inevitably causes an error due to discretization and truncation of essentially unlimited frequency spectrum. The discrete set of the time domain voltage values is defined as [40]: where n = 0, . . . , N, F is the highest frequency taken into account, N is the total number of frequency samples, ∆ f is sampling interval and ∆t is the time step.
Finally, the impulse impedance, an also essential parameter in grounding system design, is defined by the following expression [41]: where U is the voltage maximum at the discharge point and I is the injected current magnitude at the time instant when U has been reached.

Verification of the method
In this part, our model has been verified through comparison with experimental data from other published papers; the validation of our model will also be discussed.

Verification of the method
To verify the method proposed in this work, some cases solved by other authors are studied.
In the first case, from [42], which gives some grounding impedance measurement results, the grid used for measurement is a 16-mesh grid of 100 m × 100 m. The earth is modeled by a multilayered earth model. The earth model is given by (a) ρ 1 = 25Ω m, ρ 2 = 120Ω m, and thickness of the upper earth is 5 m; (b) ρ 1 = 120Ω m, ρ 2 = 25Ω m, ρ 3 = 250Ω m, thickness of the first and second layers of earth are 5 m and 9 m, respectively. The radius of the grid conductors is 0.5 cm, and they are buried at 0.5 m depth in the earth. Here, the conductivity of the copper conductors is σ Cu = 5.8 × 10 7 S/m, and the permittivity of earth is set to ε 1 = 5. The results can be seen in Table 1.
In the second case, which is from [41], a typical grounding grid configuration can be seen in Fig. 2, which was made of round copper conductors with 50 mm 2 cross section. The grounding grid was buried at 0.5 m depth in two-layer horizontal earth, whose resistivity ratio for the upper and the lower soil layers is ρ 1 /ρ 2 = 50/20, the upper layer thickness being H = 0.6 m. The inject lightning current parameter was set to T 1 = 3.5µ s, T 2 = 73µ s and I m = 12.1 A, the feed point is at the corner of the grid. For our model, the permittivities of the two layer earth model were set to ε 1 = 30ε 0 and ε 2 = 20ε 0 . The transient SEP can be seen in Fig. 3, which ultimately agreed with the measured curve in Fig. 4 (a) in [41]; meanwhile, the impulse grounding impedance was 2.12 Ω as given by [41], and it is 2.08Ω for our model.

Validation of our method
The maximum frequency of applicability of the method is limited by the quasi-stationary approximation of the electromagnetic fields, which means the propagation effect of the electromagnetic field around the grounding system can be neglected, so where γ 2 e = jωµ (jωε e + σ e ) , e = 1, . . . , N e . For most of the usual electrodes, this may be applied up to some hundreds of kHz.

Simulation result and analysis
A typical grounding system can be seen in Fig. 4. The earth is medalled by a two-layer conductive earth model, whose conductivity and permittivity is σ 1 = 500 −1 S/m, σ 2 = 900 −1 S/m, ε 1 = 10ε 0 , ε 2 = 22ε 0 , respectively, the first layer height is 5 m. The material of the grounding system conductor is Cu with conductivity σ Cu = 5.8 × 10 7 S/m. The conductor radii are 7 mm. The external excited lightning current is injected from the corner of the grounding system, which is described by a double-exponential function: I(t) = 1.29 × (e −0.019010t − e −0.292288t ) kA, which means that the parameters of the lightning current are T 1 = 10µ s, T 2 = 50µ s and I m = 1.29 kA, the lightning current has been shown in Fig. 5.
A comparison between chosen total leakage currents and injecting currents of the grounding system in frequency domain is given in Table 2. It can be seen that the total leakage current of the grounding system is close to the external injected current. All this shows the accuracy of this model.

Computational and Numerical Simulations 408
Numerical calculation for lightning response to grounding systems buried in horizontal multilayered earth model based on quasi- The quasi-static complex image in this case has two terms, the α n and β n can be seen below Table 3.  The distribution of absolute values of grounding impedance dependance |Z(jω)| on frequency can be seen in Fig. 6. This figure shows that |Z(jω)| is independent of the frequency below 100 kHz and equal to the low frequency grounding impedance, which agrees with the viewpoint of [43].
To further discuss the electromagnetic field characteristics along the surface above the grounding grid, the distribution of the electromagnetic field along the surface with three different chosen frequencies (17 kHz, 250 kHz and 800 kHz) have been given in Figs From Figs. 7-9, we know that the ground SEP rise is dependent on the magnitude of the injecting current, the ground SEP rise at 17 kHz is generated by injecting current with (-4.077,-10.964) kA, the ground SEP rise at 250 kHz is generated by injecting current with (-0.289,+ 0.143) kA, and the ground SEP rise at 800 kHz is generated by injecting current   with (-4.301,+ 3.428) kA. So the ground SEP rise at 17 kHz is the maximum, and the ground SEP rise at 250 kHz is the minimum. Meanwhile, we can see that almost an equipotential surface occurs for the low frequency case (17 kHz and 250 kHz), and at higher frequencies, the impedances of the grid conductors are no longer negligible and most of the earth currents dissipate close to the injection point. This phenomenon is well illustrated in Fig. 9 which Numerical calculation for lightning response to grounding systems buried in horizontal multilayered earth model based on quasi- represents the case of a corner current injection. Earth potentials near the injection corner present a very sharp peak, while they are quite low and flat everywhere else, with a minimum at the opposite corner.
From Figs. 10-12, we can see that as in the ground SEP rise case, the maximum value of the x-componential of the EFI E x is dependent on the magnitude of the injecting current. The distribution of the x-componential of the EFI E x is along the x-direction. Meanwhile, the distribution of the electrical field is not dependent on the current injection location at low frequencies. This conclusion no longer holds at higher frequencies, as shown in Fig. 12,  which clearly shows that the region of maximum electrical field shifts from the edge of the grid to the injection point location as the frequency increases from low to high.
The distribution of the x-componential of MFI B x is given in Figs. 13-15. We can also see that as in the ground SEP rise case, the maximum value of the x-componential of the MFI B x is dependent on the magnitude of the injecting current. Unlike the distribution of the x-componential of the EFI E x , the distribution of the x-componential of the MFI B x is along the y-direction, which can be easily explained in that the distribution of the electrical field is perpendicular to that of the magnetic field. Meanwhile, the distribution of the electrical Numerical calculation for lightning response to grounding systems buried in horizontal multilayered earth model based on quasi- field is dependent on the current injection location from low to higher frequencies, which is different from the electrical field case.

Conclusion
A novel mathematical model for accurately computing the lightning currents flowing in the grounding system of a high voltage a.c. substation, buried in multilayered earth, has been developed in this paper. Together with the FFT, not only the conducting effect of harmonic wave components of these currents, but also capacitive and inductive effects from the interface between different soil layers have been analyzed in the frequency domain. To accelerate the calculation, the QSCIM and a closed form of Green's function were introduced. With the inverse FFT, the model can calculate the distribution of lightning currents in any configuration of the grounding system. This can be used for studying the performance of transient lightning responses to grounding systems. Last, the model has been validated through some numerically simulated and experimental results from open published paper, and some numerical results have been discussed in this paper.

Computational and Numerical Simulations 414
Numerical calculation for lightning response to grounding systems buried in horizontal multilayered earth model based on quasi-