Comparing Three Approaches to the Inducing Source Setting for the Ground Electromagnetic Field Modeling due to Space Weather Events

Ground‐based technological systems, such as power grids, can be affected by geomagnetically induced currents (GIC) during geomagnetic storms and magnetospheric substorms. This motivates the necessity to numerically simulate and, ultimately, forecast GIC. The prerequisite for the GIC modeling in the region of interest is the simulation of the ground geoelectric field (GEF) in the same region. The modeling of the GEF in its turn requires spatiotemporal specification of the source which generates the GEF, as well as an adequate regional model of the Earth’s electrical conductivity. In this paper, we compare results of the GEF (and ground magnetic field) simulations using three different source models. Two models represent the source as a laterally varying sheet current flowing above the Earth. The first model is constructed using the results of a physics‐based 3‐D magnetohydrodynamic (MHD) simulation of near‐Earth space, the second one uses ground‐based magnetometers’ data and the Spherical Elementary Current Systems (SECS) method. The third model is based on a “plane wave” approximation which assumes that the source is locally laterally uniform. Fennoscandia is chosen as a study region and the simulations are performed for the September 7–8, 2017 geomagnetic storm. We conclude that ground magnetic field perturbations are reproduced more accurately using the source constructed via the SECS method compared to the source obtained on the basis of MHD simulation outputs. We also show that the difference between the GEF modeled using laterally nonuniform source and plane wave approximation is substantial in Fennoscandia.

In spite of the fact that the importance of performing simulations using fully 3-D conductivity models is currently widely recognized (Kelbert, 2020), such simulations are still rather rare in the GIC community (e.g., Liu et al., 2018;Marshalko et al., 2020;Marshall et al., 2019;Nakamura et al., 2018;Pokhrel et al., 2018;Wang et al., 2016), mostly due to the lack of the credible 3-D conduc-MARSHALKO ET AL.   tivity models of the regions of interest as well as unavailability of adequate tools to model the problem in the full complexity.
As for the source, approximating it by plane waves still prevails in the GIC-related studies (e.g., Campanya et al., 2019;Kelbert & Lucas, 2020;Kelbert et al., 2017;Lucas et al., 2018;Sokolova et al., 2019;Wang et al., 2020). This approximation seems reasonable in low and middle latitudes, where the main source of anomalous geomagnetic disturbances is a large-scale magnetospheric ring current. However, the plane wave assumption may not work in higher latitudes, where the main source of the disturbances is the auroral ionospheric current, which is extremely variable both in time and space, especially during periods of enhanced geomagnetic activity (Belakhovsky et al., 2019). Marshalko et al. (2020) provided some evidence for that by comparing ground EM fields modeled in the eastern United States using the plane wave approximation and the excitation by a laterally variable source which was constructed using outputs from 3-D magnetohydrodynamic (MHD) simulation of near-Earth space. The authors found that the difference increases toward higher latitudes where the lateral variability of the source expectedly enlarges. However, their modeling was mostly confined to the midlatitude region, thus it is still unclear how pronounced the difference between the plane wave and "laterally varying source" results could be in auroral regions. In this paper, we investigate this problem using Fennoscandia as a study region. The choice of Fennoscandia is motivated by: (a) high-latitude location of the region; (b) the availability of the 3-D ground electrical conductivity model of the region (Korja et al., 2002); (c) the existence of the regional magnetometer network (International Monitor for Auroral Geomagnetic Effect, IMAGE (Tanskanen, 2009)) allowing us to build a data-based model of a laterally variable source, which is a natural alternative to physics-based (MHD) source model in the areas with a reasonably dense net of observations. MARSHALKO ET AL.  . Snapshots of the magnitude and direction of the equivalent current computed using the SECS method at an altitude of 90 km above the surface of the Earth at 23:16 and 23:52 UT on September 7, 2017. Locations of IMAGE magnetometers (including Abisko (ABK) and Uppsala (UPS)), the data from which were used for the equivalent current construction, are marked with squares. The location of the Saint Petersburg (SPG) geomagnetic observatory is marked with a circle. Note that SPG is not a part of the IMAGE magnetometers network and its magnetic field data were not used for the equivalent current construction.
Specifically, we perform 3-D modeling of ground electric and magnetic fields in Fennoscandia using three different source models and taking September 7-8, 2017 geomagnetic storm as a space weather event. Two models approximate the source by a laterally varying sheet current flowing above the Earth's surface. One of the models is built using the results of physics-based 3-D MHD simulation of the near-Earth space, another model uses the IMAGE magnetometer data and the Spherical Elementary Current Systems (SECS) method . The third modeling is based on a "plane wave" approximation which assumes that the source is locally laterally uniform. Note that previous GIC-related studies in Fennoscandia operated with both 1-D (e.g., Myllys et al., 2014;Pulkkinen et al., 2005;Viljanen & Pirjola, 2017) and 3-D (Dimmock et al., 2019(Dimmock et al., , 2020 Earth's conductivity models, the magnetic field in most of these papers was allowed to be laterally variable, but the GEF was always calculated implicitly assuming the plane wave excitation.
We compare modeling results and discuss found differences and similarities. We also compare the results of magnetic field modeling with observations. The paper is organized as follows. The methodology used is described in Section 2 followed by the presentation of our results in Section 3. Finally, the discussion of our results and conclusions are given in Section 4. MARSHALKO ET AL.

Governing Equations and Modeling Scheme
We compute electric, E(t, r), and magnetic, B(t, r), fields for a given Earth's conductivity distribution σ(r) and a given inducing source j ext (t, r), where t and r = (x, y, z) denote time and position vector, correspondingly. These fields obey Maxwell's equations, which are written in the time domain as where μ 0 is the magnetic permeability of free space. Note that this formulation of Maxwell's equations neglects displacement currents, which are insignificant in the range of periods considered in this study. We solve Equations 1 and 2 numerically using the following three-step procedure:

Maxwell's equations in the frequency domain
are numerically solved for the corresponding angular frequencies ω = 2πf, using the scalable 3-D EM forward modeling code PGIEM2G , based on a method of volume integral equations (IE) with a contracting kernel (Pankratov & Kuvshinov, 2016).
We would like to note here that in our previous papers (Ivannikova et al., 2018;Marshalko et al., 2020) we used modeling code extrEMe (Kruglyakov et al., 2016) which is also based on the IE method. The distinction between the two codes lies in the different piece-wise bases used. PGIEM2G exploits a piecewise polynomial basis whereas extrEMe uses a piece-wise constant basis. We found that in order to properly account for the effects (in electric field) from extremely large conductivity contrasts in the Fennoscandian region, extrEMe requires significantly larger computational loads compared to the PGIEM2G. This is the reason why we used the PGIEM2G code rather than extrEMe to obtain the modeling results presented in this paper. Specifically, PGIEM2G was run with the use of first-order polynomials in lateral MARSHALKO ET AL.
10.1029/2020SW002657 7 of 18    where L is the length of the (input) times series of the inducing current j ext (t, r), and Δt is the sampling rate of this time series. In this study Δt is 1 min, and L is 8 h.
3. E(t, r) and B(t, r) are obtained with an inverse FFT of the frequency-domain fields.

3-D Conductivity Model
3-D conductivity model of the region was constructed using the SMAP (Korja et al., 2002)-a set of maps of crustal conductances (vertically integrated electrical conductivities) of the Fennoscandian Shield, surrounding seas, and continental areas. The SMAP consists of six layers of laterally variable conductance. Each layer has a thickness of 10 km. The first layer comprises contributions from the seawater, sediments, and upper crust. The other five layers describe conductivity distribution in the middle and lower crust. SMAP covers 0°E-50°E and 50°N-85°N area and has a 5′ × 5′ resolution. We converted the original SMAP database into a Cartesian 3-D conductivity model of Fennoscandia with three layers of laterally variable conductivity of 10, 20, and 30 km thicknesses (Figures 1a-1c). This vertical discretization is chosen to be compatible with that previously used by Rosenqvist and Hall (2019) and Dimmock et al. (2019Dimmock et al. ( , 2020 for GIC studies in the region. Conductivities in the second and the third layers of this model are simple averages of the conductivities in the corresponding layers of the original conductivity model with six layers. To obtain the conductivities in Cartesian coordinates, we applied the transverse Mercator map projection (latitude and longitude of the true origin are 50°N and 25°E, correspondingly) to original data and interpolated the results onto a regular lateral grid. The lateral discretization and size of the resulting conductivity model were taken as 5 × 5 km 2 and 2,550 × 2,550 km 2 , respectively. Deeper than 60 km we used a 1-D conductivity profile obtained by Grayver et al. (2017) (cf. Figure 1d).

EM Induction Source Settings
In this section, we discuss the construction of two models of a laterally variable source and also explain how the EM field is calculated in the framework of the plane wave (laterally uniform source) concept. The sources are set up for the geomagnetic storm on September 7-8, 2017, more specifically, for 8-h time period from 20:00 UT, September 7, 2017 to 03:59 UT, September 8, 2017, thus, before and during the main phase of the storm. The disturbance storm time (Dst) index during this geomagnetic storm reached −124 nT according to the World Data Center of Kyoto (http://wdc.kugi.kyoto-u.ac.jp/dstdir/). More details on the September 2017 storm can be found in Linty et al. (2018) and Dimmock et al. (2019).

Construction of the Source on the Basis of an MHD Simulation
The first source model is based on the results of the physics-based 3-D MHD simulation of near-Earth space. In this study, we employ the Space Weather Modeling Framework (SWMF, Toth et al., 2005Toth et al., , 2012. The input to this MHD model is solar wind (density, temperature, velocity) and interplanetary magnetic field parameters measured at satellites located at L1 Lagrange point, such as Advanced Composition Explorer (ACE) and Deep Space Climate Observatory (DSCOVR). The other input is the solar radio flux at F10.7 cm (2,800 MHz). The outputs are time-varying 3-D currents in the magnetosphere, horizontal currents in the ionosphere, and fieldaligned currents flowing between the magnetosphere and the ionosphere. These output data are then used to calculate (via the Biot-Savart law) external magnetic field perturbations B ext (t, r) at the ground using the Cal-cDeltaB tool (Rastätter et al., 2014).   (5) where δ(z − 0+) is Dirac's delta function, e r is the radial (outward) unit vector, and ∇ ⊥ is the surface gradient. The whole scheme of the equivalent current density calculation from the outputs of MHD simulations is discussed in Ivannikova et al. (2018).
The SWMF run, results of which are used in the current study, was performed at NASA's Community Coordinated Modeling Center (CCMC) at the Goddard Space Flight Center. The version of the SWMF is v20140611. The Rice Convection Model was used to simulate the inner magnetosphere dynamics (Toffoletto et al., 2003). The ionospheric electrodynamics is simulated using the Ridley Ionosphere Model (Ridley et al., 2004). The MHD modeling domain consists of about one million grid cells. The size of the smallest cells is 0.25 R E (where R E is the Earth radius) close to the inner boundary of the modeling domain. The size of the largest cells is 8 R E (close to the outer boundary in the distant tail). The outer boundaries are set at 32 R E in the +x upstream direction, 224 R E in the −x downstream direction, and 128 R E in the ±y and z directions (GSM coordinates). The inner boundary is located at a distance of 2.5 R E from the Earth's center. One-minute OMNI solar wind data were used as an input in this run. The F10.7 cm flux was set to 130.4. Details and results of the run are available at the CCMC website (https://ccmc.gsfc.nasa.gov, run number Naomi_Maruyama_011818_1).
We would like to note that we also performed SWMF simulations with the same input parameters as were used in the CCMC Naomi_Maruyama_011818_1 run, but with different spatial resolutions at the inner boundary of the modeling domain, namely, 0.125 and 0.0625 R E . External magnetic fields (see Figure S1 and Table S1 in the supporting information) from higher-resolution MHD simulations appeared not to differ significantly from those obtained based on Naomi_Maruyama_011818_1 run in the region of our interest. Taking into account that small differences in the external magnetic field should not notably affect modeling results, we construct the "MHD-based" source using Naomi_Maruyama_011818_1 simulation outputs.

Construction of the Source Using the SECS Method
The second model of the source was constructed using the SECS method . In this method, the elementary current systems form a set of basis functions for representing two-dimensional vector fields on a spherical surface. An important application of the SECS method, which is relevant for this study, is the estimation of the ionospheric current system from ground-based measurements of magnetic field disturbances. Note that elementary current systems, as applied to ionospheric current systems, were first introduced by Amm (1997).
With the help of the SECS technique, it is possible to separate the measured magnetic field into external and internal parts, which are represented by two equivalent sheet currents placed in the ionosphere and underground .
To construct the external sheet current, we used IMAGE 10 s vector magnetic field data from all available stations, except for Røst and Harestua, for which the baselines are not yet determined. Baselines are subtracted from variometers' measurements according to the method of van de Kamp (2013). Ionospheric current density is computed using the 2-D SECS method  with the following parameters:  Table 2 The Same Caption as in Table 1 Table 3 The Same Caption as in Table 1  Note that extrapolation of the equivalent current density up to 42°E is performed in order to cover the whole region of Fennoscandia, even though the estimates of the equivalent current far from the stations are less reliable. This applies not only to estimates in areas outside of the region covered by the stations but also to estimates inside of the region covered by the stations at locations where the distances between the nearby station are large. We further perform the equivalent current extrapolation in order to ensure its smooth decay outside the region covered by the data. This is done to avoid the occurrence of artifacts from the edges of the current sheet. We also reduce the temporal resolution of the estimated equivalent current from 10 s to 1 min in order to perform a comparison of modeling results obtained via the MHD-based and SECS-based sources. We then project the current density onto a region of interest and perform vector rotation, which is required for the results' transition from spherical to Cartesian coordinate system. After that we interpolate the current density onto a regular Cartesian grid.

Plane Wave Modeling
The scheme of the GEF calculation via the plane wave approach differs from the one described in Section 2.1. The plane wave modeling results are obtained as follows: 1. 3-D EM forward modeling is carried out via PGIEM2G code  with two (laterally uniform) plane wave sources for the SMAP conductivity model at FFT frequencies corresponding to periods from 2 min to 8 h. 3-D MT impedances Z(ω, r) (Berdichevsky & Dmitriev, 2008) that relate the surface horizontal electric field with the surface magnetic field at each grid point r are then calculated for each FFT frequency ω. 2. We then consider the magnetic field modeled using the PGIEM2G code and the SECS-based source as the "true" magnetic field, thus mimicking the actual magnetic field in the region. 3. Further, the horizontal GEF is calculated for each frequency and each grid point r as 4. Finally, an inverse FFT is performed for the frequency-domain GEF to obtain the "plane wave" GEF in the time domain.

Comparing Results at a Number of Locations in the Region
We first compare modeled and recorded magnetic field variations at the locations of the geomagnetic observa-  Table 4 The Same Caption as in Table 1

but for the Horizontal Electric Field, SECS-Based and Plane-Wave-Based Results and for Three Extra Locations
locations are shown in Figure 1. The sampling rate of the time series is 1 min. The linear trend was removed from observatory data before comparing them to modeling results.
Three upper plots in Figures 4-6 demonstrate time series of (total, i.e., external + induced) magnetic field modeled using MHD-based and SECS-based sources (hereinafter to be referred as MHD-based and SECS-based magnetic fields), as well as time series of the observed magnetic fields. We do not show in these plots "plane wave" magnetic fields because by construction they coincide with the SECS-based magnetic field (see the second item in Section 2.3.3). It is seen that the agreement between SECS-based and observed magnetic fields for ABK and UPS observatories is very good in all components. This is not very surprising because magnetic field data from these observatories were used to construct the SECS source.
As the construction is based on the least-square approach, it inevitably attempts to make predictions close to observations. In this context probably the most interesting comparison is for the SPG observatory because this observatory is not a part of the IMAGE array, and its data were not used for the SECS source construction. Remarkably, the agreement between SECS-based and observed magnetic fields for SPG is also good, except B y component. The disagreement in B y may be due to a localized geomagnetic disturbance which is not accounted for in the SECS source model. Table 1 supports quantitatively the above observations by presenting correlation coefficients between corresponding time series and the normalized root-mean-square errors which is defined as where a and b are modeled and observed time series, respectively, a i and b i are elements of these time series, and n is the number of these elements.
It is also seen from the figures and Table 2 that the agreement between MHD-based and observed magnetic field is significantly worse for all considered observatories and all components. The agreement is especially bad in the B y component. On the whole, the magnitude of MHD-based magnetic field perturbations is underestimated (compared to the SECS-based and observed magnetic field perturbations). Moreover, MHDbased magnetic field captures less of the short-period variability. These results are consistent with results of Kwagala et al. (2020), who carried out SWMF simulations for a number of space weather events and compared SWMF-based (external) magnetic fields with observed ones at a number of locations in northern Europe. According to their modeling results, the SWMF predicts the northward component of external magnetic field perturbations better than the eastward component in auroral and subauroral regions, which is also the case in our modeling of the total magnetic field. As it was mentioned by Kwagala et al. (2020), poor prediction of the eastward component of magnetic field perturbations is directly related to the northward current density in the ionosphere and may arise from the misplacement of the currents in the SWMF with respect to the magnetometer stations.
Finally, two lower plots in Figures 4-6 show plane-wave-based, SECS-based, and MHD-based horizontal GEF. Note that long-term continuous observations of GEF are absent in the region, thus only modeling results are shown in the plots.
Similar to the MHD-based magnetic field, the MHD-based GEF is underestimated compared to the SECSbased GEF. The correlation between these modeling results is very low and nRMSE are high (see Table 3).
On the contrary, SECS-based and plane-wave-based electric fields are rather close to each other, especially at locations of UPS and SPG observatories; Table 4 illustrates this quantitatively. Correlation between modeling results at ABK observatory is lower and nRMSE is higher most likely due to the fact that this observatory is situated in the region with high lateral conductivity contrasts (resistive landmass and conductive sea). To put more weight on this inference last three columns of Table 4 and Figure 7 demonstrate SECS-based and plane-wave-based results for three "sites" also located in the regions with high lateral conductivity contrasts (their locations are shown in Figure 1). Now we observe that the difference between the results is even more pronounced which is, in particular, reflected in lower correlation coefficients and higher nRMSE.

Comparing Results in the Entire Region
Contrary to the previous section where we compared modeled and observed time series of the EM field at a number of locations, in this section, we compare MHD-based, SECS-based, and plane-wave-based electric fields in the entire region for two time instants discussed earlier in the paper.
Top and middle plots in Figure 8 show magnitudes of respective SECS-based and MHD-based GEF. Bottom plots show the absolute differences between corresponding GEF magnitudes. It is seen that SECS-based GEF significantly exceeds MHD-based GEF throughout the Fennoscandian region and for both time instants. The largest differences occur in areas of strong lateral contrasts of conductivity (e.g., at the coastlines) and at higher latitudes.
In a similar manner, Figure 9 presents the comparison of SECS-based and plane-wave-based GEF. In contrast to MHD-based results, at a first glance magnitude of plane-wave-based GEF is in overall comparable with SECS-based GEF (cf. top and middle plots in the figure). However, bottom plots show that the difference is substantial but more localized (compared to the difference between SECS-based and MHD-based results), occurring, again, in areas of strong lateral contrasts of conductivity and increasing toward higher latitudes.

Conclusions and Discussion
In this work, we performed 3-D modeling of the EM field in the Fennoscandian region during the September 7-8 geomagnetic storm in 2017. The goal of this model study was to explore to what extent the resulting EM field depends on the setup of the external source which induces this field. We have used three different approaches to the EM induction source setting. The first technique is based on the retrieval of the (laterally variable) equivalent current from the dedicated MHD simulation. In the second method, the laterally variable equivalent current is constructed on the basis of IMAGE magnetometers' data using the SECS approach. The third technique exploits the plane wave concept, which implies that the source is laterally uniform locally.
Two main findings of the paper are as follows. Magnetic field perturbations in Fennoscandia are reproduced much more accurately using the SECS rather than the MHD-based source, constructed using the SWMF. The difference between the GEF modeled using the SECS-based laterally varying source and the plane wave excitation is substantial in Fennoscandia, especially in the areas of strong lateral contrasts of conductivity (e.g., at the coasts), and at higher latitudes where lateral variability of the source becomes more pronounced.
We would like to remind the reader that in order to obtain MHD-based 3-D EM modeling results presented in this paper we calculated the external magnetic field perturbations on a coarse 5° × 5° grid, which was done to reduce the computational time. Ideally, the resolution of the grid should be much higher to account for the effects of small-scale current structures. However, according to our results the external magnetic field is not reproduced accurately enough at locations of geomagnetic observatories ABK and UPS using the SMWF outputs irrespective of the resolution of the MHD modeling domain (see Figure S1 and Table S1 in the supporting information). That is why increasing the external magnetic field grid resolution most likely will not significantly improve the 3-D EM modeling results, at least in the case of the September 7-8, 2017 geomagnetic storm, the Fennoscandian region, a particular setup of the SWMF described in the current paper and three considered resolutions of the MHD modeling domain. However, it is clear that a separate study is required to investigate the influence of the spatial resolution of the MHD model on the external magnetic field perturbations at the Earth's surface.
From our study, the reader may have a biased impression that the SECS-based current system is an ideal source candidate for rigorous modeling (and eventually forecasting) ground EM field due to space weather events. However, our vision of the problem is that each source setting discussed in this study has its own advantages and drawbacks.
The MHD-based approach is the only one out of the considered three, which allows researchers to forecast the space weather impact on ground-based technological systems. This is possible due to the fact that MHD simulations are run on the basis of the satellite solar wind data collected at the L1 Lagrange point. The solar wind velocity has a typical speed of 300-500 km/s and, thus, the geomagnetic disturbance observed at the Earth's surface is usually lagged compared to the L1 point in the range of 30-90 min (Cameron & Jackel, 2019). This time is, obviously, reduced for fast CMEs; the initial speed of one of the fastest recorded CMEs, which occurred on July 23, 2012 (but was not Earth-directed), was estimated as 2,500 ± 500 km/s (Baker et al., 2013). Another advantage of the aforementioned method is the ability to compute the equivalent current and, subsequently, the EM field for any point on the Earth. It is noteworthy that this method is not dependent on ground-based geomagnetic field observations. The drawback of the approach is that it is currently the least accurate among the considered modeling techniques. Moreover, significant computational resources (in terms of CPU time and memory) are required to carry out MHD simulations. In spite of the fact that these simulations are still rather far from reproducing actual ground geomagnetic disturbances (as is shown once again in this paper) there are continuing efforts to improve the predictive power of MHD models (e.g., Zhang et al., 2019).
The SECS-based approach uses ground magnetometers' data and, thus, does not have forecasting capabilities. However, it is the most accurate among all the considered approaches, but in order to properly capture the spatiotemporal evolution of the source, it requires a dense grid of continuous geomagnetic field observations in the region of interest.
The plane wave method is most probably an optimal choice for EM modeling (due to space weather events) in low-latitude and midlatitude regions provided MT impedances are estimated in these regions on as regular and detailed grid as practicable. The plane wave approach is the least computationally expensive among the three methods considered in this study, as MT impedances can be computed/estimated in advance and then convolved with the magnetic field which, again, requires a network of continuous geomagnetic field observations in the region. However, the violation of the plane wave assumption in high latitudes leads to significant differences between GEF modeled with the use of the SECS-based laterally varying source and the plane wave approximation.

Data Availability Statement
The SWMF model is available from the University of Michigan upon acceptance of license agreement, and SWMF results used here are available at NASA's Community Coordinated Modeling Center (CCMC: https:// ccmc.gsfc.nasa.gov/results/viewrun.php?domain=GM&runnumber=Naomi_Maruyama_011818_1). OMNI solar wind data were used as an input for this run (http://omniweb.gsfc.nasa.gov). We thank Toivo Korja, Maxim Smirnov, and Lisa Rosenqvist for providing the SMAP model. The SMAP model is available via the EPOS portal (http://mt.bgs.ac.uk/EPOSMT/2019/MOD/EPOSMT201_3D.mod.json). We thank the institutes that maintain the IMAGE Magnetometer Array: Tromsø Geophysical Observatory of UiT, the Arctic University of Norway (Norway), Finnish Meteorological Institute (Finland), Institute of Geophysics Polish Academy of Sciences (Poland), GFZ German Research Center for Geosciences (Germany), Geological Survey of Sweden (Sweden), Swedish Institute of Space Physics (Sweden), Sodankylä Geophysical Observatory of the University of Oulu (Finland), and Polar Geophysical Institute (Russia). In this paper, we also used magnetic field data collected at the geomagnetic observatory Saint Petersburg. We thank the Geophysical Center of the Russian Academy of Sciences that supports it and INTERMAGNET (www.intermagnet.org) for promoting high standards of magnetic observatory practice.