Data Driven Transfer Functions and Transmission

Network Parameters for GIC Modelling

M. J. Heyns∗† , C. T. Gaunt∗ , S. I. Lotz† and P. J. Cilliers∗†
∗
Department of Electrical Engineering
University of Cape Town (UCT), Cape Town, South Africa
† South African National Space Agency (SANSA)

Hermanus, South Africa

mheyns@sansa.org.za

Abstract—Typical geomagnetically induced current (GIC) the measured geomagnetic field (B-field). The E-field and
modelling assumes the induced quasi-DC current at a node B-field can be related in the frequency domain through the
in the transmission network is linearly related to the local surface impedance. Surface impedance models have varying
geoelectric field by a pair of network parameters. Given a
limited time-series of measured geomagnetic and GIC data, degrees of complexity and of the several methods developed,
an empirical method is presented that results in a statistically a multi-layered ground conductivity is widely used due to its
significant generalised ensemble of parameter estimates with the generality, simplicity and apparent accuracy [3]. Regardless
error in the estimates identified. The method is showcased for of E-field derivation method, assumptions of DC-equivalence
different transmission networks and geomagnetic storms and, and constant network parameters remain in GIC modelling
where prior modelling exists, shows improved GIC estimation.
Furthermore, modelled networks can be locally characterised and approaches. A recent paper [4] examined the DC-equivalence
probed without any further network knowledge. Insights include assumption by comparing measured E-field and GIC data at
network parameter variation, effective network directionality a node, showing that empirical a and b parameters as defined
and response. Merging the network parameters and geoelectric by (1) in the frequency domain are frequency dependent. The
field estimation, a transfer function is derived which offers an source of the frequency dependence is difficult to pin down
alternative approach to assessing transformer exposure to GICs.
since the measured E-field at a single node is not necessarily
the same as the network effective E-field. This current work
Index Terms—Ensemble estimation, geomagnetically induced
follows the nodal formulation of (1), but differs significantly
currents (GICs), transmission network parameters
in that the network parameters are not assumed to be constant.
Relaxing this assumption allows for a simple approach of a
I. I NTRODUCTION fast, effective transform from B-field to GIC with accurate
effective network parameter estimation, while acknowledging
The effects of geomagnetically induced currents (GICs)
possible unmodelled frequency dependence and other uncer-
in communications and power systems were well known for
tainty, directly applicable to planning and operations.
many years before the first significant papers on calculating
Building from the nodal formalism of GIC modelling, the
the GICs appeared [1]. Initially, the driving near-Earth current
ensemble methodology of network parameter estimation is
system was modelled as a line or sheet, giving by first prin-
presented in Section II. Three different datasets from around
ciples different answers for the power line currents. Another
the world are used to generate network parameter ensembles
approach gave the transformer neutral currents directly - by
and test their performance in Section IV, with the data used
calculating the DC-equivalent of the voltages induced in the
described in Section III. Section IV-A looks specifically into
whole network by a uniform plane-wave [2]. For a given node,
the characteristics of GICs in the local networks derived
the traditional nodal modelling formulation is
from the network parameter ensembles. Section IV-B further
GIC(t) = aEx (t) + bEy (t), (1) expands the ensemble methodology to compute transfer func-
tions straight from B-field to GIC. Both the E-field to GIC and
where a and b are derived constants based on network
B-field to GIC ensemble methods are tested in Section IV-C.
topology and resistances assuming 1 V/km geoelectric (E-
The focus throughout this paper is on operational modelling,
field) components Ex and Ey , where x and y indicates the
with emphasis on estimating the uncertainty associated with
North and East directions respectively. These E-field com-
traditional modelling.
ponents are typically not measured, but rather derived from
II. E NSEMBLE E STIMATION
This work was supported in part by a grant from the Open Philanthropy Traditional GIC modelling recognizes three steps in cal-
Project. The authors would like to further acknowledge Eskom and the culating transformer neutral currents: derivation of local B-
EPRI Sunburst project for measured GIC data in South Africa and similarly
Powerlink Queensland in Australia and Tennessee Valley Authority in the field components from suitable measured or interpolated mag-
USA, for which CTG is a strategic partner. netic measurements; a frequency dependent transform through

21st Power Systems Computation Conference Porto, Portugal — June 29 – July 3, 2020
PSCC 2020
the surface impedance to an E-field; and network analysis. directly with accumulated modelling errors. For the ensemble
Assuming the driving disturbance B-field from near-Earth calculation an assumption is required that the network param-
current systems is spatially uniform and vertically incident at eters are constant over the time period, as is the case in a
the Earth’s surface, along with laterally homogeneous ground resistive network given no change in the network state, e.g.
conductivity, a conservative E-field would be produced through line switching. Such a change would be apparent in modelling
Faraday’s law. Given a purely resistive network, the system can as a multi-modal ensemble. Given a time-series {t1 , ..., tn },
be modelled perfectly by (1). Although these assumptions can repeated calculation of all (i, j) pairs of time instances creates
be justified as approximations, there are several challenges to an ensemble of n(n − 1)/2 ≈ n2 /2 (for large n) network
traditional GIC modelling. Ground conductivity is not laterally parameters. Relatively short time periods can by extension
homogeneous and interfaces, such as the coast, can have a result in very large, statistically meaningful ensembles for real-
significant effect on the magnitude and direction of the induced time analysis of the network state, with no actual network
E-field [5]. Adding that the driving current systems do not information needed. Figure 1 is an example of the α ensemble
produce a uniform B-field over an area the size of a network, at PAR, with the contributions to the total ensemble from
the E-field is not strictly conservative and the transmission line different GIC strengths shown as percentile range profiles.
shape is significant [3]. The network analysis needed to derive
the network parameters is not trivial and the entire network
needs to be taken into account [6]. Individual transformers can
influence each other [7], along with different voltage levels
and possible non-linear inductive network response which is
the topic of very recent research [8], [9]. There are also
variables such as the effect of rainfall on grounding resistance
that add unmodelled complexity [10]. Given the multitude of
higher-order effects, the current state of the art, which makes
use of dense electromagnetic surveys and requires detailed
network information, still does not fully model the nature
of measured GICs. In comparison the traditional modelling
framework does surprisingly well, especially when empirically
‘tuned’ via network parameters or similar scaling factors, and Fig. 1. Ensemble of α network parameter at PAR. Since the effective network
is ideal for operational application. is approximately north-south, β is negligible in comparison. Variation and
change of network parameters with GIC strength is seen in the percentile range
With more and more measured GIC data, the robustness profiles. The total ensemble parameter estimate is indicated by the dashed line,
of the traditional modelling approach can be leveraged using with the associated error signified by the interquartile range (IQR).
data driven approaches. Using measured GIC data overcomes
the fact that the entire network aggregates the driving vari- The ensemble estimation of network parameters presented
ables. The resulting integrated effective response, which may here has been shown to be stable given different ground
arise from complex interactions, can be very different from conductivity profiles [12]. The accuracy of the conductivity
modelling the actual driving variables at the single point of profile used significantly influences the spread or error of
interest. Using such an empirical methodology also allows a the ensemble, while factors mentioned above as challenges
natural error range to be assigned to each step in the modelling to traditional GIC modelling also contribute. One result from
process, or in certain cases the cumulative error of multiple previous analyses [11], [12] is that the empirical network
steps. The ensemble method described in this paper makes parameters change with GIC intensity. Possible contributions
use of actual measurements of GICs, measured B-field and to such observed parameter variation are different parts of a
derived E-field, but can be applied similarly to measured E- geomagnetic storm having different driving regimes affecting
field if available. Firstly, (1) is updated to acknowledge the the spatial homogeneity over the network area and possibly
accumulated errors involved in traditional GIC modelling and the network itself introducing some sort of response. It is
data driven models in general [11]. Now, possible to model this variation in the form of dynamic

network parameters that depend on GIC or E-field magnitude
GIC(t) + GIC(t)err =a Ex (t) + Ex (t)err + [12]. Since the empirical network parameters scale the input


b Ey (t) + Ey (t)err (2) E-field to match the measured GIC, the signal-to-noise ratio
or Γ(t) ≈αEx (t) + βEy (t), (3) (SNR) changes as well, adding to the variation of empirical
network parameter estimates; in this case the variation is away
where Xerr indicates the error made in the measurement or from a Gaussian noise profile with zero mean associated with
estimation of parameter X and Γ(t) is the GIC as measured. GIC at noise levels (as in Figure 1). Although the ensemble
Simultaneous GIC and E-field time-series are used to form method is specific only to a single node-magnetometer pair,
a large set of pairwise combinations of (3). Solving these the estimated network parameters and spread can be used to
many linear equations for α and β set up ensembles of calibrate a traditional network model, as is typically used for
parameter estimates. The spread in the ensembles is associated planning.

21st Power Systems Computation Conference Porto, Portugal — June 29 – July 3, 2020
PSCC 2020
 measured geoelectric field and sparse geomagnetic field mea-
surements (see Table I) are used as an example. Specifically,
GIC data from Tennessee, USA (PAR dataset); Queensland,
Australia (BOW dataset); and the Eastern Cape, South Africa
(GRS dataset) are used. GIC data at all these sites are
measured using a Hall-effect sensor at the transformer neutral,
with no further data cleaning or filtering besides a 1 s timing
drift corrected at PAR. The measured GIC is not the line GIC,
but rather the GIC effective to the particular transformer and
takes into account the entire network.
Corresponding B-field data provided by the nearest INTER-
Fig. 2. Standard deviation and mean of network parameter estimates given MAGNET (www.intermagnet.org) station is used to derive
different dataset lengths for 2 s (PAR) and 1 minute (GRS) cadence data.
Horizontal lines indicate best-fit network parameters, which match the total either the B-field-to-GIC transfer function or the E-field for the
dataset estimates later in Table II. ensemble method. Due to geomagnetic mid-latitude regions
being analysed, modelling can be done effectively without
spatial interpolation of the B-field to the network in ques-
A. Convergence and Statistics of Ensemble Method
tion [14]. Applying an interpolation scheme should decrease
Given that limited GIC data has been available in most the normalised ensemble spread while improving modelling
power networks, there has been a challenge to use data driven marginally. Such interpolation is sidestepped as the ensemble
modelling in the past. This has changed significantly over the method directly relates a GIC site with a magnetometer and
last few years, particularly following the FERC requirement any error being made is absorbed into the ensemble spreads.
for utilities to measure GICs and make these measurements For a utility this provides flexibility of not requiring additional
available [13]. magnetometer coverage for operational modelling. The GRS
Nevertheless, the ensemble method used here makes up for and PAR GIC data have 2 s cadence and BOW has a minimum
the lack of data by creating an ensemble of roughly n2 /2 sampling interval of 4 s during active periods. The associated
estimates for a dataset of given n. Each of these estimates is B-field data for the 2015 datasets are at the modern 1 s
a possible state of the system, with the peak of the ensemble cadence, compared to the 2003 dataset which has 60 s B-
being the most probable. As the ensembles are heavy-tailed, field data (often still the best cadence many locations have).
typical statistics may be misleading. For estimation purposes, Appropriate resampling has been done on the concurrent
the median is inherently robust and the best measure of datasets to reflect the lowest cadence for that period.
central tendency. Similarly, the interquartile range (IQR) is
an appropriate estimate for the spread in the ensemble. In TABLE I
Figure 2, the convergence of the ensemble method is shown DATASETS USED
for 2 s and 1 minute cadence data, given different dataset Dataset Type Timespan Cadence
lengths. For each case, the mean and standard deviation of GRS GIC 31/03/2001, 29-31/10/2003 [UTC] 2s
25 model runs are plotted, with the best-fit estimate shown BOW GIC 23/06/2015 [AEST]/[UTC+10] 4s
as a horizontal line. In each case the final best-fit estimate PAR GIC 22-23/06/2015 [UTC] 2s
HER MAG Same as GRS 60 s
matches the estimate given by the using the entire training set
CTA MAG Same as BOW 1s
in a single run. The 1 minute data is smoother and converges
FRD MAG Same as PAR 1s
more quickly, with variation in estimates being less than 5%
of total network parameters at around 105 estimates. The For validation, 25% of each dataset, including the largest
higher resolution 1 s data, with more variability, has similar GICs associated with the spike-like sudden storm commence-
adequate convergence at around 107 estimates. This relates ment (SSC) and main phase of intense geomagnetic storms, is
to roughly a day’s worth of data at either cadence, which kept out-of-sample. The models are trained on the remaining
means that operationally a single day of data will suffice lower amplitude GIC data with a lower SNR. Using the ‘best’
in characterising the real-time effective state of the network data available in each period to train the models on has
for GIC modelling. Defining the optimal dataset length for deliberately been avoided, with the relatively ‘weak’ model
convergence furthermore allows the ensemble method to be results being used to test the model on the validation set. Such
applied using limited resources - a 2 billion estimate ensemble an approach with limited data demonstrates conservatively
requires significant computing power for processing. If needed, what is possible in real world conditions. Given more data for
the optimal dataset lengths can be further reduced by a factor multiple storms, with more representative training throughout
of 10 if multiple runs are averaged such as in Figure 2. the storm phases, will improve modelling.
III. DATA A. E-field Derivation
In the current work three sets of results, using real world The ground conductivity (or resistivity) model is an impor-
conditions with limited GIC data, no network information, no tant contribution to the E-field. Although grounding (earthing)

21st Power Systems Computation Conference Porto, Portugal — June 29 – July 3, 2020
PSCC 2020
models are widely used in power system analysis, they are
usually restricted to power frequency models for fault analysis
or high frequency models for lightning studies. GICs are
characterized largely by the low frequency spectrum below
50 mHz, for which the skin depth is very deep, conceptually
over 50 km [15].
Although there are a few sites worldwide where E-field
measurements are recorded continuously, it is more typical
to derive the E-field from the B-field. Extending the initial
assumption that the geomagnetic disturbance field is a plane-
wave and that the conductivity of Earth solely depends on
depth, we can make use of the basic magnetotelluric (MT)
Fig. 3. Comparison of concurrent GIC measurements in the USA (blue)
equation [16]. This equation relates the horizontal components and Australia (red) with derived E-fields at the sudden storm commencement
of the B-field Bx,y , to the induced E-field Ex,y in the (SSC). The E-fields are projected onto the typical local network direction, i.e.
frequency domain, south for PAR and south-west for BOW. The scaling factor of 1/2 for the
  CTA E-field matches the E-field axis scaling with measured GIC at BOW.
~ Z(ω) ~ Zxx (ω) Zxy (ω)
E(ω) = B(ω) where, Z(ω) = .
µ0 Zyx (ω) Zyy (ω)
(4) the initial shock front of solar plasma compresses Earth’s B-
field and enhances it, followed by an intensification of the
Here, Z(ω) is the complex frequency dependent surface magnetospheric ring current which opposes the B-field and
impedance, often quoted in units of [ mV 1
km nT ], i.e. the ratio of weakens it (known as a geomagnetic storm), we also see the
electric and magnetic fields in conventional units. The most lag in B-field and hence E-field between BOW and PAR.
basic form of the MT equation assumes uniform conductivity. For this event PAR was on the Sun-ward side and felt the
This basic conductivity profile assumption can be extended brunt of the SSC before BOW. At 1 s cadence, the onset
to a layered-Earth model where conductivity in each layer is of a geomagnetic storm and GIC peak are not concurrent or
assumed to be uniform and defined solely by depth, which instantaneous at mid-latitudes globally as is often assumed. As
produces relative frequency scaling that more accurately de- a further illustration of the theoretical basis of (1), using simple
fines the Earth’s inductive filtering. For both cases, only the autoscaling of axes and projecting the E-field onto the general
off-diagonal impedance tensor components are non-zero, with network direction to mimic the role of network parameters, the
Zxy = −Zyx . When there are more complicated geological resulting measured GIC and derived general E-field profiles in
structures, such as in coastal regions there is a definite lateral Figure 3 show strikingly close correlation.
discontinuity or strike, the impedance tensor also becomes
more complicated, ultimately with each component unique for IV. R ESULTS
a measured surface impedance tensor [17].
A. Network Parameters and Characterisation
Since the ensemble method is stable regardless of the
conductivity profile, which is not typically known unless a Besides the time-series modelling capability provided by
MT survey has been done, we make use of a representative the ensemble methodology, there is the added impact of
global average layered-Earth to obtain physically relevant finding the effective network response. The resulting effective
frequency scaling [18]. This global average profile has layer response requires no prior network modelling, and can absorb
thicknesses of d = [40, 210, 160, 260, 230, 1300, 500] the inherent uncertainty in the multi-step modelling. Such
km and corresponding layer conductivities of uncertainty, seen as deviations from the traditional modelling
σ = [0.0056, 0.0095, 0.0262, 0.0776, 0.526, 1.69, 10] S/m. approach (1), takes the form of a range (the interquartile
The conductivity of the terminating half-space is 100 S/m. range in this case) about the most effective network parameters
(summarised in Table II). The network parameter range reflects
B. Data Cadence Implications an upper and lower bound for modelling purposes. Given the
Before modern 1 s cadence B-field and hence E-field data, large number of estimates, the ensembles can be broken down
the resistive approximation was not questioned at all since into percentile range profiles associated with GIC strength as
with 60 s data there is no observable time lag between GIC seen in Figure 1. These profiles show that the driving regimes
and E-field. Using the global conductivity profile described for different levels of GICs are often quite different, possibly
above, Figure 3 shows the higher resolution structure of the associated with higher-order effects, and that the empirical
profiles. An autocorrelation analysis over the geomagnetically estimates are dynamic. Practically, the percentile ranges can
active period shows BOW has 12 second delay and PAR has be used to further improve modelling, using estimated GIC
10 second delay between the E-field and the resultant GIC. to define the percentile range regime and update the GIC
By varying the Earth conductivity model, these lags may or estimation [12].
may not change, but a contribution from the power network From the empirical ensembles, the effective network direc-
cannot be ruled out. In Figure 3, which shows the SSC when tionality can also be defined. If there is a strike of sorts in

21st Power Systems Computation Conference Porto, Portugal — June 29 – July 3, 2020
PSCC 2020
 TABLE II conductivity assumption required to derive the E-field, which
S ITE SPECIFIC NETWORK PARAMETERS AND DIRECTIONALITY introduces significant error, is bypassed. The derived transfer
Site Estimates α [A/(V/km)] β [A/(V/km)] c Bearing function (TF) is the net effective transform of the Earth over
GRS ± 5 million -109.18 -1.09 0.01 181◦ the entire network area and any power network response. The
BOW ± 131 million -211.99 -167.58 0.79 218◦ B-field-to-GIC TF approach has been followed before more
PAR ± 2.1 billion -185.82 -12.97 0.07 184◦ from a traditional MT background, without decomposing the
network filtering [20]. Although the results have shown good
the ground conductivity structure, the effective directionality correspondence with known geophysical features, there are
will absorb it. When rotating a constant B-field, i.e. circular cases where such interpretations failed. When these interpre-
polarisation, a strike results in an elliptical polarised E-field tations fail, it can largely be attributed to the results of a GIC
which effectively induces larger GICs along a particular di- based TF taking into account the network. To take into account
rection. If a simple layered-Earth conductivity model is used, the effect of the network we firstly, update (3) with c = β/α.
the derived E-field would not have this characteristic and The matrix form of this equation becomes,
there would be underestimation of modelled E-field along Γ(t) = αEx (t) + αcEy (t)
with GIC if the effective network is aligned to the E-field.   
α αc Ex (t)
Using a network parameter ensemble, the additional E-field = . (5)
Ey (t)
required for accurate modelling is implicitly encoded in the
empirical network parameters and can be seen in the effective In the frequency domain, we use measured B-field data without
network direction. Other non-geophysical contributions to the the need for an assumed conductivity related impedance tensor
effective directionality, such as contributions from the network and the uncertainty that goes with it,
at large, would also be included. The effective directionality   
α αc Ex (ω)
argument comes from the arctangent of the ratio of network Γ(ω) =
Ey (ω)
parameters c = β/α [19] and is summarised with its associated    
bearing in Table II. The effective directionality can be seen α αc Zxx (ω) Zxy (ω) Bx (ω)
=
as a planning tool since alignment of the effective weighted Zyx (ω) Zyy (ω) By (ω)

T 
network contribution at a node, taking into account the entire
 
α Zxx (ω) + cZyx (ω)
Bx (ω)
network with its shape [3], [6], with the input E-field results in =
α Zxy (ω) + cZyy (ω) By (ω)
the highest GICs where measured, i.e. the transformers. The   
T (ω) Ty (ω) Bx (ω)
relation with actual local transmission line direction between = x (6)
By (ω)
nodes is not as trivial as often made out. An example of a
non-trivial directionality ensemble derived using each estimate =TFBx (ω) + TFBy (ω). (7)
from the global ensembles is shown for BOW in Figure 4. In the formalism above, Tx,y (ω) in (6) are the components of
Although it is the ‘true’ E-field that induces the measured the B-field-to-GIC TF. TFBx ,By (ω) in (7) are not TF compo-
GICs, with the ensemble method, given a consistent system nents, but rather the result of multiplying Tx,y (ω) with their
state, we have a direct association between an approximate associated B-field components. To aid in identifying effective
input E-field direction and high transformer GICs. contributions of the surface impedance tensor (4) in the TF,
the geometric scaling c is split from the network parameters.
The final TF absorbs all network filtering of α (which cannot
be separated) and can be directionally limited given certain
effective network topologies (when c or c−1 ≈ 0), severely
affecting the SNR. This is seen in Figure 5 for the GRS TF;
the north-south effective orientation of the network dampens
Tx , and by extension the contribution from Bx . In a resistive
network, α is assumed to be constant. This may not be the
case, and α may in fact have frequency dependent scaling [4].
The TF approach has a number of inherent advantages:
a. The measured B-field is directly related to measured GIC
without any further assumptions of ground or network
Fig. 4. The effective ensemble directionality at BOW for different levels response. Historically, the B-field is much more com-
of GIC strength along with the bearing of the immediate local line segment monly measured and better understood than the E-field.
indicated by the dashed lines. The long term B-field analyses allow for better extreme
value statistical modelling important in planning.
B. Empirical Transfer Functions b. The integrated ground and network effects of the entire
It is also possible to directly relate the B-field to GIC via an network are modelled for the node. Any variation in these
empirical transfer function. Through this process the ground responses are encoded in the TF spread.

21st Power Systems Computation Conference Porto, Portugal — June 29 – July 3, 2020
PSCC 2020
Fig. 5. Transfer function response at GRS, with magnitude and phase of the Tx,y components shown. Since the effective network orientation is north-south
and c ≈ 0, only the single off-diagonal magnetotelluric impedance Zxy dominates the responses, with Zyx (and hence Tx ) being swamped by noise.

c. Different regimes of driving geomagnetic disturbances variation. Given a longer dataset it can be assumed that the
have different spectral components, which are accurately TF will become significantly more accurate.
modelled in the frequency domain along with geophysical Uncertainty in the time domain reconstruction is not trivial
anomalies and network responses. The ensemble method since the frequency specific uncertainty bounds can occur in
requires further processing for similar modelling, and any combination, constructively or destructively. To obtain
failing that can only increase its uncertainty spread. an uncertainty estimate in the time domain, the time domain
d. The TF estimates are statistically and physically mean- ensemble method from earlier is reintroduced to generate time
ingful, i.e. they allow for probabilistic GIC estimation domain scaling and uncertainty bands. Taking the IFFT of the
and interpretation of the geophysical and network contri- TF result TFBx ,By (ω) in (7), we can use time domain version
butions to the TF. of the same result (denoted by TF*) analogously to the E-field-
Shortcomings of using a TF approach are the need for to-GIC ensemble methodology,
representative training data and correlation in time between
geomagnetic disturbances at the B-field measurement site and
Γ(t) = gTF∗Bx (t) + hTF∗By (t). (8)
over the effective network. Both these issues are apparent at
PAR (discussed in IV-C), in which case the more robust E-field
ensemble approach does better. The g and h parameters from (8) generate an uncertainty
In order to estimate the TF, an ensemble method in the spread and improve modelling by tuning TF* or the analogous
frequency domain is used. For a number of windows, the ‘network effective E-field’ in the time domain. Here the g and
two complex TF components are estimated at each frequency h ensembles are centred roughly around 1, which is expected
through many simultaneous equations calculations (negative since the empirical TF should correctly scale the input B-
frequencies are folded). These ensembles are then log-binned field to GIC. Any significant deviation from unity scaling
and the median and interquartile range estimated for modelling suggests the TF does not adequately link B-field to GIC, such
purposes as with the previous ensemble method. The limited as when the separation between magnetometer and substation
data in this work is a challenge since there may well be is too great. For reference, the separation of PAR-FRD pair
different spectra to consider and the SNR is poor at many is extreme at around 1000 km whereas the GRS-HER pair
frequencies. As a result only the GRS and PAR datasets were is roughly 700 km apart. When further scaling is needed,
used for TF modelling. The GRS TF used an 18 hour window, the ensemble TF can further improve modelling. In either
with 3 hour shifting, at 60 s sampling and the PAR TF used a case of an ideal TF with unity scaling or tweaked TF with
6 hour window, with 2 hour shifting at 2 s sampling. Both further scaling, the spread of the g and h parameters allow for
dataset processing schemes give a few thousand estimates accurate uncertainty estimation in the time domain and not the
per frequency, which is then significantly increased with log- frequency domain where the TF is defined. The conceptual
binning. However, at PAR the short window length may flow of the B-field-to-GIC TF and TF* ensembles and E-field
have missed certain typical lower frequencies such as daily ensembles are summarised in Figure 6.

21st Power Systems Computation Conference Porto, Portugal — June 29 – July 3, 2020
PSCC 2020
 region, suggesting the spatial distribution models for
geomagnetic disturbances need to be reviewed.
b. They are non-stationary and localised, which violates the
plane-wave assumption of the disturbance B-field. Due
to the separation between sites, during the period of
pulsations there is an offset between the measured B-
Fig. 6. Flow of the traditional ensemble updated layered-Earth (LE), transfer
function (TF) and ensemble updated transfer function (TF*) models used. field at FRD and the B-field driving the GIC at PAR.
Although the offset distorts the models, particularly the
TF’s, inspection shows the results are useful. In practice,
C. GIC Modelling with more data being accumulated and spatial interpola-
tion, the distortion seen in modelling may be minimised.
The results of the layered-Earth (LE) ensemble, TF and
TF* models discussed are summarised in Table III. Of the
measures used, the RMSE and correlation coefficient ρ are
typical. The relative error (RE) has been used in previous
studies and quantifies the percentage error made in terms of the
signal amplitude. Considering the GIC sites, only GRS [21]
and BOW [22] have previously been modelled at all and only
GRS has comparable performance metrics for the validation
set. In previous modelling the best result at GRS for the
same validation set used finite element modelling, obtaining
a RMSE of 0.98A and RE of 42% compared to 0.69A and
25.1% using the TF* model in this work.
Fig. 7. The PAR GIC time-series, with the dominant FRD By component
(c ≈ 0), shows significant driving from long-period geomagnetic pulsations
TABLE III resulting in sustained GICs of similar amplitude to the SSC spike.
M ODELLING RESULTS FOR OUT- OF - SAMPLE DATA
Site (Model) RMSE (A) ρ RE (%) Nevertheless, the presented modelling is viable from an oper-
GRS (LE) 0.91 0.86 32.3 ational outlook especially with the uncertainty band tracking
BOW (LE) 1.24 0.92 32.0 the error introduced, as shown for GRS in Figure 8.
PAR (LE) 1.23 0.89 44.6
GRS (TF) 0.81 0.86 29.4 V. C ONCLUSION AND F URTHER R ESEARCH
GRS (TF*) 0.69 0.91 25.1
PAR (TF) 1.49 0.81 57.3
The ensemble methodology has shown how the relatively
PAR (TF*) 1.32 0.86 52.5 simple governing nodal GIC equation can be leveraged given
measured data to represent much more complicated GIC mod-
elling dynamics - even with relatively few, spatially sparse,
When dealing with higher amplitude GIC and associated
measurement sites over a large area and limited temporal
higher amplitude noise, such as at PAR, the RE measure
coverage. The usefulness of the methodology is increased
can be significantly skewed. BOW and GRS have similar
by uncertainty being included, even in cases of extreme
magnitude GICs, making their RE results more comparable.
separation between B-field and GIC measurement. In a utility,
The RE measure can be skewed even further with limited
the ensemble and TF methods can be used in control centres
data used for modelling; the dynamic parameters used may
during or before a geomagnetic disturbance, giving direct
change significantly given a geomagnetic storm profile or
indication of simulated GIC exposure of several transformers
amplitude different (usually higher) to data previous trained
(providing there have been prior measurements of GIC) and
on. This is seen in the current work as an underestimation in
identifying discrepancies in state of the network, all without
modelling particularly for BOW, where the out-of-sample SSC
overloading the communication systems. In this regard the data
is significantly higher than anything in the training set.
driven methods already demonstrate the effectiveness of FERC
PAR presents a further modelling challenge since the geo- order 830 for the collection of data. Under this directive, it is
magnetic field used in training has a relatively high content of foreseeable that data driven methods become even more viable.
long-period (20 min) geomagnetic pulsations that are usually The associated B-field measurements considered do not need
associated with geomagnetic substorms and auroral regions 1 s cadence to be operationally effective, although they do al-
[23]. During a 2 hour period, the pulsations caused high am- low for high resolution GIC dynamics to be seen, such as SSC
plitude GICs, comparable to the SSC peak GIC (see Figure 7). propagation across the globe. From a planning perspective, the
These driving pulsations are interesting in two respects: ensemble method can be used to test transformer level models
a. They indicate geomagnetic activity associated more often of uneven GIC distributions between transformers in the same
with auroral regions. Here significant effects are seen substation [7] and calibrate general network models further.
in what is generally considered to be a mid-latitude The significant impact of driving geomagnetic pulsations does

21st Power Systems Computation Conference Porto, Portugal — June 29 – July 3, 2020
PSCC 2020
Fig. 8. Time-series of the resulting GIC modelling at GRS for the out-of-sample dataset, including the associated uncertainty bands. More than 80% of data
points fall within the TF* prediction band compared to 70% for the LE prediction band.

however suggest that the definition of mid-latitude GIC drivers [11] M. Wik et al., “Calculation of geomagnetically induced currents in the
and the plane-wave spatial scale be reconsidered; given the 400 kV power grid in southern Sweden,” Space Weather, vol. 6, no. 7,
pp. 1–11, 2008.
expansion of high-latitude effects, the currently understood [12] M. J. Heyns, S. I. Lotz, P. J. Cilliers and C. T. Gaunt, “Ensemble
risk to mid-latitude power systems may be underestimated. Estimation of Network Parameters: A Tool to Improve the Real-time
Estimation of GICs in the South African Power Network,” in The
R EFERENCES Proceedings of SAIP2017, the 62nd Annual Conference of the South
African Institute of Physics, edited by J. Engelbrecht, Stellenbosch, pp.
[1] V. Albertson and J. Van Baelen, “Electric and Magnetic Fields at the 270–275, 2017. Available online at http://events.saip.org.za
Earth’s Surface Due to Auroral Currents,” IEEE Transactions on Power [13] “Federal Energy Regulatory Commission: Reliability Standard for Trans-
Apparatus and Systems, vol. PAS-89, no. 4, pp. 578–584, 1970. mission System Planned Performance for Geomagnetic Disturbance
[2] M. Lehtinen and R. Pirjola, “Currents produced in earthed conductor Events.” Order 830, Sep 2016. Washington DC.
networks by geomagnetically-induced electric fields,” Annales Geophys- [14] C. M. Ngwira et al., “Limitations of the modeling of geomagnetically
icae, vol. 3, no. 4, pp. 479–484, 1985. induced currents in the South African power network,” Space Weather,
[3] R. Sun and C. Balch, “Comparison between 1-D and 3-D Geoelectric vol. 7, no. 10, pp. 1–5, 2009.
Field Methods to Calculate Geomagnetically Induced Currents: A Case [15] D. Oyedokun, M. J. Heyns, P. J. Cilliers and C. T. Gaunt,
Study,” IEEE Transactions on Power Delivery, vol. 34, no. 6, pp. 2163– “Frequency Components of Geomagnetically Induced Currents
2172, 2019. for Power System Modelling,” 2020 IEEE International
[4] R. S. Weigel and P. J. Cilliers, “An Evaluation of the Frequency SAUPEC/RobMech/PRASA Conference, pp. 1–6, 2020.
Independence Assumption of Power System Coefficients Used in Ge- doi:10.1109/SAUPEC/RobMech/PRASA48453.2020.9041021
omagnetically Induced Current Estimates,” Space Weather, vol. 17, [16] L. Cagniard, “Basic Theory Of The MagnetoTelluric Method Of Geo-
no. 12, pp. 1674–1688, 2019. physical Prospecting,” Geophysics, vol. 18, no. 3, pp. 605–635, 1953.
[5] G. M. Lucas, J. J. Love and A. Kelbert, “Calculation of Voltages [17] D. H. Boteler and R. J. Pirjola, “Modeling geomagnetically induced
in Electric Power Transmission Lines During Historic Geomagnetic currents,” Space Weather, vol. 15, no. 1, pp. 258–276, 2017.
Storms: An Investigation Using Realistic Earth Impedances,” Space [18] J. Sun, A. Kelbert and G. D. Egbert, “Ionospheric current source
Weather, vol. 16, no. 2, pp. 185–195, 2018. modeling and global geomagnetic induction using ground geomag-
[6] T. J. Overbye, K. S. Shetye, T. R. Hutchins, Q. Qiu, and J. D. Weber, netic observatory data,” Journal of Geophysical Research: Solid Earth,
“Power Grid Sensitivity Analysis of Geomagnetically Induced Currents,” vol. 120, no. 10, pp. 6771–6796, 2015.
IEEE Transactions on Power Systems, vol. 28, no. 4, pp. 4821–4828, [19] A. A. Pulkkinen, R. Pirjola, and A. Viljanen, “Determination of ground
2013. conductivity and system parameters for optimal modeling of geomag-
[7] T. Divett et al., “Transformer-Level Modeling of Geomagnetically In- netically induced current flow in technological systems,” Earth, Planets
duced Currents in New Zealand’s South Island,” Space Weather, vol. 16, and Space, vol. 59, no. 9, pp. 999–1006, 2007.
no. 6, pp. 1–18, 2018. [20] M. Ingham et al., “Assessment of GIC based on transfer function
[8] H. K. Chisepo, C. T. Gaunt and L. D. Borrill, “Measurement analysis,” Space Weather, vol. 15, no. 12, pp. 1615–1627, 2017.
and FEM analysis of DC/GIC effects on transformer magnetiza- [21] E. Matandirotya, P. J. Cilliers and R. R. Van Zyl, “Modeling geomagneti-
tion parameters,” 2019 IEEE Milan PowerTech, pp. 1–6, 2019. cally induced currents in the South African power transmission network
doi:10.1109/PTC.2019.8810423 using the finite element method,” Space Weather, vol. 13, no. 3, pp.
[9] P. Jankee et al., “Transformer models and meters in MATLAB 185–195, 2015.
and PSCAD for GIC and leakage dc studies,” 2020 IEEE In- [22] R. A. Marshall et al., “Modelling Geomagnetically Induced Currents in
ternational SAUPEC/RobMech/PRASA Conference, pp. 1–6, 2020. Australian power networks using different conductivity models,” Space
doi:10.1109/SAUPEC/RobMech/PRASA48453.2020.9041060 Weather, vol. 17, no. 5, pp. 727-756, 2019.
[10] S. P. Blake et al., “A Detailed Model of the Irish High Voltage Power [23] M. J. Heyns, S. I. Lotz and C. T. Gaunt, “Geomagnetic Pulsa-
Network for Simulating GICs,” Space Weather, vol. 16, no. 11, pp. tions Driving Geomagnetically Induced Currents,” ESSOAr (preprint on
1770–1783, 2018. https://essoar.org/), 2020. doi:10.1002/essoar.10503394.1

21st Power Systems Computation Conference Porto, Portugal — June 29 – July 3, 2020
PSCC 2020