Three‐Dimensional Modeling of the Ground Electric Field in Fennoscandia During the Halloween Geomagnetic Storm

In this study, we perform three‐dimensional (3‐D) ground electric field (GEF) modeling in Fennoscandia for three days of the Halloween geomagnetic storm (29–31 October 2003) using magnetic field data from the International Monitor for Auroral Geomagnetic Effects (IMAGE) magnetometer network and a 3‐D conductivity model of the region. To explore the influence of the inducing source model on 3‐D GEF simulations, we consider three different approaches to source approximation. Within the first two approaches, the source varies laterally, whereas in the third method, the GEF is calculated by implementing the time‐domain realization of the magnetotelluric intersite impedance method. We then compare GEF‐based geomagnetically induced current (GIC) with observations at the Mäntsälä natural gas pipeline recording point. We conclude that a high correlation between modeled and recorded GIC is observed for all considered approaches. The highest correlation is achieved when performing a 3‐D GEF simulation using a “conductivity‐based” laterally nonuniform inducing source. Our results also highlight the strong dependence of the GEF on the earth's conductivity distribution.

The primary aim of the current study is to perform three-dimensional (3-D) GEF modeling in Fennoscandia during the Halloween geomagnetic storm using available observed geomagnetic field data (Tanskanen, 2009) and a 3-D conductivity model of the region (Korja et al., 2002).As it was mentioned by Pulkkinen et al. (2017), from the engineering point of view, the spatiotemporal characteristics of the horizontal GEF provide the ideal description of a geomagnetic disturbance.Simulated GEF data of a very intense event can serve as a reference point for evaluating possible risks to ground-based technological systems in Fennoscandia from space weather, as GIC can be calculated based on the GEF data in the region, the geometry of a technological network and system design parameters (Lehtinen & Pirjola, 1985;Pirjola et al., 2022).
Worth attention, the GEF and GIC in the Fennoscandian region, and also elsewhere in Europe, have been previously modeled in several projects as described by, for example, Viljanen et al. (2014) and Myllys et al. (2014).Wei et al. (2013) followed a similar idea to model the GEF in North America.A key difference to the present work is that usually, to calculate the GEF, researchers employ the plane wave method; moreover, in most studies, only 1-D ground conductivity models have been utilized.Although different 1-D models were used for different locations, the 1-D approach cannot take properly into account effects in the GEF arising from lateral gradients in 3-D conductivity distributions (Ivannikova et al., 2018;Kelbert, 2020;Marshalko et al., 2020;Rosenqvist & Hall, 2019).
To explore the influence of the inducing source model on 3-D GEF simulations, we consider three different approaches to source approximation.Noteworthy, all methods exploit the same International Monitor for Auroral Geomagnetic Effects (IMAGE) magnetic field data to simulate the GEF, rely on the same 3-D conductivity model of the region, and use the same forward problem engine.Within the first two approaches, the source varies laterally and is factorized by spatial modes (SM) and respective expansion coefficients.In both approaches, the SM are the same and are obtained following the two-step numerical scheme introduced by Kruglyakov, Kuvshinov, and Marshalko (2022).The difference between methods lies in the calculation of the expansion coefficients.The details on these two approaches (and the reasoning to invoke the second approach) are presented in Section 2.1.Within the third approach (discussed in the same section), the GEF is calculated by implementing the time-domain realization of the magnetotelluric (MT) intersite impedance method.
Further, in Section 3, we present the results of the source recovery using the first two approaches and the results of the GEF modeling obtained using three considered methods.Besides, we explore in this section how well the observed time series of GIC at Mäntsälä natural gas pipeline recording point (Viljanen et al., 2006) during the Halloween storm are reproduced through a linear combination of the simulated horizontal GEF components at this point.We note once again that the Halloween event is chosen because it is a representative example of a big geomagnetic storm causing the largest value of GIC in the Finnish natural gas pipeline (Dimmock et al., 2019).
A summary of the results and a discussion of the possible ways forward are presented in Section 4.

Methodology
In this section, we present and discuss three approaches which we invoke to calculate the time-domain GEF.The first two methods rely on laterally varying inducing source models and the third one-on plane-wave excitation.EM modeling is performed for three days (72 hr) of the Halloween geomagnetic storm (29-31 October 2003).

Governing Equations in the Frequency Domain
We start with a discussion of the problem in the frequency domain.Maxwell's equations govern EM field variations, and, in the frequency domain, these equations read as (1) 10.1029/2022SW003370 3 of 19 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 = (x, y, z) is a position vector, in our case in the Cartesian geometry.Note that we neglected displacement currents and adopted the following Fourier convention: In problem setups, when a laterally nonuniform source is considered, we assume that the current density, j ext (r, ω), can be represented as a linear combination of SM j i (r): The form of SM 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 (2013b), Honkonen et al. (2018), Guzavina et al. (2019), Grayver et al. (2021), andKruglyakov, Kuvshinov, andNair (2022), current loops in Sun and Egbert (2012), eigenmodes from the Principal Component Analysis (PCA) of the physics-based models in Egbert et al. (2021) and Zenhausern et al. (2021), and eigenmodes from the PCA of the data-based models in Kruglyakov, Kuvshinov, and Marshalko (2022).In this paper, we will use the parameterization adopted in the latter paper.
By virtue of the linearity of Maxwell's equations with respect to the j ext (r, ω) term, we can expand electric and magnetic fields as linear combinations of individual fields E i and B i , where E i (r, ω; σ) and B i (r, ω; σ) fields are "electric" and "magnetic" solutions of the following Maxwell's equations: ∇ ×  = i. (8)

Governing Equations in the Time Domain
The transformation of Equations 5 and 6 into the time domain leads to the representation of the electric and magnetic fields as The reader is referred to Appendix A in Kruglyakov, Kuvshinov, and Marshalko (2022) for more details on the convolution integrals in the above equations.We note that we use the same notation for the fields in the time and frequency domains.Equations 9 and 10 show how the fields can be calculated provided c i (t) and conductivity model σ(r) are given.To make formulas ready for implementation, one needs to estimate the upper limits of integrals in the above equations, or, in other words, to evaluate time intervals, T E and T B , above which E i (r, τ; σ) and B i (r, τ; σ) become negligibly small.The latter will allow us to approximate Equations 9 and 10 as The details of the numerical calculation of the integrals in the above equations are presented in Kruglyakov, Kuvshinov, and Marshalko (2022) and Kruglyakov, Kuvshinov, and Nair (2022)  are more complicated and are presented in Appendix A of the current paper.A few comments on the latter equations are relevant at this point.
• T E and T B significantly differ.As shown by Kruglyakov, Kuvshinov, and Marshalko (2022) T E can be taken as small as 15 min for Fennoscandia.As for T B , it is several orders of magnitude larger than T E (Kruglyakov, Kuvshinov, & Nair, 2022); specifically T B should be taken as large as half a year.Note that in this study we only model the GEF in the region, not the magnetic field.Therefore, only T E is relevant for us.• Computation of the integrals in the right-hand side of Equations 15 and 16 is performed as follows.First, E i (r, ω; σ) and B i (r, ω; σ) are computed at 60 logarithmically spaced frequencies from 3.67 × 10 −6 to 0.054 Hz.Note that for the magnetic field, modeling at zero frequency is also required (see Equation A4 from Appendix A).Further, using cubic spline interpolation as applied to calculated E i (r, ω; σ) and B i (r, ω; σ), one can analytically compute the corresponding integrals.• Quantities      (, ; ) and      (, ; ) are time-invariant, and-for the predefined set of j i and a given conductivity model-are calculated only once, then stored and used when the calculation of E(r, t; σ) and B(r, t; σ) is required.
• One of the key ingredients to make regional EM modeling in Fennoscandia as realistic as feasible is a conductivity model of the Earth's subsurface of the region.The model adopted in this paper comprises a 3-D part (upper 60 km) and an underlying 1-D conductivity profile (Kuvshinov et al., 2021) (a part of the profile deeper than 60 km below the surface of the Earth).3-D part 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; conductivity distributions in these layers are shown in Figure 4.The lateral discretization of the model is 512 × 512 cells.Note that this model was also exploited in Marshalko et al. (2021) and Kruglyakov, Kuvshinov, and Marshalko (2022).Another note is that computations of electric and magnetic fields in this model (at a given frequency) are performed using the scalable 3-D EM forward modeling code PGIEM2G (Kruglyakov & Kuvshinov, 2018) based on a method of volume integral equation with a contracting kernel (Pankratov & Kuvshinov, 2016).
• As seen from the above equations, GEF computations require specification of SM j i (r) and estimation of time series of coefficients c i (t).We address this topic in the next two sections.

GEF Modeling Using the SECS-Based Approach
As mentioned in the Introduction, SM j i (r) (and corresponding time series of expansion coefficients c i (t)) can be obtained using the following two-step scheme (Kruglyakov, Kuvshinov, & Marshalko, 2022): 1. Spherical Elementary Current Systems (SECS) method (Vanhamäki & Juusola, 2020) is applied to 29-31 October 2003 IMAGE magnetic field data to separate the inducing and induced current systems that 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 SECS analysis.The location of IMAGE magnetometers is demonstrated in Figure 1.Note that IMAGE data for this time interval (72 hr) contain several gaps; linear interpolation was used to obtain magnetic field data in the gaps.2. The PCA is applied to the SECS-recovered inducing source resulting in the desired SM j i (r), i = 1, 2, …, L, and time series of the corresponding expansion coefficients c i (t).With L = 34, we succeeded in describing 99.9% of the inducing source variability.
Once j i (r) and c i (t) are derived, the GEF can be computed using Equation 13.Hereinafter we will refer to this two-step approach to specify the source and, consequently, compute the GEF as the SECS-based method.
It is important to discuss here one potential drawback of the SECS-based approach to isolate the inducing source (see also Juusola et al. (2020), Section 4.3).Note that this method can be viewed as a regional variant of the Gauss method-the method widely used to separate the inducing (external) and induced (internal) sources on a global scale.If the region of interest is characterized substantially by 3-D conductivity distribution (as in our case) the induced part is inevitably influenced by 3-D effects arising, in particular, from the lateral (e.g., land/ocean) conductivity contrasts.Given the deficient spatial distribution of the IMAGE sites, the SECS-based approach precludes an accurate description of the induced part affected by localized 3-D effects.Evidently, such an imperfection in the induced part description also influences the recovery of the inducing part, at least in terms of c i (t) recovery.
In the next section, we discuss an approach to circumvent this issue.

GEF Modeling Using the Conductivity-Based Inducing Source
Let us first assume that the IMAGE data analysis discussed in the previous section gives us trustworthy SM j i (r).
Assume further that the ground 3-D conductivity distribution is known to us at the inducing source construction stage; this is the reason why we call the inducing source discussed in this section the conductivity-based (CB) source.
With the above assumptions in mind, the most adequate way to obtain c i (t) at a given time instant is to reuse Equation 14.Specifically, the calculation of c i (t) at a given time instant t k = kΔt is performed as follows.Substituting coordinates of IMAGE sites into Equation 14and rearranging the terms, we obtain a system of equations to determine where , j = 1, 2, …, J, and J is the number of IMAGE sites.The expression ( 17) represents an overdetermined system of linear equations which is solved by the least-square method.This scheme was implemented and validated by Kruglyakov, Kuvshinov, and Nair (2022), who analyzed hourly-mean mid-latitude magnetic field signals of magnetospheric/ionospheric origin.As it is seen from Equation 17, computational loads to obtain c i (t k ), that is, for single time instance, are proportional to  ×   .In our scenario, L = 34,    = 180 × 24 × 60 × 6 making computational loads to be prohibitively high; recall that the value for N T = T B /Δt is obtained assuming that T B is taken as half of the year, that is, T B = 180 × 24 × 60 × 60 s, and Δt = 10 s.Note that the above-discussed approaches were developed to perform the near real-time calculations of the GEF (and magnetic field if needed).In this paper, we are interested in computing the GEF for the specific event, and, thus, we can exploit an alternative variant of the GEF modeling approach with the use of the CB source.Noteworthy, this variant has been routinely used for the last two decades to analyze the ground-based signals of magnetospheric origin (Honkonen et al., 2018;Munch et al., 2020;Olsen & Kuvshinov, 2004;Püthe & Kuvshinov, 2013a;Püthe, Kuvshinov, & Olsen, 2014).As applied to our problem setup, this variant of the method includes the following steps: 1. Magnetic field data B obs (r j , t), j = 1, 2, …, J recorded at IMAGE sites are converted from the time to frequency domain using the fast Fourier Transform (FFT).Note that in our case, J = 23; sites Ny-Ålesund, Longyearbyen, and Hornsund fall outside the modeling region.2. At each FFT frequency ω, we estimate c i (ω) by solving the over-determined system of linear equations by means of the regularized least squares method.Note that FFT frequencies range between where S is the length of the event (72 hr). 3. Time series c i (t), i = 1, …, L are then obtained by means of the inverse FFT of frequency domain coefficients c i (ω). 4. Finally, the GEF at a given time instant t k and location r is computed using Equation 13.

GEF Modeling Using the MT Intersite Impedance Method
Although, in reality, the source of the ground EM field is always laterally variable, the conventional approach to model GEF relies on the plane-wave source assumption.This assumption allows researchers to relate the frequency domain (horizontal) GEF at point r with the horizontal magnetic field at a base site r b through an intersite impedance (Kruglyakov & Kuvshinov, 2019) Note that if r coincides with r b , an intersite impedance transforms into a standard MT impedance (Berdichevsky & Dmitriev, 2008).
In order to calculate the GEF via the MT intersite impedance method, we first carry out 3-D EM forward modeling using PGIEM2G code (Kruglyakov & Kuvshinov, 2018) and the SMAP conductivity model (Korja et al., 2002) with two (laterally uniform) plane wave sources at the same 60 frequencies that are indicated in Section 2.1.2.3-D intersite impedances (Equation 19) are then computed for each of these frequencies.
The GEF at a given time instant t k and location r is calculated using a numerical scheme similar to that described in Section 2.1.2(cf., Equation 13), namely ℎ(, ; ) ≈ 1 In this study, we use the data from the nearest IMAGE magnetometer to calculate the GEF at a particular point using the MT intersite impedance method.

Original Current Versus PCA-Constrained Current
Since two approaches discussed in Sections 2.1.3and 2.1.4exploit PCA-recovered SM, we explore in this section how well the (ionospheric) equivalent current calculated using L = 34 SM identified by the PCA fits the original equivalent current obtained using the SECS technique.Figure 2 demonstrates time series of the aforementioned currents above two exemplary sites: Abisko (ABK; latitude: 68.35°N, longitude: 18.82°E) and Nurmijärvi (NUR; latitude: 60.5°N, longitude: 24.65°E); their locations are shown in Figure 3.One can hardly see the difference between the results.A perfect fit is also quantified in terms of high correlation between time series (0.9997), low normalized root-mean-square-errors (nRMSE; lower than 0.023), and low maximum absolute differences (MAD; lower than 0.041 A/m) for both x (north) and y (east) components of the equivalent current at each location.Note that nRMSE is defined as follows: where a is the ionospheric equivalent current calculated using PCA-recovered SM and b is the original ionospheric equivalent current, a i and b i are elements of these time series, and N l is the number of time instants (in our case, 6 × 60 × 72).
We conclude from this comparison that spatial structure of the equivalent current (at least for the considered 72-hr event) is very well explained by L = 34 PCA-based SM, thus supporting the usage of this SM basis in the GEF modeling approaches exploiting SECS and CB equivalent currents.

SECS-Based Current Versus Conductivity-Based Current
In this section, we compare SECS and CB equivalent currents.Figure 3a snapshots of these currents as well as current's time derivatives above Fennoscandia at 20:08:30 UT, 30 October 2003, -the moment of the largest amplification of GIC in the natural gas pipeline during the substorm event, which caused the blackout in Malmö, Sweden.It can be seen that the overall behavior of equivalent currents is similar.This is also true for the time derivatives of the equivalent currents.However, the SECS-based source and its time derivatives have a smoother spatial structure compared to the CB source.Figure 3b demonstrates the time series of the SECS and CB equivalent currents above and NUR geomagnetic observatories.It can be seen that the of variations larger in the case of the CB source for NUR.The difference between time series is especially prominent in the case of a smaller x-component of the equivalent current.The nRMSE is quite high in this case-0.8944.For the y-component, the nRMSE is 0.2826.For ABK, a good match between time series of equivalent currents is observed.
The detected difference in the recovered equivalent currents will be further assessed in the following sections by comparing GEF and GIC modeled with the use of different approaches.

Comparison of GEF Modeled by Three Methods
Figures 4a and 4b show snapshots of the magnitude and direction of the GEF in Fennoscandia modeled with the use of the SECS and CB sources, respectively, at 20:08:30 UT, 30 October 2003.Figure 4c presents the absolute difference between magnitudes of the GEF demonstrated in Figures 4a and 4b. Figure 4d shows the GEF in the area of the Finnish natural gas pipeline modeled using the CB equivalent current at 20:08:30 UT, 30 October 2003.Finally, Figures 4e-4g demonstrate conductivity distribution in three layers of the 3-D model that we use in our simulations.It is clear that the overall behavior of the GEF obtained with the use of the considered sources is very similar.However, it can be seen that differences in magnitudes of the GEF reach over 8,300 mV/km at this particular time instant at the Norwegian (see Figure 4c).Besides, the behavior of the GEF in the region is complex; the magnitude of the GEF can decrease or increase many times over short distances due to lateral variations of conductivity, especially at the ocean as it is demonstrated in Figures 4a and 4b.
Figure 5 shows the snapshots of the GEF across Fennoscandia (left figures) and in the area of the Finnish natural gas pipeline (right figures) at two moments of maximum GIC amplification at Mäntsälä, which occurred during the initial phase of the Halloween geomagnetic storm: 06:57:30 and 07:27:00 UT, 29 October 2003.
We also compare GEF modeled with the use of three approaches at several locations: ABK, NUR, MAN pipeline recording point (latitude: 60.6°N, longitude: 25.2°E) and Point X, which is located 0.5° north of MAN.Note that in the case of the intersite MT impedance method, ABK magnetic field data were used to calculate the GEF at ABK and NUR magnetic field data were used to calculate the GEF at NUR, MAN, and Point X.
Figures 6-9 demonstrate time series two time intervals: 06:00:00-08:00:00 UT, 29 October 2003 (the initial phase of the Halloween geomagnetic storm, when the largest GIC value was observed at MAN), and 19:00:00-21:00:00 UT, 30 October 2003 (the substorm event, which caused the blackout in Malmö, Sweden).It is worth mentioning that even though NUR, MAN, and Point X are located very close to each other (about 32 km apart in the case of NUR and MAN and about 56 km in the case of MAN and Point X), the magnitude of the GEF variations at these three sites is very different.The highest GEF values are observed at NUR, which is located above a resistive structure and close to a border of conductivity contrast (see Figure 4e).Point X is located above a more conductive structure.That is why GEF values are smaller at this site.
We further quantify the difference between modeled GEF in terms of correlation coefficients, nRMSE, and MAD in Table 1.Note that MAD between x-components of the GEF induced by the SECS and CB source as well as those between x-components of the GEF obtained using CB source and MT impedance method, are quite large at NUR and MAN due to the fact that GEF values at this locations are also large at the moment of maximum GEF amplification.MAD are significantly smaller in the case of Point X.However, the nRMSE between x-components of the GEF induced by SECS and CB source at NUR, MAN, and Point X are practically equal.

Comparison of Modeled and Observed GIC
GIC observations have been carried out at the MAN natural gas pipeline recording point by the Finnish Meteorological Institute (FMI) in collaboration with Gasum Oy starting from November 1998 using the differential magnetometer method.The measurements are made with the help of two magnetometers, one right above the pipe, and another at the NUR geophysical observatory (Pulkkinen et al., 2001).
It was previously demonstrated by Viljanen et al. (2006) that GIC at MAN can be reproduced accurately enough based on the horizontal electric field using the following expression: Figure 10a demonstrates the comparison of the observed GIC and GIC calculated based on the GEF simulated at MAN with the use of the CB source.
Note that modeled GIC are scaled by a factor of 4.51 in the figure.The Huber-weighted robust regression method (Holland & Welsch, 1977) was used to estimate this factor.Note that in the study of Dimmock et al. (2019) who carried out 3-D GEF and GIC modeling in the Fennoscandian region due to 7-8 September 2017 geomagnetic storm, GIC calculated based on Pulkkinen et al. (2001) parameters and GEF simulated in the SMAP conductivity model was also scaled by a factor of 4. Dimmock et al. (2019) point out that conductivities in the model adopted by Viljanen et al. (2006)  Note.In the case of MAN and Point X, MT results are obtained with the use of magnetic field data observed at NUR.   was demonstrated in Section 3.3, GEF values at Point X significantly smaller than those at due to the fact that Point X is located on a more conductive basement.Figure 10d demonstrates comparison of observed GIC and GIC calculated based on the GEF simulated at Point X with the use of the CB source.Modeled GIC are not scaled in this figure.It is clear that simulated GIC variations are of the same order of magnitude as observed ones.
Table 2 presents correlation coefficients, nRMSE, MAD between observed GIC and GIC calculated on the GEF modeled with the use of the three discussed methods.The values are demonstrated for GIC calculated both at MAN and Point X locations.It can be that the highest correlation between modeled and observed GIC with use of the CB source.GIC modeling with the use of this source outperforms GIC calculation via the MT intersite method because, first of all, the GEF modeling approach utilizing the CB source benefits from the spatial nonuniformity of the inducing source, whereas in the MT impedance method, the source is approximated by a plane wave (the inducing source is clearly nonuniform in high latitudes).Second, the GEF modeling approach with the use of the CB source accounts for magnetic field conditions better because it utilizes the data from multiple magnetometers, whereas the MT intersite impedance method relies on the data from a single site.

Concluding Remarks
In this paper, we perform 3-D GEF modeling in Fennoscandia for three days of the Halloween geomagnetic storm (29-31 October 2003).To explore the influence of the inducing source model on 3-D GEF simulations, we consider three different approaches to source approximation.Noteworthy, all methods exploit the same IMAGE magnetic field data to simulate the GEF, rely on the same high-resolution 3-D conductivity model of the region (SMAP, Korja et al. (2002)), and use the same forward problem engine (Kruglyakov & Kuvshinov, 2018).first two approaches, the source and is factorized from the original SECS-recovered source Juusola, 2020) by SM and expansion It should be noted that the insufficient density of the IMAGE magnetometer network affects the accuracy the inducing source recovery via the SECS method as it was mentioned in Section 2.1.3,which, in turn, influences the results of the sponding GEF and GIC simulations.a-c), while there is no scaling in Figure (d).GIC in all figures were calculated based on the ground electric field modeled using the CB source.
In both GEF modeling approaches with the use of a laterally nonuniform inducing source, the SM are the same and are obtained following the two-step numerical scheme introduced by Kruglyakov, Kuvshinov, and Marshalko (2022).The difference between methods lies in the calculation of the expansion coefficients; in the second approach, the expansion coefficients are obtained by taking the conductivity distribution of the earth into account.Within the third approach, the GEF calculated by implementing the time-domain realization of the MT intersite impedance method.
We modeled GIC at the MAN Finnish natural gas pipeline recording site based on the GEF obtained with the use three aforementioned modeling approaches and compared results with GIC observed there.We conclude that for all considered the correlation between and observed GIC is high.The highest correlation with GIC recordings and the lowest nRMSE is achieved with the use of the CB However, when calculating GIC based on the GEF simulated at MAN, their values appear to be overestimated 4-6 times for all modeling techniques.Similar results obtained by et al. ( 2019), who GIC MAN based on GEF computed in the SMAP model during the 7-8 September 2017 geomagnetic storm.When calculating GIC on the basis of the GEF at a point located 0.5° north of MAN (and on a significantly conductive basement), resulting GIC have the same order magnitude observed ones.As this example demonstrates, ground conductivity crucial role in estimating the GEF.Especially challenging are regions sharp gradients of near-surface conductivity.We also stress that in contrast to the spatially highly variable ground conductivity, the equivalent ionospheric currents and their time derivatives are relatively smooth.Thus, most of the lateral variation of the GEF arises from the ground conductivity.
For GIC calculation in this study, we used a simplified method, which assumes that the GEF along the pipeline is spatially uniform.One can argue that this approximation is too rough, taking into account that our modeling results demonstrate large differences between GEF values at different sites in the pipeline area (see Figures 4d  and 7-9).Calculation of the actual GIC in technological networks is an engineering task, which can be performed with a significantly higher level of accuracy by companies operating these networks and possessing all the necessary information about their configurations and parameters.Moreover, information about changes in configurations of technological systems over the years is required to model GIC properly during a particular time interval.That is why in this study, we limit ourselves to using this simplified GIC modeling method, with the help of which reasonable modeling results were previously obtained by Viljanen et al. (2006).We also share our GEF simulation results through an open-access repository (Marshalko et al., 2022).With the help of these data, companies operating technological systems in Fennoscandia will be able to assess the potential hazard to these systems from space weather.
Concerning future studies, our results of the Halloween storm serve as an explicit point of comparison.Using the same ground conductivity model, but different geomagnetic field input, we can quantify the magnitude of other events with respect to the Halloween storm.Of special interest is a recent reproduction of the Carrington storm by Blake et al. (2021) and the simulation of a sudden storm commencement due to an extreme solar wind shock (Welling et al., 2021).
Another topic for future activity is updating the 3-D conductivity model of Fennoscandia, which is to be based on a multi-scale 3-D inversion of a significant amount of new MT data collected in the region in the framework of various goal-oriented MT projects.Note that the biggest uncertainty in the GEF and GIC modeling arises due to the incompleteness of the conductivity model.In a recent study, Love et al. (2022) considered the famous magnetic storm in March 1989.They presented maps of reported GIC impacts in the contiguous United States (CONUS) power grids and compared their occurrence to the peak values of the GEF based on empirical MT impedances.There is a clear correspondence between the locations of GIC impacts and high GEF values.As Love et al. (2022) point out, geomagnetic variations tend to decrease with decreasing geomagnetic latitude.This is also seen in geoelectric hazard maps.However, the hazard across CONUS is much more prominently organized by the surface impedance, that is, the ground conductivity.This emphasizes the need for using as accurate information on the earth's conductivity as possible.In this paper, we do not compare the modeled GEF to the observed one, since, to the best of our knowledge, no GEF measurements were made in Fennoscandia during the Halloween geomagnetic storm in 2003.It should be noted that the actual GEF is subject to strong galvanic distortion (Jones, 2012) due to very local conductivity inhomogeneities, which usually cannot be accounted for during 3-D EM modeling.However, galvanic effects can be taken into account by estimating the so-called distortion matrix in a manner, which was proposed by Püthe, Manoj, and Kuvshinov (2014).In the future, we plan to analyze the available data from various MT surveys and compare observed and modeled GEF with the use of the distortion matrix correction for more recent space weather events.
Finally, the ongoing research aims to further develop the GEF modeling approach with the use of the conductivity-based inducing source to enable its real-time implementation.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).We acknowledge Gasum Oy for long-term collaboration in GIC studies of the Finnish natural gas pipeline.

Figure 1 .
Figure 1.Location of sites of the International Monitor for Auroral Geomagnetic Effects (IMAGE) magnetometer network during 29-31 October 2003.

Figure 2 .
Figure2.Time series of the original ionospheric equivalent current obtained using the Spherical Elementary Current Systems method and ionospheric equivalent current calculated using 34 spatial modes above two exemplary sites: Abisko (ABK) and Nurmijärvi (NUR).The results are in A/m.Left and right panels show x-and y-components of the currents, respectively.

Figure 3 .
Figure 3. (a) Snapshots of the magnitude and direction of the Spherical Elementary Current Systems (SECS) and CB equivalent currents and equivalent currents' time derivatives above Fennoscandia at 20:08:30 UT, 30 October 2003.International Monitor for Auroral Geomagnetic Effects sites are marked with red asterisks and white circles.(b) Time series of the SECS and CB ionospheric equivalent currents above Abisko (ABK) and Nurmijärvi (NUR) geomagnetic observatories.The left and right panels show x-and y-components of the currents, respectively.

Figure 4 .
Figure 4. (a and b) Snapshots of the magnitude and direction the ground electric field (GEF) in Fennoscandia modeled using Spherical Elementary Current Systems and CB equivalent currents, at 20:08:30 UT, 30 October 2003.(c) The absolute difference between magnitudes of the GEF in (a) and (b).(d) The magnitude and direction of the GEF in the area of the Finnish natural gas pipeline modeled with the use of the CB source at 20:08:30 UT, 30 October 2003.(e-g) Conductivity distribution in three layers of the 3-D part of the conductivity model.The Finnish natural gas pipeline network is marked in all figures.International Monitor for Auroral Geomagnetic Effects sites Abisko (ABK) and Nurmijärvi (NUR) as well as Mäntsälä (MAN) pipeline GIC recording site and Point X located 0.5° north of MAN are marked with white circles.

Figure 5 .
Figure 5. Snapshots of the magnitude and direction of the ground electric field across Fennoscandia (left) and in the Finnish natural gas pipeline area (right) modeled using the CB source at 06:57:30 UT, 29 October 2003 (top), and 07:27:00 UT, 29 October 2003 (bottom).The Finnish natural gas pipeline network is marked in all figures.International Monitor for Auroral Geomagnetic Effects sites Abisko (ABK) and Nurmijärvi (NUR) as well as Mäntsälä (MAN) pipeline GIC recording site and Point X located 0.5° north of MAN are marked with white circles.

Figure 7 .
Figure 7.The same caption as in Figure 6 but for the Nurmijärvi (NUR) geomagnetic observatory.

Figure 8 .
Figure8.The same caption as in Figure6but for the Mäntsälä (MAN) natural gas pipeline recording site.

Figure 9 .
Figure 9.The same caption as in Figure 6 but for Point X located 0.5° north of Mäntsälä natural gas pipeline GIC recording site.

Figure 10 .
Figure 10.(a) Time series of observed GIC at the Mäntsälä (MAN) natural gas pipeline recording point.and c) Time series of GIC observed and at at time intervals 06:00:00-08:00:00 UT, 29 October and 19:00:00-21:00:00 UT, 30 October 2003.(d) Time series of GIC observed and modeled at Point X located 0.5° north of MAN.Note that modeled GIC are scaled by a factor of 4.51 in Figures (a-c), while there is no scaling in Figure(d).GIC in all figures were calculated based on the ground electric field modeled using the CB source.

Table 1
Correlation Coefficients, Normalized Root Mean Square Errors, and Maximum Absolute Differences (in mV/km) Between the Ground Electric Field Obtained With the Use of Spherical Elementary Current Systems and CB Inducing Sources and Simulated Using the MT Approach at Abisko (ABK) and Nurmijärvi (NUR) Geomagnetic Observatories, Mäntsälä (MAN) GIC Recording Point, and Point X Located 0.5° North of MAN Note.Modeled GIC at MAN is scaled by factors of 4.68 (SECS source), 4.51 (CB source), and 5.39 (MT method).At Point X, GIC are not scaled.

Table 2
Correlation Coefficients, Normalized Root Mean Square Errors, and Maximum Absolute Differences (in A) Between Modeled and Observed GIC at the Mäntsälä (MAN) Pipeline Recording Point and Point X Located 0.5° North of MAN

Appendix A: Expressions for 𝑨𝑨  𝒏𝒏 𝐄𝐄 𝒊𝒊 , 𝑨𝑨  𝒏𝒏 𝐁𝐁 𝒊𝒊
Kruglyakov, Kuvshinov, and Nair (2022)and Marshalko (2022)andKruglyakov, Kuvshinov, and Nair (2022)for electric and magnetic fields, correspondingly.Even though the same technique is used in both papers, the notation and final form of expressions are different.Thus, to avoid readers' confusion, we present the expressions for    ,   the expressions are It is worth stressing here that unlike Equation A3, Equation A4 contains an additional term  Re(, ; )| =0 for n = 0 since the magnetic field is not necessarily zero for ω = 0. (A5) EM and AV were supported by Grants 314670 and 339329 from the Academy of Finland.MK was supported by the New Zealand Ministry of Business, Innovation, & Employment through Endeavour Fund Research Programme contract UOOX2002.AK was supported in the framework of Swarm DISC activities, funded by ESA contract no.4000109587, with support from EO Science for Society.