Modified GIC Estimation Using 3‐D Earth Conductivity

Geomagnetically induced currents (GICs) are quasi‐direct current (DC) electric currents that flow in technological conductors during geomagnetic storms. Extreme GICs are hazardous to man‐made infrastructure. GICs enter and exit the technological systems, such as the electric power grid, at grounding points, and their magnitudes depend on the currents that flow underground. They are, therefore, a function of the Earth's electrical conductivity, represented at ground level as Earth impedances, as well as the resistance parameters of the power network. Traditional GIC estimation practices are based on Earth impedances obtained from laterally homogeneous or piecewise layered‐Earth models. We refer to these methods, collectively, as the 1‐D approximation. However, GIC hazard mitigation can be improved with more accurate GIC modeling that takes the spatially heterogeneous Earth's conductivity into account. Here, we propose a modified approximation for GIC estimation that is very similar to the 1‐D approximation but is instead derived from empirical 3‐D Earth impedances. Our formulation sets up the computation of static, frequency‐dependent power line telluric response functions, which, once computed, may be considered part of the power grid system model. These response functions may then be used for historical scenario analysis of GIC hazards and for simplified real‐time, albeit approximate, GIC estimation in a power grid. This modest modification to the simpler local field formulation approach avoids real‐time integration of geoelectric fields along power lines while taking the realistic 3‐D Earth into account in a rigorous manner. Once implemented, the method provides a power grid operator with the benefits of convenience and computational speed for a first look real‐time operational GIC hazard assessment. We estimate that the proposed modified 3‐D GIC modeling approach produces GIC values that are well within 50% of those obtained with the full‐scale power line integration of spatially variable geoelectric fields, for storms comparable in scale to the 2003 Halloween storm, all geological structures, and power lines located in the contiguous United States and other low‐ to middle‐latitude regions.


Introduction
During geomagnetic storms, electric currents are induced in the conducting Earth. These currents propagate into power lines through the transformers that are grounded, causing quasi-direct current (DC) known as geomagnetically induced currents (GICs; e.g., Barnes et al., 1991;Lanzerotti, 1983;Love et al., 2014;Trichtchenko, 2016). As they pass through transformer windings, these currents can cause transformer half-cycle saturation, which results in adverse effects to the bulk power system, including harmonic distortion and reactive power loss. Historically, GICs have caused transformer melting (e.g., Gaunt & Coetzee, 2007) and, in several cases, a complete loss of power, such as the famous collapse of the Hydro-Québec power system in Canada during the geomagnetic storm of 13 March 1989, which left the six million residents of Québec without power for over 9 hr (e.g., Allen et al., 1989;Bolduc, 2002;Boteler, 2001Boteler, , 2019Guillon et al., 2016). Internationally, multiple other electric power transmission failures have occurred as a result of geomagnetic storms and the associated GICs (see Eastwood et al., 2017, for an overview). In the United States, the North American Electric Reliability Corporation (NERC) has developed reliability standards to address the impact of geomagnetic disturbances (GMD) on the bulk power system. Different standards have been developed for planning and for operational purposes, whereby more detailed and accurate computations are required for planning than for real-time mitigation of GIC effects. Most recently, the U.S. Federal Energy Regulatory Commission (FERC, 2018) (Order No. 851) has approved the NERC's reliability standard TPL-007-2 and the associated GMD research work plan. While the current regulatory requirements for reliability assessments only include the use of a single spatially averaged reference peak geoelectric field value of 8 V/km and another nonspatially averaged peak geoelectric field value of 12 V/km, both adjusted with regional scaling factors, much more detailed analysis is included in the GMD research work plan.
This includes the development of standardized methods for nonuniform geoelectric field modeling and promoting the adoption of such techniques in commercially available software tools (Task 3C).
Specifically, accurate methods for substation-level GIC estimation in a power grid network (e.g., Boteler & Pirjola, 2017;Overbye et al., 2013;Viljanen & Pirjola, 1994) require that ground-level geoelectric fields are integrated along power lines, to obtain line voltages. These voltages are then employed to model the induced neutral currents at transformer groundings, located at the buses and substations in the network. This is the best available proxy to transformer-level GIC measurements for most power grid systems (e.g., Divett et al., 2018). Such computations could be performed in real time, resulting in a highly accurate estimate of GIC hazard. However, due to the relative complexity of the calculation and the historically limited availability of real-time geoelectric field measurements, such calculations have only been performed for planning purposes. Such GIC risk assessment analysis based on hypothetical or historical storms has been implemented both by the power grid utilities in-house with the use of their proprietary system grid configuration models using generic, generally commercial, power grid modeling software, and in the academic literature, using approximate grid system configurations and, generally, homemade or open source modeling packages. These studies achieve different levels of modeling complexity, ranging from the use of a half-space or 1-D Earth assumption (e.g., Viljanen et al., 2014, in Europe;de Silva Barbosa et al., 2015, in Brazil;and Liu, Ganebo, et al., 2018;Senbato et al., 2018, in Ethiopia), to the use of a compilation of 1-D Earth models (e.g., Liu et al., 2014, in China;Marti et al., 2014, in Canada;Caraballo, 2016, in Uruguay;Blake et al., 2016, in Ireland;Kelly et al., 2017, in the United Kingdom and France;and Espinosa et al., 2019, in Brazil), to the laterally heterogeneous thin-sheet overlaying a 1-D profiles (Beggan, 2015, in the United Kingdom; Divett et al., 2017, in New Zealand;Nakamura et al., 2018, in Japan;and Bailey et al., 2018, in Austria), and 3-D Earth models (Rosenqvist & Hall, 2019, in Sweden;Liu, Wang, et al., 2018, in China;and Marshall et al., 2019, in Australia), the latter all recent developments (see Kelbert, 2020, for a detailed discussion). Here, a 1-D Earth refers to the assumption that the Earth's conductivity is laterally homogeneous and varies only with depth, whereas 3-D Earth models make no such assumption. These Earth conductivity models are then combined with GIC modeling on the full power grid system, to the extent of an investigator's knowledge of the network, to estimate historical or hypothetical GICs in the electric power lines. All of these approaches, in principle, allow for computation of laterally heterogeneous geoelectric fields that reflect the Earth's geology. However, the need to update spatially variable geoelectric fields, perform the integration along power lines to obtain line voltages, and to incorporate these in the power grid system model operationally in real time has resulted, to date, in very limited adoption of realistic Earth conductivity models for real-time GIC hazard monitoring in the industry.
Other methods based on empirical magnetotelluric (MT) impedances (e.g., Bonner & Schultz, 2017;Lucas et al., 2020), combined with gridded national MT surveys such as EarthScope USArray MT in the United States (Schultz et al., 2006) also compute the local geoelectric fields on a national scale with high accuracy. As is shown in Sun and Balch (2019) for one geologically complex area of the eastern United States (Virginia and the North and South Carolinas), MT impedances, smoothed and interpolated to a 0.5°regular grid, produce substantially improved estimates of GIC amplitudes compared to the alternative approaches based on 1-D and brute-force local 3-D MT impedances. All of these approaches, applied operationally, require interpolated geomagnetic field measurements to be convolved with the MT impedances in real time to obtain geoelectric fields, which are then integrated along power lines to inform the power grid system model.
On the opposite end of the complexity spectrum, simplified methods have been developed under the uniform geoelectric field assumption (Pulkkinen et al., 2007;Viljanen et al., 2006) based on a linear relationship between the local GIC and the local geoelectric field. Some applications of the approach include Ngwira et al. (2008), Bernhardi et al. (2008, and Matandirotya et al. (2015) in South Africa, Liu et al. (2009), Liu et al. (2013), and Zhang et al. (2015 in China, Pulkkinen et al. (2010) and Watari (2015) in Japan, and Wik et al. (2008) and Viljanen et al. (2012Viljanen et al. ( , 2013 in Europe. This approximation takes a local geoelectric field, computed from the local geomagnetic field, which is convolved with a 1-D model of the Earth, as input. It does not require a full system model and is instead based on empirical or calculated fitting parameters. This method yields itself readily to efficient real-time computations of GIC hazard but compromises on the accuracy of the modeling. Specifically, the 1-D Earth assumption that underlies it is only valid in a very limited number of geological scenarios (Kelbert, 2020), and the fitting parameters are only applicable at the location where they were first estimated and only until the system configuration changes.
Here, we propose a straightforward mathematical formulation that bridges the gap between the complex, multistage realistic GIC computation, and the simple linear approximation that assumes a uniform geoelectric field. Our proposed method makes full use of the 3-D MT impedances, where these measurements are available, and incorporates them as static, frequency-dependent coefficients into the power system grid configuration. This approximation provides a methodologically simple and streamlined operational GIC hazard assessment in areas with geologic complexity, for low-to middle-latitude areas where the geomagnetic field can be assumed, to first order, locally spatially homogeneous even during a magnetic storm. We set the stage for the modified GIC computation based on the commonly used Lehtinen-Pirjola method for GIC calculation in a power grid (Lehtinen & Pirjola, 1985) in section 2. While alternative formulations exist, such as the Nodal Admittance Matrix method often favored by electrical engineers, these have been shown to be mathematically equivalent (Boteler & Pirjola, 2014). Indeed, whether written in terms of the line voltages, discussed in the Lehtinen-Pirjola approach, or the nodal voltages, as in the Nodal Admittance Matrix method, both approaches are conceptually identical in how the geoelectric fields are incorporated into the computations. In section 3, we propose a modification to the Lehtinen-Pirjola algorithm which, under the locally uniform geomagnetic field approximation, allows empirical or modeled 3-D Earth impedances obtained with MT methods to be utilized when setting up the time-stationary power grid system coefficient matrices. We also discuss how our formulation may be reduced to the traditional frequency-independent local power grid approximation of, for example, Pulkkinen et al. (2007) if the Earth conductivity is assumed to be laterally homogeneous. In section 4, we provide a simple illustration of how the modified, frequency-dependent power grid system coefficients that reflect the Earth's heterogeneous 3-D electrical conductivity distribution may be derived and used in conjunction with the local geomagnetic field to improve the accuracy of GIC estimation compared with its 1-D counterpart with only a minor compromise on the algorithm's simplicity and touch upon the implications of our findings in section 5. We discuss the overall applicability and limitations of the new formulation and provide a broader context in section 6.

Review of GIC Calculations in an Electric Power Grid
GICs can flow everywhere in the power grids, entering the grid at the points of transformer groundings. 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: where 1 is the identity matrix, Y n and Z e are the network admittance and the earthing impedance matrices, respectively, and 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 Kirchoff's law, its element J i at an earthing point x i is defined as the sum of all currents I ki flowing into the node from the power grid. Using Ohm's law, where R ki are all of the power line resistances for lines, with the endpoints at x k , connected to the earthing point x i , and V ki are the induced voltages in the respective power lines. These voltages are computed as the line integrals of geoelectric field along the power line path, 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 E is not, in general, conservative.
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, I e as in Equation 1. 10.1029/2020SW002467

Space Weather
Since matrix Y n is fully defined by the resistances in the power lines, R ki , and matrix 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, its spatial structure, and the spatial structure of geoelectric field.
We would also like to note explicitly that at an arbitrary transformer j, the GIC I j (t) may be mathematically expressed as a linear combination of the line voltages in the system, V ki (t). If we define the matrix K e = (1+Y n Z e ) −1 = [K ji ], In conjunction with Equation 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. For a constant system configuration, the system parameters K ji , defined between any two grounding connections, such as transformers, and the line resistances R ki may be considered constant in time: The time dependence of a GIC comes from time variation in line voltages. From this point on, we will explicitly indicate time (or frequency) dependence, as well as continuous spatial dependence of quantities, as this is relevant to our discussion. Transformer and power line dependence will be indicated by subscripts, as before.
Here, we focus on the computation of geoelectric fields, as that is the component that requires geophysical knowledge. However, it is notable that the grounding resistances also depend upon the geology of the area. 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 (Kelbert, 2020). Since the transformer grounding systems are engineered to be within specified range of resistance so that they will safely earth the normal electrical load of the transformer, the grounding resistance is a function of both local geology and the properties of the earthing mat at the substation. It also varies over time due to corrosion, groundwater, and other time-varying factors: a complication that could be accounted for if the relevant measurements were available.
If the geoelectric field can be assumed laterally homogeneous, a simplified GIC estimation procedure holds (Pulkkinen et al., 2010;Viljanen & Pirjola, 1994). At a single transformer j, the GIC value at any point in time is where a j and b j are constants that are characteristic for each transformer 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 magnetic field produced by ionospheric currents can be approximated as a spatially uniform plane wave and the electrical conductivity of the Earth varies only with depth (1-D). An implicit assumption required for validity of this approximation is that the same 1-D Earth conductivity profile applies everywhere beneath the footprint of the power grid.

An Expression for a General GIC Calculation Using 3-D MT Impedances
As we recall from section 2, 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 Equation 2, before Equation 1 is applied. All of these are linear operations and can be reduced to Equation 5 if the geoelectric field E is uniform over the footprint of all power lines connected to transformer j. Let us, however, assume that the geoelectric field is nonuniform over the area but that the magnetic field can be assumed spatially (but not temporally) uniform over the footprint of the power grid. We shall rely on this assumption for our derivation, but will revisit it later in the discussion.
Based on Maxwell's equations, both the magnetic and electric fields depend on the Earth's electrical conductivity in a nonlinear fashion. However, we also note that the electric field is linearly related to the magnetic field at the same location. This latter relationship is traditionally expressed via the general MT response tensor, which relates the horizontal components of the inducing geomagnetic field to the horizontal components of the induced geoelectric field (e.g., Weidelt & Chave, 2012), Eðω; xÞ ¼ Zðω; xÞ Bðω; xÞ; where defines the MT response tensor as a function of angular frequency ω and location on the Earth's surface, x.
The MT response tensor tends to have the practical units of (mV/km)/nT and is related to the MT impedance μZ (measured in ohms), where μ is the magnetic permeability, via a simple rescaling of units. This full impedance tensor reflects the three-dimensional structure of the Earth, in the physical sense.
Since the voltage V ki (t) of Equation 4 is being computed at an instant in time, we prefer to allow the MT response tensor Z to assume a time domain impulse response formulation, as in Kelbert et al. (2017). By the convolution theorem (e.g., Bracewell, 2000, Chapter 3), at each location x along each of the power lines, the geoelectric field can be expressed as where Z(t,x) is the fully 3-D local Earth impedance, expressed as a continuous impulse response, and B(t, x) is the local horizontal magnetic field. Note that while Equation 8 is strictly causal, its discrete time counterpart is generally acausal (Egbert, 1992;Kelbert et al., 2017), such that Z(t,x)≠0 for t<0. However, the negative time contribution is small for practical purposes.
Under the uniform magnetic field assumption, B is no longer dependent on location x, and Equation 8 reduces to If the resistances of the power grid network components are considered constant, the local GIC measurement I j at transformer j is a linear combination of all connected line voltages in the system, which may be expressed as a convolution (Appendix A): or, equivalently, a multiplication in the frequency domain, where we define and to be the northward and eastward, frequency ω-dependent (or time t-dependent) power line telluric response functions. Here, L defines the path of the power line, and the expressions Z

Space Weather
are space and frequency ω-dependent (or time t-dependent) telluric vectors (Bahr, 1988). Example amplitudes for τ N and τ E for the high-voltage power grid lines in the United States are illustrated in Figure 1.
Here, a frequency domain representation of τ N and τ E is obtained at the frequencies at which the empirical 3-D MT impedances are available, specifically at 31 periods spanning approximately 7 to 30,000 s. The power line integrations in Equations 12 and 13 are done through a Delaunay triangulation method that interpolates between three nearby MT sites to project the impedance tensor components onto the power line, discretized as 1-km segments (Bonner & Schultz, 2017;Lucas et al., 2018). These projected impedance tensor components are then integrated to obtain the telluric vectors shown in Figure 1. As can be seen from Equations 12 and 13 above, power line telluric mapping functions are defined by the geometry of the power grid and by the 3-D Earth impedances in the area and can be viewed as stationary, frequency-dependent power grid system coefficients that are independent of the time-varying geomagnetic fields.

Substituting this expression in Equation 4
, we obtain a general 3-D expression for the GIC at transformer j: or, in the frequency domain, replacing the convolution operator with simple multiplication,

Space Weather
Note that the power line telluric mapping functions τ N ki ðωÞ and τ E ki ðωÞ (or their impulse response counterparts τ N ki ðtÞ and τ E ki ðtÞ) may be precomputed for each power line and saved as part of the system configurations. While the frequency domain representation of the power line telluric mapping functions comes out of the 3-D MT impedances most naturally, their impulse response representation may be quite easily obtained using the procedure described in Kelbert et al. (2017). They can then be multiplied (or convolved) with the local magnetic field at each transformer in real time, resulting in a highly efficient and intuitive algorithm for operational GIC hazard analysis. The only conceptual modification from preexisting approaches based on the 1-D Earth approximation would be related to the frequency dependence (or impulse-response form) of the new system parameters.
Then, in the spatial sense (though not in the time sense), the GIC at a transformer location can be considered a linear combination of the horizontal components of B in the presence of arbitrary 3-D geometry, under the assumption that the magnetic field is uniform over the area centered at the transformer location. That is, in fact, typically a reasonable assumption over a small section of power grid, connected to a single transformer. Indeed, magnetic fields at the Earth's surface are typically coherent over several hundreds of km spatial scales (e.g., Watermann et al., 2006), which is comparable with the length of a typical high-voltage power line in the United States. If the 1-D approximation of the Earth's conductivity structure is enforced, this general derivation reduces to the classical frequency-independent power grid system coefficients of Equation 5 (see Appendix B for the exact form of these coefficients). However, if the power line telluric response coefficients of Figure 1 are directly compared to those obtained under several reasonable 1-D approximations, we note that in certain geologically complex areas, the differences are comparable to the amplitudes of these coefficients ( Figures C1-C3). In all cases, the 3-D effects are more pronounced at shorter periods of tens to hundreds of seconds that are well represented in the signal power of geomagnetic storm variations of relevance to GIC estimation. This illustrates the relative importance of using the full impedance tensor versus a scalar impedance representation for line voltage computations.

Usage and Limitations of the Modified Algorithm
Here, we present a simple thought experiment inspired by Rosenqvist and Hall (2019), who short circuited a power transmission line in Sweden to Earth so that it operated in a stand-alone mode, disconnected from the rest of the power grid system. This enabled a very straightforward validation against a measured line GIC. In this simplified setting, a neutral ground GIC measurement could be expressed as where ground is the total impedance of the transmission line system, including the earthing and line resistances, measured in ohms. The inductive contribution to the transmission line system impedance may be ignored at the low frequencies relevant to this study. Rosenqvist and Hall (2019) set each of the ground and line resistances to 1 Ω.
To illustrate the utility and accuracy of our method, we analyze the estimated GIC in several hypothetical transmission lines shown in Figure 2. The hypothetical line locations were chosen to illustrate the effects of the Earth's heterogeneous electrical conductivity structure, while also being at high enough geomagnetic latitude (∼40-55°) that the locally uniform geomagnetic field approximation would also have a noticeable effect. The Lines A1-A4 run perpendicular to the Appalachian Mountains in the East Coast of the United States, while Transmission Lines B1-B4 are aligned roughly parallel to the Appalachian Mountains and the coastline. The Transmission Lines A/B 1-4 have lengths of 100, 200, 500, and 1,000 km, respectively, and a resistivity ρ L = 0.01 Ω/km. All lines intersect at (−80°longitude, 36°latitude), indicated by the yellow star in Figure 2. By analogy with Rosenqvist and Hall (2019), we set the earthing resistances at each grounding point to 1 Ω, R 1 ground ¼ R 2 ground ¼ 1. Then, for each transmission line the GIC is computed as the voltage divided by the transmission line system impedance (Equation 16), the latter calculated as Z L = ρ L ×L+2. In a realistic setting, a wide range of grounding resistances may be encountered that will substantially affect 10.1029/2020SW002467 Space Weather line GICs; moreover, GICs in a complex system do not have such a simple relationship to the voltage. Bearing that in mind, we nevertheless find that framing our illustrative discussion in terms of electric currents, obtained by scaling the voltages appropriately by the lengths of the lines, provides for a more insightful analysis compared to analyzing the voltages only.
The EarthScope USArray MT impedance tensors (Schultz et al., 2006) are used to compare the modified GIC computation approach against several different scenarios that are used in the space weather community. A summary of the methods that are discussed and compared in this section is provided in Table 1. The first option employs the full 3-D approach as described in Lucas et al. (2018) and Lucas et al. (2020) to interpolate the magnetic field using the spherical elementary current system (SECS; e.g., Amm & Viljanen, 1999;Rigler et al., 2019) method, compute the electric field at each MT site location, and then integrate the electric field along the transmission line using a Delaunay spatial integration method (e.g., Bonner & Schultz, 2017;Lucas et al., 2018) to compute an induced voltage. The second approach, newly proposed in this paper, is to integrate the MT impedance components along the transmission lines to calculate the power line telluric response functions according to Equations 12 and 13, using Equation 3 to compute the line voltage as a function of the local magnetic field. In this method, a magnetic field time series is chosen for a central location where Transmission Lines A and B intersect (yellow star in Figure 2) and assumed to be uniform over the lengths of the lines. Finally, the third method is a variant of the piecewise 1-D approach similar to, for example, Marti et al. (2014) in Canada and Blake et al. (2016) in Ireland. In this analysis, 1-D impedance tensors are estimated individually at each MT site location, based on the empirical MT impedance tensor (Appendix C) or, in an alternative formulation, on geologically motivated 1-D provinces as determined by Fernberg (2012). Then, these are used in conjunction with the SECS-interpolated geomagnetic fields to calculate the electric fields at each site location, followed by an integration of these electric fields along the transmission lines to compute the voltages. Thus, the piecewise 1-D approach makes no uniform magnetic field assumption but rather follows the same procedure as the 3-D approach, except that local 1-D impedances are computed and used. Integration along power lines was performed by first interpolating the respective values to points Figure 2. Setup of hypothetical stand-alone transmission lines. Eight hypothetical transmission lines have been set up next to the Atlantic coast: Lines A1-A4 perpendicular to the Atlantic coast and the Appalachian Mountains, and Lines B1-B4 that are parallel to the coast and the mountain range, plotted in the Mercator projection. Each line operates independently and is grounded at both ends. All lines are, independently, grounded at the yellow star where they intersect. Black dots indicate the second grounding point for each of the lines, located at 100, 200, 500, and 1,000 km from the yellow star, respectively. Note. This simple setup allows us to compute and compare line GICs as a function of the line voltages. Nominal computational times to compute the electric fields, followed by voltages on eight example transmission lines, are provided for each of the methods, as observed on a single processor of a modern desktop. Note that the methods that involve SECS computation take an additional 0.8652 s; however, it is the voltage computation part that scales with the size of the electric power grid.

10.1029/2020SW002467
Space Weather located at 1-km distances on each line (so that, e.g., the 1,000-km line has 1,000 points used for the power line integration). Computations described in Table 1 are performed in the frequency domain for consistency.
In the time domain formulation, additional efficiencies may be had in the electric field computation as well as the modified approach of Equation 14, but the voltage computation, which involves the power line integration and scales with the size of the power grid, cannot be further optimized.
The input magnetic field for the modified 3-D approach is shown in Figure 3,  Table 2, and the rootmean-square (RMS) difference between the time series is summarized in Table 3. As can be seen from the figures, the modified 3-D approach is in better agreement with the full 3-D approach than all variants of the piecewise 1-D approximation, for all transmission lines. There are several notable features of Tables 2 and 3 that are worth highlighting. In general, estimated GICs for such a stand-alone transmission line (using the full 3-D approach) are higher in the A lines that are perpendicular to the coast and across the Appalachians than in the B lines that are parallel to the coast. Indeed, the electric field in this area is predominantly polarized in the southeast-to-northwest direction (Love et al., 2019). This could be due to geological reasons: the deep crust and upper mantle beneath the B lines, extending along the Piedmont region and into the Crystalline Appalachians to the north, are reasonably resistive, while the Appalachian Valley and Ridge province further inland is a geologically complex fold and thrust belt structures, and the Appalachian Plateau to the east is characterized by upper crustal layers of electrically conducting sediments. It could, additionally, be a manifestation of the coast effect consistent with, for example, Ivannikova et al. (2018). The piecewise 1-D approach, unable to account for the electric field polarization, consistently underestimates the GICs in Lines A1-A4 and overestimates those in Lines B1-B4. For the longest A line, the 1-D methods based on the empirical impedances provide a reasonable approximation, possibly due to the averaging effect suggested by Viljanen et al. (2012). However, they do not provide a reasonable voltage or, equivalently, line GIC estimate for the B4 line of the same length, which lies parallel to the geological strike and thus perpendicular to the electric field polarization. There, the simpler Fernberg (2012) geological provinces provide a better average geoelectric field estimate and thus a closer GIC prediction.
A less obvious observation pertains to the fact that the amplitudes of the GICs do not linearly scale with the lengths of the lines. Indeed, the A-line GICs increase from A1 to A2, both of which overlie almost exclusively resistive rocks, but then decreases for the longer lines (A3 and A4), which cross over to the conductive basin  (Amm & Viljanen, 1999;Rigler et al., 2019) for the central location of our synthetic power grid, indicated by the yellow star in Figure 2. We calculate all quantities over the full time series.
in the west side of the Appalachians. Further, as is expected, the longer the transmission lines, the greater (generally) are the differences between the full 3-D and the modified 3-D approaches; however, even for the relatively strong Halloween storm, the maximum amplitude differences between the two approaches never exceed 50% except for the 1,000-km Line B4 and generally stay within 10%. This illustrates the limits of validity of the modified GIC approximation. While the dynamics of a realistic power grid would be different, we can confidently assert that the approximation is valid to well within 50% accuracy for power lines ≪1,000 km in length that are located at middle to low latitudes. In realistic power grid modeling, additional boundary effects would be observed where existing connected lines (such as those extending into the southern United States, currently uncovered with MT surveys) are not included in the model, increasing the GIC estimation errors around the edges.

Discussion
The physics of 3-D 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 obtained from 1-D modeling is a frequency-dependent complex scalar, the 3-D Earth impedance is a tensor (Equation 7) 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 direction, relative to the 1-D approximation. Because of this, the geoelectric field exhibits significant spatial variability that is not, in general, captured by Equation 5 under the uniform geoelectric field approximation. Moreover, an inherent assumption of 1-D modeling is that the layers extend to infinity, laterally. A collection of local 1-D models violates this assumption and

Space Weather
should itself be modeled in 3-D. To sum up, unless the use of Equation 5 is justified by fitting the parameters to existing GIC data at a particular transformer, the uniform geoelectric field assumption does not suffice for an accurate GIC calculation, and more detailed Earth's structure needs to be taken into account. As is shown here, a reasonably straightforward extension to Equation 5 that allows for frequency dependence of the linear regression coefficients and forgoes the use of the uniform electric fields for the uniform magnetic field time series inputs, as proposed by Weigel and Cilliers (2019), allows to better account for the 3-D Earth structure provided that the 3-D measured or modeled impedances are available for the area in question. In contrast to Equation 5 applied with a parameter fitting approach (e.g., see Pulkkinen et al., 2007;Viljanen et al., 2006, in Finland;Pulkkinen et al., 2010, in Japan;Liu et al., 2009;Zhang et al., 2015, in China;Bernhardi et al., 2008;Matandirotya et al., 2015;Ngwira et al., 2009, in South Africa), our modified method in Equations 14 and 15 can be seamlessly adapted to a changing power grid system configuration by, for example, excluding or including the relevant power line terms from the sum.
The quantities Z xx , Z yx , Z xy , and Z yy needed to fully account for the 3-D Earth in GIC modeling (Equation 15) are as defined in Equation 7. Thus, in the frequency domain, they exactly correspond to the entries of the 3-D MT response tensor (rescaled MT impedance) as available, for example, through the IRIS electromagnetic transfer functions (EMTF) online database (Kelbert et al., 2011). A challenge with setting this up for any specific power grid is the sparseness of these observations, in conjunction with very local near-surface distortion effects that are incorporated in these measurements. One solution to the problem is to use high-resolution gridded impedances such as those that are obtained from a 3-D conductivity model (see, e.g., Kelbert et al., 2019, for details and discussion). An alternative solution is obtaining a more dense MT data coverage in areas that require higher accuracy in GIC modeling than is currently achieved. This latter option has certain Note. Maximum absolute GIC value for the full 3-D calculation, compared to the maximum absolute GICs using the modified 3-D and the piecewise 1-D approaches, computed for the entirety of the Halloween storm time series as shown in Figure 3. "Modified" stands for the method proposed in this paper. "3-D no mag" stands for full 3-D computation with integration of electric fields along the power lines, but using a uniform geomagnetic field as input. "1-D" stands for Z 1 (Equation C1) and "1Dssq" stands for Z ssq (Equation C2). Figure 3 Line Modified Note. "Modified" stands for the method proposed in this paper. "3-D no mag" stands for full 3-D computation with integration of electric fields along the power lines but using a uniform geomagnetic field as input. "1-D" stands for Z 1 (Equation C1) and "1Dssq" stands for Z ssq (Equation C2). Note that as transmission lines get longer, the errors related to the modified GIC approach and the benchmark 3-D method that ignores the variation in the geomagnetic field, generally, increase but stay way within the errors of all variants of the piecewise 1-D approach.

10.1029/2020SW002467
Space Weather limitations associated with cultural noise effects. High-quality MT data can only be obtained in areas sufficiently removed from existing power lines, railways, and other infrastructure that interfere with the natural ionospheric and magnetospheric signals.
The value of accounting for the accurate 3-D Earth conductivity cannot be underestimated. Indeed, as we move away from the simplified formulation of Equation 5 to describe the Earth as a collection of 1-D models, regional 1-D representation of the Earth and even the potentially promising 1-D approximations based on empirical impedances (Appendix C) cannot, in general, recover comparable to the 3-D setting GIC amplitudes over a geologically heterogeneous region of the Earth. As our synthetic example designed in section 4 clearly illustrates, the geoelectric field polarization, which is only captured when a 3-D Earth structure is incorporated into the analysis, may play a significant role in the GIC calculation even in our simple example where line GIC is linearly related to the line voltage. In a realistic power grid, this relationship will be substantially more complex.
This effect is further illustrated by the differences in the 3-D versus 1-D power line telluric mapping functions τ N (ω) and τ E (ω) (Appendix C). The frequency-dependent τ responses determine the power line voltages for a spatially uniform external magnetic field and a given 3-D or piecewise 1-D Earth structure. Figures C1-C3 illustrate the differences in τ responses between the 3-D and various 1-D approaches and show that substantial differences exist particularly at shorter periods that do not, in general, diminish over longer lines as would be expected if the power line integration resulted in averaging over the 3-D Earth structures. However, the averaging effect exists in some areas that are primarily defined by smaller scale spatial structure of Earth's electrical conductivity variations and is more pronounced at longer periods. This may explain why it is typically harder to fit the higher-frequency GIC measurements with the 1-D approach. However, certain areas with significant large-scale geological variability, such as the Appalachians on the east coast of North America, and the midwestern region, exhibit 3-D GICs (evidenced in the non-1-D τ responses) across all the relevant temporal scales. As illustrated with the synthetic example, this effect is especially noticeable for power lines that cross a major geological discontinuity such as Lines A1-A4, or those that align with or are perpendicular to a major geologic strike such as Lines B1-B4.

Conclusions
The relative scarcity of the analysis of laterally heterogeneous geoelectric fields in the academic GIC literature is a result of practical complexity of this approach. Applied to historical geomagnetic storms, it requires knowledge of the spatially and temporally variable ground-level geoelectric fields at the time, integrated over power lines to obtain voltages, as demonstrated by Lucas et al. (2020). Similarly, performed as part of real-time hazard analysis, it requires a similar sequence of operations: a rather complex and computationally intensive real-time geoelectric field computation, followed by a real-time integration over power lines. The latter task should be made simpler, at least in the United States, by the release of the near real-time Geoelectric Field Maps data product by National Oceanic and Atmospheric Administration Space Weather Prediction Center (NOAA SWPC). While this approach is by far the most comprehensive, here we propose an approximation to geoelectric field power line integration that should fill the gap between the linear uniform geoelectric field approach and the full-scale analysis based on realistic geoelectric fields. Specifically, we derive a simple formulation that incorporates the realistic 3-D Earth MT impedances into the power grid system model, resulting in an additional set of static parameters, power line telluric response functions (τ responses), associated with each power line. If the geomagnetic field, rather than the geoelectric field, may be assumed locally uniform in space, the new set of power line parameters may be applied in lieu of real-time power line integration to obtain the line voltages as a function of time. The integration along power lines to obtain τ responses may thus be viewed as a preprocessing step that streamlines the real-time calculations.
Weigel and Cilliers (2019) have explored the validity of the frequency independence assumption of the power grid system coefficients of the linear model of Equation 5 and concluded that allowing for a frequency dependence of these coefficients provides for a better GIC predictive power even locally and that using the geomagnetic field as input in conjunction with allowing for frequency-dependent coefficients is further shown to be significantly better at estimating GIC than models that use the geoelectric field as input.

Space Weather
Indeed, as we see from Equation 15 (see also Appendix B), allowing for the frequency dependence in these coefficients partly compensates for the inadequacy of the uniform geoelectric field assumption; however, the reformatting of this equation in terms of a uniform geomagnetic field over the system grid bypasses this inaccurate approximation altogether, resulting in a preferred formulation that is, as we show here, mathematically rigorous. Weigel and Cilliers (2019) propose that this reformulation may be used when GIC and magnetic field measurements are available; however, they do not provide a path forward when the frequency-dependent coefficients cannot be obtained through a least squares fit of an existing GIC time series measurements. In contrast, this study provides an alternative approach for constructing those coefficients. Indeed, GIC measurements are not required to accurately construct the frequency-dependent power grid system coefficients, as long as the relevant MT impedances are available. This is a dramatically more practical approach since the MT impedances are static in time, in that they do not need to be continuously collected. They are also easy and relatively cheap to obtain.
Moreover, this approach allows a power grid utility to use this approximation everywhere in the power grid, rather than only locally where a GIC measurement was once obtained, and to adjust the coefficients in real time as the power grid configuration changes. We therefore suggest that our proposed formulation is a valuable step forward for operational GIC monitoring.
The validity of the uniform magnetic field assumption may itself be challenged for power grids covering large enough spatial areas. However, over each power line, a milder assumption is in fact employed: The geomagnetic field is averaged along the lengths of the power lines. Computational improvements may also be obtained if the magnetic field is instead approximated by a local value at a power grid substation. As the synthetic analysis in section 4 illustrates, the spatial variations in the magnetic field at relevant spatial scales to the power grid geometry are unlikely to cause significant inaccuracy in GIC computations, compared to that caused when the 3-D Earth is ignored, except at high latitudes and/or for very long power lines. However, this preliminary conclusion should by all means be critically considered and could be relaxed with future improvements to magnetic field monitoring and modeling. The modified algorithm proposed in this work, while computationally convenient, may need to eventually be replaced with the rigorous integration of spatially heterogeneous geoelectric field values, obtained perhaps from the operational Geoelectric Field Maps data product of NOAA SWPC, along each of the power lines. In the meantime, it may provide convenience and substantial improvement over many of the existing practices in GIC estimation in the framework of power grid system modeling.
The expressions in the square brackets are impulse responses. Equation 10 follows.

Appendix B: Reduction to the Classical Frequency-Independent Model
If the Earth's electrical conductivity structure is assumed laterally homogeneous, the space dependence of the telluric mapping functions is gone, and the MT impedance tensor is reduced to the 1-D form, which asserts Z xx = Z yy = 0 and Z xy = −Z yx = Z 1 , so that in either the frequency or time domain, and where x k and x i are the power line endpoints, and L ki is the length of the line, depend merely on the local 1-D impedance, and the network topology, yielding a 1-D approximation to Equations 14 and 15 at transformer j: Figure C1. The differences in the amplitudes of power line telluric mapping functions τ N (ω) and τ E (ω) defined by Equations 12 and 13 (3-D) and Equations B1 and B2 (1-D) on the high-voltage power grid in the United States at several sample periods (frequencies). Z 1 was employed for the 1-D approximation of the Earth impedance.

Space Weather
or, in the frequency domain, we see that B3 may be rewritten as with the fitting parameters a j and b j defined in the same way as in Equation 5. It follows that if the 1-D approximation holds, the parameters may be obtained using the expressions

Appendix C: One-Dimensional Approximations of a 3-D Impedance Tensor
Different versions of creating a 1-D approximation from a fully 3-D impedance tensor are available in the literature. An optimal 1-D approximation would by necessity be rotationally invariant and result in an accurate average 1-D conductivity profile if inverted for the Earth's structure. This matter has been a subject of much discussion, development, and even debate in the scientific literature in the past (e.g., Lilley, 1998;Marti et al., 2009;Szarka & Menvielle, 1997;Vozoff, 1991;Weaver et al., 2000). While there is not a single optimal solution to this dimensionality reduction problem that we are aware of, here we explore two Figure C3. The differences in the amplitudes of power line telluric mapping functions τ N (ω) and τ E (ω) defined by Equations 12 and 13 (3-D) and Equations B1 and B2 (1-D) on the high-voltage power grid in the United States at several sample periods (frequencies). Z ssq was employed for the 1-D approximation of the Earth impedance.
For the high-voltage power lines in the United States, the differences between the 3-D and the 1-D representations of the power line telluric mapping functions τ N (ω) and τ E (ω) are shown in Figure C1. Here, the 3-D Figure C4. The differences in the amplitudes of power line telluric mapping functions τ N (ω) and τ E (ω) defined by Equations B1 and B2 (1-D) on the high-voltage power grid in the United States at several sample periods (frequencies). Here, the difference between the responses obtained from Z 1 and from Z ssq is plotted.
values of the telluric mapping functions were obtained using the empirical MT impedances (Kelbert et al., 2011;Schultz et al., 2006) and the 1-D representation was obtained using a local 1-D approximations to the empirical 3-D impedance values, Z 1 , which, as we see here, is a better local 1-D approximation compared to the use of broad geological regions (Fernberg, 2012; Figure C2).
Very similar results are observed for the Z ssq 1-D approximation ( Figure C3). However, a more careful look suggests ( Figure C4) that the Z 1 approximation to the Earth impedance consistently underestimates the GICs compared to the Z ssq approximation. Based on Figures C1 and C3, we conclude that either of these simplifications may be appropriate in selected geologically homogeneous regions of the United States, and both are entirely inappropriate in other, geologically complex regions, consistent with the (Lucas et al., 2020) analysis.

Peer Review Disclaimer
This draft manuscript is distributed solely for purposes of scientific peer review. Its content is deliberative and predecisional, so it must not be disclosed or released by reviewers. Because the manuscript has not yet been approved for publication by the U.S. Geological Survey (USGS), it does not represent any official USGS finding or policy.