Multi‐Site Transfer Function Approach for Real‐Time Modeling of the Ground Electric Field Induced by Laterally‐Nonuniform Ionospheric Source

We propose a novel approach to model the ground electric field (GEF) induced by laterally‐nonuniform ionospheric sources in real time. The approach exploits the multi‐site transfer function concept, continuous magnetic field measurements at multiple sites in the region of interest, and spatial modes describing the ionospheric source. We compared the modeled GEFs with those measured at two locations in Fennoscandia and observed good agreement between modeled and measured GEF. Besides, we compared GEF‐based geomagnetically induced current (GIC) with that measured at the Mäntsälä natural gas pipeline recording point and again observed remarkable agreement between modeled and measured GIC.


10.1029/2023SW003621
2 of 15 (approximated by ionospheric sheet current) using the Spherical Elementary Current System (SECS) approach (Vanhamäki & Juusola, 2020) as applied to the magnetic field variations simultaneously recorded at a regional network of observations.Further, the source is factorized by spatial modes (SM) and respective time-varying expansion coefficients using time-domain Principal Component Analysis (PCA) of the SECS-recovered inducing current.
The SECS approach can be viewed as a regional variant of the Gauss method-the method used to separate the inducing and induced signals on a global scale.If the region of interest is characterized by 3-D conductivity distribution, the induced part is inevitably influenced by 3-D effects arising, particularly from the lateral (e.g., land/ ocean) conductivity contrasts.Given the usually deficient spatial distribution of the magnetic field observations, the SECS approach precludes an accurate description of the induced part affected by localized 3-D effects.Such imperfection in the induced part description also influences the recovery of the inducing part in terms of the recovery of the time series of expansion coefficients.
To circumvent this problem Marshalko et al. (2023) introduced an approach to estimate expansion coefficients more accurately.The approach also exploits the spatial modes resulting from the PCA of SECS-recovered inducing source, but an estimation of the expansion coefficients is performed using a 3-D conductivity model of the region.However, the disadvantage of their method is a high computational load to obtain expansion coefficients, thus hindering its implementation for real-time GEF modeling.This paper proposes a method that also works with spatial modes but avoids estimating the expansion coefficients.An approach is based on the multi-site (MS) transfer function (TF) concept where MS TF relates the GEF at any location with the horizontal components of the magnetic field at (fixed) multiple locations.Notably, an approach allows researchers to model the GEF in real time.The proposed technique can be considered a generalization of the inter-site TF formalism presented in Kruglyakov, Kuvshinov, and Nair (2022).In their paper, TFs relate sea-bed and single-site land-based magnetic fields, and spatial modes are obtained from spherical harmonic analysis of the magnetic data recorded at a global network of geomagnetic observatories.
Using Fennoscandia as a study region, we validate the MS TF approach by comparing modeled and measured GEF and GIC time series.The reasons for selecting this high-latitude region are manifold.First, technological systems in polar regions are especially vulnerable due to the excessive GIC.Second, the inducing source at high latitudes reveals significant lateral variability, which we are interested in accounting for in our approach.Third, the 3-D electrical conductivity model is available for Fennoscandia (Korja et al., 2002).Fourth, there is a magnetometer network in the region (International Monitor for Auroral Geomagnetic Effect, IMAGE (Tanskanen, 2009)), which provides long-term, simultaneous and continuous measurements of the magnetic field at multiple locations-the prerequisite data for successfully applying the proposed approach.Moreover, in Fennoscandia, one can find relatively long (from weeks to months) continuous GEF measurements performed during magnetotelluric surveys in the region.Finally, GIC observations have been carried out at the Mäntsälä natural gas pipeline recording point starting from November 1998 (Viljanen et al., 2006), allowing us to compare observed and modeled GIC for several space weather events.

Methodology
We assume that the horizontal GEF E h = (E x , E y ) at any location is related in the frequency domain to where x and y are directed to geographic North and East, respectively.
A set of Λ k ≡ Λ(r, r k , ω; σ) is called multi-site transfer function, where each Λ k is the following 2 × 2 matrix 10.1029/2023SW003621 3 of 15 Estimating the elements of matrices Λ k , k = 1, 2, …, N at a given frequency ω and conductivity model σ is performed as follows.First, one calculates the fields B (l) (r k , ω; σ) and E (l) (r, ω; σ) corresponding to lth spatial modes j l (r) describing the inducing (extraneous) source where L stands for a number of spatial modes.Then the elements of Λ k are estimated row-wise using the calculated   ()  ℎ and   ()  ℎ fields.Specifically, elements of Λ k are obtained as the solution of the following two systems of linear equations (SLE) and where l = 1, 2, …, L. Note, that in Equations 4 and 5 the dependency of all quantities on ω,   ()  ℎ and   ()  ℎ on σ, and elements of Λ k on σ, r and r k are omitted but implied.To resolve the TF, the number of spatial modes, L, should be equal to or larger than the doubled number of sites, N, that is, L ≥ 2N.Also, our model experiments indicate that the condition number of the system matrices in Equations 4 and 5 can be quite large, so we used the regularized ordinary least squares approach to solve the above SLEs.
Once elements of Λ k are estimated at a predefined number of frequencies, the GEF at a given time instant t i = iΔt and a given location r is then calculated using a numerical scheme similar to that described by Marshalko et al. (2023) where   obs ℎ stands for horizontal components of the observed field in corresponding locations and N t = T Λ /Δt.As discussed in Marshalko et al. (2023), computation of quantities like Λ (n) (r, r k , T Λ ; σ) can be reduced to the estimation of the following integrals Expressions for Λ (0) and  Λ (  ) are more complicated, but their calculation is similar to that presented in Appendix A of Marshalko et al. (2023) to calculate the corresponding terms for the electric field.
Computation of the integrals in Equation 7 is performed as follows.First, Λ(r, r k , ω; σ) are computed at N f logarithmically spaced frequencies (with 15 frequencies per decade) from f min to f max using Equations 4 and 5. Further, one can analytically compute the corresponding integrals by exploiting cubic spline interpolation as applied to calculated Λ.Our model experiments (not shown in the paper) indicate that N f = 71, f min = 6.13 ⋅ 10 −7 Hz, and f max = 0.054 Hz provide an accurate (with relative error of 0.001%) estimation of the integrals, in the assumption that Δt is taken as 10 s.Note that this approach to estimate similar integrals was successfully implemented in our previous studies (Kruglyakov, Kuvshinov, & Nair, 2022;Marshalko et al., 2023).
As seen from Equations 4 and 5, the estimation of Λ(r, r k , ω; σ) requires specification of j l , l = 1, 2, …, L. The form of j l (and their number, L) varies with the application.As discussed in Introduction, in this paper, we use spatial modes obtained by means of the time-domain PCA of the SECS-recovered inducing source.The reader can find further details on the derivation of these modes in the next section.
As also seen from Equations 4 and 5 (and already mentioned before), the estimation of Λ(r, r k , ω; σ) implies computations of B (l) (r k , ω; σ) and E (l) (r, ω; σ) in a given conductivity model σ of the region of interest.We performed these computations using the scalable 3-D EM forward modeling code PGIEM2G (Kruglyakov & Kuvshinov, 2018) based on a method of volume integral equation.

10.1029/2023SW003621
4 of 15 Worth noting that quantities Λ (n) (r, r k , T Λ ; σ), n = 0, 1, …, N t are time-invariant, and for a given conductivity model are calculated only once, then stored and used when the calculation of E h (r, t; σ) is required.
It is also important to note that the "memory" T Λ can be taken as short as 15 min (at least for Fennoscandia-the region we have selected to validate the presented approach).Namely, a small value for T Λ allows us to calculate E h (r, t; σ) in real time, provided long-term continuous observations of the magnetic field,   obs ℎ are available at multiple locations.Further justification of T Λ smallness is presented in the following two sections.
Ideally, Λ should relate the electric field with the magnetic field measured at as large as possible number of sites.However, there are usually gaps in the magnetic field data.At the same time, to apply the presented methodology, one needs to use continuous time series of magnetic field during the interval [t i − T Λ , t i ] to compute E h (r, t i ; σ).Thus one has to choose observatories accordingly to continuous data availability.However, one can use different (precomputed) sets of Λ for different time instances t i .Combining this option with the short T Λ means that the number of sites with magnetic field observations in use is always large enough to implement the proposed approach.
The final remark of this section is that if the spatial modes are just two polarizations of vertically indenting plane wave, then Λ degenerates to intersite impedance Z; cf.Section 2.2 of Marshalko et al. (2023).

3-D Conductivity Model
One of the critical components required for as realistic as feasible GEF modeling is a trustworthy conductivity model of the Earth's subsurface in the region of interest.The conductivity model of Fennoscandia we adopt comprises a 3-D part (upper 60 km below the surface of the Earth) and the 1-D part underneath taken as a corresponding piece (deeper than 60 km) of the 1-D profile from Kuvshinov et al. (2021).As for the 3-D part, it is based on the SMAP model (Korja et al., 2002), covers the area of 2,550 × 2,550 km 2 and consists of three layers of laterally variable conductivity of 10, 20, and 30 km thicknesses.The lateral discretization of the model is 512 × 512 cells, meaning that cells' size is about 5 × 5 km 2 .Note that this model (of the same discretization) was also used in Marshalko et al. (2021), Kruglyakov, Kuvshinov, and Marshalko (2022), and Marshalko et al. (2023).The interested reader can find details on the conductivity distribution in the model, for example, in Kruglyakov, Kuvshinov, and Marshalko (2022), cf.their Figure 7.Note also that for T Λ calculations, only observatories located in the area covered by this model are used, that is, not all IMAGE observatories.

Deriving Spatial Modes
As mentioned in the Introduction, spatial modes j l (r) are obtained using the SECS method and their derivation involves the following two steps (Kruglyakov, Kuvshinov, & Marshalko, 2022): 1. Spherical Elementary Current Systems (SECS) method (Vanhamäki & Juusola, 2020) is applied to 29-31 October 2003 (Halloween storm) ten-second IMAGE magnetic field data to separate the inducing and induced current systems that are assumed to flow 90 km above the Earth's surface and 1 m below the Earth's surface, correspondingly.The data from all 26 magnetometers were used to perform the SECS analysis.The locations of IMAGE magnetometers in use are demonstrated in Figure 1.Note that IMAGE data for this 72-hr time interval contain several gaps; linear interpolation was applied to fill these gaps.2. The time-domain PCA is applied to the SECS-recovered inducing source resulting in the desired j l (r), l = 1, 2, …, L. With L = 34, we described 99.9% of the inducing source variability.It is relevant to mention that later in the paper, we analyze the data from the Halloween storm and other space weather events.One can argue that j l obtained for a specific event could be non-adequate for other events.Kruglyakov, Kuvshinov, and Marshalko (2022) addressed this question and demonstrated 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 j l obtained from the analysis of some specific event (in our case the Halloween storm).The prerequisite to getting an adequate set of j l is that the event for their estimation should be long enough and sufficiently energetically large and spatially complex.

Comparing Modeled and Observed GEF
Even though there are multiple locations in Fennoscandia where continuous GEF measurements were performed for relatively long time intervals (from weeks to months), our analysis of the available data reveals that most of the GEF observations were conducted during periods of relatively quiet geomagnetic activity.Note that the potentially hazardous GIC are flowing in the technological systems during geomagnetic disturbances like geomagnetic storms and magnetospheric substorms, making such space weather events of particular interest in the context of this study.Fortunately, continuous (1-second) GEF measurements performed in August 2005-June 2006 at 8 of 13 sites around Joensuu, Finland, in the framework of Electromagnetic Mini Array (EMMA) Project (Smirnov et al., 2006) have successfully recorded prominent space weather event-magnetospheric substorm of 11 September 2005; during this event, amplitudes of the observed GEF in the region exceeded 2,000 mV/km.We compare modeled and measured time series of the GEF during this event at two sites (denoted below as M02 and M05) where the largest GEF amplitudes were observed and where the data quality-in terms of absence of gaps, jumps, and electrodes' drift-was the highest.The comparison is performed for the time interval from 05:15 UT to 6:15 UT.This interval includes a 30-min long event and 15-min long quiet periods before and after the event.
During this time interval, 16 IMAGE observation sites delivered uninterrupted magnetic field recordings, and thus, only the data from these sites were used to calculate the GEF using Equation 6. Figure 2 shows the location of these IMAGE sites (blue circles) and two sites (red circles) with GEF recordings selected for comparison of the observed and modeled GEF.Since GEF measurements were performed in geomagnetic coordinates, we decided to compare observed and modeled GEF components in this coordinate system.Thus, in this section, the x-and y-components of the GEF are geomagnetic North and East components, respectively.Besides, experimental 1-s time series of the GEF were re-sampled to a 10-s time series to make them compatible with IMAGE magnetic field data.
It is important to note here that depending on location, the modeled GEF might still over-or underestimate the amplitudes of the actual GEF.This is due to galvanic effects, that is, the build-up of electric charges along local near-surface heterogeneities (Jiracek, 1990) that cannot be included in the model.Galvanic effects are commonly accounted for with a 2 × 2 real-valued time-independent-unique for each location-distortion matrix, D, which linearly relates observed   obs ℎ and modeled   mod ℎ electric fields as (Püthe et al., 2014) Having observed and calculated GEF at many time instants and bearing in mind that D is time-independent, we can form highly over-determined SLE given by Equation 8to estimate elements of D using the least-square method.Specifically, these elements are estimated row-wise as the solutions of the following two SLEs and where t i runs multiple (K) time instants.Note that the dependence of all quantities in Equations 9 and 10 on r is omitted but implied.To make the comparison as fair as possible, we calculated distortion matrices using the data not from the time interval where we compare the results (05:15-06:15 September 11) but using the data far before the event, specifically, the data from one (whole; September 9) day, giving K = 24 × 3,600/10 = 8,640 equations to estimate respective two elements.
As stated in the previous section, the memory T Λ can be taken as short as 15 min in Equation 6.To explore whether this is the case, we performed GEF calculations taking T Λ as 15, 30, 60, and 90 min.Figures 3 and 4 compare observed and (corrected for the distortion) modeled E x and E y , respectively, at sites M02 and M05 for different values of T Λ .Note that the corrected modeled GEF are calculated by reusing Equation 8, that is, as where elements of D are obtained from the solution of the SLEs ( 9) and (10).As is seen, the modeled results for all values of T Λ are indistinguishably similar, meaning that T Λ can be taken as 15 min.Tables 1 and 2 confirm this inference more quantitatively by showing correlation coefficients and coefficients of determination at sites M02 and M05, respectively, for the considered 1-hr time interval.One can see that both quantities practically do not differ with respect to the value of T Λ .
What is also noticeable from the tables is that the agreement between experimental and modeled GEF (in terms of both correlation coefficients and coefficients of determination) is better in the E y component.This is most probably because the ionospheric current (source) flows predominantly in the y direction, making GEF modeling in this direction more accurate.Returning to Figures 3 and 4, we also observe that the experimental GEFs at selected sites separated by only 21 km significantly differ, more considerably than the modeled GEFs.The most probable reason for this is the existence of local 3-D conductivity heterogeneities in this region which are not accounted for in the available 3-D conductivity model of the region.

Comparing Observed and Modeled GIC
The conventional approach to simulate GIC in the pipelines is based on the following linear model (Boteler, 2013;Boteler & Pirjola, 2014) Figure 2. Location of IMAGE sites (blue-filled circles) used for ground electric field (GEF) simulation (see Equation 6).
Red circles-the location of sites where the GEF was measured and used for comparison with the simulated GEF. 10.1029/2023SW003621 7 of 15 where M is the number of pipeline nodes and ends, r i -coordinates of their location, and r is the coordinate of the pipeline point where GIC is recorded.Note that coefficients a i and b i are time-independent, but depend on pipeline physical parameters.In case when the conductivity distribution and the source are both laterally uniform, Equation 12 degenerates to
Having observed GIC and calculated GEF at many time instants, and bearing in mind that coefficients a i and b i are time-independent, we can form-as for elements of matrix D discussed in the previous section-highly over-determined SLE given by Equation 12to estimate coefficients using the least-squares method.Precisely, these coefficients are calculated as the solution of the following SLE where t i runs multiple (K) time instants.Here we omit the dependence of the coefficients and GIC obs on r because we have GIC's recordings only at one (Mäntsälä; MAN) point.Figure 5 shows the general geometry of the pipeline (gray lines) and location of 18 pipeline's nodes and ends (gold diamonds).The geometry of the pipeline is based on models presented in Pulkkinen et al. (2001) and Dimmock et al. (2019).
GIC measurements at the Mäntsälä point have been performed continuously since 1999; thus, we had an opportunity to compare observed and modeled GIC for a number of events.We have chosen events that happened in 2000, 2001, and 2003 years.It should be noted that the approach we use to determine coefficients a i and b i is implicitly based on the assumption that the configuration of the pipeline system was the same during the events in study.
Concerning the relatively short period of 2000-2003, no major modifications of the pipeline network were made.However, the present pipeline system is more extensive than in the beginning of 2000s.Note also that from 2005 the quality of GIC measurements degraded; this is why we considered earlier years' events.Also, note that the modeled GIC are calculated by reusing Equation 12 mod () = where a i and b i are obtained from the solution of the SLE ( 14).GEF in Equation 15 were calculated using T Λ = 15 min; our model experiments (not shown in the paper) with larger T Λ reveal a negligible difference in the modeled GIC.As in the case of the GEF comparison discussed in the previous section, we estimated a i and b i using GIC data not from the time interval where we compare the results but using continuous data for 27-28 October 2003, giving K = 2 × 24 × 3,600/10 = 17,280 equations to estimate respective 2 × M = 36 coefficients.
Top two panels in Figures 6-8 demonstrate time series of observed and modeled GIC for two 3-hr time intervals of corresponding years/months when large fluctuations of GIC were detected.As one can see, the agreement between observed and modeled GIC is remarkably good.The bottom right panels in the figures demonstrate the agreement differently.It shows cross-plots of the observed and modeled GIC.Ideally, the observed and modeled GIC would lie in a straight line, and this is indeed the case.Note that the bottom left panel shows IMAGE observatories, the data from which were used to calculate the GEF using Equation 6.
Finally, Table 3 demonstrates the agreement in a more quantitative way by presenting correlation coefficients and coefficients of determination for the considered events.One can see that both coefficients are close to 1, irrespective of the event.Interestingly, the agreement in GIC is better than in GEF (see previous section).We attribute better agreement in GIC to a smoothing effect of summation (see Equation 15) which effectively suppresses the potential inaccuracy of the modeled GEF at considered nodes and ends.

Conclusions
This paper presents a novel approach to the GEF calculation in real time.The approach makes it possible to take into account the 3-D Earth's conductivity distribution in the region and lateral variability of the source.It works with spatial modes describing the spatial structure of the source and exploits the multi-site (MS) transfer function (TF) concept where MS TF relates the GEF at any location with horizontal magnetic field at (fixed) multiple locations.Using Fennoscandia as an example region and the SECS method to derive the spatial modes, we compared the observed time series of the GEF and GIC available in the region with those simulated using the proposed approach.Good agreement between observed and modeled results validates the methodology.Notably, in contrast to the previous study (Marshalko et al., 2023), where GIC was calculated using the GEF at the pipeline node where GIC is measured, in this paper, we considered a more realistic scenario when GIC is calculated using simulated GEF at multiple nodes of the actual pipeline grid.
While the simulations of GICs in this paper are mainly done for validation purposes, the presented approach opens an avenue to derive a GIC activity indicator.Combining ( 6) and ( 15) one can estimate GIC in real time.Moreover, as demonstrated in this paper, the estimated GIC values correspond very accurately to the true recordings for a particular period when the conductor network can be assumed to be unchanged.Outside of this period, the resulting GIC provides a meaningful proxy.Although it may not any more give the true current, it can still be compared to other events in a relative sense.For example, in our case, we could use the Halloween storm as a benchmark.Thus this method has an immediate potential to be implemented as an operative product for space weather services.
The prerequisite for the method's application is continuous magnetic field measurements at multiple locations, as in Fennoscandia.North America and China (with existing networks of magnetic field observations), or/and New Zealand (with the network of variometers is being established now) are regions where the proposed method is worth trying.
It is worth noting that while comparing measured and simulated GEF at two selected sites separated by only 21 km, we observe that experimental GEFs at these sites differ more significantly than the modeled GEFs.Bearing in mind that the inducing source (equivalent ionospheric current) is relatively (spatially) smooth, the most probable reason for substantial lateral variability of the experimental GEF is local sharp lateral gradients of the subsurface conductivity not accounted for in the used 3-D conductivity model of the region.This prompts building a new-more accurate and detailed-3-D conductivity model of Fennoscandia via a multi-scale 3-D inversion (using new modern inverse solvers) of both old data (on which SMAP model was based) and a large amount of new MT data collected in the region in the framework of various MT projects.
Finally, we have to mention that the proposed concept of multi-site transfer functions can be easily adapted to another space weather related problem-accounting for the geomagnetic disturbances during offshore directional drilling in northern seas (like offshore of Alaska, Norway, British Isles, and North Russia), where magnetic field variations are routinely recorded at multiple adjacent land-based locations.In this scenario, MS TF will relate magnetic field variations at sea bottom with those at multiple land-based sites.

Table 3
The Agreement Between Modeled and Observed GIC in Terms of Correlation Coefficients and Coefficients of Determination for Different Events (Substorms) 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).We acknowledge Gasum Oy for long-term collaboration in GIC studies of the Finnish natural gas pipeline.Open access publishing facilitated by University of Otago, as part of the Wiley -University of Otago agreement via the Council of Australian University Librarians.

Figure 1 .
Figure 1.Location of IMAGE observatories with available data for 29-31 October 2003.

Figure 3 .
Figure 3.Comparison of the (geomagnetic) northward component of observed and modeled ground electric field at M02 and M05 sites for different values of T Λ .

Figure 4 .
Figure 4. Comparison of the (geomagnetic) eastward component of observed and modeled ground electric field at M02 and M05 sites for different values of T Λ .

Figure 5 .
Figure 5. Geometry of the pipeline (gray lines) and location of 18 pipeline's nodes and ends (gold diamonds); Mäntsälä (MAN) natural gas pipeline recording point (open red circle).For reference, the location of the IMAGE observatory, Nurmijärvi (NUR; blue-filled circle) is shown.

Figure 6 .
Figure 6.Substorms of 15/16 July 2000.Top two panels: comparison of observed and modeled GIC for two 3-hr time intervals.Bottom left panel: location of 13 IMAGE sites (blue-filled circles) used for the ground electric field simulation.Bottom right panel: cross-plots of observed and modeled GIC (cross-plots are demonstrated for the 2-day time interval: July 15-16 2000).

Figure 7 .
Figure 7. Substorms of 6 November 2001.Top two panels: comparison of observed and modeled GIC for two 3-hr time intervals.Bottom left panel: location of 18 IMAGE sites (blue-filled circles) used for the ground electric field simulation.Bottom right panel: cross-plots of observed and modeled GIC (cross-plots are demonstrated for the 2-day time interval: November 5-6 2001).

Figure 8 .
Figure 8. Substorms of 29/30 October 2003.Top two panels: comparison of observed and modeled GIC for two 3-hr time intervals.Bottom left panel: location of 17 IMAGE sites (blue-filled circles) used for the ground electric field simulation.Bottom right panel: cross-plots of observed and modeled GIC (cross-plots are demonstrated for the 3-day time interval: October 29-1 November 2003).

Table 2
The Agreement Between Modeled and Observed Ground Electric Field at Site M05 in Terms of Correlation Coefficients and Coefficients of Determination