Estimation of the 3-D geoelectric field at the Earth's surface using spherical elementary current systems

The geoelectric field drives geomagnetically induced currents (GICs) in technological conductor networks, which can affect the performance of critical ground infrastructure such as electric power transmission grids. The three-dimensional (3-D) electric field at the Earth's surface consists of an external divergence-free (DF) part due to temporally and spatially varying ionospheric and magnetospheric currents, an internal DF part due to temporally and spatially varying telluric currents, and a curl-free (CF) part due to charge accumulation at ground conductivity gradients. We have developed a new method for estimating these contributions. The external and internal parts of the DF electric field are calculated from the time derivative of the external and internal parts of the observed ground magnetic field, respectively, using DF two-dimensional (2-D) spherical elementary current systems (SECSs). The horizontal surface CF electric field is calculated from the known surface DF electric field using coefficients that linearly relate the DF electric field to the CF electric field. The coefficients were obtained from the 3-D induction model PGIEM2G (Polynomial Galerkin Integral Equation Modelling in ElectroMagnetic Geophysics). The calculations are carried out in the time domain, and only two consecutive time steps of the observed magnetic field are needed to compute the surface electric field. The external part of the DF electric field is valid at and below the ionosphere, the internal part at and above the Earth's surface, and the CF part at the Earth's surface. A dense magnetometer network is a requirement for reliable results. The external and internal parts of the DF electric field are generally oppositely directed and have comparable amplitudes, both on the ground and in the ionosphere, indicating that both contributions are significant for the total DF electric field. The largest peaks of the total DF electric field tend to occur when either the external or internal contribution is temporarily suppressed. At a given location, a DF electric field with a given amplitude can result in a total surface electric field amplitude with a difference of orders of magnitude depending on the direction of the DF electric field with respect to the locally dominant conductivity gradient structure. The electric field calculation is computationally light, facilitating operational implementation of a near-real-time 3-D surface electric field monitoring and derivation of long electric field time series.

1 Introduction

The geoelectric field at the Earth's surface drives geomagnetically induced currents (GICs) in technological conductor networks, which can affect the performance of critical ground infrastructure such as electric power transmission grids (Pulkkinen et al., 2017). In order to avoid incidents caused by space weather, a solid understanding of the possible amplitude range of the highly variable and location-dependent electric field is needed. Because electric field measurements are scarce, the surface electric field is typically modeled using magnetic field measurements and knowledge of the conductivity of the Earth. For the same reason, validation of the modeling results is difficult, which means that it would be useful to have several alternative approaches that would allow intercomparison of various modeling results.

Several modeling techniques for estimating the surface electric field exist, either based on one-dimensional (1-D) ground conductivity models (Viljanen et al., 2004, 2012, 2014), magnetotelluric transfer functions (Kelbert et al., 2017; Love et al., 2018; Lucas et al., 2020; Malone-Leigh et al., 2023), or three-dimensional (3-D) ground conductivity models (Rosenqvist and Hall, 2019; Marshalko et al., 2021; Kruglyakov et al., 2022; Marshalko et al., 2023; Kruglyakov et al., 2023). 1-D modeling only considers vertical conductivity gradients and consequently ignores lateral gradients, which are now known to be highly significant (Marshalko et al., 2021). The magnetotelluric transfer functions take lateral conductivity gradients into account but ignore spatial variations in the source, which are also known to be important (Lucas et al., 2018; Marshalko et al., 2021). 3-D electromagnetic modeling considers both vertical and lateral conductivity gradients and can also take into account spatial variations of the inducing source, but as a result of sparse magnetometer coverage and inaccuracies in 3-D conductivity models, the internal part of the magnetic field variations is not conserved (Kruglyakov et al., 2023; Marshalko et al., 2023). They can also be computationally expensive.

The modeling techniques can be broadly divided into two categories: first-principles modeling and empirical modeling. In the case of first-principles modeling, the external driver (i.e., currents in space or their magnetic field on the ground) and ground conductivity are given. The problem is fully controlled and can be applied to simulated external drivers as well as observations, providing a comprehensive way to understand the response of the conducting ground to external driving. In principle, it is possible to solve the electric field, telluric current density, and charge density everywhere, although in practice this is not currently feasible due to the computational expense. In the case of empirical modeling, ground-based magnetic field observations are used to estimate the electric field. Information on the external driver and ground conductivity is contained in the observed magnetic field, but further information, such as transfer functions, is needed to estimate the total geoelectric field. The magnetic field can be uniquely separated into an external and internal part and interpreted in terms of equivalent currents (Juusola et al., 2020). Empirical modeling is often useful in practical GIC applications.

Like any vector field on a sphere, the geoelectric field at the Earth's surface can be separated into divergence-free (DF) and curl-free (CF) parts. The DF electric field is directly associated with temporal variations of the magnetic field via Faraday's law and can be further separated into an external part due to time-varying ionospheric and magnetospheric currents in space and an internal part due to telluric currents in the conducting Earth. The CF part is created by electric charges accumulated at conductivity gradients in the ground. There can also be charges in the ionosphere and upper space. However, they are quasi-static, and the related electric field is confined between the ionosphere and Earth's surface and is perpendicular to the surface.

Vanhamäki et al. (2013) developed a method for calculating the DF surface electric field using the two-dimensional (2-D) spherical elementary current system (SECS) method (Amm, 1997; Amm and Viljanen, 1999; Pulkkinen et al., 2003a, b; McLay and Beggan, 2010; Weygand et al., 2011; Juusola et al., 2016; Marsal et al., 2017, 2020; Juusola et al., 2020; Laundal et al., 2021; Vanhamäki and Juusola, 2020; Walker et al., 2023). Their approach was to place a DF 2-D SECS layer in the ionosphere at 100 km altitude and to determine the SECS amplitudes by fitting the superposed vertical magnetic field components (B_r) of the SECSs to the measured vertical magnetic field components on the ground. The horizontal component of the vector potential (A_h) was derived by applying the equation for the vector potential of the DF 2-D SECS from Amm and Viljanen (1999). The DF electric field was calculated from the time derivative of the vector potential as

$$\begin{array}{}\text{(1)}& {\mathit{E}}_{\mathrm{DF},\mathrm{h}}=-{\displaystyle \frac{\partial {\mathit{A}}_{\mathrm{h}}}{\partial t}}.\end{array}$$

They also discussed the possibility of deriving the CF part of the surface electric field from the known ground surface conductivity and the requirement that the total current is divergence-free.

We will continue the work by Vanhamäki et al. (2013) by deriving a method for calculating the internal and external part of the DF electric field at and above the Earth's surface and the CF electric field at the Earth's surface. Similar to Vanhamäki et al. (2013), we will utilize the DF 2-D SECSs to calculate the DF electric field, but in addition to the ionospheric DF SECS layer, we will place a second layer just below the Earth's surface to represent the magnetic field of the telluric currents. A ground conductivity model is not needed for the DF electric field calculation because the conductivity distribution affects the telluric current distribution, and this is reflected in the magnetic field it produces. Using two layers instead of one allows us to fit all three components of the measured magnetic field when determining the SECS amplitudes, to derive a solution that is valid not only at the Earth's surface but also above it, and to separate the internal and external part of the DF electric field.

We will implement and test the calculation for the CF electric field suggested by Vanhamäki et al. (2013). However, we will show that the required assumption of vanishing radial currents is too restrictive and instead derive coefficients based on the first-principles 3-D induction model PGIEM2G (Polynomial Galerkin Integral Equation Modelling in ElectroMagnetic Geophysics) (Kruglyakov and Kuvshinov, 2018) that linearly relate the DF electric field to the horizontal CF electric field at the Earth's surface. PGIEM2G is a scalable 3-D electromagnetic forward modeling code based on a method of volume integral equations with a contracting kernel (Kruglyakov and Kuvshinov, 2018). The resulting geoelectric field will be compared with (a) the geoelectric field obtained via fully 3-D simulation using SECS-based ionospheric equivalent currents and PGIEM2G and (b) the geoelectric field modeled via the multi-site transfer function (MSTF) approach (Kruglyakov et al., 2023). The MSTF approach uses transfer functions that relate the horizontal geoelectric field to the horizontal magnetic field in the modeling region. The MSTFs have been calculated on the basis of electric and magnetic fields produced by PGIEM2G. The MSTF approach partly mitigates inaccuracies in the external source model, improving the accuracy of the resulting geoelectric field. We will also demonstrate how a simple linear relation, based on the same principle as the new method, can be used to relate a local magnetic field time derivative to the local geoelectric field. This application is not as flexible in terms of magnetometer configurations as the SECS method, but it can be used when the magnetometer network is too sparse for the SECS method to work reliably.

3-D induction modeling codes such as PGIEM2G determine total electric and magnetic fields by solving Maxwell's equations numerically for a given 3-D conductivity model and given inducing source. Such an approach requires relatively long time integration, during which inaccuracies can accumulate, and is sensitive to the available conductivity and source models. Our approach has the advantage that in addition to the external part of the DF electric field, the internal part of the DF electric field is also derived from measurement. Thus, no conductivity model is needed to determine the DF part of the electric field (information on the conductivity is automatically included in the measured magnetic field) and no time integration is required. However, the sensitivity to the magnetometer coverage remains, as does the dependence of the CF part of the electric field on the conductivity model.

Although first-principles 3-D induction modeling codes are very rigorous and can model ground induction as far as limited magnetic field observations and 3-D ground conductivity models allow, this does not mean that induction is perfectly understood. The geoelectric field is deterministic but highly unpredictable. In other words, theoretically speaking, the electric field can be precisely calculated from a given external source and a known ground conductivity. In practice, however, it is difficult to give any simple rules for how the ground responds to a large variety of different external drivers such as sudden impulses, temporally varying large-scale electrojets, and current vortices. As Juusola et al. (2020) show, the internal part of the time derivative of the horizontal ground magnetic field ($\mathrm{d}\mathit{H}/\mathrm{d}t$) has a much more complex spatial structure than the relatively smooth external $\mathrm{d}\mathit{H}/\mathrm{d}t$. The contrast becomes even more pronounced when comparing the external driver and the geoelectric field (Marshalko et al., 2023). Separation of the various components of the surface electric field (external and internal parts of the DF electric field and the CF electric field) provides a tool for deepening our understanding, including the combination of external driving and ground conductivity structures that pose the most severe risks for technological conductor networks. Such understanding is particularly useful for space weather services and power grid operators. Furthermore, computationally reasonable methods allow the production of long time series of surface electric field data, which can be used as material for extreme value analysis to better understand the threats of severe space weather events.

The structure of the study is as follows: the new methods are derived and old methods reviewed in Sect. 2, test results are presented and validated in Sect. 3, and further applications are discussed in Sect. 4. The conclusions are summarized in Sect. 5.

2 Methods

In this section, we will first discuss the division of the geoelectric field into DF and CF parts. In Sect. 2.2, we will derive the method for obtaining the external and internal parts of the DF electric field at and above the Earth's surface from ground magnetic field measurements using DF 2-D SECSs. In Sect. 2.3.1 we will derive the method for estimating the CF surface electric field from a known DF electric field, surface conductance, and vertical current density using CF 2-D SECS functions. In Sect. 2.3.1, we will derive the alternative method for estimating the CF surface electric field from a known DF electric field using coefficients derived from first-principles 3-D induction modeling data. In Sect. 2.4 we will describe how the principles utilized in the DF and CF electric field calculation allow a simple estimation of the total geoelectric field from the time derivative of the magnetic field. Unlike the SECS method, such an approach is rigid in terms of magnetometer station configuration, but it can be convenient when the magnetometer network is too sparse for the SECS method to work reliably. Finally, in Sect. 2.5, we will describe the existing methods, PGIEM2G modeling with the SECS-based source and the MSTF approach, which we use to derive coefficients and to compare with the new models.

2.1 Divergence-free and curl-free parts of the geoelectric field

We will approach the modeling problem by separating the different contributions to the geoelectric field. As follows from the Maxwell equations, the electric field E can be expressed as

$$\begin{array}{}\text{(2)}& \mathit{E}=-{\displaystyle \frac{\partial \mathit{A}}{\partial t}}-\mathrm{\nabla }\mathit{\varphi },\end{array}$$

where A is the vector potential and ϕ is the scalar potential. It is possible to require the vector potential to be divergence-free:

$$\begin{array}{}\text{(3)}& \mathrm{\nabla }\cdot \mathit{A}=\mathrm{0}.\end{array}$$

This is known as the Coulomb gauge condition (Jackson, 1998). It follows that the scalar potential can be expressed in terms of the charge density ρ as

$$\begin{array}{}\text{(4)}& \mathit{\varphi }(\mathit{r},t)={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }{\mathit{ϵ}}_{\mathrm{0}}}}\underset{V}{\int }{\displaystyle \frac{\mathit{\rho }({\mathit{r}}^{\prime },t)}{|\mathit{r}-{\mathit{r}}^{\prime }|}}\mathrm{d}{V}^{\prime },\end{array}$$

where ϵ₀ is the vacuum permittivity. Contrary to the often applied Lorenz gauge, the time t is not retarded, but the scalar potential is formally identical to the static case.

In the Coulomb gauge, the vector potential is given by the divergence-free (DF) part of the current density J,

$$\begin{array}{}\text{(5)}& \mathit{A}(\mathit{r},t)={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }}}\underset{V}{\int }{\displaystyle \frac{{\mathit{J}}_{\mathrm{DF}}({\mathit{r}}^{\prime },t-|\mathit{r}-{\mathit{r}}^{\prime }|/c)}{|\mathit{r}-{\mathit{r}}^{\prime }|}}\mathrm{d}{V}^{\prime },\end{array}$$

where μ₀ is the vacuum permeability and c is the speed of light. Note that here the time is retarded.

The power of the Coulomb gauge is that it automatically separates the electric field into a curl-free (CF) static part and a divergence-free inductive part:

$$\begin{array}{}\text{(6)}& \mathit{E}={\mathit{E}}_{\mathrm{CF}}+{\mathit{E}}_{\mathrm{DF}},\end{array}$$

where

$$\begin{array}{}\text{(7)}& {\mathit{E}}_{\mathrm{CF}}=-\mathrm{\nabla }\mathit{\varphi },{\mathit{E}}_{\mathrm{DF}}=-{\displaystyle \frac{\partial \mathit{A}}{\partial t}},\end{array}$$

and E_CF is produced only by charges and E_DF only by divergence-free currents as explicitly shown above.

2.2 DF electric field of a DF 2-D SECS with a time-varying amplitude

The DF part of the electric field (E) due to a time-varying magnetic field (B) can be solved from Faraday's law:

$$\begin{array}{}\text{(8)}& \mathrm{\nabla }\times {\mathit{E}}_{\mathrm{DF}}=-{\displaystyle \frac{\partial \mathit{B}}{\partial t}}.\end{array}$$

DF 2-D SECSs form a complete set of basis functions for representing ground magnetic field variations in terms of external and internal equivalent currents (Vanhamäki and Juusola, 2020). The time-varying sheet current density of the DF 2-D SECS (Amm, 1997) is

$$\begin{array}{}\text{(9)}& {\mathit{J}}_{\mathrm{DF}}\left(t\right)={\displaystyle \frac{I_{\mathrm{DF}}\left(t\right)}{\mathrm{4}\mathit{\pi }R}}\cot \left({\displaystyle \frac{{\mathit{\theta }}^{\prime }}{\mathrm{2}}}\right){\widehat{\mathit{e}}}_{{\mathit{\varphi }}^{\prime }},\end{array}$$

where I_DF(t) is the time-varying amplitude in A, R is the radius of the current sheet, and θ^′ is the co-latitude in the SECS coordinates, where ${\mathit{\theta }}^{\prime }=\mathrm{0}$ at the SECS pole. Conversion to general coordinates, e.g., geographic, follows the regular practice of the 2-D SECS method (Vanhamäki and Juusola, 2020). It should be noted that although the resulting external DF current density can only be interpreted as the DF part of the ionospheric horizontal current density at high latitudes, where the ambient magnetic field lines can be approximated as radial, it can still describe the ground magnetic field everywhere below the ionosphere (e.g., Juusola et al., 2023a, b). Thus, the applicability of the DF electric field calculation method is not limited to high latitudes.

In the SECS coordinates, the time derivative of the DF 2-D SECS magnetic field only has r and θ^′ components (Amm and Viljanen, 1999):

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\partial \mathit{B}}{\partial t}}={\displaystyle \frac{\partial B_r}{\partial t}}{\widehat{\mathit{e}}}_r+{\displaystyle \frac{\partial B_{{\mathit{\theta }}^{\prime }}}{\partial t}}{\widehat{\mathit{e}}}_{{\mathit{\theta }}^{\prime }},\end{array}$$

where

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\displaystyle \frac{\partial B_r}{\partial t}}& ={\displaystyle \frac{\partial I_{\mathrm{DF}}}{\partial t}}{\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{\mathrm{1}}{r}}\\ & \times \left({\displaystyle \frac{R}{\sqrt{r^{\mathrm{2}}-\mathrm{2}rR\cos {\mathit{\theta }}^{\prime }+R^{\mathrm{2}}}}}+\left\{\begin{array}{c}-\mathrm{1}\\ -R/r\end{array}\right\}\right)\begin{array}{c}\text{ when }r<R\\ \text{ when }r\ge R\end{array}\end{array}\end{array}$$

and

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\displaystyle \frac{\partial B_{{\mathit{\theta }}^{\prime }}}{\partial t}}& =-{\displaystyle \frac{\partial I_{\mathrm{DF}}}{\partial t}}{\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{\mathrm{1}}{r\sin {\mathit{\theta }}^{\prime }}}\\ & \times \left({\displaystyle \frac{r-R\cos {\mathit{\theta }}^{\prime }}{\sqrt{r^{\mathrm{2}}-\mathrm{2}rR\cos {\mathit{\theta }}^{\prime }+R^{\mathrm{2}}}}}+\left\{\begin{array}{c}\cos {\mathit{\theta }}^{\prime }\\ -\mathrm{1}\end{array}\right\}\right)\begin{array}{c}\text{ when }r<R\\ \text{ when }r>R\end{array}.\end{array}\end{array}$$

This follows from neglecting the displacement current, as is usual in geoelectromagnetism. However, it should be noted that the displacement current may play a role in producing ground-based magnetic field signatures that do not cause a geoelectric field on the ground (Brändlein et al., 2012). Due to the geometry of the DF SECS magnetic field, the corresponding induced electric field only has a ϕ^′ component (see also the vector potential derivation in Amm and Viljanen, 1999),

$$\begin{array}{}\text{(13)}& {\mathit{E}}_{\mathrm{DF}}=E_{\mathrm{DF},{\mathit{\varphi }}^{\prime }}{\widehat{\mathit{e}}}_{{\mathit{\varphi }}^{\prime }},\end{array}$$

which can be solved from the integral form of Eq. (8),

$$\begin{array}{}\text{(14)}& \underset{\partial A}{\oint }{\mathit{E}}_{\mathrm{DF}}\cdot \mathrm{d}\mathit{l}=-{\displaystyle \frac{\partial }{\partial t}}\underset{A}{\int }\mathit{B}\cdot \mathrm{d}\mathit{S},\end{array}$$

by calculating the change in the magnetic flux through the spherical cap surface defined by a field line of $E_{\mathrm{DF},{\mathit{\varphi }}^{\prime }}$ (Fig. 1):

$$\begin{array}{}\text{(15)}& \mathrm{2}\mathit{\pi }r\sin {\mathit{\theta }}^{\prime }{E}_{\mathrm{DF},{\mathit{\varphi }}^{\prime }}=-r^{\mathrm{2}}\underset{\mathrm{0}}{\overset{{\mathit{\theta }}^{\prime }}{\int }}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\mathrm{d}\mathit{\varphi }{\displaystyle \frac{\partial B_r}{\partial t}}.\end{array}$$

Figure 1 Geometry of the derivation of the divergence-free (DF) electric field (E_DF,ϕ) of a DF 2-D spherical elementary current system (SECS) with a time-varying amplitude. J_ϕ indicates the SECS current density, B its magnetic field, and $\mathrm{d}/\mathrm{d}t$ the time derivative. r, θ, and ϕ are spherical coordinates (radius, co-latitude, and longitude), and $\widehat{\mathit{r}}$, $\widehat{\mathit{\theta }}$, and $\widehat{\mathit{\varphi }}$ are the corresponding unit vectors.

Carrying out the calculation yields

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{E}_{\mathrm{DF},{\mathit{\varphi }}^{\prime }}& =-{\displaystyle \frac{\partial I_{\mathrm{DF}}}{\partial t}}{\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{\mathrm{1}}{r\sin {\mathit{\theta }}^{\prime }}}\\ & \times \left(\sqrt{r^{\mathrm{2}}-\mathrm{2}rR\cos {\mathit{\theta }}^{\prime }+R^{\mathrm{2}}}+\left\{\begin{array}{c}r\cos {\mathit{\theta }}^{\prime }-R\\ R\cos {\mathit{\theta }}^{\prime }-r\end{array}\right\}\right)\begin{array}{c}\text{ when }r<R\\ \text{ when }r\ge R\end{array},\end{array}\end{array}$$

which is equal to the negative time derivative of the vector potential (Eq. 1) given by Amm and Viljanen (1999) for the DF SECS when r<R. Comparison of Eq. (16) with Eq. (9) confirms that the induced electric field opposes the temporal change in the DF current, as expected from Lenz's law.

Figure 2 shows co-latitude profiles of $\partial J_{\mathit{\varphi }}/\partial t$, ground $\partial B_r/\partial t$, ground $\partial B_{\mathit{\theta }}/\partial t$, and ground E_DF,ϕ for two different cases, J_ϕ at 90 km altitude (blue curves, external contribution) and oppositely directed J_ϕ at 2 × 90 km depth (red curves, internal contribution), to mimic a simple case of ionospheric and telluric currents. Both have the same amplitude: $|\partial I_{\mathrm{DF}}/\partial t|=\mathrm{10}$ kA s⁻¹. The sum of the external and internal $\partial B_r/\partial t$, $\partial B_{\mathit{\theta }}/\partial t$, and E_ϕ on the ground is shown by the black curves. While the current density is concentrated near the SECS pole at θ=0, the magnetic field is more spread out to lower latitudes and the DF electric field even more so. On the ground, the horizontal magnetic field components B_θ from the oppositely directed ionospheric and telluric currents strengthen each other, while the vertical components B_r tend to cancel each other out. The DF electric fields also tend to cancel each other out, resulting in a much weaker total ground DF electric field than either the original ionospheric or telluric contribution.

Brändlein et al. (2012) discussed a waveguide transmission, where the wave mode on the ground has a nonzero horizontal magnetic field component but a zero horizontal electric field. In the vertical direction such a wave mode is expected to have a zero magnetic component and nonzero electric field component. The SECS reconstruction is able to reproduce such a magnetic field as a superposition of the magnetic fields of ionospheric and telluric equivalent currents. Because the vertical magnetic field is zero, E_DF would also be zero.

Although the combination of ionospheric and telluric DF current densities always produces a DF electric field that has a zero vertical component between the ionospheric and telluric equivalent current sheets, this does not necessarily mean that E_DF cannot have a vertical component in this region. This issue was investigated in detail by Pirjola and Viljanen (1998). In addition to the parts described by the DF current densities, the 3-D current distributions in the ionosphere and in the ground include a part that has a zero magnetic field between the ionosphere and ground surface. However, the corresponding vector potential A may not be zero, although ∇×A must be zero. The corresponding DF electric field is a Laplace field that has its sources above the ionosphere and inside the ground. A similar Laplace electric field could also be produced with electric charges in these regions. According to the results by Pirjola and Viljanen (1998), valid up to neglecting the displacement current, any horizontal part of this DF field is canceled by charges accumulated on the ground surface, leaving only an insignificant vertical component. Thus this part of the induction process does not drive any GIC.

Figure 2 Co-latitude profiles of $\partial J_{\mathit{\varphi }}/\partial t$, ground $\partial B_r/\partial t$, ground $\partial B_{\mathit{\theta }}/\partial t$, and ground E_DF,ϕ for two different cases: J_ϕ of a DF 2-D SECS at 90 km altitude (blue curves, external contribution) and oppositely directed J_ϕ at 2 × 90 km depth (red curves, internal contribution). Both have the same rate of change of the SECS amplitude $|\partial I_{\mathrm{DF}}/\partial t|=\mathrm{10}$ kA s⁻¹. The sum of the external and internal $\partial B_r/\partial t$, $\partial B_{\mathit{\theta }}/\partial t$, and E_DF,ϕ on the ground is shown by the black curves.

The time derivative of the SECS amplitude can be estimated from the known amplitudes I_DF as

$$\begin{array}{}\text{(17)}& {\displaystyle \frac{\partial I_{\mathrm{DF}}}{\partial t}}={\displaystyle \frac{I_{\mathrm{DF}}\left(t\right)-I_{\mathrm{DF}}(t-\mathrm{\Delta }t)}{\mathrm{\Delta }t}},\end{array}$$

where Δt is the time step of the data. When both the ionospheric equivalent current amplitudes at 90 km altitude and telluric equivalent current amplitudes at 1 m depth have been determined from ground magnetic field observations (Juusola et al., 2016, 2020; Vanhamäki and Juusola, 2020), Eq. (16) can be used to estimate the ionospheric DF current contribution to E_DF everywhere and the telluric current contribution to E_DF at and above the Earth's surface.

2.3 CF electric field at the Earth's surface driven by the DF electric field

2.3.1 Approach A: CF SECS

The DF electric field in the conducting Earth has an external and an internal component. The external component is produced by time-varying ionospheric and magnetospheric current systems and the internal component by time-varying 3-D telluric currents. This electric field drives the electric currents in the Earth. Whenever the current driven by the DF electric field crosses a conductivity gradient, there is divergence of the current. Because the total current must be divergence-free (follows from neglecting the displacement current, as usual in geoelectromagnetism), a curl-free (CF) electric field is set up by redistribution of electric charges in the Earth to satisfy this condition. The divergence of the total current can be written as

$$\begin{array}{}\text{(18)}& \begin{array}{rl}\mathrm{\nabla }\cdot \mathit{j}& =\mathrm{\nabla }\cdot (\mathit{\sigma }{\mathit{E}}_{\mathrm{DF}}+\mathit{\sigma }{\mathit{E}}_{\mathrm{CF}})=\mathrm{\nabla }\mathit{\sigma }\cdot {\mathit{E}}_{\mathrm{DF}}\\ & +\mathit{\sigma }\mathrm{\nabla }\cdot {\mathit{E}}_{\mathrm{DF}}+\mathrm{\nabla }\mathit{\sigma }\cdot {\mathit{E}}_{\mathrm{CF}}+\mathit{\sigma }\mathrm{\nabla }\cdot {\mathit{E}}_{\mathrm{CF}},\end{array}\end{array}$$

where $\mathit{\sigma }(r,\mathit{\theta },\mathit{\varphi })$ is the known, location-dependent electrical conductivity. The DF electric field is a superposition of the external and internal contribution and satisfies

$$\begin{array}{}\text{(19)}& \mathrm{\nabla }\cdot {\mathit{E}}_{\mathrm{DF}}=\mathrm{0}.\end{array}$$

We will only consider the surface layer of the conducting Earth and model it as an infinitely thin conducting sheet located at the Earth's surface. Equation (18) then becomes

$$\begin{array}{}\text{(20)}& \mathrm{\nabla }\mathrm{\Sigma }\cdot {\mathit{E}}_{\mathrm{CF}}+\mathrm{\Sigma }\mathrm{\nabla }\cdot {\mathit{E}}_{\mathrm{CF}}=-\mathrm{\nabla }\mathrm{\Sigma }\cdot {\mathit{E}}_{\mathrm{DF}}+j_r,\end{array}$$

where Σ is the conductance of the thin sheet, obtained by vertical integration of the conductivity in the surface layer, and j_r is the vertical current density.

We use CF SECS functions to obtain E_CF when Σ and E_DF are known. The electric field of a CF SECS (Vanhamäki and Juusola, 2020) is

$$\begin{array}{}\text{(21)}& {\mathit{E}}_{\mathrm{CF}}={\displaystyle \frac{Q_{\mathrm{CF}}}{\mathrm{4}\mathit{\pi }R}}{\widehat{\mathit{e}}}_{{\mathit{\theta }}^{\prime }}\left\{\begin{array}{c}\cot \left({\displaystyle \frac{{\mathit{\theta }}^{\prime }}{\mathrm{2}}}\right)\\ \tan \left(\frac{{\mathit{\theta }}^{\prime }}{\mathrm{2}}\right){\cot }^{\mathrm{2}}\left(\frac{{\mathit{\theta }}_{\mathrm{0}}}{\mathrm{2}}\right)\end{array}\right\}\begin{array}{c}\text{ when }{\mathit{\theta }}^{\prime }\ge {\mathit{\theta }}_{\mathrm{0}}\\ \text{ when }{\mathit{\theta }}^{\prime }<{\mathit{\theta }}_{\mathrm{0}}\end{array},\end{array}$$

where Q_CF is the amplitude of the CF SECS in V (i.e., Q_CF is charge divided by ϵ₀), R is the radius of the surface, and θ^′ is the co-latitude in the SECS coordinates. The source at the elementary current system's pole is spread uniformly inside a spherical cap of width θ₀. The divergence of E_CF is given by

$$\begin{array}{}\text{(22)}& \mathrm{\nabla }\cdot {\mathit{E}}_{\mathrm{CF}}={\displaystyle \frac{\mathit{\rho }}{{\mathit{ϵ}}_{\mathrm{0}}}}=-{\displaystyle \frac{Q_{\mathrm{CF}}}{\mathrm{4}\mathit{\pi }{R}^{\mathrm{2}}}}+\left\{\begin{array}{c}\mathrm{0}\\ \frac{Q_{\mathrm{CF}}}{A}\end{array}\right\}\begin{array}{c}\text{ when }{\mathit{\theta }}^{\prime }\ge {\mathit{\theta }}_{\mathrm{0}}\\ \text{ when }{\mathit{\theta }}^{\prime }<{\mathit{\theta }}_{\mathrm{0}}\end{array},\end{array}$$

where ρ is the charge density and A is the area of the spherical cap defined by θ₀:

$$\begin{array}{}\text{(23)}& A=\mathrm{2}\mathit{\pi }{R}^{\mathrm{2}}(\mathrm{1}-\cos {\mathit{\theta }}_{\mathrm{0}}).\end{array}$$

The equations above apply to a single SECS pole. The total CF electric field in the grid cell i is a superposition of the electric fields of all the CF SECSs. Equation (20) can now be written as

$$\begin{array}{}\text{(24)}& \begin{array}{rl}\sum _j& \left[\mathrm{\nabla }{\mathrm{\Sigma }}_{\mathit{\theta },i}{g}_{\mathit{\theta },i,j}+\mathrm{\nabla }{\mathrm{\Sigma }}_{\mathit{\varphi },i}{g}_{\mathit{\varphi },i,j}-{\mathrm{\Sigma }}_i{\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }{R}^{\mathrm{2}}}}\right]Q_{\mathrm{CF},j}\\ & +{\mathrm{\Sigma }}_i{\displaystyle \frac{Q_{\mathrm{CF},i}}{A_i}}=-\mathrm{\nabla }{\mathrm{\Sigma }}_{\mathit{\theta },i}{E}_{\mathrm{DF}\mathit{\theta },i}-\mathrm{\nabla }{\mathrm{\Sigma }}_{\mathit{\varphi },i}{E}_{\mathrm{DF}\mathit{\varphi },i}+j_{r,i},\end{array}\end{array}$$

where we have denoted the geometry-dependent components of E_CF by g. By assuming that the radial current either vanishes ($j_{r,i}=\mathrm{0}$) or is otherwise known, the amplitudes Q_CF,i can be determined from the matrix equation corresponding to Eq. (24). No time integration of currents is needed because the internal contribution to E_DF is obtained from measurements. Thus, only two time steps of data are needed to determine the geoelectric field at the Earth's surface for a given epoch.

2.3.2 Approach B: linear coefficients relating E_CF to E_DF

The CF 2-D SECS-based method described above assumes that the vertical current density in the surface layer is known. We will demonstrate below that this assumption is too restrictive. As an alternative approach, we assume that the CF electric field is approximately proportional to the driving DF electric field.

$$\begin{array}{}\text{(25)}& {\displaystyle }{E}_{\mathrm{CF},x}=k_{xx}{E}_{\mathrm{DF},x}+k_{xy}{E}_{\mathrm{DF},y}\text{(26)}& {\displaystyle }{E}_{\mathrm{CF},y}=k_{yy}{E}_{\mathrm{DF},y}+k_{yx}{E}_{\mathrm{DF},x}\end{array}$$

This formulation resembles that of magnetotelluric transfer functions, which define the frequency-dependent linear relationship (impedance tensors) between components of the electric field and horizontal components of the magnetic field variations measured at a given location. Our formulation does not depend on frequency but is given in the time domain. This is possible because the internal part of the DF electric field is obtained from measurements and we only need a linear relationship for describing the CF part of the electric field from the known DF part. The assumed linear relation contains a local plane wave assumption; i.e., we assume that the spatial structures of the DF electric field are larger than the range of the relevant charge accumulations that contribute to the local CF field. A more general relationship would be of the form (E_CF values) = matrix ⋅ (E_DF values), with all grid points collected in the matrices. However, the simple assumptions of Eqs. (25)–(26) appear to be sufficient, as we will demonstrate below. In principle, it should be enough to determine the time-independent coefficients from a single active interval of modeled geoelectric field data using. The most obvious problem would arise from a case where a persistent small-scale DF electric field structure biases the coefficients at some locations. However, the chaotic nature of the time derivative of the magnetic field (Kellinsalmi et al., 2022) makes such a case unlikely, as long as the event chosen is sufficiently long and active.

2.4 Relating E to $\frac{\partial \mathit{B}}{\partial t}$

The linear relation between the DF electric field and $\frac{\partial \mathit{B}}{\partial t}$ on the one hand, and the linear relation between the CF and DF electric field on the other hand, implies that the geoelectric field at location r can be estimated from the magnetic field observations at locations ${\mathit{r}}_k^{\mathrm{obs}}$ as

$$\begin{array}{}\text{(27)}& \mathit{E}(\mathit{r},t)=\sum _{k=\mathrm{1}}^{N_{\mathrm{obs}}}Q(\mathit{r},{\mathit{r}}_k^{\mathrm{obs}}){\displaystyle \frac{\partial }{\partial t}}\mathit{B}({\mathit{r}}_k^{\mathrm{obs}},t),\end{array}$$

where $Q(\mathit{r},{\mathit{r}}_k^{\mathrm{obs}})$ represents 2×3 time-independent matrices. The Q matrices can be determined using, e.g., PGIEM2G simulations and bypass the need for the SECS expansion. However, unlike the SECS expansion, the Q matrices are rigid in terms of station configuration and introduce dependence on the conductivity model to the DF electric field. Determining the matrices will be a topic for a separate study.

The SECS method requires a relatively dense regional magnetometer network. For a single location with a limited time series of electric field observations and a more extended time series of magnetic field observations, the electric field time series can be approximated by defining the coefficients q:

$$\begin{array}{}\text{(28)}& {\displaystyle }{E}_x=q_{xx}{\displaystyle \frac{\partial B_x}{\partial t}}+q_{xy}{\displaystyle \frac{\partial B_y}{\partial t}}+q_{xz}{\displaystyle \frac{\partial B_z}{\partial t}},\text{(29)}& {\displaystyle }{E}_y=q_{yx}{\displaystyle \frac{\partial B_x}{\partial t}}+q_{yy}{\displaystyle \frac{\partial B_y}{\partial t}}+q_{yz}{\displaystyle \frac{\partial B_z}{\partial t}}.\end{array}$$

An applicable case might occur when a magnetotelluric measurement is carried out close to the location of a semi-permanent magnetometer station.

2.5 Existing methods

2.5.1 PGIEM2G

PGIEM2G (Polynomial Galerkin Integral Equation Modelling in ElectroMagnetic Geophysics) (Kruglyakov and Kuvshinov, 2018) is a scalable 3-D electromagnetic forward modeling code based on a method of volume integral equations with a contracting kernel (Pankratov and Kuvshinov, 2016; Kruglyakov and Kuvshinov, 2018). The conductivity model used by PGIEM2G comprises a 3-D part (upper 60 km) and a 1-D conductivity profile (below 60 km) from Kuvshinov et al. (2021). The 3-D part is based on the SMAP model (Korja et al., 2002), covers an area of 2550×2550 km², 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. This model was also exploited in Marshalko et al. (2021), Kruglyakov et al. (2022), Marshalko et al. (2023), and Kruglyakov et al. (2023). The external source used in the PGIEM2G simulation is obtained using the SECS method (see Marshalko et al., 2021, 2023). As the integral equation approach is invoked, the model space is infinite.

Calculations in PGIEM2G code are carried out in Cartesian geometry for planar Earth. To make the transition from spherical to Cartesian geometry, a map projection (transverse mercator with the latitude and longitude of the true origin at 50° N and 25° E, respectively) and interpolation onto a regular grid (about 5×5 km) are first performed for the SMAP and inducing source data. Then an inverse transformation and interpolation onto a regular grid are carried out for the resulting electric and magnetic fields to obtain the data in spherical geometry. The final spatial resolution of electric and magnetic fields is 0.03°×0.07° in latitudinal and longitudinal directions, correspondingly. The discrepancy between the SMAP grid used in some of our calculations and the denser grid used in PGIEM2G is not considered an issue in our qualitative comparison. In the one case when we wish to insert some of the PGIEM2G data directly into our model (section 3.2), nearest-neighbor interpolation is performed.

2.5.2 MSTF

The multi-site transfer function (MSTF) approach proposed by Kruglyakov et al. (2023) is based on the use of transfer functions that linearly relate the horizontal geoelectric field ${\mathit{E}}_{\mathrm{h}}=(E_x,E_y)$ at any location in the modeling region to the horizontal magnetic field ${\mathit{B}}_{\mathrm{h}}=(B_x,B_y)$ at (fixed) N locations in frequency domain:

$$\begin{array}{}\text{(30)}& {\mathit{E}}_{\mathrm{h}}(\mathit{r},\mathit{\omega };\mathit{\sigma })=\sum _{k=\mathrm{1}}^N{\mathrm{\Lambda }}_k(\mathit{r},{\mathit{r}}_k,\mathit{\omega };\mathit{\sigma }){\mathit{B}}_{\mathrm{h}}({\mathit{r}}_k,\mathit{\omega };\mathit{\sigma }).\end{array}$$

Here r is a position vector, ω is angular frequency, σ(r) is the conductivity distribution in the modeling region, and Λ_k represents the following 2×2: matrices

$$\begin{array}{}\text{(31)}& {\mathrm{\Lambda }}_k=\left(\begin{array}{c}{\mathrm{\Lambda }}_{xx,k}{\mathrm{\Lambda }}_{xy,k}\\ {\mathrm{\Lambda }}_{yx,k}{\mathrm{\Lambda }}_{yy,k}\end{array}\right).\end{array}$$

Estimation of the elements of matrices ${\mathrm{\Lambda }}_k,k=\mathrm{1},\mathrm{2},\mathrm{\dots },N$ at a given frequency ω and conductivity model σ is performed by using the modeled fields ${\mathit{B}}^{\left(l\right)}({\mathit{r}}_k,\mathit{\omega };\mathit{\sigma })$ and ${\mathit{E}}^{\left(l\right)}(\mathit{r},\mathit{\omega };\mathit{\sigma })$ corresponding to lth spatial modes j_l(r) describing the external source,

$$\begin{array}{}\text{(32)}& {\mathit{j}}^{\text{ext}}(\mathit{r},t)=\sum _{l=\mathrm{1}}^L{c}_l\left(t\right){\mathit{j}}_l\left(\mathit{r}\right),\end{array}$$

where L stands for a number of spatial modes (see Kruglyakov et al., 2022). After these fields are calculated for a given computational model σ(r) they are substituted into Eq. (30) to form the linear system for elements of Λ_k that, in turn, are solved by means of a regularized ordinary least-squares approach.

Once elements of Λ_k are estimated at a predefined number of frequencies and assuming that time series of the observed magnetic field are provided at uniform time steps δ_t, E_h(r,t_i) for any given time instant t_i=iδ_t and a given location r is calculated using the numerical scheme described in Kruglyakov et al. (2022) and Marshalko et al. (2023):

$$\begin{array}{}\text{(33)}& {\mathit{E}}_{\mathrm{h}}(\mathit{r},t_i;\mathit{\sigma })=\sum _{n=\mathrm{0}}^{N_t}\sum _{k=\mathrm{1}}^N{\mathrm{\Lambda }}^{\left(n\right)}(\mathit{r},{\mathit{r}}_k,T_{\mathrm{\Lambda }};\mathit{\sigma }){\mathit{B}}_{\mathrm{h}}^{\mathrm{obs}}({\mathit{r}}_k,t_i-n{\mathit{\delta }}_t),\end{array}$$

where ${\mathit{B}}_{\mathrm{h}}^{\mathrm{obs}}$ stands for horizontal components of the observed field in the corresponding locations, T_Λ is the so-called “memory”, and $N_t=T_{\mathrm{\Lambda }}/{\mathit{\delta }}_t$. The detailed methodology of computing ${\mathrm{\Lambda }}^{\left(n\right)}(\mathit{r},{\mathit{r}}_k,T_{\mathrm{\Lambda }};\mathit{\sigma })$ from $\mathrm{\Lambda }(\mathit{r},{\mathit{r}}_k,\mathit{\omega };\mathit{\sigma })$ is provided in Kruglyakov et al. (2022, 2023).

The MSTFs used in this paper were calculated using the same conductivity model as in previous section and the 34 spatial modes identified previously by Marshalko et al. (2023) via the time-domain principal component analysis of the SECS ionospheric equivalent current in Fennoscandia during 3 d of the Halloween geomagnetic storm (29–31 October 2003). The corresponding forward simulations were carried out via PGIEM2G.

The calculations of geoelectric field during the 7–8 September geomagnetic storm via the MSTF approach were performed using T_Λ=15 min (see Kruglyakov et al., 2022, 2023) and the data from N=15 IMAGE magnetometer network stations (ABK, AND, DOB, HAN, IVA, KEV, KIR, NUR, OUJ, RVK, SOD, SOR, TAR, TRO, and UPS).

3 Results and validation

Marshalko et al. (2021) used PGIEM2G to model the 3-D geoelectric field during the 7–8 September 2017 storm. They show the resulting horizontal geoelectric field for two epochs, 23:16 and 23:52 UT. We use the epoch 23:16 UT to test our method. Corresponding illustrations for the epoch 23:52 UT can be found in the Appendix. Marshalko et al. (2021) ran their simulation using inducing source data with 1 min temporal resolution as an input, but we compare our results with a rerun that used 10 s data.

3.1 DF electric field

We have used 10 s ground magnetic field measurements from the International Monitor for Auroral Geomagnetic Effects (IMAGE) magnetometers (Juusola et al., 2024). After correcting the data for any erroneous spikes and jumps, a 10 d sliding median baseline was subtracted from the data. DF 2-D SECS poles were placed in the IMAGE region on uniform grids with 0.5° latitude and 1° longitude resolution at 1 m depth and at 90 km altitude, and their amplitudes were determined by fitting the superposed magnetic field of the SECSs to the three components of the measured magnetic field (Vanhamäki and Juusola, 2020). The time derivative of the SECS amplitudes was derived using Eq. (17) and the DF electric field using Eq. (16) combined with the regular coordinate manipulations of the SECS method (Vanhamäki and Juusola, 2020). Because in this application of the 2-D SECS method we are only interested in $\partial \mathit{B}/\partial t$ and not B, it would also be possible to skip the baseline subtraction and directly fit $\partial \mathit{B}/\partial t$ with the DF 2-D SECS functions.

Figure 3 shows the resulting external, internal, and total time derivative of the ground magnetic field ($\partial \mathit{B}/\partial t$) and the DF electric field (E_DF) and its curl (Eq. 8) on 7 September 2017 at 23:16 UT. E_DF is smoother than $\partial \mathit{B}/\partial t$ but peaks in the same general area as the horizontal part of $\partial \mathit{B}/\partial t$, as expected. Figure 3 also demonstrates that the internal parts have much more pronounced spatial variability than the relatively smoothly varying external parts. External sources are located more than ∼ 100 km from the ground points, whereas internal sources (telluric currents) are concentrated close to the surface and strongly modulated by conductivity gradients. The direction of the internal part of E_DF is opposite to that of the external part of E_DF and has a slightly weaker amplitude. As a result, the total E_DF has a clearly weaker amplitude and different distribution than the external and internal parts. The overall amplitude level of the DF electric field of hundreds of millivolts per kilometer (mV km⁻¹) is in agreement with the theoretical examples of Vanhamäki et al. (2013). The curl (∇×E)_r equals $-\partial B_r/\partial t$ and $\partial \mathit{B}/\partial t$ is conserved. It should be noted that this does not guarantee that the separation into internal and external parts is entirely correct. The two sources may get mixed to some degree, especially in areas where the magnetometer coverage is sparse (Juusola et al., 2020).

Figure 3 External (a, d), internal (b, e), and total (c, f) time derivative of the magnetic field ($\partial \mathit{B}/\partial t$, a–c) and DF electric field (E_DF, d–f) on the ground on 7 September 2017 at 23:16 UT. The horizontal component is shown by vectors in all plots. In (a)–(c), the vertical component of $\partial \mathit{B}/\partial t$ is shown by color and in (d)–(f) the curl of E_DF. Note that the color and arrow length scales vary between panels. Locations of the IMAGE magnetometer stations used in the analysis are indicated by the black squares. Apex coordinates (Richmond, 1995; Emmert et al., 2010; Laundal et al., 2022) are indicated with the blue grid. The north, east, and down components $(B_x,B_y,B_z)$ used in the plots correspond to $(-B_{\mathit{\theta }},B_{\mathit{\varphi }},-B_r)$.

3.2 CF electric field

Similar to Marshalko et al. (2021), we have used the SMAP model (Korja et al., 2002) to estimate the ground conductivity σ. SMAP provides the conductivity with a 5^′ latitude and 5^′ longitude resolution. The conductivity distribution in the top 0–10 km layer is displayed in Fig. 4.

Figure 4 Conductivity in the 0–10 km layer (i.e., the top layer, which is 10 km thick) according to the SMAP model (Korja et al., 2002). The color bar maximum value and minimum value correspond to the largest and smallest conductivity value, respectively, placing the white color of the diverging color map in the middle of these values.

We have determined the CF SECS amplitudes on a grid where the conductance, its gradient, and grid cell area are obtained with the finite-difference approach.

$$\begin{array}{}\text{(34)}& {\displaystyle }{\mathrm{\Sigma }}_{i,j}={\mathit{\sigma }}_{i,j}\cdot \mathrm{10}\mathrm{km}\text{(35)}& {\displaystyle }{\mathrm{\Sigma }}_{i+\frac{\mathrm{1}}{\mathrm{2}},j+\frac{\mathrm{1}}{\mathrm{2}}}={\displaystyle \frac{{\mathrm{\Sigma }}_{i,j}+{\mathrm{\Sigma }}_{i+\mathrm{1},j}+{\mathrm{\Sigma }}_{i,j+\mathrm{1}}+{\mathrm{\Sigma }}_{i+\mathrm{1},j+\mathrm{1}}}{\mathrm{4}}}\text{(36)}& {\displaystyle }\mathrm{\nabla }{\mathrm{\Sigma }}_{\mathit{\theta },i+\frac{\mathrm{1}}{\mathrm{2}},j+\frac{\mathrm{1}}{\mathrm{2}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \left({\displaystyle \frac{{\mathrm{\Sigma }}_{i+\mathrm{1},j}-{\mathrm{\Sigma }}_{i,j}}{R_{\mathrm{E}}\mathrm{\Delta }\mathit{\theta }}}+{\displaystyle \frac{{\mathrm{\Sigma }}_{i+\mathrm{1},j+\mathrm{1}}-{\mathrm{\Sigma }}_{i,j+\mathrm{1}}}{R_{\mathrm{E}}\mathrm{\Delta }\mathit{\theta }}}\right)\text{(37)}& {\displaystyle }\begin{array}{rl}& \mathrm{\nabla }{\mathrm{\Sigma }}_{\mathit{\varphi },i+\frac{\mathrm{1}}{\mathrm{2}},j+\frac{\mathrm{1}}{\mathrm{2}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\\ & \cdot \left({\displaystyle \frac{{\mathrm{\Sigma }}_{i,j+\mathrm{1}}-{\mathrm{\Sigma }}_{i,j}}{R_{\mathrm{E}}\sin {\mathit{\theta }}_{\mathit{\varphi },i+\frac{\mathrm{1}}{\mathrm{2}},j+\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Delta }\mathit{\varphi }}}+{\displaystyle \frac{{\mathrm{\Sigma }}_{i+\mathrm{1},j+\mathrm{1}}-{\mathrm{\Sigma }}_{i+\mathrm{1},j}}{R_{\mathrm{E}}\sin {\mathit{\theta }}_{\mathit{\varphi },i+\frac{\mathrm{1}}{\mathrm{2}},j+\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Delta }\mathit{\varphi }}}\right)\end{array}\text{(38)}& {\displaystyle }{A}_{i+\frac{\mathrm{1}}{\mathrm{2}},j+\frac{\mathrm{1}}{\mathrm{2}}}=\mathrm{2}{R}_{\mathrm{E}}^{\mathrm{2}}\mathrm{\Delta }\mathit{\varphi }\sin {\mathit{\theta }}_{i+\frac{\mathrm{1}}{\mathrm{2}},j+\frac{\mathrm{1}}{\mathrm{2}}}\sin {\displaystyle \frac{\mathrm{\Delta }\mathit{\theta }}{\mathrm{2}}}\text{(39)}& {\displaystyle }{j}_{r,i+\frac{\mathrm{1}}{\mathrm{2}},j+\frac{\mathrm{1}}{\mathrm{2}}}=\mathrm{0}.\end{array}$$

Here, R_E is the Earth's radius, Δθ and Δϕ are the co-latitude and longitude resolution of the original conductivity grid, i is an index of the co-latitude grid, and j is an index of the longitude grid. As we do not have any information on the vertical current density, it is assumed to vanish. The SECS polar cap angle θ₀ was calculated from the grid cell area using Eq. (23). The final CF electric field was evaluated at the original conductance grid, i.e., at grid points (i,j).

Because of the computational limitations of a regular laptop, some optimization of the CF electric field calculation was needed to process the full IMAGE area: we split the conductivity model area into 240 pieces with an extent of 2° in latitude with ±2° padding and 1.9° in longitude with ±4° padding. The full area was assembled by fitting these 240 subareas together. This approximation did not extend to the DF electric field calculation, for which the full area was utilized. Naturally, the best approach would be to utilize more powerful computing resources so that no optimization would be needed at all. For our simple test the approximation suffices.

Figure 5 Divergence of the horizontal electric field (a), horizontal CF electric field amplitude (b), and horizontal total electric field amplitude (c) calculated using the CF SECS method with the DF electric field from Fig. 3f and vanishing vertical current density. (d) Divergence of the electric field calculated using the CF SECS method when the DF electric field and radial current density are from PGIEM2G (Fig. 6b, e). (e–f) CF electric field (e) and total electric field (f) approximated using Eqs. (25)–(26) with the coefficients from Fig. 7 and the DF electric field from Fig. 3f. For the electric field amplitude plots, the color bar minimum value and maximum value have been set to 10⁰ and 10⁴, respectively, following Marshalko et al. (2021). This places the white color of the diverging color map at 10².

Figure 5 displays the resulting charge density (Fig. 5a), amplitude of the horizontal part of the CF electric field E_CF (Fig. 5b), and total electric field E_DF+E_CF (Fig. 5c) for the same epoch as in Fig. 3. The edges of the computational subareas can be detected in some regions, indicating that larger areas would be better. As expected, charges accumulate at conductivity gradients, producing a total electric field that is significantly modulated compared to the DF electric field of Fig. 3f. The tendency for higher electric field in areas of lower conductivity and lower electric field in areas of higher conductivity is especially more pronounced.

3.3 Comparison with PGIEM2G

In this section, we compare the electric field modeled using the DF and CF SECS methods with the electric field modeled using PGIEM2G (Kruglyakov and Kuvshinov, 2018). Figure 6a shows the amplitude of the horizontal surface electric field modeled using PGIEM2G for the same epoch as in Fig. 5c. The rest of the panels show the DF part of the electric field and its curl (Fig. 6b), CF part of the electric field (Fig. 6c), divergence of the electric field (Fig. 6d), and divergence of the horizontal current density divided by conductance (Fig. 6e). We will come back to Fig. 6f later.

The 2-D SECS method can also be used to separate any vector field on a sphere into its DF and CF parts (Amm, 1997). We have utilized it to separate the total electric field produced by PGIEM2G. The DF part of the electric field was obtained using the modeled curl $(\mathrm{\nabla }\times \mathit{E}{)}_r=-\frac{\partial B_r}{\partial t}$ to estimate DF 2-D SECS amplitudes I_DF on the PGIEM2G grid:

$$\begin{array}{}\text{(40)}& {\displaystyle }(\mathrm{\nabla }\times \mathit{E}{)}_{r,i}={\displaystyle \frac{I_{\mathrm{DF},i}}{A_i}}-\sum _j{\displaystyle \frac{I_{\mathrm{DF},j}}{\mathrm{4}\mathit{\pi }{R}^{\mathrm{2}}}},\text{(41)}& {\displaystyle }{I}_{\mathrm{DF},i+\frac{\mathrm{1}}{\mathrm{2}},j+\frac{\mathrm{1}}{\mathrm{2}}}\approx A_{i+\frac{\mathrm{1}}{\mathrm{2}},j+\frac{\mathrm{1}}{\mathrm{2}}}(\mathrm{\nabla }\times \mathit{E}{)}_{r,i+\frac{\mathrm{1}}{\mathrm{2}},j+\frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

The CF part of the electric field was obtained by subtracting the DF part of the electric field from the total electric field. The divergence of the electric field (∇⋅E) was estimated using finite-difference calculation.

Comparison of the PGIEM2G (Fig. 6a) and SECS (Fig. 5c) results reveals that, unlike PGIEM2G, the SECS method does not produce a narrow band of intense CF electric field on the side of the low conductance at conductance gradients but enhances the electric field across the entire area of lower conductance. The SECS method assumes that the vertical current density is zero, but clearly this is not a valid assumption, as can be seen by examining the divergence of the horizontal current density in Fig. 6e. Comparison of Fig. 6e with the divergence of the electric field in Fig. 6d confirms that vertical currents are clearly responsible for producing the second layer of charges near conductivity gradients that largely confine the CF electric field near the gradients.

Another possible source of discrepancy between PGIEM2G and the SECS method is the internal part of the DF electric field. Whereas the SECS method reconstructs it from ground-based geomagnetic measurements, PGIEM2G models it based on the 3-D conductivity of SMAP. While sparse magnetometer coverage affects the DF electric field of the SECS method, inaccuracies in the conductivity model affect PGIEM2G. The difference can be seen by comparing the respective DF electric fields in Figs. 3f and Fig. 6b. PGIEM2G produces small-scale structures that cannot be resolved with the operational magnetometers. On the other hand, the curl of the electric field from the SECS method corresponds to the measured $-\partial B_r/\partial t$ at all IMAGE stations, whereas for PGIEM2G this is not true, as can be seen by comparing with $\partial B_r/\partial t$ in Fig. 3c.

Figure 6 Horizontal electric field amplitude modeled using PGIEM2G (Marshalko et al., 2021) (a), DF part of the electric field and its curl (arrows and color in panel b), CF part of the electric field (c), divergence of the electric field (d), and divergence of the horizontal current density divided by conductance (e). Panel (f) shows the divergence of the electric field derived using the SECS method when the DF electric field and radial current density are from PGIEM2G.

We test the above conclusions by running the SECS CF electric field calculation with PGIEM2G's DF electric field and vertical current density. The resulting divergence of the electric field is given in Fig. 5d for comparison with the CF SECS method and repeated in Fig. 6f for comparison with PGIEM2G. Now there is a much better correspondence between the CF SECS result in Fig. 6f and PGIEM2G result in Fig. 6d. Thus, we conclude that without knowledge of the vertical current density, the SECS method cannot reliably estimate the CF part of the electric field. The DF part of the electric field, however, can be considered reliable, at least in areas where the magnetometer coverage is good.

3.4 Coefficients relating E_CF to E_DF

As we cannot use the CF SECS method to estimate the CF part of the electric field, we try an alternative approach: the method described in Sect. 2.4, which assumes that the CF electric field is linearly proportional to the driving DF electric field. Figure 7 shows the coefficients derived using the CF and DF part of PGIEM2G's horizontal surface electric field between 7 September 2017 at 23:00 UT and 8 September 2017 at 01:00 UT. The standard deviation (SD) errors of the coefficients are displayed in Fig. 8. These are small compared to the values of the coefficients.

Figure 7 Proportionality coefficients (Eqs. 25–26) derived using the CF and DF part of PGIEM2G's horizontal surface electric field between 7 September 2017 at 23:00 UT and 8 September 2017 at 01:00 UT.

Figure 8 The same as Fig. 7 except that instead of the coefficients their standard deviation (SD) errors are shown. The color bar minimum value and maximum value have been set to 0 and 0.25, respectively, placing the white color of the diverging color map in the middle of these values.

In order to examine the validity of Eqs. (25)–(26), Fig. 9 shows a time series of the CF electric field at SOD between 7 September 2017 at 23:00 UT and 8 September 2017 at 01:00 UT. The original PGIEM2G-modeled value is shown in black and the approximation based on Eqs. (25)–(26) is shown in red. There is good correspondence between the original curve and approximation. Furthermore, Fig. 10 shows maps of the CF part of the horizontal electric field modeled by PGIEM2G on 7 September 2017 at 23:16:00 UT (Fig. 10a), the CF electric field approximated using Eqs. (25)–(26) with the coefficients from Fig. 7 and DF electric field from PGIEM2G as shown in Fig. 6b (Fig. 10b), and the difference between the original and approximated CF field (Fig. 10c). The differences are generally not very large, indicating that this approximation could be used to estimate the CF electric field from a known DF electric field. However, it is worth noting that the results are obtained for Fennoscandia and may be not valid in the case of different conductivity distributions.

Figure 9 CF electric field at SOD between 7 September 2017 at 23:00 UT and 8 September 2017 at 01:00 UT. The original PGIEM2G-modeled value is shown in black and an approximation based on Eqs. (25)–(26) is shown in red.

Figure 10 CF part of the horizontal electric field modeled by PGIEM2G on 7 September 2017 at 23:16:00 UT (a), CF electric field approximated using Eqs. (25)–(26) with the coefficients from Fig. 7 and DF electric field from PGIEM2G (Fig. 6b) (b), and the difference between the original and approximated CF field (c).

Figure 5e shows a map of the CF electric field approximated using Eqs. (25)–(26) with the coefficients from Fig. 7 and the IMAGE-based DF electric field from Fig. 3f. Figure 5f shows the corresponding total electric field, i.e., the sum of the DF electric field in Fig. 3f and the CF electric field in Fig. 5e. Now the CF and total electric field structures are similar to those produced by PGIEM2G (Fig. 6c), although there are some differences due to the DF electric field (Fig. 3f and Fig. 6b). Furthermore, Fig. 11 shows E_x and E_y at SOD between 7 September 2017 at 23:00 UT and 8 September 2017 at 01:00 UT. PGIEM2G modeling is shown in black, DF SECS modeling based on the IMAGE magnetic field with the CF part obtained using the coefficients from Fig. 7 in red, and DF SECS modeling based on the PGIEM2G magnetic field with the CF part obtained using the coefficients from Fig. 7 in green. We have also added estimates obtained via the MSTF (Kruglyakov et al., 2023) approach (blue curves). In Fig. 11, the SECS modeling based on IMAGE data shows peaks of higher amplitude than either PGIEM2G or the MSTF approach. Because the SECS modeling based on PGIEM2G magnetic field data has much better agreement with PGIEM2G, the peaks in the IMAGE-based SECS are clearly due to the difference in the internal DF electric field.

Figure 11 E_x and E_y at SOD between 7 September 2017 at 23:00 UT and 8 September 2017 at 01:00 UT. PGIEM2G modeling is shown in black and SECS (DF electric field from SECS, CF electric field from Eqs. 25–26) modeling in red (calculated using the magnetic field from IMAGE) and green (calculated using the magnetic field from PGIEM2G). The multi-site transfer function (MSTF) approach (Kruglyakov et al., 2023) is shown in blue.

Finally, Fig. 12 shows a time series of the electric field as a sum of the external DF, internal DF, and CF parts at the location of the magnetometer station SOD between 7 September 2017 at 23:00 UT and 8 September 2017 at 01:00 UT. The DF electric field has been modeled using the DF SECS method and IMAGE data, and the CF electric field has been obtained using the coefficients from Fig. 7. As can be seen from Fig. 4, SOD is located close to a conductivity gradient. The two top panels of Fig. 12 show the north (x) component of the electric field, the two middle panels show the east (y) component, and the two bottom panels show the amplitude. The external, internal, and total DF electric fields are shown in the top, third, and fifth panel and the total DF, CF, and electric field in the second, fourth, and sixth panel. Comparison of the external, internal, and total DF electric field curves reveals that the total DF electric field is sometimes dominated by the external part and sometimes by the internal part, and sometimes large, oppositely directed peaks in both parts cancel each other out, producing a small total DF electric field. Comparison of the DF and CF electric field curves show that in the x direction, the CF electric field is directed in the same way as the DF electric field and has a larger amplitude than the DF electric field, resulting in a total electric field x component which is clearly stronger than the DF electric field amplitude. In the y direction, the CF electric field also has a stronger amplitude than the DF electric field but the CF and DF parts are oppositely directed, resulting in a total electric field that is stronger and oppositely directed to the DF electric field.

Figure 12 Time series of the north (x) and east (y) components and amplitudes of the external (E_DF,ext) and internal (E_DF,int) part of the DF electric field, total (${\mathit{E}}_{\mathrm{DF}}={\mathit{E}}_{\mathrm{DF},\mathrm{ext}}+{\mathit{E}}_{\mathrm{DF},\mathrm{int}}$) DF electric field, CF electric field (E_CF), and total (E_DF+E_DF) electric field at the location of the magnetometer station SOD between 7 September 2017 at 23:00 UT and 8 September 2017 at 01:00 UT. The DF electric field has been modeled using the DF SECS method and IMAGE data, and the CF electric field has been obtained using the coefficients from Fig. 7. The dashed vertical lines indicate the epochs 23:16 and 23:52 UT displayed in Figs. 3, 5, A1, and A2.

3.5 Comparison with observations

3.5.1 GIC observations

Recordings of geomagnetically induced currents (GICs) in the Finnish natural gas pipeline were carried out close to the Mäntsälä (MAN) compressor station in southern Finland (60.6° N, 25.2° E) from November 1998 to October 2023. The measurements were based on the differential magnetometer method, utilizing two magnetometers: one right above the pipe and the other at the nearby IMAGE station NUR (60.50° N, 24.65° E) (Pulkkinen et al., 2001). The quality of the MAN magnetometer data is variable, but good data are available for the Halloween storm on 29–31 October 2003. GIC at MAN can also be estimated from the horizontal geoelectric field using the empirical expression (Pulkkinen et al., 2001; Viljanen et al., 2006)

$$\begin{array}{}\text{(42)}& \mathrm{GIC}\left(t\right)=-\mathrm{70}\mathrm{A}\mathrm{km}{\mathrm{V}}^{-\mathrm{1}}{E}_x\left(t\right)+\mathrm{88}\mathrm{A}\mathrm{km}{\mathrm{V}}^{-\mathrm{1}}{E}_y\left(t\right),\end{array}$$

which assumes a spatially uniform field in the region of the pipeline. In 2003, the magnetometer network around MAN was very sparse and the conductivity model in this region is known to be inaccurate (Marshalko et al., 2023), making direct comparison with the electric field modeled using the SECS method less than ideal. However, we can utilize the simplified version of Sect. 2.4.

For MAN, the contribution from the nearby magnetometer station NUR is likely to be much larger than that from the other, much more distant, stations. Combined with the linear approximation of Eq. (42), this allows us to estimate the GIC at MAN based on NUR $\frac{\partial \mathit{B}}{\partial t}$ alone:

$$\begin{array}{}\text{(43)}& \mathrm{GIC}\left(t\right)=a_x{\displaystyle \frac{\partial B_x}{\partial t}}+a_y{\displaystyle \frac{\partial B_y}{\partial t}}+a_z{\displaystyle \frac{\partial B_z}{\partial t}}.\end{array}$$

To test this approach, we have used the period from 29 October 2003 at 08:00:00 UT to 31 October 2003 at 23:59:50 UT to determine the coefficients. The resulting values are provided in Table 1. Observed and modeled GICs for the immediately preceding period on 29 October 2003 from 06:00:00 to 07:59:50 UT are shown in Fig. 13 and show good agreement in terms of time development and amplitude. The correlation coefficient is CC=0.80 and the coefficient of determination,

$$\begin{array}{}\text{(44)}& R^{\mathrm{2}}=\mathrm{1}-{\displaystyle \frac{{\sum }_i({\mathrm{GIC}}_i^{\mathrm{obs}}-{\mathrm{GIC}}_i^{\mathrm{mod}}{)}^{\mathrm{2}}}{{\sum }_i({\mathrm{GIC}}_i^{\mathrm{obs}}-\mathrm{mean}({\mathrm{GIC}}^{\mathrm{obs}}){)}^{\mathrm{2}}}},\end{array}$$

is R²=0.64.

Table 1 Coefficients for the linear relation of Eq. (43) between MAN GIC (in A) and NUR $\frac{\partial \mathit{B}}{\partial t}$ (in nT s⁻¹) obtained as a fit to the observed GIC from 29 October 2003 at 08:00:00 UT to 31 October 2003 at 23:59:50 UT. For comparison, the bottom set of coefficients is determined as a fit to the longer period from 29 October 2003 at 00:00:00 UT to 31 October 2003 at 23:59:50 UT.

Figure 13 GIC observed at MAN on 29 October 2003 from 06:00:00 UT to 07:59:50 UT (black) and GIC modeled based on NUR magnetic field observations using Eq. (43) and the coefficients provided in Table 1.

3.5.2 Electric field observations

Electric field observations suitable for comparison with modeling results are rare. However, there is an interval on 11 September 2005 from 05:15:00 to 06:15:00 UT with observations (Smirnov et al., 2006) from two sites near the IMAGE station MEK that has been used by Kruglyakov et al. (2023). Kruglyakov et al. (2023) called these sites M02 (63.043740° N, 30.657030° E) and M05 (62.938890° N, 30.993910° E). Similar to MAN, M02 and M05 are located in a region where the magnetometer coverage is far from ideal. Thus, we will again utilize Eq. (27) and approximate the electric field components by fits to $\frac{\partial \mathit{B}}{\partial t}$ from the nearby station MEK (62.77° N, 30.97° E) using Eqs. (28)–(29).

Because the electric field components provided by Kruglyakov and Marshalko (2023) are in geomagnetic coordinates and the magnetic field observation from IMAGE in geographic coordinates, the coefficients in Eqs. (28)–(29) also include rotation of the coordinate system. Geoelectric field observations contain effects due to charge build-up at local small-scale conductivity structures (see, e.g., Kruglyakov et al., 2023). These effects are also included in the coefficients. The coefficients resulting from fitting the period from 05:15:00–05:29:50 UT are provided in Table 2, and the observed and modeled electric field components are shown in Fig. 14. The period used for fitting is shaded. There is good agreement between the observed electric field components and the fits to MEK $\frac{\partial \mathit{B}}{\partial t}$ observations. For the period not used in the fitting, i.e., 05:30:00–06:15:00 UT, the correlation coefficients are 0.91 for M02 E_x, 0.95 for M02 E_y, 0.95 for M05 E_x, and 0.94 for M05 E_y. The coefficients of determination are 0.84 for M02 E_x, 0.89 for M02 E_y, 0.88 for M05 E_x, and 0.86 for M05 E_y. The good agreement between observations and modeling results here and in Sect. 3.5.1 supports the assumptions of linear dependence between the parameters.

Table 2 Coefficients for the linear relations of Eqs. (28)–(29) between M02 and M05 geomagnetic north (E_x) and south (E_y) components of the geoelectric field (in mV km⁻¹) and MEK geographic $\frac{\partial \mathit{B}}{\partial t}$ (in nT s⁻¹) obtained as a fit to the observed geoelectric field on 11 September 2005 from 05:15:00 to 05:29:50 UT.

Figure 14 Geomagnetic north (E_x) and east (E_y) components of the geoelectric field observed at the sites M02 (63.043740° N, 30.657030° E) and M05 (62.938890° N, 30.993910° E) on 11 September 2005 from 05:15:00 to 06:15:00 UT (black) and E_x and E_y modeled based on nearby MEK (62.77° N, 30.97° E) geographic magnetic field observations using Eqs. (28)–(29) and the coefficients provided in Table 2. The interval used to determine the coefficients is shaded.

4 Discussion

In the previous section we have shown that obtaining the DF part of the geoelectric field from ground-based magnetometer network data using the 2-D SECS method and estimating the CF part from the DF part using Eqs. (25)–(26) can be used to estimate the geoelectric field at the Earth's surface. In cases when the magnetometer network is too sparse for the 2-D SECS method, the same principle of linear relation can be used to estimate the geoelectric field from local $\partial \mathit{B}/\partial t$. In this section we will discuss further applications of the new method.

4.1 Maximum amplitudes of the event

The temporal variation in the DF and CF electric field behavior in Fig. 12 is caused by the combination of the spatiotemporally complex external source and the dynamic response of the fixed 3-D ground conductivity. Although induction in the 3-D ground is deterministic, it is unpredictable and difficult to understand, as mentioned in the Introduction. Ideally, we could determine the electric field everywhere in the Earth, calculate the telluric current density and charge density, and visualize them in 3-D to understand the complete induction process. However, this is not feasible with present computational resources. Our new method can be used to analyze the spatiotemporal variations of the different geoelectric field contributions: time variations of ionospheric currents, time variations of telluric currents, and accumulated charges in the ground. This can help clarify the complicated interaction between the ionosphere and the conducting ground.

As a demonstration, the top part of Table 3 provides the times, locations, and amplitudes of the maximum electric field contributions on the ground in the IMAGE area during a geomagnetic storm between 7 September 2017 at 23:00:00 UT and 8 September 2017 at 01:00:00 UT (Dimmock et al., 2019; Juusola et al., 2023a). The DF part of the electric has been calculated using the SECS method with IMAGE magnetic field as input, and the CF electric field has been estimated using Eqs. (25)–(26). The maximum of the external DF electric field took place on 7 September 2017 at 23:59:50 UT near stations ABK and AND. At this time, the total DF electric field also had a maximum in the same area, but a slightly larger value had already occurred earlier at 23:15:30 UT, near station KIR, due to a temporary suppression of the internal part of the DF electric field around that location (not shown). The largest internal DF electric field occurred near station DON on 8 September at 00:30:00 UT. The CF electric field maximum occurred at the same time as the maximum of the external DF field, on 7 September 2017 at 23:59:50 UT, at a location a little southeast of RVK. Due to the very large amplitude of the CF field, this was also the time and location of the total electric field maximum.

The total electric field maximum value exceeding 72 V km⁻¹ may sound high compared to the 20 V km⁻¹ value cited by Pulkkinen et al. (2012) as an extreme (one in 100 years) geoelectric field event or to the 8.5 V km⁻¹ value cited by Lanabere et al. (2024) as an extreme for Sweden. However, both Pulkkinen et al. (2012) and Lanabere et al. (2024) obtained their values using 1-D modeling, whereas our results are based on 3-D modeling. Because 1-D modeling of the electric field does not consider charges, it produces an estimate for the DF electric field. Our maximum DF electric field value of 3753 mV km⁻¹ is clearly below the extreme values cited by Pulkkinen et al. (2012) and Lanabere et al. (2024). In our modeling, the largest contribution to the maximum electric field comes from the CF electric field, highlighting the importance of 3-D effects when estimating extreme values. We are currently preparing a study on the estimation of the 3-D geoelectric field during extreme geomagnetic storms and another study on a detailed comparison of 3-D and 1-D geoelectric field modeling.

Table 3 Time, location, and amplitude of the maximum horizontal electric field contributions and time derivative of the horizontal magnetic field ($|\partial \mathit{H}/\partial t|$) contributions in the IMAGE area during a geomagnetic storm between 7 September 2017 at 23:00:00 UT and 8 September 2017 at 01:00:00 UT. The DF electric field maxima are given both on the ground (0 km altitude) and in the ionosphere (90 km altitude) and the CF and total electric field peaks as well as $|\partial \mathit{H}/\partial t|$ peaks only on the ground. The CF and total electric field maxima have been determined using the denser (0.03°×0.07°) PGIEM2G grid and the rest using the sparser (0.5°×1°) SECS grid.

The middle rows of Table 3 provide the times, locations, and amplitudes of the maximum time derivative of the horizontal magnetic field,

$$\begin{array}{}\text{(45)}& \left|{\displaystyle \frac{\partial \mathit{H}\left(t\right)}{\partial t}}\right|={\displaystyle \frac{\sqrt{\left[B_x\right(t)-B_x(t-\mathrm{\Delta }t){]}^{\mathrm{2}}+[B_y\left(t\right)-B_y(t-\mathrm{\Delta }t){]}^{\mathrm{2}}}}{\mathrm{\Delta }t}},\end{array}$$

on the ground, which is often used as a proxy for the geoelectric field (Viljanen, 1998; Viljanen et al., 2001). This proxy is based on an approximation of Faraday's law (Eq. 8) at the Earth's surface,

$$\begin{array}{}\text{(46)}& {\displaystyle \frac{\partial B_{\mathit{\theta }}}{\partial t}}={\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \left(r{E}_{\mathrm{DF},\mathit{\varphi }}\right)}{\partial r}}={\displaystyle \frac{E_{\mathrm{DF},\mathit{\varphi }}}{r}}+{\displaystyle \frac{\partial E_{\mathrm{DF},\mathit{\varphi }}}{\partial r}}\approx {\displaystyle \frac{E_{\mathrm{DF},\mathit{\varphi }}}{R_{\mathrm{E}}}},\text{(47)}& {\displaystyle \frac{\partial B_{\mathit{\varphi }}}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \left(r{E}_{\mathrm{DF},\mathit{\theta }}\right)}{\partial r}}=-{\displaystyle \frac{E_{\mathrm{DF},\mathit{\theta }}}{r}}-{\displaystyle \frac{\partial E_{\mathrm{DF},\mathit{\theta }}}{\partial r}}\approx -{\displaystyle \frac{E_{\mathrm{DF},\mathit{\theta }}}{R_{\mathrm{E}}}},\end{array}$$

where the first half of the equations is based on the results of Sect. 2.2 that at the Earth's surface $E_{\mathrm{DF},r}=\mathrm{0}$. The times and locations of the peak external, internal, or total DF electric field and $|\partial \mathit{H}/\partial t|$ values in Table 3 do not match. This demonstrates the complexity of geomagnetic induction as mentioned in Introduction. It also indicates that the drivers of the most intense geoelectric field peaks may not be exactly the same as the drivers of rapid geomagnetic variations (Juusola et al., 2023a).

4.2 Directionality of the geoelectric field

Because of the assumed linear dependence of the CF electric field on the DF electric field, it is possible to determine the direction of the driving DF field for which the electric field amplitude at a given location maximizes.

$$\begin{array}{}\text{(48)}& {\displaystyle }{E}_x=k_{xx}{E}_{\mathrm{DF},x}+k_{xy}{E}_{\mathrm{DF},y}+E_{\mathrm{DF},x}\text{(49)}& {\displaystyle }{E}_y=k_{yx}{E}_{\mathrm{DF},x}+k_{yy}{E}_{\mathrm{DF},y}+E_{\mathrm{DF},y}\text{(50)}& {\displaystyle }{E}_{\mathrm{DF},x}=\left|{\mathit{E}}_{\mathrm{DF}}\right|\cos \mathit{\alpha }\text{(51)}& {\displaystyle }{E}_{\mathrm{DF},y}=\left|{\mathit{E}}_{\mathrm{DF}}\right|\sin \mathit{\alpha }\text{(52)}& {\displaystyle }\left|\mathit{E}\right|=\left|{\mathit{E}}_{\mathrm{DF}}\right|\sqrt{\begin{array}{c}(k_{xx}\cos \mathit{\alpha }+k_{xy}\sin \mathit{\alpha }+\cos \mathit{\alpha }{)}^{\mathrm{2}}\\ +(k_{yx}\cos \mathit{\alpha }+k_{yy}\sin \mathit{\alpha }+\sin \mathit{\alpha }{)}^{\mathrm{2}}\end{array}}\end{array}$$

Figure 15 shows maps of the maximum (Fig. 15a) and minimum $\left|\mathit{E}\right|/\left|{\mathit{E}}_{\mathrm{DF}}\right|$ (Fig. 15b) and the angle α of maximum $\left|\mathit{E}\right|/\left|{\mathit{E}}_{\mathrm{DF}}\right|$ (Fig. 15c). The value of $\left|\mathit{E}\right|/\left|{\mathit{E}}_{\mathrm{DF}}\right|$ is the same for ±E_DF, which is why there is repetition in the color bar of the bottom panel. In some areas the maximum $\left|\mathit{E}\right|/\left|{\mathit{E}}_{\mathrm{DF}}\right|$ is greater than 1, indicating that the CF field enhances the DF field, and in some areas that ratio is smaller than 1, indicating that the CF field always weakens the DF field. The absolute maximum value of $\left|\mathit{E}\right|/\left|{\mathit{E}}_{\mathrm{DF}}\right|=\mathrm{40.66}$ in the area takes place at 64.48° and 11.77° longitude and occurs when the DF electric field angle is $\mathit{\alpha }=-\mathrm{52}\mathit{°}$; i.e., the field is directed perpendicular to the nearby coastline, where there are large conductance gradients. The condition that the derivative of Eq. (52) with respect to α is zero is repeated at intervals of $\mathit{\pi }/\mathrm{2}$, indicating that the minimum $\left|\mathit{E}\right|/\left|{\mathit{E}}_{\mathrm{DF}}\right|$ occurs when the DF field is perpendicular to the direction of the maximum enhancement. In most areas the minimum $\left|\mathit{E}\right|/\left|{\mathit{E}}_{\mathrm{DF}}\right|$ is close to 1; i.e., the CF part is close to zero and the DF part determines the electric field amplitude. Figure 15c emphasizes the active role of the passive Earth in generating geoelectric field peaks.

The above discussion on the dependence of the geoelectric field amplitude on its direction is related to the standard magnetotelluric concept of electric field polarization, typically described using polarization ellipses (Love et al., 2022). The importance of electric field polarization in GIC studies has been highlighted in several papers (Cordell et al., 2021; Love et al., 2019, 2022; Malone-Leigh et al., 2024; Murphy et al., 2021). Exploring the similarities and differences of the polarization as described by polarization ellipses and our method could be yet another interesting topic for a future study.

Figure 15 Maximum (a) and minimum $\left|\mathit{E}\right|/\left|{\mathit{E}}_{\mathrm{DF}}\right|$ (b) and the DF electric field direction where maximum $\left|\mathit{E}\right|/\left|{\mathit{E}}_{\mathrm{DF}}\right|$ is reached (c). For the panels in the top row, the color bar maximum value and minimum value correspond to the largest and smallest value of the shown parameter, respectively, placing the white color of the diverging color map in the middle of these values.

4.3 DF electric field in the ionosphere

The DF electric field in the ionosphere is typically ignored, although in some dynamical situations inductive effects are not negligible and the ionospheric electric field is not a pure CF field but has a significant DF part (Vanhamäki et al., 2007; Madelaire et al., 2024). The method used to estimate the DF electric field is valid at and above the Earth's surface. Hence, it can be used to estimate the DF electric field due to telluric and DF horizontal ionospheric currents in the ionosphere as well. The resulting DF electric field is not the total DF electric field above the ionospheric horizontal current sheet because the contribution from the CF horizontal ionospheric currents and field-aligned currents is missing. Below and at the ionospheric horizontal current sheet, this part of the DF electric field is zero because the combined magnetic field from the horizontal CF currents and (radial) field-aligned currents is zero (Fukushima, 1976; Amm, 1997). Moreover, Vanhamäki et al. (2007) argue that the DF electric field due to the CF horizontal and field-aligned currents should be very small.

Figure 16 shows the external (Fig. 16a), internal (Fig. 16b), and total (Fig. 16c) DF electric field in the ionosphere at 90 km altitude on 7 September 2017. Similar to the ground (Fig. 3), the external and internal parts of the DF electric field are more or less oppositely directed, but whereas on the ground the amplitudes of the external and internal parts were almost equal, in the ionosphere the external part has a clearly stronger amplitude than the internal part. The internal part is also clearly smoother in the ionosphere than on the ground. Nonetheless, the internal part significantly modifies both the pattern and amplitude of the resulting total DF electric field, indicating that ground induction should be included when ionospheric induction is considered. The maximum amplitudes of the total DF electric field, indicated below the scale arrows in Fig. 16, have similar values as the commonly observed ionospheric situations modeled by Vanhamäki et al. (2007).

Vanhamäki et al. (2005) used the 1-D complex image method (CIM) to estimate that the electric field caused by the Earth's induction is relatively small (at most 400 mV km⁻¹) and smooth at the ionospheric altitude. This is not in agreement with our result and indicates that the simple 1-D CIM modeling may not be sufficient for the task.

Figure 16 External (a), internal (b), and total (c) DF electric field (arrows) and its curl (color) in the ionosphere at 90 km altitude on 7 September 2017 at 23:16:00 UT.

The bottom part of Table 3 shows the times, locations, and amplitudes of the maximum ionospheric DF electric field contributions in the IMAGE area during a geomagnetic storm between 7 September 2017 at 23:00:00 UT and 8 September 2017 at 01:00:00 UT. The largest external and total DF electric fields both occurred at the same time and location, on 7 September 2017 at 23:15:30 UT. This was the same time when the DF total electric field maximum occurred on the ground as well. The internal DF electric field had its maximum almost at the same time both in the ionosphere and on the ground.

The method for deriving the CF electric field from known conductance and vertical current distributions derived in Sect. 2.3 can also be applied to the ionosphere. In that case, the ground conductance Σ is replaced by Pedersen conductance Σ_P and terms for the Hall conductance need to be added:

$$\begin{array}{}\text{(53)}& {\displaystyle }\mathit{J}={\mathrm{\Sigma }}_P\mathit{E}+{\mathrm{\Sigma }}_H\mathit{E}\times {\widehat{\mathit{e}}}_{\mathit{r}},\text{(54)}& {\displaystyle }\begin{array}{rl}\mathrm{\nabla }\cdot \mathit{J}& =\mathrm{\nabla }{\mathrm{\Sigma }}_P\cdot ({\mathit{E}}_{\mathrm{CF}}+{\mathit{E}}_{\mathrm{DF}})+{\mathrm{\Sigma }}_P\mathrm{\nabla }\cdot {\mathit{E}}_{\mathrm{CF}}+\mathrm{\nabla }{\mathrm{\Sigma }}_H\\ & \cdot \left[\right({\mathit{E}}_{\mathrm{CF}}+{\mathit{E}}_{\mathrm{DF}})\times {\widehat{\mathit{e}}}_{\mathit{r}}]+{\mathrm{\Sigma }}_H(\mathrm{\nabla }\times {\mathit{E}}_{\mathrm{DF}}{)}_r,\end{array}\text{(55)}& {\displaystyle }\begin{array}{rl}(\mathrm{\nabla }\times \mathit{J}{)}_r& =[\mathrm{\nabla }{\mathrm{\Sigma }}_P\times ({\mathit{E}}_{\mathrm{CF}}+{\mathit{E}}_{\mathrm{DF}}){]}_r+{\mathrm{\Sigma }}_P(\mathrm{\nabla }\times {\mathit{E}}_{\mathrm{DF}}{)}_r\\ & +[\mathrm{\nabla }{\mathrm{\Sigma }}_H\times (({\mathit{E}}_{\mathrm{CF}}+{\mathit{E}}_{\mathrm{DF}})\times {\widehat{\mathit{e}}}_{\mathit{r}}){]}_r\\ & -{\mathrm{\Sigma }}_H(\mathrm{\nabla }\cdot {\mathit{E}}_{\mathrm{CF}}).\end{array}\end{array}$$

Equation (53) is Ohm's law in a thin-sheet ionosphere, Eq. (54) is the ionospheric equivalent of Eq. (20), and Eq. (55) is a similar expression for the curl of the horizontal current. If E_DF is estimated to be insignificant compared to E_CF, Eq. (54) becomes the traditional problem of the electrostatic ionosphere of global magnetosphere–ionosphere simulations, where E_CF is solved from known conductance and vertical current distributions. Although E_DF is typically ignored, it would be possible to include it in the calculation by first using Ohm's law to derive an expression for I_DF(t) as a function of the SECS amplitudes Q_CF(t), I_CF(t), and I_DF(t−Δt) and inserting this into Eq. (54) to obtain an equation for Q_CF(t) in terms of the known I_CF(t) and I_DF(t−Δt). This approach resembles the inductive ionosphere solver presented by Vanhamäki (2011). Ground-based magnetometer networks can be used to determine (∇×J)_r and E_DF, but obtaining Σ_P and Σ_H is challenging. If they could be estimated, for example from all-sky camera images, Eq. (55) would yield E_CF. This approach otherwise resembles the corresponding solver suggested by Vanhamäki (2011), except that E_DF would be directly obtained from ground-based magnetometer data, and thus no time integration would be needed.

A thorough study on the significance of the internal part of the DF electric field on the total ionospheric DF electric field and on the significance of the ionospheric DF electric field to the ionospheric total electric field is a topic for further study. The results may be of interest for global simulations, which typically ignore induction (e.g., Mukhopadhyay et al., 2022; Sorathia et al., 2023; Ganse et al., 2025).

5 Conclusions

We have developed a new method for estimating various contributions to the 3-D geoelectric field at the Earth's surface. The surface electric field consists of an external DF electric field due to time-varying ionospheric and magnetospheric currents, an internal DF electric field due to time-varying telluric currents, and a CF electric field due to charge accumulation at ground conductivity gradients.

1. The external part of the DF electric field is calculated from the time derivative of the external part of the observed ground magnetic field using DF 2-D SECSs.

2. The internal part of the DF electric field is calculated from the time derivative of internal part of the observed ground magnetic field using DF 2-D SECSs.

3. The surface CF electric field is calculated from the known surface DF electric field and coefficients that linearly relate the DF electric field to the CF electric field. The coefficients were obtained from the PGIEM2G model.

4. The calculations are carried out in the time domain, and only two consecutive time steps of the observed magnetic field are needed to compute the surface electric field. This makes the method ideal for near-real-time applications. The external part of the DF electric field is valid everywhere, the internal part at and above the Earth's surface, and the CF part at the Earth's surface. A dense magnetometer network is required for good results.

5. The external and internal parts of the DF electric field are generally oppositely directed and have comparable amplitudes both on the ground and in the ionosphere, indicating that both contributions are significant for the total DF electric field. The largest peaks in the total DF electric field tend to occur when either the external or internal contribution is temporarily suppressed at the location of interest.

6. At a given location, a DF electric field with a given amplitude can result in a total surface electric field amplitude with a difference of orders of magnitude depending on the direction of the DF electric field with respect to the locally dominant conductivity gradient structure.

7. Peak amplitudes of the various electric field contributions did not occur at the same time or at the same location as the peak amplitudes of the time derivative of the horizontal magnetic field for our example event. This indicates that analysis of rapid magnetic field variations may not describe all relevant aspects of the electric field behavior.

8. The linear dependence of the DF electric field on $\partial \mathit{B}/\partial t$ observations on the one hand, and of the CF electric field on the DF electric field on the other hand, makes it possible to estimate the total geoelectric field directly from nearby magnetic field observations. As an example, we have determined coefficients that relate the geoelectric-field-driven GIC at MAN to NUR $\partial \mathit{B}/\partial t$.

Analyzing the separated contributions from currents (DF electric field) and charges (CF electric field) to the geoelectric field can help clarify the complicated interaction between the ionosphere and the conducting ground. As we have shown, the DF field is generally spatially clearly smoother than the CF field. A significant point is that a 3-D ground conductivity leads to charge accumulation at conductivity gradients and to large localized increases in the CF field. Consequently, although the ground conductivity does not depend on time, the response of the ground to external driving is highly dynamical and can produce complicated patterns of the geoelectric field.

Appendix A: Figures for the epoch on 7 September 2017 at 23:52:00 UT

In this Appendix we show Figs. A1, A2, A3, A4, and A5. They correspond to Figs. 3, 5, 6, 10, and 16, respectively, but illustrate the epoch on 7 September 2017 at 23:52:00 UT instead of the epoch on 7 September 2017 at 23:16:00 UT.

Figure A1 The same as Fig. 3 except for 23:52:00 UT instead of 23:16:00 UT.

Figure A2 The same as Fig. 5 except for 23:52:00 UT instead of 23:16:00 UT.

Figure A3 The same as Fig. 6 except for 23:52:00 UT instead of 23:16:00 UT.

Figure A4 The same as Fig. 10 except for 23:52:00 UT instead of 23:16:00 UT.

Figure A5 The same as Fig. 16 except for 23:52:00 UT instead of 23:16:00 UT.