This article has been accepted for publication in a future issue of this journal, but has not been

fully edited. Content may change prior to final publication. Citation information: DOI
10.1109/TIA.2015.2409801, IEEE Transactions on Industry Applications

Inclusion of Field Solver-Based Tower Footing

Grounding Models in Electromagnetic Transients
Programs

Andrzej Holdyk Bjorn Gustavsen

SINTEF Energy Research SINTEF Energy Research
Norway Norway
andrzej.holdyk@sintef.no bjorn.gustavsen@sintef.no

Abstract— This paper presents a new approach to wide-band EMTPs have no dedicated models capable of representing
modeling of tower footing grounding impedances. Appropriate grounding structures over a wide band of frequencies. It is
representation of transmission tower footing impedances is common practice to use a fixed resistor representation
essential in electromagnetic transient simulations for proper calculated at power frequency, or using a lumped RLC model
evaluation of overvoltages and flashovers in the system during
with typical values [4]. Much more accurate results can be
fast-front transient conditions such as lightning strikes. In this
paper, a field solver is used to precisely model the most popular obtained using methods based on electromagnetic field
types of tower footing grounding systems, i.e. rods, counterpoises calculations but such methods are difficult to integrate into
and ring electrodes, for a variety of dimensions and soil EMTPs. Recently, a number of studies have appeared
parameters. Results are stored in a lookup table and used for presenting methods for integration of electromagnetic field
generating frequency-dependent models that are exported to calculation models of complex grounding systems into popular
electromagnetic (EM) transient simulators. The table and model EMTPs [5-7]. These works characterize the grounding
generation tool can be used by utility companies to model tower systems in terms of n-terminal impedance or admittance
footings while taking into account local variations along the matrices in the frequency domain, which are fitted with
right-of-way, such as local variations in soil resistivity. The
rational functions for inclusion in EMTPs. It follows that the
presented examples show close agreement of results of lightning
studies using developed electrode models included in EM tools: field solver must be reapplied every time a new configuration
PSCAD, EMTP-RV and ATP-EMTP. A sensitivity study shows is to be studied.
that the maximum voltage at the tower top during a direct This paper presents a different approach compared to [5-7]
lightning strike can be up to 19% higher when using the adopted as the field solver is no longer needed by the end user of the
modeling procedure, compared to a standard approach of using a model. The main idea is to store the results obtained by the
fixed resistor equal to the power frequency resistance. field solver in a multi-dimensional lookup array for later use.
Our focus is on precise modeling of the most popular types of
Keywords—tower footing impedance; wide-band model; tower footing grounding systems (rods, counterpoises and ring
lightning studies; electromagnetic transients
electrodes). The results (single-terminal admittance) are stored
I. INTRODUCTION in a multi-dimensional table which covers a large range of
grounding system parameters, and the information in the table
The tower footing grounding system is of major importance
is later used for generating frequency-dependent models that
for the lightning performance of overhead lines [1]. During a
are exported to EMTP simulators. The approach permits utility
direct lightning strike to the tower or the overhead ground
companies to accurately represent transmission tower footing
wires, the amplitude of the overvoltage at the tower top is
impedances in EMTP simulations without having access to the
strongly influenced by the impedance of the tower footing
field solver. The utility company simply chooses the class of
grounding system due to wave reflections from the bottom of
tower grounding systems and specifies its parameters and the
the tower [2, 3], thereby being a decisive factor for the
program automatically generates an EMTP-compatible model
probability of a back flashover to the phase conductors. The
based on the information. That way, the utility company can
lightning performance of overhead lines is usually assessed
easily model the tower footings while taking into account local
using electromagnetic transients simulation software
variations along the right-of-way, e.g. local variations in soil
(EMTPs), e.g. PSCAD, EMTP-RV, and ATP-EMTP.
resistivity.
The impedance of a grounding system is strongly
The paper describes the generation of the look-up table in
frequency-dependent and depends also on the characteristics
terms of electrode types and parameters, and the procedure for
of the ground electrode and the soil. So far, the available
model generation for various EMTP type tools. We also show
The work was supported by the Norwegian Research Council (RENERGI the functionalities of a graphical user interface (GUI) for
Programme) with additional support from DONG Energy, EdF, EirGrid, streamlining the model generation step for use with PSCAD,
Hafslund Nett, National Grid, Nexans Norway, RTE and Siemens

0093-9994 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See
http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI
10.1109/TIA.2015.2409801, IEEE Transactions on Industry Applications

EMTP-RV and ATP-EMTP. Finally, the methodology is responsible for interpolation of the admittance between
applied to the modeling of an overhead line which is struck by parameter values, fitting the admittance with rational
lightning. The obtained results are compared with those functions, enforcing passivity of the model, and exporting the
obtained using a fixed-resistor representation of the tower model parameters into formats that are compatible with the
footing grounding. most popular EMTPs, i.e. PSCAD, EMTP-RV and ATP-
EMTP.
II. PROGRAM STRUCTURE
The program consists of three main parts: a graphical user III. ELECTRODE MODELING
interface (GUI), a look-up table containing pre-computed A. Electrodes and Parameter Variations
driving-point admittances of tower footing electrodes, and a Each electrode has a single point of connection (current
set of additional routines for selecting geometry and extracting
injection point) and is described by its specific dimensions and
a model. Fig. 1 shows the general structure of the program and
soil parameters, i.e. relative permittivity 𝜖𝑟 and resistivity 𝜌,
the type of data exchanged between its parts. All program
as indicated at Fig. 2 and Fig. 3.
parts were developed in Matlab.

Fig. 3. Grounding electrodes

A range of values of the parameters was specified according to

common practice used by local utilities [8, 9] with the
boundary values shown in Table 1. The soil parameters cover
Fig. 1. General structure of the program and communication between its parts. a wide range of soil conditions, i.e. from rich, agricultural,
highly moist ground to dry and rocky ground. All conductors
The look-up table is conveniently represented by a structure are copper. The cross-section area for counterpoise and ring
as shown in Fig. 2. The structure format allows for easy electrodes is 50 𝑚𝑚2 and for rod electrodes 120 𝑚𝑚2 .
alteration in the future, e.g. by adding more electrode types or
parameters. Table 1. Range of electrode parameters
Electrode type Parameters Range
Soil resistivity [Ω𝑚] 100 – 50000
All electrodes
Relative permittivity 1 – 100
Rod Length [𝑚] 1 – 20
Length of radial [𝑚] 1 – 80
Counterpoise:
Depth [𝑚] 0.5 – 1
Radius [𝑚] 1 – 10
Ring
Depth [𝑚] 0.5 – 1

B. Field Solver
CDEGS is a software package specifically designed for
analysis and design of power system grounding using field
calculations. It calculates conductor currents and
electromagnetic fields generated by an arbitrary network of
energized conductors, both above and below ground [10]. A
HIFREQ module is used for wide band modeling of electrodes
as the module is designed for precise electromagnetic field
analysis. The module also allows user-defined specification of
soil resistivity, permittivity and permeability.
C. Calculation of Frequency Dependent Impedance
The driving point impedances are calculated in CDEGS by
Fig. 2. Storage of look-up table in a Matlab structure injecting a unit current at discrete values of frequency between
1 Hz and 10 MHz through a 0.1 m long copper conductor into
The GUI is responsible for communication with the user, an electrode (current injection point at Fig. 3) and measuring
displaying and exporting the results. Additional routines are the ground potential rise (GPR) of the conductor. The GPR is

0093-9994 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See
http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI
10.1109/TIA.2015.2409801, IEEE Transactions on Industry Applications

equal to the driving-point impedance of the electrode due to to file as branch cards for ATP-EMTP [18]. The circuit is
the unit current used. The electrode is placed in a homogenous imported using the $Include statement.
soil, i.e. there is a single soil layer. Most of the simulations are
performed automatically by employing scripting and batch V. EXAMPLE
processing capabilities of the software. High accuracy of the A. System
simulations, especially at high frequencies, is ensured by
We consider overvoltages caused by lightning strikes to the
conductor segmentation which assured accurate computation
system in Fig. 4.
of current distribution within conductors. Any mutual
The towers are modeled as lossless transmission line
electromagnetic coupling between adjacent electrodes is
segments with constant parameters: surge impedance 𝑍 =
neglected.
150 Ω and velocity of 80 % of speed of light in vacuum. The
IV. INTERFACING WITH EMTP earth wires are represented using the same model and
parameters, but with 𝑍 = 280 Ω [3]. The phase conductors are
The interfacing with EMTPs is done in four steps: 1) The
not included as the study investigates the maximum
driving point impedance calculated by CDEGS is extracted
overvoltages in the tower top. Only three towers are included
from the look-up table and converted into admittance 𝑌(𝜔) =
in the model: the tower struck by lightning and the two
𝑍(𝜔)−1 ; 2) 𝑌(𝜔) is subjected to fitting with rational functions
adjacent ones. The length of earth wires outgoing from the
using Vector Fitting; 3) Resulting model is subjected to
adjacent towers is increased to 40 km to prevent reflections
passivity enforcement; 4) Final model is saved into proper
from the far ends within the investigated time frame. The
format for inclusion into alternative EMTPs. Steps 2-4 are
lightning stroke is modeled using CIGRE concave lightning
described in more detail below.
current source [1] with parameters: peak value 𝐼𝑚𝑚𝑚 = 31 𝑘𝑘,
A. Vector Fitting maximum steepness 𝑆𝑚 = 35 𝑘𝑘/𝜇𝜇, front time 𝑡𝑓 = 1.1 𝜇𝜇,
Vector Fitting [11-13] is used for approximating a selected and time to half-value 𝑡ℎ = 60 𝜇𝜇. Each tower is grounded
driving-point admittance Y (ω ) with a pole-residue model, assuming a counterpoise electrode with identical parameters:
soil relative permittivity µ r = 1 , soil resistivity ρ = 700 Ωm ,
N rj
y (ω ) ≅ d + ∑ (1) counterpoise length of 25 𝑚 and depth of 1 𝑚 (see Fig. 3).
j =1 jω − a j

The model is physical in the sense that all poles {a j } are

stable and all poles and residues {rj } are real or come in
complex conjugate pairs. Also, 𝑑 is real.
B. Passivity enforcement
In order to ensure a stable time domain simulation, passivity
of the model must also be enforced. This is achieved by
solving a constrained linear least squares problem where the Fig. 4. Lightning stroke to a tower. The modeled system includes one ground
model's residues {rj } are perturbed such that the model wire and three towers with associated grounding electrodes.
satisfies the passivity criterion,
B. Voltage Waveforms Calculated by Alternative
Re{ y (ω )} > 0 ∀ω (2) Electromagnetic Transients Simulation Programs
while minimizing the change to the model's behavior [14]. The Fig. 5 shows the simulated voltage waveform at the top of
passivity enforcement is applied iteratively together with an the tower struck by lightning and at the top of the adjacent
algebraic test [15] for passivity assessment, until all passivity towers, as obtained by the three alternative EMT simulation
violations have been removed. programs. The results agree very closely as expected.
However, a small difference exists between ATP-EMTP and
C. Model Extraction to EMTPs the two other programs, being smaller than 3% and thus not
The model is exported to the desired EMTP programs by visible in the figure above. This deviation is caused by the
saving it in a text file with appropriate format. For PSCAD, a finite precision of the file format used for representing the
Norton equivalent is used, which is based on recursive pole-residue model by an equivalent electrical network [19].
convolution [16, 17] and realized in a user-defined module. The investigated system is symmetric which causes the
However, starting from version 4.6 it will be possible to use a voltages at both adjacent towers to be identical; therefore only
built-in Frequency Dependent Network Equivalent (FDNE) results from one tower are shown.
block instead. For EMTP-RV, the built-in state space block is
used and the model is imported directly by the equivalent
state-space form of (1). Internally, recursive convolution is
used also in EMTP-RV. For ATP, the pole-residue model (1)
is converted into an equivalent electrical circuit which is saved

0093-9994 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See
http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI
10.1109/TIA.2015.2409801, IEEE Transactions on Industry Applications

Fig. 5. Voltage at the tower top during lightning strike. Comparison of results Fig. 6. Voltage at tower top during lightning for two models of ground
from three simulation programs: PSCAD, ATP-EMTP and EMTP-RV. electrodes: wide-band and lumped resistance.

C. Comparison with lumped model based on power frequency D. Sensitivity Study

resistance The use of the wide-band grounding model in EMT
Reference [4] recommends using a simplified model (RLC programs is as easy as any standard lumped component.
circuit) for representation of footing impedances in fast-front Therefore, the wide-band grounding electrode models might
transient studies, or alternatively a transmission line in case of be used for a variety of optimization or sensitivity studies. The
counterpoises. In practice, it is difficult to obtain parameters following example assesses the simultaneous influence of the
for both models and often a fixed resistance calculated at lightning current front steepness and the soil resistivity on the
power frequency is used instead. voltage peak value in the struck tower. Here, EMTP-RV is
The power frequency resistance of counterpoise buried in used to calculate a maximum value of voltage at the stroke
homogenous soil can be obtained by [20]: tower for a large range of values of the two parameters.
 π ⋅m  All grounding electrodes have the same parameters and the
n −1 1 + sin
ρ  2⋅l n  resistivity of soil is varied between 100 and 5000 Ω𝑚 in steps
=R0 ⋅  ln −1+ ∑  (3) of 100 Ω𝑚.
n ⋅π ⋅ l  2⋅a⋅d π ⋅m 
 m =1
sin  The parameters of the circuit are the same as in the previous
 n 
example except for the following:
where ρ - resistivity of soil, 𝑛-number of radial wires o 𝑆𝑚 is linearly distributed between 5 and 100 𝑘𝑘/𝜇𝜇,
(𝑛 = 4), 𝑙 - length of counterpoise wire, 𝑎 – radius of the wire, which corresponds approximately to lower 1% and upper
𝑑 – depth of counterpoise. Equation (3) assumes the same 99% of cumulative distribution function (PDF) of 𝑆𝑚 for
angle between any two adjacent wires, which corresponds impulse current crest of 31 𝑘𝑘 [1];
with our case. The resistance at 50 Hz calculated from the o Front time depends on the value of 𝑆𝑚 , i.e. t f = I F / S m ;
wide-band model is 14.76 Ω while (3) results in 14.78 Ω.
Reference [9] compares power frequency resistance calculated o Soil permittivity is 3.
with (3) and CDEGS for different length of counterpoise wire. Fig. 7 and 8 show, respectively, the peak voltage obtained
Deviations up to 23 % were shown for short counterpoise using the wide-band electrode model and the model based on
arms but decreasing with increasing counterpoise arm length. (3). The peak voltage is clearly dependent on both frequency
This result might be due to differences in the mutual coupling and soil resistivity. Fig. 9 shows the difference in percent
effect between arms when the arm length changes. Much between the two results. Positive value means the crest voltage
better agreement was found for ring and rod electrodes. obtained with wide-band electrode model is greater than the
Fig. 6 shows a comparison of voltages at the top of the corresponding result obtained using constant value of
struck tower for two models of footing impedance: 1) the resistance.
wide-band model used in this work, and 2) the lumped It is often assumed that modeling the ground impedance by
resistance (3) at power frequency, as recommended in [12]. a constant resistance at power frequency guarantees
The crest value of voltage at the tower is observed to be 11 % conservative results; however, in this example conservative
higher when the wide-band model is used. All calculations results were obtained only for large values of soil resistivity
were performed in PSCAD. and impulse front steepness. In all other cases the voltages
obtained with constant resistance are underestimated, with a
maximum difference of 19 %.

0093-9994 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See
http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI
10.1109/TIA.2015.2409801, IEEE Transactions on Industry Applications

VI. DISCUSSION
A. Soil Ionization Effects
The proposed approach does not take into account possible
soil ionization effects which in case of large currents might
play an important role in reducing the effective impedance
[21]. The use of the presented wide-band models can therefore
be considered as a conservative approach by overestimating
the grounding impedance.
B. Electrode Parameter Sampling Density
Simulations using a field solver are very time consuming;
therefore, the number of parameter values should be as low as
possible without decreasing the accuracy of the results and
thus the usability of the software. Interpolation techniques
Fig. 7. Maximum value of voltage at the top of a tower struck by lightning for might be used to obtain any value of parameters within the
different values of 1) lightning current front steepness 𝑆𝑚 and 2) soil range given in Table 1; however, this might decrease the
resistivity. Tower footing impedance represented by wide-band model. accuracy of the model or, depending on the parameters, even
cause the model to become unphysical. As the impedance of
the tower footing electrodes depends strongly on the
frequency, soil parameters and electrode's dimensions, the
proper choice of values for each of the parameter is not trivial.
As an example, Fig. 10 shows dependency of earthing rod's
resistance on length and frequency of excitation for a case of
high soil permittivity. Several resonance frequencies are now
visible for electrodes of length below 10 m. Therefore, a
denser sampling density should in this case be used for short-
length electrodes.
With proper parameter sampling density we have obtained
acceptable accuracy using linear interpolation for extracting
intermediate values. This however, has not been tested for all
combinations of parameters.
Fig. 8. Maximum value of voltage at the top of a tower struck by lightning for
different values of 1) lightning current front steepness 𝑆𝑚 and 2) soil
resistivity. Tower footing impedance represented by constant resistance
obtained by (3).

Fig. 10. Resistance of earthing rod electrode in frequency depending on the

electrode's length. Soil resistivity is 1 and permittivity is 100.

VII. SUMMARY AND CONCLUSION

Fig. 9. Difference in percentage between results shown at Fig. 8 and Fig. 7.

This paper presents a new approach to wide-band modeling of
tower footing ground impedance, especially useful for utility
companies. The three most popular types of tower footing
electrodes are modeled in a field solver for a variety of
dimensions and soil parameters. This information is collected
in a multi-dimensional database and used to create models
compatible with electromagnetic transients simulators. The

0093-9994 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See
http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI
10.1109/TIA.2015.2409801, IEEE Transactions on Industry Applications

presented examples show excellent agreement between results [9] I. Tikhomirova, "Project report in TET 5500: Modeling of grounding
structures," Department of Electrical Power Engineering, Norwegian
obtained by PSCAD, EMTP-RV and ATP-EMTP. University of Science and Technology, Trondheim 2012.
A sensitivity study was used to analyze the peak voltage at [10] SES. (2015). Web Page of CDEGS software package. Available:
a tower top struck by lightning. It was found that the wide- http://www.sestech.com/Products/SoftPackages/CDEGS.htm
band electrode model predicts a maximum voltage which [11] B. Gustavsen and A. Semlyen, "Rational approximation of frequency
domain responses by vector fitting," IEEE Transactions on Power
differs between −11 % and +19 % from that obtained by a Delivery, vol. 14, pp. 1052-1061, 1999.
constant resistance model evaluated at power frequency. The [12] B. Gustavsen, "Improving the pole relocating properties of vector
wide-band model should therefore be used as the constant fitting," Power Delivery, IEEE Transactions on, vol. 21, pp. 1587-1592,
resistance model may result in a significant underestimate of 2006.
[13] D. Deschrijver, M. Mrozowski, T. Dhaene and D. De Zutter,
the overvoltages. "Macromodeling of Multiport Systems Using a Fast Implementation of
the Vector Fitting Method," Microwave and Wireless Components
ACKNOWLEDGMENT Letters, IEEE, vol. 18, pp. 383-385, 2008.
The authors gratefully acknowledge the contribution of M.Sc [14] A. Semlyen and B. Gustavsen, "A Half-Size Singularity Test Matrix for
Fast and Reliable Passivity Assessment of Rational Models," Power
student Irina Tikhomirova, who performed field solver Delivery, IEEE Transactions on, vol. 24, pp. 345-351, 2009.
calculations. We also extend our gratitude to professor Hans [15] B. Gustavsen, "Fast Passivity Enforcement for Pole-Residue Models by
Kristian Høidalen for his advice. Perturbation of Residue Matrix Eigenvalues," Power Delivery, IEEE
Transactions on, vol. 23, pp. 2278-2285, 2008.
REFERENCES [16] B. Gustavsen and O. Mo, "Interfacing convolution based linear models
to an electromagnetic transients program," in Proc. Int. Conf. Power
[1] "CIGRÉ WG 33-01 (convener A. Eriksson): Guide to Procedures for Systems Transients (IPST), 2007, pp. 4-7.
Estimating the Lightning Performance of Transmission Lines," in [17] B. Gustavsen and H. M. J. De Silva, "Inclusion of Rational Models in an
CIGRÉ Brochure vol. 063, ed: CIGRÉ 1991. Electromagnetic Transients Program: Y-Parameters, Z-Parameters, S-
[2] "Modeling guidelines for fast front transients," Power Delivery, IEEE Parameters, Transfer Functions," IEEE Transactions on Power Delivery,
Transactions on, vol. 11, pp. 493-506, 1996. vol. 28, pp. 1164-1174, 2013.
[3] J. Martinez-Velasco, Power system transients. Parameter determination. [18] B. Gustavsen, "Computer code for rational approximation of frequency
Boca Raton, Fla. u.a.: CRC Press, 2010. dependent admittance matrices," Power Delivery, IEEE Transactions on,
[4] IEC/TR, "Part 4: Computational guide to insulation co-ordination and vol. 17, pp. 1093-1098, 2002.
modelling of electrical networks, Ed.1," in 60071-4 Insulation co- [19] A. C. S. Lima, B. Gustavsen and A. B. Fernandes, "Inaccuracies in
ordination, ed: International Electrotechnical Commission, 2004. network realization of rational models due to finite precision of RLC
[5] M. Heimbach and L. D. Grcev, "Grounding system analysis in transients branches," in International Conference on Power System Transients,
programs applying electromagnetic field approach," IEEE Transactions IPST, Lyon, France, 2007.
on Power Delivery, vol. 12, pp. 186-193, Jan 1997. [20] E. D. Sunde, Earth conduction effects in transmission systems. New
[6] P. Sarajcev, S. Vujevic and D. Lovric, "Interfacing harmonic York: D. Van Nostrad Company, Inc., 1949.
electromagnetic models of grounding systems with the EMTP-ATP [21] M. Heimbach and L. Grcev, "Simulation of grounding structures within
software package," Renewable Energy, vol. 68, pp. 163-170, Aug 2014. EMTP," in Proc. of 10th International Symposium on High Voltage
[7] K. Sheshyekani, M. Akbari, B. Tabei and R. Kazemi, "Wideband Engineering in Montreal, 1997, pp. 131-135.
Modeling of Large Grounding Systems to Interface With
Electromagnetic Transient Solvers," Power Delivery, IEEE Transactions
on, vol. PP, pp. 1-1, 2014.
[8] I. Tikhomirova, "MSc Thesis: Electromagnetic Transients Modelling of
Grounding Structures," Master of Science in Electric Power
Engineering, Department of Electric Power Engineering, Norwegian
University of Science and Technology, Trondheim, 2013.

0093-9994 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See
http://www.ieee.org/publications_standards/publications/rights/index.html for more information.