$E_{0}(g,s,\Phi)=\Phi(g,s)+\mathbb{B}(g,s,\Phi)+W_{0}(g,s,\Phi)$ (308)

where

$\mathbb{B}(g,s,\Phi)=\sum_{\gamma\in\Gamma_{\infty}\backslash\mathrm{SL}_{2}(\mathbb{Z})}B(m(\gamma)g,s,\Phi)$
(309)

with

$B(g,s,\Phi)=W_{0}\left(e,s+\tfrac{1}{2},\left(\eta^{*}\circ
r(g)\right)\Phi\right).$ (310)

A similar argument to the proof of the Proposition 4.19 shows that

$\mathbb{B}(g_{v},0,\Phi^{(\ell)})=\mathbb{B}(g_{v},0,\Phi^{sp,(\ell)})=0$
(311)

and

$\mathrm{ord}_{s=0}\mathbb{B}(g_{v},s,\Phi^{L})\geq 2.$ (312)

Moreover, if $\Phi^{*}_{2}$ is any of the sections $\Phi^{\mathcal{L}}_{2}$,
$\Phi^{L}_{2}$, $\Phi^{(\ell)}_{2}$ or $\Phi^{sp,(\ell)}_{2}$, the definitions
imply the formula

$\Phi^{*}_{2}(g_{v},s)=\det(v)^{\frac{s}{2}+\frac{3}{4}}.$ (313)

Thus, we find

$E^{\prime}_{0}(g_{Nv},0,\Phi^{\mathcal{L}}_{2})=N^{\frac{3}{2}}\det(v)^{\frac{3}{4}}\left(\log
N+\frac{1}{2}\log\det
v\right)+W^{\prime}_{0}(g_{Nv},0,\Phi_{2}^{\mathcal{L}}),$ (314)

and

$E^{\prime}_{0}(g_{v},0,\Phi^{L}_{2})=\det(v)^{\frac{3}{4}}\left(\frac{1}{2}\log\det
v\right)+W^{\prime}_{0}(g_{v},0,\Phi_{2}^{L});$ (315)

in addition, since $\Phi^{L}_{2}$ and $\Phi^{sp,(\ell)}$ are incoherent, we
have

$E_{0}(g_{v},0,\Phi_{2}^{L})=E_{0}(g_{v},0,\Phi_{2}^{sp,(\ell)})=0$ (316)

so

$W_{0}(g_{v},0,\Phi^{L})=W_{0}(g_{v},0,\Phi^{sp,(\ell)})=-\det(v)^{\frac{3}{4}}.$
(317)

Thus, Proposition 5.1 implies that

$\displaystyle N^{-\frac{3}{2}}$ $\displaystyle
E^{\prime}_{0}(g_{Nv},0,\Phi_{2}^{\mathcal{L}})-E^{\prime}_{0}(g_{v},0,\Phi_{2}^{L})-\sum_{\ell|N}\frac{\ell-1}{2(\ell+1)}E_{0}(g_{v},0,\Phi^{(\ell)}_{2})\log\ell$
(318) $\displaystyle=\det(v)^{\frac{3}{4}}\log
N+\left(\sum_{\ell|N}\frac{\ell+1}{2(\ell-1)}\log\ell\right)W_{0}(g_{v},0,\Phi_{2}^{L})$
$\displaystyle\qquad-\sum_{\ell}\frac{2}{\ell^{2}-1}W_{0}(g_{v},0,\Phi_{2}^{sp,(\ell)})\log\ell-\sum_{\ell}\frac{\ell-1}{2(\ell+1)}\log\ell\cdot\det(v)^{\frac{3}{4}}$
(319)
$\displaystyle=\det(v)^{\frac{3}{4}}\left(\sum_{\ell}1-\frac{\ell+1}{2(\ell-1)}+\frac{2}{\ell^{2}-1}-\frac{\ell-1}{2(\ell+1)}\right)\log\ell$
(320) $\displaystyle=0.$

This implies the proposition for the case $T=0$. ∎

######  Remark 5.4.

The preceding proposition can be phrased in more classical language, as
follows, cf. Section 2.6. If $\Phi^{*}_{2}$ is one of the sections
$\Phi^{\mathcal{L}}_{2},\Phi^{L}_{2}$ or $\Phi^{(\ell)}_{2}$, and
$\tau\in\mathbb{H}_{2}$, we set

$E(\tau,s,\Phi^{*}_{2})=\det(v)^{-\frac{3}{4}}E(g_{\tau},s,\Phi^{*}_{2}).$
(321)

Let

$\Gamma_{0}(4N)=\left\\{\begin{pmatrix}A&B\\\
C&D\end{pmatrix}\in\operatorname{Sp}_{r}(\mathbb{Z})\ |\ C\equiv
0\pmod{4N}\right\\}$ (322)

denote the usual congruence subgroup. Then it can be verified that
$E(\tau,s,\Phi^{*}_{2})$ transforms like a Siegel modular form of scalar
weight $3/2$ and level $\Gamma_{0}(4N)$ with character
$\chi(\gamma)=\mathrm{sgn}\det(d)$, see [KRY06, §8.5.6] for a more precise
formulation.

Let $M_{\frac{3}{2}}(\Gamma_{0}(4N))$ denote the space of (non-holomorphic)
Siegel modular forms of scalar weight $3/2$ and level $\Gamma_{0}(4N)$ and
character $\chi$, and consider the Hecke operator

$U_{N}\colon M_{\frac{3}{2}}(\Gamma_{0}(4N))\to
M_{\frac{3}{2}}(\Gamma_{0}(4N))$ (323)

given by

$U_{N}(F)(\tau)=\sum_{u\in\operatorname{Sym}_{2}(\mathbb{Z}/N\mathbb{Z})}F|_{3/2}\left(\begin{smallmatrix}1&u\\\
&N\end{smallmatrix}\right)=N^{\frac{3}{2}}\sum_{u}F(N^{-1}\tau+N^{-1}u)$ (324)

Then the previous proposition can be rewritten in classical terms as the
identity

$E^{\prime}(\tau,0,\Phi^{\mathcal{L}}_{2})=U_{N}\left(E^{\prime}(\tau,0,\Phi^{L}_{2})+\sum_{\ell|N}\frac{\ell-1}{2(\ell+1)}E(\tau,0,\Phi^{(\ell)})\log\ell\right).$
(325)

## References

  * [BGKK07] J. I. Burgos Gil, J. Kramer, and U. Kühn.  Cohomological arithmetic Chow rings.  J. Inst. Math. Jussieu, 6(1):1–172, 2007.
  * [BY09] Jan Hendrik Bruinier and Tonghai Yang.  Faltings heights of CM cycles and derivatives of $L$-functions.  Invent. Math., 177(3):631–681, 2009.
  * [DY19] Tuoping Du and Tonghai Yang.  Arithmetic Siegel-Weil formula on $X_{0}(N)$.  Adv. Math., 345:702–755, 2019.
  * [GK93] Benedict H Gross and Kevin Keating.  On the intersection of modular correspondences.  Inventiones mathematicae, 112(1):225–245, 1993.
  * [GQT14] Wee Teck Gan, Yannan Qiu, and Shuichiro Takeda.  The regularized Siegel-Weil formula (the second term identity) and the Rallis inner product formula.  Invent. Math., 198(3):739–831, 2014.
  * [GS19] Luis E. Garcia and Siddarth Sankaran.  Green forms and the arithmetic Siegel-Weil formula.  Invent. Math., 215(3):863–975, 2019.
  * [How20] Benjamin Howard.  Arithmetic volumes of unitary Shimura curves.  arXiv e-prints, page arXiv:2010.07362, October 2020.
  * [Kit83] Yoshiyuki Kitaoka.  A note on local densities of quadratic forms.  Nagoya Math. J., 92:145–152, 1983.
  * [KM85] Nicholas M. Katz and Barry Mazur.  Arithmetic moduli of elliptic curves, volume 108 of Annals of Mathematics Studies.  Princeton University Press, Princeton, NJ, 1985.
  * [KR88] Stephen S. Kudla and Stephen Rallis.  On the Weil-Siegel formula.  J. Reine Angew. Math., 387:1–68, 1988.
  * [KR94] Stephen S. Kudla and Stephen Rallis.  A regularized Siegel-Weil formula: the first term identity.  Ann. of Math. (2), 140(1):1–80, 1994.
  * [KRY04] Stephen S. Kudla, Michael Rapoport, and Tonghai Yang.  Derivatives of Eisenstein series and Faltings heights.  Compos. Math., 140(4):887–951, 2004.
  * [KRY06] Stephen S. Kudla, Michael Rapoport, and Tonghai Yang.  Modular forms and special cycles on Shimura curves, volume 161 of Annals of Mathematics Studies.  Princeton University Press, Princeton, NJ, 2006.
  * [Kud97] Stephen S. Kudla.  Central derivatives of Eisenstein series and height pairings.  Ann. of Math. (2), 146(3):545–646, 1997.
  * [Küh01] Ulf Kühn.  Generalized arithmetic intersection numbers.  J. Reine Angew. Math., 534:209–236, 2001.
  * [Rib90] K. A. Ribet.  On modular representations of ${\rm Gal}(\overline{\bf Q}/{\bf Q})$ arising from modular forms.  Invent. Math., 100(2):431–476, 1990.
  * [RR93] R. Ranga Rao.  On some explicit formulas in the theory of Weil representation.  Pacific J. Math., 157(2):335–371, 1993.
  * [Sie80] Carl Ludwig Siegel.  Advanced analytic number theory, volume 9 of Tata Institute of Fundamental Research Studies in Mathematics.  Tata Institute of Fundamental Research, Bombay, second edition, 1980.
  * [Sou92] C. Soulé.  Lectures on Arakelov geometry, volume 33 of Cambridge Studies in Advanced Mathematics.  Cambridge University Press, Cambridge, 1992.  With the collaboration of D. Abramovich, J.-F. Burnol and J. Kramer.
  * [Yan98] Tonghai Yang.  An explicit formula for local densities of quadratic forms.  J. Number Theory, 72(2):309–356, 1998.
  * [Yan04] Tonghai Yang.  Local densities of 2-adic quadratic forms.  J. Number Theory, 108(2):287–345, 2004.

Department of Mathematics, University of Manitoba
E-mail address<EMAIL_ADDRESS>
Department of Mathematics, University of Wisconsin, Madison
E-mail address<EMAIL_ADDRESS>
Department of Mathematics, University of Wisconsin, Madison
E-mail address<EMAIL_ADDRESS>


# Global transformer overheating from geomagnetic storms

Morgan Rivers  0000-0001-9367-8280 Alliance to Feed the Earth in Disasters
(ALLFED), 603 S. Public Rd #57 Lafayette, CO 80026, USA<EMAIL_ADDRESS>Łukasz G. Gajewski  0000-0003-3097-0131 Alliance to Feed the Earth in
Disasters (ALLFED), 603 S. Public Rd #57 Lafayette, CO 80026, USA  David
Denkenberger  0000-0002-6773-6405 Alliance to Feed the Earth in Disasters
(ALLFED), 603 S. Public Rd #57 Lafayette, CO 80026, USA  Department of
Mechanical Engineering, University of Canterbury, Christchurch, Canterbury
8041, NZ

###### Abstract

Geomagnetic storms occurring due to sustained, high-speed solar winds are
known to induce currents in power distribution networks. These geomagnetically
induced currents (GICs) can cause high voltage transformers (HVT) to overheat,
thus resulting in a catastrophic electricity loss (CEL). Since significant
portions of infrastructures around the world rely heavily on access to
electric power, it is essential to estimate the risks associated with GICs on
a global scale. We assemble multiple methodologies across various scientific
disciplines to develop a framework assessing the probability of a severe
geomagnetic storm causing a long-term, widespread power outage. Our model
incorporates thermal models of HVT tie bar hot spots, historical geoelectric
field estimates, and a global conductivity model to estimate the risk of long-
term power outage for regions between $-70^{\circ}$ and $80^{\circ}$
geomagnetic latitude due to transformer overheating failure. Assuming a
uniform $33\%$ HVT spare capacity, our analysis indicates that a 1 in 10,000
year storm would result in approximately $1\%$ of the population in Europe and
North America experiencing a long-term (months to years) electricity loss.

geomagnetic storms, power grid resilience, catastrophic electricity loss,
critical infrastructure, interdisciplinary systems modelling, global
catastrophic risk

††preprint: APS/123-QED

##  I Introduction

The phenomenon known as the geomagnetic storms has been studied for hundreds
of years[1], but only recently, due to advances in electrification, some are
starting to worry about the potential risks they can pose to our society [2,
3]. A society that is increasingly dependent on electricity and related
technologies [2, 4].

Geoelectric Fields for Severe Storm • Historical $B$ field levels. • Apparent
conductivity and resistivity at TF sites. • $E$ field and rates of storm
recurrence. • Log-normal fit to adjusted $E$ field levels to an equivalent
field at a reference location. • Global $E$ field levels. GIC and Transformer
Overheating in the High Voltage Power Distribution Network • GICs and
transformer properties. • Probability of transformer failure due to
overheating. Effects of Transformer Failure on the Population and Electricity
Grid • Regions supplied by transformer nodes. • Electricity Loss for each
region. Final results.

Figure 1: Simplified flowchart of our proposed framework. The general model
consists of three main stages and results in concrete estimates of regional
power outages for arbitrary locations on Earth and how those outages affect
the population. A detailed flowchart is available along with the model’s
source code at https://github.com/allfed/GeomagneticModel.

Geomagnetic storms result from an energy exchange from the solar wind into
Earth’s magnetosphere. They are triggered by sustained periods of high-speed
solar wind, and the most significant storms involve solar coronal mass
ejections (CMEs) [5]. CMEs are the largest-scale eruptive phenomenon in the
solar system – the bulk of plasma with a mass of $10^{13}$ kg is hauled up out
to the interplanetary space with a velocity of more than 1000 km/s [6]. These
storms induce intense currents in Earth’s magnetosphere, affecting the
radiation belts, ionosphere, and thermosphere. We typically measure them by
indices such as the disturbance storm time – Dst – and a planetary geomagnetic
disturbance index – Kp [7]. The consequences include increased atmospheric
density affecting satellite orbits, disruptions to navigation systems like
GPS, and the generation of harmful geomagnetic induced currents (GICs) in
power grids. These extreme events can also be visually stunning in the form of
aurorae, such as the commonly known northern lights, despite their potentially
disruptive nature [8].

GICs are what they sound like – electric currents that appear (are induced)
in, e.g., power lines as a consequence of rapid changes in Earth’s magnetic
field, and those currents can be a cause of worry for the stability and health
of power distribution networks [9].

This worry about GICs mainly consists of two phenomena: overheating components
(predominantly transformers) and tripping safety switches causing cascade
failures [10]. The jury on whether these constitute a significant societal
risk is still out. Scientific studies on this issue are mostly localised
(i.e., consider one particular region instead of all of the Earth), and
severity claims vary from no risk at all [11, 12] to catastrophic, such as
trillions of dollars in damages and years of recovery [13, 14, 15], with most
ending somewhere in between [16, 17, 18, 10, 19, 20]. It is worth noting that
while the global risk of GICs is still undetermined, there are documented
cases of damage to power grids from GICs [19, 20], including severe cases like
the Quebec power failure of $1989$ [21]. The problem with assessing the risk
here is that strong enough events to cause damage are exceedingly rare, and a
replay of the largest documented case from $1859$, commonly referred to as the
Carrington event[22], is yet to be seen.

For further reading on this topic, we refer the Reader to [18, 10, 21] and
references therein.

Since there seems to be a gap in the body of literature dealing with
widespread events reaching global scales, we attempt to partially fill this
gap in this paper. Even though our analysis here shall be limited to Europe
and North America due to limits in data availability, we develop a methodology
to assess the consequences of various strengths of geomagnetic storms on a
global scale. The framework can be applied to any geographical region where
the power grid is known with sufficient accuracy and completeness. The
software implementation of our model is open-
source111https://github.com/allfed/GeomagneticModel, written in Python, and
can be run on a modern PC without the need for a computing cluster or a super-
computer.

In this article, we specifically tackle the issue of overheating transformers
in power grids due to GICs on an international scale, utilising data from
global magnetic field observation (MT) stations [24], electromagnetic transfer
function (EMTF) repository [25, 26, 27], crowd-sourced power distribution
network [28] as well as OpenStreetMap [29], North American Electric
Reliability Corporation (NERC) [30, 31], and Electric Power Research Institute
(EPRI) [32, 33, 34, 35]. In order to assess the effect on population, we used
data from the U.S. Energy Information Administration, and the Center for
International Earth Science Information Network, Columbia University [36, 37,
38].

Our primary target is estimating the population numbers experiencing a
catastrophic electricity loss (CEL) given a recurrence class of a geomagnetic
storm (e.g., a $1$ in $100$ year event).

##  II Methods

In Fig. 1, we present a simplified flowchart of our model (a full picture is
available in the source code repository:
https://github.com/allfed/GeomagneticModel). We partition the model into three
distinct stages, described in detail as follows.

###  II.1 Stage I

This stage is based on the methodology developed by Love et al., [39] and
consists of estimating historical magnetic field – $B$ – levels, determining
apparent conductivity and resistivity at electromagnetic transfer function
sites (EMTF or TF, for short), calculating electric field – $E$ – levels at
magnetotelluric (MT) stations and recurrence rates for those levels, adjusting
$E$ field levels to a reference location, and finally computing $E$ fields
around the globe from the adjusted $E$ fields.

MT stations have been measuring the $B$ field regularly for many decades.
Through these data, we are able to determine the statistics on the recurrence
of $E$ fields over a long period and thus forecast the repeat rates of the
less frequent, larger amplitude events. We obtained the data at the Fresno,
California, USA MT site (FRN) from the global ground-based magnetometer
initiative (SuperMAG) [24]. These data contain measurements of the magnetic
field every minute for $37$ years ending in $2019$. It is worth noting that
the results presented here depend on the choice of the site due to
inaccuracies in the apparent conductivity model. The FRN dataset was picked
because of its proximity to an EMTF site (RET06), sufficient data size, and
time span. For more precise results, data from all existing MT and EMTF sites
would need to be combined into one dataset or a better apparent conductivity
model established, ideally both; however, this falls out of the scope of this
paper.

We then process this dataset as described in [39] Sec. $5$: the time series of
$B$ is detrended by subtracting a second-order polynomial fitted to the entire
duration of the time series, and any missing data are filled by linear
interpolation. Then, the series is put through a fast Fourier transform (FFT),
multiplied by an appropriate transfer tensor for the survey site, and finally,
a reverse FFT is performed, resulting in a geoelectric – $E$ – series. This
procedure is rooted in the relationship between the geomagnetic and
geoelectric fields at the Earth’s surface: [40, 41]

$\mathbf{E}(f,x,y)=\frac{1}{\mu}\mathbf{Z}(f,x,y)\cdot\mathbf{B}(f,x,y),$ (1)

where $f$ is the frequency of sinusoidal variation, $x$ and $y$ are horizontal
Cartesian components, $\mu$ is permeability, and $\mathbf{Z}$ is the impedance
tensor specifying the amplitude, polarisation and phase of the geoelectric
field. The transfer tensor mentioned earlier is $\mathbf{Z}/\mu$.

In order to compute $E$ fields using Eq.(1) we need to determine the apparent
conductivity and resistivity at the corresponding TF site. We acquire the
necessary data from the SPUD EMTF repository [25, 26, 27] at the Dog Creek,
California, USA site (RET06) corresponding to the chosen MT site in Fresno.
The apparent conductivity can be expressed as: [39]

$\sigma_{A}(f)=\frac{2\pi\mu f}{|\mathbf{Z}(f)|^{2}},$ (2)

where $|\cdot|$ is the Frobenius norm[40], and the apparent resistivity is
given by: [42, 40]

$\bm{\rho}(f)=\frac{|\mathbf{Z}(f)|_{\odot}^{2}}{2\pi\mu f},$ (3)

where $|\cdot|_{\odot}$ represents an element-wise norm.

With $E$ fields established, we now consider an adjustment coefficient that
describes the relationship between fields’ strengths and durations (SI Fig.
S1). The peak electric field levels that occur once every ten years change
with the duration over which they are sustained; namely, longer duration
corresponds to lower peak field levels. Thus, we compare the mean field level
lasting sixty seconds to longer-lasting fields. We shall use this adjustment
later to predict the strength of GICs in transformers to ensure we do not
overestimate the current levels of long-duration GICs.

Our ultimate goal of this stage is to compute $E$ fields (their strength and
occurrence rates) at an arbitrary location on Earth, extrapolating from the
measurement site. For this purpose, we adjust all measured fields to an
arbitrarily chosen reference with conductivity
$\sigma_{ref}=8.7*10^{-6}[\frac{\text{S}}{\text{m}}]$, located at a magnetic
latitude $L_{M}=72^{\circ}$ N so that we can model $E$ fields globally with a
single function regardless of the choice of the MT site or even ensemble of
sites. To this end, we incorporate an auroral boundary shift using an
estimated auroral boundary movement of $0.004$ degrees threshold shift south
per [nT] Dst as per the recommendation in the report from EPRI [32], with the
Dst estimates taken from Oughton et al., [18]. Subsequently, we use the
geomagnetic latitude distribution data from Ngwira et al.,[43] to compute a
magnetic latitude adjustment $\hat{L_{M}}$, and we establish the apparent
conductivity adjustment as $\hat{\sigma_{A}}=\sigma_{A}/\sigma_{A}^{ref}$.

Thus, we obtain a formula for the reference $E$ field:

$|E_{ref}|=\hat{L_{M}}\hat{\sigma_{A}}|E|.$ (4)

We can now fit a log-normal function to all $E$ field data (SI Fig. S2) and
then use this model to estimate geoelectric fields at an arbitrary point on
Earth [44, 39].

This arbitrary location $E$ field estimation is conducted with the help of
Earth’s ground conductivity at a $15$ arc-minute basis from Alekseev et al.,
[45]. It is a layered model, so we calculate the impedance recursively as
prescribed in the NERC report [31]. Finally, we adjust the E field at the
reference site back to the expected field levels for all points on Earth at
$15$ arc-minute resolution via:

$|E|=\frac{1}{\hat{L_{M}}\hat{\sigma_{A}}}|E_{ref}|.$ (5)

This concludes stage one, and in Fig. 2, we present our prediction for the
expected peak magnitude of a sixty-second $E$ field in a $1$-in-$100$ year
storm for Europe and North America.

Figure 2: Peak $60$-second geoelectric field magnitude for a $1$-in-$100$ year
storm in Europe (top) and North America (bottom). The colour indicates the
field level [V/km]. Figure 3: European substations at risk of electricity loss
after a $1$-in-$10000$ year storm. Colour indicates the fraction of HVTs we
expect to fail within a given Voronoi tile.

###  II.2 Stage II

In this stage, we determine the GICs in power grids and, thus, the transformer
failure due to overheating. We gathered the worldwide power distribution
network data from OpenGridMap [28]. It is an open-source initiative that
utilises crowd-sourced data to automatically generate power grid models in a
Common Information Model (CIM) format, which is the power industry standard
[46]. Since we focus in this paper on high voltage networks and high voltage
transformers (HVT), we selected all lines above $100$ [kV]. As pointed out in
the introduction, we will be restricting ourselves to only two continents, and
the primary reason for this is data availability (or lack thereof) in the
OpenGridMap. We present plots of the selected networks in SI Fig. S3.

On these grids, the GICs were computed with the help of the open-source
package – GEOMAGICA [47], which implements the Lehtinen-Pirjola method [48]
for the Austrian grid. We modified this code such that it can be applied to an
arbitrary location on Earth. Each power line length was preserved accurately,
and the full length of the line was used to estimate overall line resistance;
however, the GEOMAGICA package only allows for straight line estimates of
induced GIC. The resistance per kilometre, as well as HVT winding conductance,
were estimated using linear interpolation of line conductance as a function of
AC voltage using values provided in a NERC report [30], with the ground
resistance being set to $0.5[\Omega]$.

Final GIC per phase values for transformers were estimated as the expected GIC
per phase flowing through the ground grid divided by $1.3$. This factor brings
the final values closer to those obtained with more complex transformer
circuit configurations and is discussed in further detail in [18]. The
absolute values of geoelectric fields were applied uniformly in the East-West
direction, as this is the primary field orientation [49, 50]. If the GIC
duration is longer than a minute, we scale its magnitude down to account for
the decreasing statistical likelihood of an $E$ field sustained for a longer
duration (as explained in stage one, also see, SI Fig. S1), and the fact that
GICs are proportional to the electric field $I_{GI}\propto E$.

With GIC levels and their recurrence rates established, we move on to
determine whether they cause sufficient heat increase in transformers to
result in transformer failure. We source the steady state temperature rise
$T^{*}(I_{GI})$ for $42$ different types of transformers and two tie bar
geometries from an EPRI report [34] and linearly interpolate the data where
necessary, as well as the temperature response for both designs – $\tau$.
These data were extracted by hand from the said report. Both $T^{*}$ and
$\tau$ are multiplied by $\frac{2}{3}$ to accommodate the fact that the
studied tie bar geometry designs represent the extreme “best” and “worst”
cases; thus, we are considering the midpoint between them (this scaling factor
can be found via inspecting the EPRI report). The HVT temperature can be found
as: [33]

$T_{HVT}=T^{*}(I_{GI})\left(1-e^{-\frac{t}{\tau}}\right)+T_{oil},$ (6)

where $t$ is the duration, and $T_{oil}=90^{\circ}C$ is the surface oil
temperature assumed in a severe geomagnetic storm [34].

To determine whether the temperature in Eq. (6) makes a transformer fail, we
look at the structural temperature limit in [33], which is a function of the
transformer’s age. We estimate the transformer population age distribution
based on the statistics in the same report. Each transformer has an associated
voltage class ($V_{c}$), which we assign based on the proximity to one of
$\\{230,345,500,765\\}$ [kV]. These values are chosen specifically as these
have associated statistics on the transformer phase demographics in the
aforementioned report and prevalence in the US grid per class [35]. The US
dataset was the easiest to obtain, and we assume grids around the world would
not have significantly different characteristics in these particular features.
However, should that not be the case, it is fairly straightforward to
substitute this number for a region-specific one if needed (for those who have
access to the data). We use the cumulative fraction of transformer population
exceeding safe temperatures as the proxy for the probability of failure for
the transformer in a given class; therefore, we can write:

$F(V_{c})=f(V_{c})\sum_{A}f(T_{HVT}>T_{max})f(A),$ (7)

where $F(V_{c})$ is the fraction of transformers overheating in class $V_{c}$,
$A$ is the age category, with $f(A)$ the fraction of transformers in given age
category, and $f(T_{HVT}>T_{max})$ is the fraction of transformers above safe
temperature for the age, and $f(V_{c})$ is the fraction associated with the
transformer type. For each of the $42$ transformer types, we associate $V_{c}$
and a phase $\phi$ with them, and now, each transformer is a member of a
group. E.g., there are seven types in a group defined by
$\\{V_{c}=500,\phi=1\\}$, so for each transformer in this group,
$f(V_{c})=1/7$. All the values for transformer types, temperature limits, and
transformer population are presented in SI Table S1, S2, and S3.

Ultimately, we want to assess whether a power network node $n$ fails, and we
achieve that by estimating the percentage of transformers failing at any given
node as:

$F(n)=\sum_{\\{\phi,V_{c}\\}}\left(f(\phi)f_{grid}(V_{c})\sum_{V_{c}}F(V_{c})\right).$
(8)

$f_{grid}(V_{c})$ is the prevalence of transformers in class $V_{c}$ in the
power grid, and $f(\phi)$ is the fraction of transformers with phase $\phi$.
If $F(n)>0.33$, then we conclude that node $n$ fails longer term. This cutoff
value is based on the report for Oak Ridge National Laboratory[51], in which
we can find that “the standard approach to spares has been to purchase an
extra single phase transformer for a three phase bank”. Thus, we assume that
if the failure rate is more than $33\%$ in a given region, the region suffers
an extended power outage.

###  II.3 Stage III

We conclude how transformer failures affect population access to electricity.
Our consideration here is based on the tessellation of the continents. It
could also be performed in countries and alike, but we focus on continents in
the analysis. Unfortunately, the region served by each substation node is not
publicly available, so we use a proxy approach – a Voronoi diagram. Each power
grid node in the power distribution data serves as the centre for a Voronoi
cell. This way, we have a good approximation of real-world electric districts
since each power substation serves its nearest surroundings. Naturally, we
limit the Voronoi regions to the extent of the geographic boundaries of the
networks. We use the $15$ arc-minute grid population density data to estimate
the population in each cell [37]. Specifically, in each grid cell, the land
area was multiplied by the population density to obtain the absolute
population estimates. The total population in each Voronoi region was
calculated as the total number of centres of $15$ arc-minute population grid
cells within the Voronoi polygon represented by a power node.

We assume the power outages to be long-term as the grid is often run at near-
maximum capacity, and it is estimated that a region with insufficient
transformers would be unable to serve electricity for many months until those
transformers were replaced, which is a slow and expensive process. This is a
lower bound estimate, however, and could well take much longer [51].

In addition to the number of people being cut off from power, we provide
estimates on the total electric power loss in [GW]:

$P_{loss}=\frac{P_{cons}}{N}\sum_{n}C(n)[F(n)>0.33].$ (9)

$P_{loss}$ is power lost, $P_{cons}$ is total power consumption (for
continent/country, etc.) [38], $N$ is total population [36], $C(n)$ is the
population in a given Voronoi region $n$, $F(n)$ is from Eq. (8), and
$[\dots]$ is the Iverson bracket notation.

##  III Results and Discussion

Figure 4: North American substations at risk of electricity loss after a
$1$-in-$10000$ year storm. Colour indicates the fraction of HVTs we expect to
fail within a given Voronoi tile.

To recap, In Fig. 1, we presented a simplified flowchart of our model. The
data-driven framework we propose runs in three stages estimating: 1)
geoelectric fields for severe storms, 2) GICs and transformer overheating in
high voltage power distribution networks, and 3) effects of transformer
failure on the electric grid and population’s access to electricity. This
three-stage process is an agglomeration of multiple smaller models we
assembled across various scientific fields and modified or expanded where
appropriate to suit our needs.

We utilise magnetotelluric (MT) [24] and EMTF [25, 26, 27] site data to find
geoelectric fields as described by Love et al., [39]. Then, an adjustment
process is applied to all fields to a reference location, with an intermediary
step incorporating auroral boundary movement. The reference adjustment is done
by using the geomagnetic latitude data from Ngwira et al., [43]. To this data,
we fit a log-normal function that serves as our model for the recurrence of
the geoelectric fields [44, 39]. Finally, an impedance estimate is computed
using Earth’s ground conductivity model from Alekseev et al., [45]. With all
those elements, we can compute $E$ field levels at any location on Earth for
any storm rarity and duration.

This general model can be used for an arbitrary geographical region, including
the whole globe, as long as power distribution network data are available. The
software implementation is written in Python and the source code is open and
available for anyone to use and modify. It is also worth noting that our
implementation is additionally accessible to a wide audience as it does not
require an expensive computer cluster, and can be run on a reasonably modern
machine in less than an hour on the scale of a continent or smaller.

In order to see how particular regions are affected by geomagnetic storms, we
include the Voronoi diagram plots for both continents depicting a specific
number of substations at risk in a granular manner; see Fig. 3 for Europe, and
Fig. 4 for North America in a $1$-in-$10000$ year scenario. The colour
indicates the ratio of HVTs within a given Voronoi cell that are expected to
fail. We acknowledge that Sweden has in fact hardened their grid to address
the risks discussed here [52]; however, for simplicity’s sake, we include it
in our results.

Figure 5: Population volume experiencing CEL as a function of storm rarity in
$[\%]$ (left), the number of people (centre), and electricity loss in
gigawatts (right).

We determine a region to be expected to suffer a power loss by considering the
transformer failures exceeding available spare parts [51]. Thus, we assume
that if the failure rate is more than $33\%$ in a given region, the region
suffers an extended power outage. In Fig. 5, we present the estimates of the
population numbers experiencing a CEL as a result of a geomagnetic storm, and
total electricity loss in gigawatts. This plot illustrates how many people
will lose power for months or perhaps even years when an extreme storm
arrives. This extreme event’s strength is the highest expected to occur once
per certain amount of years. For instance, for an event of a strength that is
expected to happen once in $10000$ years, in Europe, we expect $0.99\%$ (over
$7.4$ million people) of the population to experience CEL, whereas in North
America, it would be $0.86\%$ (over $4.1$ million people). The total expected
loss is $49.71$ [GW] and $56.12$ [GW], respectively. We limited ourselves here
to these particular metrics, but of course, this analysis can be expanded to
include the expected loss of GDP, industry production, and so on.

Geomagnetically induced currents are a serious cause for concern, considering
how reliant many of humanity’s critical infrastructure systems are on electric
power [53]. While plenty of papers tackle this issue in a localised and narrow
manner, to our knowledge, no attempt at a global risk analysis has been
conducted. We tackle the problem of filling this knowledge gap and develop a
general model that combines all adequate specific methods from thermal models
of HVT tie bar hot spots through historical geoelectric field estimates to a
global conductivity model. It is our belief that the methodology presented in
this paper can be applied on a global scale and inform the decisions of
policymakers in regulating power distribution network resilience.

Our main analysis was limited to only two continents – Europe and North
America – due to a lack of quality data in other regions; however, these two
serve as very good proofs of concept. The key results show that a severe
geomagnetic storm can severely affect millions of people. We believe that
rather than approaching our results as final, they should be considered as a
starting point for further analysis and a powerful tool and framework for
future work, where more accurate grid connections, transformer designs and
voltages, grounding resistances, and other specifics on a global scale can be
used to obtain realistic results for a given region. The unpredictable and
rare nature of coronal mass ejections impacting Earth makes the risk
assessment much more difficult since, at the time of writing this manuscript,
it is virtually impossible to say what precisely the chances of a
$1$-in-$10000$ year storm occurring within the next decade are. While it seems
prudent and most likely cost-effective to prepare the grid for severe storms
instead of waiting for failures to happen, a global-scale complete power
outage does not seem possible with the assumptions we have made in our model.

With that said, the data we have used to predict HVT overheating is limited.
Our grid model relies on potentially incomplete or incorrect infrastructure
data crowd-sourced by OpenStreetMaps, due to the lack of official public high
voltage infrastructure data. However, this does not represent a failure of the
model, as electrical utilities and government agencies can readily replace our
grid map with their own private data for enhanced prediction accuracy.

One could potentially adapt our model to use a characteristic power line
length or a synthetic network instead of real data, but this leads to the
problem of losing regional precision in predicted transformer outages. We
would still have a global estimate of expected power loss, but the exact
locations of predicted outages would be unreliable.

Similarly, the magnetotelluric data are still very limited, and while our
method of normalising these data to a reference location and, from there,
extrapolating worldwide has the advantage of universality, it does come with a
drawback of potentially over-generalisation. This can lead to underestimating
storm recurrence in some specific regions.

Most importantly, we only consider transformer failure due to overheating;
however, the safety switch cascade failure scenario is also a very important
avenue of investigation that should be explored in the future, as well as
other more complex interactions with the increasingly complex digital
infrastructure automatically routing power through the electrical grid.

Overall, to our knowledge, the approach proposed in this paper is the first of
its kind in assessing power distribution network resilience to GICs and can be
broadly applied at municipality, country, continent or even a global level to
inform governing bodies and private entities alike of the dangers associated
with transformer overheating during geomagnetic storms.

## Disclosure statement

Authors report that there are no competing interests to declare.

## Data availability statement

All data used in this research are publicly available. The program source code
and instructions on how to procure the data can be found in the github
repository of the project: https://github.com/allfed/GeomagneticModel.

## References

  * Lakhina and Tsurutani [2016] G. S. Lakhina and B. T. Tsurutani, Geomagnetic storms: historical perspective to modern view, Geoscience Letters 3, 1 (2016).
  * Organek [2017] J. F. Organek, Black sky hazards: Resilience planning, in  _Congress on Technical Advancement 2017_ (2017) pp. 29–40.
  * Schulz [2024] K. Schulz, What a Major Solar Storm Could Do to Our Planet, https://www.newyorker.com/magazine/2024/03/04/what-a-major-solar-storm-could-do-to-our-planet (2024), accessed: 2024-03-05.
  * Löfquist [2020] L. Löfquist, Is there a universal human right to electricity?, The International Journal of Human Rights 24, 711 (2020).
  * Gosling [1993] J. T. Gosling, The solar flare myth, Journal of Geophysical Research: Space Physics 98, 18937 (1993).
  * Chen [2011] P. Chen, Coronal mass ejections: models and their observational basis, Living Reviews in Solar Physics 8, 1 (2011).
  * Gonzalez  _et al._ [1994] W. Gonzalez, J.-A. Joselyn, Y. Kamide, H. W. Kroehl, G. Rostoker, B. T. Tsurutani, and V. Vasyliunas, What is a geomagnetic storm?, Journal of Geophysical Research: Space Physics 99, 5771 (1994).
  * noa [2014] Geomagnetic Storms, https://www.swpc.noaa.gov/phenomena/geomagnetic-storms (2014), accessed: 2023-11-24.
  * Pulkkinen  _et al._ [2017] A. Pulkkinen, E. Bernabeu, A. Thomson, A. Viljanen, R. Pirjola, D. Boteler, J. Eichner, P. Cilliers, D. Welling, N. Savani, _et al._ , Geomagnetically induced currents: Science, engineering, and applications readiness, Space weather 15, 828 (2017).
  * Roodman [2015] D. Roodman, The risk of geomagnetic storms to the grid: A preliminary review, Available at SSRN 3635861 (2015).
  * Švanda  _et al._ [2021] M. Švanda, A. Smičková, and T. Vỳbošt’oková, Modelling of geomagnetically induced currents in the Czech transmission grid, Earth, Planets and Space 73, 1 (2021).
  * Girgis and Vedante [2012] R. Girgis and K. Vedante, Effects of GIC on power transformers and power systems, in  _PES T &D 2012_ (IEEE, 2012) pp. 1–8.
  * National Research Council [2008] National Research Council,  _Severe space weather events: understanding societal and economic impacts: a workshop report_ (The National Academies Press, Washington, DC, 2008) https://nap.nationalacademies.org/catalog/12507/severe-space-weather-events-understanding-societal-and-economic-impacts-a.
  * Said [2020] C. Said, Everything I’ve learned about solar storm risk and EMP attacks, https://chris-said.io/2020/06/18/everything-ive-learned-about-solar-storm-risk-and-emp-attacks/ (2020), accessed: 2023-11-24.
  * Weiss and Weiss [2019] M. Weiss and M. Weiss, An assessment of threats to the American power grid, Energy, Sustainability and Society 9, 1 (2019).
  * Torta  _et al._ [2014] J. M. Torta, S. Marsal, and M. Quintana, Assessing the hazard from geomagnetically induced currents to the entire high-voltage power network in Spain, Earth, Planets and Space 66, 1 (2014).
  * Gil  _et al._ [2021] A. Gil, M. Berendt-Marchel, R. Modzelewska, S. Moskwa, A. Siluszyk, M. Siluszyk, L. Tomasik, A. Wawrzaszek, and A. Wawrzynczak, Evaluating the relationship between strong geomagnetic storms and electric grid failures in Poland using the geoelectric field as a GIC proxy, Journal of Space Weather and Space Climate 11, 30 (2021).
  * Oughton  _et al._ [2019] E. J. Oughton, M. Hapgood, G. S. Richardson, C. D. Beggan, A. W. P. Thomson, M. Gibbs, C. Burnett, C. T. Gaunt, M. Trichas, R. Dada, and R. B. Horne, A Risk Assessment Framework for the Socioeconomic Impacts of Electricity Transmission Infrastructure Failure Due to Space Weather: An Application to the United Kingdom, Risk Analysis 39, 1022 (2019), https://onlinelibrary.wiley.com/doi/pdf/10.1111/risa.13229 .
  * Gaunt and Coetzee [2007] C. Gaunt and G. Coetzee, Transformer failures in regions incorrectly considered to have low GIC-risk, in  _2007 IEEE Lausanne power tech_ (IEEE, 2007) pp. 807–812.
  * Moodley and Gaunt [2012] N. Moodley and C. Gaunt, Developing a power transformer low energy degradation assessment triangle, in  _IEEE Power and Energy Society Conference and Exposition in Africa: Intelligent Grid Integration of Renewable Energy Resources (PowerAfrica)_ (IEEE, 2012) pp. 1–6.
  * Gaunt [2016] C. Gaunt, Why space weather is relevant to electrical power systems, Space Weather 14, 2 (2016).
  * Tsurutani  _et al._ [2003] B. T. Tsurutani, W. D. Gonzalez, G. Lakhina, and S. Alex, The extreme magnetic storm of 1–2 September 1859, Journal of Geophysical Research: Space Physics 108 (2003).
  * Note [1] https://github.com/allfed/GeomagneticModel.
  * Gjerloev [2009] J. Gjerloev, A global ground-based magnetometer initiative, Eos, Transactions American Geophysical Union 90, 230 (2009), https://supermag.jhuapl.edu/. Accessed: 2023-12-11.
  * Kelbert  _et al._ [2011] A. Kelbert, G. Egbert, and A. Schultz, IRIS DMC data services products: EMTF, the magnetotelluric transfer functions, National Geoelectromagnetic Facility Technical Report (2011).
  * Kelbert [2019] A. Kelbert, Taking magnetotelluric data out of the drawer, in  _AGU Fall Meeting Abstracts_ , Vol. 2019 (2019) pp. IN21A–01, https://ds.iris.edu/spud/emtf. Accessed: 2023-12-11.
  * Kelbert  _et al._ [2019] A. Kelbert, P. A. Bedrosian, and B. S. Murphy, The first 3D conductivity model of the contiguous United States: Reflections on geologic structure and application to induction hazards, Geomagnetically induced currents from the Sun to the power grid , 127 (2019).
  * Rivera  _et al._ [2017] J. Rivera, J. Leimhofer, and H.-A. Jacobsen, OpenGridMap: towards automatic power grid simulation model generation from crowdsourced data, Computer Science-Research and Development 32, 13 (2017), https://github.com/OpenGridMap/transnet-models. Accessed: 2023-12-11.
  * OpenStreetMap contributors [2017] OpenStreetMap contributors, Region extent data retrieved from https://download.geofabrik.de/, https://www.openstreetmap.org (2017), accessed: 2023-12-13.
  * North American Electric Reliability Corporation (2013a) [NERC] North American Electric Reliability Corporation (NERC), Network applicability, project 2013-03, EOP-010-1, https://www.nerc.com/pa/Stand/Project201303GeomagneticDisturbanceMitigation/ApplicableNetwork.pdf (2013a), accessed: 2023-12-11.
  * North American Electric Reliability Corporation (2013b) [NERC] North American Electric Reliability Corporation (NERC), Computing geomagnetically-induced current in the bulk-power system, https://www.nerc.com/comm/PC/Geomagnetic%20Disturbance%20Task%20Force%20GMDTF%202013/GIC%20Application%20Guide%202013_approved.pdf (2013b), accessed: 2023-12-12.
  * Electric Power Research Institute , Palo Alto, CA: 2020. 3002017952.(2020) [EPRI] Electric Power Research Institute (EPRI), Palo Alto, CA: 2020. 3002017952., Magnetohydrodynamic (MHD) modeling for the further understanding of geoelectric field enhancements and auroral behavior during geomagnetic disturbance events, https://www.epri.com/research/products/000000003002017952 (2020), accessed: 2023-12-12.
  * Electric Power Research Institute , EPRI, Palo Alto, CA: 2017. 3002009001(2017) [EPRI] Electric Power Research Institute (EPRI), EPRI, Palo Alto, CA: 2017. 3002009001, Magnetohydrodynamic electromagnetic pulse assessment of the continental U.S. electric grid: Geomagnetically induced current and transformer thermal analysis, https://ed-thelen.org/pics/EPRI-HEMP-2017.pdf (2017), accessed: 2023-12-12.
  * Electric Power Research Institute , Palo Alto, CA: 2019. 3002017708.(2019) [EPRI] Electric Power Research Institute (EPRI), Palo Alto, CA: 2019. 3002017708., Transformer thermal impact assessments for dc withstand capability: Examining the impacts of geomagnetically induced current (GIC) on transformer thermal performance, https://www.epri.com/research/products/000000003002017708 (2019), accessed: 2023-12-12.
  * Kappenman and Mcconnell [2010] J. G. Kappenman and B. W. Mcconnell, Meta-R-322 Low-Frequency Protection Concepts for the Electric Power Grid: Geomagnetically Induced Current (GIC) and E 3 HEMP Mitigation, Metatech Corporation, 358 S. Fairview Ave., Suite E, Goleta, CA 93117 https://api.semanticscholar.org/CorpusID:62781689 (2010).
  * U.S. Energy Information Administration [2020a] U.S. Energy Information Administration, Population by country, https://www.eia.gov/international/data/world/other-statistics/population-by-country (2020a), accessed: 2023-12-13.
  * Center for International Earth Science Information Network - CIESIN - Columbia University [2020] Center for International Earth Science Information Network - CIESIN - Columbia University, Gridded Population of the World, Version 4 (GPWv4): Population Density, Palisades, New York: NASA Socioeconomic Data and Applications Center (SEDAC). (2020), https://doi.org/10.7927/H49C6VHW. Accessed: 2023-12-13.
  * U.S. Energy Information Administration [2020b] U.S. Energy Information Administration, Electricity consumption by country, https://www.eia.gov/international/data/world/electricity/electricity-consumption (2020b), accessed: 2023-12-13.
  * Love  _et al._ [2018] J. J. Love, G. M. Lucas, A. Kelbert, and P. A. Bedrosian, Geoelectric hazard maps for the Mid-Atlantic United States: 100 year extreme values and the 1989 magnetic storm, Geophysical Research Letters 45, 5 (2018).
  * Berdichevsky and Dmitriev [2008] M. N. Berdichevsky and V. I. Dmitriev,  _Models and methods of magnetotellurics_ (Springer Science & Business Media, 2008).
  * Chave and Jones [2012] A. D. Chave and A. G. Jones,  _The magnetotelluric method: Theory and practice_ (Cambridge University Press, 2012).
  * Love  _et al._ [2019] J. J. Love, G. M. Lucas, P. A. Bedrosian, and A. Kelbert, Extreme-value geoelectric amplitude and polarization across the Northeast United States, Space weather 17, 379 (2019).
  * Ngwira  _et al._ [2013] C. M. Ngwira, A. Pulkkinen, F. D. Wilder, and G. Crowley, Extended study of extreme geoelectric field event scenarios for geomagnetically induced current applications, Space Weather 11, 121 (2013).
  * Love [2020] J. J. Love, Some experiments in extreme-value statistical modeling of magnetic superstorm intensities, Space Weather 18, e2019SW002255 (2020).
  * Alekseev  _et al._ [2015] D. Alekseev, A. Kuvshinov, and N. Palshin, Compilation of 3D global conductivity model of the Earth for space weather applications, Earth, Planets and Space 67, 1 (2015), https://globalconductivity.ocean.ru/. Accessed: 2023-12-11.
  * Simmins [2011] J. J. Simmins, The impact of PAP 8 on the Common Information Model (CIM), in  _2011 IEEE/PES Power Systems Conference and Exposition_ (IEEE, 2011) pp. 1–2.
  * Weidelt, P. and Beggan, C. and Turnbull, A. and Mckay, K. and Bailey, R. [2018] Weidelt, P. and Beggan, C. and Turnbull, A. and Mckay, K. and Bailey, R., Geomagnetically induced currents in Austria (GEOMAGICA), https://github.com/geomagpy/GEOMAGICA (2018), accessed: 2023-12-11.
  * Lehtinen and Pirjola [1985] M. Lehtinen and R. J. Pirjola, Currents produced in earthed conductor networks by geomagnetically-induced electric fields, Annales Geophysicae 3, 479 (1985).
  * Kikuchi  _et al._ [2008] T. Kikuchi, K. K. Hashimoto, and K. Nozaki, Penetration of magnetospheric electric fields to the equator during a geomagnetic storm, Journal of Geophysical Research: Space Physics 113 (2008).
  * Marshall  _et al._ [2013] R. Marshall, H. Gorniak, T. Van Der Walt, C. Waters, M. Sciffer, M. Miller, M. Dalzell, T. Daly, G. Pouferis, G. Hesse, _et al._ , Observations of geomagnetically induced currents in the Australian power network, Space Weather 11, 6 (2013).
  * Kappenman [2010] J. Kappenman, Meta-R-319 Geomagnetic storms and their impacts on the US power grid, Metatech Corporation, 358 S. Fairview Ave., Suite E, Goleta, CA 93117 https://www.ferc.gov/sites/default/files/2020-05/ferc_meta-r-319.pdf (2010).
  * E. Krausmann, E. Andersson , W. Murtagh and N. Mitchison [2013] E. Krausmann, E. Andersson , W. Murtagh and N. Mitchison,  _Space Weather and Power Grids: Findings and Outlook. An event co-organised by the European Commission’s Joint Research Centre, the Swedish Civil Contingencies Agency and the NOAA Space Weather Prediction Centre 29-30 October, 2013, Ispra, Italy_ , Tech. Rep. LB-NA-26370-EN-N (Luxembourg, 2013) https://publications.jrc.ec.europa.eu/repository/handle/JRC86658.
  * Foster  _et al._ [2008] J. S. Foster, E. Gjelde, W. R. Graham, R. J. Hermann, H. M. Kluepfel, R. L. Lawson, G. K. Soper, L. L. Wood, and J. B. Woodard, Report of the commission to assess the threat to the United States from electromagnetic pulse (EMP) attack, critical national infrastructures, Critical National Infrastructures (2008), https://apps.dtic.mil/sti/pdfs/ADA484672.pdf.