**2. 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

$$\mathbf{x}(t) = F^{-1}\mathcal{W}(j\omega) \cdot F[i(t)] \tag{1}$$

where *<sup>F</sup>* and *<sup>F</sup>*−<sup>1</sup> are the Fourier and inverse Fourier transforms, respectively, *<sup>W</sup>*(*jω*) is the transfer function, and *ω* is the angular frequency.

The physical model is based on the following assumptions.

2 Computational and Numerical Simulationse

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

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

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

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

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.

from a grounding system buried in multilayered earth model.

of Maclaurin's series to quickly calculate the Green's function.

coefficient have analytical formulae so as to avoid numerical integration.

up to some hundreds of kHz.

**2. Frequency domain analysis**

