Transient Overvoltage Analysis in the Medium Voltage Substations Based on Full-Wave Modeling of Two-layer grounding System

In this paper, a comprehensive investigation of substation grounding soil structure effects on switching/lightning transients is presented. The full-wave computational solution Method of Moments (MoM) is adopted to solve Maxwell’s equations. The significant advantage of this technique is its ability to deal effectively with substation grounding systems (SGSs) at high frequency. This evaluation is carried out considering two soil structures: (1) a uniform soil and (2) soil with a layered structure. The impacts of both soil structures on the peak transient overvoltage values are compared. To begin with, the harmonic impedance matrix (HIM) of the layered multi-terminal SGS is computed. Then, a vector fitting approach is employed to incorporate the accurate SGS behavior by relying on the obtained HIM into time domain simulations. Finally, the transient overvoltages and maximum ground potential rise (GPR) are computed. The chief contribution of this paper is the transient assessment of the two-layer multi-terminal soil structure of the large SGS taking into account the frequency dependence model of the electrical parameters of soil, which is not generally investigated and instead, treated simplistically. According to the results of the study, the lack of an exact grounding system models may lead to underestimations/overestimations of the GPR and the overvoltage values. The results demonstrate the need for the grounding modeling of SGS considering layered soil structure to be incorporated into power system transient analysis in a wide range of frequencies.


Introduction
A well-designed substation grounding systems is a remarkably effective operational factor in the field of protection, human safety, insulation coordination, and electromagnetic compatibility.The grounding system of substations is composed of vertical and horizontal electrodes inserted in either uniform or layered soil to achieve a suitable grounding impedance [1][2][3].Understanding how soil layers, depths, and traverses affect ground impedance are important.The SGS can reduce the risk to acceptable levels posed by source disturbance (e.g., lightning, switching, and ground fault) in the electrical equipment [4].It has a dynamic (.., inductive or capacitive) behavior and can affect the electrical system components since all are electrically connected to the ground terminal through earth terminations [5,6].While the SGS is being excited by lightning with high-frequency content or fast front waveform [4], power transformers (PTs), supporting insulators between ground and active parts, and other substation equipment may experience dangerous overvoltages due to the dynamic behavior of the SGS.Based on these known risks, understanding the soil structure of the SGS can help engineers to predict the dynamic behaviors of * Corresponding author.E-mail addresses: mhg@et.aau.dk(M.Ghomi), ffs@et.aau.dk(F.Faria da Silva), shayegani@ut.ac.ir (A.A. Shayegani Akmal), clb@et.aau.dk(C.Leth Bak).
the grounding system and induced voltage values.The used methods provide results with insufficient accuracy in high-frequency situations.Traditionally, the most widely known time-domain simulation tools are applied for transient assessment.Although, in these software packages, the frequency-dependent response of the SGS is not taken into account.It is modeled as a lumped parameter.The quasi-static theory is used to calculate the impedance of SGSs.These models fail to yield reliable predictions if used to predict overvoltages at high-frequency events.With considering a uniform soil structure, Alipio et al. are analyzed the multi-terminal response of the large SGS, including the frequency dependence of soil electrical parameters, and are presented the differences between the low frequency and high-frequency behavior of grounding systems [7,8].Although, there is no proper equivalent circuit for modeling the multi-terminal two-layer SGS considering soil frequency dependency effects [9][10][11].
When constructing substations, most horizontal electrodes are placed in the upper layer, whereas the vertical electrodes are placed in a combination of two or more layers of the soil.This justifies the need for a multi-layer soil grounding system approach.There is no https://doi.org/10.1016/j.epsr.2022.108139Received 6 September 2021; Received in revised form 19 April 2022; Accepted 24 May 2022 uniformity in the soil in practice.Different numerical models have been developed to model a nonuniform soil using horizontal layers characterized by different soil electrical parameters.Based on IEEE standard recommendation, two-layer soil configurations should be accounted for when estimating the transient behavior of grounding systems in the non-homogeneous soil structures with more accuracy is required [1,2].Hence, an accurate assessment of the transient performance of an integrated power system needs its corresponding grounding system to be suitably included in the investigation as well.Nevertheless, in the transient overvoltages and insulation coordination studies, the quasistatic assumption of the grounding system can lead to unreasonable results while inexact soil structure has been used (i.e., using the uniform soil when soil is layered) at high-frequency ranges [12,13].
To model the GS and evaluation of the performance of the grounding grid, the TLM method has been discussed in [12,22,23].The TLM technique has a lower accuracy for lightning overvoltage studies in comparison to the full-wave approach [6].The full-wave methods have been proposed in [24][25][26][27][28][29].They did not analyze the performance of the grounding grid, including the layered structure, model grounding system considering the frequency-dependent impact of the soil electrical parameter, and making the link between the developed grounding system model and existing time-domain simulators.In the high-frequency regime, several approximations have been applied in [14,[30][31][32] to calculate the admittance matrix using the impedance matrix inversion.Nevertheless, the employed approximations can decrease the accuracy of the predicted overvoltage values because of the accuracy reduction in the GS model.To this aim and due to the frequency content of a transient impulse, a full-wave approach up to several MHz is highly recommended for grounding grids modeling [6], considering the soil dispersive property (frequency dependence of electrical parameters, (.., resistivity and permittivity)) [33][34][35][36] and the corresponding overvoltages.Using these models, the impact of frequency dependence of soil electrical parameters on the lightning response of SGS can be evaluated.Simulations are carried out in the frequency range of DC to 10 MHz considering the frequency dependence of soil electrical parameters using the causal model presented by  and   [37,38].This modeling should be well suited for time-domain simulation to analyze transients of power grids.
This paper is divided into three parts.
(1) A full-wave appraisal of SGSs up to several MHz: to evaluate the SGS, a full-wave approach entails solving Maxwell's equations numerically.Among the most extensively studied and improved electromagnetic methods in recent years, the MoM solution solves Maxwell's equations for complex geometries, propagation of electromagnetic waves, layered, and lossy soil structures.This method has accurate results over a wide frequency range and computational advantages when used properly.A thin wire approximation of the SGS can be applied to obtain the harmonic impedance.It is assumed that soil parameters depend on frequency [33][34][35][36]39]; (2) Based on frequency dependence representation, a time-domain representation of layered SGS is derived.The calculated rationally fitted HIM is used to produce a multi-terminal SGS based on a layered soil formation in the time domain.The vector fitting approach (VFA) is then used for calculating a rational expansion of self, and mutual impedances [29,40,41].Once the rational fitting is computed, it could be easily transformed into a time-domain state-space representation that is, in turn, suitable to be implemented in EMT-tools [42][43][44].The number of iterations required by VFA was significantly decreased by selecting an appropriate set of starting poles.
(3) computing transient overvoltages on substation equipment: in this part, it is vital to use high-frequency model of the substation equipment in the time domain.Lightning and switching studies are carried out to predict maximum overvoltage and GPR values.The studies are performed in the EMTP-RV platform [43].
A myriad of approaches have been put forward by researchers for grounding system simulating in recent years [12,45,46].Refs.[45,46] present studies that explain the full-wave modeling of a large grounding system.However, the presented results do not apply to the layered soil structure.It is adopted just for two-terminal grounding systems, which cannot represent the transient overvoltage for multi-terminal grounding systems [45], as demonstrated in Section 2 of this paper.By employing the interpolation procedure described in [12], the Sommerfeld integrals [47] can be derived quickly via reducing the number of direct computations while preserving the exactness of the full-wave approach.It has to be noted that significant research has been dedicated to obtaining transient overvoltage of large SGS considering the dynamic behavior of the grounding system.To the extent of our knowledge, the influence of soil frequency-dependency on transient overvoltages has not been addressed in these works for the complicated, layered grounding systems.There were research [46] regarding transient overvoltages caused by lightning in the MV substation; however, it: • did not take into account the impact of layered soil structure on HIM; • did not incorporate the frequency-dependent of soil parameters effect and by not considering it, yields conservative results; • did not take into account multi-terminal SGS.
As far as we know, no comprehensive, accurate method, except recent works presented by Markovski, Grcev [12], Popov [10], and Manunza [48] has been proposed yet for the inclusion of the fullwave model of SGSs into time-domain simulation tools considering multi-layer soil structure precisely.The proposed approach assumes that the soil geometry is identified.This means that the soil electrical parameters and depths of layers can be determined or measured.
The following is the structure of this paper.The used SGS modeling procedure for layered soil structure is described curtly in Section 2. The state-space model development is presented by incorporating the frequency response of the SGS in time-domain simulation in Section 3. Section 4 presents the process of system study and modeling, including problem specifications, substation equipment modeling, HIM calculation, and the frequency-dependent soil model.The study system performed is based on the realistic 63/20 kV substation.Sections 5-7 deal with the obtained results, analytical discussion, and conclusion.Finally, in Section 8, the future works are explained.

Modeling procedure of substation grounding systems
In this part, the SGS procedure of the layered grounding system is reviewed briefly.More details on MoM solutions can be found in [12][13][14][15].Fig. 2 illustrates the full-wave method to calculate the direct characterization HIM and transient overvoltage of SGSs by using the MoM.
Generally speaking, modeling of grounding systems as an electromagnetic problem can be explained by using the well-known Maxwell equations with scalar and magnetic vector potentials.A full-wave approach (electromagnetic field theory) can investigate the performance of grounding systems over the desired frequency interval.To model the grounding system, the same procedure presented in [47,49,50] is adopted.
To solve the mixed-potential integral equation (MPIE) numerically, the method of moment developed by Harrington [15] is employed.Several studies have proved that MPIE is well-adapted for investigating multi-layer media [12].Based on the MoM solution, the grounding conductors of SGS are considered as a thin wire.All electrodes should be divided into fictitious elements.In this procedure, a widely used excitation model, namely the impressed-current technique presented in [51] is employed.A half sub-sectional current is supposed to be impressed at each excitation terminal.It is presumed that SGS is a multi-terminal system.The MPIE is achieved by meeting the boundary condition of the conductors of the SGS on their surfaces given by ( 1) where   and  i are scattered and incident fields by an external source, respectively.
To present the expression of Green's function, the scalar and magnetic vector potentials are used [49,50,52].The current distribution through the numerical solution determines the electromagnetic interactions between each segment.The radius of a conductor is supposed to be shorter than the wavelength.Generally, the current passes in the direction of the axis of the conductor.The endpoint of the electrode segment, the triangular shape of basis functions, should be defined as a zero value.With this in mind, solving the MPIE electrodes buried in the layered soil structure calculates the current distribution of the electrode segments of SGS.The axial current creates an electric field around the medium in this condition.It is computed using the current times Green's function, and the integral of an obtained function is solved over the length of the conductors.To fulfill the boundary condition, assuming a perfect electric conductor, the summation of the tangential electric field must be zero.After Green's function calculation of the SGS, Sommerfeld's integrals are employed to consider the effect of the soil layers.They are solved numerically to assess the spatial domain of Green's functions.Finally, the harmonic impedances of the SGSs are computed from the MoM's harmonic impedance matrix [12].

Grounding system interface
In the transient softwares such as PSCAD/EMTDC, ATP-EMTP, and EMTP-RV, there is no comprehensive model for describing the multilayer multi-terminal SGS over the wide range of frequency interest.After computing the HIM of the SGS, to obtain the time-domain overvoltage at the ground terminal, across the bushing of PTs, or to do insulation coordination, it is necessary to convert the frequency responses so they should be linked with electromagnetic transient tools.
To perform this, the VFA [40,41] which can be used for each element of HIM in an iterative manner, is employed to make an interface between obtained SGS impedance/admittance matrix with the transient software packages.
First, consider the model of pole-residue of the HIM that is expressed in (2).
where   is passive poles.  is the matrix of residues.Their values are real or complex conjugate pairs of HIM elements.The element value of matrix  and  are real. poles are employed for fitting the vectors in the VFA.The size of matrix , , and  are respectively,  × ,  × 1, and 1 × .The poles (M), multiplying by the terminals (T), represent the total number of grounding system state variables (N = M×T).The passivity enforcement is performed to modify all or some of the parameters like   ,   , and  ( is initially set to zero) of () to guarantee that all eigenvalues are positive.The impedance in the frequency domain, given at  frequency sample points, can be explained in the state space expression given by (3).

𝐙(𝑠)
By using a recursive convolution method, the voltage v(t) from the injected current i(t) at each terminal of the grounding system can be obtained based on a fixed time-step  [30].The state-space model of the multi-terminal SGS is linked with an electromagnetic transient tools like EMTP-RV.More details of the method used can be found in [40].
In the time domain, () can be expressed by a state-space model.It is a first-order matrix differential equation and is shown in (4).
where () is defined as a state vector.Respectively, the () and () are input (current) and output (voltage) vectors.The first and second equations are called state and output equations, respectively.The improved matrix of Jordan-canonical equations of the state-space framework is expressed by (5).
It is supposed that the matrix of   ,   , and   consists of only the real poles and their residues values.The matrix of   ,   , and   have residues and their complex poles.Thus, the parameters are precisely identified, so the fitting is nearly perfect (RMS error is 2E-14).The SGS state-space model can then be developed in the time domain software (for example, state-space block in EMTP-RV or FDNE block in PSCAD/EMTDC).

Table 1
Low frequency values of the soil resistivity.

Study system
To assess the impression of modeling of multiterminal layered SGS on the transient performance of power systems, a realistic data of the   substation, which is a 63/20 kV operated by the Iran's TSO (see Fig. 1(a)).The conductors of the grounding grid are constructed from copper with a 15-mm diameter buried in soil with different soil configurations at a depth of 0.6 m.The thickness of first layer is set to 1-m.The grounding system is simulated as an eight-terminal network whose main elements of the substation are jointed to these terminals (see Fig. 1(b)).Appropriate models of PT and lightning arrester (highfrequency models) are used in the study as well [53].To this aim, we focus on the transient generated overvoltages at both sides of PT (T1) as well as its neutral terminal.The effects of the transformer's neutral opening on the transient voltage at the neutral terminal at 63 kV are analyzed.
As presented in Fig. 1, the 63 kV double circuit PTL is connected to the two PTs (Ynd11) of MV substation.The total distance of the power transmission line (PTL) is 50 km.The incoming tower is connected to terminal 1 of the SGS as shown in Fig. 1(a).PTs' high voltage (HV) side is connected solidly to the ground terminals in 2 and 6.Both PTs' low voltage (LV) side is connected to the earthing transformers (ETs).The ETs are connected to ground terminals (GTs) 3 and 7 respectively for T1 and T2.Fig. 1(b) shows the grounding layout substation, which covers an area of 4800 m 2 .The cross-section of all electrodes is 150 mm 2 .For the SGS model, two different models, low frequency (LF) and the full-wave, are considered for both uniform (Case#1) and two-layer (Case#2,3) soil structures.The value of low-frequency (LF) resistivity of soil for two configurations (three cases) is shown in Table 1.

Frequency dependent model of the electrical parameters of soil
Generally, SGSs are ordinarily stimulated by a very simple or nonlinear resistor in transient studies.It should be pointed out that the ground impedance pertains significantly to soil resistivity and permittivity.Furthermore, neglecting this effect may lead to meaningful errors in the harmonic impedance and GPR values, especially for soil with high resistivity.For modeling of the soil frequency dependent (FD) effect, Alipio and Visacro have developed a causal model to explain the frequency dependence of soil electrical parameters [37].To consider this impact on the simulations, according to the recommendation of  [38], these causal models of resistivity and permittivity are utilized here that is given by ( 6) and ( 7) as below: where a frequency is  and   is the resistivity of soil at 100 Hz (LF range).  and  can be computed at each frequency.
The FD impact decreases the amplitude of the GS harmonic impedance and consequently, the generated overvoltages caused by lightning stroke [54].It means that the soil electrical parameter frequency dependence is in charge for reducing the GPR value of substation grounding system [34,36,39].

Harmonic impedance matrix
The integral equation can be changed into linear equations, whose solution can be provided by the MoM.The current is injected into the terminals of SGS, and the voltages between each injected point and a remote terminal are obtained by integrating the electric field on the specified path.Having obtained I  and V  at each terminal, the HIM is computed for a given grounding system.The computed HIM is, in fact, the frequency domain response of the SGS.It is worth mentioning that in the quasi-static methods, the voltage between two terminals is pathindependent [55].The computed voltage depends on the integration path at the high-frequency range due to the lossy ground situation.In this paper, the scalar potential is employed to deal with the issue of the path dependence effect on the voltage and consequently calculated HIM.The HIM elements (self and mutual) can be calculated by where V  ( ) and I  ( ) are electric potential phasors at the injected terminal in reference to the remote terminal point and impressed current at each frequency, respectively.Fig. 3 exhibits the multi-terminal grounding system, which is based on the layered soil structure of SGS.Once each terminal is excited by a current source, the other terminals are assumed as an open circuit.HIM dimensions for two cases according to the number of terminals (T) is 8 × 8.The absolute value and phase angle of the self and manual impedances for a few terminals are presented (see Figs.

𝐙(𝑠)
where  is the Laplace variable, it is also shown that the self impedance is dependent on the soil resistivity at the lower frequencies.In contrast, the self impedance's well-known inductive behavior for both cases increases at higher frequencies up to 10 MHz.This behavior is remarkably dissimilar from the grounding system's LF model, which misses presenting an accurate frequency response at high frequencies.The value of mutual impedance at high frequency depends very much on the relative positioning of the GTs and soil resistivity.It can also be observed that the self-impedances are different at each terminal.Fig. 4 illustrates the resulting frequency-domain self and mutual impedance (absolute values and the phase angles) considering two different soil structures for the SGS.Assume that the soil in Case#1 is uniform with a resistivity of  1 =  2 =100 Ω m.In Case#2 and 3, as shown in Fig. 1, the grounding grid is buried in two layers of stratified soil.In Case#2, the LF soil resistivity in the upper and lower soil layers are, respectively,  1 =100 Ω m and  2 =10 Ω m, respectively.In Case#3, soil is characterized by resistivities of  1 =100 Ω m. and  2 =1000 Ω m, respectively, for upper and lower layers.The resistivity values are used in ( 6) and (7) as   when the MoM approach is performed.The LF value of the soil relative permittivity is assumed 10 for all cases.
It can be seen that in the early-frequency region up to the frequency of about 100-200 kHz, the harmonic impedance is not frequencydependent and equal to the low-frequency resistance to ground.It depends on the soil resistivity and configuration of the grounding grid.In the high-frequency (HF) range, it is frequency-dependent, and its value becomes higher than LF ones.The harmonic impedance demonstrates inductive behavior at higher frequencies.
As seen from Fig. 4, to verify the presented approach, the self and mutual impedance computed by the MoM are in acceptable agreement with those calculated using FEM [56].In this simulation, the amount of memory needed for both methods reached of about 3 TB.The number of segments being analyzed in the MoM is 21,793.A singlefrequency simulation was completed in about 13 min on an Intel-Xeon E2286 G Processor (12M Cache, 4.00 GHz) desktop computer.The computational time in the FEM is about 54 min at each frequency.The used frequency samples are 129.

Substation equipment modeling
The substation's single line diagram (SLD) of the installed main types of equipment is illustrated in Fig. 5.A simple capacitive model is selected for modeling the capacitance voltage transformer (CVT), disconnect switch (DS), and circuit breaker (CB).The values of    ,   , and   are 6.6 pF, 60 nF, and 60 nF respectively based of manufacture data sheets [57].

Power transformer modeling
For the PT modeling, the  model with the parameters obtained from the datasheet of manufacture is used [58].It requires very detailed data about the transformer due to inductive and capacitive impacts from the windings, core, and tank.As well as this, a linear model proposed by Gustavsson can be used in EMTP-type programs to calculate a black-box model based on measurements at the transformer terminals [59].It is noted that this model can be used for the HF range of interest [53].The amounts of the parallel capacitance between the ground and HV terminals and (  ), the ground LV terminals (  ), and series capacitance between the HV and LV terminals of the power transformer (  ) are listed in Table 2.

Lightning arrester modeling
It is commonly known that lightning arresters (LAs) are among the most important equipment in the protection and insulation coordination of substations.Two sets of LA are installed at the substation entrance (GTs 1 and 5).Also, both sides of PTs are protected by LAs connected to the ground terminals 2 and 3 for T1 and 6 and 7 for T2, respectively.Due to the frequency-dependent behavior of the LAs, an accurate representation is necessary to duplicate arrester efficiency in high-frequency transients.To model a LA in EMTP-RV, the proposed procedure by Martine is used [60].  and   are nonlinear resistors.They can be selected based on LA's impulse residual voltages value with the rise time of 45 μs.A Typical zinc oxide (ZnO) lightning arrester is modeled.The maximum continuous operating voltages (MCOV) for HV and LV sides are 60 kV and 14.4 kV, respectively.Fig. 6 shows the model of LA implemented in the time domain simulation.The physical data (overall height, block diameter, number of columns) and its V-I characteristic are used to determine the model parameters obtained from the manufacturer's datasheet.These parameters are calculated from the given formulas by ( where  are  are the LA height (in meters) and the number of columns of ZnO plates installed in parallel.In this simulation,  is 2,   and   are 1-m and 1.4-m respectively.

Tower modeling
For the modeling of the PTL tower, the proposed multi-story model is used in this analysis according to the well-known analytical formula presented in [61].The tower height is 24.45-m.The surge impedance can be calculated in the time domain, having reached the tower geometry.Tower consists of four sections.Each section modeled by a distributed parameter of transmission line   .It is connected series with a parallel  circuit model.The parallel  circuit represents distortion and traveling-wave attenuation of the tower.The   1 ,   2 ,   3 , and   4 are the surge impedance from top of tower to the higher phase arm ( 1 ), higher to middle ( 2 ), middle to lower ( 3 ), and lower phase arm to tower bottom ( 4 ), respectively.It is assumed that the   1 −   3 are the same and set to 220 Ω and   4 is set to 150 Ω [61]. ) where  is defined as a tower height, the  = ∕ and =0.89 are traveling time and attenuation along with the tower, respectively.=300 m∕μs is a free space speed of light.The average value of measured tower-footing impedance at low frequency for all towers is around 8 Ω.

Power transmission line modeling
The modeling of PTLs is implemented by distributed lines with the FD parameters in the EMTP-RV platform.The incoming 63-kV line is modeled as a multi-phase untransposed model, implemented by three spans (236, 276, 260 m).To reduce the impacts of wave reflection from the line termination at other substations, a long span of 49.3 km is added.

Time domain model verification at the power frequency
To further assess the exactness of the presented model, short-circuit generated overvoltages are also investigated in the time domain simulation.To this aim, a phase to ground fault (L-G) has happened in the LV side of substation (phase-a) at =10 ms.The fault duration time is 20 ms.The 3-phase voltage waveforms of both LV (left) and HV (right) sides due to a fault (L-G) computed using the full-wave (solid line), and low frequency (dotted line ) approaches are illustrated in Fig. 7.As shown, the voltages of the other two non-fault phases (phases B and C) increase by a factor of k= √ 3 due to the type of transformer's connection.The results show that the developed full-wave model of the SGS is in good agreement with the low-frequency model in the time domain simulation at the power frequency.It is also worth mentioning that the calculated results are identical with EMT tools (LF model), and the proposed model is also verified to analyze SGS buried in two-layer soil.This is due to the L-G fault's frequency content for which the two SGS modeling approaches take the same values leading to identical overvoltages, and it is predictable.Indeed, the frequency content of the short circuit is placed in the low-frequency range (power frequency).The grounding substation impedance response (self and mutual impedance) of the full-wave and LF model are the same at this range of frequency based on obtained results in the time domain (see Fig. 7).A comparison of the resulting 3-phase waveforms in time domain simulation validates the model's accuracy.

Simulation study results
This part of the paper presents the results of the transient overvoltages for lightning and switching impulses.It is noted that the SGS is simulated following the procedures shown in Section 2. Different cases are illustrated with different soil resistance values and soil structures (uniform and two-layer).
Case#1 refers to a uniform soil, whereas Case#2 and 3 correspond to the two-layer soil structure with a depth of 1-m for the first layer.The transient overvoltage for the uniform and two-layer soil structure is compared in different scenarios.The simulation results are evaluated when (1) the neutral of T1 is solidly grounded and ( 2) it is opened.The neutral power transformer T2 is connected to grounding terminal 6 permanently.For each ground structure (uniform and two-layer soil), the maximum transient overvoltage value on each phase and GPR value on the grounding terminal of PT T1 are computed.

Model of the lightning current
To compute the lightning overvoltage and GPR, the parameters of Heidler's functions are adopted to approximate the subsequent stroke (SS) currents from the obtained empirical data by Berger et al. [62][63][64].
In accordance with Heidler's functions, the SS has a front time, steepness, and peak value of 0.8 μs, 40 kA/, and 12 kA.More precisely, the fastest front time of the lightning waveform in the time domain produces in high-frequency contents (more than 200 kHz) in the frequency domain.Having enough high-frequency content of subsequent lightning strokes might cause a large transient overvoltage at the injection point due to the inductive behavior of SGS at higher frequencies [6].This is the main reason for considering the SS in this work.

Substation ground potential rise
Lightning current leads to generate a GPR in the MV and LV equipment and it is one of the major source of overvoltages in substations.The GPR is computed by integrating the electric field along a predetermined path.The soil structure can be analyzed to determine its impact on GPR and better explain the significance of the SGS soil dispersive behavior.Fig. 8(a) shows the comparison of the developed GPR using full-wave and LF (or resistive) models for uniform soil structure of SGS.The obtained results show about 30 kV voltage difference between fullwave and LF models, which should be considered in the insulation coordination of PT.
For Case#2, two-layer soil, the peak value of GPR at the neutral point of the HV side of PT T1 (terminal 2) is presented in Fig. 8(b).It observes that the peak of GPR, which is calculated using the fullwave approach, is higher than when the model of grounding system is assumed as a low frequency (LF) model.This is owing to the fact that the GPR takes different peak values at higher frequencies than in comparison with the LF model.These values are about 10 kV when the neutral of PT T1 is grounded solidly or opened in the post-fault condition around  = 2 μs.It can be seen from Fig. 8(c) that the peak values obtained from LF and Full-wave model are about 70 kV and 120 kV when the neutral of PT T1 is grounded solidly and opened, respectively (for Case#3).
It is evident from the results that because the dynamic behavior of the SGS at high frequencies is ignored in the LF model, the GPR might be underestimated.For instance, when the neutral of PT T1 is connected to terminal 2 of the substation grounding grid, the differences of the voltage peak are approximately 70% in (Case#1), 80% in the (Case#2), and 40% in the (Case#3).Once the ground terminal opens, the voltage peak value differences are 66%, 20%, and 44% for (Case#1), (Case#2), and (Case#3) respectively.The results demonstrate that the impact of the full-wave model of SGS at high frequencies is the source of the considerable variation in the GPRs peak values.It can be observed that the effect of the full-wave method and layered soil structure, determined from the proposed solution, should be considered in the GPR calculations.Based on the full-wave model of SGS used in this study, GPR peak values are substantially higher than LF ones, which significantly impacts simulation results.

FD impact of the soil electrical parameters on the GPR
To further investigate, the computed GPR curves for the identical SGS subjected to the subsequent lightning stroke are exhibited in Fig. 9 that highlights the influence of the FD of the soil electrical parameters on the GS lightning performance.As we discussed earlier in IV-A, the soil electrical parameters are a function of frequency.For assessing the impact of soil electrical parameters frequency-dependency on the GPR value, the soil electrical parameters are assumed to be a constant (CP) and FD.
Fig. 9 illustrates the GPR values of SGS for Case#1 and Case#3.This effect is clearly more pronounced for GSs buried in the high resistive soil when exposed to a subsequent lightning strike.It is obtained that the effect of soil electrical parameter frequency dependence is in charge of reducing the GPR of the SGS.The outcomes demonstrate that the influence of the frequency dependence of soil on the GPR is ignoble in Case#1 with low resistivity of the soil of about 100 Ω m.In this regard, neglecting the FD effect cannot cause an error in the GPR values.If the soil resistivity is low, the constant parameter model is applicable to reduce computation time (see Fig. 9(a)).It can be a significant error of the GPR value in Case#3 considering the CP model.In fact, ignoring the frequency dependence effect results in conservative GPR values due to the higher resistivity of soil (see Fig. 9(b)).Moreover, a comparative study is performed to demonstrate the importance of including frequency-dependent soil parameters.This is important in terms of the insulation strength of transformer winding since the winding close to the neutral point is generally more prone to overvoltages.By calculating GPR and electrical stresses on the substation PTs, engineers can perform parametric analysis, determining whether to use FD or CP models.1).

Lightning overvoltages
The proposed technique is utilized for the calculation of lightning overvoltage.It is assumed that the SS occurs on phase A of the incoming 63 kV of the double-circuit tower at a 50-m distance from the transformer.
Fig. 10 shows the computed overvoltages at transformer HV (right) and LV (left) sides for all phases when the full-wave and LF (resistive) modeling approaches of the uniform SGS are accomplished (Case#1).It is clearly seen from this figure that the full-wave modeling approach generates more severe overvoltages than the LF approach.
Figs. 11 and 12 show the generated overvoltage of SGS buried in two layer-soil structure as given in the (Case#2) and (Case#3).In this part, the lightning studies of the substation in the proposed model revealed the importance of high-frequency modeling of two-layer SGSs LV side (Left) and HV side (right), (see (Case#2) in Table 1).Fig. 12. Generated overvoltage when subsequent lightning strikes the phase conductor (phase A) at shielding failure situation considering two-layer soil structure.(a) The neutral of T1 is grounded solidly (state 1), (b) the neutral of T1 is opened (state 0).LV side (Left) and HV side (right), (see (Case#3) in Table 1).
for the fast transient overvoltage.This is owing to the fact that the lightning performance of a two-layer grounding system is entirely different from when SGS is assumed uniform as a simplified model.It is clearly shown from this figure that the full-wave approach predicts more severe overvoltages than the LF method.In the LF model, ignoring the inductive behavior of SGS results in underestimating values of the fast transient overvoltage at PT's bushings.The influence of the soil layer on the generated overvoltage is demonstrated.In general, the soil layer assumption in the model of SGS can change the final results of overvoltage values.

Switching overvoltages
Switching overvoltages in power systems is usually a matter of concern in insulation coordination studies since they can damage insulation or cause insulation flashover.This transient can occur when PTLs  1).
are energized and de-energized or during fault clearing, with a wave traveling along with the PTL following switching.This phenomenon of a switching overvoltage is, thus, a matter of a few microseconds with a range of frequency lower than lightning (5-150 kHz).It is vital to limit the amplitude of the transient overvoltage to avoid exceeding the surge rating on this substation equipment which is not specified for IEC in ranges below 220 kV [65].In this section, to better understand the impact of the full-wave modeling of the layered SGS, the generated switching transient overvoltages considering uniform and two-layer soil configuration of SGSs are compared.For the modeling of the switching situation, the switching operation is executed at the time,  = 10 μs.Fig. 13 shows switching overvoltages at all phases (both LV and HV sides) when SGS is considered a uniform soil structure.
To evaluate the effects of the type of grounding connection of PT, the simulation results are analyzed while the neutral of power transformer T1 is solidly grounded (see Fig. 13(a)).Then, overvoltage values are computed when the neutral of PT T1 is opened (see Fig. 13(b)).
Figs. 14 and 15 show the switching overvoltage for all phases (both LV and HV sides) when the two-layer soil structure of the substation is taken into account.These are also calculated with (see Figs. 14(a) and 15(a)) and without (see Figs. 14(b) and 15(b)) jointing the neutral point of T1 to the grounding terminal.These figures demonstrate the overvoltage values when the model of SGS is full-wave, differs from LF ones but less than for lightning which is expected given the lower frequency.Indeed, the maximum frequency content of the switching transients is about 150 kHz.Grounding substation harmonic impedance responses do not have considerable differences at low-frequency range up to 150 kHz for all cases.

Discussion
This section evaluates the impact of full-wave modeling of layered SGS on the switching and lightning performance of MV substations.The effectiveness of soil layers and structures in the MV SGS on HIM has been studied.Transient overvoltages have been investigated, considering frequency-dependent influences from layered soil.The choice of a proper model of SGS is an essential part of transient studies.The grounding system has been modeled based on several approaches from quasi-static (LF) to exact (full wave) method.The level of soil structure details should be considered according to the best agreement between the complexity and correctness of the case studies.The factors like the  1).Fig. 15.Generated overvoltage considering two-layer soil structure when switching occurs the phase conductor (phase A).(a) The neutral of T1 is grounded solidly (state 1), (b) the neutral of T1 is opened (state 0).LV side (Left) and HV side (right), (see (Case#3) in Table 1).impact of soil structure of SGSs and connection state of power transformer neutral on transient overvoltages have been examined.Two soil structures of the grounding system were used to investigate those effects on transient overvoltage estimations.In the studied cases, simulation results indicate that LF models may underestimate/overestimate potentially transient overvoltage on the power transformer (T1), and equipment of substation may damage or increase the cost of insulation due to these underestimation/overestimations.To extra assess the impression of the frequency dependence of SGS responses, an analytical comparison is performed on the calculated overvoltage values based on Table 3.This table shows the voltage peak values resulting from lightning and switching events with different states of the neutral connection to the ground for transformer T1.Table 3 indicates that the full-wave modeling of the SGS could significantly affect overvoltage values compared to the LF model.The results provide information that can be considered in studies of overvoltage protection of substations  considering soil structure and frequency-dependent model of SGS.The present results are not general but specific for the case study analyzed in this research, and further works and cases will be needed.For these purposes, the following two indexes are defined.

Low frequency error index (LFEI)
The LFEI criteria are computed for all installed power transformer T1 phases in the substation.It is assumed that the SGS buried in the soil with the same structure.As shown in Table 4, the full-wave modeling of SGS impacts generated lightning and switching overvoltages.This table illustrates the calculated overvoltages and demonstrates how ignoring the frequency dependence of the SGS model results in mostly an underestimation, with some cases of overestimation for lightning and switching values.

Soil structure impact index (SSII)
The SSII for all phases of the installed power transformer T1 with the SGS buried in the different soils structures (uniform and two-layer) are calculated.

𝑆𝑆𝐼𝐼 =
max{Overvoltage value (uniform soil)} max{Overvoltage value (two layer soil)} According to Tables 5 and 6, ignoring the soil structure and using the LF model of SGS might lead to some errors in the estimated overvoltage.LF model or simple soil structure assumptions of SGS in the high frequency events may lead to a remarkable underestimation/overestimation in the value of transient overvoltage on the phase conductors of PT as well as of the GPR.

Conclusion
This paper can be regarded as an extension of the researches proposed by Markovski, Popov, Sheshyekani, and Legrand, in which the same approach was adopted [10,12,14,45], and [49].The impact of full-wave SGS modeling on the lightning/switching performance of the MV substation has been examined.The chief contribution of this paper is the transient evaluation of the two-layer multi-terminal soil structure of the large SGS, taking into account the frequency dependence model of the electrical parameters of soil, which is not generally investigated simultaneously for the real MV substation.Accordingly, the ignoring FD effect results in conservative GPR values due to the higher resistivity of soil.Moreover, a comparative study has been performed to demonstrate the importance of including frequency-dependent soil parameters.However, in these situations, the FD of the soil electrical parameters ought to be considered to accurately assess the SGS performance:1) for GSs buried in soils with high resistivity and 2) to investigate the precise behavior of GSs exposed to a fast front waveform with high-frequency content.
It has been demonstrated that utilizing the full-wave GS modeling provides a more accurate assessment of the lightning overvoltages than the quasi-static methods.Accordingly, applying a simple (LF) model of the GS may yield an inaccurate estimation of the lightning overvoltages on the apparatus bushings joined to the phase conductors or neutral points of the installations.The MoM solution was selected to govern MPIE.Then, the HIM of SGS was computed over the required frequency interval.The proposed approach analyzed the MV substation with a minimum 350 kV BIL.The subsequent stroke was modeled.In order to use an exact model of the SGS that is matchable with an EMT software, a state-space model was employed to make an interface between the proposed approach of the SGS and the time domain simulation.This model defined the correlation between voltages and currents at each terminal of the multi-layer multi-terminal grounding systems.To this aim, the pole-residue model provided by VFA was obtained for two cases.The precision of the full-wave technique was substantiated by comparison with FEM.The accuracy is validated both in the frequency and time domain.
From the results shown in this paper, it has been demonstrated that the harmonic impedance of the SGS depends on the soil structure and resistivity plays a vital role in the definition of the accurate transient overvoltage values and GPRs.It was shown in particular the effects of the full-wave method and layered soil structure, determined from the proposed solution, should be considered in the GPR calculations.Disregarding the dynamic behavior of grounding grids or approximating it by LF values could cause erroneous results despite that the substation equipment above the ground has been simulated within a wide range of frequencies.
The comparative results of the computed transient (lightning and switching) overvoltages at the transformer bushings for different soil structures and the type of transformer's connection have been summarized in Tables 3-6.It has been observed that the soil structure and resistivity of the soil may lead to influence the overvoltage values compared to the LF model.

Future work
The authors believe that further work is required to reduce the calculation time for the MoM, which is time-consuming, especially on large SGS with several electrodes and terminals considering multi-layer soil structures.Also, more accurate models for equipment like power transformers with more details will be applied to compute results based on the desired precision.The effective length of the first layer depth of soil in multi-layer structures on the input impedance will be researched in future works.It can help reduce the SGS complexity and computation time from an engineering point of view.In addition, the FS current will be taken into account to compute transient overvoltage, which has a higher magnitude than the SS.

Fig. 5 .
Fig. 5. Single line diagram of the MV substation and equipment arrangement.

Fig. 7 .
Fig. 7. Comparison the LF and MoM models on the over voltage generated at the power transformer bushings when single phase to ground fault occurs at LV outgoing distribution feeder (Phase a).[(a): LV side and (b): HV side (right)]; PT (T1) neutral point is connected to ground solidly.

Table 2
-circuit capacitance values for PTs adapted from the data sheet of the supplier.
Fig. 6.Fast front model of lightning arrester.

Table 3
Overvoltage peak values due to lightning (subsequent stroke current) and switching events.

Table 4
LFEI for the overvoltage values shown in Table3.