Frequency Dependence of Multilayer Soil Electrical Parameters: Effects on Ground Potential Rise

—The paper presents a comprehensive study on the effect of frequency dependence of multilayer soil electrical parameters (SEPs) on ground potential rise (GPR). The lightning current has a large effect on the GPR. Typical analysis has often been based on quasistatic approximation that disregards the exact model of the grounding systems in the transient analysis of electrical networks, input impedance and GPR. To compute GPR in this paper, the full-wave approach based on the MoM solution to Maxwell equations is employed. To model the transient behavior of grounding systems, this method uses minimal approximations to find out GPR. The frequency dependence of SEPs combined with a multilayer soil are two complex issues that are modeled together. In the analysis of the GPR, special attention is given to the influences emanating from the soil multilayer structure and frequency dependency of SEPs. The GPR is calculated for three different soil configuration namely: 1) single layer, 2) two layer and 3) multilayer. From the simulations, the peak value of the GPR reduces when the frequency dependence of SEPs are taken into account. It is shown that considering the frequency dependent of multilayer SEPs can affect on the input impedance of grounding systems and GPR calculations.

packages, FD behavior of the grounding systems and precise model of grounding systems are disregarded.Thereby, the grounding system is modeled either as a lumped resistor or as an equivalent circuit of the lumped inductor, resistor and capacitor [12] whose values are calculated based upon the quasi-static assumptions.Thus, these models fail in providing accurate results at high frequencies for the calculation of lightning generated GPR.
More recently, other methods based on finite-difference time domain (FDTD) [13], method of moments (MoM) solution to Maxwell's equations [14], [15], [31] finite-element method (FEM) [16] have been also applied for analyzing the input Abstract-The paper presents a comprehensive study on the effect of frequency dependence of multilayer soil electrical parameters (SEPs) on ground potential rise (GPR).The lightning current has a large effect on the GPR.Typical analysis has often been based on quasistatic approximation that disregards the exact model of the grounding systems in the transient analysis of electrical networks, input impedance and GPR.To compute GPR in this paper, the full-wave approach based on the MoM solution to Maxwell equations is employed.To model the transient behavior of grounding systems, this method uses minimal approximations to find out GPR.The frequency dependence of SEPs combined with a multilayer soil are two complex issues that are modeled together.In the analysis of the GPR, special attention is given to the influences emanating from the soil multilayer structure and frequency dependency of SEPs.The GPR is calculated for three different soil configuration namely: 1) single layer, 2) two layer and 3) multilayer.From the simulations, the peak value of the GPR reduces when the frequency dependence of SEPs are taken into account.It is shown that considering the frequency dependent of multilayer SEPs can affect on the input impedance of grounding systems and GPR calculations.
Index Terms-Frequency dependence, Method of moment, Multi-layered soil, Ground potential rise, Input impedance.

I. INTRODUCTION
Ground potential rise (GPR) is defined in IEEE Std.367-2012 as the product of a ground electrode impedance, referenced to remote earth, and the current that flows through that electrode impedance.In other words, GPR is expressed by the total ground input impedance at the power system fundamental frequency to a remote earthing point.The GPR occurs when a large amount of current enters the soil originated by a fault or by lightning.It should be noted that the potential rise of the grounding systems of power transformers or lightning protective devices under the impact of lightning stroke is one of the main overvoltages sources on MV and LV equipment.
An exact analysis of GPR requires enough knowledge on the source modeling (e.g.lightning) [2], input impedance modeling [3], [4], and nature of the soil (e.g., electrical parameters frequency dependence [5], layers [6], and nonlinearity [7]).Therefore, to use an accurate model of soil is essential to calculate GPR precisely.By having a better 978-1-6654-2346-5/21/$31.00 ©2021 European Union impedance of grounding systems and GPR calculation.The evaluation of GPR has been usually carried out assuming the soil has an uniform medium characterized by constant electrical parameters.In general, soil frequency dependence effect is usually ignored in the calculation of the GPR.The constant SEPs assumption leads to considerable errors: the magnitude of GPR is around 60% higher than the measured one in the analyzed case in [17].Such large errors in the GPR show that the frequency dependence of SEPs should not be neglected in lightning-protection applications and GPR calculation.Soil frequency dependency is also disregarded in the input impedance of grounding systems electrodes buried in multilayer soil.In fact, the soil resistivity is assumed to be constant and equal to the value measured at low frequencies.Owing to these facts, accurate modeling of grounding system as a FD in a multilayer medium might help reducing the error GPR estimation [18].
Within this context in this paper, the GPR of a vertical electrode buried in multilayer soil subjected ligthning current is reduced by using MoM approach when considering FD effect.
This paper is organized as follows.Section II discusses the full wave approach for modeling of grounding systems.In this section, the curve-fit expression of FD soil model with application to analysis of grounding system proposed by Longmire and Smith is discussed.Section III presents the GPR calculation procedure based on MoM approach.Section IV describes the problem modeling and simulation results and relevant discussions are presented in Section V. Finally, general conclusions are provided in Section VI.

II. FULL-WAVE MOM APPROACH TO CALCULATE OF
GROUNDING SYSTEM INPUT IMPEDANCE The input impedance of grounding system is dependent only on the electrode geometry and electrical properties of the soil.
With reference to Fig. 1, a two-dimensional structure of the multilayer soil is shown.The resistivity and permittivity of lth layer is defined as ρ l and ε rl respectively.A vertical electrode buried in multilayer soil with FD resitivity and permittivity is considered.The thickness of each layer is characterized by h l .The grounding electrodes are divided into several small segments.The space is divided into several medias: air and the multi layer soil.For the input impedance calculation, we have applied MoM solution to Maxwell's equations to the calculation of grounding system impedance matrix in the frequency range of dc to 10 MHz.In this approach, the thin wire approximation is used.
It should be noted that the method used in this paper adopts a frequency-domain approach.Therefore, at the excited point, the time-domain lightning return stroke is identified by its spectral contents.A proper fast Fourier transform algorithm must be used for calculating in frequency domain.Considering the steps showed in [6] and [19], the spectral domain Green's functions are solved for obtaining the magnetic vector potential and the scalar potential.Then, by using the inverse Fourier transform, the spatial domain Green's function is calculated [20].
where GA,qe and G A,qe show the spectral domain and spatial domain Green's functions, respectively.The magnetic vector potential and the scalar potential are defined by A and qe superscripts respectively.The ρ-component of the wave vector for the lth layer is u ρ .The J 0 is the zero-order Bessel's function of the first kind.The radial distance in the cylindrical coordinate system is ρ.
+∞ 0 where u z is the z-component of the wave vector for the lth layer, r is the observation vector, and J 1 is the first order Bessel's function of the first kind.The more description of spectral domain Green's functions deriving and the appraisal of Sommerfeld's integral can be found in [21].
The mixed potential integral equation is caculated by enforcing the boundary condition for the electric current density, I, on the surface S placed in the multilayer soil.The relation between scattered and incident fields is explained in (4) where E s and E i are scattered and incident fields by an external source, respectively.To find the current distribution along the segments, electric field integral equation (EFIE) obtained by making use of the MoM.Equation ( 5) is presented EFIE for the distributed currents along the segments. n.
where µ is the magnetic permeability, G (r, r ) is the electric field Green's function at r due to a current element at r , I i (r ) is the unknown induced current along the electrode, and n is the unit vector tangential to the wire path x.The scattered field can be defined as (6).
where A, ω, and φ are the magnetic vector potential, angular frequency, and the electric scalar potential, respectively.The equations ( 7) and ( 8) are determined the magnetic vector potential and the electric scalar potential respectively.
where r is the source vector.The G A and G qe , respectively, are the spatial domain and spectral domain Green's functions, the spatial domain Green's function for magnetic vector potential, and spatial domain Green's function for scalar potential.Equation ( 9) is expressed the relation between the electric current density I(r) and the electric charge density ρ s (r).
In the MoM approach, the distribution of current on the electrodes of grounding system is expanded in a finite series as ( 10) where I n is the vector of unknown coefficients to be determined.F n (r) is basis function (triangular curve), which is utilized to distribute the electric current density in the vertical electrode.By substituting from ( 6) to (10) in (5) and multiplying the resulting equation by F n (r), and then integrating over the space, I n can be calculated from the resultant matrix equations (11).
where [A], [X], and [B] define as the impedance matrix, the unknown vector, and the excitation vector, respectively.

A. Frequency Dependent of Multilayer Soil
Frequency dependence effects of SEPs have been considered during the last few decades to evaluate the input impedance of electrodes buried in uniform soil [5] and [22].One of the first studies on FD of SEPs have been proposed by Smith-Rose [23] who presented experimental measurements of a variety of soil samples with different moisture contents in various frequencies ranging from 1 kHz to 10 MHz.At low frequencies range, the soil resistivity can be obtained by the socalled four-potential measuring technique (Wenner method).At high frequency range, measuring the attenuation of the electromagnetic wave is used to find the SEPs.The assumption of constant resistivity and permittivity leads to significant errors: the magnitude of the simulated GPR is around 60% higher than the measured one in the analyzed case [17].Such large errors in the GPR reveal that the dispersive properties of soil should not be ignored in lightning-protection applications.
The resistivity and related pemitivity of air are ρ = ∞ and ε r = 1 respectively.The resitivity and permittivity of soil are function of frequency.So, the electrical parameters of each soil layer should be defined a s f unctions o f frequency.
To this aim, analytical formulae that have been proposed by Longmire and Smith [24] based on experimental data [25] are used.They allow an analytical representation a curve-fit model for frequencies ranging from dc to 1MHz of FD of SEPs.
where the high frequency limit of dielectric constant is ε ∞ and set to 5. ρ LF is the low-frequency soil resistivity at 100 Hz; f is the frequency; The electrical parameters of each soil layer as functions of frequency ε r (f ) and ρ(f ) are defined, respectively, the soil resistivity and permittivity at each frequency; M is moisture content of soil; and, the coefficients a n are given in Table I for different values of n.
The overall process of the proposed method is presented in Fig. 2.

III. GROUND POTENTIAL RISE
A ground potential rise (GPR) occurs when a large current such as lightning flows to soil through an earth grid impedance.The potential relative to a distant point on the earth is highest at the point where current enters the soil, and declines with distance from the source.An engineering analysis of the power system under fault conditions can be used to determine whether or not hazardous step and touch voltages will develop.The result of this analysis can show the need for protective measures and can guide the selection of appropriate precautions.GPR is a concern in the design of electrical substations because the high potential may be a hazard to people or equipment.Grounding systems should guarantee the safety of equipment and human beings.Furthermore, recent experimental research demonstrates that the maximum value of GPR in the modeling at the injection point is considerably larger than the measured GPR value considering constant electrical parameters of soil in response to the impression of lightning currents to electrodes [26].Thus, SEP frequency dependence and layers also have a key role to reduce conservative estimation in GPR calculations.

A. Input Impedance of Grounding Systems
The grounding system consists of a vertical and horizontal electrodes, and it is a significant part of electrical system such as electrified railway systems, towers, substation, wind turbines, and large buildings.Often it is required to calculate the effect of the grounding system on the spread of surge voltages and currents within the electrical system during lightning strikes.For this aim, many researchers have developed different models for analyzing the transient behavior of grounding systems under lighting strikes.The precise modeling of grounding systems can help estimating the GPR with more accuracy.In this paper, the input impedance of grounding system is calculated by MoM in frequency domain.In MoM modeling, we define the SEPs as functions of frequency.In the simulations, a current of 1A is injected into the grounding system at each frequency.The excitation current spreads to an impressed half subsectional basis function located at top position feed point and the input impedance Z(ω) of the grounding system is computed as a function of frequency.
where V (ω), I(ω), and i(t) denote the electric potential phasor at the feed point in reference to the remote earth, the injected current phasor, and lightning current in time domain , respectively.The L denotes Laplace transform.The electric potential is calculated by integrating the electric field along a straight path starting from the electrode conductor surface to a point very far from the injection point in which the current approaches zero and could be physically considered as a voltage reference point.

B. Lightning Model
In the analysis, for the excitation current, a lightning current waveform is used to strike the tower.A Heidler's functions [27] corresponding to the typical first and subsequent return strokes are presented in Table II.Heidler's function parameters are used because of its convex shape similar to the real lightning and its zero time derivative at t = 0 is appropriated to represent the current [2].
η = e −(τ1/τ2)(n(τ2/τ1)) 1/n In ( 16), I 0 is the amplitude of lightning current; τ 1 is the front time constant; τ 2 is the decay time constant; n is exponent and set to 2; and η is the amplitude correction factor.
Note that the first stroke (FS) curve is reproduced by Heidler's functions with parameters presented in Table II.The sum of two Heidler's functions is used to represent the subsequent stroke (SS) current.The first stroke current has larger intensity, while the subsequent stroke current has a larger rate in its front.The constants in (16) for these two lightning stroke currents are shown in Table II.The subsequent stroke is greatly affected by the propagation behavior of a lightning current [28].The frequency content of FS lightning current is below 100 kHz, and SS current is below 1 MHz.Thus, SS suffers more from the propagation effect, which behaves as a higher transient ground impedance.In this paper, for further investigating of the impact of the frequency dependent, SS is considered for GPR calculation.Fig. 3 shows the curves affiliated with lightning first and subsequent stroke currents, obtained using Heidler's functions parameters [27] shown in Table II.

C. GPR Calculation
The input impedance of grounding system, Z(ω), is multiplied by I(ω) to obtain the frequency response of the grounding electrode (17).
where V (ω) is the GPR in frequency domain at the feed point.Then, using (18), the transient response is calculated by means of the inverse Laplace transform.v(t) is the GPR in time domain at injected point.
All simulated results comprising the input grounding impedance and GPR are presented in the section V, along with the analyses of some derived parameters, such as layer effect, FD effect, and the length of electrodes.
IV. PROBLEM MODELING Fig. 4 shows the different ground layer structures considered in this paper.For constant electrical parameters, the soil is assumed to have a soil resistivity ρ LF and a permittivity ε .These parameters are also utilized in ( 12) and (13) for the case in which the soil is considered to be FD.The evaluation is carried out for three structures of grounding system (single layer, two layer, and multilayer).After finding the input impedance of various geometry of grounding system, The GPR calculates based on input impedance and SS lightning impulse considering FD effects as explained in section II.The SEPs of each case and configuration are presented in Table III.
The cases 1 and 2 refer to a single layer soil, cases 3 and 4 represent two layer soil, whereas cases 5 and 6 represent multilayer soil for which the depth of the first and second layer are set to h=1m for two and multilayer structures.The burial depth is set to 0.6 m.It is assumed that the permittivity of grounding system considering constant parameters is equal to 10.The electrode cross section radius is 1.5 cm.The resistivity of each layer considering with and without FD of soil ia adopted from Table III at low frequency (100 Hz).The electrical parameters of all soil layers are defined as functions of frequency.We use analytical formulae that have been proposed by Longmire and Smith based on experimental measurements.Ionization of soil is disregarded in all simulations.

V. SIMULATION RESULTS
In order to evaluate the effect of frequency dependence of multilayer soil of electrical parameters on GPR, vertical grounding electrodes with different lengths buried in single layer and multilayer soil are studied considering both constant and FD parameters of soil.For the numerical simulation, the full-wave approach based on MoM is employed.In our simulations, Heidler's function [27] is adopted for modeling the lightning subsequent stroke currents.This study is done for three soil configuration: 1) single layer; 2)two layer; and 3) multilayer.show that the frequency-dependence of electrical parameters of soil can affect the high frequency behavior of ground electrodes.The frequency dependence of SEPs has a reduction effect on the amplitude of the grounding system input impedance and maximum GPR [5].As it can be seen from Figs. 5(a) and (b), for highly resistive soils, the input impedance of the vertical electrode at high frequencies show a capacitive behavior for both soils with constant and      FD parameters (i.e., the phase angle decreases to -90 and -50 degrees for constant and FD SEPs, respectively).Fig. 6 shows the simulation results of the GPR of the vertical electrode subjected to lightning SS current, under the constant and FD soil parameter assumptions.Some differences in the voltage peak and voltage waveform are observed, which could be caused by the frequency dependent effects of SEPs on the input impedance are underestimated which are dominant at the high frequencies.In these fast transient (e.g., subsequent lightning current), the frequency can reach up to 10 MHz that a slight change of grounding impedance would lead to large errors in GPR calculation.

B. Two layer soil
In simulation cases 3 and 4, the input impedance of vertical electrode for different electrode length of l= 3-m and 15-m is evaluated.These electrodes are buried in a two layer soil.In case 3, the upper and lower soil layers are characterized by electric resistivity of ρ 1 = 1000 Ωm, and the lower layer resistivity is ρ 2 = 10 Ωm, respectively.In the case 4, the upper and lower soil layers are characterized by electric resistivity of ρ 1 = 100 Ωm, and the lower layer resistivity of ρ 2 = 1000 Ωm,  respectively.The magnitude and the phase of these electrodes are shown in Fig. 7(a) and (b).
It is worth noting that the GPR shown in Fig. 8 are obtained while the vertical electrode subjected to subsequent lightning current under the constant and FD soil electrical parameter assumptions.From the results presented here, it is found that the vertical electrodes input impedance and GPR are affected by different parameters such as length of electrode, depth and SEP of various layers.By comparing the GPR waveform of FD and CP results in Fig. 8, we can find t hat G PR i s barely influenced b y t he i nput i mpedance, t he l ength o f electrode, and the soil resistivity at low frequency.As shown in Fig. 10, the peak difference of GPR between CP and FD models of soil electrical parameters for case 5 and 6 is 27 and 18 kV, respectively.The current propagation in the soil is related to soil layer characteristics, length of electrode, and parameters of the lightening current.The simulated results denote interesting aspects related to the impact of the frequency dependence of soil parameters on GPR.

D. Discussion
The previous simulations show that the assumption of constant soil electrical parameters gives significant errors when estimating the GPR [22].Additionally, all soil layers significantly contribute to the determination of the input impedance of grounding system and GPR.To further analyze the effect of multilayer soil and its frequency dependency, a reduction factor (RF) is used ( 19) where GPR(FD) and GPR(CP) are, respectively, the peak of the GPR considering FD soil parameters and the GPR considering constant electrical soil parameters [29].
The RF values presented in Tables IV and are computed for all cases subjected to subsequent lightning current.It is clear that the frequency dependence effect of SEPs on the GPR is mainly marked by the length of the electrode, the soil resistivity of layer, geometry of grounding system, and the lightning current properties.It could be also valuable to avoid overestimations of impulse impedance of grounding system.The errors are found for GPR under constant electrical parameters of soil assumption.The impulse impedance of grounding system is explained by ( 20) where V m the peak value of the GPR and I m is the peak value of the lightning current pulse at injected point [30].It is seen that the impulse impedance has correlation with the frequency dependence of SEPs.The impulse impedance is quite different from the low-frequency resistance, due to transient behavior of grounding system (capacitive or inductive), to the frequency dependence of soil parameters effects.

VI. CONCLUSION
In this paper, a precise full-wave MoM-based solution of Maxwell's equations for the evaluate of GPR considering the frequency dependence of multilayer SEPs has been used.The importance of the frequency dependence of multi layer soil parameters on the GPR subject to lightning currents has been investigated.To this aim, the GPR was computed for single layer, two layer, and multilayer soil structure with different grounding system geometries (vertical electrodes with different lengths) subjected to lightning subsequent current.It should be note that the peak value of the GPR reduces when the frequency dependence of SEPs are taken into account.It has been shown that considering the FD of multilayer SEPs could affect the impulse impedance of grounding systems.The assumption of constant soil resistivity and permittivity yields simulated GPRs much higher than the simulated ones considering frequency dependency effect of (a)

Fig. 8 .
Fig. 8. GPR of a vertical electrode subjected to lightning subsequent strokes current.(a) The vertical electrode of length l = 3 m, the upper layer resistivity is ρ 1 = 1000 Ωm, and the lower layer resistivity is ρ 2 = 10 Ωm; (b) the vertical electrode of length l = 15 m, the upper layer resistivity is ρ 1 = 100 Ωm, and the lower layer resistivity is ρ 2 = 1000 under the constant (purple-solid line) and FD soil parameter (green-dashed line) assumptions.The permittivity of both layers are ε 1 = ε 2 =10 at low frequency(100 Hz).

Fig. 9 (
Fig. 9(a) shows the simulation results to evaluate the effect of frequency dependent multilayer soil on the input impedance [absolute and phase values] of vertical electrodes.The results are obtained for a vertical electrode of length l = 3 m, the upper layer resistivity is ρ 1 = 10 Ωm, middle layer resistivity is ρ 2 = 100 Ωm and the lower

TABLE I COEFFICIENTS
an APPLIED IN LONGMIRE AND SMITH MODEL FOR FD OF SOIL MODELING

TABLE II HEIDLER
[27]UNCTIONS PARAMETERS USED TO REPRODUCE THE CURVE OF LIGHTNING CURRENT[27]

TABLE III LOW
FREQUENCY SOIL PROPERTIES OF GROUNDING SYSTEM CONFIGURATION SL:Single layer, TL:Two layer; ML:Multi layer; M:Percentage of soil moisture

TABLE IV RF
FOR THE LIGHTNING SUBSEQUENT STROKE CURRENT