Estimating the Geoelectric Field and Electric Power Transmission Line Voltage During a Geomagnetic Storm in Alberta, Canada Using Measured Magnetotelluric Impedance Data: The Influence of Three‐Dimensional Electrical Structures in the Lithosphere

Estimating the effect of geomagnetic disturbances on power grid infrastructure is an important problem since they can induce damaging currents in electric power transmission lines. In this study, an array of magnetotelluric (MT) impedance measurements in Alberta and southeastern British Columbia are used to estimate the geoelectric field resulting from a magnetic storm on September 8, 2017. The resulting geoelectric field is compared to the geoelectric field modeled using the more common method that uses a piecewise‐continuous 1‐D conductivity model. The 1‐D model assumes horizontal layers, which result in orthogonal induced electric fields while the measured MT impedance data can account for fully 3‐D conductivity structure. The geoelectric field derived from measured MT impedance data is partially polarized in southern Alberta, and the geoelectric field magnitude is largest in northeastern Alberta where the resistive Canadian Shield outcrops. The induced voltage in the Alberta transmission network is estimated to be ∼120 V larger in northeastern Alberta when using the measured MT impedances compared to the piecewise‐continuous 1‐D model. Estimated voltages on transmission lines oriented NW‐SE in southern Alberta are 10%–20% larger when using the MT impedances due to the polarized geoelectric field. As shown with forward modeling tests, the polarization is due to a feature in the lower crust (20–30 km depth) called the Southern Alberta British Columbia conductor that is associated with a Proterozoic tectonic suture zone. This forms an important link between ancient tectonic processes and modern‐day geoelectric hazards that cannot be modeled with a 1‐D analysis.

Previous work in e.g., the United States, Europe, China, Australia and New Zealand has shown that there can be significant discrepancies between the predicted and observed geoelectric fields, power transmission line voltages, and/or GICs (e.g., Bedrosian & Love, 2015;Dimmock et al., 2019;Divett et al., 2020;Gannon et al., 2017;Kazerooni et al., 2017;Liu et al., 2018;Lucas et al., 2018;Marshall et al., 2019;Shetye et al., 2018;Simpson & Bahr, 2020aSun & Balch, 2019;Torta et al., 2017). The discrepancies have been hypothesized to be due to the fact that deep Earth conductivity structure is 2-D or 3-D which violates the widely used 1-D modeling assumptions (e.g., Gannon et al., 2017;Kazerooni et al., 2017;Kelbert & Lucas, 2020;Lucas et al., 2018;Shetye et al., 2018;Simpson & Bahr, 2021) but these discrepancies have not been explicitly linked to specific geological structures. Other cases have noted the effects of deep oceans around islands and near coastlines (e.g., Liu et al., 2018;Marshall et al., 2019) or used coarse-grained, global 3-D conductivity models (e.g., Honkonen et al., 2018). A similar analysis of the robustness of the 1-D approximation has not yet been applied in Canada which is proximal to the auroral zone and thus susceptible to damaging GMDs (Boteler, 2001). It is important for power transmission companies and regulators to know under what circumstances the 1-D assumption is valid and, if not, what mitigation may be required to avoid unexpected GICs. This paper focuses on the estimation of the geoelectric field in Alberta and southeastern British Columbia and compares two different approaches: (a) "The 1-D Model-Space Method" and (b) "The Data-Space Method". The 1-D Model-Space Method uses zone-based, piecewise continuous, 1-D models of the subsurface conductivity from Trichtchenko et al. (2019) to compute synthetic surface impedance data. These models are a priori 1-D and do not take into account lateral variations or the effects of 3-D electromagnetic induction. In contrast, the Data-Space Method uses measured magnetotelluric (MT) impedance tensors to provide information about the subsurface conductivity without any a priori assumptions (e.g., Lucas et al., 2018;Simpson & Bahr, 2021). The data are determined by geological structures and thus incorporate the complex 3-D electromagnetic effects of irregular geological boundaries, dipping structures, topography, coastlines, etc. which cannot be incorporated into an a priori 1-D model. The comparison between the methods is made by first using a frequency-domain synthetic example and then using a real, time-domain GMD from a large geomagnetic storm on September 8, 2017 (see e.g., Wu et al., 2019). It is worth noting that there are various other methods for estimating the geoelectric field using e.g., uniform conductivity models, singe n-layered conductivity models, thin sheet conductance models, and 3-D conductivity models (Marshall et al., 2019). All these methods can be classified as "Model-Space Methods" of varying degree in that they all use some type of a priori conductivity model to compute synthetic surface impedance data (Kelbert, 2020). The "Data-Space Method" is unique in that it uses measured MT impedance tensor data to calculate geoelectric fields directly. Model-Space methods generally suffer from having to make some simplifying assumptions about the unknown conductivity of the Earth while Data-Space methods suffer from the effects of noisy data and sparse data coverage (Kelbert, 2020). The majority of Model-Space methods discussed in the literature do not incorporate the effects of 3-D electromagnetic induction and instead rely on the 1-D assumption.
Alberta and southeastern British Columbia present an interesting case study to test the 1-D approximation and examine the influence of 3-D Earth structure and lateral variations on geoelectric hazards due to the complex subsurface structure that has been mapped by extensive geophysical and geological exploration. The region is characterized by the geological boundary between the Canadian Cordillera and the Western Canada Sedimentary Basin (WCSB) (Figure 1a). The Cordillera consists of several different NW-SE trending mountain belts of Phanerozoic age (e.g., <0.5 Ga), while the WCSB is a thick (0-5 km) basin of sedimentary rocks of Phanerozoic age (Dickinson, 2004;Porter et al., 1982). The WCSB is underlain by 0.5-2.5 Ga Paleoproterozoic rocks to the north and 2.5 to 4 Ga Archean rocks in the south (Figure 1a). The Archean rocks to the south are known as the Hearne craton which is one the oldest parts of North America while the Paleoproterozoic rocks to the north represent terranes which were joined to the Hearne craton in a series Figure 1. Maps of geoelectric field vectors calculated using a synthetic north-polarized sinusoidal magnetic field with a frequency of 0.01 Hz and magnitude of 1,000 nT. (a) The geoelectric field calculated in each regional zone using the 1-D Model Space Method. Blue squares are population centers mentioned in the text. Green shades approximately correspond to Archean (2.5-4 Ga) and Proterozoic (0.5-2.5 Ga) basement units overlain by the Phanerozoic (<0.5 Ga) Western Canada Sedimentary Basin (WCSB). Blue shades approximately correspond to Cordillera units, and yellow shades are where Canadian Shield outcrops. (b) The geoelectric field calculated using the Data Space Method with measured magnetotelluric impedance data. Yellow squares show the representative site locations of ABT175, SAB060, and ABT272. of collisions during the formation of North America (Whitmeyer & Karlstrom, 2007). The primary suture which links the Archean craton to the Paleoproterozoic terranes is known as the Snowbird Tectonic Zone (STZ). The STZ also forms an important suture in joining the Archean Hearne and Rae cratons further to the northeast (Whitmeyer & Karlstrom, 2007). The boundary between the WCSB and the Cordillera approximately coincides with the geographical division between the provinces of Alberta and British Columbia. This geological boundary is evidenced by clear changes in surficial geology and topography between the relatively flat sedimentary units of the WCSB and the highly deformed, glacially eroded units of the Canadian Cordillera. A simple hypothesis is that the relatively uniform stratigraphy of the WCSB and flat topography of Alberta would be amenable to a 1-D approximation whereas the complex geology of the Cordillera and variable topography of British Columbia would violate the 1-D assumption and have significant 3-D induction effects.
Previous GIC studies on portions of the Alberta Interconnected Electric System (AIES) using a 1-D Model-Space approach found that 7 transformers in Alberta were at high risk and 60 transformers were at medium risk to damage by GICs (Haque et al., 2017). An outstanding question is whether the effects of 3-D electromagnetic induction due to lateral variations in the Earth's conductivity would alter the risk assessment. This paper focuses primarily on the first step in the modeling workflow: the comparison between the geoelectric field and transmission line voltages computed using the 1-D Model-Space Method with those computed from the Data-Space Method. An important question is how much the direction of the induced geoelectric field depends on the direction of the inducing magnetic field versus how much the direction is influenced by Earth's conductivity structure. This is important because if the geoelectric field direction differs over a large regional area when comparing the two methods, it will impact the line integral calculation along transmission lines in the time-domain (Kelbert, 2020). In particular, if the geoelectric field is partially polarized (i.e., its direction is relatively constant and independent of the GMD direction in the time domain) due to interactions with the conductive subsurface, then transmission lines parallel to the polarization direction will have larger voltages than lines perpendicular to the polarization direction. This effect would not be captured by a 1-D analysis and requires methods such as the Data-Space Method which can take into account fully 3-D electromagnetic induction.

Methodology
Geoelectric fields are most often calculated using the 1-D Model-Space Method with an assumption of layered conductivity over a large region, or several 1-D models combined into a single piecewise-continuous model (Kelbert, 2020 and references therein). For any of these 1-D Earth models, the horizontal geoelectric field components can be calculated by multiplying the horizontal magnetic field components (e.g., the GMD) by the impedance tensor in the frequency domain (Bedrosian & Love, 2015):  unless otherwise stated. Note that in the 1-D case, the impedance tensor is anti-diagonal and anti-symmetric and thus the geoelectric field will always be orthogonal to the GMD in the frequency domain.
Power transmission companies and regulators generally calculate the geoelectric field and transmission-line voltages using the 1-D Model-Space method because detailed knowledge about the 3-D subsurface conductivity is often not available (Canadian 1-D examples : Boteler, 1998: Boteler, , 2001Haque et al., 2017;Jacobson et al., 2014;Nikitina & Trichtchenko, 2015;Trichtchenko et al., 2019). In order to perform an analysis which takes into account fully 3-D electromagnetic induction, either an a priori 3-D model must be constructed to produce synthetic impedance data, or the Earth's frequency-dependent impedance must be measured at many locations over a large area. The a priori 3-D model approach is limited in that it requires some simplifying assumptions to be made about the 3-D conductivity structure but it allows for noise-free impedance data to be computed anywhere on the surface. In contrast, measured impedance data make no assumptions about the subsurface and are instead determined by the structure. However, measured data are noisy and sparsely distributed which requires interpolation in both frequency and spatial domains (Campanyà et al., 2019). There are various methods to deal with uncertainties in measured impedance data and interpolation errors including the use of interstation transfer functions but these are beyond the scope of the current study (see Campanyà et al., 2019 for a detailed discussion).
In Alberta and southeastern British Columbia, a large regional array of long-period MT data has been collected at 530 locations using 10-20 km average site spacing south of 55°N and sparser coverage to the north (Figure 1b; Hanneson et al., 2021;Wang, 2019). The MT method is a well-established geophysical exploration method and parts of this array have previously been used to study the geology of Western Canada in the search for hydrocarbon, mineral and geothermal resources (Hanneson et al., 2021;Liddell et al., 2016;Nieuwenhuis et al., 2014;Wang et al., 2018;Xiao & Unsworth, 2006). MT data were collected using several types of instruments. The first data were collected using the Long-period Intelligent Magnetotelluric Systems (LIMS) during the Lithoprobe project with a 4 s sample rate (Boerner et al., 1999). More recent data were collected using the Narod Intelligent Magnetotelluric System (NIMS) data logger with an 8 Hz sample rate (Wang, 2019). Both instruments used a fluxgate magnetometer and non-polarizing electrodes. Several sites in the array are located in the United States and were collected as part of the USArray (Bedrosian & Feucht, 2014). The MT data set was edited to remove outlier sites that had much larger or smaller impedance values than surrounding sites. These outliers may have been influenced by galvanic distortion caused by local heterogeneities which are much smaller than the station spacing (Bonner & Schultz, 2017). For the remaining stations, it is unlikely that galvanic distortion on the length-scale of interest would be present given that the impedance varies smoothly from station-to-station. In total, 513 sites were used in this study, and 17 were removed. This unique regional data set in Western Canada is denser than other large MT arrays such as USArray (70 km site spacing; Lucas et al., 2018) or AusLamp arrays (55 km site spacing; Marshall et al., 2019). As a result, there is higher spatial resolution and fewer interpolation uncertainties.
The MT time series were converted to the frequency-domain to calculate the frequency-dependent wave impedance tensor at measurement points at the surface of the Earth. The wave impedance captures the full 3-D electromagnetic induction processes related to subsurface conductivity structure. For a 3-D conductivity distribution, the impedance in Equation 2 becomes a full non-symmetric matrix that depends on position: Here, , denotes that all variables are dependent on horizontal position on the surface of the Earth. In this case, the frequency-domain geoelectric field computed using Equation 4 would not necessarily be orthogonal to the GMD.
In this study of Alberta and British Columbia, the 1-D Model-Space Method and the Data-Space Method for computing the geoelectric field are compared. For the 1-D approach, we follow the method of Trichtchenko et al. (2019) which divided Alberta into 10 zones and British Columbia into 13 zones based on known geological domains (Figure 1a). The zone numbers for Alberta match the numbering of Trichtchenko et al. (2019). The zone numbers for British Columbia (Zones 11, 12, 13, and 14 in Figure 1a) correspond to British Columbia Zones 3S, 4S, 5S, and 6 in Trichtchenko et al. (2019). In this study, 9 of the zones in British Columbia were outside the region of interest and were not included. In total, 14 zones were used and each was given a different 1-D conductivity model based on previously published models (Trichtchenko et al., 2019) to create a piecewise-continuous 1-D conductivity model over the regional area. The synthetic, noise-free frequency-dependent impedance, E Z , was calculated for each zone using an analytical recursion relation (Trichtchenko et al., 2019;Wait, 1954). The 1-D models and computed data are shown in Figures S1 and S2 in Supporting Information S1. For a given GMD, the synthetic 1-D impedance data can then be used to calculate the geoelectric field using Equation 2. For the Data-Space Method, the measured MT impedance tensor data at each site were used directly to compute the geoelectric field using Equation 4 without any reference to any conductivity model. Thus, the Data-Space Method makes no a priori assumptions about the Earth's structure.
Once the geoelectric field is computed at a particular location (using either the impedance derived from the 1-D piecewise continuous model or measured MT impedance tensor data), the field can be converted to the time-domain using the inverse Fourier transform (Bedrosian & Love, 2015;Kelbert et al., 2017;Lucas et al., 2018). The time-domain geoelectric field,   , E t e r , can then be used to calculate the line integral along transmission lines to estimate the voltage in the power transmission line at a given time, .
In order to examine whether there are significant differences between the 1-D Model-Space Method and the Data-Space Method when modeling the geoelectric field induced during a GMD, both methods are used to compute the geoelectric field for two scenarios. The first is a synthetic, idealized GMD with a frequency of 0.01 Hz and the second uses a real GMD from September 8, 2017. For the real GMD scenario, the computed frequency-domain geoelectric field is converted back to the time domain to examine how the two modeling methodologies compare in the time-domain. Finally, geoelectric field line integrals along >240 kV power transmission lines in Alberta are computed to assess the line voltages.

Synthetic Geomagnetic Disturbance
Under the 1-D assumption, the induced geoelectric field will always be orthogonal to the inducing geomagnetic field at a particular frequency due to the symmetry of the impedance tensor as per Equation 2. A synthetic example is shown in Figure 1 using an idealized sinusoidal GMD that is magnetically polarized in the north-south direction with a frequency of 0.01 Hz and a magnitude of 1,000 nT. The frequency is chosen because spectral power of real GMDs is often concentrated around 0.01 Hz (Bedrosian & Love, 2015). The magnitude is chosen as representative of a large event with a time derivative ( dB dt fB /  2 0  ) of around 3,600 nT/min. Thomson et al. (2011) estimated that a 1-in 100 year GMD event could have a magnitude as high as 5,000 nT with dB dt / between 1,000-4,000 nT/min. Of course, a real GMD would not have a simple sinusoidal signal but this simplistic synthetic signal is primarily used for ease of illustration.
As expected for a north-south polarized GMD, the calculated geoelectric field using the 1-D method is oriented in the east-west direction ( Figure 1a). A single geoelectric field vector is shown for each regional zone. The strength of the geoelectric field is scaled by the magnitude of the impedance at 0.01 Hz as per Equation 2 which, for a 1-D regional assumption, is independent of the location within a zone (Trichtchenko et al., 2019). Note that the orthogonality is maintained even when crossing zone boundaries for 1-D models; it is only the magnitude of the calculated geoelectric field that changes when crossing zone boundaries.
Using the same idealized 1,000 nT, 0.01 Hz GMD as an input, the geoelectric field can be calculated using the measured impedances at 0.01 Hz for each MT site location. This calculation is shown in Figure 1b where, instead of only one geoelectric field vector for each regional zone in the piecewise continuous model, there are now 513 geoelectric field vectors, one at each MT site. It can be seen that there are deviations from the 1-D assumption and this results in a much more complex geoelectric field that is not always orthogonal to the input geomagnetic field. What is most notable is that the geoelectric field vectors are uniformly deflected to the northwest over much of the survey area and this is most extreme in southern Alberta near Calgary. The magnitude of the geoelectric field also varies across the province and does not align cleanly with zone boundaries. For example, the eastern areas of Zones 6 and 7 near Fort McMurray have a larger magnitude geoelectric field than the western areas of the same zones. The 1-D Model-Space Method cannot account for this intra-zone variability. In this simple example, the GMD is spatially homogeneous so the heterogeneity in the geoelectric field is most likely due to Earth structure. This analysis is also completed for a 0.1 and 0.001 Hz GMD and the results are shown in Figures S3 and S4 in Supporting Information S1.
The above analysis was done in the frequency domain using a single frequency (0.01 Hz) and an idealized sinusoidal polarized GMD source. This was done to simplify the problem to explicitly examine the differences between the 1-D Model-Space Method and Data-Space Method. However, the spatial-and temporal-dependence of a real GMD means that there will be broadband frequency content in the signal, and the impedance data themselves are also frequency dependent. The next section explores the additional complexity of a real GMD.

Real Geomagnetic Disturbance
To examine the effects of using the 1-D Model-Space Method versus the Data-Space Method during a real GMD, the magnetic storm that occurred on September 8, 2017 was used as a case study because it is one of the largest geomagnetic storms of the last solar cycle for which there is abundant 1 Hz magnetic data available (Blagoveshchensky & Sergeeva, 2019; Clilverd et al., 2021;Wu et al., 2019). The magnetic storm had two distinct peaks in intensity around 02:00 UT and 14:00 UT and here we analyze the second peak (from 06:00 UT to 23:59 UT) which had a larger amplitude over the area of study. This storm was associated with measurable and unexpectedly large GICs in transmission networks in several mid-latitude and high-latitude countries (e.g., Clilverd et al., 2021;Dimmock et al., 2019) but no data has been publicly published for Alberta. Horizontal magnetic field data with a 1 Hz sample rate were downloaded from the NRCAN, USGS, and CARISMA (Mann et al., 2008) magnetometer networks for 22 magnetic field observatories over the region of study (see Figure S5 in Supporting Information S1). An example time series and spectrum from the Meanook observatory in Alberta is shown in Figure S6 in Supporting Information S1. The magnetic field was interpolated onto a 0.5° by 0.5° grid at each time step and then converted to the frequency domain with a 10,000-sample padding following Kelbert et al. (2017). The magnetic field observatories are separated by much greater distances than the grid size and thus small-scale fluctuations in the magnetic field are not captured in this analysis.

Frequency-Domain Analysis
For the 1-D Model-Space Method, the geoelectric field was computed in the frequency domain at interpolated grid points using the 1-D synthetic impedance for the zone in which the grid point was located. The resulting geoelectric field is shown in Figure 2a at 0.01 Hz (the inducing magnetic field at the same frequency is shown in Figure S7 in Supporting Information S1). For the Data-Space Method shown in Figure 2b, the measured MT impedance at each MT site location along with the nearest interpolated magnetic field grid point were used to compute the geoelectric field at each MT site location using Equation 4. Maps of the geoelectric field at 0.1 and 0.001 Hz are shown in Figures S8 and S9 in Supporting Information S1 for the 1-D and 3-D methods, respectively.
As shown in Figure 2, when using a real GMD, the input geomagnetic field varies spatially (see Figure S7 in Supporting Information S1) which complicates the overall patterns for both the 1-D synthetic impedance and the measured impedance calculations and may obscure the influence of regional geology. However, even given the spatial heterogeneity of the source, it can be seen most clearly in southern Alberta that the geoelectric field is deflected to the NW-SE direction when using the measured MT impedances while the 1-D Model-Space Method results in a west-oriented field. In northern Alberta there is a broader agreement in the direction of the geoelectric field whether using the 1-D synthetic or the measured impedances.
There are also some notable discrepancies in the magnitude between the two methods at 0.01 Hz. In northeastern Alberta near Fort McMurray, the estimated geoelectric field amplitude using the measured impedances is larger than when using the 1-D synthetic impedances while in southeastern Alberta, the estimated magnitude when using the measured impedances is nearly half the estimate when using 1-D synthetic impedances. However, at lower frequencies (e.g., 0.001 Hz; Figure S9 in Supporting Information S1), the magnitude of the geoelectric field in northeastern Alberta is smaller when using the Data-Space Method which may indicate more conductive rocks at greater depths in the measured MT data compared with the 1-D zone-based models.  Figure S10 in Supporting Information S1 for measured MT data at each site). For these sites, the horizontal electric field vector direction is computed at all frequencies (from 0.5 Hz to 10 −5 Hz) using both the 1-D Model-Space Method and Data Space Method. This is shown as histograms in Figure 3. At Site ABT175, the difference in the histograms when using the 1-D synthetic or the measured impedances is relatively small and the histogram is relatively flat indicating that both the source GMD and the resulting geoelectric field are randomly oriented in the frequency domain. In contrast, there is a significant difference in the methods at site SAB060 in southern Alberta. Here, the estimate using the 1-D Model-Space Method results in a relatively flat histogram indicating a randomly oriented GMD source field, but the Data-Space Method using measured impedances shows that at nearly all frequencies, the geoelectric field points NW-SE. The effect at ABT272 in northeastern Alberta (where the magnitude is larger) lies between the extremes of ABT175 and SAB060. At ABT272, there is some evidence that the geoelectric field is preferentially oriented in the northeast-southwest direction at some frequencies, although it is not as extreme as in southern Alberta at SAB060. This suggests that significant 2-D or 3-D subsurface structures or lateral variations are influencing the measured MT data, especially at SAB060, which is not captured by the 1-D conductivity models.

Time-Domain Analysis
To examine the time-domain geoelectric field behavior, the inverse Fourier transform is applied to the frequency-domain geoelectric field. A map of the geoelectric field calculated using the 1-D Model-Space Method and the Data-Space Method is shown in Figure 4 at 14:02:22 UT. This time was chosen because it is approximately halfway between the average peak geoelectric field calculated with the 1-D Model-Space Method (14:02:29) and the average peak geoelectric field calculated with the Data-Space Method (14:02:11). As in the frequency domain, it can be seen that in the time domain, the geoelectric field is relatively spatially smooth when using the 1-D method reflecting the spatial smoothness of the interpolated GMD. At this particular time step, the geoelectric field calculated using the 1-D method points roughly to the north or northeast over much of Alberta. The magnitude is largest in north-central Alberta and relatively small in northwestern and southern Alberta. Note that a different GMD (or a different time step) could have a larger magnitude in some other area of the province as the magnitude varies both spatially and temporally. The geoelectric field calculated using the Data-Space Method has a similar pattern of magnitude but the direction of the geoelectric field in southern Alberta is deflected to the northwest.
As in the frequency-domain maps in Figure 2 which only show a single frequency, the maps shown in Figure 4 only show a single time step and the NW-SE deflection of the geoelectric field may not be significant at all time steps. Histograms of the direction of the geoelectric field at all times (06:00 UT to 23:59 UT) for representative sites ABT175, SAB060, and ABT272 are shown in Figure 5. The estimate of the direction of the geoelectric field at ABT175 in northern Alberta when using the 1-D Model-Space Method is effectively random with a relatively flat histogram showing no preferential direction. When using the Data-Space Method, there is some directionality to the geoelectric field in a NW-SE oriented direction perhaps related to low-frequencies which have NW-SE oriented electric fields across much of the province (see Figure S4 in Supporting Information S1). In southern Alberta at SAB060, the direction of the geoelectric field when using the 1-D Model-Space Method has some subtle preference in the east-west direction perhaps related to the source GMD preferentially being oriented north-south. However, when using the measured MT impedances in the Data-Space Method at SAB060 the geoelectric field points approximately NW-SE at nearly all times. This indicates that the geoelectric field is partially polarized by 3-D geological structure as measured by the measured MT impedances (i.e., the geoelectric field nearly exclusively oscillates in the NW-SE direction through time). Site ABT272 in northeastern Alberta lies between these two extremes. When using the Data-Space Method at Site ABT272, the geoelectric field points in the northeast-southwest direction at most times.
The geoelectric field time series for all three sites is shown in Figure 6 from 12:00 UT to 15:00 UT. It is difficult to visualize the direction of the geoelectric field by looking at the time series, but the correlation between the 1-D Model Space Method and the Data-Space Method is greatest for ABT175 (e.g., the site which best agrees with the 1-D approximation) and least correlated for SAB060 (which shows the greatest polarization in the geoelectric field derived from measured MT impedances and thus the greatest violation of the 1-D assumption). In fact, the y-component (i.e., east-west) of SAB060 is partially anti-correlated with a negative correlation coefficient. The correlation for ABT272 is quite well-correlated relative to SAB060. Also note that at site ABT272, the magnitude of the geoelectric field is much greater when using the Data-Space Method than when using the 1-D Model-Space Method whereas at the other sites, the magnitude is comparable (ABT175) or much smaller (SAB060) when comparing the two methods. The magnitude of the fields is also comparable to geoelectric fields estimated in other studies with similar storms (e.g., Love et al., 2018;Lucas et al., 2018;Simpson & Bahr, 2020a)

Calculating Line Voltages on Transmission Lines
The fact that the geoelectric field estimated with the measured MT impedances can be partially polarized has implications for voltage calculations along transmission lines since the magnitude of the line voltage will vary depending on whether the transmission line is parallel or perpendicular to the geoelectric field. In the case of southern Alberta, transmission lines oriented NW-SE may be more susceptible to GICs than NE-SW oriented transmission lines. This additional directional information would not be captured by an analysis which assumes 1-D conductivity models and 1-D electromagnetic induction. It is also important to note that, when comparing methods, one method does not always result in larger magnitudes than another method. In southern Alberta, the 1-D Model-Space Method results in larger magnitudes than the Data-Space method whereas in northeastern Alberta, the opposite is true.
To investigate the effects of the variations in magnitude and direction of the geoelectric field between methods, line voltages along >240 kV transmission lines in a simplified version of the AIES were calculated at times between 13:45 UT and 14:15 UT. The line voltage was calculated along long transmission lines between substations but not every substation or line in the network was included for simplification purposes. Only line voltages are calculated since currents (in amperes) require additional information about  transformer groundings and line resistances. The line voltage was calculated at each time step using a nearest neighbor interpolation of the geoelectric fields and numerical line integration . The difference in the line voltage calculated when using the 1-D Model-Space Method and the Data-Space Method is shown at 14:02:22 UT in map view (Figure 7). Maps of the absolute voltages for each method at that time step are shown in Figure S11 in Supporting Information S1 and data for this time step are tabulated in Table S14 in Supporting Information S1. As can be seen in Figure 7, the difference in the computed line voltage is greatest in northeastern Alberta. Lines near Fort McMurray (Lines 6, 7, 8) had maximum voltages of 139 V, 180 V, and 223 V, respectively, when using the 1-D Model-Space method. When using the Data-Space method with measured MT impedances, the voltages were estimated to be 207 V, 301 V, and 311 V, respectively. This is an increase of 40%-70% compared to the voltages anticipated using a 1-D analysis (see Table S14 in Supporting Information S1). These voltages are similar to estimates made in similar studies . These three lines near Fort McMurray not only have a larger line voltage at the particular time step shown in Figure 7, but also have a larger line voltage at nearly all time steps evaluated during the peak of the geomagnetic storm (83%, 92%, and 97% of the time steps, respectively; Table S14 in Supporting Information S1).
Over much of the AIES, the line voltage calculated using the Data-Space Method is less than the 1-D Model-Space Method because the measured impedances (and resulting geoelectric fields) have smaller magnitudes than the 1-D synthetic impedance (see Figure 4 and Figures 6a-6f for representative sites in southern and northwestern Alberta). However, the partial polarization of the geoelectric field calculated using measured MT impedances appears to play a role because transmission lines which have segments oriented approximately NW-SE in southern Alberta (e.g., Lines 18, 27, and 33) all have 10%-20% larger line voltages when using the Data-Space Method compared to the 1-D Model-Space Method, whereas all other transmission lines outside of northeastern Alberta have smaller line voltages when using the Data-Space Method. Lines 18 and 33 also have larger line voltages when using the Data-Space Method more often (61% and 55% of the time, respectively) compared with the average for other lines (39%), although this does not hold for Line 27 (34%) (see Table S14 in Supporting Information S1). The impact of the NW-SE polarized geoelectric field is also shown using a hypothetical line oriented NW-SE from Calgary to Nordegg, which crosses the area of greatest polarization (Line 36). Such a transmission line would have the fourth largest line voltage (after Lines 6, 7 and 8) when using the Data-Space method at 14:02:22 UT, with a voltage over 50% greater than the 1-D Model-Space Method (31 vs. 20 V). This hypothetical line also has a larger line voltage using the Data-Space Method over 95% of the time reflecting the fact that the Data-Space geoelectric field aligns with the transmission line orientation at nearly all times. These discrepancies in Lines 18, 27, 33 and 36 are despite the fact that the absolute magnitude of the geoelectric field is often larger when using the 1-D Model-Space Method (e.g., SAB060 in Figure 6). Thus, the above analysis suggests that regional partial polarization of the geoelectric field may impact the final assessment of line voltages and associated GICs on transmission lines and assessing the magnitude only is not sufficient.

Discussion
The 1-D assumption of layered conductivity without lateral variations would result in geoelectric fields oriented perpendicular to the driving GMD in the frequency domain. When using the measured MT impedances, the departure from the 1-D assumption is especially pronounced near Calgary (e.g., Zones 2 and 3) whereas areas further from the Cordillera (e.g., Zone 7 and eastern parts of Zone 6) tend to have a geoelectric field which more closely resembles the geoelectric field computed using the 1-D synthetic impedances. An important result is that the partial polarization is widespread over much of southern Alberta when using the measured MT impedances which influences the line integrals when computing the voltage along transmission lines. The Calgary-Edmonton corridor in southern Alberta is the most densely populated area of the province and has more infrastructure than elsewhere in Alberta including several 500 kV transmission lines. This is also the region that shows MT data with the greatest departure from the 1-D assumption making it all the more important to ensure the best modeling methods are used given the geology. The particular GMD which was investigated in this paper happens to not have such a large peak magnitude in southern Alberta, but a larger GMD could impact southern Alberta in the future. Discrepancies between the 1-D Model-Space Method and the Data-Space Method related to geoelectric field polarization would only be magnified during a larger GMD. It is also worth noting that we do not consider pipeline infrastructure in this study, although such infrastructure is also susceptible to GMDs (e.g., Trichtchenko & Boteler, 2002). There are several major gas pipelines currently undergoing expansion that are oriented NW-SE in the Calgary-Nordegg corridor where the geoelectric field is most strongly polarized.

The Cause of Geoelectric Field Polarization in Southern Alberta
The geoelectric field polarization in southern Alberta estimated using the measured MT impedances is likely due to the unique subsurface geological structure and tectonic history of this region. A number of inver- sions of the MT data have been undertaken to produce 3-D resistivity models of Alberta, including a recent 3-D inversion of the entire data set by Wang (2019) and this model is shown in Figure S12 in Supporting Information S1. Seismic, magnetic and gravity studies have also provided information about lithospheric structure (Chen et al., 2019;Ross, 2002).
Together these studies show that the Western Canada Sedimentary Basin (WCSB) is thickest in western and southwestern Alberta with thicknesses between 3-4 km near Nordegg, and that the WCSB pinches out to the northeast near Fort McMurray. As a result, there is an abrupt change from the thick, relatively conductive sediments of the WCSB to the relatively resistive deformed units of the Canadian Cordillera. The varying WCSB thickness, and the conductivity contrast between WCSB and Canadian Cordillera are relatively shallow but may influence geoelectric fields in southern Alberta. In addition, the WCSB is underlain by the thick (∼250 km), resistive (>1,000 Ωm) lithosphere of the Precambrian North American craton, while the Canadian Cordillera is underlain by an anomalously thin (∼150 km) lithosphere and a shallow, conductive (<60 Ωm) asthenosphere as evidenced in previous geophysical studies (Bao et al., 2014;Chen et al., 2019;Gu et al., 2011;Rippe et al., 2013). Lithospheric thickness has been previously shown to affect induced geoelectric fields (Simpson & Bahr, 2020b). Finally, the resistive basement of the North American craton beneath the WCSB is punctuated by several lower crustal conductors (<10 Ωm), most prominently the Southern Alberta British Columbia (SABC) conductor (also called the Loverna Conductor by some authors) (Gough et al., 1982;Nieuwenhuis et al., 2014;Wang, 2019). The SABC is associated with the Paleoproterozoic Snowbird Tectonic Zone, a major crustal suture where a Proterozoic terrane collided with the Archean Hearne craton (Figure 1a; Ross, 2002). The SABC conductor trends northeast-southwest in the lower crust (20-30 km depth) and may be due to enrichment and/or graphite associated with ancient Proterozoic subduction events during the assembly of the North American craton (Gough et al., 1982;Nieuwenhuis et al., 2014;Ross & Eaton, 2002;Wang, 2019). It is worth noting that the 1-D piecewise-continuous model of Trichtchenko et al. (2019) was developed using knowledge from previous geophysical studies. Zones 1 through 10 and Zone 14 include a conductive upper layer corresponding to the WCSB in Alberta and northeastern British Columbia, while Zones 11 through 13 include resistive upper layers corresponding to the resistive Cordillera (see Figure S1 in Supporting Information S1). Zone 5 approximately corresponds to the location of the SABC and the 1-D model for Zone 5 has a more conductive middle and lower crust compared to the surrounding zones. However, the impedance for each zone is calculated independently, and does not incorporate the 3-D electromagnetic inductive effects of surrounding zones.
To systematically test whether the WCSB thickness, the change in the lithospheric thickness, and/or the lower crustal conductors are responsible for the distribution of geoelectric fields estimated from the measured MT impedances, a "3-D Model-Space Method" can be used whereby 3-D conductivity models are constructed and synthetic MT impedances are calculated from the 3-D model. This process of generating synthetic MT impedance data is similar to the 1-D Model-Space Method, except it now uses finite difference numerical methods to compute electric and magnetic fields which take into account the effects of fully 3-D electromagnetic induction. This process is far more computationally intensive than the 1-D method. The impedance data from the piecewise-continuous 1-D model can be computed analytically on the order of milliseconds using a single processor, while the numerical calculation of impedance data using a 3-D model is orders of magnitude slower and takes ∼30 min using tens of parallel processors. Two separate simulations were run using two different three-dimensional models of conductivity over the region of study to test which features are responsible for the patterns of geoelectric fields observed in the measured MT data (see Figure S13 in Supporting Information S1 and Figure 8). The first model included a 50 Ωm layer to represent the WCSB at the surface which thickens to the southwest, and a 60 Ωm layer to represent the asthenosphere which is at 100 km depth beneath British Columbia and 250 km depth beneath the North American craton based on Bao et al. (2014), Rippe et al. (2013), and Wang (2019) (see Figure S13 in Supporting Information S1). The lithosphere was modeled at 1,000 Ωm (Trichtchenko et al., 2019). The second model, shown in Figure 8, was identical to the first, but also included a northeast-southwest trending 1 Ωm conductor to represent the SABC in south-central Alberta between 20 and 30 km depth based on Gough et al. (1982), Nieuwenhuis et al. (2014) and Wang (2019) (Figure 8). Note, that in both models, ocean bathymetry was included in padding cells with ocean model cells set to 0.3 Ωm. Synthetic MT impedance data were computed for each conductivity model using the ModEM finite difference modeling algorithm (Kelbert et al., 2014). The synthetic MT impedances were then used to estimate the geoelectric field given a synthetic north-polarized geomagnetic source field (see Section 3.1). The first model showed no deflection or polarization in the geoelectric field in Alberta indicating that neither (a) the lithospheric thickness nor (b) the sedimentary basin thickness was the cause (see Figure S13 in Supporting Information S1). The geoelectric fields computed from this model resembled the 1-D Model-Space Method shown in Figure 1a. However, the second model clearly showed a strong deflection of the geoelectric field in the NW-SE direction across southern Alberta (Figure 8d) confirming that the lower crustal SABC conductor is the likely cause of the geoelectric field polarization. The resulting geoelectric field resembles the field calculated using the measured MT data shown in Figure 1b. Note that in Figure S4 in Supporting Information S1, the geoelectric field over much of northern Alberta is also oriented NW-SE at low-frequencies (0.001 Hz) for a north-polarized magnetic field. Neither the SABC nor the change in lithospheric thickness can explain this widespread deflection in northern Alberta and perhaps some other much deeper structure is the cause. Speculatively, it may be related to the Snowbird Tectonic Zone or Birch Hills conductor at depths of 100-200 km (Wang, 2019). Regardless, these low frequencies do not seem to have as much impact on the time-domain geoelectric field or line voltage estimates.
The direction of the geoelectric field is primarily oriented NW-SE because of the complicated 3-D electromagnetic induction associated with NE-SW trending SABC. Figure 8 shows the effect for a north-polarized magnetic source field whereas the geoelectric strike of the SABC trends NE-SW. If the inducing magnetic field is oriented NW-SE (i.e., perpendicular to the geoelectric strike of the SABC), then the transverse electric mode associated with 2-D induction is approximated resulting in near-zero geoelectric fields inside the conductor and stronger geoelectric fields outside the conductor (oriented NE-SW). If the inducing magnetic field is oriented NE-SW (i.e., parallel to the geoelectric strike of the SABC), then the transverse magnetic (TM) mode associated with 2-D induction is approximated resulting in NW-SE oriented geoelectric fields which are relatively unaffected by the SABC. As a result, regardless of the incoming magnetic field direction, the geoelectric field is either near-zero, or is oriented NW-SE. Thus, regardless of the incoming GMD, transmission lines in southern Alberta oriented NW-SE will generally see much higher quasi-DC voltages than transmission lines in other orientations. This 3-D modeling exercise confirms that the SABC is the primary cause of the geoelectric field polarization in southern Alberta and suggests that lower crustal conductors related to ancient tectonics can play a significant role in geoelectric field modeling even in relatively undeformed sedimentary basins and cratonic geologic settings.

Cause of Stronger Geoelectric Fields in Northeastern Alberta
The magnitude of the calculated geoelectric field and estimated line voltage varies significantly depending on whether the 1-D Model-Space Method or Data-Space Method is used. In general, the measured impedances result in smaller geoelectric field estimates over much of the province but an important exception is the area around Fort McMurray in northeastern Alberta which tends to have a stronger geoelectric field when using the Data-Space approach (e.g., Figure 6) as well as the largest difference in line voltage (>100 V) compared to the 1-D Model-Space method (Figure 7). This area has a very resistive surface layer as it marks the edge of the exposed Paleoproterozoic and Archean basement rocks of the Canadian Shield (Figure 1a; Ross & Eaton, 2002) and previous studies indicate that the subsurface structure may be complicated and potentially electrically anisotropic (Liddell et al., 2016). The 1-D model for Zone 7 from Trichtchenko et al. (2019) has a 1.4 km thick 10 Ωm shallow layer representing sedimentary rocks related to the WCSB. The limitation of the one-dimensional piecewise continuous model is that this thickness is assigned for all of Zone 7 even though the WCSB varies in thickness from 0 m in the northeast of Zone 7 to more than 2,000 m in the southwest of Zone 7 near Peace River. This lateral variation in thickness is reliably represented by the measured MT impedance calculation that predicts large geoelectric fields in the northeast and smaller geoelectric fields in the southwest where the conductive sedimentary rocks of the WCSB are thickest (e.g., Figure 1b). This is also apparent in the 3-D modeling exercise described in Section 4.1 (see Figure 8d) with larger geoelectric fields in northeastern Alberta compared to west-central Alberta although the effect is not as dramatic as in the real data. This discrepancy between the 3-D model and the measured MT data suggests that perhaps the conductivity contrast between the Precambrian basement and the WCSB is larger than what was used for the 3-D model in Figure 8. It is also worth noting that the transmission lines with the smallest line voltages calculated using the Data-Space method relative to the 1-D Model-Space method (e.g., Lines 1, 2, 4, 5, and 11) are in the southwestern portion of Zone 7. This is where the thickness of the conductive sedimentary basin is underestimated by the piecewise-continuous 1-D model. Thus, the geoelectric fields are much weaker when using the measured MT impedances which are sensitive to the lateral variations in the thickness of this shallow conductor. The area around Fort McMurray contains significant electrical and pipeline infrastructure related to oil sands mining. It is an area which might be most susceptible to GIC impacts that also sees the greatest differences between the 1-D Model-Space Method and the Data-Space Method.

The Geoelectric Field in British Columbia
The accuracy of the absolute magnitude and direction of the geoelectric field may be susceptible to poor MT data quality and/or distortion effects (Bonner & Schultz, 2017). However, the fact that trends in the geoelectric field are consistent across many sites in Alberta would suggest that the differences in magnitude are due to real structures on the length-scale of interest and are not due to small-scale distorting heterogeneities as such features would be expected to be site-specific. Southeastern British Columbia is not discussed in detail in this paper primarily because there is a significant amount of variability in data quality and coverage in that region without clear trends in the direction or magnitude of the geoelectric field when using the measured impedances. Many MT sites may be influenced by small-scale distorting heterogeneities associated with the complex geology of the Rocky Mountain Trench (Finley, 2020) and the complicated topography of southeastern British Columbia. This would make interpolation between sites prone to significant errors. However, it is interesting to note that the edge of the SABC, which is hypothesized to terminate in southeastern British Columbia (Rippe et al., 2013), may be partially responsible for the complex pattern of scattered geoelectric fields due to edge effects. Broadly, the differences between British Columbia and Alberta illustrate the influence of the thin lithosphere beneath British Columbia. In particular, the magnitude of the geoelectric field is smaller at long periods in British Columbia reflecting the influence of the more conductive shallow asthenosphere whereas in Alberta, at similar periods, the geoelectric field is larger since it is being influenced by the thicker, resistive lithosphere (see Figure 8 and Figures S3 and S4 in Supporting Information S1). This is somewhat contrary to the expectation that geoelectric fields are smaller in sedimentary basins and larger in resistive areas; the sedimentary basin thickness may be negligible at long-periods which are sensitive to much deeper lithospheric structures. Whether this influences time-domain geoelectric fields in British Columbia is beyond the scope of this study and a denser MT array would be required.

Conclusion
This study has shown that estimates of the geoelectric field caused by a GMD will depend on the assumptions made about Earth's conductivity structure even in sedimentary basins and cratonic environments located far from oceans and associated coast effects. The 1-D Model-Space Method assumes 1-D structure a priori, whereas the Data-Space method makes no such assumptions and allows for the possibility of 2-D or 3-D structures to influence the geoelectric field. The results of the Data-Space Method were compared to the 1-D Model-Space Method to highlight differences between the methods. This was first shown using an idealized geomagnetic disturbance to isolate any spatio-temporal heterogeneity in the source field such that all variation in the estimated geoelectric field was due to the Earth's conductivity structure (Figure 1). It was then shown that the measured MT impedances result in a geoelectric field which can be preferentially oriented in a particular direction in the frequency domain even for a real GMD over large regional areas (Figures 2 and 3). This suggests that the measured MT impedance data are sensitive to 3-D geological structures which are not reproduced by the 1-D approximation. This partial polarization can be observed even when the estimated geoelectric field is transformed into the time domain whereby the geoelectric field primarily oscillates in one direction (Figures 4 and 5). Finally, differences in the magnitude and direction of the geoelectric field when using the 1-D Model-Space Method and the Data-Space approaches can ultimately have an impact on the estimated voltage in a transmission line (Figures 6 and 7). The strength of the estimated voltage on a particular transmission line could be 10%-20% larger when using the Data-Space Method depending on the orientation of the line relative to the geoelectric field polarization. The prefer-ential direction of geoelectric field polarization in southern Alberta would not be identified by an analysis which uses 1-D conductivity models.
The cause of the geoelectric field behavior in Alberta and the reasons for the departure from the 1-D assumption were evaluated using 3-D modeling tests (Figure 8). Conductivity models were based on inversions of 3-D MT datasets and indicated that the primary cause of the geoelectric field distortion in southern Alberta was neither (a) changes in lithospheric thickness variations from the Cordillera to the North American craton nor (b) changes in the thickness of sedimentary rocks in the WCSB. Instead, it was shown that the primary cause was the presence of a previously reported lower crustal conductor called the SABC which trends northeast-southwest across Alberta and into southeastern British Columbia. In northeastern Alberta, the large geoelectric field was attributed to the surface contact between the high resistivity Precambrian basement rocks of the Canadian shield and the low resistivity sedimentary rocks of the WCSB that thickens to the southwest. These results suggest that measured MT impedance data may provide a more reliable way to incorporate geological information (and the resulting effects of 3-D electromagnetic induction) into the geoelectric field calculation than the 1-D Model-Space Method. Ultimately, this suggests that it may be beneficial to include 3-D crustal conductivity information from MT impedances in the GIC modeling workflow and this study provides a link between tectonic history and modern-day geoelectric hazards.
Future work could expand upon this study by calculating the estimated GIC in the AIES using the 3-D modeling approach and a power network model which incorporates transformer grounding and line resistances (e.g., the Nodal Admittance Matrix method; Marshall et al., 2019). This would allow for a complete transformer risk assessment which could be compared to e.g., Haque et al. (2017). The primary purpose of our analysis is to show that: (a) the geoelectric field and resulting line voltages can be influenced directly by 3-D crustal-scale structures; (b) there are important differences between the results of the widely used 1-D Model-Space Method and Data-Space Method which lead to discrepancies in estimated line voltages even in sedimentary basin and cratonic environments; and (c) that measured MT impedances provide important additional information about lateral variations which are not obtained by an analysis which assumes a 1-D conductivity structure. Measuring the GIC in real-time during a GMD at a transformer or by using remote fluxgate magnetometer measurements along transmission lines (e.g., Hübert et al., 2020) would provide an important validation to the modeling methodologies presented. We would like to thank David Boteler and Larisa Trichtchenko for providing information about the one-dimensional piecewise continuous model. Thank you to Greg Lucas for early conversations about the project. We would like to thank the many researchers, students and field technicians associated with the University of Alberta who have collected magnetotelluric data in Alberta and British Columbia over many years. Thank you to Chris Thomas at the Alberta Electric System Operator (AESO) for providing information on transmission line routings. We also acknowledge the National Science and Engineering Research Council of Canada (NSERC) support through Discovery Grants to Martyn Unsworth and Ian Mann, and support for data collection from the Canada Foundation for Innovation (CFI), Alberta Geological Survey, Geological Survey of Canada (GEMS) and the Helmholtz Alberta Initiative. CARISMA is operated by the University of Alberta, funded by the Canadian Space Agency. Special thanks to two anonymous reviewers and editor Jennifer Gannon for helpful comments in improving the manuscript during the review process.