A decade of the fast-varying ionospheric and magnetospheric magnetic fields from ground and multisatellite observations

SUMMARY

The time-varying geomagnetic field is a superposition of contributions from multiple internal and external current systems. A major source of geomagnetic variations at periods less than a few years are current systems external to the solid Earth, namely the ionospheric and magnetospheric currents, as well as associated induced currents. The separation of these three sources is mathematically underdetermined using either ground or satellite measurements alone, but becomes tractable when the two data sets are combined. Based on this concept, we developed a new geomagnetic field modelling approach that allows us to simultaneously characterize the mid-latitude ionospheric, magnetospheric and the internal induced magnetic fields using ground and satellite observations for all local times and magnetic conditions, and without prescribing any harmonic behaviour on these current systems in time, as is typical in other models. By applying this new method to a 10-yr data set of ground observatory and multisatellite measurements from 2014 to 2023, we obtained the time-series of the spherical harmonic coefficients of the ionospheric, magnetospheric and induced fields. These new time-series allow the study of complex non-periodic dynamics of the external magnetic fields during global geomagnetic storms, as well as periodicities in the magnetospheric coefficients linked to solar activities and periodic ionospheric magnetic fields linked to lunar daily variations, contributing to a more complete picture of the dynamics of the external currents and magnetosphere–ionosphere interactions, and facilitating more accurate space weather nowcast and forecast. Finally, the new approach allows for a better characterization of internal induced field sources, leading to higher quality electromagnetic transfer functions.

Geomagnetic induction, Magnetic field variations through time, Satellite magnetics, Fourier analysis

1 INTRODUCTION

Time variations of the magnetic field observed on the ground or at a spacecraft represent a superposition of contributions from multiple time-varying current systems. The electric current generated by the dynamo action in the fluid outer core is responsible for sustaining the main field over geological timescales, and is also the major contributor to magnetic field variations at long periods, known as secular variations (Backus et al. 1996). Meanwhile, the extraneous current systems in the ionosphere and the magnetosphere generate time-varying magnetic fields, which then induce currents in the electrically conductive solid Earth (Schuster 1889; Olsen 1999a). Characterizing the ionospheric and magnetospheric contributions to the observed magnetic field is of interest for several reasons. First, magnetic fields of ionospheric and magnetospheric origins are used to quantitatively study respective electric current systems, and serve as a proxy for understanding the underlying dynamics in these regions (Yamazaki & Maute 2017; Tsyganenko 2019; Laundal et al. 2021). On the other hand, the relation between these magnetic fields of external origin and their induced counterparts in the Earth interior can be used to probe the electrical conductivity in the Earth’s subsurface (Olsen 1999b; Kelbert et al. 2009; Kuvshinov 2012; Grayver 2024). More accurate space physics and subsurface conductivity models are then essential for a practically useful space weather hazard evaluation, as the latter is driven by the dynamics in the extraneous current systems (Tóth et al. 2005; Pulkkinen et al. 2013) but sensitive to the subsurface electrical conductivity as well (Juusola et al. 2020; Kelbert & Lucas 2020). Last but not least, better characterization of the ionospheric and magnetospheric systems also improves the core field modelling (Finlay et al. 2017), and indirectly leads to improvements in crust and ocean field models (Olsen et al. 2017; Grayver et al. 2024).

Ionospheric and magnetospheric contributions in the geomagnetic field models have long been constrained by using the magnetic field measurements from ground observatories (Sugiura 1964; Pulkkinen et al. 2003; Olsen et al. 2014). While this conventional source of geomagnetic observations provides high-quality validated long time-series (Love & Chulliat 2013; Macmillan & Olsen 2013), there are some shortcomings. First, their distribution over the globe is very non-uniform and the number of observatories does not grow. This prohibits the accurate retrieval of geomagnetic field features at higher spherical harmonic degrees, and poses a limit to the spatial resolution of the geomagnetic field models, especially over the oceans. Another fundamental limitation is related to separation of the external field sources. Since all magnetic field contributions of external origin are external to the ground observatories, it is mathematically impossible to separate the ionospheric signal from the magnetospheric signal, unless restrictive prior constraints are imposed. Previous attempts to study the individual extraneous magnetic signals mostly rely on temporal intervals, frequency bands or spatial modes where the contribution of one of the current systems is assumed to be dominant a priori. For instance, the magnetospheric contribution is often studied by working with the Dst index (Sugiura 1964) or its analogues (e.g. the RC-index; Finlay et al. 2015), which is assumed to reflect the axisymmetric part of the magnetospheric ring current (Maus & Weidelt 2004; Olsen et al. 2005), although it most certainly contains part of the ionospheric field as well. The large-scale ionospheric Sq fields were studied by selecting the data during magnetic quiet times and parametrizing the temporal evolution using a few space and time harmonics (Winch 1981; Schmucker 1999; Yamazaki & Maute 2017). In studies by Egbert et al. (2021) and Zenhäusern et al. (2021), authors adopted a more advanced spatial basis derived from a physics-driven model of the ionospheric dynamics called TIE-GCM. While these models provide a more detailed spatial description of magnetic fields, and allow one to perform physics-based interpolation in regions with lack of observations, their temporal content is limited and mostly focuses on modelling dynamics within the daily band. In addition, it is not trivial to extend these models to include largely non-periodic magnetospheric fields. Further, more work is needed to be able to reconcile this approach with satellite observations.

Noting the problems and limitations with geomagnetic ground observatories, recent studies increasingly incorporate satellite measurements in geomagnetic field modelling. Continuous magnetic field observations from Low-Earth-Orbit (LEO) satellites have enabled significant advances in geomagnetic field modelling, leading to improved spatial resolution of various magnetic field models, including those for the slow-varying or static core (Alken et al. 2021), crustal (Maus et al. 2007; Olsen et al. 2017) and oceanic (Sabaka et al. 2020; Grayver et al. 2024) fields. This also holds for models of ionospheric and magnetospheric fields (Lesur et al. 2010; Chulliat et al. 2016; Sabaka et al. 2018; Baerenzung et al. 2020; Finlay et al. 2020). However, the second problem regarding separation of external field sources is still mostly at large. Most external field models attempt to achieve this by reducing the spatial and/or temporal complexity of the external field model parametrizations. For instance, by taking only the dark-side magnetic observations, these models attribute the magnetic field external to observatories and satellites to the magnetosphere, assuming ionospheric magnetic fields can be neglected at night (Lesur et al. 2010; Baerenzung et al. 2020; Finlay et al. 2020). This approach limits the usable data to only the dark-side field measurements, typically discarding more than half of the available data. In addition, even the dark-side data may contain ionospheric contributions. Apart from the potential dynamics in the F region during the night, the transient behaviour of electromagnetic induction also means that the dark-side data contain internal signals that are induced by the day-time ionosphere and persist into the night (Maus & Weidelt 2004; Grayver et al. 2021). Another common practice is to use only observations taken during quiet magnetic conditions. For instance, a dedicated magnetic field model of the ionosphere and equatorial electrojet (EEJ) by Chulliat et al. (2016) delivers a model of the mid-latitude ionosphere during quiet phases. Although the model has a high spatial resolution, the temporal dynamics of the ionosphere is approximated by a small set of daily and seasonal time harmonics. Other temporal components are approximated by modulating the predicted field with an external time-varying proxy (e.g. F10.7 flux). More elaborate comprehensive inversion (CI) approach tries to parametrize the ionospheric magnetic field together with other components, including a lowest degree magnetospheric field (Sabaka et al. 2018). While the parametrization of ionospheric currents in CI is quite elaborate and fits observations at mid and equatorial latitudes well, a small set of discrete points in the frequency domain at annual, semiannual and diurnal bands also means that other temporal variations are not captured. Similarly, the model is constructed using only quiet magnetic periods. Further, these models are difficult to update on the fly as small volumes of new observations are added in time, limiting their applications for operational space-weather or for studying individual events, such as geomagnetic storms.

To overcome the limitations of current models, and to better exploit the increasing number of satellite observations delivered by dedicated missions and platform magnetometers, we present here a novel approach to simultaneously characterize the ionospheric, magnetospheric and internally induced magnetic fields using all available ground and space observations and high temporal resolution. Our approach utilizes the geometrical configuration of ground and low Earth orbit (LEO) geomagnetic observations. Indeed, the ground observatories alone cannot distinguish between the ionospheric and magnetospheric contributions, as both are external to the observations. The introduction of satellite observations in the region between the ionosphere and magnetospheric currents, however, renders the problem feasible. We construct model in short time-bins, allowing for a continuous Swarm-era model, including all local times and magnetic activity conditions. In addition to ionospheric and magnetospheric components, we also estimate time-series of the spherical harmonic coefficients of the internally induced field, rendering the model output also relevant for mantle induction studies.

The rest of the article is organized as follows. We introduce the model setup in Section 2, focusing on how our formulation would allow simultaneous co-estimation of external sources using short time bins instead of pre-defined time harmonics. This is followed by a description of the data pre-processing (Section 3) and model parametrization (Section 4). In Section 5 we present the results from processing the 10-yr time-series with our designated methodology. The results provide some new inputs for electromagnetic induction studies, and insights into ionospheric and magnetospheric dynamics, which we will discuss in Section 6.

2 METHODOLOGY

The time-dependent magnetic field at location |$\mathbf {r}$| and time t is described by the Maxwell’s equations, which read

$$\begin{eqnarray} \begin{aligned} \nabla \cdot \mathbf {B}(\mathbf {r}, t) &= 0, \\ \nabla \times \mathbf {B}(\mathbf {r}, t) &= \mu \left(\mathbf {j}(\mathbf {r}, t) + \varepsilon \frac{\partial \mathbf {E}(\mathbf {r}, t)}{\partial t}\right), \end{aligned} \end{eqnarray}$$

(1)

where |$\mathbf {j}$| is the total current density, |$\mu$| is the magnetic permeability, |$\varepsilon$| is the electric permittivity and |$\mathbf {B}$| and |$\mathbf {E}$| are the magnetic and electric fields, respectively. Unless otherwise stated, the physical variables are given in SI units. In an electrically insulating medium under magneto-quasi-static approximation (MQS) (Larsson 2007), the right hand side of the second equation vanishes identically, yielding an irrotational magnetic field, in which case the magnetic field is a potential field, expressed through a scalar potential V,

$$\begin{eqnarray} \mathbf {B}(\mathbf {r}, t) = - \nabla V(\mathbf {r}, t). \end{eqnarray}$$

(2)

Combined with the solenoidal property of the magnetic field, the scalar potential is shown to satisfy the Laplace equation

$$\begin{eqnarray} \nabla ^2 V = 0. \end{eqnarray}$$

(3)

In a spherical shell, this leads to the following form of general solutions

$$\begin{eqnarray} V(\mathbf {r}, t) &=& a \sum _{n=1}^{+\infty } \sum _{m=0}^n \left [\left(g_{nm}(t) \cos m\phi + h_{nm}(t) \sin m\phi \right) \left(\frac{r}{a}\right)^{-(n+1)}\right. \\ &&\left. + \left(q_{nm}(t) \cos m\phi + s_{nm}(t) \sin m\phi \right) \left(\frac{r}{a}\right)^n \right ] P_n^m(\cos \theta ), \end{eqnarray}$$

(4)

where |$(r,\theta ,\phi )$| are the spherical radius, the colatitude and the azimuth, |$P_n^m(\cos \theta )$| is the associated Legendre polynomial of degree n and order m, and a is a scaling factor for radius. The time-dependent coefficients |$g_{nm}$| and |$h_{nm}$| represent the contribution to the magnetic field of internal origin (i.e. from electric current sources within the inner boundary of the shell), and |$q_{nm}$| and |$s_{nm}$| represent the contributions of external origin (i.e. from electric current sources outside the outer boundary of the shell).

In geomagnetism, the spherical harmonic (SH) representation (4) applies to any observations made in the electric current free region surrounding the Earth. In this case, the radius scaling factor a is conveniently chosen to be the mean radius of the Earth, and |$P_n^m(\cos \theta )$| is chosen to be Schmidt quasi-normalized associated Legendre polynomial (Winch et al. 2005). Coefficients |$g_n^m$|⁠, |$h_n^m$|⁠, |$q_n^m$| and |$s_n^m$| are referred to as the real Gauss coefficients.

The SH representation has long been used for constructing geomagnetic field models from ground-based observations (Gauss 1839). Since the ground observations are made in the electrically insulating neutral atmosphere, the potential formulation of the magnetic field in a spherical shell is valid. In recent decades, an increasing number of studies also included satellite observations in the same SH representation (Langel et al. 1980; Bloxham & Jackson 1989; Langlais & Mandea 2000; Olsen et al. 2006; Baerenzung et al. 2020). This approach is justified since satellite observations taken at 400–800 km altitude are away from major source regions: the magnetospheric current systems (Ganushkina et al. 2018), and the E-region dynamo in the ionosphere (Yamazaki & Maute 2017).

Note, however, because of the presence of ionospheric currents, for instance those generated in the E-region dynamo, the ground observatories and the satellites do not ‘see’ the same set of source configurations. Taking this into account, we introduce two sets of Gauss coefficients. For the ground observations, the SH representation reads

$$\begin{eqnarray} V_\mathrm{obs}(\mathbf {r}, t) &=& a \sum _{n=1}^{+\infty } \sum _{m=0}^n \left [\left(g_{nm}^\mathrm{obs}(t) \cos m\phi + h_{nm}^\mathrm{obs}(t) \sin m\phi \right) \left(\frac{r}{a}\right)^{-(n+1)} \right. \\ && \left.+ \left(q_{nm}^\mathrm{obs} (t)\cos m\phi + s_{nm}^\mathrm{obs}(t) \sin m\phi \right) \left(\frac{r}{a}\right)^n \right ] P_n^m(\cos \theta ), \\ \end{eqnarray}$$

(5)

while for the satellite observations we have

$$\begin{eqnarray} V_\mathrm{sat}(\mathbf {r}, t) &=& a \sum _{n=1}^{+\infty } \sum _{m=0}^n \left [\left(g_{nm}^\mathrm{sat}(t) \cos m\phi + h_{nm}^\mathrm{sat}(t) \sin m\phi \right) \left(\frac{r}{a}\right)^{-(n+1)}\right. \\ &&\left. + \left(q_{nm}^\mathrm{sat}(t) \cos m\phi + s_{nm}^\mathrm{sat}(t) \sin m\phi \right) \left(\frac{r}{a}\right)^n \right ] P_n^m(\cos \theta ). \\ \end{eqnarray}$$

(6)

In this study, we are interested in the magnetic fields that are ultimately attributed to extraneous current systems. For this purpose, we consider three major electrical current systems contributing to the magnetic field: (1) the magnetospheric currents, the most dominant of which is the ring current, (2) the ionospheric currents, mostly from the E-region dynamo and (3) the corresponding induced signals originating in the solid Earth or in the ocean, and hence internal to all observations. As we are most interested in the large-scale features in the magnetic field, we shall restrict the ionospheric currents to the mid-latitude currents. This excludes the equatorial and polar current systems, which require dedicated parametrizations. Meanwhile, the core field with its secular variations and the lithospheric field are considered removed from the system. The geometric setup of our model is illustrated in Fig. 1.

Figure 1.

Geometric configuration of the model, showing time-dependent spatially varying external and internal electric currents at middle geomagnetic latitudes (i.e. excluding the polar regions).

Open in new tab Download slide

Under these considerations, the observatory Gauss coefficients of internal origin represent contributions from electric currents induced in the Earth’s interior or in the ocean by nature of electromagnetic induction. We write these coefficients alternatively as

$$\begin{eqnarray} g_{nm}^\mathrm{obs} = g_{nm}^\mathrm{int}, \quad h_{nm}^\mathrm{obs} = h_{nm}^\mathrm{int}. \end{eqnarray}$$

(7)

The satellite Gauss coefficients of external origin represent contributions from magnetospheric currents, and are written as

$$\begin{eqnarray} q_{nm}^\mathrm{sat} = q_{nm}^\mathrm{mag}, \quad s_{nm}^\mathrm{sat} = s_{nm}^\mathrm{mag}. \end{eqnarray}$$

(8)

The observatory Gauss coefficients of external origin contain contributions from the magnetospheric currents as well as from the ionospheric currents. They can thus be decomposed as

$$\begin{eqnarray} q_{nm}^\mathrm{obs} = q_{nm}^\mathrm{ion} + q_{nm}^\mathrm{mag}, \quad s_{nm}^\mathrm{obs} = s_{nm}^\mathrm{ion} + s_{nm}^\mathrm{mag}. \end{eqnarray}$$

(9)

The same holds for the satellite Gauss coefficients of internal origin. They describe the field generated by both ionospheric currents and induced currents, and can be decomposed as

$$\begin{eqnarray} g_{nm}^\mathrm{sat} = g_{nm}^\mathrm{ion} + g_{nm}^\mathrm{int}, \quad h_{nm}^\mathrm{sat} = h_{nm}^\mathrm{ion} + h_{nm}^\mathrm{int}. \end{eqnarray}$$

(10)

Note that the coefficients |$g_{nm}^\mathrm{int}$| and |$h_{nm}^\mathrm{int}$| in eq. (10) are the same as in (7). Since the induced currents reside within the Earth, and are internal to both the ground and the satellite observations, the field generated at the surface can be upward continued as the same potential field to a satellite altitude. Similarly, the coefficients |$q_{nm}^\mathrm{mag}$| and |$h_{nm}^\mathrm{mag}$| in eqs (9) and (8) are also the same. The only component different is the ionospheric component, described by |$q_{nm}^\mathrm{ion}$| and |$s_{nm}^\mathrm{ion}$| in the ground setting, and |$g_{nm}^\mathrm{ion}$| and |$h_{nm}^\mathrm{ion}$| in the satellite setting. It is generally not possible to establish a one-to-one relation between the fields generated by the ionosphere within the ionosphere and outside the ionosphere, without knowing a priori the spatial distribution of the ionospheric currents. However, assuming that the electric current in the ionosphere occupies a domain with limited radial extent, one can adopt the thin-sheet approximation, whereby the ionospheric current is approximated as a spherical sheet of radius |$a+h$|⁠. This approximation does not only associate a unique sheet current to an internal or an external magnetic field, but also provides the continuity of the radial magnetic field across the thin sheet:

$$\begin{eqnarray} B_r|_{r\rightarrow (a+h)^+} = B_r|_{r\rightarrow (a+h)^-}. \end{eqnarray}$$

(11)

Using |$B_r = -\partial _r V_\mathrm{obs}$| for |$r < a+h$| and |$B_r = -\partial _r V_\mathrm{sat}$| for |$r>a+h$|⁠, a one-to-one relation between the ionospheric Gauss coefficients can be established (see also Sabaka et al. 2002)

$$\begin{eqnarray} \begin{aligned} n\left(\frac{a+h}{a}\right)^{n-1} q_{nm}^\mathrm{ion} &= - (n+1)\left(\frac{a+h}{a}\right)^{-(n+2)} g_{nm}^\mathrm{ion},\\ n\left(\frac{a+h}{a}\right)^{n-1} s_{nm}^\mathrm{ion} &= - (n+1)\left(\frac{a+h}{a}\right)^{-(n+2)} h_{nm}^\mathrm{ion}. \end{aligned} \end{eqnarray}$$

(12)

Combining eqs (7)–(10) with (12), we can rewrite the SH representation for ground observations (5) as

$$\begin{eqnarray} &V_\mathrm{obs}(\mathbf {r},t) = a \sum _{n=1}^{+\infty } \sum _{m=0}^n \Big [\left(g_{nm}^\mathrm{int}(t) \cos m\phi + h_{nm}^\mathrm{int}(t) \sin m\phi \right) \left(\frac{r}{a}\right)^{-(n+1)} \\ &+ \left(\left(q_{nm}^\mathrm{ion}(t) + q_{nm}^\mathrm{mag}(t)\right) \cos m\phi + \left(s_{nm}^\mathrm{ion}(t) + s_{nm}^\mathrm{mag}(t)\right) \sin m\phi \right) \left(\frac{r}{a}\right)^n \Big ] P_n^m(\cos \theta ), \end{eqnarray}$$

(13)

and the SH representation for satellites (6) as

$$\begin{eqnarray} \begin{aligned} V_\mathrm{sat}(\mathbf {r}, t) = a &\sum _{n=1}^{+\infty } \sum _{m=0}^n \Bigg [\Bigg (\left(g_{nm}^\mathrm{int}(t) - \frac{n}{n+1} \left(\frac{a+h}{a}\right)^{2n+1} q_{nm}^\mathrm{ion}(t)\right) \cos m\phi \\ + &\left(h_{nm}^\mathrm{int}(t) - \frac{n}{n+1} \left(\frac{a+h}{a}\right)^{2n+1} s_{nm}^\mathrm{ion}(t) \right) \sin m\phi \Bigg ) \left(\frac{r}{a}\right)^{-(n+1)} \\ + &\left(q_{nm}^\mathrm{mag}(t) \cos m\phi + s_{nm}^\mathrm{mag}(t) \sin m\phi \right) \left(\frac{r}{a}\right)^n \Bigg ] P_n^m(\cos \theta ). \end{aligned} \end{eqnarray}$$

(14)

These representations reduce the magnetic potential to the source-specific Gauss coefficients |$(g_{nm}^\mathrm{int}, h_{nm}^\mathrm{int})$|⁠, |$(q_{nm}^\mathrm{ion}, s_{nm}^\mathrm{ion})$|⁠, |$(q_{nm}^\mathrm{mag}, s _{nm}^\mathrm{mag})$|⁠, hereinafter referred to as internally induced coefficients, ionospheric coefficients and magnetospheric coefficients, respectively. These Gauss coefficients can be estimated by minimizing the data misfit of the combined satellite and ground observatory data set. The compilation and pre-processing of the data set and the model estimation process are detailed in Sections 3 and 4, respectively.

3 DATA

We compiled a joint data set of ground observatory and satellite observations from data provided by the European Space Agency’s (ESA) VirES for Swarm service.¹ The data set spans a total length of nearly 10 yr, from January 2014 to December 2023. The compilation and pre-processing of the data is largely facilitated by the Python interface viresclient (Smith et al. 2024).

While eqs (13) and (14) provide magnetic field models that are continuous in both space and time, we implement a time-discrete version by estimating the Gauss coefficients within non-overlapping time bins. For this purpose, the ground and satellite observations must be both arranged into separate time bins, within which the corresponding SH coefficients are assumed to be constant. This inevitably introduces a trade-off between spatial coverage and the temporal resolution of the satellite data. We found that a time bin of 3 hr provides a fair balance between the spatial and temporal resolution. The model parametrization can be easily adjusted depending on future applications and availability of new satellite missions and observatories.

3.1 Satellite data

We used vector magnetic field measurements from the Swarm satellites (Swarm data product SW_MAGx_LR_1B where |$\tt {x} \in \lbrace \tt {A,B,C}\rbrace$|⁠) within the designated time period (Olsen et al. 2013). For better local time coverage, the Swarm data are further supplemented by the vector magnetic field measurements from CryoSat-2 (Swarm data product CS_OPER_MAG) from 2014 to 2023 (Olsen et al. 2020), as well as those from Grace-FO satellites (Swarm data product GFx_OPER_FGM_ACAL_CORR where |$\tt {x} \in \lbrace \tt {1,2}\rbrace$|⁠) from 2018 to 2023 (Stolle et al. 2021). All satellite measurements are downsampled to a one-minute sampling interval. As we are interested in the magnetic fields that are ultimately attributed to the extraneous sources, the core field (product SW_MCO_SHA_2C) and the static lithospheric field (product SW_MLI_SHA_2C) as given by the comprehensive inversion chain (Sabaka et al. 2020) are subtracted from all magnetic field measurements.

We then filter out satellite observations taken at locations with an absolute quasi-dipole (QD) latitude (Richmond 1995) below |$5^\circ$| or above |$56^\circ$|⁠. This is to avoid the effects of equatorial and polar electrojets, which are localized current systems requiring dedicated parametrizations and cannot be well described by low-degree spherical harmonics. Contrary to approaches deriving full geomagnetic models with a focus on core models, such as the CHAOS model (Finlay et al. 2020) or the Comprehensive Model (Sabaka et al. 2020), we do not filter observations based on solar zenith angle or magnetic indices, nor do we parametrize the ionospheric field using specific diurnal/seasonal time harmonics. Since the spatial and temporal behaviour of the external field is of prime interest here, we impose no constraint on the magnetic conditions and inject no prior knowledge on the periodicity of reconstructed sources. Instead, the parameter estimation is carried out independently for each time window, and the temporal behaviour will be analysed in a completely a posteriori fashion. A significant advantage of such model is the ease of its update, which can be performed on the fly in real-time (limited only by the geomagnetic data latency).

3.2 Ground observatory data

We compiled the vector magnetic field data measured by ground observatories from the INTERMAGNET and World Data Centre for Geomagnetism (WDC) within the same time period (INTERMAGNET et al. 2024). The hourly mean time-series are binned in 3-hr time bins for consistency with the satellite data. The same filtering based on QD latitude (⁠|$5^\circ - 56^\circ$| north and south) has been applied. The distribution of ground observatories used in this study is shown in Fig. 2.

$Distribution of ground observatories within $5^\circ$N–$56^\circ$N and $5^\circ$S–$56^\circ$S in magnetic latitude. The colour of the markers indicates the percentage of time bins where the observatory data are available in the period of 2013–2023.$

Figure 2.

Open in new tab Download slide

3.2.1 Observatory bias

One of the difficulties in dealing with ground observatory data lies in the observatory biases. They appear as large offsets in the components of the vector magnetic field, which are not explained by the model, and at the same time are not well correlated even between relatively nearby stations. They are interpreted to represent short-wavelength, localized anomaly crustal fields that are not captured by a lithospheric field model (Langel et al. 1982; Sabaka et al. 2002). In addition, observatory-specific baseline errors introduce biases that cannot be represented in magnetic field models (Lesur et al. 2017).

While these observatory biases are not described by the magnetic field model to be derived, they nevertheless appear as intercepts in the regression and hence affect the parameter estimation process. For this reason, they can be seen as nuisance parameters in our model. Provided an estimate for such biases already exists, one can circumvent this issue by using the time differential ground observations for data fitting (Finlay et al. 2015). Determining the observatory bias from observations, on the other hand, is more challenging. One simple solution is to take the arithmetic mean of each vector component at a ground observatory over magnetic quiet times (Olsen et al. 2014). This serves as a pragmatic approach, effectively removing the means and leaving only zero-mean time-series from which external and the induced fields are estimated, but lacks physical reasoning. From the data perspective, it has long been recognized that this bias cannot be determined unless satellite data unaffected by such bias are used (Langel et al. 1982; for an explanation, see Appendix A). Even when the satellite data are used, the estimation only works when the region between satellite and ground is void of electric currents, implying that no field originates in the ionosphere below satellites (Langel et al. 1982), or when the ionosphere is parametrized using time harmonics (Sabaka et al. 2002). Relaxing the temporal prior knowledge on the ionosphere as we do in this study results in a non-trivial nullspace when the bias is included (Appendix B).

Since the observatory bias cannot be co-estimated in our model with a piece-wise constant temporal parametrization of the ionosphere, we derive these biases in a separate step. We took ground observatory and satellite data from magnetic quiet times (defined as times when |$K_p \le 2$| and |$|dDst/dt| \le 3\mathrm{nT/h}$|⁠) with the Sun at least |$10^\circ$| below the horizon. Data fulfilling these conditions are referred to as quiet night data. The quiet night data are then fitted using a magnetic field model with only magnetospheric field and the internally induced field, and the observatory bias consisting of one time independent vector per observatory. Additional vector biases are introduced for observatories where large offsets are seen in different segments of observation. Due to the limited spatial coverage of the filtered data, the magnetospheric field and the internally induced field are parametrized only up to SH degree and order one. Although the dark-side data are not fully free of ionospheric and its corresponding induced field contributions, such effects are minimal due to much weaker ionospheric currents (Price & Wilkins 1963). Using a parametrization without the ionosphere on the quiet night data is a necessary compromise since the addition of the ionospheric source would lead to an intrinsic non-uniqueness in the bias estimation problem (Appendix B). In addition, note that the nullspace of the forward problem with the ionosphere is characterized by a static field. Therefore, the omission of the ionosphere in the bias estimation can only introduce a small static offset for the entire data span (10 yr in our model), which does not interfere with the temporal variability of the field estimates.

The bias values are then determined using a computation- and memory-efficient linear least-squares algorithm for block-diagonal matrix with intercept (Appendix C). The estimated observatory biases are shown in Supplementary Fig. S1. The observatory biases derived in this study agree well with the CI-derived biases, and both are close to the quiet night mean values (Supplementary Fig. S2, S3 and Table S1). The reason to carry out this estimation instead of using the CI-derived biases is to retain the consistency with our data set, and with our model which does not prescribe a periodic behaviour of the ionospheric field, whereas CI does. These vector biases are then removed from the ground observatory measurements.

The pre-processing of the data yields a total number of 7.7 million ground observatory measurements, and 15.0 million satellite measurements of the three-component vector magnetic field. An average time bin contains 800 magnetic field vector component data sampled by (on average) 89 ground observatories, and 1558 vector component data sampled at 519 satellite locations. These numbers vary in time depending on the availability of satellites and observatories, as is shown in Fig. 3.

Figure 3.

Average number of data points per time bin as a function of time. Each vector component is treated as an independent datum.

Open in new tab Download slide

4 MODEL CONSTRUCTION AND SELECTION

The magnetic field model is parametrized in three sets of source-specific Gauss coefficients: |$g_{nm}^\mathrm{int}(t)$| and |$h_{nm}^\mathrm{int}(t)$| for the internally induced field, |$q_{nm}^\mathrm{ion}(t)$| and |$s_{nm}^\mathrm{ion}(t)$| for the ionospheric field, and |$q_{nm}^\mathrm{mag}(t)$| and |$s_{nm}^\mathrm{mag}(t)$| for the magnetospheric field (eqs 13 and 14). Denoting these Gauss coefficients in the i-th time bin jointly as |$\mathbf {m}_i$|⁠, and the observatory and satellite data jointly as |$\mathbf {d}_i$|⁠, the model estimation procedure in each time bin can be formulated as solving the following optimization problem

$$\begin{eqnarray} \hat{\mathbf {m}}_i = \arg \min _{\mathbf {m}_i} L \left(\mathbf {W}_i (\mathbf {d}_i - \mathbf {G}_i \mathbf {m}_i) \right) , \end{eqnarray}$$

(15)

where |$\mathbf {W}_i$| is the data weight matrix for the data in the i-th time bin, |$\mathbf {G}_i$| is the Gauss matrix and |$L(\cdot )$| is a generic misfit function. We assume independent sampling of each data point and assign uniform variance to each observation, hence |$\mathbf {W} = \mathbf {I}$|⁠. The Gauss matrix |$\mathbf {G}_i$| is the design matrix for the regression problem, and can be obtained by evaluating the gradients of the potentials in eqs (13) and (14) at observation locations.

Since no temporal behaviour is imposed on any of the SH modes, all the modes and their coefficients with a spherical degree no greater than the truncation degree are treated equally in eqs (13) and (14). Since spherical harmonics form an irreducible representation of the 3-D rotation group, there is no difference between using different coordinate systems in our formulation. We chose to work in geocentric geographic coordinates to simplify integration with the existing data set, as well as with other ground-based observations that are naturally provided in this coordinate system.

The truncation degree for the spherical harmonics is selected based on two arguments. First, the local time coverage of the satellites gives a theoretical upper bound on the maximum azimuthal wavenumber. Our data set contains data from six satellites in polar orbits, which, with an approximate orbital period of 90–100 min, cover 24 local times within a 3-hr time bin (each satellite covers two local times per orbit, and completes two orbits per time bin). However, among these satellites, two pairs always share the same local times as they are very close in orbits (Swarm A and C, Grace-FO 1 and 2), and Swarm B satellite also evolves to the same orbital plane as Swarm A and C around 2021. Therefore, the satellites at best sample the magnetic field at 12 local times within one time bin. According to the Nyquist sampling theorem, the field that can be reconstructed without aliasing is limited to azimuthal wavenumber |$m < 6$|⁠, even with uniform sampling. Beyond this azimuthal wavenumber, even if the combined external field can be estimated from the observatory data set, it can no longer be convincingly separated into ionospheric and magnetospheric contributions without being aliased by the satellites. As we are using |$N=M$|⁠, this limits the truncation SH degree to |$N < 6$|⁠.

Given |$N < 6$|⁠, the specific parametrization is further selected using K-fold cross-validation (CV), a technique that provides an estimate on the generalization error of the model. In a K-fold CV, the full data set to be used for parameter estimation is split into K parts. Model parameters are estimated in each fold using |$K-1$| parts of the data sets, while the remaining part serves as the validation set for testing how the estimated model performs on unseen data. The average of the K validation errors then yields the proxy for the generalization error (Arlot & Celisse 2010). While a pure extension of the model parametrization will always reduce the total data misfit (provided global minimum is reached), the generalization error will stagnate or deteriorate when the model is being increasingly fitted to the noise rather than the true signal, and hence provides means to detect overfitting.

We performed a 5-fold CV on the combined data set in each time bin, for SH truncation degrees ranging from 3 to 7. The performance of the model is evaluated using the coefficient of determination, defined as

$$\begin{eqnarray} R^2 = 1 - \frac{\Vert \mathbf {d} - \hat{\mathbf {d}}\Vert _2^2}{\Vert \mathbf {d}\Vert _2^2 - N_d \overline{d}^2} , \end{eqnarray}$$

(16)

where |$\mathbf {d}$| is the observed data vector, |$\hat{\mathbf {d}}$| is data predicted by a model, |$\overline{d}$| is the arithmetic mean of |$\mathbf {d}$| and |$N_d = \dim \mathbf {d}$| the dimension of the data. For each of the model parametrizations, the coefficients of determination |$R^2$| are computed on the validation set of the ground observatory and satellite data sets separately and averaged over the CV folds. Statistics of the CV |$R^2$| score shows that while the estimated model seems to perform better as the truncation degree increases on the satellite data set (blue line in Fig. 5), the CV score for the ground observatory data set stagnates at |$N=5$| and then decreases for |$N=7$| (blue line in Fig. 4). This observation is consistent with the theoretical argument that SH modes with degrees |$>5$| cannot be resolved with the current data. A more careful investigation shows that a model with induced and magnetospheric fields parametrized to |$N=4$|⁠, but ionospheric field extended to |$N=5$| already explains most of the variance in the data and is located on the CV score plateau of |$N=5$|⁠. Compared to this model, the 25% more free parameters in the |$N=5$| model only yields a marginal increase of 0.002 in the CV score. This implies that the extra degrees of freedom in the |$N=5$| model does not yield a significant gain in explaining the data, while increasing the risk of overfitting. Based on these observations, we adopted the parametrization in which the induced and magnetospheric fields are expanded to |$N=4$|⁠, and the ionospheric field is expanded to |$N=5$|⁠. For each time bin, Gauss coefficients |$\mathbf {m}_i$| are estimated by solving eq. (15) using a robust Huber loss as the loss function |$L(\cdot )$|⁠.

Figure 4.

The coefficient of determination |$R^2$| from fitting the whole ground observatory data set and from CV for models estimated using SH truncation degrees from 3 to 7. The square markers show the median of the |$R^2$| scores, while the upper and lower error bar limits show the 1st and 3rd quantiles, respectively.

Open in new tab Download slide

Figure 5.

Same as Fig. 4, but for the satellite data set.

Open in new tab Download slide

5 RESULTS

5.1 The role of ionosphere in the model

By estimating coefficients in each time bin, we obtained the 10-yr time-series of Gauss coefficients describing the magnetospheric field, the ionospheric field and the induced counterpart. The introduction of the ionospheric field in our approach, as described in Section 2, is the major reason all three fields can be characterized individually. While the inclusion of the ionosphere is dictated by reality, it results in additional model parameters and requires the sheet current approximation. Hence, a statistical experiment is needed to justify that adding the ionosphere indeed results in a better model.

To this end, we constructed another model with the same model parametrization, but without the ionosphere layer. This model effectively attributes all satellite and observatory data variations to a common internal source (internally induced currents) and a common external source (magnetospheric currents). Comparing the results from the reduced model (model without the ionosphere) to the results we obtained using the complete model with the ionosphere provides a data-based criterion for justifying the addition of the ionospheric component into our model.

Naturally, introducing an ionospheric field is able to explain more variance in the data simply by expanding the parameter space. The data fit alone is, hence, not a good indicator of the quality of a model, as the model can also achieve lower misfits by fitting noise. To objectively compare the quality of the model and rule out the possibility of overfitting, we once again employ the CV approach. Using a 5-fold CV, we see that there is a drastic improvement in the CV averaged |$R^2$| score (Fig. 6) when an ionosphere layer is included. The parametrization with an ionosphere is thus not only better at explaining the data variance overall, but can explain the unseen observations better. This strongly indicates that the model parametrization with an ionosphere provides better descriptions of the geomagnetic field.

Figure 6.

Histogram of the coefficient of determination |$R^2$| based on the 5-fold cross-validation.

Open in new tab Download slide

As is evident from different |$R^2$| values, introducing the ionosphere in our model parametrization certainly changes the data fit. It is however unclear how the introduction of the ionospheric field affects the existing field estimates—that is, how the estimates of the induced and magnetospheric magnetic fields are affected by the introduction of the ionosphere. To further quantify the impact of the ionospheric layer on the magnetospheric field coefficient estimates, we turn to two statistical overviews of our field estimates. First, as an overview of the difference in the Gauss coefficients, we compute the root-mean-square (RMS) differences between the magnetospheric field coefficient time-series estimated with and without the ionosphere, defined as

$$\begin{eqnarray} \mathrm{RMSD}_{nm}^\mathrm{mag} = \left\lbrace \begin{aligned} \sqrt{\frac{1}{N_t} \sum _{i=1}^{N_t} \left[q_{nm}^\mathrm{mag}(t_i) - \widetilde{q}_{nm}^\mathrm{mag}(t_i)\right]^2},\quad m \ge 0 \\ \sqrt{\frac{1}{N_t} \sum _{i=1}^{N_t} \left[s_{n|m|}^\mathrm{mag}(t_i) - \widetilde{s}_{n|m|}^\mathrm{mag}(t_i)\right]^2},\quad m < 0 \end{aligned}\right. \end{eqnarray}$$

(17)

where |$N_t$| is the total number of time bins. The absolute RMS differences, together with their relative values using the full model estimates as a reference, are shown in Fig. 7 for different spherical harmonics. The notation |$\widetilde{f}$| here denotes the field coefficients f estimated in the reduced model, without an ionosphere. In the notation for the RMS differences, a negative m is used to denote sine modes, whereas the non-negative m is reserved for cosine modes. For |$q_{10}^\mathrm{mag}$|⁠, the RMS difference with and without an ionosphere is 2.4 nT during storm time (⁠|$Kp \ge 5$|⁠), and as low as 0.86 nT during quiet time (⁠|$Kp \le 2$|⁠). Due to the strong magnetospheric |$P_1^0$| signal, these differences translate to merely |$3.8\, \rm per\,cent$| and |$4.7\, \rm per\,cent$| during storm time and quiet time, respectively. The similarity of the two model estimates also extends to other spherical harmonics, such as the |$q_{21}^\mathrm{mag}$| component. For this mode, the estimates with and without an ionosphere have an RMS difference of 2.1 nT during storm time and 1.0 nT during quiet time. Although the RMS differences are always smaller during quiet times, we observe that the relative differences are typically higher (Fig. 7), due to the overall weaker magnetospheric field and greater contributions from the ionospheric field.

$RMS difference for magnetospheric Gauss coefficients estimated using the models with and without the ionosphere. The rows show absolute (upper) and relative (lower) RMS differences, respectively. The columns show the RMS difference for the complete time-series (left), the $Kp\le 2$ subset (middle) and the $Kp \ge 5$ subset (right), respectively.$

Figure 7.

RMS difference for magnetospheric Gauss coefficients estimated using the models with and without the ionosphere. The rows show absolute (upper) and relative (lower) RMS differences, respectively. The columns show the RMS difference for the complete time-series (left), the |$Kp\le 2$| subset (middle) and the |$Kp \ge 5$| subset (right), respectively.

Open in new tab Download slide

Next, we turn to the distribution of the Gauss coefficient pair estimated with and without an ionosphere. In other words, we are interested in the distribution of vectors |$(f_{nm}, \tilde{f}_{nm})$|⁠, where f denotes a Gauss coefficient estimated using our complete model (13–14), and |$\tilde{f}$| again denotes the same quantity estimated using the reduced model without the ionosphere.

If the ionospheric field does not strongly affect the estimate of the Gauss coefficient f, then the distribution of |$(f_{nm}, \tilde{f}_{nm})$| would closely follow the straight line |$\tilde{f}_{nm} = f_{nm}$|⁠. In other words, a linear model in |$f_{nm}$| with a slope close to 1 should explain most of the variance in |$\tilde{f}_{nm}$|⁠. We see that such a linear model is indeed observed for |$q_{10}^\mathrm{mag}$|⁠, the first zonal harmonic of the magnetospheric field estimate (Fig. 8 left panel). The correlation between the two magnetospheric field estimates deteriorates as the SH degree goes higher, for example, for |$q_{21}^\mathrm{mag}$|⁠, as shown in the right panel of Fig. 8. However, while there is a considerably larger variance and deviation from the |$\tilde{q}_{21}^\mathrm{mag} = q_{21}^\mathrm{mag}$| line, the distribution is overall largely unimodal and follows a weaker linear trend. In summary, for the joint ground observatory–satellite data set, the introduction of the ionosphere to the model does not appear to affect strongly the estimates for the field of magnetospheric origin.

Figure 8.

2-D histogram showing the distribution of the magnetospheric field Gauss coefficients estimated with and without the ionosphere layer in the model. The left panel shows the distribution for the |$P_1^0$| mode, and the right panel shows the distribution for the |$P_2^1$| mode. Lines show a linear fit to the distribution of points.

Open in new tab Download slide

The induced field estimates, on the other hand, are much more susceptible to whether or not the ionosphere is co-modelled. This sensitivity is clearly shown in the RMS differences between the estimated made with and without the ionosphere, summarized in Fig. 9. Although the absolute differences are only slightly greater than in the magnetospheric coefficients, the relatively differences are much larger. During magnetic quiet time, where the ionospheric signal is significant in proportion to the magnetosphere, the induced field coefficients are typically over |$80\, \rm per\,cent$| difference for degrees |$n\ge 2$|⁠, rendering the estimates with and without an ionosphere drastically different. Even the most coherent mode, |$g_{10}^\mathrm{int}$|⁠, shows |$> 40\, \rm per\,cent$| difference between the models with and without an ionosphere during quiet times.

Figure 9.

RMS difference for Gauss coefficients of the internally induced field estimated using the models with and without the ionosphere. The rows and columns are the same as in Fig. 7.

Open in new tab Download slide

In the |$g_{10}^\mathrm{int}$| mode, the difference between the two model estimates can be seen as the increased variability of data points around a linear trend (Fig. 10 left panel), compared to the same spatial mode for the magnetospheric field (Fig. 8 left panel). The linear model completely fails for spatial modes such as |$g_{21}^\mathrm{int}$|⁠, where the |$R^2$| score decreases to a mere 0.34. Visualization of its distribution shows that the density is at least bimodal in the |$(g_{21}^\mathrm{int}, \tilde{g}_{21}^\mathrm{int})$| space. This is suggestive of a complex contribution to the difference in the estimated time-series that cannot be described by a simple scaling in the amplitude, but must also include other effects, such as phase shifts. The phase content of the spectrum would therefore be substantially different from one estimate to another, which in turn yields considerable differences in the transfer function estimates, as will be seen in Section 6. In summary, we observe that the introduction of the ionosphere to the model has a considerable effect on the induced field estimate in both the amplitude and phase, especially when the ionospheric field is strong.

Figure 10.

2-D histogram showing the distribution of the induced field Gauss coefficients estimated with and without ionosphere. The left panel shows the distribution for the |$P_1^0$| mode, and the right panel shows the distribution for the |$P_2^1$| mode. Lines show the linear fit to the distribution of points.

Open in new tab Download slide

5.2 Time-series of external and induced coefficients

We have seen in the previous section the time-cumulative distribution of the Gauss coefficients estimated with and without the ionosphere layer. Now we investigate these coefficients further by disentangling the distribution into time-series. The time-series of the estimated magnetospheric field coefficients generally show relatively small differences between different models (Figs 11–12), consistent with our observation on their distributions. Fig. 11 shows time-series of the estimated first zonal harmonic of the magnetospheric field (⁠|$q_{10}^\mathrm{mag}$|⁠) for a period with a geomagnetic storm and a quiet period. The estimates from the complete and reduced models almost coincide with each other, consistent with the low RMS difference (Fig. 7) and almost the perfect correlation between the two (Fig. 8). The close match between the model estimates with and without the ionosphere can also be seen in the |$P_2^1$| mode (Fig. 12), where the two models exhibit nearly coincident peaks during the storm, and in-phase diurnal oscillations at quiet times. Discrepancies are however visible in the amplitudes of the |$P_2^1$| estimates, especially during quiet times, which, combined with the lower energy, contribute to the aforementioned increased relative RMS difference (Fig. 7).

$Time-series of the magnetospheric Gauss coefficient $q_{10}^\mathrm{mag}(t)$ around the 2017 September geomagnetic storm (upper panel) and during magnetic quiet time, in January / February 2018 (lower panel). The four curves in each plot show our model estimates, the estimates with no ionosphere parametrization, the estimates from the CI model and estimates from the CHAOS-7 external field model.$

Figure 11.

Time-series of the magnetospheric Gauss coefficient |$q_{10}^\mathrm{mag}(t)$| around the 2017 September geomagnetic storm (upper panel) and during magnetic quiet time, in January / February 2018 (lower panel). The four curves in each plot show our model estimates, the estimates with no ionosphere parametrization, the estimates from the CI model and estimates from the CHAOS-7 external field model.

Open in new tab Download slide

$Time-series of the magnetospheric Gauss coefficient $q_{21}^\mathrm{mag}(t)$ around the 2017 September geomagnetic storm (upper panel) and during magnetic quiet time, in January / February 2018 (lower panel).$

Figure 12.

Time-series of the magnetospheric Gauss coefficient |$q_{21}^\mathrm{mag}(t)$| around the 2017 September geomagnetic storm (upper panel) and during magnetic quiet time, in January / February 2018 (lower panel).

Open in new tab Download slide

In addition to the same model without the ionosphere, we also compare our estimates with the external field model from the CI model (Sabaka et al. 2020) and the CHAOS-7 model (Finlay et al. 2020). All model estimates for |$q_{10}^\mathrm{mag}$| capture well the two successive geomagnetic storms (as two peaks in the time-series) on 2017 September 8 (Tassev et al. 2017; Zhang et al. 2019), and are in overall good agreement during recovery phases of the geomagnetic storm and during magnetic quiet periods (Fig. 11). Discrepancies can be observed in the peak amplitude of the field coefficient in the geomagnetic storm, as well as in the oscillations in magnetic quiet times. This is anticipated considering the difference in the data used, and the different parametrization of the external field in different models. We note that CHAOS model estimate in particular shows quasi-diurnal oscillations in the geomagnetic quiet times, which are not visible in our model estimates or CI estimates. The discrepancies are much large for higher degrees, such as in the |$q_{21}^\mathrm{mag}$| estimates (Fig. 12). The CHAOS external field estimation is primarily intended for quiet time on the dark side, and only allows for modes that are stationary in either the solar magnetic coordinates (SM) or the geocentric solar magnetospheric coordinates (GSM) beyond spherical harmonic degree one. As a result, it does not describe the fluctuations in the external field beyond degree one, such as the high-amplitude fluctuations for |$q_{21}^\mathrm{mag}$| during a magnetic storm, which are clearly visible in our model and the CI model estimates. Compared to the CI model estimates, our model yields estimates of |$q_{21}^\mathrm{mag}$| that have a larger dynamic range, especially during geomagnetic storms.

In addition to the magnetospheric field, our model also simultaneously estimates the ionospheric field Gauss coefficients. The |$q_{21}^\mathrm{ion}$| spatial mode (Fig. 13) exhibits a large amplitude in the diurnal band in the geographic frame, reflecting a standing mode in the solar-fixed frame, an important constituent of the Sq variation (Malin 1973; Schmucker 1999). The diurnal variation of |$q_{21}^\mathrm{ion}$| is very clear in the reconstructed time-series not only for magnetic quiet times (Fig. 13 lower panel), but around the geomagnetic storms as well (Fig. 13 upper panel). This suggests that the Sq variation is a robust feature that is modulated by the magnetosphere. For instance, Sq dominant structure quickly recovers after the storm’s peak, although its amplitude is dwarfed by a much larger |$q_{10}^\mathrm{mag}$| magnetospheric signal (Fig. 11). However, we observe a considerable deviation from the quiet time daily variation during the storm’s main phase (Fig. 13 upper panel). Interestingly, this deviation is anticorrelated with the magnetospheric mode |$q_{21}^\mathrm{mag}$| (Fig. 12). This anticorrelation may represent an interaction between the ionosphere and magnetosphere current systems during the main phase of a geomagnetic storm.

$Time-series of the ionospheric Gauss coefficient $q_{21}^\mathrm{ion}(t)$ around the 2017 September geomagnetic storm (upper panel) and during magnetic quiet time, in January / February 2018 (lower panel).$

Figure 13.

Time-series of the ionospheric Gauss coefficient |$q_{21}^\mathrm{ion}(t)$| around the 2017 September geomagnetic storm (upper panel) and during magnetic quiet time, in January / February 2018 (lower panel).

Open in new tab Download slide

As has been observed in Figs 9–10, the induced field estimate is more susceptible to the presence of the ionospheric layer to the model. This is confirmed in the time-series in Figs 14 and 15. The discrepancy is especially pronounced in spatial modes where the ionosphere is comparatively energetic, for example, the |$P_2^1$| SH mode. During the magnetic quiet time (cf. Fig. 15, lower panel), the amplitude of the induced field estimate in the reduced model is less than 50% of that estimated using the full model with an ionosphere. We interpret this systematically smaller amplitude as an imprint of the missing ionospheric field on the induced field. When the ionosphere is neglected, the internal contribution to the satellite data which should have been attributed to the ionosphere is forced into the induced coefficients. This, in turn, may weaken the induction signal when the induced field is correlated to the ionospheric signal (eq. 14). The difference in the model estimates with and without the ionosphere is not only in the amplitude, but also in the phase. The phase change can already be observed in |$g_{21}^\mathrm{int}$|⁠, but it is even more pronounced in higher SH degrees. The semidiurnal oscillation of |$h_{32}^\mathrm{int}$|⁠, for instance, appears to be |$\sim 3$| hr ahead in the model without an ionosphere compared to the full model estimate (Fig. 16).

$Time-series of the internally induced Gauss coefficient $g_{10}^\mathrm{int}(t)$ around the 2017 September geomagnetic storm (upper panel) and during magnetic quiet time, in January / February 2018 (lower panel). The two curves in each plot show our full model estimates and the estimates with no ionosphere parametrization.$

Figure 14.

Time-series of the internally induced Gauss coefficient |$g_{10}^\mathrm{int}(t)$| around the 2017 September geomagnetic storm (upper panel) and during magnetic quiet time, in January / February 2018 (lower panel). The two curves in each plot show our full model estimates and the estimates with no ionosphere parametrization.

Open in new tab Download slide

$Time-series of the internally induced Gauss coefficient $g_{21}^\mathrm{int}(t)$ around the 2017 September geomagnetic storm (upper panel) and during magnetic quiet time, in January / February 2018 (lower panel).$

Figure 15.

Time-series of the internally induced Gauss coefficient |$g_{21}^\mathrm{int}(t)$| around the 2017 September geomagnetic storm (upper panel) and during magnetic quiet time, in January / February 2018 (lower panel).

Open in new tab Download slide

$Time-series of the internally induced Gauss coefficient $h_{32}^\mathrm{int}(t)$ during magnetic quiet time, in February 2018.$

Figure 16.

Time-series of the internally induced Gauss coefficient |$h_{32}^\mathrm{int}(t)$| during magnetic quiet time, in February 2018.

Open in new tab Download slide

5.3 Fourier spectrum of external and induced fields

Our 10-yr time-series of the SH coefficients has a theoretical bandwidth of over four decades, that is, |$f_\mathrm{Nyq}/\Delta f > 10^4$|⁠, where |$f_\mathrm{Nyq}$| is the Nyquist frequency and |$\Delta f$| is the frequency resolution. As such, it covers the entire period range from subdiurnal to interannual periods. In this section, we explore the frequency content of the estimated fields by deriving their amplitude spectra. The spectra are computed directly as discrete Fourier transforms of the time-series of the Gauss coefficients, using a forward normalization that preserves the time-domain amplitude of a harmonic signal in the amplitude spectrum (i.e. a sinusoidal signal with 1 nT amplitude will result in a 1 nT peak at the corresponding frequency in the amplitude spectrum).

The first zonal harmonic of the magnetospheric magnetic field, largely attributed to the ring current, is the most energetic contribution to the magnetic field of the external origin at periods longer than one day. The power of the spectrum for the first zonal harmonic quickly rolls off at periods shorter than one day, diurnal frequency bands, while exhibiting strong amplitude at longer periods (Fig. 17, upper panel). At longest periods, a prominent spectral peak can be observed between 2.0 and 2.5 yr. This biennial oscillation has been previously observed in the solar magnetic field, and has been attributed to a harmonic of the solar cycle (Knaack & Stenflo 2005). This notable feature can be similarly observed in the ionospheric signal (Fig. 17 lower panel) and in other degree-one modes (Supplementary Fig. S4). The second, rather sharp peak in the magnetospheric field spectrum is at the semi-annual period (2 cycles per year). Since our estimation is conducted in the Earth-fixed frame, the semi-annual peak is likely associated with the axial tilt of the Earth and the associated seasonal variability in the Earth’s frame (Maus & Lühr 2005; Olsen et al. 2005; Balasis & Egbert 2006). A secondary, wider peak can be observed at the annual period (1 cycle per year). This can be a signature of the tail current in the magnetosphere (Maus & Lühr 2005). As the |$q_{10}^\mathrm{mag}$| component experiences semi-annual changes in the solar radiation, it is further modulated by the change in the interplanetary magnetic field. Both annual and semi-annual peaks are observed not only in the first zonal harmonic component, but also in other dipole components of the magnetospheric field (Supplementary Fig. S4 upper panel). This can result from the transformation of the quasi-static magnetospheric field in GSM/SM frames to the Earth-fixed frame. In contrast, these peaks are weaker or not visible in the ionospheric spectrum (Figs 17 and S4, lower panels), suggesting that ionospheric dynamics at |$n=1$| is not as strongly modulated by the orbit-related periodicities. Note that this observation only holds for the |$n=1$| spherical harmonic signals. The annual, semi-annual even quarterly peaks are clearly visible in both magnetospheric and ionospheric Gauss coefficients of |$n>1$| (Supplementary Figs S5 and S6).

$Amplitude spectra of the magnetospheric coefficient $q_{10}^\mathrm{mag}$ (upper panel), the ionospheric coefficient $q_{10}^\mathrm{ion}$ coefficient (middle panel) and the internally induced coefficient $g_{10}^\mathrm{int}$ (lower panel). Spectral peaks observed in the spectra are marked with triangles, with their peak periods annotated. Selected theoretical periods of the system are plotted as solid vertical lines for reference.$

Figure 17.

Amplitude spectra of the magnetospheric coefficient |$q_{10}^\mathrm{mag}$| (upper panel), the ionospheric coefficient |$q_{10}^\mathrm{ion}$| coefficient (middle panel) and the internally induced coefficient |$g_{10}^\mathrm{int}$| (lower panel). Spectral peaks observed in the spectra are marked with triangles, with their peak periods annotated. Selected theoretical periods of the system are plotted as solid vertical lines for reference.

Open in new tab Download slide

Further, we observe spectral peaks at quasi-monthly periods in the |$q_{10}^\mathrm{mag}$| spectra (Fig. 17 upper panel). These include peaks at |$29.7\pm 0.2$| d and |$25.5\pm 0.2$| d, and a smaller peak at |$27.0\pm 0.2$| d. The uncertainty reported here is the resolution at the given peaks. This set of peaks is a robust feature that also appears in other degree-one mode spectra of the magnetospheric coefficients (Supplementary Fig. S4 upper panel). In contrast, these peaks are weaker (Fig. 17 lower panel) or absent (Supplementary Fig. S4 lower panel) in ionospheric coefficients. We associate these with the well-known |$\sim 27$|-d harmonic in terrestrial magnetic field (Bartels 1932), stemming from mechanisms associated with the solar rotation, which has a synodic Carrington period of |$\sim 27.7$| d (Lockwood & Owens 2024). The strong peaks at 25.5 and 29.7 d are likely due to the variability in the solar rotation, and are consistent with the analysis of the interplanetary magnetic field (IMF) in the solar cycle 24, which shows strong spectral peaks at 26 and 30 d (Chowdhury et al. 2015).

Additional spectral peaks are also observed at multiples or submultiples of the Carrington period. The quasi-fortnightly oscillation, observed at 14.4 d in |$q_{10}^{\mathrm{mag}}$| and |$q_{11}^{\mathrm{mag}}$| and at 14.8 d in |$q_{10}^\mathrm{ion}$|⁠, is likely a 2nd harmonic of the |$\sim 27$| d signal. At periods longer than the Carrington period, spectral peaks can be observed at |$54.9\pm 0.8$| d (Fig. 17 |$q_{10}^\mathrm{mag}$|⁠, Supplementary Fig. S4 |$q_{11}^\mathrm{mag}$|⁠), |$134 \pm 5$| d (Supplementary Fig. S4, |$q_{11}^\mathrm{mag}$|⁠, |$q_{11}^\mathrm{ion}$|⁠), |$157\pm 7$| d (⁠|$q_{10}^\mathrm{mag}$|⁠, ibid.) and |$241\pm 17$| d (⁠|$q_{11}^\mathrm{mag}$|⁠, |$q_{11}^\mathrm{ion}$|⁠, ibid.). This series of spectral peaks, termed the Rieger-type periods, was first observed by Rieger et al. (1984) in the 154-d periodic occurrence of solar flares. The periodicity was later extended to shorter periods of |$\sim 130$| and |$\sim 50$| d (Bai & Sturrock 1991), and the series of periods have since also been observed in the IMF (Cane et al. 1998), coronal mass ejections (Lou et al. 2003) and sunspot occurrence (Chowdhury et al. 2015). Since these peaks are located close to multiples of the solar rotation period, it has long been postulated that the signals associated with these periods are sub-harmonics of a fundamental period at the solar rotation (Bai & Sturrock 1991). It has been proposed that the mechanism underlying this signal is linked to r-modes or equatorial trapped Rossby-type waves in the Sun (Lou 2000). Although the conclusion regarding the mechanism of the Rieger-type periodicity is presently not fully established, the current evidence suggests that these spectral peaks observed in our reconstruction of the magnetospheric and ionospheric fields is likely modulated by periodicities in the Sun’s activity.

In contrast to variations at long periods, the variations in diurnal period band are dominated by the ionospheric sources. One of the major features in the ionospheric signal in this period band is the diurnal mode and its harmonics. These are clearly visible as discrete spectral peaks at one, two, three and four cycles per day (cpd) (Fig. 18, lower panel). For the |$P_2^1$| spatial mode, the fundamental time harmonic with a period of one day is the strongest, with an amplitude two orders of magnitudes above the average spectrum level. This mode, which has an azimuthal wavenumber |$m=1$|⁠, and has a 1-cpd periodicity observed in the Earth-fixed geographic frame, represents a stationary |$n=2$| mode in a solar-fixed frame (local-time term). This mode is predominantly influenced by the day-side structure of the Sq current system (Malin 1973; Winch 1981; Schmucker 1999). The same holds also for the SH coefficients |$s_{32}^\mathrm{ion}$| and |$s_{43}^\mathrm{ion}$| with the azimuthal wavenumbers of |$m=2$| and |$m=3$|⁠, respectively. They exhibit strongest spectral peaks at 2 cpd (Supplementary Fig. S7 lower panel) and 3 cpd (Supplementary Fig. S8 lower panel), respectively. These correspond to the |$n=3$| and |$n=4$| modes that are stationary in a solar-fixed frame. The diurnal periods are also visible in the spectra of magnetospheric coefficients (Fig. 18, Supplementary Figs S7 and S8 upper panels), although they have a lower signal-to-noise ratio compared to the ionosphere.

$Amplitude spectra of the magnetospheric coefficient $q_{21}^\mathrm{mag}$ (upper panel), the ionospheric coefficient $q_{21}^\mathrm{ion}$ coefficient (middle panel) and the internally induced coefficient $g_{21}^\mathrm{int}$ (lower panel) around the diurnal band. Spectral peaks observed in the spectra are marked with triangles, with their peak periods annotated. Selected theoretical periods of the system, including the diurnal period and its harmonics, and the periods for lunar daily variations $L_p$ with $p=1,2,3,4$ are plotted as solid vertical lines for reference.$

Figure 18.

Amplitude spectra of the magnetospheric coefficient |$q_{21}^\mathrm{mag}$| (upper panel), the ionospheric coefficient |$q_{21}^\mathrm{ion}$| coefficient (middle panel) and the internally induced coefficient |$g_{21}^\mathrm{int}$| (lower panel) around the diurnal band. Spectral peaks observed in the spectra are marked with triangles, with their peak periods annotated. Selected theoretical periods of the system, including the diurnal period and its harmonics, and the periods for lunar daily variations |$L_p$| with |$p=1,2,3,4$| are plotted as solid vertical lines for reference.

Open in new tab Download slide

In addition to dominant solar daily variations, we also observe another set of weaker peaks at nearby periods. These periodic signals are due to lunar daily variations (Chapman 1919; Chapman & Miller 1940), hereinafter denoted by |$L_p$| following the original notation by Malin & Chapman (1970). |$L_p$| variations can be well described as harmonic signals that obey the Chapman’s phase law (Chapman & Miller 1940; Malin & Chapman 1970)

$$\begin{eqnarray} L_p(t,\nu ) = l_p\sin {\left( pt - 2\nu +\lambda _p\right)}, \end{eqnarray}$$

(18)

where |$p = 1, 2, 3 \cdots$| corresponds to diurnal, semidiurnal, terdiurnal, etc., terms with t solar time in hours and |$\nu$| the moon phase, measured by the hour angle between the Moon and the Sun. |$l_p$| and |$\lambda _p$| are the corresponding amplitude and phase values. With the Lunar synodic period of 29.53059 d, the periods of the first four lunar daily harmonics |$L_p$| are 25.7435, 12.4206, 8.1847 and 6.1033 hr for |$p = 1\dots 4$|⁠, respectively. These periods are marked in the power spectra of ionospheric SH coefficients in Fig. 18 and Supplementary Figs S7–S8. Note that these peaks are not equally strong across SH coefficients owing to a non-uniform (and daily modulated) spatial structure of the E-region conductivity (Winch 1981). As anticipated, peaks are strongest in the SH coefficients that correspond to local-time terms with degree |$n = p + 1$| and order |$m = p$| (Schmucker 1999).

Although lunar daily variations are apparent in all SH coefficients describing the field of the ionospheric origin, they are absent in any of the SH coefficients describing the field of the magnetospheric origin (Figs 18, Supplementary Figs S7 and S8). This confirms their ionospheric origin and illustrates high quality of the field separation in the present model.

The estimated internal field contains signals that originate from both the magnetospheric and ionospheric field variations. For the |$P_1^0$| component, for instance, the internal coefficient |$g_{10}^{\mathrm{int}}$| inherits the strong semi-annual and Rieger-type spectral peaks as well as the peak at 29.7 and 25.5 d from the magnetospheric signal, whereas these peaks are absent in the ionospheric signal (Fig. 17). For |$n>1$| in diurnal bands, the induced coefficients show clear spectral peaks at |$L_p$| tidal periods, which are only observed in the ionosphere but not in the reconstructed magnetospheric signal (Fig. 18). Additionally, our induced coefficients also contain oceanic tidal magnetic signals (Grayver et al. 2024). However, since we use both day-side and dark-side data and adopt a low-degree spherical harmonic representation of the field, the ocean tidal contributions are superposed by signals resulting from the ionospheric lunar daily variations (Malin & Chapman 1970).

6 DISCUSSION AND CONCLUSIONS

Compared to the previous geomagnetic field models, our new modelling approach with the ionosphere parametrization has several distinct features. First, the magnetic fields with ionospheric and magnetospheric origins are separated. This sets our model apart from the modelling approaches where the ionospheric field is neglected (Lesur et al. 2010; Finlay et al. 2020) or the magnetospheric field is removed by baseline removal and ring current corrections (Yamazaki & Maute 2017; Guzavina et al. 2019). Secondly, no temporal behaviour is prescribed on either of the external fields. This is in contrast to many other models. For instance, the temporal behaviour for magnetospheric fields beyond |$P_1^0$| is prescribed to be standing modes in SM or GSM coordinates in the CHAOS model (Finlay et al. 2020), or the ionospheric field is assumed to consist of daily and seasonal time-harmonics in the comprehensive model (Sabaka et al. 2018). Furthermore, the secondary induced field is co-estimated alongside the external fields, rather than modelled using an existing conductivity model, as is the case in Chulliat et al. (2016).

These new features make our model complementary to other existing model and in some cases more suitable for dedicated external field or mantle induction studies. First, since we do not assume a subsurface electrical conductivity a priori for modelling the induction effects, the co-estimated induced field coefficients combined with the external field coefficients can be used to probe the electrical conductivity in the interior of the Earth. Previous attempts to image the interior conductivity using joint satellite-observatory data set tried to avoid the inconsistency between the two data sets by removing CI-based ionospheric field from the data (see e.g. Kuvshinov et al. 2021; Velímský & Knopp 2021). This approach is limited since the ionospheric model is based on quiet-time data, and is parametrized using daily and seasonal time harmonics. For instance, we have seen in Section 5.2 that during geomagnetic storms, the ionospheric field exhibits strong non-periodic behaviour.

When the geomagnetic field model fails to describe the ionospheric sources properly, the remnant ionospheric signal will be compensated by biasing the induced field and the magnetospheric field estimates. To test whether the induced and the magnetospheric fields are indeed biased if the ionospheric field is neglected, we need a physics-based metric to examine how the model with an ionosphere compares with a model without the ionosphere. For this purpose, we estimated the transfer functions (TFs) relating the induced and the external parts of the geomagnetic potential. In the frequency domain, assuming a radially symmetric electrical conductivity of the Earth, the Gauss coefficients of the induced field are related to those of the external field as dic occurrence of solar flares. The periodici

$$\begin{eqnarray} \iota _{nm}(\omega ) = Q_n(\omega ; \sigma ) \, \varepsilon _{nm}(\omega ), \end{eqnarray}$$

(19)

where |$\omega$| is the angular frequency, and |$Q_n$| is the so-called Q-response (Schmucker 1985), a functional of the subsurface conductivity |$\sigma$|⁠. For a general 3-D subsurface conductivity, the relation above is not complete since a single external mode would induce many internal modes, all of which are coupled through a so-called Q-matrix (Olsen 1999a). As the subsurface conductivity can be approximated to the first order as a radial conductivity profile, the diagonal elements of this matrix dominate, and the relation between the induced and external coefficients can to first order be approximated by the expression above. Therefore, the equation above is still a valid test since any realistic 3-D conductivity effects will be dwarfed by much stronger source effects.

The complex Gauss coefficients |$\iota _{nm}$| and |$\varepsilon _{nm}$| are the expansion coefficients of the geomagnetic potential V using the complex SH basis |$Y_n^m(\theta , \phi ) \propto P_n^{|m|}(\cos \theta ) e^{im\phi }$| and are related to their real counterparts as

$$\begin{eqnarray} \varepsilon _{nm} = \left\lbrace \begin{aligned} &q_{n0},& &m=0 \\ &\frac{1}{2} \left(q_{n|m|} - i s_{n|m|}\right),& &m > 0, \\ &\frac{1}{2} \left(q_{n|m|} + i s_{n|m|}\right),& &m < 0, \end{aligned}\right.\qquad \iota _{nm} = \left\lbrace \begin{aligned} &g_{n0},& &m=0 \\ &\frac{1}{2} \left(g_{n|m|} - i h_{n|m|}\right),& &m >0, \\ &\frac{1}{2} \left(g_{n|m|} + i h_{n|m|}\right),& &m < 0. \end{aligned}\right. \end{eqnarray}$$

(20)

For the model with the ionosphere layer, both ionospheric and magnetospheric fields are external, and hence we use |$\varepsilon _{nm} := \varepsilon _{nm}^\mathrm{ion} + \varepsilon _{nm}^\mathrm{mag}$| for |$\varepsilon _{nm}$| in (19). The internal coefficients are simply the Gauss coefficients for the induced field, thus we have |$\iota _{nm}:=\iota _{nm}^\mathrm{int}$|⁠. For the model without ionosphere, the magnetospheric and the induced fields serve as the external and the internal fields, respectively.

For a radially symmetric electrical conductivity of the Earth, the Q-response is only dependent on the SH degree n, but is independent of the order m. Therefore, an estimation of |$Q_n$| can be conducted using any spherical harmonic mode of the same degree. The results are further converted to the global C-response (Schmucker 1985; Olsen 1999a) as

$$\begin{eqnarray} C_n(\omega ) = \frac{a}{n+1} \frac{1 - \frac{n+1}{n} Q_n(\omega )}{1 + Q_n(\omega )}. \end{eqnarray}$$

(21)

The real part of the |$C_n$|-response is a depth of the maximum induced current, thus approximating the sounding depth at a given frequency (Weidelt 1972). The value of |$C_n$| at |$\omega$| is estimated following the general method outlined in Olsen (1998). First, the time-series of the complex Gauss coefficients |$\varepsilon _{nm}$| and |$\iota _{nm}$| (obtained from real counterparts using eq. 20) are split into segments with length |$3T$|⁠, where |$T = 2\pi /\omega$| is the period. Secondly, the short-time Fourier transform of the time-series with Hamming time window is computed for every segment |$\tau _j$| (⁠|$j=1,2,\cdots N_\tau$|⁠), generating the spectrum points within the segment, denoted as |$\iota _{nm}(\tau _j, \omega )$| and |$\varepsilon _{nm}(\tau _j, \omega )$|⁠. The final transfer function is estimated by performing a robust regression on the linear system (19), with |$\varepsilon _{nm}(\tau _j, \omega )$| and |$\iota _{nm}(\tau _j, \omega )$| estimated in the previous step. The resulting |$C_n$|-responses from the field estimates of the model with and without the ionosphere are visualized in Fig. 19 at 9 log-spaced periods between 1.5 and 10 d, as well as at the diurnal (24 hr), semidiurnal (12 hr) and terdiurnal (8 hr) periods. It is also possible to estimate these transfer functions at longer periods, but the effect of the ionosphere becomes smaller at periods longer than 20–30 d. To evaluate the linear correlation between the internal and external fields, we also computed the squared coherence (⁠|$\mathrm{Coh}^2$|⁠) for each spherical harmonic mode and each frequency, defined as

$$\begin{eqnarray} \mathrm{Coh}^2 = \frac{\left|\sum _{j=1}^{N_\tau } \iota ^{*}_{nm}(\tau _j, \omega ) \, \varepsilon _{nm}(\tau _j, \omega )\right|^2}{\sum _{j=1}^{N_\tau }\left|\iota _{nm}(\tau _j, \omega )\right|^2 \sum _{j=1}^{N_\tau }\left|\varepsilon _{nm}(\tau _j, \omega )\right|^2}. \end{eqnarray}$$

(22)

Figure 19.

|$C_n$|-responses (upper panel) and their squared coherences (lower panel) at periods between 8 hr and 10 d. The columns show the |$C_1$|⁠, |$C_2$|⁠, |$C_3$| and |$C_4$| responses estimated from the SH coefficients described by the |$P_1^0$|⁠, |$P_2^1$|⁠, |$P_3^2$| and |$P_4^3$| modes, respectively.

Open in new tab Download slide

The TFs for the model with the ionosphere have generally higher squared coherence, especially for higher SH degress where the contribution of the ionosphere is expected to be larger. For instance, |$C_3$| estimated from the |$P_3^2$| mode boasts a squared coherence of |$\mathrm{Coh}^2 = 0.95$| at the dominant period of 12 hr, a drastic improvement over the mere |$\mathrm{Coh}^2 = 0.42$| for the model without the ionosphere; |$C_4$| estimated from the |$P_4^3$| mode has |$\mathrm{Coh}^2 = 0.78$| at the period of 8 hr, compared to |$\mathrm{Coh}^2=0.43$| for the model that neglects ionosphere (Fig. 19). Similar but somewhat smaller advantages can also be observed for |$C_2$| estimated from the |$P_2^1$| mode around the diurnal period band, and across almost the entire period range for |$C_3$| and |$C_4$|⁠, albeit both have lower coherences at longer periods. For the |$C_1$| estimated from the first zonal SH |$P_1^0$|⁠, we note that the transfer functions from both models deteriorate as they approach the diurnal band compared to longer periods, where the two estimates are nearly identical in terms of coherence.

Besides improved coherences, |$C_n$|-response values appear more physical when the ionosphere is modelled. In higher degree responses (⁠|$C_2$|⁠, |$C_3$| and |$C_4$|⁠), the C-responses derived using the reduced model without the ionosphere feature spurious fluctuations and non-monotonic behaviour, occasionally changing the sign in the real parts (e.g. |$C_3$| and |$C_4$| at diurnal periods). This behaviour is not physical since global |$C_n$|-responses are strictly monotonic for a 1-D model and can deviate from this in 3-D settings, but such deviations are smooth and mild for realistic 3-D models (Püthe & Kuvshinov 2014; Grayver et al. 2021). In contrast, these strong fluctuations and unrealistic values of the |$C_n$|-responses are completely absent in the TFs computed from the model with the ionosphere.

The difference in the estimated transfer functions and their coherence shows that a lack of proper treatment of the ionosphere in the model results in transfer functions that are in some cases significantly incoherent and physically inconsistent. Fortunately, this does not affect the |$C_1$|-response computed from the first zonal harmonic at periods longer than one day, rendering previous global reference models (e.g. Grayver et al. 2017; Kuvshinov et al. 2021) still valid. The situation is different for the 3-D mantle inversions. Previous studies (Kuvshinov et al. 2021; Velímský & Knopp 2021) used external and internal coefficients derived from satellite and ground observations whereby the ionospheric effect was not explicitly modelled. The authors subtracted a model of the ionosphere (together with the corresponding induced field) with a prescribed periodic behaviour, constructed using quiet-time magnetic data as a part of the CI model (Sabaka et al. 2018). Failing to co-model ionospheric sources or extrapolating quiet-time ionosphere models to the entire data set will likely result in source effects being propagated to 3-D conductivity models. To illustrate this, we also first subtracted a CI-derived ionospheric model (with its induced field) from our data set and then carried out the internal-external field separation. The resulting magnetospheric spectra show lunar tidal peaks (Supplementary Fig. S9), implying the ionospheric signal may have leaked into the magnetospheric estimates. At the same time, there is almost no improvement and even occasional deterioration in the coherence of the estimated TFs (Supplementary Fig. S1) compared to the internal-external field separation without the ionosphere on the original data set (Fig. 19). Therefore, explicitly co-estimating the ionosphere is necessary to achieve coherent and physical transfer functions for probing the 3-D subsurface conductivity of the Earth with satellite and ground data.

Besides new opportunities for induction studies, our model also enables a more flexible and versatile characterization of the external magnetic field sources. Thanks to the relaxed assumptions on temporal behaviour of the fields and the environmental setting, our model is capable of exploiting data across all local times and during all magnetic conditions. To get a glimpse of the spacial structure and the temporal behaviour of the estimated external field, we visualize the results in terms of the streamfunction of an equivalent current, defined as

$$\begin{eqnarray} \Psi (\mathbf {r}, t) &=& - \frac{a}{\mu _0} \sum _{nm} \frac{2n+1}{n+1} \left(\frac{a + h}{a}\right)^n \left[q_{nm}(t) \cos m\phi + s_{nm}(t)\sin \phi \right] \\ && P_n^m(\cos \theta ). \end{eqnarray}$$

(23)

The electric current is then linked to the streamfunction via

$$\begin{eqnarray} \mathbf {j}(\mathbf {r}, t) = - \delta (r - (a+h)) \, \hat{\mathbf {e}}_r \times \nabla _H \Psi (\mathbf {r}, t) , \end{eqnarray}$$

(24)

where |$\hat{\mathbf {e}}_r$| is the unit vector in the radial direction, and |$\nabla _H$| is the surface gradient. On the surface of a sphere with radius b, the surface gradient can be explicitly written as |$\nabla _H = \hat{\mathbf {e}}_\theta b^{-1} \partial _\theta + \hat{\mathbf {e}}_\phi (b\sin \theta )^{-1} \partial _\phi$|⁠. This formulation, describing purely toroidal electric currents concentrated in an infinitely thin layer at a radius |$a+h$|⁠, is consistent with the current sheet approximation for the ionospheric currents (Section 2). Therefore, the ionospheric equivalent current streamfunction |$\Psi ^{\mathrm{ion}}$| can be calculated using the expansion (23) and |$q_{nm}^\mathrm{ion}(t)$|⁠, |$s_{nm}^\mathrm{ion}(t)$|⁠.

During geomagnetic quiet times, our model recovers the dominant Sq current system in the ionosphere, as is shown in Fig. 20. The equivalent current of the Sq current system features one large-scale counterclockwise vortex in the Northern Hemisphere, and one large-scale clockwise vortex in the Southern Hemisphere. Near the summer solstice, the foci of the Sq currents in the Northern Hemisphere are observed to consistently lag behind in local time compared to the southern foci, consistent with the previous observations on the seasonal variation in the local times of the focus positions (Yamazaki & Maute 2017). Near the spring equinox, the northern and southern foci of the Sq currents alternate to lead in local times on a semidiurnal basis. These seasonal changes of the Sq system are also reflected by the annual, semi-annual and seasonal periodicities in the ionospheric Gauss coefficient spectra at spherical harmonic degree |$n > 1$| (Supplementary Figs S5 and S6, lower panels).

$Streamfunctions of the equivalent current for the ionospheric field ($\Psi ^{\mathrm{ion}}$) on magnetic quiet days close to a spring equinox (left column) and a summer solstice (right column). The streamfunction is cut off above an absolute latitude of 70 degrees to filter out the values in high latitudes where our data set and parametrization have limited resolution. The central longitude for each subplot is at the local noon.$

Figure 20.

Streamfunctions of the equivalent current for the ionospheric field (⁠|$\Psi ^{\mathrm{ion}}$|⁠) on magnetic quiet days close to a spring equinox (left column) and a summer solstice (right column). The streamfunction is cut off above an absolute latitude of 70 degrees to filter out the values in high latitudes where our data set and parametrization have limited resolution. The central longitude for each subplot is at the local noon.

Open in new tab Download slide

During geomagnetic storms, the estimated ionospheric current deviates considerably from a quiet Sq current, exemplified by the ionospheric equivalent current reconstructed during the geomagnetic storm on 2017 September 8 (Fig. 21). Compared to the Sq current during quiet times, the equivalent sheet shows two distinct features during geomagnetic storms. First, the day–night dichotomy, ubiquitous during geomagnetic quiet times, is no longer present. Instead, the sheet current shows strong current amplitude across all local times (see UT=1:30 and 13:30 on 8 September). Secondly, the structurally robust vortices centred at mid-latitudes in the Northern and Southern Hemispheres give way to multiple transient, rapid changing vortices centred at high latitudes (see UT = 22:30, 1:30, 13:30, 16:30 on 8 September, when |$Kp \ge 7$|⁠). This may imply the energy input from the magnetosphere in high latitudes, which then propagates to lower latitudes by virtue of electric field penetration (Yamazaki & Maute 2017). Interestingly, we see that in between the time bins with |$Kp \ge 7$|⁠, the Sq current system still seems to persist, observable as the two vortices in the Northern and Southern Hemispheres.

$Streamfunction of the equivalent current for the ionospheric field ($\Psi ^{\mathrm{ion}}$) during a geomagnetic storm on 2017 September 8. The same cutoff latitude of 70 degrees is used. Note that the maximum value of the colour scale is extended to twice of that for the quiet day plot (Fig. 20) to accommodate higher amplitudes of the streamfunction during the storm. Nevertheless, the strong equivalent currents saturate the colour scales when $Kp > 7$.$

Figure 21.

Streamfunction of the equivalent current for the ionospheric field (⁠|$\Psi ^{\mathrm{ion}}$|⁠) during a geomagnetic storm on 2017 September 8. The same cutoff latitude of 70 degrees is used. Note that the maximum value of the colour scale is extended to twice of that for the quiet day plot (Fig. 20) to accommodate higher amplitudes of the streamfunction during the storm. Nevertheless, the strong equivalent currents saturate the colour scales when |$Kp > 7$|⁠.

Open in new tab Download slide

Our model provides a continuous description of the mid-latitude ionospheric and magnetospheric fields for the period 2014–2023, covering all local times and magnetic conditions. The temporal resolution of the presented model is 3 hr, indicating that dynamics on a timescale that is similar to or less than 3 hr cannot be well captured. This is especially pronounced in phases of geomagnetic storms with fast Dst variations, where we observe higher misfits of the model (⁠|$R^2 \sim 0.7$| for UT = 22:30, 01:30, 13:30, Fig. 21). In comparison, the misfits during other times tend to be lower (⁠|$R^2 > 0.85$| for all other snapshots). In addition, the spatial parametrization of the current model is up to SH degree and order 4 for induced and magnetospheric fields, and up to degree and order 5 for ionospheric field. Localized features, such as equatorial electrojets and the polar current system, cannot be described within our parametrization and hence have been avoided in data selection. The (in-)ability to capture certain temporal and spatial variations of the magnetic fields is controlled by a trade-off between space and time coverage of the data, which is ultimately dictated by the distribution of the observations.

In light of these limitations, our framework can benefit both from the incorporation of additional data, as well as dedicated parametrization for localized current systems. The new low-orbit geomagnetic satellites, such as the Macau Science Satellites (MSS) and a planned ESA NanoMagSat mission, will significantly improve the spatial (especially local-time) coverage, allowing us to increase the spatiotemporal resolution of the model. Dedicated parametrization of polar electrojets based on loop or point currents can be integrated into our model to improve the fit at polar regions. Thanks to its flexibility, the model can be updated on the fly in real time, thus making it particularly suitable for working with low-time-lag data sets, such as the Swarm FAST data, in support of space weather operations and hazard assessment.

AUTHOR CONTRIBUTIONS

Jingtao Min (Conceptualization, Data curation, Investigation, Methodology, Resources, Software, Visualization, Writing—original draft, Writing—review and editing), and Alexander Grayver (Conceptualization, Investigation, Methodology, Project administration, Resources, Visualization, Writing—original draft, Writing—review and editing).

ACKNOWLEDGMENTS

JM is grateful for funding from the European Research Council (agreement No. 833848-UEMHP) under the Horizon 2020 program, and the funding from Swiss National Science Foundation (grant No. 219247) under the MINT 2023 call. AG was supported by the Heisenberg Grant from the German Research Foundation, Deutsche Forschungsgemeinschaft (Project No. 465486300) and European Space Agency Swarm DISC project No. 4000109587.

DATA AVAILABILITY

The input Swarm, CryoSat-2 and Grace-FO data are part of the Swarm product and are available at https://swarmhandbook.earth.esa.int/catalogue/. The input ground observatory data are from INTERMAGNET and WDC data, accessed via the VirES platform (https://vires.services/). The CI core and lithospheric models subtracted from the data and the CI magnetospheric model used for comparison are available at https://swarmhandbook.earth.esa.int/catalogue/. The CHAOS-7.18 model used for comparison is available at https://spacecenter.dk/files/magnetic-models/CHAOS-7. Time series of estimated ionospheric and magnetospheric Gauss coefficients are available on Zenodo, at https://dx.doi.org/10.5281/zenodo.14787587.

Footnotes

VirES for Swarm service, https://vires.services/

REFERENCES

Alken

et al. ,

2021

International Geomagnetic Reference Field: the thirteenth generation

Earth Planets Space

(

10.1186/s40623-020-01288-x

Arlot

Celisse

2010

A survey of cross-validation procedures for model selection

Stat. Surv.

–

Backus

Parker

R.L.

Constable

1996

Found. Geomagn.

Cambridge Univ. Press

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Baerenzung

Holschneider

Wicht

Lesur

Sanchez

2020

The Kalmag model as a candidate for IGRF-13

Earth Planets Space

(

163

10.1186/s40623-020-01295-y

Bai

Sturrock

P.A.

1991

The 154-day and related periodicities of solar activity as subharmonics of a fundamental period

Nature

350

(

6314

141

–

143

Balasis

Egbert

G.D.

2006

Empirical orthogonal function analysis of magnetic observatory data: Further evidence for non‐axisymmetric magnetospheric sources for satellite induction studies

Geophys. Res. Lett.

2006GL025721

doi:10.1029/2006GL025721

Bartels

1932

Terrestrial-magnetic activity and its relations to solar phenomena

Terr. Magnet. Atmos. Electr.

(

–

10.1029/TE037i001p00001

Bloxham

Jackson

1989

Simultaneous stochastic inversion for geomagnetic main field and secular variation: 2. 1820–1980

J. geophys. Res.: Solid Earth

(

B11

15753

–

15769

10.1029/JB094iB11p15753

Cane

H.V.

Richardson

I.G.

von Rosenvinge

T.T.

1998

Interplanetary magnetic field periodicity of |$\sim$|153 days

Geophys. Res. Lett.

(

4437

–

4440

Chapman

1919

The solar and lunar diurnal variations of terrestrial magnetism

Phil. Trans. R. Soc. Lond.. Series A, Containing Papers of a Mathematical or Physical Character

218

–

118

Google Scholar

OpenURL Placeholder Text

WorldCat

Chapman

Miller

1940

The statistical determination of lunar daily variations in geomagnetic and meteorological elements

Geophys. Suppl. Mon. Not. Roy. Astron. Soc.

(

649

–

669

Google Scholar

OpenURL Placeholder Text

WorldCat

Chowdhury

Choudhary

D.P.

Gosain

Moon

Y.-J.

2015

Short-term periodicities in interplanetary, geomagnetic and solar phenomena during solar cycle 24

Astrophys. Space Sci.

356

(

–

10.1007/s10509-014-2188-0

Chulliat

Vigneron

Hulot

2016

First results from the Swarm Dedicated Ionospheric Field Inversion chain

Earth Planets Space

(

104

doi:10.1186/s40623-016-0481-6

10.1186/s40623-016-0481-6

Egbert

G.D.

Alken

Maute

Zhang

2021

Modelling diurnal variation magnetic fields due to ionospheric currents

Geophys. J. Int.

225

(

1086

–

1109

Finlay

C.C.

Olsen

Tøffner-Clausen

2015

DTU candidate field models for IGRF-12 and the CHAOS-5 geomagnetic field model International Geomagnetic Reference Field - the twelfth generation

Earth Planets Space

–

Finlay

C.C.

Lesur

Thébault

Vervelidou

Morschhauser

Shore

2017

Challenges handling magnetospheric and ionospheric signals in internal geomagnetic field modelling

Space Sci. Rev.

206

(

157

–

189

10.1007/s11214-016-0285-9

Finlay

C.C.

Kloss

Olsen

Hammer

M.D.

Tøffner-Clausen

Grayver

Kuvshinov

2020

The CHAOS-7 geomagnetic field model and observed changes in the South Atlantic Anomaly

Earth Planets Space

(

156

doi:10.1186/s40623-020-01252-9

10.1186/s40623-020-01252-9

Ganushkina

N.Y.

Liemohn

M.W.

Dubyagin

2018

Current systems in the Earth’s magnetosphere

Rev. Geophys.

(

309

–

332

Gauss

C.F.

1839

Allgemeine Theorie des Erdmagnetismus

, in

Resultate aus den Beobachtungen des magnetischen Vereins im Jahre 1838

, pp.

–

Weidmannsche Buchhandlung

Grayver

2024

Unravelling the electrical conductivity of earth and planets

Surv. Geophys.

187

–

238

Grayver

Finlay

C.C.

Olsen

2024

Magnetic signals from oceanic tides: new satellite observations and applications

Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci.

382

(

2286

20240078

Grayver

A.V.

Munch

F.D.

Kuvshinov

A.V.

Khan

Sabaka

T.J.

Tøffner‐Clausen

2017

Joint inversion of satellite‐detected tidal and magnetospheric signals constrains electrical conductivity and water content of the upper mantle and transition zone

Geophys. Res. Lett.

6074

–

6081

Grayver

A.V.

Kuvshinov

Werthmüller

2021

Time-domain modeling of three-dimensional earth’s and planetary electromagnetic induction effect in ground and satellite observations

J. geophys. Res.: Space Phys.

126

e2020JA028672

doi:10.1029/2020JA028672

Guzavina

Grayver

Kuvshinov

2019

Probing upper mantle electrical conductivity with daily magnetic variations using global-to-local transfer functions

Geophys. J. Int.

219

2125

–

2147

Google Scholar

OpenURL Placeholder Text

WorldCat

INTERMAGNET

et al. .,

2024

Intermagnet Reference Data Set (IRDS) 2020–Definitive Magnetic Observatory Data

GFZ Data Services

Google Scholar

Juusola

Vanhamäki

Viljanen

Smirnov

2020

Induced currents due to 3D ground conductivity play a major role in the interpretation of geomagnetic variations

Ann. Geophys.

983

–

998

10.5194/angeo-38-983-2020

Kelbert

Lucas

G.M.

2020

Modified GIC estimation using 3-D earth conductivity

Space Weather

(

e2020SW002467

doi:10.1029/2020SW002467

Kelbert

Schultz

Egbert

2009

Global electromagnetic induction constraints on transition-zone water content variations

Nature

460

1003

–

1006

Knaack

Stenflo

J.O.

2005

Spherical harmonic decomposition of solar magnetic fields

Astron. Astrophys.

438

(

349

–

363

10.1051/0004-6361:20052765

Kuvshinov

2012

Deep electromagnetic studies from land, sea, and space: progress status in the past 10 years

Surv. Geophys.

169

–

209

10.1007/s10712-011-9118-2

Kuvshinov

Grayver

Tøffner-Clausen

Olsen

2021

Probing 3-D electrical conductivity of the mantle using 6 years of Swarm, CryoSat-2 and observatory magnetic data and exploiting matrix Q-responses approach

Earth Planets Space

(

10.1186/s40623-020-01341-9

Langel

R.A.

Estes

R.H.

Mead

G.D.

Fabiano

E.B.

Lancaster

E.R.

1980

Initial geomagnetic field model from Magsat vector data

Geophys. Res. Lett.

(

793

–

796

10.1029/GL007i010p00793

Langel

R.A.

Estes

R.H.

Mead

G.D.

1982

Some new methods in geomagnetic field modeling applied to the 1960-1980 epoch

J. Geomagn. Geoelectr.

(

327

–

349

Langlais

Mandea

2000

An IGRF candidate main geomagnetic field model for epoch 2000 and a secular variation model for 2000–2005

Earth Planets Space

(

1137

–

1148

Larsson

2007

Electromagnetics from a quasistatic perspective

Am. J. Phys.

(

230

–

239

Laundal

K.M.

et al. ,

2021

Electrojet estimates from mesospheric magnetic field measurements

J. geophys. Res.: Space Phys.

126

(

e2020JA028644

doi:10.1029/2020JA028644

Lesur

Wardinski

Hamoudi

Rother

2010

The second generation of the GFZ Reference Internal Magnetic Model: GRIMM-2

Earth Planets Space

(

doi:10.5047/eps.2010.07.007

10.5047/eps.2010.07.007

Lesur

Heumez

Telali

Lalanne

Soloviev

2017

Estimating error statistics for Chambon-la-Forêt observatory definitive data

Ann. Geophys.

(

939

–

952

10.5194/angeo-35-939-2017

Lockwood

Owens

2024

Reconstruction of Carrington rotation means of open solar flux over the past 154 years

Sol. Phys.

299

(

doi:10.1007/s11207-024-02268-0

10.1007/s11207-024-02268-0

Lou

Y.-Q.

2000

Rossby-type wave-induced periodicities in flare activities and sunspot areas or groups during solar maxima

Astrophys. J.

540

(

1102

doi:10.1086/309387

Lou

Y.-Q.

Wang

Y.-M.

Fan

Wang

J.X.

2003

Periodicities in solar coronal mass ejections

Mon. Not. Roy. Astron. Soc.

345

(

809

–

818

10.1046/j.1365-8711.2003.06993.x

Love

J.J.

Chulliat

2013

An international network of magnetic observatories

EoS, Trans. Am. Geophys. Un.

(

373

–

374

Macmillan

Olsen

2013

Observatory data and the Swarm mission

Earth Planets Space

1355

–

1362

10.5047/eps.2013.07.011

Malin

S. R.C.

1973

Worldwide distribution of geomagnetic tides

Phil. Trans. R. Soc. Lond. Ser. A, Math. Phys. Sci.

274

(

1243

551

–

594

Google Scholar

OpenURL Placeholder Text

WorldCat

Malin

S. R.C.

Chapman

1970

The determination of lunar daily geophysical variations by the Chapman-Miller method

Geophys. J. Int.

(

–

10.1111/j.1365-246X.1970.tb06738.x

Maus

Lühr

2005

Signature of the quiet-time magnetospheric magnetic field and its electromagnetic induction in the rotating Earth

Geophys. J. Int.

162

(

755

–

763

10.1111/j.1365-246X.2005.02691.x

Maus

Weidelt

2004

Separating the magnetospheric disturbance magnetic field into external and transient internal contributions using a 1D conductivity model of the Earth

Geophys. Res. Lett.

(

doi:10.1029/2006GC001521

10.1029/2004GL020232

Google Scholar

OpenURL Placeholder Text

WorldCat

Crossref

Maus

Lühr

Rother

Hemant

Balasis

Ritter

Stolle

2007

Fifth-generation lithospheric magnetic field model from CHAMP satellite measurements

Geochem. Geophys. Geosyst.

(

doi:10.1029/2006GC001521

10.1029/2006GC001521

Google Scholar

OpenURL Placeholder Text

WorldCat

Crossref

Olsen

1998

The electrical conductivity of the mantle beneath Europe derived from C-responses from 3 to 720 hr

Geophys. J. Int.

133

298

–

308

10.1046/j.1365-246X.1998.00503.x

Olsen

1999a

Induction studies with satellite data

Surv. Geophys.

309

–

340

10.1023/A:1006611303582

Olsen

1999b

Long-period (30 days-1 year) electromagnetic sounding and the electrical conductivity of the lower mantle beneath Europe

Geophys. J. Int.

138

179

–

187

10.1046/j.1365-246x.1999.00854.x

Olsen

Sabaka

T.J.

Lowes

2005

New parameterization of external and induced fields in geomagnetic field modeling, and a candidate model for IGRF 2005

Earth Planets Space

doi:10.1186/BF03352512

10.1186/BF03352512

Google Scholar

OpenURL Placeholder Text

WorldCat

Crossref

Olsen

Lühr

Sabaka

T.J.

Mandea

Rother

Tøffner-Clausen

Choi

2006

CHAOS - a model of the earth’s magnetic field derived from CHAMP, Ørsted, and SAC-C magnetic satellite data

Geophys. J. Int.

166

–

10.1111/j.1365-246X.2006.02959.x

Olsen

et al. ,

2013

The Swarm Satellite Constellation Application and Research Facility (SCARF) and Swarm data products

Earth Planets Space

(

1189

–

1200

10.5047/eps.2013.07.001

Olsen

Lühr

Finlay

C.C.

Sabaka

T.J.

Michaelis

Rauberg

Tøffner-Clausen

2014

The CHAOS-4 geomagnetic field model

Geophys. J. Int.

197

(

815

–

827

Olsen

Ravat

Finlay

C.C.

Kother

L.K.

2017

LCS-1: a high-resolution global model of the lithospheric magnetic field derived from CHAMP and Swarm satellite observations

Geophys. J. Int.

211

(

1461

–

1477

Olsen

Albini

Bouffard

Parrinello

Tøffner-Clausen

2020

Magnetic observations from CryoSat-2: calibration and processing of satellite platform magnetometer data

Earth Planets Space

(

doi:10.1186/s40623-020-01171-9

10.1186/s40623-020-01171-9

Price

A.T.

Wilkins

G.A.

1963

New methods for the analysis of geomagnetic fields and their application to the Sq field of 1932-3

Phil. Trans. R. Soc. Lond.. Ser. A, Math. Phys. Sci.

256

(

1066

–

Google Scholar

OpenURL Placeholder Text

WorldCat

Püthe

Kuvshinov

2014

Mapping 3-D mantle electrical conductivity from space: a new 3-D inversion scheme based on analysis of matrix Q-responses

Geophys. J. Int.

197

(

768

–

784

Pulkkinen

Amm

Viljanen

2003

Ionospheric equivalent current distributions determined with the method of spherical elementary current systems

J. geophys. Res.: Space Phys.

108

(

doi:10.1029/2001JA005085

10.1029/2001JA005085

Google Scholar

OpenURL Placeholder Text

WorldCat

Crossref

Pulkkinen

et al. ,

2013

Community-wide validation of geospace model ground magnetic field perturbation predictions to support model transition to operations

Space Weather

(

369

–

385

Richmond

A.D.

1995

Ionospheric electrodynamics using magnetic apex coordinates

J. Geomagn. Geoelectr.

(

191

–

212

Rieger

G.H.

Forrest

D.J.

Kanbach

Reppin

Chupp

E.L.

1984

A 154-day periodicity in the occurrence of hard solar flares?

Nature

312

(

5995

623

–

625

Sabaka

T.J.

Olsen

Langel

R.A.

2002

A comprehensive model of the quiet-time, near-Earth magnetic field: phase 3

Geophys. J. Int.

151

(

–

10.1046/j.1365-246X.2002.01774.x

Sabaka

T.J.

Tøffner-Clausen

Olsen

Finlay

C.C.

2018

A comprehensive model of Earth’s magnetic field determined from 4 years of Swarm satellite observations

Earth Planets Space

(

130

doi:10.1186/s40623-018-0896-3

10.1186/s40623-018-0896-3

Sabaka

T.J.

Tøffner-Clausen

Olsen

Finlay

C.C.

2020

CM6: a comprehensive geomagnetic field model derived from both CHAMP and Swarm satellite observations

Earth Planets Space

(

doi:10.1186/s40623-020-01210-5

10.1186/s40623-020-01210-5

Schmucker

1985

Magnetic and electric fields due to electromagnetic induction by external sources

, in

Geophysics of the Solid Earth, the Moon and the Planets, Vol. 2B of Landolt-Börnstein, New Series - Group V Geophysics

Springer

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Schmucker

1999

A spherical harmonic analysis of solar daily variations in the years 1964-1965: Response estimates and source fields for global induction-I. methods

Geophys. J. Int.

136

439

–

454

10.1046/j.1365-246X.1999.00742.x

Schuster

1889

The diurnal variation of terrestrial magnetism

Phil. Trans. R. Soc. A

180

467

–

518

Google Scholar

OpenURL Placeholder Text

WorldCat

Smith

A.R.A.

Pačes

Santillan

2024

. ESA-VirES/VirES-Python-Client,

Zenodo

doi:10.5281/zenodo.14900299

10.5281/zenodo.14900299

Stolle

Michaelis

Xiong

Rother

Usbeck

Yamazaki

Rauberg

Styp-Rekowski

2021

Observing Earth’s magnetic environment with the GRACE-FO mission

Earth Planets Space

(

doi:10.1186/s40623-021-01364-w

10.1186/s40623-021-01364-w

Sugiura

1964

Hourly values of equatorial Dst for the IGY

Ann. Int. Geophys. Yr.

–

Google Scholar

OpenURL Placeholder Text

WorldCat

Tassev

Velinov

P.I.Y.

Tomova

Mateev

2017

Analysis of extreme solar activity in early September 2017: G4-severe geomagnetic storm (07-08.09) and GLE72 (10.09) in solar minimum

C. R. Acad. Bulg. Sci.

(

1445

–

1456

Google Scholar

OpenURL Placeholder Text

WorldCat

Tsyganenko

N.A.

2019

Secular drift of the auroral ovals: How fast do they actually move?

Geophys. Res. Lett.

3017

–

3023

Tóth

et al. ,

2005

Space Weather Modeling Framework: a new tool for the space science community

J. geophys. Res.: Space Phys.

110

(

A12

doi:10.1029/2005JA011126

10.1029/2005JA011126

Google Scholar

OpenURL Placeholder Text

WorldCat

Crossref

Velímský

Knopp

2021

Lateral variations of electrical conductivity in the lower mantle constrained by Swarm and CryoSat-2 missions

Earth Planets Space

Weidelt

1972

The inverse problem of geomagnetic induction

J. Geophys.

257

–

289

Google Scholar

OpenURL Placeholder Text

WorldCat

Winch

D.E.

1981

Spherical harmonic analysis of geomagnetic tides, 1964-1965

Philos. Trans. R. Soc. Lond.. Ser. A, Math. Phys. Sci.

303

(

1473

–

104

10.1098/rsta.1981.0193

Winch

D.E.

Ivers

D.J.

Turner

J. P.R.

Stening

R.J.

2005

Geomagnetism and Schmidt quasi-normalization

Geophys. J. Int.

160

(

487

–

504

10.1111/j.1365-246X.2004.02472.x

Yamazaki

Maute

2017

Sq and EEJ—a review on the daily variation of the geomagnetic field caused by ionospheric dynamo currents

Space Sci. Rev.

206

(

1-4

299

–

405

10.1007/s11214-016-0282-z

Zenhäusern

Kuvshinov

Guzavina

Maute

2021

Towards probing earth’s upper mantle with daily magnetic field variations: exploring a physics-based parametrization of the source

Earth Planets Space

136

doi:10.1186/s40623-021-01455-8

10.1186/s40623-021-01455-8

Zhang

Xiong

Meng

Yuan

Zhang

2019

The influence of geomagnetic storm of 7-8 September 2017 on the Swarm precise orbit determination

J. geophys. Res.: Space Phys.

124

(

6971

–

6984

APPENDIX A: NON-UNIQUE LEAST-SQUARES ESTIMATION OF A PIECE-WISE CONSTANT FIELD AND OBSERVATORY BIAS FROM OBSERVATORY DATA

Given only observatory data, a piece-wise constant field along with the observatory bias cannot be uniquely determined in the least-squares sense. This has been noted already in Langel et al. (1982). Although not complicated, there is no reference to our knowledge where a clear mathematical proof or physical explanation is given. We thus document our mathematical proof here, followed by a physical explanation.

In the i-th time bin (⁠|$i=1,2,\cdots N_t$| where |$N_t$| is the total number of time bins), the data |$\mathbf {d}_i\in \mathbb {R}^{N_i}$| measured at ground observatories are linearly related to the field coefficients |$\mathbf {x}_i \in \mathbb {R}^{M_0}$| via the design matrix |$\mathbf {G}_i \in \mathbb {R}^{N_i\times M_0}$|⁠, combined with a contribution from the observatory bias

$$\begin{eqnarray} \mathbf {d}_i \simeq \mathbf {G}_i \mathbf {x}_i + \mathbf {b}_i. \end{eqnarray}$$

(A1)

Here |$\mathbf {b}_i \in \mathbb {R}^{N_i}$| represents the observatory bias corresponding to the measurements in this time bin. We use |$N_i$| and |$M_0$| to denote the dimensionality of the data and field coefficients, respectively. The dimensionality of the field coefficients |$M_0$| is assumed to be uniform for all time bins. Since the observatory biases are static biases, the bias term in the i-th time bin is simply sampled from the collective observatory biases according to the measurements available in this time bin, and hence can be related to the vector representing the collective observatory biases |$\mathbf {b} \in \mathbb {R}^{M_b}$| by a sampling matrix |$\mathbf {S}_i \in \mathbb {R}^{N_i\times M_b}$|⁠. Therefore, we have

$$\begin{eqnarray} \mathbf {d}_i \simeq \mathbf {G}_i \mathbf {x}_i + \mathbf {S}_i \mathbf {b} = {\begin{pmatrix}\mathbf {G}_i & \mathbf {S}_i \end{pmatrix}} {\begin{pmatrix}\mathbf {x}_i \\ \mathbf {b} \end{pmatrix}}. \end{eqnarray}$$

(A2)

Since the bias |$\mathbf {b}$| is shared across all time bins, one cannot solve each time window separately, but has to solve the full system, written as

$$\begin{eqnarray} \mathbf {d} = {\begin{pmatrix}\mathbf {d}_1 \\ \vdots \\ \mathbf {d}_{N_t} \end{pmatrix}} \simeq {\begin{pmatrix}\mathbf {G}_1 & & \mathbf {0} & \mathbf {S}_1 \\ & \ddots & & \vdots \\ \mathbf {0} & & \mathbf {G}_{N_t} & \mathbf {S}_{N_t} \\ \end{pmatrix}} {\begin{pmatrix}\mathbf {x}_1 \\ \vdots \\ \mathbf {x}_{N_t} \\ \mathbf {b} \end{pmatrix}} = \widetilde{\mathbf {G}} \widetilde{\mathbf {x}}, \end{eqnarray}$$

(A3)

where |$\widetilde{\mathbf {G}} \in \mathbb {R}^{N\times M}$| is the full design matrix for the problem, with |$N = \sum _{i=1}^{N_t} N_i$| and |$M = {N_t} M_0 + M_b$|⁠. |$\widetilde{\mathbf {x}} \in \mathbb {R}^M$| is the augmented model vector. Solving the least-squares problem above is mathematically equivalent to solving the normal equation

$$\begin{eqnarray} \widetilde{\mathbf {G}}^\top \widetilde{\mathbf {G}} \widetilde{\mathbf {x}} = {\begin{pmatrix}\mathbf {G}_1^\top \mathbf {G}_1 & & \mathbf {0} & \mathbf {G}_1^\top \mathbf {S}_1 \\ & \ddots & & \vdots \\ \mathbf {0} & & \mathbf {G}_{N_t}^\top \mathbf {G}_{N_t} & \mathbf {G}_{N_t}^\top \mathbf {S}_{N_t} \\ \\ \mathbf {S}_1^\top \mathbf {G}_1 & \cdots & \mathbf {S}_{N_t}^\top \mathbf {G}_{N_t} & \sum _i \mathbf {S}_i^\top \mathbf {S}_i \end{pmatrix}} \widetilde{\mathbf {x}} = \widetilde{\mathbf {G}}^\top \mathbf {d}. \end{eqnarray}$$

(A4)

We now want to calculate the rank or the nullity of the matrix |$\widetilde{\mathbf {G}}^\top \widetilde{\mathbf {G}} \in \mathbb {R}^{M\times M}$|⁠. We assume that the submatrices |$\mathbf {G}_i^\top \mathbf {G}_i$| are all invertible (i.e. the field coefficients can be estimated time-bin-wise if there is no bias on the measurements). In this case, we can apply the Guttman rank additivity formula

$$\begin{eqnarray} \mathrm{rank} {\begin{pmatrix}\mathbf {A} & \mathbf {B} \\ \mathbf {C} & \mathbf {D} \end{pmatrix}} = \mathrm{rank}(\mathbf {A}) + \mathrm{rank}(\mathbf {D} - \mathbf {C} \mathbf {A}^{-1} \mathbf {B}) \end{eqnarray}$$

(A5)

iteratively to the |$\widetilde{\mathbf {G}}^\top \widetilde{\mathbf {G}}$| matrix. The rank reduces to

$$\begin{eqnarray} \mathrm{rank}\left(\widetilde{\mathbf {G}}^\top \widetilde{\mathbf {G}}\right) &=& \mathrm{rank}\left(\mathbf {G}_1^\top \mathbf {G}_1\right) \\ &&+ \mathrm{rank} {\begin{pmatrix}\mathbf {G}_2^\top \mathbf {G}_2 & & \mathbf {0} & \mathbf {G}_2^\top \mathbf {S}_2 \\ & \ddots & & \vdots \\ \mathbf {0} & & \mathbf {G}_{N_t}^\top \mathbf {G}_{N_t} & \mathbf {G}_{N_t}^\top \mathbf {S}_{N_t} \\ \\ \mathbf {S}_2^\top \mathbf {G}_2 & \cdots & \mathbf {S}_{N_t}^\top \mathbf {G}_{N_t} & \sum _i \mathbf {S}_i^\top \mathbf {S}_i - \mathbf {S}_1^\top \mathbf {G}_1 \mathbf {G}_1^+ \mathbf {S}_1 \end{pmatrix}} \\ &=& \sum _{i=1}^{N_t} \mathrm{rank}\left(\mathbf {G}_i^\top \mathbf {G}_i\right) + \mathrm{rank} \left(\sum _{i=1}^{N_t} \mathbf {S}_i^\top \mathbf {P}_{\mathbf {G}_i}^\perp \mathbf {S}_i\right) ,\end{eqnarray}$$

(A6)

where |$\mathbf {A}^+ = (\mathbf {A}^\top \mathbf {A})^{-1} \mathbf {A}^\top$| is the Moore–Penrose pseudo-inverse of matrix |$\mathbf {A}$|⁠, and |$\mathbf {P}_{\mathbf {A}}^\perp = \mathbf {I} - \mathbf {A} \mathbf {A}^+$| is the projection operator onto the orthogonal complement of the image of |$\mathbf {A}$|⁠. Under our assumption, |$\mathbf {G}_i$| is full column rank and |$\mathrm{rank}(\mathbf {G}_i^\top \mathbf {G}_i) = M_0$|⁠. Applying the rank-nullity theorem, we have

$$\begin{eqnarray} \mathrm{dim} \, \mathrm{ker}\left(\widetilde{\mathbf {G}}^\top \widetilde{\mathbf {G}}\right) &=& \mathrm{dim}\left(\widetilde{\mathbf {G}}^\top \widetilde{\mathbf {G}}\right) - \mathrm{rank} \left(\widetilde{\mathbf {G}}^\top \widetilde{\mathbf {G}}\right) \\ &=& \mathrm{dim}\left(\widetilde{\mathbf {G}}^\top \widetilde{\mathbf {G}}\right) - \sum _{i=1}^{N_t} \mathrm{rank}\left(\mathbf {G}_i^\top \mathbf {G}_i\right) - \mathrm{dim} \left(\sum _{i=1}^{N_t} \mathbf {S}_i^\top \mathbf {P}_{\mathbf {G}_i}^\perp \mathbf {S}_i\right) + \mathrm{dim}\, \mathrm{ker} \left(\sum _{i=1}^{N_t} \mathbf {S}_i^\top \mathbf {P}_{\mathbf {G}_i}^\perp \mathbf {S}_i\right)\\ &=& M - {N_t} M_0 - M_b + \mathrm{dim}\, \mathrm{ker} \left(\sum _{i=1}^{N_t} \mathbf {S}_i^\top \mathbf {P}_{\mathbf {G}_i}^\perp \mathbf {S}_i\right) = \mathrm{dim}\, \mathrm{ker} \left(\sum _{i=1}^{N_t} \mathbf {S}_i^\top \mathbf {P}_{\mathbf {G}_i}^\perp \mathbf {S}_i\right) . \end{eqnarray}$$

(A7)

Therefore, the dimension of the nullspace of the original linear system can be identified with that of the matrix |$\sum _i \mathbf {S}_i^\top \mathbf {P}_{\mathbf {G}_i}^\perp \mathbf {S}_i$|⁠. In a piece-wise constant field model with observatory data, not only does |$\mathbf {S}_i$| sample the biases in each time bin, but the same sampling matrix also samples the design matrix. In other words, the design matrices in each time bin are sampled from the rows of a design matrix |$\mathbf {G}_0 \in \mathbb {R}^{M_b\times M_0}$|⁠, whose rows contain the collective design vectors for all observatories, hence |$\mathbf {G}_i = \mathbf {S}_i \mathbf {G}_0$|⁠. The whole column space of |$\mathbf {G}_0$| hence maps to zero under the matrix |$\sum _i \mathbf {S}_i^\top \mathbf {P}_{\mathbf {G}_i}^\perp \mathbf {S}_i$|

$$\begin{eqnarray} \left(\sum _{i=1}^{N_t} \mathbf {S}_i^\top \mathbf {P}_{\mathbf {G}_i}^\perp \mathbf {S}_i\right) \mathbf {G}_0 = \sum _{i=1}^{N_t} \mathbf {S}_i^\top \mathbf {P}_{\mathbf {G}_i}^\perp (\mathbf {S}_i \mathbf {G}_0) = \sum _{i=1}^{N_t} \mathbf {S}_i^\top \mathbf {P}_{\mathbf {G}_i}^\perp \mathbf {G}_i = \mathbf {0} \end{eqnarray}$$

(A8)

since |$\mathbf {P}_{\mathbf {A}}^\perp \mathbf {A}\equiv \mathbf {0}$|⁠, |$\forall \mathbf {A}$|⁠. Therefore, we arrive at

$$\begin{eqnarray} \mathrm{dim} \, \mathrm{ker}\left(\widetilde{\mathbf {G}}^\top \widetilde{\mathbf {G}}\right) = \mathrm{dim}\, \mathrm{ker} \left(\sum _{i=1}^{N_t} \mathbf {S}_i^\top \mathbf {P}_{\mathbf {G}_i}^\perp \mathbf {S}_i\right) \ge \mathrm{rank}(\mathbf {G}_0) = M_0. \\ \end{eqnarray}$$

(A9)

The full linear system thus is singular, with a nullspace of dimension |$M_0$|⁠, the same dimension as the field coefficients. There is no unique least-squares solution to this problem.

From a physical point of view, introducing the observatory bias in the regular spherical harmonic analysis allows one to embed an arbitrary time-invariant field that can be represented with the field coefficients in the bias. If a set of field coefficients solves the original least-squares problem, adding an arbitrary constant vector |$\delta \mathbf {x}$| to the field coefficients in all time bins uniformly also solves the problem, as long as the observatory biases are adjusted accordingly to compensate for the static field. Therefore there is no uniqueness in the solution.

APPENDIX B: NON-UNIQUE LEAST-SQUARES ESTIMATION OF A PIECE-WISE CONSTANT FIELD AND OBSERVATORY BIAS FROM GROUND AND SATELLITE DATA

The observatory bias along with the field coefficients can be uniquely determined in a least-squares sense if both satellite and ground observatory data are used, provided that the two data sets share the same internal and external contributions in the geomagnetic field parametrization. However, since we introduced an extra ionospheric field in our complete model which contributes to the fields measured by observatories and satellites differently, the solution remains non-unique. Here we also provide a mathematical and a physical explanation for this non-uniqueness.

The least-squares problem in this scenario can be formulated in the same way as in the last section. However, the setup of the submatrices are different from the previous section due to the different data set and field parametrization. In this scenario, we have both the data from ground observatories, and the data from satellites, which are unaffected by the observatory bias. Therefore, the data vectors, the design matrices and the sampling matrices have row divisions

$$\begin{eqnarray} \mathbf {d}_i = {\begin{pmatrix}\mathbf {d}_{i,\mathrm{obs}} \\ \mathbf {d}_{i,\mathrm{sat}} \end{pmatrix}}, \quad \mathbf {G}_i = {\begin{pmatrix}\mathbf {G}_{i,\mathrm{obs}} \\ \mathbf {G}_{i,\mathrm{sat}} \end{pmatrix}}, \quad \mathbf {S}_i = {\begin{pmatrix}\mathbf {S}^{\prime }_{i} \\ \mathbf {0} \end{pmatrix}}. \end{eqnarray}$$

(B1)

On the other hand, in our complete model, the field coefficients are separated into internal, ionospheric and magnetospheric parts. The field coefficient vector and the design matrix has the following row and column divisions,

$$\begin{eqnarray} \mathbf {x}_i = {\begin{pmatrix}\mathbf {x}_i^{\mathrm{int}} \\ \mathbf {x}_i^{\mathrm{ion}} \\ \mathbf {x}_i^{\mathrm{mag}} \end{pmatrix}}, \quad \mathbf {G}_i = {\begin{pmatrix}\mathbf {G}_{i,\mathrm{obs}} \\ \mathbf {G}_{i,\mathrm{sat}} \end{pmatrix}} = {\begin{pmatrix}\mathbf {G}_{i,\mathrm{obs}}^{\mathrm{int}} & \mathbf {G}_{i,\mathrm{obs}}^{\mathrm{ion}} & \mathbf {G}_{i,\mathrm{obs}}^{\mathrm{mag}} \\ \\ \mathbf {G}_{i,\mathrm{sat}}^{\mathrm{int}} & \mathbf {G}_{i,\mathrm{sat}}^{\mathrm{ion}} & \mathbf {G}_{i,\mathrm{sat}}^{\mathrm{mag}} \end{pmatrix}}. \end{eqnarray}$$

(B2)

For simplicity, we shall only discuss the case where the three fields are parametrized in the same set of spherical harmonics, hence |$\mathbf {x}_i^{\mathrm{int}}, \mathbf {x}_i^{\mathrm{ion}}, \mathbf {x}_i^{\mathrm{mag}} \in \mathbb {R}^{M_0/3}$|⁠. A similar conclusion holds for a non-uniform truncation of the fields, but the notation will be much more complicated.

Following eq. (13), we have |$\mathbf {G}_{i,\mathrm{obs}}^{\mathrm{ion}} = \mathbf {G}_{i,\mathrm{obs}}^{\mathrm{mag}} = \mathbf {G}_{i,\mathrm{obs}}^{\mathrm{ext}}$|⁠, that is, the design matrix linking the ionospheric and magnetospheric coefficients to the ground observations are the same, and are simply the design matrix for external coefficients |$\mathbf {G}_{i,\mathrm{obs}}^{\mathrm{ext}}$| in conventional internal–external field separation. These submatrices corresponding to the ground observations are also sampled by the sampling matrix from the conglomerate design matrices, that is,

$$\begin{eqnarray} \mathbf {G}_{i,\mathrm{obs}}^\mathrm{int} = \mathbf {S}_i^{\prime } \mathbf {G}_{0,\mathrm{obs}}^\mathrm{int},\quad \mathbf {G}_{i,\mathrm{obs}}^\mathrm{ext} = \mathbf {S}_i^{\prime } \mathbf {G}_{0,\mathrm{obs}}^\mathrm{ext}. \end{eqnarray}$$

(B3)

On the other hand, following eq. (14), we have |$\mathbf {G}_{i,\mathrm{sat}}^\mathrm{mag} = \mathbf {G}_{i,\mathrm{sat}}^\mathrm{ext}$|⁠, but |$\mathbf {G}_{i,\mathrm{sat}}^\mathrm{ion} = \mathbf {G}_{i,\mathrm{sat}}^\mathrm{int} \boldsymbol {\Lambda }$|⁠. |$\boldsymbol {\Lambda }$| is diagonal matrix responsible for changing the internal coefficients with ionospheric origins observed at satellite altitudes into the external coefficients with ionospheric origins observed on the ground, as described by eq. (12). In the end, the submatrix |$\mathbf {G}_i$| takes the form

$$\begin{eqnarray} \mathbf {G}_i = {\begin{pmatrix}\mathbf {S}_i^{\prime } \mathbf {G}_{0,\mathrm{obs}}^{\mathrm{int}} & \mathbf {S}_i^{\prime } \mathbf {G}_{0,\mathrm{obs}}^{\mathrm{ext}} & \mathbf {S}_i^{\prime } \mathbf {G}_{0,\mathrm{obs}}^{\mathrm{ext}} \\ \\ \mathbf {G}_{k,\mathrm{sat}}^{\mathrm{int}} & \mathbf {G}_{k,\mathrm{sat}}^{\mathrm{int}} \boldsymbol {\Lambda } & \mathbf {G}_{k,\mathrm{sat}}^{\mathrm{ext}} \end{pmatrix}}. \end{eqnarray}$$

(B4)

The mathematical formulation of this proof follows the same line as that in Appendix A. The dimension of the nullspace of the original linear system can still be identified with the that of the matrix |$\sum _i \mathbf {S}_i^\top \mathbf {P}_{\mathbf {G}_i}^\perp \mathbf {S}_i$| (eq. A7). The underlying assumption here is that the field coefficients can still be determined uniquely in a least-squares sense in the absence of observatory bias, that is, |$\mathbf {G}_i^\top \mathbf {G}_i$| invertible. A corollary is that |$\mathbf {G}_{0,\mathrm{obs}}^\mathrm{int}$| and |$\mathbf {G}_{0,\mathrm{obs}}^\mathrm{ext}$| must both be full column rank, and their column spaces should be mutually linearly independent. Therefore, we can introduce a full column rank matrix |$\mathbf {G}^{\prime } = \mathbf {G}_{0,\mathrm{obs}}^\mathrm{int} \boldsymbol {\Lambda } - \mathbf {G}_{0,\mathrm{obs}}^\mathrm{ext}$|⁠. The action of the sampling matrix |$\mathbf {S}_i$| on |$\mathbf {G}^{\prime }$| lies within the image of |$\mathbf {G}_i$|⁠, as

$$\begin{eqnarray} \mathbf {S}_i \mathbf {G}^{\prime } = {\begin{pmatrix}\mathbf {S}_i^{\prime } \\ \mathbf {0} \end{pmatrix}} ( \mathbf {G}_{0,\mathrm{obs}}^\mathrm{int} \boldsymbol {\Lambda } - \mathbf {G}_{0,\mathrm{obs}}^\mathrm{ext}) = {\begin{pmatrix}\mathbf {S}_i^{\prime } \mathbf {G}_{0,\mathrm{obs}}^\mathrm{int} \boldsymbol {\Lambda } - \mathbf {S}^{\prime }_i \mathbf {G}_{0,\mathrm{obs}}^\mathrm{ext} \\ \mathbf {0} \end{pmatrix}} = \mathbf {G}_i {\begin{pmatrix}\boldsymbol {\Lambda } \\ - \mathbf {I} \\ \mathbf {0} \end{pmatrix}} ,\end{eqnarray}$$

(B5)

where we used eq. (B4) in the last step. As a result, the image of |$\mathbf {G}^{\prime }$| lies in the nullspace of |$\sum _i \mathbf {S}_i^\top \mathbf {P}_{\mathbf {G}_i}^\perp \mathbf {S}_i$|⁠,

$$\begin{eqnarray} \left(\sum _{i=1}^{N_t} \mathbf {S}_i^\top \mathbf {P}_{\mathbf {G}_i}^\perp \mathbf {S}_i\right) \mathbf {G}^{\prime } = \sum _{i=1}^{N_t} \mathbf {S}_i^\top \mathbf {P}_{\mathbf {G}_i}^\perp \mathbf {G}_i {\begin{pmatrix}\boldsymbol {\Lambda } \\ - \mathbf {I} \\ \mathbf {0} \end{pmatrix}} = \mathbf {0}. \end{eqnarray}$$

(B6)

The nullity of the full system can be calculated via

(B7)

which concludes our proof that the problem has no unique least-squares solution.

From a physical point of view, when the complete model is augmented with observatory bias, an arbitrary time-invariant field can be embedded in the induced field, which can always be compensated by the ionospheric field and the observatory bias. If a set of field coefficients solves the original least-squares problem, we can add an arbitrary vector |$\delta \mathbf {x}^\mathrm{int}$| to the induced coefficients in every time bin. An additional ionospheric field vector |$\delta \mathbf {x}^\mathrm{ion} = - \boldsymbol {\Lambda }^{-1} \delta \mathbf {x}^\mathrm{int}$| can be added so that they cancel out in the satellite data. The residual in the observatory data can be compensated by adjusting the observatory biases accordingly. The adjusted model vectors will provide the exact same data fit as the original model vectors, and hence the solution is non-unique.

APPENDIX C: SOLVING BLOCK-DIAGONAL LINEAR SYSTEM WITH INTERCEPT

In spherical harmonic analysis of a piece-wise constant geomagnetic field given vector magnetic field measurements, the least-squares problem that arises takes the form

$$\begin{eqnarray} {\begin{pmatrix}\mathbf {G}_1 & & \mathbf {0} \\ & \ddots & \\ \mathbf {0} & & \mathbf {G}_{N_t} \end{pmatrix}} {\begin{pmatrix}\mathbf {x}_1 \\ \vdots \\ \mathbf {x}_{N_t} \end{pmatrix}} \simeq {\begin{pmatrix}\mathbf {d}_1 \\ \vdots \\ \mathbf {d}_{N_t} \end{pmatrix}} , \end{eqnarray}$$

(C1)

which we denote using the shorthand notation |$\mathbf {G} \mathbf {x} \simeq \mathbf {d}$|⁠. Since the matrix |$\mathbf {G}$| is block diagonal, the field coefficients for different time bins naturally decouple, and the least-squares problem can be solved for each time bin separately. The computational and memory cost is therefore bounded by the summation of the costs for solving subsystems |$\mathbf {G}_i \mathbf {x}_i \simeq \mathbf {d}_i$|⁠, without involving the action of the full matrix |$\mathbf {G}$|⁠.

In our estimation of the observatory bias, the least-squares problem that arises takes the form

$$\begin{eqnarray} {\begin{pmatrix}\mathbf {G} & \mathbf {S} \end{pmatrix}} {\begin{pmatrix}\mathbf {x} \\ \mathbf {b} \end{pmatrix}} = {\begin{pmatrix}\mathbf {G}_1 & & \mathbf {0} & \mathbf {S}_1 \\ & \ddots & & \vdots \\ \mathbf {0} & & \mathbf {G}_{N_t} & \mathbf {S}_{N_t} \\ \end{pmatrix}} {\begin{pmatrix}\mathbf {x}_1 \\ \vdots \\ \mathbf {x}_{N_t} \\ \mathbf {b} \end{pmatrix}} \simeq {\begin{pmatrix}\mathbf {d}_1 \\ \vdots \\ \mathbf {d}_{N_t} \end{pmatrix}}. \end{eqnarray}$$

(C2)

This linear system has been described in Appendix A, and is denoted |$\widetilde{\mathbf {G}} \widetilde{\mathbf {x}} \simeq \mathbf {d}$|⁠. Although the left part of the matrix is still block diagonal, the right submatrix couples all the subsystems together by sampling the contribution from the observatory biases in every subsystem, and prohibits solving each subsystem separately. Naive implementation might result in solving the full system, incurring large computational and memory costs.

We present here an algorithm that reduces the full system into solving a series of smaller subsystems by virtue of matrix partitioning. The algorithm consists of two steps. In the first step, the bias |$\mathbf {b}$| is estimated. After removing the bias contribution from the data, the |$\mathbf {x}_i$| vectors can then be estimated using the block-diagonal linear least-squares solver in the second step. To begin with, solving the aforementioned linear least squares is equivalent mathematically to solving the normal equation

$$\begin{eqnarray} \widetilde{\mathbf {G}}^\top \widetilde{\mathbf {G}} \widetilde{\mathbf {x}} = {\begin{pmatrix}\mathbf {G}^\top \mathbf {G} & \mathbf {G}^\top \mathbf {S} \\ \mathbf {S}^\top \mathbf {G} & \mathbf {S}^\top \mathbf {S} \end{pmatrix}} {\begin{pmatrix}\mathbf {x} \\ \mathbf {b} \end{pmatrix}} = {\begin{pmatrix}\mathbf {G}^\top \mathbf {d} \\ \mathbf {S}^\top \mathbf {d} \end{pmatrix}} = \widetilde{\mathbf {G}}^\top \mathbf {d}. \end{eqnarray}$$

(C3)

The solution can be formally written as |$\widetilde{\mathbf {x}} = (\widetilde{\mathbf {G}}^\top \widetilde{\mathbf {G}})^{-1} \widetilde{\mathbf {G}}^\top \mathbf {d}$|⁠, where we assume that the matrix |$\widetilde{\mathbf {G}}^\top \widetilde{\mathbf {G}}$| is invertible. Assuming |$\mathbf {G}^\top \mathbf {G}$| is also invertible, we can use the first row block to express |$\mathbf {x}$| as

$$\begin{eqnarray} \mathbf {x} = (\mathbf {G}^\top \mathbf {G})^{-1} \mathbf {G}^\top (\mathbf {d} - \mathbf {S} \mathbf {b}) = \mathbf {G}^+ (\mathbf {d} - \mathbf {S} \mathbf {b}) , \end{eqnarray}$$

(C4)

which when inserted into the second row block, yields

$$\begin{eqnarray} {\begin{array}{c}\left(\mathbf {S}^\top \mathbf {S} - \mathbf {S}^\top \mathbf {G} \mathbf {G}^+ \mathbf {S}\right)\mathbf {b} = \mathbf {S}^\top \mathbf {d} - \mathbf {S}^\top \mathbf {G} \mathbf {G}^+ \mathbf {d} \\ \left(\mathbf {S}^\top \mathbf {P}_{\mathbf {G}}^\perp \mathbf {S}\right) \mathbf {b} = \mathbf {S}^\top \mathbf {P}_{\mathbf {G}}^\perp \mathbf {d}. \end{array}} \end{eqnarray}$$

(C5)

It can be shown using Guttman rank additivity theorem that if |$\widetilde{\mathbf {G}}^\top \widetilde{\mathbf {G}}$| and |$\mathbf {G}^\top \mathbf {G}$| are both invertible, then |$\mathbf {S}^\top \mathbf {P}_{\mathbf {G}}^\perp \mathbf {S}$| must be invertible as well. Therefore, we arrive at the solution of |$\mathbf {b}$|

$$\begin{eqnarray} \mathbf {b} = \left(\mathbf {S}^\top \mathbf {P}_{\mathbf {G}}^\perp \mathbf {S}\right)^{-1} \mathbf {S}^\top \mathbf {P}_{\mathbf {G}}^\perp \mathbf {d}. \end{eqnarray}$$

(C6)

Thanks to the block-diagonality of the |$\mathbf {G}$| matrix, the matrix–matrix and matrix–vector multiplications in the equation above can be drastically simplified. The seemingly high-dimensional multiplications turn out to be merely a summation of low-dimensional products; the solution of |$\mathbf {b}$| therefore reads

$$\begin{eqnarray} \mathbf {b} = \left(\sum _{i=1}^{N_t} \mathbf {S}_i^\top \mathbf {P}_{\mathbf {G}_i}^\perp \mathbf {S}_i\right)^{-1} \sum _{i=1}^{N_t} \mathbf {S}_i^\top \mathbf {P}_{\mathbf {G}_i}^\perp \mathbf {d}_i. \end{eqnarray}$$

(C7)

The bias estimation problem is reduced to collecting the matrix–matrix and matrix–vector products in each subsystem, and then performing one low-dimensional linear solve. This completes the first step.

In the second step, the contribution from |$\mathbf {b}$| is removed from each subsystem. The remaining subsystems can be solved separately as they are now decoupled.

This algorithm does not require the action or the decomposition of the full matrix, but merely requires matrix multiplications of low-dimensional submatrices. As such, the algorithm provides a computation- and memory-efficient way to solve such block-diagonal systems with bias contribution.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Download all slides

Month:	Total Views:
February 2025	48
March 2025	190
April 2025	113

Article Contents

A decade of the fast-varying ionospheric and magnetospheric magnetic fields from ground and multisatellite observations

SUMMARY

1 INTRODUCTION

2 METHODOLOGY

3 DATA

3.1 Satellite data

3.2 Ground observatory data

3.2.1 Observatory bias

4 MODEL CONSTRUCTION AND SELECTION

5 RESULTS

5.1 The role of ionosphere in the model

5.2 Time-series of external and induced coefficients

5.3 Fourier spectrum of external and induced fields

6 DISCUSSION AND CONCLUSIONS

AUTHOR CONTRIBUTIONS

ACKNOWLEDGMENTS

DATA AVAILABILITY

Footnotes

REFERENCES

APPENDIX A: NON-UNIQUE LEAST-SQUARES ESTIMATION OF A PIECE-WISE CONSTANT FIELD AND OBSERVATORY BIAS FROM OBSERVATORY DATA

APPENDIX B: NON-UNIQUE LEAST-SQUARES ESTIMATION OF A PIECE-WISE CONSTANT FIELD AND OBSERVATORY BIAS FROM GROUND AND SATELLITE DATA

APPENDIX C: SOLVING BLOCK-DIAGONAL LINEAR SYSTEM WITH INTERCEPT

Supplementary data

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

Article Contents

A decade of the fast-varying ionospheric and magnetospheric magnetic fields from ground and multisatellite observations

SUMMARY

1 INTRODUCTION

2 METHODOLOGY

3 DATA

3.1 Satellite data

3.2 Ground observatory data

3.2.1 Observatory bias

4 MODEL CONSTRUCTION AND SELECTION

5 RESULTS

5.1 The role of ionosphere in the model

5.2 Time-series of external and induced coefficients

5.3 Fourier spectrum of external and induced fields

6 DISCUSSION AND CONCLUSIONS

AUTHOR CONTRIBUTIONS

ACKNOWLEDGMENTS

DATA AVAILABILITY

Footnotes

REFERENCES

APPENDIX A: NON-UNIQUE LEAST-SQUARES ESTIMATION OF A PIECE-WISE CONSTANT FIELD AND OBSERVATORY BIAS FROM OBSERVATORY DATA

APPENDIX B: NON-UNIQUE LEAST-SQUARES ESTIMATION OF A PIECE-WISE CONSTANT FIELD AND OBSERVATORY BIAS FROM GROUND AND SATELLITE DATA

APPENDIX C: SOLVING BLOCK-DIAGONAL LINEAR SYSTEM WITH INTERCEPT

Supplementary data

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only