Real-Time 3-D Modeling of the Ground Electric Field Due To Space Weather Events. A Concept and Its Validation

We present a methodology that allows researchers to simulate in real time the spatiotemporal dynamics of the ground electric field (GEF) in a given 3-D conductivity model of the Earth based on continuously augmented data on the spatiotemporal evolution of the inducing source. The formalism relies on the factorization of the source by spatial modes (SM) and time series of respective expansion coefficients and exploits precomputed GEF kernels generated by corresponding SM. To validate the formalism, we invoke a high-resolution 3-D conductivity model of Fennoscandia and consider a realistic source built using the Spherical Elementary Current Systems (SECS) method as applied to magnetic field data from the International Monitor for Auroral Geomagnetic Effect network of observations. The factorization of the SECS-recovered source is then performed using the principal component analysis. Eventually, we show that the GEF computation at a given time instant on a 512 × 512 grid requires less than 0.025 s provided that GEF kernels due to pre-selected SM are computed in advance. Taking the 7–8 September 2017 geomagnetic storm as a space weather event, we show that real-time high-resolution 3-D modeling of the GEF is feasible. This opens a practical opportunity for GEF (and eventually geomagnetically induced currents) nowcasting and forecasting.


Governing Equations in the Frequency Domain
We start with the discussion of the problem in the frequency domain. Maxwell's equations govern electromagnetic (EM) field variations, and in the frequency domain, these equations are read as where μ 0 is the magnetic permeability of free space; ω is angular frequency; j ext (r, ω) is the extraneous (inducing) electric current density; B(r, ω; σ) and E(r, ω; σ) are magnetic and electric fields, respectively. σ(r) is the spatial distribution of electrical conductivity, r is a position vector, either in the spherical or Cartesian coordinates. Note that we neglected displacement currents and adopt the following Fourier convention Note that we will use the same notation for the fields in the time and frequency domain. We also assume that the current density, j ext (r, ω), can be represented as a linear combination of spatial modes (SM) j i (r), Note also that the form of j i (r) (and their number, L) varies with application. For example, j ext (r, ω) is parameterized via spherical harmonics in Püthe and Kuvshinov (2013), Honkonen et al. (2018), Guzavina et al. (2019), and Grayver et al. (2021), current loops in Sun et al. (2015), or eigenmodes from the Principal component analysis (PCA) of the physics-based models in Egbert et al. (2021) and Zenhausern et al. (2021).
By virtue of the linearity of Maxwell's equations with respect to the j ext (r, ω), we can expand the total (i.e., inducing plus induced) electric field as a linear combination of individual fields E i , where the E i (r, ω; σ) is the "electric" solution of the following Maxwell's equations:

Governing Equations in the Time Domain
The transformation of Equation 5 into the time domain leads to the representation of the time-varying GEF as a sum of convolution integrals or equivalently where r s stands for the position vector at the surface of the Earth. The reader is referred to Appendix A for more details on the convolution integrals in Equations 8 and 9.
Since the radial component of the GEF is negligibly small (due to insulating air) and is not used in GIC calculations (Kelbert, 2020), we will confine ourselves to modeling the horizontal electric field solely; thus, hereinafter, E i will stand for E i = (E x,i E y,i ).

Real-Time Modeling of the GEF. A Concept
Equation 9 shows how the GEF can be modeled using continuously augmented data on the temporal evolution of the nowcasted or forecasted c i (the reader is referred to Section 4.2 where nowcasting/forecasting of c i is briefly discussed). To make the formula ready for implementation, one needs: (a) to specify a set of SM, j i , i = 1, 2, …, L in the region, where GIC nowcasting/forecasting is required; we will discuss the construction of j i in Section 3.1; (b) to set up a 3-D conductivity model in this region; and (c) to estimate an upper limit of integrals in Equation 9, or, in other words, to estimate the integration time interval, T, above which E i (r s , τ; σ) becomes negligibly small. The latter will allow us to rewrite Equation 9 as Note that the upper limit in the integrals could be different for different SM, different components of the field, and different locations. However, one can choose a conservative approach, taking a single T irrespective of modes/ components/locations as a maximum from all individual upper limit estimates. We will discuss the estimation of T in Sections 3.3 and 3.4.
The details of numerical calculation of the integrals in Equation 10 are presented in Appendix B. In short, assuming that c i (t), i = 1, 2, …, L are time series with the sampling interval Δt and T = N t Δt, we approximate E(r s , t k ; σ) at t k = kΔt as where d i is defined as The reasoning to represent the time-dependent part in Equation 11 in this form is given in Appendix B. Note also that quantities  ( ; ) and  ( , ; ) are time-invariant, and for the given j i , i = 1, 2, …, L and 3-D conductivity model are calculated only once, then stored and used, when the calculation of E(r s , t k ; σ) is required. Actual form and estimation of kernels  ( ; ) and  ( , ; ) are also discussed in Appendix B. Equation 11 is the essence of the real-time GEF calculation, showing that  ( × × ) summations and multiplications are required to compute the GEF at a (current) time instant t k plus some overhead to read the precomputed  ( ; ) and  ( , ; ) from the disc. Note that N g is a number of points r s , at which the GEF is computed.

Real-Time Modeling of the GEF. Validation of the Concept
The validation of the presented concept will be performed using Fennoscandia as a test region. The choice of Fennoscandia is motivated by several reasons. First, it is a high-latitude region, where GIC are expected to be especially large. Second, there exists a 3-D electrical conductivity model of the region (Korja et al., 2002). Third, the regional magnetometer network (International Monitor for Auroral Geomagnetic Effect (IMAGE; Tanskanen, 2009)), allows us to build a realistic model of the source. Finally, the last but not the least consideration to choose this region is the fact that we have already performed a comprehensive 3-D EM model study in this region (Marshalko et al., 2021).

Building a Model of the Source
First, let us rewrite Equation 4 in the time domain We will further assume that the extraneous current j ext (r, t) is divergence-free, it flows in a thin layer at the altitude of h = 90 km, and this layer is separated from the Earth by the insulating atmosphere. Following the Spherical Elementary Current Systems (SECS) method (Vanhamäki & Juusola, 2020), this current is represented as where δ is Dirac's delta function, e ϑ and e φ are unit vectors of the spherical coordinate system, r = (R, ϑ, φ), r m = (R, ϑ m , φ m ), R = a + h, a is a mean radius of the Earth, r m is the location of the pole of the mth spherical elementary current system and S m is the so-called scalar factor associated with the mth pole. Expressions for P(r, r m ) and Q(r, r m ) are presented in Appendix D. Note that in practice r and r m are usually taken as the nodes of two (similar) grids, which are slightly shifted with respect to each other (the reason for the shift is explained in Appendix D). Once S m (t), m = 1…M are obtained by means of the SECS method as applied to some real data for some event, one can perform the PCA of S m (t) expecting that the spatial structure of S m (t) will be well approximated with a small number of principal components (PC) , i = 1, 2, …L allowing to represent j i as The aim of this section is to obtain and, consequently, j i (using Equation 15). To this end, we apply the SECS method to 10-s vector magnetic field data from all available (  The PCA of S m (t) is performed in a similar manner as it was done, for example, in Alken et al. (2017), Egbert et al. (2021), and Zenhausern et al. (2021). Specifically, we construct a matrix F as where is S m (t) estimated at the nth time instant at the mth grid point. Further, according to the PCA concept, and apply an eigenvalue decomposition to R where Λ is a diagonal matrix containing the eigenvalues λ i , i = 1, 2, …, M of R. The at M grid points is represented by ith column vector of V which is in its turn the eigenvector of R corresponding to the eigenvalue λ i . Both PC are usually sorted in order from the largest to the smallest eigenvalues. The PC corresponding to the largest eigenvalue will explain the most variance, followed by the second, third PC, etc… In practice, the PC corresponding to a few of the largest eigenvalues explain most of the analyzed fields' variance. The cumulative variance of L PC can be calculated as (Alken et al., 2017 Figure 2 presents the cumulative variance for the first 30 PC. Horizontal dashed line allows us to estimate the number of modes needed to explain 99% of the spatial variability of S m (t). It is seen from the figure that one needs L = 21 PC to explain most (99%) of the variance. This is a dramatic reduction from the original M = 3,456 SECS poles. These 21 PC will be used in the further discussion of the real-time calculation of the GEF. Figure 3 shows j i corresponding to PC of different i, illustrating the fact that the modes with larger i capture smaller spatial structures of the source. The respective time series c i are presented in Figure 4. Figure 5 compares the maps of the original and the PCA-based sources (i.e., external equivalent currents) for two snapshots of the enhanced geomagnetic activity. The original source is built using the SECS method (Equation 14), whereas PCA-based source is calculated using Equations 13 and 15. It is seen that the agreement between the original and PCA-based sources is very good both in terms of the amplitude and spatial pattern. In addition, Figure 6 demonstrates the comparison of the time series of these sources for two exemplary sites (shown in Figure 5 as white circles): one is located in the region where the significant source current is observed (Jäckvik (JCK)), another-aside from this region (Tartu (TAR)). Again, we observe good agreement between the two sources, especially for the site above which the source current is large.

3-D Conductivity Model of Fennoscandia
We 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 the original data, and then performed the interpolation to a laterally regular grid. Note that a similar procedure was invoked to convert j i from spherical to Cartesian coordinates. The lateral discretization and the size of the resulting 3-D part of the conductivity model of Fennoscandia were taken as 5 × 5 km 2 and 2,550 × 2,550 km 2 , respectively. Deeper than 60 km, we used the 1-D conductivity profile obtained by Kuvshinov et al. (2021) (Figure 7d), which is an updated version of the 1-D profile from Grayver et al. (2017).
Note that the lateral discretization and the size of the conductivity model of Fennoscandia imply that the GEF is calculated at a grid comprising N g = 512 × 512 points.

Computation of E i (r s , ω; σ)
To implement the real-time modeling concept one needs-as is seen from Equations B13 and C2-to compute E i (r s , ω; σ) at a number of frequencies, or, in other words, to solve Maxwell's equations (Equations 6 and 7). These equations are numerically solved using the 3-D EM forward modeling code PGIEM2G , which is based on a method of volume integral equations with a contracting kernel (Pankratov & Kuvshinov, 2016). PGIEM2G exploits a piece-wise polynomial basis; in this study, PGIEM2G was run using the first-order polynomials in lateral directions and third-order polynomials in the vertical direction.  Finally, it is important to note that E i decrease-irrespective of the mode and location-as frequency decreases; specifically, the magnitude of E i drops down more than two orders of magnitude as frequency decreases from 1 Hz down to 10 −3 Hz. These plots suggest a value for T in Equation 10; recall, that useful rule of thumb is that the value for T corresponds to the inverse of frequency at which the field becomes small compared to the higher frequencies. Following this rule, T = 1,000 s seems to be a reasonable choice which will be further justified in the next section.

Model Study to Justify a Value for T
To justify the value for T we first calculate the reference ("true") time-domain electric field for the chosen 8-hr event. This reference field was computed using the procedure presented in Ivannikova et al. (2018) and Marshalko et al. (2020Marshalko et al. ( , 2021. Specifically, we calculate j ext (t, r) using Equations 13 and 15 and taking 21 terms in Equation 13. Further, according to Marshalko et al. (2021), we calculate the reference electric field as follows: 1. The source j ext (t, r) is transformed from the time to the frequency domain with a fast Fourier transform (FFT). 2. Frequency domain Maxwell's equations 1 and 2 are numerically solved using PGIEM2G at FFT frequencies between 1 and 1 2Δ where K is the length of the event, and Δt is the sampling rate of the considered time series. In this study, Δt is 10 s, and K is 8 hr.
3. E(t, r) is obtained with an inverse FFT of the frequency domain field. Then we calculate electric fields using Equation 11 with T = 900 s (15 min) and with T = 3,600 s (1 hr) and compare them with the reference field. Figures 11-13 show comparison of electric field time series, again, at locations of ABK, UPS, and SPG observatories. It is seen that both "real-time" (either calculated using T = 15 min or T = 1 hr) electric fields agree well with the reference electric field. Tables 1 and 2 confirm this quantitatively by presenting correlation coefficients between corresponding time series and the normalized root-mean-square errors; the latter are defined as where a and b are the GEF time series calculated exploiting real-time scheme and the reference GEF time series, respectively, a i and b i are elements of these time series, and N is the number of time instants. Since results for T = 15 min and T = 1 hr appear to be very similar, we present in the next section the estimates of computational loads for the case when T is taken as 15 min.

Computational Loads for the Real-Time GEF Calculation
Once  ( ; ) and  ( , ; ) are computed and stored on the disc, GEF at a grid = × and time instant t k is computed using Equation 11. In accordance with this equation, the GEF calculation requires forecasting/nowcasting the L × N t array c, reading the L × N t × N g array  and L × N g array  , and performing  ( × × ) summations and multiplications. For our problem setup with N g = 512 × 512, N t = 90 and L = 21, the calculation of E(r s , t k ; σ) takes from 0.00625 to 0.025 s, depending on the computational environment. Note that to store arrays for this setup one needs 7.25 Gigabytes of disc space.

Further Justification of the Concept
So far, we have demonstrated the concept's validity on an example of a single space weather event. However, one can argue that j i (r) obtained for a specific event could be non-adequate for other events. To address this question, we performed the following modeling experiment. First, we built the "sources" (i.e. , where ( ) are taken from the PCA of 7-8 September 2017 data, but S m (t) are formed using new event data. To ensure that the spatiotemporal structure of the source for new events is different from that of the (reference) 7-8 September 2017 event, we took the new events' lengths as 48 and 72 hr, respectively; recall that the duration of 7-8 September 2017 event was taken as 8 hr. Figure 14 shows snapshots of the original and SM-based sources for St. Patrick's Day (top panels) and Halloween (bottom panels) geomagnetic storms. It is seen from the figure that the SM-based source (with SM obtained from another event) approximates very well the source of the other two events. Figure 15 confirms this inference by showing the agreement between the time series of the original and SM-based sources above exemplary site JCK, again, for Halloween (top panels) and St. Patrick's Day (bottom panels) storms. These results suggest that irrespective of the event (which corresponds to sources of different geometry), the spatial structure of these sources is well approximated by a finite number of SM obtained from the analysis of some specific event. The prerequisite to getting adequate set of SM is that the event to be used for SM estimation should be long enough and sufficiently energetically large and spatially complex.
The linked question we also address is whether T = 15 min is a valid choice for the real-time modeling of the GEF during the above-discussed events. As in Section 3.4, we calculate electric fields using Equation 11 with T = 15 min and with T = 1 hr and compare them with the reference fields. Two top panels in Figures 16 and 17 show the comparison of electric field time series at the location of the ABK observatory for Halloween and St. Patrick's Day events, respectively. Similar to the 7-8 September 2017 event results, both "real-time" electric fields agree well with the reference electric field. Besides, bottom panels in Figures 16 and 17 show the difference between "real-time" electric fields calculated using T = 15 min and T = 1 hr. It is seen that this difference is small compared to signals themselves (cf., top panels in the same figures), once again confirming the fact that T = 15 min can be used for the real-time GEF modeling.
Tables 1 and 2 quantify the agreement for these two events at three geomagnetic observatories, as it was already done for 7-8 September 2017 event. It is seen that the agreement between results for two new events is as good as for the 7-8 September 2017 geomagnetic storm. Notably, throughout all three events, the correlation coefficient is lower, and nRMSE is larger for ABK than for UPS and SPG. This is probably because ABK is located in the region with a large lateral contrast of conductivity (Figure 7), where modeling results are expectedly less accurate.

Nowcasting and Forecasting GEF Using the Proposed Concept
In this section, we discuss how the proposed concept could be implemented for nowcasting/forecasting of the GEF. Specifically, a scheme to nowcast GEF could work as follows: 1. Using magnetic field data collected at an observational network for historical (past) event/several events one obtains at r m . This is done by exploiting the procedure described in Section 3.1. These allow us to represent the source at any time instant t and at any position r via Equations 13-15. In this paper, we used IMAGE network of magnetic field data to obtain L = 21 ( ) and further j i (r). Using IMAGE data, we confine ourselves to Scandinavian region. If Canada, for example, is a region of interest, one would use the data from the Canadian networks of magnetic field observations, like CARISMA (Mann et al., 2008) and AUTUMNX (Connors et al., 2016). 2. Once ( ) and, subsequently, j i (r) are obtained and stored, one estimates electric field at the current time instant, t k , using Equation 11. This, in particular, requires knowledge of coefficients c i at time instant t k and at a number of time instants in the past, t k − Δt, t k − 2Δt, …, t k − N t Δt. The coefficients at these instants are obtained by reusing Equations 16 and 19, namely.
where S m are computed from the available ground magnetic field data. Note that N t Δt = T, where T is either 15 min or 1 hr in our example. It is also important to stress that to model the GEF at the next time instant, t k + Δt, one needs to update only c i (t k + Δt).
As for forecasting of the GEF, the scheme could include the following steps: 1. One obtains ( ) and, subsequently, j i (r) in a similar manner as it is done in case of the GEF nowcasting. 2. One trains the neural network (NN) using as input data the time series of solar wind parameters collected by satellite(s) at the L1 Lagrangian point and c i (t) as output data. Time series c i (t) for the training period are obtained using Equations 16 and 19. There is a common understanding that the longer time series are used for the training phase, the better the quality of the forecasted results. Therefore, this period preferably should include multiple years of the L1 and ground magnetic field data; recall that c i (t) during the training phase are obtained from the ground magnetic field data. 3. One forecasts GEF using the trained NN. Ideally, one has to forecast well ahead. However, given observations made at the L1 point, a geomagnetic disturbance is seen on the ground as fast as an hour ahead. This time latency can be further shrunk to half an hour or so, depending on the solar wind speed. This, in particular, advocates real-time modeling of the GEF which is a topic of this paper.

Conclusions
In this paper, we presented a formalism for real-time computation of the GEF in a given 3-D Earth's conductivity model excited by a continuously augmented spatially and temporally varying source responsible for a space weather event.
The formalism relies on the factorization of the source by SM and time series of respective expansion coefficients, and exploits precomputed GEF kernels generated by corresponding SM.
To validate the formalism, we invoked a high-resolution 3-D conductivity model of Fennoscandia and considered a realistic source built with the use of the SECS method as applied to magnetic field data from the IMAGE network of observations. Factorization of the SECS-recovered source is then performed using the PCA. Eventually, we show that the GEF computation at a given time instant on a 512 × 512 grid requires at most 0.025 s provided that GEF kernels due to the pre-selected SM are computed in advance. This opens a practical opportunity for GEF nowcasting, using ground magnetic field data, or even forecasting, using both ground magnetic field and L1 data.
We illustrate the concept on a Cartesian geometry problem setup. Global-scale implementation is rather straightforward; for this scenario, the source could be obtained either using magnetic field data from the global network of geomagnetic observatories or exploiting the results of the first-principle modeling of the global magnetosphere-ionosphere system.

Appendix A: Properties of Transfer Functions and Impulse Responses
The convolution integrals in Equation 9 represent the response of the medium to a time-varying extraneous current. These relations follow from the properties of a physical system we consider. We list these properties below and discuss implications. The presentation closely follows a more detailed analysis by Svetov (1991). Note that for the sake of clarity, we discuss the properties on an example of abstract scalar quantities and omit their dependence on spatial variables and electrical conductivity pertinent to our application.
1. Linearity allows us to define the response, ζ(t), of the medium at time t to an extraneous forcing as where χ is the extraneous forcing that depends on time t′ and  ( ′ ) is the medium Green's function.
2. Stationarity implies that the response of the medium does not depend on the time of occurrence of the excitation. In this case  ( ′ ) ≡ ( − ′ ) and Equation A1 is rewritten as a convolution integral where f(t) represents the impulse response of the medium. In the frequency domain, the convolution integral degenerates tõ( Figure 13. The same caption as in Figure 11 but for the Saint Petersburg (SPG) geomagnetic observatory.
where ̃( ) is called the transfer function and we use tilde sign (⋅) to denote complex-valued quantities. Equations A2 and A3 are related through the Fourier transform̃( 3. Since we work in the time domain with a real-valued forcing, the impulse response is also real. To see implications of this, let us define the inverse Fourier transform of ̃( ) = ( ) + i ( ) as For an impulse response to be real, the last term in the integral (Equation A5) has to vanish. This is possible only if f R (ω) and f I (ω) are even and odd functions of ω, respectively. Therefore, Equation A5 reduces to 4. Impulse response is causal. This property implies that Under this assumption, the convolution integral (Equation A2) is recast to Also, due to causality (Equation A7), and exploiting Equation A6, one can write for t < 0 Further, using the fact that cos(ωt) and sin(ωt) are odd and even functions with respect of t, one obtains for t > 0 The integrals  ( , ; ) can be computed using the digital filter technique (see Appendix C), whereas first term in the RHS of Equation B5 is estimated as follows.
Taking into account that we have c i (t) at discrete time instants, t = nΔt, n = 0, 1, … , we approximate d i (t, τ; T) using the Whittaker-Shannon (sinc) interpolation formula Recall that sinc interpolation is a method to construct a continuous band-limited function from a sequence of real numbers, in our case time series d i at time instants t = nΔt, n = 0, 1, …. Note that in our context, the term "band-limited function" means that non-zero values of a Fourier transform of this function are confined to the frequencies Using the approximation (Equation B7) and taking into account that E i (r s , τ; σ) = 0, τ < 0 (Appendix A), one obtains Further, following the properties of the Fourier transform as applied to sinc function, we obtain that Integrals in Equation C2 can be efficiently estimated using the digital filter technique. Specifically, one needs to construct a digital filter for the following integral transform