The Role of Global/Regional Earth Conductivity Models in Natural Geomagnetic Hazard Mitigation

Abstract

Geomagnetic disturbances cause perturbations in the Earth’s magnetic field which, by the principle of electromagnetic induction, in turn cause electric currents to flow in the Earth. These geomagnetically induced currents (GICs) also enter man-made technological conductors that are grounded; notably, telegraph systems, submarine cables and pipelines, and, perhaps most significantly, electric power grids, where transformer groundings at power grid substations serve as entry points for GICs. The strength of the GICs that flow through a transformer depends on multiple factors, including the spatiotemporal signature of the geomagnetic disturbance, the geometry and specifications of the power grid, and the electrical conductivity structure of the Earth’s subsurface. Strong GICs are hazardous to power grids and other infrastructure; for example, they can severely damage transformers and thereby cause extensive blackouts. Extreme space weather is therefore hazardous to man-made technologies. The phenomena of extreme geomagnetic disturbances, including storms and substorms, and their effects on human activity are commonly referred to as geomagnetic hazards. Here, we provide a review of relevant GIC studies from around the world and describe their common and unique features, while focusing especially on the effects that the Earth’s electrical conductivity has on the GICs flowing in the electric power grids.

1 Introduction

The phenomena of space weather (e.g., Song et al. 2001; Lanzerotti 2001; Moldwin 2008; Love et al. 2014; Fujita et al. 2016) have been observed to cause disruption in man-made technological systems for well over a century. The first paper on the impact of a geomagnetic disturbance on technological infrastructure, known to the author, appeared in 1849 in relation to space weather effects on telegraph systems (Barlow 1849). In fact, as Boteler (2003) points out, Bain’s “chemical telegraph” used specially prepared paper; current from a stylus caused a chemical reaction leaving a colored mark on the paper. During the geomagnetic storm of February 19, 1852, the current increased so much that a “flame of fire” followed the pen and set fire to the paper (Prescott 1866, p. 318). A series of similar reports followed, of which the most dramatic occurrences were associated with the geomagnetic storm now known as the Carrington event, observed on September 1–2, 1859 (Carrington 1859), that set telegraph wires ablaze (Loomis 1861; Boteler 2006; Shea and Smart 2006). Throughout the twentieth century, expanding phone and radio communication networks were also repeatedly disturbed by space weather events (e.g., Boteler 2003).

Apart from loss of communications, another long-standing risk to infrastructure from geomagnetic disturbances is pipeline corrosion, which can result due to a cumulative effect of many space weather events driving the pipeline potentials outside of their optimum voltage range (e.g., Peabody 1979; Gummow and Eng 2002), and thus accelerating corrosion. Finally, the second half of the twentieth century has seen space weather effects interfering with the operations of electric power transmission systems. Severe geomagnetic disturbances can permanently damage transformers, which are difficult and expensive to replace; they can also cause blackouts. The March 13, 1989, geomagnetic storm famously caused a complete collapse of the Hydro-Québec power system in Canada, leaving the 6 million residents of Québec without power for over 9 h (e.g., Allen et al. 1989; Boteler 2001; Bolduc 2002). It is less widely known that the voltage oscillations caused tripping of protective equipment, nearly bringing most of the northeastern USA down into a cascading collapse (e.g., Homeier and Wei 2013; see also North American Electric Reliability Corporation 1990 for a compilation of effects). To list a few failure reports from the early 2000s, the geomagnetic storm of November 6, 2001, caused the failure of electric power transmission systems in New Zealand (Béland and Small 2005), and the Halloween geomagnetic storm of October 29–31, 2003, caused operational failures in Swedish power grid systems (Pulkkinen 2005) (see Eastwood et al. 2017 for a list of such examples all over the world). Railway automatic system failures have also been observed in Russia during geomagnetic storms (Eroshenko et al. 2010 and references therein).

Indeed, the hazards of space weather to the modern technological infrastructure are varied (e.g., Lanzerotti 2001; Cannon et al. 2013; Schrijver 2015), ranging from the impacts of ionizing radiation on avionics and ground systems, to radiation impacts on satellites, to impacting satellite communication and the global navigation satellite system (GNSS) positioning, to power grid failures, with potentially enormous societal impact and economic ramifications (e.g., Baker et al. 2008; Eastwood et al. 2017). Because of our increasing reliance on electrical power, communications, and navigation systems in our daily lives, a modern day repeat of an extreme, Carrington-class geomagnetic storm (such as, for example, the 2012 solar eruptive event that, thankfully, missed the Earth; Ngwira et al. 2013) could result in significant damage and disruption. We should also note that probability estimates for a Carrington-class storm show a real possibility of recurrence (Riley 2012; Kataoka 2013; Love et al. 2015). In fact, a compilation of evidence from solar, stellar, terrestrial, and lunar observations and residual proxies led Schrijver and Beer (2014) to conclude that very much more energetic events could occur on the Sun, albeit with a significantly lower probability of occurrence. The possibility of a solar storm of Carrington-class or higher intensity hitting the Earth has caused the governments of countries as geographically distributed as the USA, UK, China, Australia, and others to develop national strategies and risk assessment procedures with respect to space weather hazard mitigation.

For the purposes of this review, we are specifically focused on geomagnetic hazards to infrastructure located at or near the Earth’s surface. Through the phenomenon of electromagnetic induction, perturbations of the magnetosphere and the related ionospheric currents caused by a space weather event induce electric currents to flow in the conducting Earth. These currents then enter technological conductors through grounding points, causing damage. The term geomagnetically induced currents (GICs) most commonly refers to the quasi-DC electric current flowing in man-made technological conductors during geomagnetic disturbances. The most direct measurement proxy for GICs is the neutral to ground current, measured at the grounding points of man-made conductors. The modern era of investigation of geomagnetic hazards came in the 1970s and 1980s, when the effects of geomagnetic variations on pipelines and power grids were first recorded and formally analyzed. GICs were first recorded and considered analytically on the trans-Alaska pipeline (Campbell 1978) and in a Bavarian pipeline (Brasse and Junge 1984). Rigorous GIC modeling techniques for pipelines were developed by Pulkkinen et al. (2001). Initial studies on GICs in the power grids appeared at about the same time (e.g., Albertson and Van Baelen 1970; Anderson III et al. 1974; Lanzerotti 1983; Pirjola and Lehtinen 1985), and the analytical groundwork for modeling the GICs in power lines was laid out in the 1980s (Lehtinen and Pirjola 1985). Excellent GIC monitoring conducted in Finland since 1977 allowed Viljanen and Pirjola (1989) and related studies, for the first time, to develop and validate complete power grid GIC modeling procedures. Slight variations on these procedures are used to this date by many of the research groups, and power grid system operators, all over the world.

In general, the Earth’s electrical conductivity plays a critical role in the amplitudes of GICs that enter man-made conductors. However, GICs in grounded infrastructure are distinctly different depending on whether the infrastructure is continuously grounded (such as the pipelines, submarine cables, or railways), or discretely grounded such as an electric power grid. These two types of infrastructure call for different analysis methods. While few analyses of GICs in pipelines attempt a quantitative model to data comparison (most studies, e.g., Hejda and Bochníček 2005 are satisfied with observed correlations between quantities), the few relevant studies that we could find, most notably Pulkkinen et al. (2001) and Viljanen et al. (2006), use a very simple, 1D depth-dependent Earth conductivity for GIC modeling in pipelines with reasonable success. For example, Pulkkinen et al. (2001) GIC estimates, when compared against measurements in a Finnish pipeline, exhibit high correlation, while the peak amplitudes are only of some 30–60%; this can be viewed as a success considering the Earth was assumed uniformly conductive. Viljanen et al. (2006) show a nearly perfect match to data with a simple two-layer Earth model. This relative insensitivity of pipeline GICs to the true Earth’s structure could be explained by the lengths of typical pipelines, combined with the nature of their grounding. We find, however, that the influence of the Earth electrical conductivity model on GIC estimates in discretely grounded conductors such as the power lines is rather dramatic. We shall therefore focus primarily on the electric power grid studies for this review.

A comprehensive review of the GIC estimation in Finland by Viljanen and Pirjola (1994) has outlined the central aspects, and a sequence of analysis steps necessary for a GIC estimation study. In this and prior works, the authors have been very careful to delineate the practical aspects of the modeling, which were then dramatically limited by the computing power and data availability, from the physics of the process. The authors have accurately proposed that the geophysical phenomena responsible for GICs are very complex and include both the ionospheric and magnetospheric current systems as they are affected by the solar wind and the Earth’s electrical conductivity structure. At the time, however, only the simplest models of the ionospheric–magnetospheric current system were tractable and were therefore explored; namely, the source structure due to auroral current systems was modeled as a localized electrojet, while the current systems away from the auroral regions were modeled as a homogeneous plane wave. Oftentimes, an effective model involved a combination of these sources, with some free parameters to fit to geomagnetic field measurements. Similarly, on the Earth’s conductivity side, the most sophisticated conductivity model that was considered involved several one-dimensional (1D) profiles and an explicit assumption that their finite lateral extent does not affect the 1D approximation and the corresponding Earth’s impedance. The overall system allowed the authors to fit suitable parameter values to measured GIC data.

Both the computational capabilities and data availability have dramatically improved in the past several decades, but the overall layout still stands strong. Fueled by wide availability of modern computing capabilities, a wide range of GIC estimation techniques and applications have flourished. Some researchers turned their attention to improving the power system modeling, and others to enhanced electric field modeling capabilities. Even though physics-based source current modeling does not yet incorporate all of the necessary physics and is active work in progress, non-uniform ionospheric source modeling has been successfully applied over the years (e.g., Viljanen et al. 2004) and the use of a physics-based parameterization of the ionospheric currents is now a common practice. Accurate incorporation of Earth’s electrical conductivity in these calculations has remained, to this date, the least explored aspect of the problem.

Here, we aim to provide an overview of the existing techniques, organized both by method and by country. We discuss the variety of GIC studies qualitatively in Sect. 2. In Sect. 3, we go into sufficient details on the GIC estimation method to provide the background for the review that follows. In Sect. 4, we briefly discuss the geoelectric field variability and address the question of its relevance to GIC estimation. We use Sect. 5 to discuss the different ways in which the Earth’s electrical conductivity enters these calculations. In Sect. 6, we describe the various flavors of Earth conductivity models used for geomagnetic hazards and discuss their applicability and limits of validity. Sections 7 and 8 cover, respectively, the global and regional modeling efforts of relevance to GIC studies and assume the form of a literature review. Some final thoughts are provided in Sect. 9.

2 GIC Studies: Flavors and Expectations

The physics of GIC modeling is most naturally decomposed into three distinct steps. The first involves measuring or reconstructing the geomagnetic field at the Earth’s surface. Geomagnetic field measurements at ground-level global geomagnetic observatories are geographically sparse (INTERMAGNET Geomagnetic Data 2013). They can be (and often are) supplemented with variometer measurements (e.g., Gjerloev 2009) which, however, have temporal continuity and reliability challenges. To supplement the sparse measurements and/or for scenario analysis and forecasting purposes, a model for input geomagnetic fields is often constructed. These models provide a means for physics-based geomagnetic field interpolation (e.g., Spherical Elementary Current Systems (SECS); Amm and Viljanen 1999), or for physics-based Sun-to-Earth modeling from first principles (e.g., Pulkkinen et al. 2013; Honkonen et al. 2018), the latter still very much work in progress. Indeed, the most general and forward-looking approach to this problem involves an estimation of ionospheric and magnetospheric sources, from which the induced geomagnetic fields are then obtained. Alternatively, a conceptual magnetic field approximation (based on simplified ionospheric current models such as plane wave or a line current) is sometimes set up. An ultimate goal of this step of GIC modeling is successful validation of the interpolated or modeled geomagnetic fields against observatory or variometer magnetic field measurements.

The second step toward obtaining GIC estimates in a power grid requires geophysical knowledge and involves estimation of geoelectric fields at the Earth’s surface. Specifically, the Earth impedances may be obtained through construction of a 1D (e.g., Ádám et al. 2012; Fernberg 2012) or 3D (e.g., Kelbert et al. 2019) Earth conductivity model, or, alternatively, empirical 3D Earth impedances and interstation transfer functions may be used (e.g., Bonner and Schultz 2017; Campanyà et al. 2019). These different options for choosing the Earth impedances are further discussed in Sect. 5. These are then used, in conjunction with the local geomagnetic field values, to construct local geoelectric fields as a function of both space and time. This may be done in either the time domain (through converting the Earth impedance to an impulse response and convolving it with the geomagnetic field at the location of interest), or via a more straightforward multiplication in the frequency domain. The latter involves a Fourier transform of the magnetic field and the inverse Fourier transform of the output geoelectric field back to the time domain. A conceptual comparison of these two techniques, as applied especially to a near real-time geoelectric field estimation, may be found, for example, in Kelbert et al. (2017). A GIC study that computes geoelectric fields may be reasonably expected to validate their outputs against measured geoelectric fields from magnetotelluric surveys, and/or cross-validate several interpolation or convolution techniques.

Finally, the third step of the calculation, computation of GIC time series in power lines or transformers, involves a power grid network model, which may or may not be accurate and detailed, depending on the goals of the study. Some studies, which focus on comparing different techniques or conductivity models, compute the transmission-line voltages (e.g., Lucas et al. 2018 in the USA). For that, only the power-line geometry is needed. Others go as far as computing accurate transformer-level GICs (e.g., Divett et al. 2018 in New Zealand) and require detailed knowledge of grounding resistances and transformer characteristics. Most studies use a simplified model of the power grid to estimate GICs at the substation level, while also making multiple approximations of the details of the power grid system parameters which are rarely precisely known. Studies that go as far as computing GICs in the substations and/or transformers may be reasonably expected to validate their analysis against one or several GIC measurements in the network. Such measurements are either obtained at the transformer ground connections at substations and buses of the network (i.e., nodes), where the neutral current measurements are made, or via a magnetometer which is placed directly beneath a power line, to produce a proxy for the GIC flowing in the line. Both types of measurements may be used to validate estimated GICs. As we shall see throughout this review, correlations between estimated and measured values are extremely easy to obtain (see Sect. 6 for a discussion as to why), and no sophisticated techniques are needed for that; a uniform plane-wave source and a uniform half-space model for the Earth are sufficient. In fact, such correlations even exist with GIC proxies, which are computed just that way (Marshall et al. 2012). On the other hand, getting an accurate estimate of GIC amplitudes requires a much more sophisticated model for ground-level magnetic fields and the Earth’s conductivity, but is also considerably more useful.

Depending on the ultimate goal of the study, any of these approaches, or a combination, may be used. For the purposes of this review, we are mostly concerned with the role of the Earth electrical conductivity structure (as encapsulated within the Earth’s impedances; step two above) and its influence on the overall outcome of the analysis. To provide sufficient context for this discussion, here we briefly summarize the various conceptual goals pursued by GIC studies around the world and provide some examples.

Goal 1: method development and scenario analysis Most GIC studies belong to this category; Sects. 7 and 8 cover many such works in some detail. While method development may focus on any subset of the operations discussed above, scenario analysis may be considered a critical aspect of the former and is tightly linked to validation. It most commonly refers to computing geoelectric fields, line voltages and/or GICs for a specific historic geomagnetic storm, and comparing with measured data. Once the methodology is validated and deemed sufficiently accurate, such calculations are also made for historic storms of interest where measured GIC data are not available, to calculate the GICs that would have occurred in a particular power network during such a geomagnetic disturbance. Sufficient accuracy for such a historical analysis has only been obtained to date using various transfer function methods that rely on fitting parameters to existing GIC data and are therefore only applicable at a few locations (Ingham et al. 2017 is an example of a rigorous application of such a technique).

Goal 2: extreme-value hazard analysis Extreme-value hazard analysis refers to a statistical technique that compiles a database of prior geomagnetic storms and substorms and uses extreme-value theory or other statistical approaches to predict the geomagnetic and/or geoelectric field amplitudes (and, sometimes, also a power grid network response) of a future extreme but rare event, such as a once in a 100 yrs geomagnetic storm. Such an analysis has been performed at several locations around the world and included work based on a 1D Earth approximation (Canada, e.g., Boteler 2001; Japan, Watari 2015), thin-sheet approximation (UK, Beggan et al. 2013), and for a 3D Earth, in the USA (Love et al. 2016, 2019, 2018a, b utilizing the empirical MT impedances from USArray MT; Schultz et al. 2006–2018). Behavior of extreme geomagnetic disturbances was also studied statistically in Europe by Thomson et al. (2011).

Goal 3: real-time estimation Real-time or near real-time analysis is a goal of all operational geoelectric field and GIC estimation. Indeed, local geoelectric field information supplied in near real-time provides power grid operators in the control room with situational awareness and informs decision-making. It comes with its own unique set of challenges. For example, criteria such as computational efficiency and the speed of data transmission become important. Algorithms need to be robust enough to correctly deal with gaps in real-time data streams, data transmission delays, and site fallouts. Finally, certain conceptual considerations to do with causality come into play. Discussions of these issues may be found in Marti et al. (2014a), Boteler and Pirjola (2017) and Kelbert et al. (2017).

Goal 4: Forecast Short-term, minutes to hours, forecast of GICs is critical for natural geomagnetic hazard mitigation and is a long-standing operational goal in the USA and the UK, among other countries. In fact, developing and validating such a capability is specifically called out in the 2015 US National Space Weather Strategy and Action Plan. However, the methods are far from mature. Forecast of GICs starts with a forecast of geomagnetic fields; if that problem is solved, existing and developing real-time geoelectric field estimation algorithms and efficient power grid system models may be applied to provide GIC forecasts. Space weather forecasts are a hard problem, which is, however, not unsolvable. Recent model validation efforts have shown that both empirical modeling and magnetohydrodynamic (MHD) modeling of the magnetosphere–ionosphere system show promise in predicting surface ${\text {d}}B/{\hbox {d}}t$ (e.g., Pulkkinen et al. 2013; Glocer et al. 2016). Global MHD is especially important because it captures global dynamics using only upstream solar wind inputs and is fast enough to run in an operational setting. Upstream observations of solar wind values are a critical input. Data coming in from NASA’s Advanced Composition Explorer (ACE) spacecraft and the newer Deep Space Climate Observatory (DSCOVR) satellite, which orbit the L1 point between the Earth and Sun, could give the Sun-to-Earth modeling anywhere between 20 min and 1 h of lead time before the storm hits the Earth. (Indeed, even a very slow solar wind (300 km/s) reaches the magnetosphere within about 1 h from L1.) In an operational environment, prompt action such as power load redistribution, possibly involving a decision to take particularly critical or vulnerable grid segments offline, could be used to mitigate long-term damage from a geomagnetic storm. Thus, forecasting is the way forward for operational geomagnetic hazard mitigation, and efficient real-time algorithms will facilitate the adoption of such methods once they mature. Dramatic improvements in regional accuracy of geomagnetic field forecasts would be enabled by real-time data assimilation that takes advantage of both satellite and ground-level magnetic field measurements. However, before data assimilation may be approached, certain deficiencies in the full-physics Sun-to-Earth modeling need to be eliminated. This is work in progress (see, e.g., Sect. 7). Another promising avenue for prediction of coronal mass ejections on the Sun and related geomagnetic field disturbances involves machine learning techniques (e.g., Bobra and Ilonidis 2016; Gruet et al. 2018; Camporeale et al. 2018), providing an alternative approach for future GIC forecasts.

3 Geomagnetically Induced Current Modeling in a Power Grid

GICs flow everywhere in the power grids, entering the grid at the points of transformer groundings, but they are most commonly measured and modeled as the quasi-DC earthing currents flowing in or out of the system at substations. As defined in Lehtinen and Pirjola (1985) and Viljanen and Pirjola (1994), these transformer neutral currents are given by a column matrix

$$\begin{aligned} {\mathbf {I}}^e = ({\mathbf {1}} + {\mathbf {Y}}^n {\mathbf {Z}}^e)^{-1} {\mathbf {J}}^e, \end{aligned}$$

(1)

where ${\mathbf {1}}$ is the unit matrix, ${\mathbf {Y}}^n$ and ${\mathbf {Z}}^e$ are the network admittance and the earthing impedance matrices, respectively, and ${\mathbf {J}}^e$ is a column matrix of induced currents that flow into the ground at earthing points, also known as the “perfect earthing” currents. By Kirchhoff’s law, which applies because the frequencies of the GIC time variation are low enough, its element $J_i$ at an earthing point ${\mathbf {x}}_i$ is defined as the sum of all currents $I_{ki}$ flowing into the node from the power grid. Using Ohm’s law,

$$\begin{aligned} J_i = \sum _{k\ne i} I_{ki} = \sum _{k\ne i} \frac{V_{ki}}{R_{ki}}, \end{aligned}$$

(2)

where $R_{ki}$ are all of the power-line resistances for lines connected to the earthing point ${\mathbf {x}}_i$, $V_{ki} = V_{ki}(t)$ are the induced voltages in the respective power lines, and t is time. These voltages are computed as the line integrals of the time-dependent geoelectric field along the power-line path,

$$\begin{aligned} V_{ki}(t) = \int _{{\mathbf {x}}_k}^{{\mathbf {x}}_i} {\mathbf {E}}(t) \cdot {\mathbf {dl}}. \end{aligned}$$

(3)

It is notable that the geometry of the line may be arbitrary, and an accurate calculation of the voltage would depend on the integration path because the ground-level electric field ${\mathbf {E}}$ is not, in general, conservative. Indeed, the assumption that $\oint {\mathbf {E}} \cdot {\mathbf {dl}} = 0$ over any closed loop would imply that ${\mathbf {E}}$ is irrotational, i.e., $\nabla \times {\mathbf {E}} = 0$, which is not universally true.

Whenever GIC measurements are discussed in the published literature, this most typically refers to the transformer neutral currents at grounding connections at a power grid substation, ${\mathbf {I}}^e$ as in Eq. 1. Since matrix ${\mathbf {Y}}^n$ is fully defined by the resistances in the power lines, $R_{ki}$, and matrix ${\mathbf {Z}}^e$ is defined by the earthing resistances at grounding connections, we see that the GICs at a given point in time are determined by the resistances in the system, the network topology, and the spatial structure of the geoelectric field.

We should also like to note explicitly that, at an arbitrary transformer j, the GIC $I_j(t)$ is a linear combination of the line voltages in the system, $V_{ki}(t)$. If we define the matrix ${\mathbf {K}}^e = ({\mathbf {1}} + {\mathbf {Y}}^n {\mathbf {Z}}^e)^{-1} = [K_{ji}]$,

$$\begin{aligned} I_j(t) = \sum _{i} K_{ji} J_i(t) = \sum _{i} \sum _{k\ne i} \left[ \frac{K_{ji}}{R_{ki}}\right] V_{ki}(t). \end{aligned}$$

(4)

In conjunction with Eq. 3, this expression makes it clear that GICs are a function of the spatially and temporally variable geoelectric fields, integrated over a set of complicated paths. Moreover, we also see that GIC at a given point depends on line voltages within the whole grid. In practice, regions far away from a specific location have only a small effect (e.g., Pirjola 2005). This makes it easier, for example, to study national power grids while applying crude approximations to the (often unknown) details of the power grid systems in the neighboring countries.

Here, we focus on the computation of geoelectric fields, as that is the component that requires geophysical knowledge. Specifically, the Earth’s electrical conductivity (or its inverse, the resistivity of the subsurface) and the grounding resistances at substations are two factors that play a critical role in these computations. The resistivity of the Earth’s crust and mantle is determined by the geology of the region. However, it is notable that the grounding resistances also depend upon the local geology. In fact, while it is intuitively clear that higher GICs would be observed in less conductive regions (giving rise to higher geoelectric fields), if the near-surface is also a poor conductor then the effect may be compensated by higher grounding resistances. Based on these considerations, the highest GICs would be observed in regions with low grounding resistances overlaying a highly resistive structure, such as a sedimentary basin over a deep cratonic root. Indeed, the low grounding resistances would allow the currents induced in the Earth to enter the power grid, while the high ground resistivity would also direct these currents, preferentially, into the power grid.

Until a few years ago, power grid system modeling software was utilizing the approach of Viljanen and Pirjola (1994) under the uniform electric field assumption. Specifically, the engineering software would take in a homogeneous electric field of a certain amplitude. Since the direction of the electric field would be undetermined, the software would rotate the field through all possible azimuths and compute the maximum possible GIC due to an electric field of this amplitude, at each of the grid nodes. A modest modification of this approach allowed the user to define a plane-wave or an electrojet model at a certain height and of a certain amplitude and compute the spatially (but not temporally) variable geoelectric field from that model. In this case, a 1D electrical conductivity profile could also be input to the software, to determine the Earth impedance needed for geoelectric field computation. More advanced variants of the same approach allowed the input of several 1D conductivity profiles beneath different regions of the power grid, for a piecewise 1D approach to the Earth impedance computation, to match the capabilities described by Viljanen and Pirjola (1994).

While there is some utility to this approach for preventative scenario analysis, this utility is limited. As we shall show, the GIC at every substation at every point in time during a space weather event is actually dependent on a combination of spatiotemporal variations of geomagnetic fields and spatial variations of the electrical conductivity of the subsurface, which jointly determine both the amplitude and the direction of geoelectric fields. Earth’s 3D electrical conductivity structure, which determines the Earth impedance, plays a primary role in these calculations, modulating the amplitude, polarization, and phase of the geoelectric field. While some bounds on the GICs can be explored with these simplistic models, using, for example, a maximum observed rate of change of geomagnetic fields for a given geomagnetic latitude, these estimates run the risk of being misleading unless the true Earth impedance is taken into account.

In the past 5 years or so, efforts have flourished to upgrade these tools to incorporate a user-supplied spatially and temporally variable geoelectric field. This new set of methods, while still very raw, has great potential for improved accuracy in scenario analysis. Moreover, this new approach can also be successfully utilized for real-time estimation and forecasting. However, the few software tools that support this capability are of a commercial nature and are therefore not accessible to the academic community worldwide. Neither are the detailed power grid system configurations, which, with a few exceptions, are generally considered proprietary and sensitive information. This situation resulted in widespread use of alternative approaches, as we shall see in Sect. 8, which are, nevertheless, all based on a variant of Eq. (1).

4 Spatial Structure of Geoelectric Fields and the “Smoothing Effect” of Power-Line Integration

Spatial variability of geoelectric fields is determined by non-uniformity of the geomagnetic field, combined with a nonlinear dependence on the 3D Earth’s electrical conductivity structure. Even in the case of a 1D ground model, spatial variations of the magnetic field can produce a strongly non-uniform electric field (e.g., Viljanen et al. 1999; Viljanen and Pirjola 2017). Vice versa, a realistic 3D ground conductivity model makes the electric field strongly non-uniform even when a plane-wave magnetic field is assumed (e.g., Bedrosian and Love 2015). The relative contributions of these two effects have not, to the author’s knowledge, been rigorously explored; however, simultaneous magnetic and electric field measurements collected as part of USArray MT surveys worldwide (see Sect. 6) exhibit orders of magnitude greater site-to-site variability in the electric field versus that in the magnetic field measurements, suggesting that the bulk of the variability found in the former is governed by the Earth’s structure. Spatiotemporal structure of geoelectric fields, when applied to a power grid, determines the GIC hazard. It is therefore important to better understand geoelectric field variability and its causes.

High-quality, long-term geoelectric field measurements are notoriously difficult to obtain. The three geomagnetic observatories in Japan, Kakioka (KAK), Kanoya (KNY), and Memambetsu (MMB) have pioneered long-term high-quality geoelectric field data collection; geoelectric field data recordings at these three observatories commenced, respectively, in the years 1932, 1948, and 1949, and are still ongoing. Similarly, geomagnetic observatory Nagycenk (NCK) in Hungary has been recording geoelectric field data since 1957, and in digital form since 1996 (Kis et al. 2007). Other geomagnetic observatories that collected geoelectric field measurements over the years include Chambon-la-Foret (CLF; France), Gan (GAN, Maldives), Eskdalemuir (ESK), Hartland (HAD), Lerwick (LER) in UK, Valentia (VAL, Ireland), and Niemegk (NGK, Germany). Recently, the U.S. Geological Survey initiated 1 s geoelectric field recordings at Boulder (BOU) geomagnetic observatory in Colorado, USA (Blum et al. 2017).

In spite of these efforts, such long-term measurements are rare. However, geoelectric field data are also recorded over short periods of time, generally, days to months, every time an MT site is installed. Such geoelectric field measurements therefore sample thousands of locations worldwide. These time series can be used to analyze qualities of the real geoelectric field variations, such as variations in amplitude and direction. We know well from magnetotelluric measurements that geoelectric fields vary dramatically over short length scales, often within the footprint of a power grid (e.g., see Schultz et al. 2006–2018).

Other methods to analyze geoelectric field variations include modeling, which can be performed using real magnetic field measurements (see Sect. 5 for details), or by studying the response of geoelectric fields to synthetic magnetic field waveforms. As is seen in Bedrosian and Love (2015) in the simplified example of a synthetic sinusoidal, spatially uniform geomagnetic field variation, both the amplitude and direction of the geoelectric field vary dramatically (by 2 orders of magnitude) in space and in time entirely due to the spatial variations in Earth’s electrical conductivity. Realistic geoelectric fields obtained with real geomagnetic field measurements and estimated using magnetotelluric impedances also exhibit pronounced variability (e.g., Cuttler et al. 2018; Love et al. 2019, 2018a). Depending on the configuration of the power grid, certain directions of geoelectric field cause significantly higher transmission-line voltages (and therefore, GICs) than others; moreover, the geoelectric field direction that would cause a high GIC value varies between substations (e.g., Viljanen and Pirjola 1994, their Figure 6). Similar conclusions have been reached by a Gannon et al. (2017) study which compared simulated GICs using Fernberg (2012) 1D models vs Schultz et al. (2006–2018) in the Pacific NW of the USA. This work is in some disagreement with other published and unpublished reports for reasons discussed in Love et al. (2018b), but it does emphasize the role of geoelectric field orientation in GIC estimation.

In spite of clear evidence for significant spatial heterogeneity of geoelectric fields, it has been repeatedly argued in the literature (e.g., Viljanen et al. 2012) that heterogeneity of the Earth, as reflected in the electric fields, is not a primary factor for GIC estimation because of the smoothing effects of the power-line integration in Eq. 3. Of course, in reality (recall Eq. 3), GICs depend on geoelectric fields not at a single point, but along each of the power-line paths. Here, the assertion is that realistic 3D electric fields, when integrated, produce line voltages that are approximately the same as the average 1D estimate of the electric field would provide; thus, a good 1D approximation is sufficient. We have also seen discussed as early as Viljanen and Pirjola (1994) that significant random perturbations of up to 50% of the average electric field are effectively smoothed out in the power-line integration, producing line voltages that are almost indistinguishable from those for the unperturbed electric field. Viljanen et al. (2012) also propose that this averaging effect, in general, justifies a 1D approximation for the electrical conductivity, even though different 1D models might need to be used in different geological regions.

While it is certainly true that a synthetic electric field with random variations against a known average will average itself out in the process of power-line integration, it is naive to assume that this will always happen in reality, where conductive and resistive structures that matter are often of comparable spatial scales to the power grid footprint, are not distributed randomly, and may be aligned, or not, with any of the power lines. See Appendix A for a discussion of the fortunate, but unlikely, situations when 3D and 1D Earth models produce identical GICs. That this is a rare event at least in geologically complicated regions is well illustrated by the analysis presented in Lucas et al. (2018).

Fig. 1

a Geomagnetic field at the USGS Fredericksburg geomagnetic observatory during the 1989 geomagnetic storm. b Snapshot at 21:30 13 March 1989 of the calculated voltages across the transmission lines in the mid-Atlantic region of the USA with the use of empirical 3D impedances (Schultz et al. 2006–2018). c Calculated voltages across the transmission lines using regional 1D impedances (Fernberg 2012). Note that the colormap is logarithmically scaled, indicating the wide range of voltages that are generated across transmission lines. The circled region is an area where the use of empirical 3-D impedances produces higher voltages compared to the voltages calculated with 1D impedances. Also, a few transmission lines have smaller voltages across them when 3D impedances are used, for example, just to the west and southwest of the circled region (Lucas et al. 2018)

Lucas et al. (2018) estimated the voltages across just over 500 transmission lines exceeding 150 kV in the mid-Atlantic region of the USA for several historic geomagnetic storms (March 1989 and July 2000), using for simplicity the data stream from a single geomagnetic observatory, at Fredericksburg, VA, USA, as input to geoelectric-field estimation algorithms. This study compared the voltages computed using the commonly used (Fernberg 2012) 1D electrical conductivity models, relative to 64 real Earth impedances ${\mathbf {Z}}$ measured through the magnetotelluric component of the NSF’s EarthScope USArray project (Schultz et al. 2006–2018). Figure 1 shows the voltage distributions in the power lines in these two cases. While there is certainly substantial similarity between these two sets of computations, the differences (e.g., as circled) are just as substantial and are determined by regional, small-scale geologic structures as they align with or cross the power lines. We should also note that the specific case of 3D effect caused by the land/ocean contrast cannot be eliminated by power-line integration, affecting particularly those lines that are close to and run parallel to the coast (see, e.g., Ivannikova et al. 2018; Liu et al. 2019a for relevant discussions).

While it could be argued that the 1D models that are used here for comparison have room for improvement, there are deeper reasons for these differences (see, e.g., Bedrosian and Love 2015; Lucas et al. 2018 for discussion). The physics of 3D induction is distinct from that of induction in a layered setting, even if a large collection of layered models is provided. While the Earth impedance ${Z}^{1}$ obtained from 1D modeling is a frequency-dependent scalar, the 3D Earth impedance ${\mathbf {Z}}$ is a tensor and has four complex components at each frequency. This structure of the impedance allows for all kinds of distortions of the electric fields, both in amplitude and in direction, relative to the 1D approximation. Because of this, the geoelectric field exhibits significant spatial variability that is not, in general, captured by the uniform electric field approximation. Moreover, an inherent assumption of 1D modeling is that the layers extend to infinity, laterally. A collection of local 1D models violates this assumption and should itself be modeled in 3D (see, e.g., Liu et al. 2019a for a recent analysis of this effect).

5 Data and Model-Space Approaches to Geoelectric Field Estimation

Now that we have established the need for an accurate geoelectric field estimation algorithm, let us consider several alternative approaches to this problem.

Fig. 2

Two alternative workflows for estimation of gridded geoelectric fields (purple) are illustrated in this flowchart. Starting with the empirical geomagnetic and magnetotelluric data (green), one could either use interpolation methods to estimate real-time gridded geoelectric fields (“data-space workflow,” shown in blue; e.g., Bonner and Schultz 2017; Campanyà et al. 2019), or a conductivity model could be used for a physics-based variant of the interpolation (“model-space workflow,” shown in red, e.g., Marshall et al. 2019). Method 1: interpolation using multi-station transfer functions; not based on physics. Method 2: spherical elementary currents (SECs; Amm and Viljanen 1999) or other physics-based interpolation (Rigler et al. 2019). Method 3: multiplication in frequency domain or convolution in the time domain. Method 4: MT inversion. Method 5: numerical modeling based on Maxwell’s equations. Method 6: interpolation, e.g., using Delaunay triangulation; not based on physics

The flowchart presented in Fig. 2 [modified from Kelbert et al. (2019)] illustrates two alternative workflows for estimation of gridded geoelectric fields. For example, estimation of real-time, ground-level geoelectric fields starts with two empirical data types: a time series of geomagnetic field measurements (such as the U.S. geomagnetic observatory data; USGS 1901) and the magnetotelluric (MT) transfer functions (such as the USArray MT data; Schultz et al. 2006–2018). Post-analysis and forecast approaches could use a very similar framework but, possibly, a different magnetic field data stream. Using the geomagnetic field time series at a limited set of locations, some version of interpolation may be used to estimate geomagnetic fields at the locations of the MT data sites. These estimates may be multiplied (in the frequency domain) or convolved (in the time domain) with the MT transfer functions to obtain geoelectric field time series at the locations of MT sites. Then, another interpolation technique may be applied to estimate geoelectric fields on a grid. This method does not require an electrical conductivity model, but may suffer from possible drawbacks due to the fact that the geoelectric field interpolation it requires is not based on any underlying physics.

Alternatively, using the geomagnetic field time series measured at geomagnetic observatories, one could interpolate directly to a grid, to obtain gridded magnetic fields at ground level. Independently, one would use the MT transfer functions to obtain a 3D electrical conductivity model via the method of MT inversion (e.g., Avdeev 2005; Egbert and Kelbert 2012) and use the numerical MT forward modeling technique to obtain MT transfer functions directly on the grid. This would amount to a rigorous, physics-based interpolation method. The gridded magnetic fields and gridded MT data would then be multiplied (in the frequency domain) or convolved (in the time domain) to estimate ground-level geoelectric fields directly on the grid.

Both pathways provide real-time gridded geoelectric fields to the power grid industry, an operational goal identified, for example, by the interagency Space Weather Operations, Research, and Mitigation (SWORM) Task Force in support of the power grid industry in the USA. Ground-level geoelectric fields could be used as the input to power grid flow models, both operationally, for situational awareness in the control room, as well as for historical scenario analysis that allows a power grid to investigate potential points of failure in their grid in advance of a space weather event. The use of 3D MT transfer functions as part of the workflow has been shown to provide an improvement in ground-level geoelectric fields relative to the use of 1D conductivity models (e.g., Weigel 2017; Bonner and Schultz 2017; Kelbert et al. 2017; Campanyà et al. 2019). We suggest that interpolation based on a numerical solution of Maxwell’s equations using a 3D conductivity model such as that of the contiguous USA (CONUS; Kelbert et al. 2019) could provide a further improvement against the utilization of empirical MT transfer functions interpolated to the locations of interest. Cross-validation of these two alternative techniques is a subject of ongoing space weather research.

Each of the studies we discuss in Sects. 7 and 8 adopts either the data-space or the model-space approach to geoelectric fields and/or GIC estimation problem. In the past, it has been most common to use a variant of a model-space approach, oftentimes with simplifying approximations (see Sect. 6). In recent years, however, several nationwide gridded MT data surveys, such as USArray MT in the USA (Schultz et al. 2006–2018) and AusLAMP in Australia (Stolz 2013), and the wider availability of MT data in general (Kelbert et al. 2011) have made the alternative data-space approach more appealing. These developments have enabled derivation of fully 3D electrical conductivity models that are now becoming available for geomagnetic hazard applications. In fact, while the plane-wave approximation inherent to MT impedances could be a cause of inaccuracy in geoelectric field estimates during a major magnetic storm when this assumption is violated (e.g., up to 15% during the peak of the Halloween storm, Kelbert et al. 2017), direct use of 3D Earth conductivity models in conjunction with non-uniform ionospheric sources suffers from no such drawbacks (Fig. 3).

Fig. 3

A generalized model-space formulation that bypasses the MT impedances and the plane-wave source assumption to model the gridded geoelectric fields with non-uniform ionospheric and magnetospheric sources. Two alternative workflows (Method 7 + Method 9 and Method 8 + Method 9) are shown [Alexey Kuvshinov 2019, pers. comm.] Method 7 assumes some parameterization of the source (using spherical harmonics, electric dipoles, elementary currents, or alternative basis functions) and some type of regularization and requires the knowledge of Earth’s conductivity model. Method 8 assumes that radial external field is known on a grid. Method 9 involves numerical solution of Maxwell’s equations with a given source (of any complexity) and a given 3D conductivity model. Method 7 + Method 9 workflow was employed by Püthe and Kuvshinov (2013a) and Püthe et al. (2014) to compute geoelectric fields on a global scale, and Method 8 + Method 9 combination was used by Honkonen et al. (2018) and Ivannikova et al. (2018) to compute geoelectric fields on global and regional scales, respectively

6 Flavors of Earth Conductivity Modeling Approaches

Electrical conductivity of the Earth’s lithosphere is the primary factor affecting the frequency content of GIC variations. Geomagnetic disturbances typically last between several hours to several days, causing most intense GICs in the period range between 10 s and several thousand seconds. However, as suggested by Kappenman (2003), signal power applicable to geomagnetic hazards spans a wider range of periods that extends to tens of thousands of seconds. This translates to electrical conductivity variations that extend from the near-surface to about the bottom of the mantle transition zone (at 670 km depth). Lateral scales of relevance are harder to constrain, but at a minimum they are commensurate with a typical footprint of a power line and could therefore range between several kilometers to several hundreds of kilometers. As is known from MT inversions, these spatial scales exhibit lateral variability of electrical conductivity of up to 3–5 orders of magnitude. Additional discussion of spatial and temporal scales of relevance to the GIC problem may be found, e.g., in Love et al. (2018c) and Kelbert et al. (2019). Many regional electrical conductivity models have been compiled for geomagnetic hazard mitigation. However, since the depths of interest are greater than any MT models can provide, global models are useful to fill in any missing information at depths, and for global-scale geomagnetic hazard analysis.

Here, we briefly discuss the types of Earth conductivity models used in GIC studies and expand on their applicability and limitations. Usage examples for many of these methodologies are provided in Fig. 4 and are further discussed in Chapter 8.

Fig. 4

Regional GIC studies categorized conceptually by methodology and the approach to Earth’s electrical conductivity modeling. A detailed review of these analyses is provided in Sect. 8. Shaded areas correspond to formulations that would not be self-consistent or otherwise optimal. Only studies that involve modeling and/or comparison with GIC data are included in this summary. Rather than being exhaustive, this list is intended as a useful reference for the reader

1D approximation is a very common approach to GIC modeling which involves an assumption that the Earth’s electrical conductivity may be considered, for all practical purposes, laterally homogeneous. If the lateral geoelectric field is assumed spatially constant, a simplified GIC estimation procedure holds (Viljanen and Pirjola 1994; Pulkkinen et al. 2007). At a single transformer j, the GIC value at any point in time is

$$\begin{aligned} I_j(t) = a_j E_x(t) + b_j E_y(t), \end{aligned}$$

(5)

where $a_j$ and $b_j$ are coefficients that characterize a single transformer or a power line (hereafter referred to as the (a, b) parameters) and depend on the resistances and geometry of the power system, and $E_x$ and $E_y$ are the northward and eastward components of the geoelectric field at that location. This assumption is valid if the ionospheric current can be approximated as a spatially uniform plane wave, and the electrical conductivity of the Earth is one-dimensional and the same everywhere.

A range of possible approaches may be applied to obtain the (a, b) parameters. The most common approach is to estimate these linear coefficients based on the power grid model (e.g., Viljanen et al. 1999; Beggan et al. 2013; Beggan 2015; Viljanen and Pirjola 2017). Other methods, particularly applicable in the absence of such information, but when selected GIC measurements are available, involve least-squares parameter estimation (e.g., Pulkkinen et al. 2007; Ngwira et al. 2008). In spite of the considerations put forward in Sect. 4, the 1D GIC approximation in Eq. (5) often gives extremely satisfactory results, fitting local GIC measurements with a high level of accuracy. Let us examine why this simple approximation sometimes works well, and explore its limitations.

For example, Pulkkinen et al. (2010) employed an empirical parameter estimation method of Pulkkinen et al. (2007) to very accurately model the GIC in a highly inhomogeneous area of Hokkaido, Japan. Pulkkinen et al. (2010) used a 1D conductivity model for the area and the magnetic field measurements from a neighboring Memanbetsu geomagnetic observatory, to compute a 1D version of the geoelectric field at this location. They proceeded to fit the (a, b) parameters to a measured GIC data stream. They then used Eq. 5 to compute the GIC values at the same location for a different time period. We should note that both the GIC measurements and the computed values were downsampled from 1 s temporal resolution to 60 s temporal resolution, thus removing both the unphysical spikes and some of the signal, somewhat reducing the complexity of the problem. Still, the parameter fitting algorithm worked incredibly well and resulted in GIC estimates with the misfit centered on zero and the maximum mismatch of less that 0.3 A. Based on these appealing results, it is tempting to conclude that the 1D approach may be applied in any, even the most 3D, geological setting. Let us, however, consider more carefully the physics of this approximation.

As we recall from Sect. 3, a physics-based computation of the GIC at a transformer requires that the geoelectric field is integrated along all of the incoming power lines and summed up as in Eq. 2, before Eq. 1 is applied. All of these are linear operations and can be reduced to Eq. 5 if geoelectric field ${\mathbf {E}}$ is uniform over the footprint of all power lines connected to transformer j (Appendix A). In the 1D approach, the space dependence is gone, resulting in a simplified Earth impedance: $Z_{xx} = Z_{yy} = 0$ and $Z_{yx} = - Z_{xy} = Z^1$. Let us, however, assume that the geoelectric field is non-uniform over the area, but that the magnetic field can be assumed uniform over the same area. As is shown in Appendix A, this results in a very similar expression to that of the 1D approximation, in that in both cases, in the spatial sense (though not in the temporal sense), the GIC at a transformer location can be considered a linear combination of the horizontal components of ${\mathbf {B}}$. In addition to the network topology and configuration, in the presence of 3D geometry, these coefficients depend of the 3D MT impedances in the area, while for the 1D approximation they depend only on $Z^1$. Both the electric fields ${\mathbf {E}}$ and the voltages V will be dramatically different under the 1D Earth assumption. However, since in both cases the GICs have the same mathematical form, parameter fitting easily compensates for the inadequacy of this approximation. This explains the great success of this approach in GIC modeling even in highly non-1D settings.

Unfortunately, the hard limits of this simple and powerful parameter fitting approach make it impractical for use in geomagnetic hazard mitigation. Since the (a, b) parameters are fit to existing data, rather than computed from first principles using realistic 3D Earth impedances, system geometry, and grounding resistances, this method needs GIC measurements at all locations of interest. Indeed, parameter fitting can only be applied in selected locations where prior GIC measurements have been taken; moreover, the parameters are only valid until the grid configuration changes. It is a common but unfortunate misconception that since the parameter fitting method works well in this very narrow range of circumstances, it must also be true that the 1D approximation holds and may be applied away from a location of prior GIC measurements—a leap of reasoning exhibited, e.g., by Liu et al. (2013). That turns out not to be the case, since anywhere else in the grid, the application of Eqs. 1–3 cannot be avoided and requires an accurate estimate of geoelectric fields. This accurate estimate cannot, in general, be obtained with a 1D Earth approximation. On the contrary, when derived from a given power grid model, the (a, b) parameters are exact under the assumption of a spatially uniform electric field. However, as the uniform electric field assumption is commonly broken in realistic geological settings, this approach has to be distinguished from the parameter fitting approach in both the procedure and the modeling outcomes.

A complementary parameter fitting approach bypasses the (a, b) parameters and instead uses the full power grid system model in conjunction with local GIC measurements to linearly invert for the local electric field (e.g., Kazerooni et al. 2013). This approach requires no electrical conductivity model but is implicitly assuming a 1D Earth, since an electric field obtained at one location is then used in a regional sense.

In light of the above, extra care needs to be taken when power grid modeling software is adapted from the use of a uniform geoelectric field, to using more realistic, both spatially and temporally variable geoelectric field estimates. Prior parameter calibrations no longer hold when the full physics of the problem is modeled. Thus, any system configuration model that employed the parameter fitting technique would need to be redesigned to adapt the generalized 3D approach discussed above. For example, a recent work by Butala et al. (2017) carefully compared the use of 1D to measured 3D impedances at several locations in the USA and validated these approaches against observed GICs. They found no significant improvement when the 3D Earth impedances were used. However, in their GIC modeling, the authors utilized a version of the (a, b) parameter model of Eq. 5, which they recognized as a limitation. For the 3D version of their approach, they used that same equation, while replacing the local 1D impedance ${\mathbf {Z}}$ with the local 3D impedance. Since that equation does not involve integration of variable 3D impedances along the power-line paths like Eq. 10 does, there is no reason that this modified 1D approach would work any better than the original.

A similar, but rigorous and explicit, variant of this approximation has been illustrated in Scotland by Thomson et al. (2005), and more recently in New Zealand by Ingham et al. (2017). Thomson et al. (2005) defined transfer functions between nearby geomagnetic observatory data and measured GICs (see also McKay 2003) and applied them quite successfully to reproduce the GICs recorded at the same transformers a couple years later, in spite of the system configuration changes that occurred over this time frame. In Ingham et al. (2017), transfer functions were developed between the rate of change of local geomagnetic field and measured GICs. These were then employed to estimate GICs at the same transformers for other times when the measurements were not available. Just like with the more traditional, (a, b) parameter fitting approach, this technique cannot be used in places other than these few transformers, but it is a great tool for validation of other approaches.

Fig. 5

Europe (Ádám et al. 2012, openly available from http://real.mtak.hu/id/eprint/2957) (left) and the USA (Fernberg 2012, as adapted by the U.S. Geological Survey) (right) are subdivided into regions, with boundaries defined by the outlines. Each region is associated with a distinct 1D electrical conductivity model, resulting in 1D compilations which are used for GIC modeling

Another common approach found in the literature is the use of 1D compilations, also known as piecewise layered Earth models (e.g., Marti et al. 2014b). These compilations are collections of 1D models of Earth’s conductivity, derived from a variety of sources; notable examples are Ádám et al. (2012) in Europe and Fernberg (2012) for the USA (Fig. 5). For example, the Fernberg (2012) 1D model compilation has been used by Wei et al. (2013) to estimate ground-level geoelectric fields for several historical geomagnetic storms. A similar, but finer resolution, MT-based 1D compilation has been obtained by Blake et al. (2016) for Ireland. 1D compilations have the virtue that GICs can be locally computed in the 1D sense. Thus, while covering a wide area on a map without assuming identical 1D structure everywhere, computationally this approach is no different than the more traditional 1D approach. The drawbacks of this approach, however, are notable. As discussed in Sect. 4, this assumption cannot account for the tensor form of a realistic Earth impedance and is inherently non-physical.

The approach taken by Viljanen et al. (2012) for Europe can be taken as a guide of the validity of Eq. (5) for system-wide GIC estimation. There, the authors use a collection of 1D models (Ádám et al. 2012) and the local Eq. (5) approximation, computing the (a, b) parameters from the grid system configurations, rather than fitting them to existing GIC data. The authors validate their 1D approach first at the geoelectric field level, comparing against measured E-field data at Nagycenk, Hungary, and second, to GIC measurements at a substation in Vykhodnoy, Russia. They find that their geoelectric field estimate is extremely good, validating the three-layer 1D conductivity model in Hungary, at least at periods relevant to 1 min temporal resolution of the data. Their GIC estimates, however, are far from perfect. The authors accurately suggest that there may be multiple reasons for it, including the unknowns in the system parameters, the (two- to five-layer) electrical conductivity, and data noise. Finally, a better result is obtained for GICs estimated in a pipeline in Finland; a problem that, as discussed in the Introduction, seems to lend itself to the 1D Earth approximation well.

Fig. 6

Seawater surrounding an island creates a significant electrical conductivity variation that is causing a complicated coast effect. These example subsurface models are inspired by bathymetry, sediment thickness, and surface and bedrock geology. Fujita et al. (2018, Reuse permission obtained from Kakioka Magnetic Observatory) (left) created a fully 3D electrical conductivity model using bathymetry and sediment thickness in and around Japan. The top 1 km is shown. Beggan et al. (2013, Reuse permission purchased from Wiley) (center) defined a thin-sheet model for the UK based on geological information. Surface conductance is shown. Divett et al. (2017, Modified from their Fig. 1) (right) developed a block-wise thin-sheet model of New Zealand that is similarly inspired by geology and bathymetry

Thin sheet approach uses an infinitely thin layer of laterally variable conductance on top of a 1D model (Price 1949). This method was developed in Cartesian coordinates by Vasseur and Weidelt (1977) and McKirdy et al. (1985), as well as a multi-sheet Cartesian solution by Fainberg et al. (1993). Spherical coordinates, global thin-sheet EM modeling approach was first presented by Fainberg et al. (1990) and further developed by Kuvshinov et al. (1999). In the context of geomagnetic hazards, the Vasseur and Weidelt (1977) technique was applied in the UK (e.g., Thomson et al. 2005; Beggan et al. 2013; Fig. 6, center) and New Zealand (e.g., Divett et al. 2017; Fig. 6, right). Independently, a global thin-sheet approximation was also applied to geomagnetic hazards (e.g., Püthe and Kuvshinov 2013a; see Sect. 7) and a variation with distortion parameter fitting was developed for Japan by Püthe et al. (2014). We should note that both of these studies used a fully 3D global modeling code as described in Kuvshinov (2008), but applied it to a thin-sheet model of the Earth. While the thin sheet approach is not a fully 3D formulation in the sense that 3D Earth conductivity variations are not supported other than through an infinitely thin heterogeneous layer at the top, it is similar in the sense that it provides a means to account for the lateral inhomogeneity of the Earth’s conductivity.

Finally, fully 3D modeling is rather rare in the modern GIC estimation landscape, but is coming into fashion as new MT surveys are being completed and the regional electrical conductivity models are obtained. Some examples include the recent work in Japan (Fujita et al. 2018; Nakamura et al. 2018; Fig. 6, left) and the modeling study by Ivannikova et al. (2018) in the UK. In the USA, a fully 3D conductivity model compilation exists (Kelbert et al. 2019), to be used for geomagnetic hazards.

Coast effect analysis defines a distinct class of studies focused on specifically analyzing the effect of the coastline conductivity contrast on geoelectric fields. Classical works on this topic include Swift and Wescott (1964), Wescott (1967), and Parkinson and Jones (1979). More recently, Gilbert (2005) and Goto (2015) addressed this topic in the context of geomagnetic hazards. Liu et al. (2019a, b) discussed the implications of the coast effect and similar sharp lateral conductivity contrasts to geoelectric fields in the immediate proximity of the contrast, using conceptual synthetic modeling. While these and similar studies are invaluable for improving our fundamental understanding of the behavior of geoelectric fields near seawater to land interfaces, these are not intended for practical geoelectric field estimation. Instead, it is a way to isolate, numerically and analytically, the coast effect from the rest of the 3D structure affecting these fields. For realistic geoelectric field computation, the coast effect should be considered jointly with the rest of heterogeneity in 3D electrical conductivity of the Earth’s crust and mantle. Specifically, the ocean coastlines, bathymetry, and conductivity can be modeled as part of the thin-sheet and 3D conductivity models, as discussed in Sect. 8 [and addressed in more detail in Püthe and Kuvshinov (2013a), Ivannikova et al. (2018), Liu et al. (2018b), and Pokhrel et al. (2018)].

7 Global Geoelectric Field Estimation

7.1 Models

Global 3D electrical conductivity models in existence to date (Kelbert et al. 2009; Semenov and Kuvshinov 2012; Sun et al. 2015) resolve electrical conductivity structures of the Earth’s mantle, globally, approximately between 410 km (upper transition zone boundary) and 1200 km depth. They have all been obtained by rigorous 3D inversions of INTERMAGNET Geomagnetic Data (2013) ground geomagnetic observatory data and parameterized by layered spherical harmonics. Due to the period ranges in which the simplifying source approximations are applicable (several days to 100 days), they have limited to no resolution in Earth’s crust and upper mantle.

Several well-developed inversion strategies that can be applied to invert global geomagnetic observatory data (Kelbert et al. 2008; Kuvshinov and Semenov 2012) and several other recent and in development global EM inversion methods applicable to satellite data (Püthe and Kuvshinov 2013b is the only such study published to date) will, no doubt, allow us to dramatically improve on our imaging of Earth’s electrical conductivity, globally, in the upcoming decade. As our understanding of ionospheric source structure is better incorporated into the modeling, we also envision significant improvements in the resolution of global electrical conductivity of the upper mantle.

In contrast, local and regional 3D electrical conductivity models of interest span the depths of several kilometers to several hundred kilometers ($\sim 410$ km depth at the maximum). These models have been obtained through a variety of methods, but primarily with 3D magnetotelluric inversion techniques (e.g., Siripunvaraporn et al. 2005; Egbert and Kelbert 2012). Recent nationwide MT surveys such as USArray MT in the USA, Schultz et al. (2006–2018); AusLAMP in Australia, Stolz (2013); SinoProbe in China, Dong et al. (2013); and a number of European and Canadian initiatives that are just coming into fruition provide unprecedented potential for dramatic improvements in regional 3D electrical conductivity models. Unfortunately, of these, only the USArray survey data are publicly available to the international academic community to date. Efforts to invert such data sets in a regional sense make it possible to compile these regional electrical conductivity models into a coherent framework that could be used globally for geomagnetic hazard mitigation purposes. The only published effort of this kind that is relevant to global GIC modeling is Alekseev et al. (2015). This work combined a global bathymetry distribution, sedimentary cover thicknesses, and several inverse models from MT in USA (Meqbel et al. 2014) and Europe (Korja et al. 2002), as well as a regional geomagnetic observatory and magnetometer array data from Australia (Wang et al. 2014) into a fully 3D, 100-km-deep global conductivity model's compilation at $0.25^\circ \times 0.25^\circ$ resolution. Maintaining the same approach with more region-specific defaults and a more comprehensive compilation of inverse models would provide a comprehensive global 3D map of electrical conductivity in the Earth’s crust and uppermost mantle. This map could also be expanded into the upper mantle at 100–410 km depths using MT models at regional to continental scales, and global to regional magnetometer data inversions.

7.2 Methods

Prediction of space weather hazards involves full physics, high-resolution geospace modeling, and data assimilation. Development of comprehensive methods for global geoelectric field modeling in the framework of geomagnetic hazard mitigation is very much an ongoing effort. While GIC estimation efforts are often limited to nationwide borders because of the nature of human-made infrastructure (Sect. 8) and the related data sensitivities, geomagnetic disturbances are a global phenomenon and need to be modeled globally to enable forecasts. Space weather modeling of the magnetosphere–ionosphere system is driven primarily by solar wind observations (Sect. 2) and provides forecasts of the global ionospheric current systems, as well as the ground-level geomagnetic fields. These forecasts are a work in progress and do not currently capture regional dynamics accurately enough to be of practical use. Methods that could ingest these forecasts and generate a full-physics global geoelectric field in real time are also being developed. Here, we outline several parallel efforts.

The most comprehensive work on this, to date, belongs to the research group at ETH Zurich in Switzerland, most notably Püthe and Kuvshinov (2013a), Püthe et al. (2014), and Honkonen et al. (2018). Global magnetic and geoelectric field modeling of Püthe and Kuvshinov (2013a) and Püthe et al. (2014) is (somewhat misleadingly) referred to as 3D modeling. Indeed, the software that is used in this approach (a variant of Kuvshinov and Semenov 2012) is capable of modeling fully 3D Earth structures, albeit at the expense of significant computational time. However, in this particular set of studies, the Earth model that is employed, instead, follows the thin sheet approach, resulting in a formulation similar to Vasseur and Weidelt (1977) for regional studies (used also by, e.g., Thomson et al. (2005) and Beggan et al. (2013) in the UK and by Divett et al. (2017) in New Zealand), except in this case it is a global modeling technique. A major difference between this and other thin sheet approaches is its flexibility with respect to the source representation. Since the magnetic and electric fields are both linear in terms of the source, it makes sense to represent the global ionospheric sources as a combination of spherical harmonics and pre-compute the global magnetic and electric fields distribution, in the frequency domain, for each spherical harmonic coefficient. An appropriate linear combination of the output “unit” magnetic and electric field representations, then, allows for a fast computation of the fields with arbitrary ionospheric sources, making this a powerful technique for practical modeling.

Püthe and Kuvshinov (2013a) successfully validate this method against magnetic field measurements at several coastal geomagnetic observatories. They also discuss the difference between the 1D and the thin sheet + 1D (with 3D modeling) approaches. Then, Püthe et al. (2014) go further by reproducing the geoelectric field measurements at the three geomagnetic observatories in Japan (Fujii et al. 2015). Notably, the coast of Japan exhibits rather extreme electrical conductivity variations, while the conductivity model employed for the Püthe et al. (2014) study is a $1^\circ \times 1^\circ$ thin sheet, underlain by a 1D profile. This resolution is clearly not sufficient to reproduce the effect of 3D electrical conductivity variations in Japan. The authors address this issue by introducing a frequency-independent (and, therefore, time-independent) distortion matrix at each of the observatory locations:

$$\begin{aligned} {\mathbf {E}}^\mathrm{{observed}}(t) = {\mathbf {G}} {\mathbf {E}}^\mathrm{{modeled}}(t), \end{aligned}$$

(6)

where ${\mathbf {G}}$ is a $2\times 2$ matrix that is unique to each geographic point of measurement. The data are fit with a least-squares method to obtain these distortion matrices ${\mathbf {G}}$, which is then used to correct the modeled geoelectric field data at that location. The logic behind this approach includes the traditional consideration of galvanic and inductive effects; the distortion matrix is meant to separate out the galvanic effects (due to accumulation of electric charges along conductivity contrasts), which often happen at small enough scales that they are corrected for in many traditional MT studies. The authors demonstrate that this method allows them to fit the observed (detrended) geoelectric fields quite well. Even though the amplitudes of measured geoelectric field perturbations are not perfectly reproduced, the overall signal waveform, and even the amplitudes of the distortions in geoelectric fields, are roughly recovered (Fig. 7).

The authors proceed to apply their approach in forecast mode to assess its applicability for accurate geoelectric field prediction based on the $D_{st}$ forecast from the Advanced Composition Explorer (ACE) satellite for the October 2003 “Halloween” storm. They demonstrate that at present the (10 min) $D_{st}$ forecast lacks temporal resolution and energy, to reproduce real geoelectric field variations.

We find the approach of Püthe et al. (2014) very promising. We would, however, like to note that the strategy of compensating for high-resolution, fully 3D electrical conductivity by distortion matrix fitting may not be sufficient in many cases where real induction effects cannot be warped into a galvanic distortion matrix (an example most clearly illustrated by the method’s inability to reproduce geoelectric field amplitudes at MMB; Fig. 7). Specifically, the frequency-independent galvanic distortion cannot, in general, compensate for the frequency-dependent inductive effects. The approach, as presented, may be directly applicable on the seafloor and in other select locations where a coarse thin-sheet conductance model is a good approximation to Earth conductivity. For general practical use, however, the distortion matrices would need to be calculated at each location of interest, which would in theory mean everywhere on land, in a gridded manner. They can only of course be computed in places where geoelectric field measurements have been previously obtained. In some sense, this makes the approach a close cousin of, e.g., Pulkkinen et al. (2010), which requires GIC measurements at a location in order to accurately model new GICs there, by fitting the GICs first to proxy geoelectric field equivalents. Having said that, geoelectric time series already exist everywhere at MT site locations, so the distortion matrices ${\mathbf {G}}$ may be computed at each modern MT site, globally, except that the measured time series and not the MT impedances are then required. This possibility makes the approach of Püthe et al. (2014) much more appealing, especially if the baseline 3D conductivity model can be dramatically improved (and the corresponding unit magnetic and electric fields recomputed). Unfortunately, incorporation of model improvements is, by far, the most expensive part of the approach.

Finally, Honkonen et al. (2018) couple the approach of Püthe et al. (2014) to the Space Weather Modeling Framework (SWMF; Tóth et al. 2005) running at NASA’s Community Coordinated Modeling Center (CCMC; Pulkkinen et al. 2013). SWMF is used to compute geomagnetic fields at the Earth’s surface. These are then converted to equivalent sheet ionospheric current representation in terms of spherical harmonics, as in Kuvshinov and Semenov (2012), their Appendix G. The global $1^\circ \times 1^\circ$ thin-sheet model configuration of Püthe et al. (2014) is forced by these equivalent currents; frequency-domain geoelectric fields are computed and inverse-Fourier-transformed into time domain. A geomagnetic storm of December 2006 is run through the system and compared against both geomagnetic and geoelectric measured data at several geomagnetic observatory locations around the world. It is encouraging that SWMF can produce ground magnetic fields whose amplitudes and time derivatives have realistic magnitudes. However, while the concept is well developed and comprehensive, the results clearly show that predictive capabilities of SWMF as used in this study are not yet suitable for regional forecasts. The Earth conductivity model probably also needs to be substantially refined, as discussed above. A similar approach was applied to regional-scale modeling by Ivannikova et al. (2018).

Fig. 7

Computed electric field data at geomagnetic observatories KAK (left) and MMB (right) for the Halloween storm of Oct 2003 are compared to the measurements. Computations are performed using the local geomagnetic field measurements, the global thin-sheet conductance model, and local distortion matrices. Reuse permission obtained from Püthe et al. (2014), their Figs. 4 and 6

There are comparable ongoing efforts in the USA, though none of them are ready for prime time yet (e.g., Welling et al. 2017). In 2016, the U.S. Department of Energy funded the Shielding Society from Space Weather (SHIELDS) initiative through the Los Alamos National Laboratory (LANL) Laboratory Directed Research and Development (LDRD) program. The SHIELDS framework is being developed by a cross-disciplinary team of scientists in the fields of space science and computational plasma physics. SHIELDS global Sun-to-Earth modeling framework at LANL now includes a global EM modeling component; for proof of concept, they currently work with the Alekseev et al. (2015) model compilation. For a description of other ongoing efforts, please refer to Sect. 8.17.

8 Regional Geoelectric Field and GIC Estimation

Because of the national nature of power grids, it is typical for GIC studies to be focused on a specific country, rather than on a broader geographic region. We therefore structure our discussion geopolitically, sort country names alphabetically, and within each subsection, we describe multiple methods, if they exist. Within Europe, we order our discussion roughly geographically, North to South. As outlined in Sect. 2, some of these studies go all the way to GIC estimation, some stop short of GICs and estimate the voltages in the lines, and others yet focus of geoelectric field estimation. Here, we are most interested in the effect that the choice of electrical conductivity model plays in geoelectric field and GIC estimation. We also note that many extremely interesting studies that do not explicitly discuss Earth’s electrical conductivity are, regrettably, outside the scope for this paper. Such studies are not brought up except where they are directly relevant to our discussion.

8.1 Australia

Geomagnetically induced current studies in Australia in relation to power grids commenced relatively recently, but there are a few studies of significance that warrant discussion. Relevant work on the continent to date was performed in the frequency domain and was geared more toward method development, fundamental understanding, and validation than toward realistic real-time estimation. The analysis by Marshall et al. (2013) is a pilot study, reporting the first GIC observations in the Australian power network. GIC monitoring was initiated at eight strategic locations in eastern Australia, in collaboration with the Australian Energy Market Operator Limited. In spite of their proximity to the coastline, in this particular analysis the Earth was approximated as a uniform half-space, and the uniform field (a, b) approximation was utilized (see Sect. 6). For this simple approximation, transformer neutral currents were found to be very well correlated with the west-to-east component of the geoelectric field estimate, $E_y$, demonstrating a direct relationship of these currents to space weather. The scaling of the GIC modeling part would need to be adjusted with more sophisticated techniques in order to be useful. However, these valuable measurements demonstrate that GICs produced by space weather are readily observable in parts of the Australian power grid and that their amplitudes are sufficient to suggest the need for establishment of geomagnetic hazard mitigation procedures in Australia.

Independently, in a rare study involving a full 3D electrical conductivity model, Wang et al. (2016) computed frequency-domain geoelectric fields in Australia for the south-to-north and west-to-east plane-wave sources (corresponding to the two polarizations in 3D magnetotellurics). This was accomplished by using a pre-computed 3D electrical conductivity model of Australia (Wang et al. 2014) based on the AWAGS magnetometer array (Chamalaun and Barton 1990) to forward model 3D magnetotelluric impedances at the quasi-regular AWAGS locations covering the continent of Australia. In a similar manner, a 1D model and a simpler 3D model were explored; the latter only comprised the land/ocean contrasts. It was shown that inland electrical conductivity structures in northern Australia, which reduce the conductivity contrast between coast and ocean, also notably reduce geoelectric fields relative to models that do not account for inland conductivity. Modeled fields were also reduced in central Australia. In contrast, modeled geoelectric field was significantly increased on the resistive Western Australia craton, throughout the period range of the study.

This work, however, should still be considered qualitative because of the inherent plane-wave approximation for the magnetic field that was employed. The next logical stage of GIC modeling in Australia would incorporate these 3D impedances into operational geoelectric field and, subsequently, GIC estimation, fueled by real-time geomagnetic measurements on the continent. Validation against measured geoelectric and GIC data would form an important part of this development. This path forward seems especially attractive in view of the comprehensive ongoing long-period MT data collection initiative in Australia, AusLAMP (e.g., Stolz 2013; Robertson et al. 2016, 2017), which upon completion should cover the whole of the Australian continent with approximately 2,800 estimated 3D MT impedances on a quasi-regular 50 km-by-50 km grid spacing. These impedances could either be used directly for geomagnetic hazard mitigation purposes, or the evolving high-resolution electrical conductivity models (e.g., Robertson et al. 2018) could be used to forward model the impedances, in a manner similar to Kelbert et al. (2019) in the USA (using the data-space and model-space workflows of Fig. 2, respectively).

Indeed, a state-of-the-art validation analysis of GIC modeling using multiple different conductivity model formulations in Australia was recently published by Marshall et al. (2019). Four electrical conductivity models were set up with varying degrees of complexity, ranging from a uniform conductivity model, to a couple of 1D models, to a fully 3D model, derived from AWAGS magnetometer data as in Wang et al. (2014), from which gridded impedance data were then produced (as described in Fig. 2). These were compared with AusLAMP MT data (Duan and Kyi 2018) at several locations, and a good correspondence was observed. The four different geoelectric field models were then applied to the same elaborate power network model within the east Australian power system, and the power-line integration is performed as described by Lehtinen and Pirjola (1985). The modeling for four different geomagnetic storms is validated against observed GIC time series at up to eight different locations within the network. While some storms and locations were reproduced very well, the unusually thorough (for a GIC modeling study) validation allowed the authors to observe that none of the models did consistently well everywhere and at all times. Overall, the 3D model was found to perform better than the simpler conductivity models, but with some caveats. According to Marshall et al. (2019), the relatively good performance of all models could be attributed to the reduced coast effect and otherwise relatively high conductivity of the Earth beneath the footprint of this power grid. The 3D conductivity model resolution could also be substantially improved by applying community MT inversion software (e.g., Kelbert et al. 2014) to the newly obtained AusLAMP MT data; this higher-resolution model could then be combined with the coarser AWAGS-derived model at depths.

8.2 Azerbaijan

We were happy to find rich literature in the Republic of Azerbaijan in the South Caucasus pertaining to space weather effects on power grids (Babayev et al. 2006), and also, as a surprising aside, on human brain activity, cardiovascular and psychological health (e.g., Babayev 2003; Babayev and Allahverdiyeva 2007 and multiple others). Azerbaijani studies show mature yet qualitative understanding of the role of Earth’s electrical conductivity in natural geomagnetic hazards.

8.3 Brazil

One of the important demonstrations of GICs in a low latitude region, Trivedi et al. (2007) presents the first GIC measurements in Brazil that were obtained during the November 2004 geomagnetic storm, in cooperation with the electric power company in Brazil, FURNAS, and the federal government national research institute, Instituto Nacional de Pesquisas Espaciais (INPE). Direct GIC monitoring was initiated at a later time at one of the transformers on the network; however, these data do not cover November 2004. Instead, GIC measurements presented by Trivedi et al. (2007) are indirect, calculated from measured magnetic field variations using sensors installed under two transmission lines in southeastern Brazil. Geoelectric fields were also estimated for the same time period, using a 1D conductivity model that was informed by a 2D MT profile in between the lines, and the corresponding inverse 2D conductivity structure. Just like Marshall et al. (2013) in Australia, similarity between geoelectric field estimates and the GIC measurements was observed, but no scaling was yet attempted.

More recently, the same GIC measurements were used for validation of a more rigorous GIC estimation technique in de Silva Barbosa et al. (2015). Even though existing localized magnetotelluric studies show significant 3D structure in the northeast corner of the Brazilian grid (Padilha et al. 2014), as well as some deep anisotropy in the southeast (Padilha et al. 2006), a comprehensive MT data set covering all of Brazil does not yet exist. At present, 1D conductivity approximation is utilized for these studies. Specifically, de Silva Barbosa et al. (2015) created a four-layer 1D conductivity model and validated it against the GIC measurements described in Trivedi et al. (2007), using a formal system grid model (Viljanen and Pirjola 1994). As usual with such analyses, the timing of the GIC extremes was estimated very accurately, while the amplitudes are somewhat off as compared to the values presented in Trivedi et al. (2007), especially as far as reproducing all the peaks. This model, however, fares better than a half-space. As is pointed out by the authors, another unexplored complication stems from the presence of the South Atlantic Magnetic Anomaly (SAMA) above Brazil, linked with an increase in the electron precipitation in this area. Nevertheless, based on the existing model and approximations, the maximum GIC amplitude for the October 2003 “Halloween” storm is estimated at 25 A. Based on this publication, we understand that the current work is focused on method development and is not intended to provide operational quality estimates, until the sensitivity of GIC estimates to the geologic setting and power network characteristics is better explored.

Finally, a recent analysis by Espinosa et al. (2019) focused on the effects of both ionospheric and ground electrical conductivity on GICs at the geomagnetic equator in Brazil. To circumvent the relative lack of MT data in Brazil, four magnetometers from the Embrace MagNet network (Denardini et al. 2018) were chosen at such locations that the underground geological structure was amenable to 1D interpretation. Four corresponding 1D electrical conductivity models were set up at these locations. In an approach similar to Viljanen et al. (2013), a power grid model was set up based on a 500 kV network operated in central Brazil, with some of the network parameters adopted from de Silva Barbosa et al. (2015), and the others set to default values that provided the best fit to measured GICs. By analogy, we assume that power-line integration of the electric fields was performed; however, the boundaries of the four conductivity models are not indicated in Espinosa et al. (2019), and a detailed discussion of the procedure is absent. The single modeled-to-measured GIC comparison presented in the paper suggests that this methodology is on a par with other international results based on a similar approach and can sometimes reproduce up to 60—80% of the variations. The effect of the ground conductivity was also isolated by using the same input magnetic field at all four locations. As expected, the relative GIC magnitudes are controlled by the electrical conductivity profile (surface impedance) and the sampling rate of the geomagnetic field (1 sample per minute, in this case). Espinosa et al. (2019) suggest that crustal variability, particularly that of the upper crust, is primarily responsible for the observed variations due to the Earth’s impedance.

8.4 Canada

Located at a higher geomagnetic latitude and having experienced a major GIC-related blackout (Hydro-Québec in 1989; Allen et al. 1989; Boteler 2001; Bolduc 2002; Guillon et al. 2016), Canada is very cognizant of geomagnetic hazards. In Canada, the Natural Resources Canada (NRCan) is the organization responsible for communicating forecasts and warnings on geomagnetic storms based on models and observed data from the Canadian Automatic Magnetometer Observatory System (CANMOS), NOAA satellites, and other space-based and ground-based observations (e.g., Canadian riometer and GNSS networks for monitoring the ionosphere and radio communication; Larisa Trichtchenko, 2018, written pers. comm.). NRCan’s Canadian Geomagnetic Forecast Service has been in operation since 1974 (Trichtchenko et al. 2009). Over the years, NRCan scientists have authored a most comprehensive body of work on a number of GIC-related topics, such as: (1) historical damage from GICs, e.g., Boteler (2003) and Boteler (2006), (2) modeling of GIC impacts on pipelines and power grids, e.g., Trichtchenko et al. (2010) and Boteler and Pirjola (2017), and (3) statistical geomagnetic hazard assessments (including the extreme-value statistics) in Canada and worldwide, e.g., Nikitina et al. (2016).

NRCan telluric hazard assessments for power grids and pipelines, as well as telluric hazard maps (NRCan 2006) and online geoelectric field tool (since 2006), are based on the 1D Earth’s electrical conductivity models inferred from a number of magnetotelluric studies available in the open literature at each particular period of time. For example, a large collection of 1D models obtained from 400 MT sites in the Cordillera using the Niblett–Bostick transformation were summarized and used by Boteler (1997), Boteler et al. (1998), Boteler (2001), and Boteler (2003) to analyze geomagnetic hazards to the five major power systems in Canada. For example, Hydro One in Ontario, Canada, has collaborated with NRCan to implement their own, rather sophisticated, control-room real-time GIC estimation procedures, all based on 1D Earth model compilations (Marti et al. 2014a, b). Subsequent statistical assessments of the geoelectric field performed since 2017 were based on the updated 1D Earth models (Trichtchenko and Fernberg 2012; Trichtchenko et al. 2015 and references therein). More recently, Trichtchenko (2016) authored a review on the status of geoelectric field modeling in Canada used by the NRCAN online tool, which employs a detailed collection of 1D models to compute geoelectric fields at or near the CANMOS observatories. Merging and updating 1D models currently used for telluric hazard assessment in specific areas of interest (such as the footprints of power grids and pipelines) into one comprehensive model across Canada is work in progress. One suggestion for possible improvement on this approach would be to use the available MT data from the Lithoprobe project (Jones et al. 1984–2003, 2014), possibly supplemented with additional gridded surveys, to attempt 3D geoelectric field modeling, via either the data-space or the model-space approaches (Fig. 2).

8.5 China

As is evident from the GICs observed in many other low- to mid-latitude countries (e.g., Australia, New Zealand, South Africa, Brazil, and others), GICs are a real problem at latitudes relevant to Chinese infrastructure, even if the nature of the currents is somewhat different from that at higher latitudes—for example, they may more often be observed in conjunction with the sudden commencement phase of the storm. At first glance, it may seem unexpected that, for a large country with well-developed infrastructure, China has comparatively few publications covering the topic in the international literature. However, as is explained in Liu et al. (2013), the impact of GICs on power grids in China is a relatively new problem, due to the rapidly increasing scale of the power grid in China in recent years. High-voltage power lines tend to be more susceptible to space weather effects, and in China, a 500 kV system had become the backbone of the large-scale power grid by the year 2000. Moreover, power demand and distribution are highly unbalanced in China, where both coal and hydropower are focused in the West, while the load is primarily in the East, and, to mitigate the situation, a 1000 kV system is being built from the West to the East of China by the year 2020 (Liu et al. 2013). Higher-voltage transmission lines also have lower DC resistances per unit length, making them more vulnerable to GICs.

Consequently, in 2012 China initiated the National High Technology Research and Development Program (also known as the 863 Program), which started its space weather operations at the end of 2014, launching a multi-point GIC monitoring campaign which was able to record GIC signals for a March 2015 geomagnetic storm at two substations on the coast of China (Zhang et al. 2015). Even before the new monitoring campaign, as early as 2001, neutral current monitors were installed at several transformers in China to evaluate some abnormal vibration and transformer noise that, at the time, did not have a known cause. One of these monitors was operating during the November 2004 magnetic storms and recorded GICs as high as 75.5 A at Ling’ao power plant (over three transformer phases; Liu et al. 2014). Liu et al. (2013) estimated how dramatically these GIC measurements would scale to the Carrington-sized event (i.e., similar to September 1859 geomagnetic storm), in the new 1000 kV system.

Thus, GIC measurements in China exist, even though they are rather sparse at present. Comprehensive MHD modeling codes have also been developed and validated in China (e.g., Zhang et al. 2015). However, methods related to geoelectric field estimation still have room to develop. Several GIC method development papers have explored the 1D (a, b) parameter fitting approach (see Sect. 6) in China (Liu et al. 2009; Zhang et al. 2015). The original method development paper, Liu et al. (2009), implemented the 1D approximation to the Lehtinen and Pirjola (1985) approach and obtained very reasonable GIC data fits by fitting the (a, b) parameters directly to measured GICs. However, it is a mistake to consider that this method can, therefore, be reliably used with geoelectric fields computed from geologically inspired 1D conductivity models as in Liu et al. (2013, 2014). In fact, Liu et al. (2013) make important GIC amplitude estimates based on an understanding that the methodology “has been proved correct” by Liu et al. (2009). This is a dangerous misunderstanding; indeed, that the parameter fitting method works well for GIC amplitude estimation is no indication that the same 1D approach will work without the benefit of parameter fitting. See Sect. 6 for a discussion of this common misunderstanding. Unfortunately, the Liu et al. (2014) study that uses this method to compute GICs in the Chinese Northwest has no measured GICs to compare to, as all GIC measurements to date seem to have been obtained on the Southeast coast of China. A recent study by Liu et al. (2018b) used these data to cross-validate a time-domain finite element 3D Earth modeling approach against the previously implemented 1D formulation and found that accounting for the coast effect noticeably improves the GIC estimates compared to the use of a geological 1D model. They also note that China’s complex tectonics would likely cause substantial lateral inhomogeneities in the Earth’s conductivity inland as well as in the coastal regions, and conclude that 3D Earth conductivity models should be adopted for GIC estimation efforts in China.

Finally, the ongoing SinoProbe project (Dong et al. 2013, 2014) is collecting comprehensive MT measurements covering all of China. These data, if released, would allow rapid development of advanced new GIC estimation methods in China that would have the benefit of the 3D MT impedances and models, similar to those that are covering the USA (e.g., Bonner and Schultz 2017; Kelbert et al. 2019). Unfortunately, the Chinese MT data are not openly available to the international research community.

8.6 Ethiopia

Senbato et al. (2018) and Liu et al. (2018a) developed a full node modeling of GICs in Ethiopia based on the Ethiopian power grid system model. A subset of 11 geomagnetic storms between the years 2003 and 2017 has been analyzed statistically, and the high- and low-probability extreme GIC values were estimated, with the low-probability extreme-value GIC exceeding 20 A at 89% of substations. This comprehensive approach, however, makes a rather simplified assumption on the electrical conductivity structure of the Ethiopian lithosphere. A single 1D conductivity model is used for Ethiopia based on De Villiers et al. (2016), which is, in turn, loosely adopted from the local profile at Addis Ababa out of the Alekseev et al. (2015) compilation at depths above 100 km and a global profile Civet et al. (2015) constrained by satellite data below 100 km. Notwithstanding the consideration that a 1D profile may not be an appropriate electrical conductivity model for all of Ethiopia, another concern might stem from the fact that the electrical resistivity of the upper crust (0.1–8 km depths) is set to 2e5 Ohm·m, a number substantially higher than is reasonable for the East African Rift. This is likely a misread from Alekseev et al. (2015) which suggests more like 2e3 Ohm·m at these depths. This factor of 100 difference could have substantially increased the estimates of GIC hazards in this area.

We should also note that the electrical conductivity of Africa in Alekseev et al. (2015) is not inferred from magnetotelluric data [(as suggested by De Villiers et al. (2016)] but rather combined from the distribution of sediments and a default lithospheric electrical conductivity value. By contrast, geothermal fields in Ethiopia have been extensively studied with magnetotellurics, see e.g., Whaler and Hautot (2006), Desissa et al. (2013), Didana et al. (2015), Samrock et al. (2015), Hübert et al. (2018) among others, consistently recovering three orders of magnitude variations in electrical conductivity throughout the Earth’s crust (and sometimes lithosphere).

8.7 Europe: The Nordic Countries

The Nordic countries, specifically Finland, Norway, and Sweden, have a long history of GIC modeling, thanks to unusually strong cooperation with the power grid industry. GIC recordings in the earthing wires of Finnish high-voltage power system stations have been actively studied and monitored in the framework of a tight collaboration between the Fingrid power company (formerly Imatran Voima Oy) and the Finnish Meteorological Institute (FMI) (Viljanen and Pirjola 1994). Recordings were initiated in 1977, and mostly included analogue and digital 10 s recordings at one to several stations. In the years 1991–1992, a more focused study was performed, where digital 10 s recordings were measured at five sites of the 400 kV power grid. There have been some sporadic recordings also after 1992. This wealth of information allowed for pioneering method development and validation (e.g., Lehtinen and Pirjola 1985; Wik et al. 2008; Viljanen et al. 2012), primarily based on the 1D conductivity approach (Sect. 6).

A landmark study by Viljanen et al. (2004) demonstrated the use of the uniform field, 1D Earth (a, b) approximation to the power grid system parameters (see Sect. 6) and concluded that this fast and simple method gave good and consistent results for both geoelectric field and GIC estimation. Indeed, at the time of the publication of this work, it was a breakthrough in practical GIC estimation, and the methodology is widely used to this day. In a later work, Viljanen et al. (2012) set up a prototype model of high-voltage power grids covering all of Europe. This study performs the end-to-end modeling of both ground geoelectric fields, and the GICs in the prototype power grids in all of Europe, using the 1D collection of ground conductivity models from EURHOM (Ádám et al. 2012), and the (a, b) approximation to the system (see Sect. 6). However, some limitations of this approach can also be pointed out.

Consistent with other validation and method comparison studies that do not involve direct parameter fitting (e.g., Wik et al. 2008; Viljanen et al. 2012), Viljanen et al. (2004) illustrate that the quality of geoelectric field estimates obtained using this 1D approach is variable. Specifically, while at some locations this approximation happens to work adequately well (e.g., site B22 of the BEAR array), in others (e.g., B13, B16) the geoelectric field is dramatically overestimated by about a factor of 4, while in others yet (e.g., B24) it is underestimated by about the same factor. Similar results are observed for GIC estimates in the power grid, but GIC estimates in the pipelines are nearly perfect (see also Pulkkinen et al. 2001 and Viljanen et al. 2006 for examples of successful application of the 1D GIC estimation approach in pipelines).

The latest developments include full 3D modeling of GICs in Sweden using three heterogeneous layers to represent upper crust, mid- to lower crust, and the upper mantle, overlaying a uniform Earth (Rosenqvist and Hall 2019).

8.8 Europe: The British Isles

The British Isles are composed of the UK, which includes mainland UK (Great Britain) and Northern Ireland, and the Republic of Ireland, which is politically independent. While the geology and geophysics of Ireland and the UK are closely interlinked, the two countries have independent national power transmission grids and independent sources and strategies of science funding. Even though we fully recognize and respect the sovereignty of both Ireland and the UK, for the purposes of this review, which is focused on geophysics, we discuss these two countries, as well as several smaller islands nearby, jointly as the British Isles. Where relevant, such as for the historical overview, we resort to separate discussions.

Great Britain and Ireland are both islands, surrounded by the shallow sea shelf on the East, and by the deeper Atlantic Ocean to the West. It was therefore clear from the onset that the GIC estimates in the UK would be influenced by the so-called coast effect (Sect. 6) in a manner that could not be accommodated by anything other than 3D modeling. GIC studies in the UK commenced in the late 1990s when several confidential technical reports, on 3D electrical conductivity and on GIC risk assessment, were commissioned to the British Geological Survey (BGS) by the National Grid in the UK (David Beamish 2018, written pers. comm.). This work culminated in the earlier National Grid models and in the so-called thin-sheet conductance model of the British Isles (partly informed by MT data, which were then available mainly in Scotland) that were publicly reported in Beamish et al. (2002) and McKay (2003). While full 3D analysis was computationally prohibitive at the time, the thin-sheet conductance approach based on Vasseur and Weidelt (1977) at 20 km resolution was used for geoelectric field modeling. A variant of this model was also used in the thorough development and validation study of GICs in the Scottish Power grid by Thomson et al. (2005). The thin-sheet conductance map, at 10 km resolution, is referred to by the name of AMK2003 in later literature (e.g., Beggan 2015).

Both the model of the National Grid and the thin-sheet conductance model were substantially improved in the following decade. The newer thin-sheet conductivity model of the British Isles (referred to as BGS2012, Beggan 2015) was derived by the BGS from the analysis of the resistivity properties of bedrock material based on the BGS’s 1:625,000 geological map of the UK and Northern Ireland (Fig. 6, center; Beggan et al. 2013). We should note that while BGS’s detailed tectonic map of the British Isles, and the corresponding conductance map, includes Ireland as well as mainland UK, until recently all published GIC studies have covered mainland UK only.

The first published GIC study in Ireland (Blake et al. 2016) did not make use of this conductance map, instead follow an entirely different paradigm of using a collection of MT-inspired 1D models for their work. For a couple of smaller storms discussed in this work, GIC estimates obtained with this approach were substantially improved compared to using the Ádám et al. (2012) model collection or a uniform half-space. Since then, two additional comprehensive analyses of GICs in the Irish power network appeared in the press. Blake et al. (2018) designed an elaborate power grid model for Ireland and tested its sensitivity to various critical parameters. This included a thorough investigation of the effects of the grounding resistances on GIC estimates using the network’s response to a uniform geoelectric field. Earth’s electrical conductivity was beyond the scope of this study. Most recently, a novel geoelectric field modeling approach based on local and interstation transfer functions was developed and validated in Ireland and the UK by Campanyà et al. (2019).

In addition to work that focused specifically on the UK, larger-scale European studies also included modeling of the islands, albeit with a broad brush. The 1D model compilation by Ádám et al. (2012) includes a low-resolution set of 1D conductivity models in the UK, used later by Viljanen et al. (2013) and by Kelly et al. (2017) for large-scale European analysis of GICs in the UK and elsewhere. Finally, Ivannikova et al. (2018) published a more conceptual comparison of 3D, thin sheet, and 1D modeled geoelectric fields in the UK, with the emphasis on the coast effect.

To conclude, we would like to mention exciting recent developments in the UK and Ireland that should, over time, dramatically improve GIC modeling in the isles, as well as allow a transition to a fully 3D approach. In May 2017, the Natural Environment Research Council (NERC) has funded a major 4-year effort in the UK by a consortium of 10 institutes, led by the Geomagnetism Team at the BGS, to assess Space Weather Impacts on Ground-based Systems (SWIGS). Within the framework of this project, both MT and GIC measurements are now being made to improve the UK conductivity model and to validate BGS’s UK-wide GIC network model. BGS also works in partnership with National Grid plc (the system operator) and EDF Energy (a major electrical utility in the UK) and has worked with them for several years now, allowing them access to some company-measured GIC data. BGS’s GIC network model is being used to monitor potential “hot spots” in the grid. In Ireland, there is a parallel project under development funded by the Geological Survey of Ireland with the aim of covering the island of Ireland with MT data for Space Weather applications and for modeling Ireland’s lithosphere (Joan Campanya 2018, written pers. comm.).

8.9 Europe: Spain

Spanish GIC research commenced in 2012 (Torta et al. 2012) with a focus on northeastern Spain and has since been expanded to all of Spain by Torta et al. (2014). Specifically, a model of the entire 400 kV Spanish power network was set up under certain reasonable assumptions about the earthing and transformer winding resistances.

A couple of layered models were assumed for the Earth and compared to a uniformly resistive Earth model at 0.001 S/m conductivity. GICs predicted under these assumptions were compared with GIC measurements at one location, and Torta et al. (2014) decided that the uniform Earth model fared somewhat better. That uniform Earth model was then used to estimate GICs at each substation for two major storms, October 29–31, 2003 (“Halloween storm”), and March 24, 1991. This yielded maximum amplitudes of 70 and more than 100 A, respectively, for these storms in Spain. Unfortunately, as has been argued by, e.g., Bedrosian and Love (2015) as well as here in Sect. 6, while uniform or layered Earth assumptions allow computation of geoelectric fields and GICs that are well correlated with measurements, the amplitudes may be off by up to several orders of magnitude, depending on the local electrical conductivity structure of the true Earth.

This modeling approach was significantly improved by Torta et al. (2017). There, a single MT impedance tensor was acquired in the proximity of the Vandellòs substation of the Spanish national power grid, which is located in a coastal area. The authors then suggest that this measured impedance can be considered as the representative response of the regional geoelectrical structure for long periods and use it in a regional sense to reflect the 3D structure in the area. While there is no further elaboration or verification of this proposition, the local 3D geoelectric field is obtained, and the formula provided in Eq. 5 is used to compute the GICs at the substation. This formula, as we discuss in Sect. 6, is strictly only applicable to 1D Earth situations. However, it provides a very satisfactory fit to measured data, that is far superior to the estimates obtained with the best 1D approximation of this impedance tensor. Additionally, a 2D approximation of the impedance is also considered, after a formal distortion correction, and is found to provide a fit comparable to the 3D approach.

We suggest that, in this coastal area, where the geoelectric strike is well known, a 2D impedance tensor that is corrected for near-surface distortion may indeed provide a good representative regional response. In this very special case, the 1D approximation for the GICs, Eq. 5, may also be directly applicable to the 2D case. The 3D approach applied under these same restrictions does not provide any additional benefits compared to the 2D approximation, but it may do so once the fully 3D approach to GIC computation (Eq. 10) is implemented, using more than a single MT impedance to describe the Earth’s structure. Overall, this work constitutes highly valuable validation and sensitivity analysis applying 1D, 2D, and 3D MT impedance to compute GICs in a coastal region. In the future, additional MT data acquisition, combined with shedding the 1D GIC approximation of Eq. 5 in favor of Eq. 10, may provide improved results everywhere in Spain.

8.10 Europe: Austria

Several studies of interest have been performed in Austria in recent years. The earlier studies in Austria established basic models of GICs in the Austrian power network managed by the Austrian Power Grid (APG) and analyzed the first sets of measurements (Halbedl et al. 2014; Bailey et al. 2017; Halbedl et al. 2017). Measurements of direct currents flowing through the neutral points of transformers have been carried out across the network by the APG and Technical University of Graz, starting in 2014 with one specific node and adding four more stations in 2016 and 2017.

The most recent study, Bailey et al. (2018), is a comprehensive analysis of GICs in the Austrian power network based on modeling using the thin-sheet approach. The 2D surface conductivity model utilized for the thin-sheet was obtained by the Austrian Geological Survey (GBA) using airborne EM data and has been described in detail in Bailey et al. (2017) and Schattauer et al. (2017). Bailey et al. (2018) isolate the Austrian power network from the rest of Europe for GIC modeling purposes and use a careful integration of the electric fields over power lines to obtain GICs. The comparisons between the modeled and measured data at the three stations in Austria over a day’s worth of measurements in 2017 show an impressive correspondence.

8.11 Europe: Russia

The author has been pointed to a wealth of scientific publications by Russian authors on the potential vulnerability of electric power grid equipment to geomagnetic storms. Unfortunately, the vast majority of these works is published in national journals and does not make its way to the international scientific discourse. This speaks to the unfortunate reality of very limited scientific communication between Russia and the West in the past three or four decades. To the best of the authors’ knowledge, all of these publications have been focused on the engineering aspects of the issue and are therefore not directly relevant to the discussion at hand.

Having said that, a few recent exceptions to the rule certainly deserve a mention. For example, Belakhovsky et al. (2018) analyze the time and spatial variability of geomagnetic fields at high latitudes, as compared to the GIC measurements at nearby locations. They explore the equivalent current systems causing this variability and conclude that their spatial complexity implies that electric power lines stretching not only in the longitudinal but also in the latitudinal directions are affected by even a moderate geomagnetic storm at these latitudes. We should note that substantial GIC values of up to 40 A in magnitude are recorded for this storm (St Patrick’s Day of 2015) in a power grid that extends entirely latitudinally. Kozyreva et al. (2018) further suggest that the key factor of GIC excitation is the medium-scale and fast-burst component of substorm evolution. Vorobev et al. (2019) explore the statistical relationships between these GIC measurements and the time variability of local geomagnetic fields using the log-normal and generalized Pareto distributions. As might be expected from the physics of the problem, they find high correlations between these values and, indeed, estimate that a linear regression model between the two parameters can infer the local GIC values with a very modest error of a few Amperes, on average, for a moderate-sized storm. This is consistent with the observation of Pulkkinen et al. (2006) on the linear relationship between the GIC values and the local time variations of the geomagnetic field. Independently, Sokolova et al. (2018) performed a preliminary risk analysis of the Siberian unified power grid, based primarily on the network configuration and equipment condition. These authors do not currently invoke electric field estimates or electrical conductivity models in their analyses.

8.12 Europe: Italy

A preliminary GIC risk analysis for Italy has been recently performed by Tozzi et al. (2019). This work is based on the estimation and study of the so-called GIC index (Marshall et al. 2011, see also Sect. 8.15), which serves as a proxy for geoelectric fields. “Low” to “moderate” risk levels were estimated over most of the Italian power grid, with a possibility of significant damage during a future large geomagnetic storm. Earth’s structure has not been explicitly considered in this analysis.

8.13 Europe: Turkey

Kalafatoglu Eiguler and Kaymaz (2017) present measured magnetic and geoelectric field time series from a single MT station that they installed in Turkey during several geomagnetic disturbance periods and discuss the relationships between the time variations of the magnetic field, geoelectric fields, and certain solar wind parameters. A qualitative interpretation in terms of 1D Earth conductivity structure is provided.

8.14 Japan

Japan, a country located at a low geomagnetic latitude, was not considered to be at risk for GICs until recently. GIC measurements in Japan are scarce and proprietary, and the only GIC recording available for scientific analysis was obtained at a transformer in central Hokkaido (an area in the north of Japan) and discussed thoroughly by Pulkkinen et al. (2010) and Watari (2015). The data showed a small GIC signal of several Amperes even in the middle of a moderate geomagnetic disturbance. Furthermore, there have never been any reports of serious GIC-induced power-line failures in Japan (Fujita et al. 2018).

However, the susceptibility of the Japanese power grid to an extreme geomagnetic event has never been studied. The Japanese government is acutely aware of the importance of generating resilience to extreme natural disaster scenarios. In 2015, the Ministry of Economy, Trade and Industry in Japan announced their support for basic research of GICs. Since, at present, GIC and geoelectric field measurements only exist at a few locations in Japan, and the MT data and models only cover several smaller areas, some synthetic analysis is appropriate.

In congruence with the GIC discussions of much of Europe, especially the British Isles (Sect. 8.8), and of New Zealand (Sect. 8.15), coast effect is a primary consideration in Japan (e.g., see Uyeshima et al. 2001; Goto 2015, and Sect. 6). This raises some issues for 1D GIC modeling as in Sect. 6 (Pulkkinen et al. 2007) that involves the (a, b) parameter fitting. For example, Watari (2015) fit the (a, b) parameters to the GIC measurement in Hokkaido, Japan, and then apply these same parameters to geoelectric field measurements at Memanbetsu observatory, which is also located in Hokkaido, Japan, to estimate the GICs for the March 1989 and October 2003 geomagnetic storms. These were found to be logistically insignificant in magnitude. These results, however, cannot be considered conclusive since, as discussed in Sect. 6, the (a, b) parameter fitting approach is only good at the location where GIC data were previously acquired, and only for the same power grid system configuration. While these computations provide a good estimate for that transformer, applying these conclusions to a different location can be misleading in terms of the modeled GIC amplitudes. This is especially true for areas of significantly heterogeneous electrical conductivity.

For this reason, Fujita et al. (2018) address the nation’s GIC needs by creating an approximate 3D conductivity map of Japan and using that to model geoelectric fields in Japan in the frequency domain. They do so by utilizing the bathymetry maps from ETOPO1 (Amante and Eakins 2009) and a global digital map of sediment thickness (Laske and Masters 1997), interpolating these to accommodate the higher resolution of their model. They generate a 3D near-surface resistivity distribution for the uppermost 1 km model layer by assigning typical resistivity values to each of seawater, sediments, and rocks (Fig. 6, left). Even though, conceptually, Fujita et al. (2018) follow the thin-sheet approximation approach of the Püthe and Kuvshinov (2013a) global study, Thomson et al. (2005) in the UK, and Divett et al. (2017) in New Zealand, in practice the conductivity model they construct for Japan, even if approximate, is fully 3D, and therefore allows the authors to employ the full-physics approach of 3D EM.

A fully 3D, global geomagnetic field modeling code based on the staggered-grid finite-difference approach (Uyeshima and Schultz 2000) was used for the modeling. A $0.125^\circ \times 0.125^\circ$ lateral grid resolution was used in and around Japan, with the mesh getting gradually coarser away from the area. Ground geoelectric fields are then approximated from fully 3D Earth currents, which the code also computes. For the source, the approach of Püthe and Kuvshinov (2013a) (see Sect. 7) was followed, where the unit electric fields were computed for each spherical harmonic of the source, and a linear combination of these was taken to reflect realistic sources. However, in this particular study, the authors employed a simple plane-wave like equatorial current source known as the $P_1^0$; the source was rotated by an inclination angle $\alpha$ to model magnetic sources with arbitrary polarizations. The coast effect in geoelectric fields from the West to East equatorial ring current in Japan is evident in Fujita et al. (2018) their Figure 5, which shows the ratio of 3D modeled vs 1D modeled geoelectric fields for the same source current. Unfortunately, we disagree with the authors in that the comparison of modeled MT impedances (telluric vectors) with the measurements of Uyeshima et al. (2001) shows any obvious similarity. It is clear that smaller-scale structure is not yet captured by this approximate model, but there are possibly also other sources of error. Overall, higher resolution modeling better informed by data is needed. Full GIC modeling in Japan is the next step, if the power grid system parameters can be acquired or approximated (Nakamura et al. 2018).

For completeness, we would like to also touch upon several localized studies in Japan. Because of the availability of unusually high-quality, long-term geoelectric field measurements in Japan (Fujii et al. 2015), Japan’s geological complexity and ocean effects, and some GIC data availability, several international studies set up their geoelectric field modeling algorithms in Japan to enable validation. A study by Pulkkinen et al. (2010) used a plane-wave, 1D parameter fitting approach to reproduce measured GICs in Hokkaido, Japan (see Sect. 6). This approach is only applicable in locations where GIC data exist. The global study by Püthe et al. (2014) reproduced geoelectric field measurements at KAK, KNY and MMB geomagnetic observatories in Japan using rigorous global thin-sheet at $1^\circ \times 1^\circ$ resolution + 1D modeling, with local adjustments (see Sect. 7). This approach is only applicable where geoelectric field data exist, or in select 1D Earth regions such as on the seafloor. Finally, Kelbert et al. (2017) illustrated a time-domain approach for local geoelectric field modeling using an empirical 3D MT impedance, and reproducing the geoelectric field measurements in Japan at MMB to within 15%. This approach is only applicable where empirical 3D MT impedances have been obtained. It is also applicable to modeled 3D MT impedances, but that would require additional validation to ensure that the 3D model is accurate enough.

8.15 New Zealand

New Zealand has an exceptional network of GIC observations, collected by the nation’s power grid operator, Transpower, since 2001 in a near-continuous manner. As is clearly explained by Mac Manus et al. (2017), the data collection was initiated to assist with New Zealand’s comprehensive renewable energy program. In the South Island, which is blessed with an extensive river network, 98% of the electricity generation comes from hydroelectricity. On the other hand, during dry years, the South Island tends to benefit from the thermal generation capacity of the North. Since 75% of the population of New Zealand lives on the North Island, in 1965, a 611-km-long high-voltage direct current (HVDC) link was established between the two islands to deliver the electricity to where it is needed most (Fig. 8). The New Zealand HVDC system usually operates as a bipole with equal currents traveling on the transmission conductors between the two islands. However, it is also not uncommon for the system to operate in a monopolar mode, i.e., using a single wire and the ground. To enable this Earth return operation, grounded electrodes have been installed at the two ends of the HVDC link. In this situation, the return currents mostly travel through the Earth between the two electrodes, but a small fraction (called “stray currents”) travel instead to the earthing points of network transformers in much of the South Island and complete the circuit across the AC transmission lines.

The neutral current measurements, then, were installed to monitor these stray currents. In 2001, monitors at 12 different substations were installed, often monitoring multiple transformers within a substation. In 2009, the network was expanded to 17 substations. Overall between 2001 and 2015, up to 58 individual transformers were simultaneously monitored to measure the stray current returns. However, since the magnitudes of these currents are related to the total HVDC Earth return current in a linear fashion, Mac Manus et al. (2017) were able to correct these long-term measurements for network originating signals, resulting in a comprehensive GIC database.

Fig. 8

On the left, a map of the South Island of New Zealand showing the Transpower New Zealand electrical transmission network, with different colors representing different line voltages. The heavy purple line is the HVDC transmission line linking the South Island and North Island electrical networks. On the right, an illustration of the mechanism through which the HVDC return current travels through the Earth between the two electrodes (red triangles). A small fraction of the current finds paths traveling from the northern electrode to the earthing points of network transformers in much of the South Island and then completing the circuit across the AC transmission lines (“stray currents”). Transpower has installed DC current measuring devices on all the transformers for which significant stray HVDC return current is expected. Mac Manus et al. (2017, Reuse permission purchased from Wiley) their Figs. 1 and 2

Initial analysis of these data by Marshall et al. (2012) involved an illustration that these were, indeed, natural signals associated with space weather. The authors convolved geomagnetic field measurements with an impedance from a uniform half-space unity conductivity model, resulting in a proxy for the geoelectric field they called a GIC index. While the true geoelectric field would be a linear combination of the two GIC indices, and the actual GIC would depend not only on the regional geoelectric field but also on the system configurations, they were able to use the index to show qualitative, highly convincing correlation between the (proxy) natural geoelectric field signal, and the corrected neutral current measurements in the system. We should also like to note that recent GIC measurements from this network for the September 2017 geomagnetic storm exceed 45 A at certain locations (Clilverd et al. 2018), and that an accurate transfer function based estimate for what the GICs in New Zealand would have been during the Carrington event of 1859 nearly reaches the value of 1000 A (Ingham et al. 2017). While the extensive monitoring network is invaluable, accurate modeling of GICs in New Zealand is the path forward for geomagnetic hazard mitigation.

Not unlike the British Isles and Japan, New Zealand is a moderately sized island, which of course means a land mass surrounded by seawater, exposing any GIC studies in New Zealand to the challenges of the coast effect (Sect. 6). This is best captured through a 3D Earth conductivity model. Short of that, a thin-sheet conductance map approach has previously been developed for the UK (e.g., Thomson et al. 2005; Beggan et al. 2013) and has now been adopted in New Zealand by Divett et al. (2017). A synthetic geologically inspired thin-sheet conductance model has been developed for both islands of New Zealand, including the ocean bathymetry and some on-land variations, which are primarily aligned with the main southwest to the northeast axis of New Zealand’s mountain backbone. The 2D conductance map lies on top of a three-layer 1D conductivity model. Justification for the validity of this approximation is given by the authors, although a more careful analysis of skin depth variability in the area might be warranted. A comparison of 34 modeled and measured (Chamalaun and McKnight 1993) induction vectors was performed over the islands of New Zealand, and adjustments to the conductance map were made, resulting in the model and data fits as shown in Divett et al. (2017) their Figure 3. Indeed, the modeled induction vectors are mostly of comparable size to those measured, with the average offset in their direction of about $20^\circ$.

This model was then used to compute estimated geoelectric fields for a variety of orientations of synthetic (plane-wave) horizontal geomagnetic field variations, and then combined with a well-known power grid system model, provided by Transpower as a snapshot from 2015 configuration, to estimate substation-level GICs. Note that of the 63 substations on the high-voltage network, Divett et al. (2017) showed that the highest geoelectric fields, and the highest GICs, were observed when the magnetic field was oriented parallel to the southwest to the northeast axis of the islands. This makes perfect sense, since this orientation of the inducing field would correspond to primarily southeast to northwest orientation of the Earth currents, precisely the direction in which the Earth currents could not flow easily, causing buildup of galvanic charges in the across-island (TM) mode, and would therefore more readily enter the power grid AC network. Striking features of the calculated GIC distribution (Divett et al. 2017 their Figure 8) included the high spatial variability in the GIC magnitude, and the small number of affected substations, quite in a similar manner to those observed in the UK (Beggan et al. 2013 their Figure 8). A related analysis by Divett et al. (2018) expanded this effort to consider transformer-level (rather than substation-level) GICs, thus eliminating one more uncertainty in the interpretation of real, transformer-level GIC measurements.

The wealth of multifaceted ground measurements in New Zealand, including the induction vectors from a magnetometer array covering all of New Zealand (Chamalaun and McKnight 1993), multiple regional MT surveys discussed in Divett et al. (2017), and of course the GIC observations mentioned above, is unique on the global stage of space weather hazard mitigation efforts. In conjunction with the manageable size of New Zealand’s network and modeling domain, these measurements make it an ideal natural laboratory for examining the uses of Earth conductivity modeling techniques in GIC estimation. Indeed, making full use of the modern computational capabilities, a technique as straightforward as a Monte Carlo inversion would enable dramatic improvement in the thin-sheet conductance + 1D model of New Zealand, while also exploring the limits of this approximation.

The strong collaboration between New Zealand’s primary power provider and academia, therefore, allows New Zealand’s researchers to explore GIC measurements in unique ways, and use their access to the neutral line current measurements to develop accurate and efficient GIC estimation techniques that would inform the power grid of ongoing space weather events. While this is not yet a solved problem, New Zealand’s space weather program has dramatically advanced in recent years compared to most of the rest of the world thanks to this unique data availability paradigm.

8.16 South Africa

Ever since transformer damages were reported on the South African power grid (e.g., Kappenman 2005), multiple studies addressed the GIC estimation problem in South Africa, all focused on the 1D Earth approximation approach to date. For example, Ngwira et al. (2009) report a comparison of modeled and measured GICs at one location for the Halloween storm, using three different streams of magnetic observatory data. The modeled GICs are obtained using a 1D electrical conductivity model at a single location and the (a, b) 1D method (Sect. 6). These improved (a, b) coefficients and the 1D electrical conductivity model have been obtained by Ngwira et al. (2008) by least-squares fitting to an existing GIC time series which, as accurately noted by the authors, are likely only valid for the specific magnetic observatory and GIC station pair and do not otherwise reflect the system configuration and the true conductivity structure of the Earth, respectively. Depending on the magnetic field data input, the modeled GICs in the Ngwira et al. (2009) study either fit the observations well at about 80% of the signal, or recover only about 20% of the variation. The authors also suggest that, confirming the results of Bernhardi et al. (2008), the geoelectric fields in southern Africa are indeed spatially variable and need to be modeled as such.

A less straightforward approach to GIC modeling in South Africa has been applied by Matandirotya et al. (2015). There, the authors employ 2D finite element modeling of the Earth. Unfortunately, this powerful technique was applied to the usual 1D model of the Earth and supplemented by the usual (a, b) parameter estimation method, although the actual method by which the parameters were estimated was modified. Comparison against measured GIC data at a single station is highly satisfactory, as should be expected when a version of parameter estimation is utilized.

Comprehensive MT data coverage (Jones et al. 2003–2008, 2009; Evans et al. 2011) in South Africa makes more advanced methods, see e.g., Sect. 5, directly applicable in this region.

8.17 USA

Academic research on geomagnetic hazards is flourishing in the USA. Thanks in no small part to the nationwide MT survey, Schultz et al. (2006–2018), substantial progress has been made on the methods of geoelectric field estimation (e.g., Wei et al. 2013; Bonner and Schultz 2017; Kelbert et al. 2017; Weigel 2017), geomagnetic hazard maps (Love et al. 2016, 2018a, b, c, analysis of power grid voltages using alternative electrical conductivity approaches (Lucas et al. 2018), and electrical conductivity models (Kelbert et al. 2019), to name only those areas of research that are of direct relevance to the topic of this review. In fact, mitigation of space weather hazards is a national priority for the USA. In the USA, geomagnetic disturbances have been officially recognized at a threat to national infrastructure by the Federal Energy Reliability Commission (FERC 2016), Federal Emergency Management Agency (FEMA) (MacAlester and Murtagh 2014), and the White House (e.g., NSTC 2015, Goals 1.1 and 5.5). Along with the US Air Force, the Space Weather Prediction Center (SWPC) of the National Oceanic and Atmospheric Administration (NOAA) is the primary agency that is responsible for monitoring and reporting on space weather activities. NOAA has also developed a Space Weather Scale based on severity levels similar to hurricanes. Geomagnetic Storms are rated G1 (Minor) to G5 (Extreme), Solar Radiation storms are rated S1 to S5, and Radio Blackouts are rated R1 to R5. In addition to this crude warning indicator, NOAA is working with the U.S. Geological Survey (USGS) to develop detailed real-time geoelectric field maps for the contiguous USA. The initial operational product (released in September 2019) that is based on 1D regions of Fernberg (2012) will soon be supplemented with a corresponding geoelectric field map product based on the USArray MT data (Schultz et al. 2006–2018) directly, and/or with an alternative version based on the conductivity model compilation by Kelbert et al. (2019).

Full-physics modeling, a prerequisite for actionable GIC prediction, is another active area of research in the USA. An extensive selection of space weather modeling codes is hosted at NASA’s Community Coordinated Modeling Center (CCMC). One of these community models, a global Geospace model (e.g., Welling and Ridley 2010) based on the Space Weather Modeling Framework (SWMF), also runs operationally at NOAA SWPC in the USA, but this model does not yet reproduce the data well enough to be useful for practical geomagnetic hazard mitigation. Several projects funded by the NSF’s Prediction of and Resilience Against Extreme Events (PREEVENTS) program in 2017, aim to improve these modeling capabilities. One of these funded projects, the PREEVENTS Comprehensive Hazard Analysis for Resilience to Geomagnetic Extreme Disturbances (CHARGED), aims to improve (Welling et al. 2017) SWMF with realistic particle precipitation and ionospheric conductance physics, and to integrate it with a global 3D finite-difference time-domain (FDTD) Earth model (Simpson 2009). Collaboration with the USGS has been established to incorporate up-to-date 3D Earth conductivity models (Kelbert et al. 2019) (Fig. 9) into the Earth modeling to estimate geoelectric fields in the time domain. The other funded PREEVENTS project, under the name of Integrated Modeling of Extreme Space Weather Events from Electron to Global Scales, is focused on integrating magnetohydrodynamic modeling such as Geospace with small-scale kinetic modeling that accounts for magnetic field reconnection. Both of these projects address the development of better Earth system modeling for magnetic fields, which is a critical stepping stone toward forecasting. However, these projects are not yet intended to provide operational and predictive space weather capabilities. The latter would involve real-time and multifaceted data assimilation. Although focused ionospheric data assimilation codes exist (e.g., Richmond et al. 1992; Bust et al. 2004), no prior or current work on Sun-to-Earth modeling has yet reached that stage of maturity.

One aspect where the USA is lacking is validation against GIC measurements, which are rarely available due to sensitivity and proprietary concerns. Even though other countries have had successful collaborations between the academic researchers, the government, and the power grid industry on geomagnetically induced current estimation and validation, the USA unfortunately has a historical tradition of maintaining all of their measurements, including the transformer neutral currents and system configurations, highly sensitive and proprietary. Development of a useful hazard mitigation system would be much accelerated were these data more easily accessible. Progress is being gradually made on that front, but, at present, the academic and government researchers are limited in the validation procedures they can accomplish. A notable recent exception to this general rule has resulted in an important government and industry collaboration (Sun and Balch 2019), which has enabled a rigorous validation of alternative 1D and 3D approaches to geoelectric field computation against real GIC data at least in one area of the USA (Virginia and the North and South Carolinas). Based on the three available GIC time series collected during the September 2017 geomagnetic storm, this study concluded that MT impedances, smoothed and interpolated to a $0.5^{\circ }$ regular grid, produce substantially improved estimates of GIC amplitudes compared to the alternative approaches based on 1D and brute-force local 3D MT impedances. Several additional studies have had access to GIC data for validation, but these have, to date, been preliminary and rather inconclusive (e.g., Butala et al. 2017; Shetye et al. 2018). An ongoing NSF-funded smart GIC-grid project, led by Oregon State University, that aims to bring real-time GIC estimates into the operating rooms of nationwide power grids, might help to further bridge this gap.

Fig. 9

A depth slice at 30 km depth from a fully 3D compilation of electrical conductivity models in North America, obtained by multiple authors using USArray MT (Schultz et al. 2006–2018) and other magnetotelluric data in the USA and ModEM MT inversion software (Kelbert et al. 2014). BREP, Basin and Range extensional province; CAC, Cascades arc conductor; GS, Grenville suture; KBML, Klamath-Blue Mountains lineament; NACP, North American Central Plains anomaly; PT, Piedmont terrane (Kelbert et al. 2019)

8.18 Uruguay

Significant progress on modeling the GICs in the Uruguayan high-voltage power grid has been made by Caraballo et al. (2013) and Caraballo (2016). The former of these studies sets up the comprehensive power grid modeling for Uruguay, using the Lehtinen and Pirjola (1985) approach and integration of electric fields along power lines, as well as SECS (Amm and Viljanen 1999) for geomagnetic field interpolation. The follow-up work, Caraballo (2016), employs this setup to estimate GICs in the power grid. Even though the authors refer to a 2D conductivity model, this appears to be a compilation of two 1D profiles, one 4-layer and overall more conductive, the other 3-layer and overall more resistive, which are used to represent Uruguay’s electrical conductivity. Under this approximation, GICs are estimated using the Lehtinen and Pirjola (1985) approach in the 500 kV Uruguayan power network. Under certain additional assumptions on the line resistances, as well as the grounding resistances of the network, modeled GICs at two substations are obtained under different system configurations and conductivity model assumptions, and compared.

9 Conclusions

This paper aims to provide a comprehensive review on the status of international natural geomagnetic hazard mitigation efforts. As we have seen, hazard mitigation efforts include gaining a quantitative understanding of GIC generation in the affected technological systems and taking measures to reduce the related risk, either in a preventative fashion, or in real time during a geomagnetic disturbance. Geomagnetic hazards are similar to most other natural hazards, such as earthquakes, volcanic eruptions, and landslides, in that mitigation strategies and system hardening have a critical up-front component; however, prediction of such events is extremely hard. To date, most efforts have been directed at investigation of past events in order to increase the robustness of the infrastructure, and (most recently) on real-time geoelectric field estimation that would provide situational awareness, for example, to the power grid operators in a control room. These and other efforts would benefit from incorporating the realistic Earth’s electrical conductivity into their procedures. Improvements along these lines are an active area of research in the international space weather and magnetotelluric communities.

While we focus on the role of the solid Earth for the purposes of this review, we should explicitly note the inherent complexity of the natural geomagnetic hazard mitigation problem, which extends far beyond geophysics. Quantifying GIC hazards requires a highly multidisciplinary approach, as it follows the physics and the observations of solar wind through the magnetosphere, its interaction with the ionosphere, the electromagnetic induction causing currents to flow in the Earth and, finally, the subtle engineering details of the power grid configuration. All these steps are necessary to provide an accurate GIC estimate in a power line or at a transformer grounding. To provide a real-time estimate or an actionable prediction of this value, much of this complexity must be considered. While this is a challenge, we also invite the electromagnetic geophysics community to view this complexity as an opportunity to expand our collaborations beyond the domain of geophysics: a move that will, no doubt, also benefit basic science of the Earth system in the long term.

Appendix A

Equation (4) can also be expressed as

$$\begin{aligned} I_j(t) = \left[ \sum _{i} \sum _{k\ne i} \left[ \frac{K_{ji}}{R_{ki}}\right] \tau ^N_{ki}(t) * B^N_{ki}(t)\right] + \left[ \sum _{i} \sum _{k\ne i} \left[ \frac{K_{ji}}{R_{ki}}\right] \tau ^E_{ki}(t) * B^E_{ki}(t)\right] , \end{aligned}$$

(7)

where $B^N_{ki}(t)$ and $B^E_{ki}(t)$ are the northward and eastward magnetic fields averaged over each of the transmission lines, $*$ stands for the convolution operator, and

$$\begin{aligned} \tau ^N_{ki} = \int _{{\mathbf {x}}_k}^{{\mathbf {x}}_i} Z_{xx} \mathrm{d}x + \int _{{\mathbf {x}}_k}^{{\mathbf {x}}_i} Z_{yx} \mathrm{d}y = \int _L {\mathbf {Z}}_x \cdot {\mathbf {dl}}, \end{aligned}$$

(8)

and

$$\begin{aligned} \tau ^E_{ki} = \int _{{\mathbf {x}}_k}^{{\mathbf {x}}_i} Z_{xy} \mathrm{d}x + \int _{{\mathbf {x}}_k}^{{\mathbf {x}}_i} Z_{yy} \mathrm{d}y = \int _L {\mathbf {Z}}_y \cdot {\mathbf {dl}} \end{aligned}$$

(9)

are the northward and eastward, frequency $\omega$ (or time t)-dependent power-line telluric response functions, obtained by integration of the local MT impedance components over the electric power-line paths. Here, L defines the path of the power line, and the expressions ${\mathbf {Z}}_x = \left[ \begin{array}{cc} Z_{xx}&Z_{yx} \end{array}\right]$ and ${\mathbf {Z}}_y = \left[ \begin{array}{cc} Z_{xy}&Z_{yy} \end{array}\right]$ are space and frequency $\omega$ (or time t) -dependent telluric vectors (Bahr 1988).

If the magnetic field is assumed to be locally uniform in the area of the power grid of direct relevance to transformer j, Eq. 7 is reduced to

$$\begin{aligned} I_j(t) = \left[ \sum _{i} \sum _{k\ne i} \left[ \frac{K_{ji}}{R_{ki}}\right] \tau ^N_{ki}(t)\right] * B_x(t) + \left[ \sum _{i} \sum _{k\ne i} \left[ \frac{K_{ji}}{R_{ki}}\right] \tau ^E_{ki}(t)\right] * B_y(t), \end{aligned}$$

(10)

an expression which, in spite of its relative complexity, can be recognized as a straightforward linear combination of the two local horizontal magnetic field components. Expression (11) easily follows from Eq. (10) under the 1D Earth assumption:

$$\begin{aligned} I_j(t) = - \left[ \sum _{i} \sum _{k\ne i} \frac{K_{ji} L^y_{ki}}{R_{ki}} \right] Z^1(t) * B_x(t) + \left[ \sum _{i} \sum _{k\ne i} \frac{K_{ji} L^x_{ki}}{R_{ki}}\right] Z^1(t) * B_y(t), \end{aligned}$$

(11)

where $L^x_{ki} = \int _{{\mathbf {x}}_k}^{{\mathbf {x}}_i} \mathrm{d}x$, $L^y_{ki} = \int _{{\mathbf {x}}_k}^{{\mathbf {x}}_i} \mathrm{d}y$ are the line integrals, $Z^1(t)$ is the impulse response representation of the 1D Earth impedance. Of course, in the 1D Earth approximation, $E_x(t) = Z^1(t) * B_y(t)$ and $E_y(t) = - Z^1(t) * B_x(t)$. Thus, Eq. 11 is just another representation of the 1D GIC formulation Eq. 5 of Sect. 6.

Note that the fortunate, but unlikely, turn of events when a 1D approximation is sufficient corresponds to a situation in which the $\tau ^N_{ki}$ and $\tau ^E_{ki}$ coefficients of Eq. 10 are exactly equal to their 1D counterparts in Eq. 11.