Grounding System Modeling and Evaluation Using Integrated Circuit Based Fast Relaxed Vector Fitting Approach, Considering Soil Ionization

: Since high voltage transmission line towers or wind turbines structures are installed in high-altitude areas, it is essential to achieve a high overvoltage protection system against direct and indirect lightning strikes collisions. The lightning current must be discharged quickly into the protective earth, to prevent the dangerous over-voltages formation and deﬁne a reference voltage node. This paper presents a novel model to assess the behavior of the grounding system, based on Pocklington integral equations under lightning magnetic ﬁelds and variations in soil ionization, in which an explicit circuit-based vector ﬁtting RLC admittance branches are proposed. The frequency-dependent behavior of grounding system frequency response and soil ionization e ﬀ ect is modeled in time domain, straightly to implement into the electro-magnetic transient program (EMTP). The model veriﬁcation contains horizontal, vertical, and their combinations of grounding grids to represent the complete investigations under lightning strikes. The harmonic impedance mathematical formulations and principles are derived based on a rational function, that could be applicable on ground potential rise (GPR) investigation.


Introduction
Grounding systems, or the buried electrodes under the soil (either individually or in combination and interconnected with regular arrangements), are essential and significant to provide the health and safety for the staff occupied at high-voltage substations, power plants and industrial areas. The main purpose of grounding systems is to prevent unauthorized step voltage and touch voltage that faces human health at risk. Their other task is to provide the same common reference potential for all electrical and electronic components, especially protective equipment connected to the power grid.
In almost a few years, the development and installation of impeccable lightning protection system for high voltage towers and wind turbine structures, due to their high altitudes and vulnerability to potential lightning strikes, has become one of the fields of interest and study for researchers [1,2].
In all types of electrical grounding systems, these two requirements have to be considered: (1) to provide the safety and health for the people who are in contact with electrical equipment; (2) to help to maintain the integrity of electrical equipment and increase their toughness to ensure performance, regardless of safety issues [3,4]. In the first view, the principle of safety is for people in the community who are the end consumers of electrical energy, to prevent firing, equipment short circuit (SC) breakdown, and to stopover SC current. This kind of strategy is called protective earth, the accuracy of computations. In some literature, the circuit-based models have been proposed without considering the effect of soil ionization, or some negative resistance parameters have been calculated which are physically unacceptable [27,28].
In high frequency analysis, the soil ionization causes non-linearization of the transient behavior of the buried grounding system in the Earth. Currently, if these equipments are buried in multilayer soils with different special resistances and with different geometries such as cylindrical or spherical soils, the modeling of its electrodes, in order to analyze power networks in software packages, is practically impossible. Therefore, this paper presents a novel technique to investigate all these phenomena, in order to investigate the electromagnetic transient states of lightning strikes by fast relaxed vector fitting (FRVF) algorithm. Using this method, each arbitrary grounding system can be implemented and evaluated with just a few simple RLC admittance branches with a low-order state-space transfer function and frequency response in a very wide frequency range. Another feature of this process is its high precision and computational speed.

Problem Formulation
When the lightning strikes collide to a transmission line tower, the current flows to the grounding grid directly. The Maxwell equations integrated with Pocklington equations can model the lightning current flowed to any grounding electrodes or rods [29,30]. These sets of Pocklington-Maxwell equations will be solved and evaluate the high-frequency grounding grid behaviors. The excitation function (E exc sm (s)) in the loss-full medium is derived by the following equation. (2) I n (s ) is the current induced along a horizontal electrode and a vertical rod as shown in Figure 1, and g o (s, s ) represents the half-power green function, while the g i (s, s ) increases according to image theory. Thus, they are found by (3) and (4).
g 0 (s, s * ) = e −jk 1 R 1 R 1 (4) where R 0 and R 1 are the distance between the reference point and the image point to the observer. Furthermore, k 0 and k 1 represent the propagation constant in air and loss-full ground, respectively, in which they are defined with the following formulations.
Appl. Sci. 2020, 10, 5632 4 of 18 where r and σ are Earth permeability and conductivity coefficients. Since ω = 2πf represents the angular frequency, the → G s (s, s ) could be extracted through electric dipoles using Sumerfeld integral law.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 4 of 18 G ⃗ (s, s ) = x . s × G . ρ ⃗ + G . φ ⃗ + G . z ⃗ + z . s × G . ρ ⃗ + G . z ⃗ The unknown current I (ϑ) flowing in n th section of the grounding grid, can be approximated by the sum of finite and linear independent vectors f with weight coefficients as (8).
where the iso-parametric variables result in: Using the equations written in (10) will give us a proper linearization and approximation at border elements.
Then, the impedance matrix is calculated in (11).

Lightning Current Formulation
There are many current functions with mathematical formulations that have been proposed to model the lightning behavior, such as 'Double Exponential', 'Jones Modification', 'Heidler', and The unknown current I e n (ϑ) flowing in nth section of the grounding grid, can be approximated by the sum of finite and linear independent vectors f ni with weight coefficients as (8).
where the iso-parametric variables result in: Using the equations written in (10) will give us a proper linearization and approximation at border elements.
Then, the impedance matrix is calculated in (11). The matrixes {f} and f consist of image functions, while {D} and D represent the derivative operators. The abovementioned equations are solved using numerical analysis.

Lightning Current Formulation
There are many current functions with mathematical formulations that have been proposed to model the lightning behavior, such as 'Double Exponential', 'Jones Modification', 'Heidler', and 'Pulse' functions [31,32]. However, a proper function that can accurately model the lightning strike must have the following characteristics [33]: • Shape the lightning waveform with a very good approximation, as the base shape described in [34]. • It must be able to provide the typical parameters of the lightning waveform, such as the maximum current peak, the maximum current slope, rise time, and settling time.

•
No discontinuity should occur in the first and second derivatives of the current lightning model over time (especially at the instant t = 0, due to magnetic and electric field calculations).

•
The current lightning model must be differentiable, in order to calculate the electric and magnetic fields generated.

•
It shouldn't be complicated to avoid appearing "Dirac delta function" or "doublet function" and their higher-order derivative through mathematical computations.
Therefore, the CIGRE model is used to represent the lightning current behavior in this paper. This single-phase model uses two portions as follows; one is to model the front-wave and the other for the tail-wave current model, where these two parts are fully continuous at the interconnect point.

Front-Wave Model
The front-wave model is expressed with Equation (12). where, The above unknown variables are explained in (15)- (17).
t f and S m represent the front time of the waveform and the maximum slope, respectively.

Tail-Wave Model
The tail-wave lightning current is modeled using (18).
Appl. Sci. 2020, 10, 5632 6 of 18 In Equation (18), the coefficients are defined as: In (19), t h is the time duration that the lightning current, reaches 50% of its maximum value. It is also worth noting that (18) is used when t ≥ t n + t start . The current waveform described is shown in Figure 2.

Parameter Derivation
The grounding system equivalent impedance spectrum, as extracted in (11), cannot be utilized straightly to ATP/EMTP software packages. Therefore, it is essential to convert this transfer function to the passive parameters within all high-frequency soil behaviors.

Fast Relaxed Vector Fitting
If the frequency response of a section of the power system is available as (23), formerly, it could be approximated with several finite rational functions, as in (24).
where 's' is the Laplace operator and and coefficients represent the residuals and fitted poles respectively. Moreover, and ℎ are real numbers which show the constant system impedance. The FRVF method will estimate all these coefficients with high accuracy to determine the rational functions. This approach includes two steps to approximate stable poles and proper residuals. In the first step, a set of stable preliminary poles is considered as an initial guess. Then, the matrix-fitting presented in [35,36] will be used to obtain the deterministic real and complex conjugate fitted poles. At the second step, the residuals are defined by solving a set of linear equations as . = obtained from (24). However, the solutions attained from this approach are acceptable, but when they are compared with the real impedances, some of the resistors may get negative values. Therefore, a modification is needed to apply this criterion to the fitting process. As the matrix fitting is initiated with a set of stable preliminary poles, the round off errors might result in a high variance between the approximated poles and the real poles. Then, it could be solved by using the imaginary part of the poles, as in (26). a = −p ± j q (25) q = 100 p (26)

Parameter Derivation
The grounding system equivalent impedance spectrum, as extracted in (11), cannot be utilized straightly to ATP/EMTP software packages. Therefore, it is essential to convert this transfer function to the passive parameters within all high-frequency soil behaviors.

Fast Relaxed Vector Fitting
If the frequency response of a section of the power system is available as (23), formerly, it could be approximated with several finite rational functions, as in (24).
c n s − a n + d + sh (24) where 's' is the Laplace operator and c n and a n coefficients represent the residuals and fitted poles respectively. Moreover, d and h are real numbers which show the constant system impedance. The FRVF method will estimate all these coefficients with high accuracy to determine the rational functions. This approach includes two steps to approximate stable poles and proper residuals. In the first step, a set of stable preliminary poles is considered as an initial guess. Then, the matrix-fitting presented in [35,36] will be used to obtain the deterministic real and complex conjugate fitted poles. At the second step, the residuals are defined by solving a set of linear equations as A.x = B obtained from (24). However, the solutions attained from this approach are acceptable, but when they are compared with the real impedances, some of the resistors may get negative values. Therefore, a modification is needed to apply this criterion to the fitting process. As the matrix fitting is initiated with a set of stable preliminary poles, the round off errors might result in a high variance between the approximated poles and the real poles. Then, it could be solved by using the imaginary part of the poles, as in (26).
Then, the low value imaginary parts of fitted poles are resolved, and the FRVF method diagram is shown in Figure 3.
Appl. Sci. 2020, 10  However, there are many air-filled pores among the soil components; as shown in Figure 4, the soil model ionization has to be modified to represent the are transitions in air-filled pores, additionally. Therefore, the frequency-dependent consideration of soil conductivity and relative permittivity is not just enough. The model that could represent the full soil behavior is shown in Figure 5, that consists of a paralleled arc resistance and its capacitance series to the CIGRE resistance. The CIGRE model is mathematically formulated in (29) and (30).
where E is the critical electric field intensity of the soil, and R is the ground resistance with low current at low frequency. Consequently, the capacitance of air-filled pores is presented in (31) [39].

Soil Ionization
To illustrate the behavior of high-frequency soil conductivity and relative soil permittivity, the Scott model is habitually considered [37,38]. To determine the mathematical formulation, Equations (27) and (28) should have been substituted with conventional soil conductivity and relative permittivity coefficients.
However, there are many air-filled pores among the soil components; as shown in Figure 4, the soil model ionization has to be modified to represent the are transitions in air-filled pores, additionally.
Appl. Sci. 2020, 10, 5632 8 of 18 Therefore, the frequency-dependent consideration of soil conductivity and relative permittivity is not just enough. The model that could represent the full soil behavior is shown in Figure 5, that consists of a paralleled arc resistance and its capacitance series to the CIGRE resistance. The CIGRE model is mathematically formulated in (29) and (30).
where E 0 is the critical electric field intensity of the soil, and R 0 is the ground resistance with low current at low frequency. Consequently, the capacitance of air-filled pores is presented in (31) [39].
where R arc is the arc resistance that created in the air-filled pores when the magnetic field could not be strengthened by the air. R arc is formulated with Browne's combined theory [40].
where and represent the electrode diameter and length, respectively.
where and represent the electrode diameter and length, respectively.
Appl. Sci. 2020, 10, 5632 where a and l represent the electrode diameter and length, respectively.
where and represent the electrode diameter and length, respectively. Currently, in order to evaluate the high-frequency behavior of the grounding system, along with the effect of soil ionization, Figure 7 shows the circuit-based diagram of the proposed strategy. As is observed, there are some series impedance branches with their real or complex conjugate poles; those are enforced to be stable. It means that the real part of all fitted poles will have a negative value to be located at the left side of the imaginary axis in the complex plane. In (24), however, all the and Currently, in order to evaluate the high-frequency behavior of the grounding system, along with the effect of soil ionization, Figure 7 shows the circuit-based diagram of the proposed strategy. As is observed, there are some series impedance branches with their real or complex conjugate poles; those are enforced to be stable. It means that the real part of all fitted poles will have a negative value to be located at the left side of the imaginary axis in the complex plane. In (24), however, all the c n and a n coefficients could be real or complex, but d and h are certainly real if they exist. The passive parameters of Figure 7 will be calculated with the following algorithm: Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 18 coefficients could be real or complex, but d and h are certainly real if they exist. The passive parameters of Figure 7 will be calculated with the following algorithm: Proposed high frequency grounding system behavior considering soil ionization.

•
The d and h coefficients represent the resistance and the inductance of grounding grid in the rest mode, which are determined by and , respectively.
• If the partial transfer function is defined as (37), then the poles are real and stable, thus we have: Other criteria must be applied to the parameter determination algorithm which prevent the resistance to be negative.
• If in (39) is negative in transfer functions, it could be considered as (40). Then, while the poles are stable, this rational function is simplified to two partial ones. Therefore, the RL branch parameters will be computed with (41) and the newer rest impedance will be calculated as (42). • The d and h coefficients represent the resistance and the inductance of grounding grid in the rest mode, which are determined by R 0 and L 0 respectively.
• If the partial transfer function is defined as (37), then the poles are real and stable, thus we have: Other criteria must be applied to the parameter determination algorithm which prevent the resistance to be negative.

•
If c j in (39) is negative in transfer functions, it could be considered as (40). Then, while the poles are stable, this rational function is simplified to two partial ones. Therefore, the RL branch parameters will be computed with (41) and the newer rest impedance will be calculated as (42).

Model Verification
The proposed formulations and modeling strategy are examined on three different grounding systems, to evaluate their impedance magnitudes and phase angles. The proposed model is also compared and validated by the Grcev model, which is one of the most valuable structures in this study field. In the first scenario, a typical horizontal electrode is considered as a grounding system, and the frequency response of ground impedance is evaluated. In the second scenario, a high-capacitive capability earthling system will be analyzed, and the GPR will be checked. In the third scenario, we evaluate a real grounding grid, in which the effects of soil ionization parameters have been considered. Figure 8 represents a horizontal grounding electrode that is excited by a 20 kA lightning current. As could be observed, the grounding system frequency response has a nearly constant magnitude that occurs at low frequencies. However, since the impedance of the soil increases with increasing the frequency, the frequency response will similarly be proportional to those variations. For low-soil resistivity (100 to 1000 Ω·m), the impedance magnitude increases at high frequencies, which can be justified by inductive behavior modeling of the structure, as shown in Figure 9. Figure 10 represents the phase angles correspondingly that relies on inductive behavior of grounding grids, due to their reachability to 90 • mostly. As the soil resistivity increases, impedance amplitude growth begins at lower frequencies. As a result, in more fertile soils, the frequency response angle becomes faster to achieve to 90 • . Figure 8 represents a horizontal grounding electrode that is excited by a 20 kA lightning current. As could be observed, the grounding system frequency response has a nearly constant magnitude that occurs at low frequencies. However, since the impedance of the soil increases with increasing the frequency, the frequency response will similarly be proportional to those variations. For low-soil resistivity (100 to 1000 Ω. m), the impedance magnitude increases at high frequencies, which can be justified by inductive behavior modeling of the structure, as shown in Figure 9. Figure 10 represents the phase angles correspondingly that relies on inductive behavior of grounding grids, due to their reachability to 90° mostly. As the soil resistivity increases, impedance amplitude growth begins at lower frequencies. As a result, in more fertile soils, the frequency response angle becomes faster to achieve to 90°.    Figure 8 represents a horizontal grounding electrode that is excited by a 20 kA lightning current. As could be observed, the grounding system frequency response has a nearly constant magnitude that occurs at low frequencies. However, since the impedance of the soil increases with increasing the frequency, the frequency response will similarly be proportional to those variations. For low-soil resistivity (100 to 1000 Ω. m), the impedance magnitude increases at high frequencies, which can be justified by inductive behavior modeling of the structure, as shown in Figure 9. Figure 10 represents the phase angles correspondingly that relies on inductive behavior of grounding grids, due to their reachability to 90° mostly. As the soil resistivity increases, impedance amplitude growth begins at lower frequencies. As a result, in more fertile soils, the frequency response angle becomes faster to achieve to 90°.   In the investigation of the optimal length of the horizontal electrodes of the grounding system, the impedance frequency response is calculated for a length of 2 m to 100 m at the same condition. When the same current amplitude is injected into one of the two ends of the electrode, the impedance changes as the soil resistivity changes, which is shown in Figure 11. This plot diagram, which is drawn at a frequency of 50 Hz, reflects the electrode impedance, decreasing since the length is increased. However, as the frequency of the injected current increases, the impedance reduction will stop after a while, and its magnitude will be fixed at a definite constant value. The length of the electrode, when the impedance amplitude begins to be constant, is called the optimal length, and should be used in the grounding grid designs. These results indicate that the lower the frequency, the effective length will be less for a specific reduction. For example, at a frequency of 1 kHz, the effective length for a specific resistance of 10,000 ohms-m is about 1000 m long. However, this effective length can be reduced to 40 m in the soil resistivity of 10 ohms. The point that was previously expected is that, as the frequency increases, the effective length of the horizontal electrodes decreases. In the investigation of the optimal length of the horizontal electrodes of the grounding system, the impedance frequency response is calculated for a length of 2 m to 100 m at the same condition. When the same current amplitude is injected into one of the two ends of the electrode, the impedance changes as the soil resistivity changes, which is shown in Figure 11. This plot diagram, which is drawn at a frequency of 50 Hz, reflects the electrode impedance, decreasing since the length is increased. However, as the frequency of the injected current increases, the impedance reduction will stop after a while, and its magnitude will be fixed at a definite constant value. The length of the electrode, when the impedance amplitude begins to be constant, is called the optimal length, and should be used in the grounding grid designs. These results indicate that the lower the frequency, the effective length will be less for a specific reduction. For example, at a frequency of 1 kHz, the effective length for a specific resistance of 10,000 ohms-m is about 1000 m long. However, this effective length can be reduced to 40 m in the soil resistivity of 10 ohms. The point that was previously expected is that, as the frequency increases, the effective length of the horizontal electrodes decreases. electrode, when the impedance amplitude begins to be constant, is called the optimal length, and should be used in the grounding grid designs. These results indicate that the lower the frequency, the effective length will be less for a specific reduction. For example, at a frequency of 1 kHz, the effective length for a specific resistance of 10,000 ohms-m is about 1000 m long. However, this effective length can be reduced to 40 m in the soil resistivity of 10 ohms. The point that was previously expected is that, as the frequency increases, the effective length of the horizontal electrodes decreases.  Figure 12 shows a relatively simple grounding system that consists of several horizontal and vertical electrodes. The evaluation of this grounding system as a second case study is due to the revision of soil capacitive behavior. The Earth capacitive behavior in the theory of the transmission line is fully justified [17].   Figure 12 shows a relatively simple grounding system that consists of several horizontal and vertical electrodes. The evaluation of this grounding system as a second case study is due to the revision of soil capacitive behavior. The Earth capacitive behavior in the theory of the transmission line is fully justified [17]. electrode, when the impedance amplitude begins to be constant, is called the optimal length, and should be used in the grounding grid designs. These results indicate that the lower the frequency, the effective length will be less for a specific reduction. For example, at a frequency of 1 kHz, the effective length for a specific resistance of 10,000 ohms-m is about 1000 m long. However, this effective length can be reduced to 40 m in the soil resistivity of 10 ohms. The point that was previously expected is that, as the frequency increases, the effective length of the horizontal electrodes decreases.  Figure 12 shows a relatively simple grounding system that consists of several horizontal and vertical electrodes. The evaluation of this grounding system as a second case study is due to the revision of soil capacitive behavior. The Earth capacitive behavior in the theory of the transmission line is fully justified [17].  Currently, the frequency response diagrams of the grounding grid represented in Figures 13 and 14, show that initially, at low frequencies, the system is considered to be capacitive, and the ground impedance is high due to this capacitive property. The phase angle diagram of ground frequency response is also distributed between −90 • to 0 • . As the frequency increases, the grounding system decreases the capacitive property, and enhances the inductive characteristics. This reduces the impedance magnitude, which pushes the phase angle to zero. Then, the higher the frequency, the more indicative of its inductive property, thus, the impedance magnitude of the grounding system will increase. As expected, the frequency response phase angle will reach around 90 degrees.

A Simple Grounding Grid
Appl. Sci. 2020, 10, x FOR PEER REVIEW 12 of 18 Currently, the frequency response diagrams of the grounding grid represented in Figures 13 and 14, show that initially, at low frequencies, the system is considered to be capacitive, and the ground impedance is high due to this capacitive property. The phase angle diagram of ground frequency response is also distributed between −90° to 0°. As the frequency increases, the grounding system decreases the capacitive property, and enhances the inductive characteristics. This reduces the impedance magnitude, which pushes the phase angle to zero. Then, the higher the frequency, the more indicative of its inductive property, thus, the impedance magnitude of the grounding system will increase. As expected, the frequency response phase angle will reach around 90 degrees.
As explained in [18], the Grcev model uses the electromagnetic formulation to derive the semiprecise equivalent circuit of horizontal and vertical model. Since this model is approximately matched with the real data up to one MHz, we use a comparison of proposed approach to Grcev model it this range of frequency. The formulations founded in [18], despite the accurate solution, are hired, to derive the equivalent circuit of each horizontal and vertical electrodes, individually, while the proposed approach presented here uses only one equivalent circuit for all grounding grid. Therefore, our strategy can be implemented much more easily than the Grcev model.   The ground potential rise (GPR), or Earth potential rise, is a phenomenon that occurs when large amounts of electricity enter the earth. This is typically caused when substations or high-voltage towers fault, or when lightning strikes occur. When currents of large magnitude enter the earth from a grounding system, not only will the grounding system rise in electrical potential, but so will the surrounding soil as well. The voltages produced by a GPR or Earth potential rise event can be hazardous to both personnel and equipment. As described earlier, the soil has resistance known as soil resistivity, which will allow an electrical potential gradient or voltage drop to occur along the path of the fault current in the soil.
The resulting potential differences will cause currents to flow into any and all nearby grounded conductive bodies, including concrete, pipes, copper wires and people [25][26][27][28][29]. The GPR calculated in different soil resistivity is visible in Figure 15. As can be seen, the CIGRE function as an excitation As explained in [18], the Grcev model uses the electromagnetic formulation to derive the semi-precise equivalent circuit of horizontal and vertical model. Since this model is approximately matched with the real data up to one MHz, we use a comparison of proposed approach to Grcev model it this range of frequency. The formulations founded in [18], despite the accurate solution, are hired, to derive the equivalent circuit of each horizontal and vertical electrodes, individually, while the proposed approach presented here uses only one equivalent circuit for all grounding grid. Therefore, our strategy can be implemented much more easily than the Grcev model.
The ground potential rise (GPR), or Earth potential rise, is a phenomenon that occurs when large amounts of electricity enter the earth. This is typically caused when substations or high-voltage towers fault, or when lightning strikes occur. When currents of large magnitude enter the earth from a grounding system, not only will the grounding system rise in electrical potential, but so will the surrounding soil as well. The voltages produced by a GPR or Earth potential rise event can be hazardous to both personnel and equipment. As described earlier, the soil has resistance known as soil resistivity, which will allow an electrical potential gradient or voltage drop to occur along the path of the fault current in the soil.
The resulting potential differences will cause currents to flow into any and all nearby grounded conductive bodies, including concrete, pipes, copper wires and people [25][26][27][28][29]. The GPR calculated in different soil resistivity is visible in Figure 15. As can be seen, the CIGRE function as an excitation source with a 400 ohm lightning channel is multiplied to the proposed impedance diagram, and GPR is calculated compared to the ideal node within zero potential relatives. The higher the soil resistance, the higher the GPR production will be.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 18 source with a 400 ohm lightning channel is multiplied to the proposed impedance diagram, and GPR is calculated compared to the ideal node within zero potential relatives. The higher the soil resistance, the higher the GPR production will be. Figure 15. GPR of a simple grounding grid.

A Complicated Grounding Grid
The grounding system under consideration in this section, as shown in Figure 16, includes several circular electrodes, a number of vertical rods, and numerous horizontal electrodes that are interconnected. Their apparent characteristics, in terms of length and positioning, have been adjusted so that the grounding system can be used as a wind turbine system or a high voltage electric tower grid. Thus, the effect of soil ionization dependence on frequency is also evaluated, and the proposed modeling effect is compared with the frequency-independent conditions.

A Complicated Grounding Grid
The grounding system under consideration in this section, as shown in Figure 16, includes several circular electrodes, a number of vertical rods, and numerous horizontal electrodes that are interconnected. Their apparent characteristics, in terms of length and positioning, have been adjusted so that the grounding system can be used as a wind turbine system or a high voltage electric tower grid. Thus, the effect of soil ionization dependence on frequency is also evaluated, and the proposed modeling effect is compared with the frequency-independent conditions. Figure 15. GPR of a simple grounding grid.

A Complicated Grounding Grid
The grounding system under consideration in this section, as shown in Figure 16, includes several circular electrodes, a number of vertical rods, and numerous horizontal electrodes that are interconnected. Their apparent characteristics, in terms of length and positioning, have been adjusted so that the grounding system can be used as a wind turbine system or a high voltage electric tower grid. Thus, the effect of soil ionization dependence on frequency is also evaluated, and the proposed modeling effect is compared with the frequency-independent conditions.  Figure 17 show the impedance magnitude and its phase angle in several soil resistivity with the effect of frequency-dependent soil ionization. As is observed, the proposed method is well suited to the data extracted from the CDEGS software, which suggests its high accuracy and effectiveness. Conductivity and permittivity are two important electrical properties of soil. Several factors, such as soil compaction, temperature, moisture, and grain size, affect soil conductivity and permittivity. Soil ionization occurs because of the electric field enhancement in soil when subjected to high impulse current discharge, which in turn causes the electric breakdown of air voids enclosed between soil particles. The air breakdown causes current to flow in the voids, and this is usually in the form of an arc. The growth of the arc within the whole soil structure is further described as containing many discrete breakdown channels and branching in nature. Soil medium is a mixture of fragmented  Figure 17 show the impedance magnitude and its phase angle in several soil resistivity with the effect of frequency-dependent soil ionization. As is observed, the proposed method is well suited to the data extracted from the CDEGS software, which suggests its high accuracy and effectiveness. Conductivity and permittivity are two important electrical properties of soil. Several factors, such as soil compaction, temperature, moisture, and grain size, affect soil conductivity and permittivity. Soil ionization occurs because of the electric field enhancement in soil when subjected to high impulse current discharge, which in turn causes the electric breakdown of air voids enclosed between soil particles. The air breakdown causes current to flow in the voids, and this is usually in the form of an arc. The growth of the arc within the whole soil structure is further described as containing many discrete breakdown channels and branching in nature. Soil medium is a mixture of fragmented mineral materials, organic matter, and gases. The size and distribution of both the soil particles and the air-filled pores are not uniform. To facilitate the computations for engineering purposes, the soil can be represented by an ordered network of soil particles and air voids, as shown in Figure 4. The impulse current mostly flows through the soil particles when there is no ionization, because of the much lower soil resistance. Therefore, the soil ionization phenomena are very important to be considered. If soil ionization is assumed to be frequency independent, Figure 18 are obtained, and the outcomes of the proposed algorithm will still match the data extracted from the CDEGS finite element software. It can be concluded that the impedance of the typical grounding grids buried in soil changes as a function of frequency. The characteristic of the grounding electrode is assessed depending on the soil resistivity. The characteristics of the grounding electrode in low resistivity and high-resistivity soils are dominantly capacitive and inductive, respectively. In both cases of frequency-independent and frequency-dependent, the magnitude and phase angle in all figures varied correspondingly as the frequency increases. In the case of frequency-dependent, the increase of magnitude is a bit higher than the frequency-independent case. On the contrary, the phase angle, in the case of frequency-dependent, increases, but not in the same way as the angle in the case of frequency-independent. There are some spikes in frequency dependent diagrams which rely on the evacuation of the magnetic field at some susceptible points in the soil. frequency-independent and frequency-dependent, the magnitude and phase angle in all figures varied correspondingly as the frequency increases. In the case of frequency-dependent, the increase of magnitude is a bit higher than the frequency-independent case. On the contrary, the phase angle, in the case of frequency-dependent, increases, but not in the same way as the angle in the case of frequency-independent. There are some spikes in frequency dependent diagrams which rely on the evacuation of the magnetic field at some susceptible points in the soil. varied correspondingly as the frequency increases. In the case of frequency-dependent, the increase of magnitude is a bit higher than the frequency-independent case. On the contrary, the phase angle, in the case of frequency-dependent, increases, but not in the same way as the angle in the case of frequency-independent. There are some spikes in frequency dependent diagrams which rely on the evacuation of the magnetic field at some susceptible points in the soil.

Conclusions
When the high-frequency lightning strokes collide to the ground as a shock wave and excite it, the Earth high-frequency components emerge from the ground, causing soil ionization. Afterwards, the soil particles properties cause the non-linearization of the transient state of buried grounding grid under the soil. If the grounding system is buried in multi-layered soils with different special resistivity through different geometries, such as cylindrical or spherical soils components, the modeling of its electrodes to analyze power networks in software packages is practically impossible. This paper presents a novel method to evaluate all these phenomena, in order to analyze the transient

Conclusions
When the high-frequency lightning strokes collide to the ground as a shock wave and excite it, the Earth high-frequency components emerge from the ground, causing soil ionization. Afterwards, the soil particles properties cause the non-linearization of the transient state of buried grounding grid under the soil. If the grounding system is buried in multi-layered soils with different special resistivity through different geometries, such as cylindrical or spherical soils components, the modeling of its electrodes to analyze power networks in software packages is practically impossible. This paper presents a novel method to evaluate all these phenomena, in order to analyze the transient state of electromagnetic waves due to the lightning strikes. Considering this technique, each grounding grid can be modeled and appraised using a simple multi-circuit model, with low-order characteristic approximation equations in a wide frequency range. Another feature of the proposed strategy is the high accuracy compared with some existing methods. Regarding those approaches, there is no possibility to model and predict the wave propagation delay. In addition, the calculation speed is more accurate than that of the electromagnetic waves theory. On the other hand, one can analyze the most complex grounding systems by any conventional processor. However, in many previous methods, parallel processors are required with very high power, due to the high-usage of current and voltage dependent sources. Finally, the results and output diagrams obtained from the proposed technique are validated with the Grcev formulation. Table 1 compares all existing approaches.