**1. Introduction**

A grounding system, which is essential for proper, reliable and safe operation of a power system, can be accomplished by providing a true reference to the electrical system that controls the discharge path of high energy faults. The performance of the grounding system during power frequency faults is a present key driver in general design of substation grounding system [1]. The transient behavior of the grounding system during a lightning discharge is characterized by a steep front that induces inductive effects in the grounding system. Thus, a large short-term voltage rises within the grounding system region, close to the injection point. The uneven voltage distribution can be explained in terms of current and voltage waves traveling along the grounding grid conductors that can be modeled by the telegrapher's equation [2]. Problematic lightning surge behavior has long been recognized by the industry and several models to describe these transient events. Related topics have been proposed in the relevant research literatures [3–7]. The reviewed literatures tend to treat high energy faults originating from direct lightning strikes discharged trough the grounding system. A method of implementation to evaluate the effect associated with lightning transient in the transmission system and the corresponding effect on the grounding system is not publicly available.

Lightning surges on overhead transmission lines may introduce travelling waves which have the potential to penetrate deep into a substation and present hazards ElectroMagnetic Interference (EMI) for sensitive equipment. To ensure reliable operation in the design of transmission systems, the grounding system should rather be considered as an integrated part. This paper presents an approach whereby the combined simulation of the grounding and transmission system is performed based on software implementation. Since the substation surge arrester is evaluated as the injection point to the grounding system, the surge arrester performance is considered. Simultaneously, the time and spatial distribution of the potential rise in the grounding itself is presented to consider EMI preventative action in the substation area.

The grounding system is implemented in Matlab/Simulink in combination with the specialized commercial software for power system transients studies, EMTP-RV. Detailed software specification is found in Appendix A.

#### **2. Parameters of the Grounding System**

Based on the classic work by Sunde, the grounding system in this work has been modelled as a lossy transmission line [2]. The soil medium with electrical resistivity and permittivity surrounds the grounding wires that are characterized by their electrical parameters thus forming a unified system. The electrical parameters in per-unit length are defined through Equation (1):

$$\mathbf{Z} = \mathbf{j}\omega \mathbf{L} \tag{1a}$$

$$
\mathcal{Y} = \mathcal{G} + j\omega \mathcal{C} \tag{1b}
$$

where the per-unit length *G* is defined as the grounding system conductance (S), *C* is the capacitance (F) and *L* is the inductance (H).

With relatively short conductor length and large cross section of grounding wires, the internal resistance (including the skin effect) and inductance are significantly smaller than the external self-inductance. With this consideration, a simplification that only includes the self-inductance is made [8]. The connection between the soil properties and the grounding system elements defined above is represented in Equations (2)–(4) [2] for a horizontal buried grounding wire.

$$G\_i = G\_j = \frac{\pi}{\rho\_{soil} \left[ \ln \left( \frac{2I}{\sqrt{2ad}} - 1 \right) \right]} \tag{2}$$

$$\mathbf{C}\_{i} = \mathbf{C}\_{j} = G\_{i} \rho\_{soil} \left( \epsilon\_{0} \times \epsilon r\_{soil} \right) \tag{3}$$

$$L\_i = L\_j = \frac{\mu\_0}{2\pi} \left( \ln \frac{2l}{a} - 1 \right) \tag{4}$$

where *ρsoil* is the soil resistivity ( Ωm),  *rsoil* is the soil relative permittivity (-), *a* is the conductor radius (m) and *d* is the buried depth (m).

A grounding wire conducting a lightning impulse current will exert a time-varying electric field, outwards through the wire and into the surrounding soil. The soil itself, depending on soil resistivity and properties, will conduct a current from the grounding wire, dissipated into the soil. Depending on the dissipated current density , the electric field and the current density are given in Equation (5a). The surface current density of a round wire is found in Equation (5b), which exerts the electric field on the soil [9] (p. 1586). The linear behaviour between current density and electric field are valid below the soil critical breakdown value, *Ec*.

$$E\_{\rm soil} = J\_{\rm soil} \rho\_{\rm soil} \tag{5a}$$

$$J\_{\rm soil} = \frac{I\_{\rm soil}}{2\pi al} \tag{5b}$$

#### **3. Modeling the Grounding System**

The element length in per-unit model of the grounding wire is implemented in Matlab as a T-section of the transmission line model. The grounding wire implementation is illustrated with fundamental electrical elements and measurement nodes in Figure 1 which forms a horizontal buried grounding wire. From the grounding wire t-section development, the grounding system is implemented by an appropriate number of T-sections and additional nodes using Simulink graphical block elements. The grounding system layout with formation strategy is illustrated Figure 2 which gives references to the simulation database is found in Appendix B. To assure a smooth representation of the interconnection between elements, and to possibly extend the model functionality to include mutual couplings, the T-section represent one meter of ground wire (in per-unit) [4]. Consequently, the model can be extended to arbitrary lengths and configurations.

**Figure 1.** Element lenght in per-unit grounding wire illustration of implementation.

**Figure 2.** The grounding system layout with formation strategy.

This consideration gives the total number of required logged variables for the formed horizontal buried grounding grid square meshes as an indication for two different mesh sizes and total area in Table 1.


**Table 1.** Required number of logged variables for grounding grid of two different square mesh sizes and total area values.

As described trough Section 2 the grounding system is modelled as a transmission line, connecting the soil and grounding conductor properties. The matrices of the grounding system may be expressed trough the coupled telegraphers equations in frequency domain trough Equation (6). Where *<sup>V</sup>*(*z*) and *<sup>I</sup>*(*z*) are the line phasor of voltage and current [10]:

$$\frac{d^2}{dz^2}\hat{\mathcal{V}}(z) = \mathcal{Z}\hat{\mathcal{V}}(z)\tag{6a}$$

$$\frac{d^2}{dz^2}\mathcal{I}(z) = \hat{Y}\mathcal{Z}\mathcal{I}(z)\tag{6b}$$

In transient analysis, considering the spatial distribution of voltage potential and current flowing in the grounding system, it is mandatory to simulate each element in time domain. With the layout of grounding grids, as is illustrated Equation (2), several connection points exists. Consequently, the time domain model is required to account for the coupling between elements as feeding points. For *J* feeding points the impact of the grounding system, as an electrical network exited by an injected lightning current, is global. Each element in the grounding system matrices is numerically solved independently by the Matlab ODE23t solver, which implements the trapezoidal rule using a "free" interpolation for the time domain solution. The implementation account for the injected current characteristics and corresponding frequency response of the grounding system.

#### **4. Grounding Model Verification**

The grounding system model is validated by the work based on the ElectroMagnetic Field (EMF) theory first performed by Grcev [11] and later by Jardines et al. [7] who introduced a variant of the Multi-conductor Transmission Line (MTL) approach. In these cases of model verifications, the current source is implemented using the double exponential waveform, *<sup>i</sup>*(*t*) = ˆ*<sup>I</sup>*(*e*<sup>−</sup>*α<sup>t</sup>* − *<sup>e</sup>*<sup>−</sup>*β<sup>t</sup>*), and the source parameters where adjusted to fit the given stroke function. The referenced work is reproduced from manual reading.

#### *4.1. Grounding Rod of 15 m Length*

A grounding rod of 15 m was simulated and measured by Grcev [11] and later evaluated by Jardines et al. [7] as shown in Figure 3a. Using two independent modeling approaches and relative similarities in results may justify the present modeling accuracy. The grounding rod is horizontally buried in soil of *ρsoil* = 70 Ωm and  *rsoil* = 15, with a wire radius of *a* = 0.012 m at depth *d* = 0.6 m. The current source was set with the amplitude of ˆ *I* = 36 A and the stroke of 0.36/12 μs which leads to the new parameter values *α* = 32 × 10<sup>3</sup> and *β* = 7.6 × 106. The results from the implemented model are shown in Figure 3b. As it can be observed from Figure 3a in the referenced work, the peak values at the injection point are approximately 560 V. While using the implemented model, Figure 3b, the peak value shows 569 V. The second point for comparison is at 7 m from the injection point where Grcev simulations show approximate 200 V while the results in this work gives 172 V after 0.7 μs. Besides the point value readings, the surge wave propagation along the grounding rod and voltage distribution are of similar character.

**Figure 3.** Voltage distribution along an horizontal copper wire in ground (length *l* = 15 m, conductor radius *a* = 0.012 m, at soil depth *d* = 0.6 m) in soil of *ρsoil* = 70 Ωm and  *rsoil* = 15: (**a**) the reproduced results for comparison from the research given in [11] (p. 818); and (**b**) results from implemented model with excitation current ˆ *I* = 36 A of double exponential waveform (*α* = 32 × 103, *β* = 7.6 × 106).

#### *4.2. Grounding Grid of 10 m Meshes with Total Size 3600 m*<sup>2</sup>

A grounding system of 10 × 10 m meshes size and a total square area of 3600 m<sup>2</sup> was simulated by Jardines et al. [7] using a variant of the MTL approach (see Figure 4a). The grounding grid was buried in soil of *ρsoil* = 100 Ωm and  *rsoil* = 36, with a wire radius of American Wire Gauge (AWG) 2/0 (*a* ≈ 0.004126 m) at depth *d* = 0.6 m. The current source was set with amplitude of ˆ*I* = 1 kA and a 1/20 μs which gave adjusted parameters to *α* = 38 × 10<sup>3</sup> and *β* = 2.54 × 106. The results from the implemented model are given in Figure 4b. Since the injected current was not given in the referenced work (see Figure 4a) it is worth noting that a small deviation in the injected current will have a large impact on the grounding system response, especially for the region contributing to limiting the peak voltage. As it can be observed, both the voltage distribution and the propagation characteristics correspond well. When comparing the simulation to Grcev [11] which was based on EMF against the simulation in the presented work, the peak value and the propagation characteristics are of comparable values. However, the distributed voltage at the given nodal points is more conservative in the transmission line approach, both in the implemented model and in Jardines et al. work [7]. This unveils a model accuracy difference compared to the EMF. As reviewed, Liu [4] developed the non-uniform transmission line approach to compensate for the inaccuracy of this method and are treated in details in her work.

**Figure 4.** Voltage distribution in a grounding grid consisting of 6 × 6 meshes of 10 m size. The grounding grid consist of copper conductors of AWG 2/0, buried at *d* = 0.6 m in soil of *ρsoil* = 100 Ωm and  *rsoil* = 36: (**a**) the reproduced results for comparison from research given in [7] (p. 31); and (**b**) results from implemented model with excitation current ˆ *I* = 1 kA of double exponential waveform (*α* = 38 × 103, *β* = 2.54 × 106).

#### **5. Integration of Grounding and Transmission System**

The Matlab/Simulink grounding grid are integrated with EMTP-RV through a newly developed FMI software, which was released by Powersys Solutions in early 2018 [12,13]. The FMI package gives possibilities for co-simulation with information exchange at a per simulation time-step interval (sequentially processed). The specific FMI interface used in this study is found in Appendix A.2. From the transmission system surge arrester, the injection point impedance describes the grounding system response trough Equation (7).

$$Z\_{InjectPoint}(j\omega) = \frac{\mathcal{U}\_{InjectPoint}(j\omega)}{I\_{arrester}(j\omega)}\tag{7}$$

When the transmission system surge arrester reaches the breakdown voltage a current surge is injected into the grounding system. The current injection value is exchanged from EMTP-RV to Matlab, which simulates the grounding system response. The grounding system initial impedance is set to the power frequency value corresponding to the grounding system resistance [14]. From the injection point current and the induced voltage potential rise, the impulse impedance is calculated and exchanged from Matlab to EMTP-RV, which connects the dynamic response of the grounding system to the transmission system. With a large number of logged variables required by the grounding system (see Table 1) and presented implementation strategy, the additional measured values of the transmission system were exchanged from EMTP-RV to Matlab to provide a common simulation log for pre-processing. The advanced functionality offered by the Matlab/Simulink modeling of the

grounding system lies in the preprocessing of large data-sets. A schematic overview of the model integration is illustrated in Figure 5.

**Figure 5.** Schematic overview of the integration between grounding and transmission system.

#### **6. An Example of Results: Case Study**

A simplified transmission system network is shown in Figure 6 with the grounding system model interfaced in Matlab. In this model the transmission line, the cable, the transformer and the surge arrester are selected from the standard EMTP-RV software library. At 10 km distance from a substation, a lightning strikes the 300 kV overhead transmission line (*Zcl* = 400 Ω). A shielding failure causes an injected current with a magnitude of 1 kA and 1.2/50 μs of CIGRE waveform stroke to flow towards the substation. In the substation, the cable (*Zcc*) between the surge arrester (Appendix C) and transformer are 10 m. Figure 7 shows the simulation results of the transmission system when the grounding system is ignored, thus giving a peak voltage of 171 kV at the transformer.

**Figure 6.** Simplified illustration of implemented transmission and grounding system.

A grounding system with 10 × 10 m mesh size is added with a total grid area of 3600 m<sup>2</sup> in soil of *ρsoil* = 2000 Ωm,  *rsoil* = 16, with a wire radius of *a* = 0.04126 m at depth *d* = 0.6 m is obtained. The transmission system conditions are similar. The surge arrester is connected in center of the grid. The simulation results of the transmission system are given in Figure 8. As it can be observed, the surge arrester performance is reduced to give a peak transformer voltage of 190 kV.

**Figure 7.** Ignoring grounding system: transmission system nodal voltages and CIGRE 1.2/50 μs injected current stroke in far-end (*ll* = 10 km).

**Figure 8.** Application case when the grounding system is added: (**a**) the effect in the transmission system; and (**b**) the results at the surge arrester injection point.

Then, the EMI analysis of the substation could be performed based on the overall voltage distributions in the grounding system as shown in Figure 9. Figure 9a shows measurments of the grounding system, with nodal points selected diagonaly outwards from the center injection point. For the selected events an overall voltage distribution is presented. Figure 9b shows the distribution at center peak voltage and Figure 9c when the voltage potential in the grid corners are at peak.

With the comprehensive log dataset further analysis is exemplified in Figures 10 and 11. The impulse effective area of the grounding grid defined in [15] is shown in Figure 10. This is the total area of the grounding grid which limits the peak voltage in the grid. There exists several definitions and empirical formulas for estimating the effective length when optimizing the grounding grid design that was recently evaluated [5]. However, these approaches have not taken the transmission line network itself as an integrated element into consideration. In addition, the electric field exerted on the soil close to the grounding wires that was based on the current density leaked to the soil is shown in Figure 11.

Institute of Electrical and Electronics Engineers (IEEE) indicates (quoted in standard) a critical breakdown value of *Ec* = 1000 kV/m [16] (p. 1263). This definition gives reference to the above mentioned experimental relation between soil resistivity and *Ec* [17]. Further evaluation of IEEE standards gives a value of *Ec* = 400 kV/m, a level which are referred without reference in [18] (p. 38) (refereed standard is currently under review for update). If the ionization level is reached the electric field in the soil has pronounced influence on the impulse peak voltage. Moreover, ionization has a positive effect by lowering the peak voltage due to arcing or puncturing in the soil. Soil ionization can be included in dynamic simulations as it has been proposed in [9]. This would modify the apparent grounding conductor radius in Equations (2) and (3). For this application case when considering strokes in the transmission system, soil ionization phenomena has been evaluated to have minor deviation of the results due to relative small currents injected trough the surge arrester.

**Figure 9.** Application case when the grounding system is added: (**a**) nodal voltage measurement values for selected points in the grounding grid; (**b**) the overall voltage distribution in the grounding grid at peak; and (**c**) the overall voltage distribution in the grounding grid at corner peak.

**Figure 10.** Application case: used area.

**Figure 11.** Application case: peak field.
