Introduction

The main part of the Earth’s magnetic field is generated by motions in the electrical conducting liquid outer core, in a process known as the geodynamo. This magnetic field, termed the core field, exhibits both spatial and temporal changes over a broad range of scales. Magnetic measurements from satellites have increased the recovery of small-scale features of this field and revealed rapid changes in its temporal behavior (e.g., Alken et al. 2020a; Finlay et al. 2020; Baerenzung et al. 2020; Ropp et al. 2020; Sabaka et al. 2020). From ground and space magnetic observations, variations in the first and second time derivatives of the field, termed the secular variation (SV) and acceleration (SA), respectively, may now be resolved down to periods of about 1 to 2 years (Lesur et al. 2010, 2017; Ropp et al. 2020).

Studies have revealed oscillating SA pulse-like field features at the core–mantle boundary (CMB) focused in the region around the geographical equator (Chulliat et al. 2010; Chulliat and Maus 2014; Chulliat et al. 2015; Sabaka et al. 2018; Alken et al. 2020a). The interpretation and geophysical mechanisms responsible for driving such distinctive behavior in the SA signal remains under debate (e.g., Gillet 2019; Buffett and Matsui 2019; Aubert and Finlay 2019; Gerick et al. 2020), as is the connection to abrupt changes in the SV observed at ground observatories (Mandea et al. 2010). The secular acceleration must be characterized with care, paying attention to those spatial and temporal scales that are well resolved, as its observed spatial spectra at the core surface is blue, showing increasing power with spherical harmonic degree, and its observed temporal spectra seems to be rather flat, meaning that there could be important unresolved fast variations (Christensen et al. 2012; Bouligand et al. 2016; Lesur et al. 2017; Gillet 2019). In this respect, assessing the limitations of the information obtained from measurements by analyzing their resolving power is crucial when aiming to investigate the SA signal at small length scales and short timescales.

Magnetic field measurements from low-Earth orbiting (LEO) satellites provide global field observations, which have proved important for mapping the spatial structures of the core field signal (e.g. Olsen and Stolle 2012). Beginning with the launch of the Danish Ørsted (1999–2014) satellite, the German CHAMP (2000–2010) satellite and the European Swarm (2013-) satellite trio, satellites have provided high quality magnetic field measurements, enabling global investigations into the spatio-temporal variations of field. Unfortunately, the CHAMP mission ended in September 2010 and since reliable vector measurements from the Ørsted satellite extend only up to 2006, there is a gap from 2010 and 2014 in the satellite magnetic records (Finlay et al. 2016). However, other satellite missions, not dedicated to measuring the magnetic field, offer a possibility to fill in this gap adding information about the field. In particular, the CryoSat-2 (2010-) mission, intended for measuring polar ice thickness, carries three platform magnetometers for navigational purposes. Calibrated CryoSat-2 measurements, where vector fluxgate magnetometer readings have been transformed into reliable magnetic field vector outputs, from August 2010 to December 2018 have recently become available (Olsen et al. 2020), such that there are now in total 20 years of continuous satellite measurements. Regarding the data from CryoSat-2, it is clearly important to assess the quality and limitations of the calibrated platform magnetometer measurements and to test what contribution they can make to the study of core field variations.

The standard approach to using satellite magnetometer data for core field studies is to construct spherical harmonic (SH) field models by least-squares inversion methods (see e.g., Langel 1987). In such global models B-splines are often used to parameterize the model time-dependence. Typically this necessitates temporal regularization which modifies the time-dependence of the harmonics in a non-uniform manner. An undesirable consequence is that, for the higher harmonics, the first time derivative effectively becomes a time average over an increasingly long interval, rather than an estimate of the instantaneous secular variation (Olsen et al. 2009). CryoSat-2 data has already been used in the construction of such time-dependent spherical harmonic field models by Alken et al. (2020a), Finlay et al. (2020) and Kloss et al. (2021). From the viewpoint of accessing more detailed information on the spatial and temporal structure of the field, it is also of interest to look at alternative techniques for studying secular variation that can complement the traditional SH approach. In this paper, we focus on two local methods for studying core field variations as recorded in satellite measurements, with a focus on assessing the quality and resolving ability of CryoSat-2 magnetic data.

In a first assessment of the quality of the CryoSat-2 data and its ability to map the SV field at satellite altitude together with CHAMP and Swarm we use the Geomagnetic Virtual Observatory (GVO) technique proposed by Mandea and Olsen (2006) and Olsen and Mandea (2007). This technique involves computing time series of field estimates at specified target locations at satellite altitude, from satellite measurements taken nearby. We apply the processing algorithm recently developed to derive 4-monthly Swarm GVO data series (Hammer et al. 2021) and derive time series on a global grid of 300 GVOs. This network of GVOs allows field changes at satellite altitude to be investigated globally at fixed locations. The GVOs provide a useful compression of satellite magnetic measurements and are a convenient dataset for workers wishing to use constraints from satellite measurements for studies of core dynamics. The GVO dataset involves a series of independent local constraints that can be separately assessed rather than the inherently global constraints provided by SH models. In addition the GVOs have well understood error covariances that can be assigned by methods similar to those used with ground observatories. GVOs have already been used by a number of groups for studies of core dynamics (e.g., Whaler and Beggan 2015; Barrois et al. 2018; Domingos et al. 2019).

In the second part of this study, we directly map the radial field SV at the core–mantle boundary, using the technique of Subtractive Optimally Localized Averages (SOLA) that was adapted to geomagnetism by Hammer and Finlay (2019). The SOLA technique can be used to compute estimates of the radial field SV directly at the CMB, based on local spatial averages the SV field centered on target locations of interest and time-averages over chosen time windows. By collecting many individual SOLA estimates on a grid at the CMB, the SV field can be mapped on regional or global scales. An important foundation of the SOLA technique is that, for noise-free data and a linear forward problem, any linear combination of the data provides a specific average of the true model. With noisy data, a variance is ascribed to this spatial average value such that a trade-off between resolution and variance arises (Oldenburg 1984; Parker 1994). The SOLA technique readily provides information on the resolution offered by a given set of magnetic field observations, in the form of averaging kernels, as well as estimates of the variance of the locally averaged field. We compare SOLA-based maps of SV and SA estimates derived from CryoSat-2 and Swarm data, in order to asses the suitability of the CryoSat-2 data for mapping these fields at the CMB. Demonstrating the usefulness of the CryoSat-2 data, we then take advantage of this data to map the time evolution of the SA field along the geographic equator at the CMB from 2001 to 2019.

We wish to emphasize that both the GVO and the SOLA methods can result in patterns of secular variation different from those seen in the CHAOS field model, despite the fact they use similar data selection schemes and processing steps that involve the same magnetospheric field model. In the GVO and SOLA methods data close to a location of interest are effectively used to determine localized field or SV estimates. In contrast, estimation of the coefficients of truncated spherical harmonic expansions, in models such as CHAOS, is an inherently global procedure that aims at finding the best possible global model. In the GVO method only data from within a 700-km radius of a target location is used to determine the local potential. In the SOLA method measurements far from the target location have essentially no influence on the estimated CMB SV because the data kernels, based on Green’s functions for Laplace’s equation under Neumann boundary conditions, have decreasing sensitivity far away from the target location.

Moreover, the CHAOS model involves temporal regularization whereby one minimizes global norms based on time derivatives of the CMB radial field, integrated over the entire timespan of the model. In the GVO method there is no temporal smoothing beyond the choice of 4-month data windows and use of annual differences to produce SV series. In the SOLA method we use 1- or 2-year windows to estimate the SV and then annual differences to estimate SA. The SOLA method involves a trade-off parameter specifying the balance between the spatial resolution and the variance of each local estimate, each estimate being time-averaged over a 1- or 2-year time window; there is no global regularization over longer time spans. Both methods presented here therefore constitute a localized compression of information contained within the satellite data on the potential field near the location of interest within the specified time window. The GVOs and SOLA can thus give a different picture of the core field evolution compared to the spherical-harmonic based field model that involve global temporal regularization, especially regarding rapid changes at short wavelengths when the data quality is high (see Sect. 4.2).

Section 2 describes the satellite measurements used in this study, including the CryoSat-2 platform magnetometer data, including how these have been selected and processed. Section 3 presents the GVO technique and results concerning global GVO time series. Section 4 presents the SOLA technique and results of applying this to estimate SV and SA at the CMB. Conclusions and perspectives are given in Sect. 5.

Data selection

We use satellite vector magnetic field measurements from the Ørsted satellite between July 2000 and December 2005, from CHAMP taking L3 magnetic data between July 2000 and September 2010, and from the Swarm trio taking Level 1b MAG-L data, version 0505/0506, between January 2014 and April 2020. Most importantly for this study we make use of platform magnetometer data from the CryoSat-2 mission, taking calibrated vector measurements with a sampling rate of 4s from the FGM1 magnetometer dataset, version 3, between August 2010 to December 2018. This dataset has had extensive corrections applied for disturbances fields, and was calibrated using a reference field model—for full details see Olsen et al. (2020).

From the Ørsted, CHAMP and Swarm measurements we produced two data sets: dataset #1 used in the GVO application taking a 15 s subsampling of the vector field measurements from Ørsted, CHAMP and Swarm, while taking every 4th measurement from the CryoSat-2 dataset (i.e., a 16-s subsampling); dataset #2 used in the SOLA application taking a 5-s subsampling of the vector field measurements from Ørsted, CHAMP and Swarm, while taking every element of the CryoSat-2 dataset with its 4 s sampling rate. Field measurements having gross data outliers for which the vector field components deviated more than 500 nT from the CHAOS-7.2 internal field model predictions (Finlay et al. 2020) were rejected. For both data sets we apply a dark geomagnetically quiet-time selection criteria scheme. See Table 1 for full details of the selection requirements for the datasets used in the GVO and SOLA applications. In both cases we required the sun to be at least \(10^{\circ }\) below the horizon, adding restrictions on the geomagnetic activity index \(K_p\) and the change in the ring current index (see Olsen et al. 2014), as well as constraints on the merging electric field at magnetopause, and on the magnitudes of the \(B_Y\) and \(B_Z\) components of the interplanetary magnetic field (e.g., Finlay et al. 2020; Ritter et al. 2004). We used minute values of the IMF components and solar wind speed from the OMNI database, http://omniweb.gsfc.nasa.gov, computing two-hourly means prior to the time of the considered datum (Finlay et al. 2016).

Since our focus is the core field, we have applied corrections to the data for the lithospheric and external fields. For the lithospheric field model, we used the LCS-1 model (Olsen et al. 2017); the precise choice of lithospheric field is not crucial for studies of the SV and SA. For the solar-quiet ionospheric field and associated induced fields we used the CIY4 model (Sabaka et al. 2018). For the magnetospheric field and related induced fields we used the CHAOS-7 model (Finlay et al. 2020). These models were chosen as they are well established and compatible with the data selection criteria described above.

Table 1 Selection criteria, and model corrections applied for the GVO and SOLA datasets used in this study
Full size table

As noted in Table 1, for the GVO application we use the sums and differences of the magnetic field measurements. It has been shown by Olsen (2015) and Sabaka et al. (2018), that taking differences of the satellite measurements along-track and East–West (between Swarm satellites A and C) helps the recovery of the small-scale core field, as this reduces the impact of correlated errors caused by unmodeled large-scale external fields. Here we follow such an approach by taking differences of the measurements, but we also include along-track and East–West sums of the measurements in order to ensure sufficient constraint on the larger wavelengths of the field (Sabaka et al. 2013; Hammer 2018). We denote the magnetic vector measurements by \(B_k(\mathbf {r})\), where k is any of the three given vector component of the field, \(\Delta d_k\) and \(\Sigma d_k\) denote measurement differences and sums of this particular component, respectively. Here the along-track (AT) and East–West (EW) data differences are denoted by \(\Delta d_k=(\Delta d_k^{\mathrm{AT}},\Delta d_k^{\mathrm{EW}})\), and the data sums by \(\Sigma d_k=(\Sigma d_k^{\mathrm{AT}},\Sigma d_k^{\mathrm{EW}})\). The along-track data differences are calculated using the 15 s differences \(\Delta d_k^{\mathrm{AT}} = [B_k(\mathbf {r},t) - B_k(\mathbf {r}+\delta \mathbf {r},t+15s)]\). The along-track sums were calculated as \(\Sigma d_k^{AT} = [B_k(\mathbf {r},t) + B_k(\mathbf {r}+\delta \mathbf {r},t+15s)]/2\). For Swarm, East–West differences were calculated as \(\Delta d_k^{\mathrm{EW}} = [B_k^{\mathrm{SWA}}(\mathbf {r}_1,t_1) - B_k^{\mathrm{SWC}}(\mathbf {r}_2,t_2)]\) having an East–West orbit separation between the Swarm Alpha (SWA) and Charlie (SWC) satellites of \(\approx 1.4^{\circ }\) corresponding to 155 km at the equator (Olsen 2015). The East–West sums were calculated as \(\Sigma d_k^{\mathrm{EW}} = [B_k^{\mathrm{SWA}}(\mathbf {r}_1,t_1) + B_k^{\mathrm{SWC}}(\mathbf {r}_2,t_2)]/2\). For a particular orbit of Swarm Alpha the corresponding Swarm Charlie data were selected to be those closest in colatitude with the condition that \(\vert \Delta t\vert =\vert t_1-t_2\vert <50s\) (Olsen 2015).

Application I: geomagnetic virtual observatories

4-monthly core field GVOs and secular variation estimates

In the first application, we compute Geomagnetic Virtual Observatory time series derived from dataset #1. Of particular interest is the quality of the GVO series obtained from CryoSat-2 data compared with similar series obtained from CHAMP and Swarm data. The time series consist of estimates of the geocentric spherical polar vector components of the magnetic field at specified target points, referred to as GVOs (Mandea and Olsen 2006; Olsen and Mandea 2007). Here we use the same algorithm described in detail by Hammer et al. (2021) (see also http://www.spacecenter.dk/files/magnetic-models/GVO/GVO_Product_Algorithm.pdf) to produce the Swarm GVO product, and derive global grids of 300 uniformly distributed GVO time series each having 4-month cadence. The GVOs are located in an approximately equal area grid based on the sphere computed using the algorithm of Leopardi (2006).

For each GVO in the grid we take data from within a cylinder of horizontal radius \(r_{cyl}=700\,\hbox {km}\). The GVOs have the spherical polar coordinates \(\mathbf {r}_{GVO}=(r,\theta ,\phi )\), and are placed at fixed altitudes \(r=r_a+h_{GVO}\) where \(h_{GVO}\) is the height above the Earth’s mean spherical radius \(r_a=6371.2\,\hbox {km}\). For the CHAMP, CryoSat-2 and Swarm missions the GVO altitudes were chosen as \(h_{GVO}=370\,\hbox {km}\), \(h_{GVO}=727\,\hbox {km}\) and \(h_{GVO}=490\,\hbox {km}\), respectively, such that the GVOs are located at approximately the mean orbital altitude for each mission during the time interval considered.

The input measurements of dataset #1 are provided in an Earth-Centered-Earth-Fixed (ECEF) coordinate frame by the spherical polar components \(\mathbf {B}^{obs}=(B_r,B_{\theta },B_{\phi })\). From the extracted vector field measurements surrounding each GVO, within a radius of 700 km and within a 4-month time window, a residual magnetic field, \(\delta \mathbf {B}\), is first computed by subtracting off estimates of the main field and non-core fields:

$$\begin{aligned} \delta \mathbf {B} = \mathbf {B}^{obs} - \mathbf {B}^{MF} -\mathbf {B}^{lith}- \mathbf {B}^{mag}- \mathbf {B}^{iono}, \end{aligned}$$
(1)

where the field estimates removed are: a) \(\mathbf {B}^{MF}\) the internal field for SH degrees \(n \in [1,13]\) as given by IGRF-13 (Alken et al. 2020b), b) \(\mathbf {B}^{lith}\) the static internal field for SH degrees \(n \in [14,185]\) as given by the LCS-1 model (Olsen et al. 2017), c) \(\mathbf {B}^{mag}\) the magnetospheric and associated induced field as given by the CHAOS-7.2, model (Finlay et al. 2020), d) \(\mathbf {B}^{iono}\) the ionospheric and associated induced field as given by the CIY4 model (Sabaka et al. 2018). Note that estimates of the main field from IGRF-13 (with linear time dependence over 5-year intervals) are subtracted here. At a later stage, main field estimates, again from IGRF-13, are added back for the specified GVO times and positions. Removal of a main field at this stage allows for a more effective pre-whitening of the data, such that Huber weights, used in the robust GVO estimation scheme, can be well determined. We emphasize that this approach still allows us to capture departures from the subtracted main field when required by the satellite data.

Of the remaining residual field, we are interested in the core field part of that signal. Although we have removed estimates of the external fields and their associated Earth-induced counterpart, contributions from non-core sources remain in the residual field. A particular concern is contamination from fields caused by rapidly varying ionospheric currents, for example polar electrojet currents at high latitudes in the E-layer, and possibly also signatures of F-layer currents at mid and low latitudes. Further work is needed on these aspects. We attempt to mitigate leakage of field-aligned currents by removing toroidal field estimates (Sabaka et al. 2010) obtained by performing an epoch-by-epoch spherical harmonic analysis performed on the global network of GVOs (Hammer et al. 2021), see below for further details.

Next, the residual magnetic field data and their positions are transformed from the spherical system to a right-handed local topocentric Cartesian system (x, y, z) with origin at the GVO target location. At the GVO location and only at this location, x points towards geographic south, y points towards east and z points radially upwards (Hammer et al. 2021). Assuming that the satellite measurements are made in a source free region, the residual magnetic field, \(\delta \mathbf {B}\), is a Laplacian potential field. In the local Cartesian coordinate system the magnetic scalar potential, V, can be expanded as a sum of polynomials of the form \(C_{abc} x^a y^b z^c\) (Backus et al. 1996). Here we use this expansion out to cubic terms:

$$\begin{aligned} V(x,y,z)&= C_{100}x + C_{010}y + C_{001}z + C_{200}x^2 + C_{020}y^2 + C_{002}z^2 \nonumber \\ &\quad + C_{110}xy + C_{101}xz + C_{011}yz + C_{300}x^3 + C_{030}y^3 + C_{003}z^3 \nonumber \\&\quad + C_{210}x^2y+ C_{201}x^2z + C_{120}y^2x + C_{021}y^2z + C_{102}z^2x + C_{012}z^2y \nonumber \\ &\quad + C_{111}xyz. \end{aligned}$$
(2)

The forward problem linking the vector of GVO model coefficients, \(\mathbf {m}={[C_{100},C_{010},..., C_{111}]}^T\), with the data vector \(\mathbf {d}^{vec}\) containing the residual field components, \(\delta \mathbf {B}\), can be written

$$\begin{aligned} \mathbf {d}^{vec} = \underline{\underline{\mathbf {G}}}^{vec} \mathbf {m}, \end{aligned}$$
(3)

where \(\underline{\underline{\mathbf {G}}}^{vec}\) is a design matrix derived from the spatial derivatives of Eq. (2). As noted above, instead of using residual vector field components to compute the potential, we use sums and differences of the residual field vector components such that the data vector is \(\mathbf {d}=[\Delta d_{x}^{vec},\Delta d_{y}^{vec},\Delta d_{z}^{vec},\Sigma d_{x}^{vec},\Sigma d_{y}^{vec},\Sigma d_{z}^{vec}]^T\), where \(\Delta\) and \(\Sigma\) denotes the differences and sums of the computed residual field as described in Sect. 2. The relevant design matrix is then constructed as \(\underline{\underline{\mathbf {G}}}=[\Delta G_x^{vec};\Delta G_y^{vec};\Delta G_z^{vec};\Sigma G_x^{vec};\Sigma G_y^{vec};\Sigma G_z^{vec})]\) where \(\Delta G_{k}^{vec} =[G_{k}^{vec}(\mathbf {r}_1) - G_{k}^{vec}(\mathbf {r}_2)]\) and \(\Sigma G_{k}^{vec} = [G_{k}^{vec}(\mathbf {r}_1) +G_{k}^{vec}(\mathbf {r}_2)]/2\), where \(k=(x,y,z)\).

To determine the GVO model coefficients, we use a robust iteratively-reweighted least-squares inversion scheme, based on a diagonal weight matrix consisting of robust (Huber) weights for each entry in the data vector (e.g., Constable 1988). We also include a down-weighting factor of 1/2 for the Swarm satellites Alpha and Charlie accounting for the fact that these two satellites fly side-by-side and therefore do not provide completely independent measurements. A minimum number of 30 data points were required to computing the inversion. Using the resulting coefficients of the potential for a given GVO target location and time, derived from the associated sums and differences satellite data, a prediction for the mean residual field at the GVO target point and epoch can be computed as \(\delta \mathbf {B}_{GVO}(x,y,z)=-\nabla V(0,0,0)=-(C_{100},C_{010},C_{001})\).

Moving back to the vector components in spherical polar coordinates, \(\delta B_{GVO,r}=\delta B_{GVO,z}\), \(\delta B_{GVO,\theta }=\delta B_{GVO,x}\), \(\delta B_{GVO,\phi }=\delta B_{GVO,y}\). We then add back the IGRF-13 main field predictions for the given target point and epoch, \(\mathbf {B}_{GVO}^{MF}(\mathbf {r}_{GVO},t)\) to obtain

$$\begin{aligned} \mathbf {B}_{GVO}(\mathbf {r}_{GVO},t) = \delta \mathbf {B}_{GVO}(\mathbf {r}_{GVO},t) + \mathbf {B}_{GVO}^{MF}(\mathbf {r}_{GVO},t). \end{aligned}$$
(4)

The above procedure is then repeated for each epoch and each component to obtain time series of GVO estimates of the vector magnetic field at the GVO target locations.

The GVO method assumes that the residual field in Eq.(1), is a potential field. However, because the satellite measurements are made in the ionospheric F-region, in situ currents can cause non-potential fields to leak into the estimated potential (Olsen and Mandea 2007). Therefore, in a final post-processing step, we carry out an epoch-by epoch spherical harmonic analysis of our 4-monthly GVOs, estimating external and toroidal field contributions to SH degree 13 or to degree 6 if fewer than 300 GVOs are available. For epochs having an insufficient number of GVOs available to ensure a stable solution, the external and toroidal coefficients were obtained via a linear interpolation between nearby epochs. The external and toroidal field estimates, reaching a level of \(\pm 15\)nT at high latitudes, are then removed to obtain the estimated Core Field GVO time series (Hammer et al. 2021).

Secular variation at a given GVO location is computed using annual differences between values at time \(t+6\) months and at time \(t-6\) months. We have chosen to take annual differences, as this helps to avoid annual non-core signals that may persists in the GVO series despite our best efforts in reducing such contamination.

Figure 1 presents the number of 4-monthly Core Field GVO estimates during the past 20 years. The maximum possible number of GVOs per epoch is 300. A strong dip in the number of GVOs is seen during 2002-2004 due to increased solar activity that meant there were fewer data meeting our selection criteria. As noted above, if fewer than 30 measurements are available within a GVO target cylinder during a given 4-month window, we were unable to reliably determine GVO estimates. The remaining epochs from 2004-2020 are well covered with only few epochs having less than 250 GVOs available. We find that using CryoSat-2 data, we are able to provide between 200 and 300 GVO estimates at all times during the gap between the end of the CHAMP mission and the start of the Swarm mission.

Fig. 1
figure 1

The available number of GVOs for each epoch during CHAMP (purple), CryoSat-2 (pink) and Swarm (green)

Full size image

GVO results: global time series of secular variation from 2002 to 2020

Table 2 presents the root-mean-square (rms) and mean of the residuals between the contributing satellite data (sums and differences) and GVO model predictions, summing over all GVOs for a given vector component and a given region (polar or non-polar). Here we defined polar to be polward of \(\pm 54^{\circ }\) geographic latitude. The polar rms values for both sums and differences are higher than the non-polar, and the CHAMP values are slightly higher than the Swarm values. The CryoSat-2 values are seen to be higher for all components but not unreasonable, given they are derived from platform magnetometer data. The non-polar rms values for all components are below 2 nT during both CHAMP and Swarm times. The CryoSat-2 GVO’s rms values are as expected larger, and especially for the along-track differences, indicating that along-track correlated noise is less dominant, due to the presence of other noise sources in platform magnetometer data.

Table 2 GVO model rms misfit statistics between GVO estimates and the contributing data for the global grid of 300 GVOs during CHAMP, CryoSat-2 and Swarm. Here \(\scriptsize \sum\) and \(\Delta\) represent data sums and data differences, respectively, split into the North-South (NS) and East–West (EW) components
Full size table

The CryoSat-2 magnetometer data have been cleaned from known platform signals and calibrated as described by Olsen et al. (2020). This calibration relies on computing residuals with respect to a reference magnetic field which was taken from the CHAOS-6-x9 field model (Finlay et al. 2016). Here, we carried out an experiment to verify that the GVO secular variation signals obtain from CryoSat-2 data are independent of the main field model used for data calibration. To do this we computed CryoSat-2 GVO estimates using two different datasets, the first being the official CryoSat-2 dataset calibrated using the CHAOS-6-x9 field model and the second a test version of the CryoSat-2 data calibrated instead using IGRF-13 (Alken et al. 2020b). Having estimated global grids of GVO series for each dataset, to each series we fit cubic smoothing splines, with a knot spacing at every 4 months, and with the smoothing parameter determined using a GCV (generalized cross-validation) approach (Green and Silverman 1993). Figure 2 presents SV series for the three field components at three example GVOs, with colatitude/longitude \((36^{\circ },77^{\circ })\) (top), \((72^{\circ },110^{\circ })\) (Center) and \((84^{\circ },-137^{\circ })\) (bottom). The SV time series derived from the CHAOS-6x9 calibrated dataset are shown with red dots (spline fit in the red line) and similar series from the IGRF-13 calibrated dataset are shown with blue dots (spline fit in the blue line). SV model predictions from the CHAOS-6-x9 model up to SH degree 16 (in green), and IGRF-13 (in black), are also shown for reference. There are certainly some differences in the GVO SV estimates derived using the CHAOS-6x9 and IGRF-13 calibrated data. Since the CryoSat-2 calibration relies on computing residuals with respect to a reference field model (here either CHAOS or IGRF), we do expect such differences, especially since the IGRF model assumes a crude piecewise constant SV field. Nevertheless, the IGRF calibrated GVO SV estimates are clearly closer to the CHAOS-6x9 SV predictions than to the IGRF-13 SV predictions, by an rms of 1.2nT/yr for the horizontal components, and 3.9nT/yr for the radial component. This indicates that the CryoSat-2 data does possess an SV signal regardless of the model used to calibrate them.

We find that the CHAOS-6x9 and IGRF-13 calibrated CryoSat-2 GVO SV series show similar sub-decadal changes, neither of which exactly match those seen in CHAOS-6x9. The overall similarity of the IGRF and CHAOS calibrated SV series give us confidence that the sub-decadal secular variation seen in CryoSat-2 data is independent of the field model used in its calibration. In addition, we note a strong acceleration in the radial SV component from 2014 to 2018 in both sets of GVOs (top plot) located above Northwest Siberia, a change which is only partially captured by the piecewise constant IGRF-13 SV model.

Fig. 2
figure 2

Time series of annual differences of the spherical polar field components of CryoSat-2 GVO time series derived using CHAOS-6-x9 (red) and IGRF-13 (blue) calibration along with GCV spline fits at the three example grid locations: 26 (top), 100 (center) and 140 (bottom). Also shown are IGRF-13 (black) and CHAOS6-x9 predictions (green)

Full size image

In Table 3, we report the mean rms differences between the GVO SV estimates and GCV spline fits, separated by component (\(B_r,B_\theta ,B_\phi\) and also the intensity F) and by polar and non-polar regions. These numbers give an indication of the scatter in the SV datasets and allow the quality of the GVO SV series obtained from CHAMP, Swarm and CryoSat-2 (both CHAOS-6x9 and IGRF-13 calibrated versions) to be compared. Similar results are seen for the CHAOS-6-x9 and IGRF-13 model calibrated GVO SV series, with rms differences between GCV fits and intensity SV data, taken the over all GVOs, of \(3.5\,\hbox {nT}/\hbox {yr}\) and \(3.4\,\hbox {nT}/\hbox {yr}\), respectively. As expected, similar numbers for the CHAMP and especially the Swarm derived GVOs, are much smaller being \(1.8\,\hbox {nT}/\hbox {yr}\) and \(0.9\,\hbox {nT}/\hbox {yr}\), respectively. This indicates the lower scatter in the CHAMP and Swarm GVO series.

Table 3 Mean of the rms differences (in nT/year) between GVO SV series and GCV cubic spline fits. Results are shown for both the CHAOS-6-x9 and IGRF-13 calibrated CryoSat-2 datasets, and for GVO SV series derived using Swarm and CHAMP data
Full size table
Fig. 3
figure 3

Time series of annual differences of the radial field component for GVOs derived from CryoSat-2 data calibrated using CHAOS-6-x9 (red dots) and IGRF-13 (blue dots), and relocated to a common altitude of 700 km. GVO locations are marked with a black cross

Full size image

Figure 3 goes beyond detailed comparisons at a few example locations and presents a global map of CryoSat-2 GVO time series for the radial field SV, showing both the official CHAOS-6x9 calibrated dataset (red dots) and the IGRF-13 calibrated test dataset (blue dots). The scale is shown in the bottom left corner, with the y-axis being 20 nT/yr and the x-axis going from 2011 to 2018. Visual inspection of Fig. 3 supports the results of Table 3, and shows similar scatter levels with rms differences of less than \(1.0\hbox {nT}/\hbox {yr}\) in all three components of the SV series derived from the CHAOS-6-x9 and IGRF-13 calibrations.

With the availability of CryoSat-2 magnetic data, it is now possible to use GVOs to study secular variation globally over the past 20 years. In order to illustrate the quality of information they can provide, in Fig. 4 we first present comparisons of GVO SV series, from CHAMP, CryoSat-2 and Swarm, to annual differences of Revised Monthly Means (rmm) (Olsen et al. 2014) from example high quality ground observatories, Kourou in South America (top plots), from Novosibirsk observatory in Siberia (middle plots) and from Honolulu ground observatory in the central Pacific (bottom plots). Each plot shows the spherical polar components of the annual differences of revised monthly means (black dots) computed from the ground observatory data (Olsen et al. 2014), and GVO SV time series derived from CHAMP (purple dots), CHAOS-6x9 calibrated CryoSat-2 (blue dots) and Swarm (red dots) data mapped to ground level. Here the GVO time series for each satellite mission have been mapped to ground by subtracting the field difference as given by the CHAOS-7.2 model for SH degrees 1-20, between ground and GVO altitudes which are close to the mean orbital altitude for each mission. At all three locations the GVOs derived from CryoSat-2 data have more scatter – this is particularly noticeable in the \(\theta\)- and \(\phi\)-components for the examples shown, whereas the scatter in the radial component is closer to the level seen in the CHAMP and Swarm derived GVOs. Table 3 indicates that at non-polar latitudes the scatter in the r- and \(\theta\)-components is generally similar, 3.4 nT/year compared to 3.8 nT/year. Notice that the scatter in the ground observatory rmm’s is also enhanced in the horizontal components. The good agreement between the independently estimated CryoSat-2 and Swarm-derived GVOs at overlapping epochs from 2014 to 2018, is particularly evident in the r- and \(\phi\)-components. SV variations are coherent in both phase and amplitude between the CryoSat-2 GVO time series and the ground observatory records, thus confirming that the CryoSat-2 GVOs are able to track the same field changes as observed by ground observatory records on timescales of 1 year and longer.

Fig. 4
figure 4

SV times series obtained from annual differences of revised monthly means (black dots) (Olsen et al. 2014), and 4-monthly Core Field GVOs derived from CHAMP (purple dots), CryoSat-2 (blue dots) and Swarm (red dots) data. Note the GVO estimates have been mapped from satellite altitude to ground in order to aid the comparison

Full size image

Next, in Fig. 5 we present a global map of the SV of the radial field component over the past 18 years. Here, to ease visualization, GVO SV series from CHAMP (covering 2002-2010), CHAOS-6x9 calibrated CryoSat-2 (covering 2010-2014) and Swarm (covering 2014-2020) have been mapped to a common altitude of 700 km, again using the CHAOS-7.2 field model, and combined into one composite time series. This allows for easy visual investigation of global patterns of sub-decadal SV. We find that regions at low latitudes display strong sub-decadal variations, however, not simultaneously at all longitudes. For instance, we observe a change of slope in the radial SV field occurring over the south Atlantic region around 2007 and again in 2014, over Indonesia around 2014, and in the Pacific region centered in 2017. Some of these variations are characterized by distinct “\(\Lambda\)” and “V”-shaped behavior occurring over time spans of 5–10 years and locally confined to specific regions, but otherwise reminiscence of the often discussed phenomenon of geomagnetic jerks (Mandea et al. 2010). The availability of CryoSat-2 magnetic field data plays a key role in permitting a continuous coverage from satellite-based time series without a gap between 2010 and 2013, thus allowing the study of global patterns in the time-varying core field secular variation over the past 20 years.

Fig. 5
figure 5

Composite GVO SV time series for the radial field component, computed using annual differences of 4-monthly core field GVOs (red dots) which have been mapped to a common altitude of 700 km. GVO locations are marked with a black cross

Full size image

Application II: SOLA

SOLA method for local estimation of CMB radial field SV

We now move on to investigate the behavior of the core field not at satellite altitude but down at the core–mantle boundary (CMB), on the edge of the region where it originates. To do this we use the Green’s functions of the Neumann boundary value problem that links the magnetic field at satellite altitude to the radial field at the CMB. Following Hammer and Finlay (2019), a localized estimate, \(\widehat{B}_r\), of the radial magnetic field at a target location and time, \((\mathbf {r}_0,t_0)\), at the CMB can be computed as a localized spatial average around the target location time-averaged over a specified interval. Because the CMB radial magnetic field is linearly related to the spherical polar components of the vector field at satellite altitude, we can write \(\widehat{B}_r\) as a weighted linear combination of the satellite magnetic measurements with weights \(q_n\) (Backus and Gilbert 1968, 1970; Hammer and Finlay 2019)

$$\begin{aligned} \widehat{B}_r(\mathbf {r}_0,t_0) = \sum _n^N q_n(\mathbf {r}_0,t_0)\, d_n(\mathbf {r}_n,t_n), \end{aligned}$$
(5)

where \(d_n\) are satellite magnetic measurements \((n=1,...,N)\) within a specified time window and \(q_n\) are weighting coefficients to be determined. Here the data \(d_{n}\), at positions \(\mathbf {r}_n\) and times \(t_n\), are taken from dataset #2 and for simplicity we consider using only observations of the radial component of the field. Corrections for the lithospheric field for SH degrees \(n \in [14,185]\) as given by the LCS-1 model (Olsen et al. 2017), for the magnetospheric and associated induced fields as given by the CHAOS-7.2 model (Finlay et al. 2020), and for the ionospheric and associated induced fields as given by the CIY4 model, (Sabaka et al. 2018) are removed from the observations in a pre-processing step. The radial magnetic field measurements, \(d_{n}(\mathbf {r}_n,t_n)\), are then related to the radial magnetic field \(B_r(\mathbf {r}',t_n)\), integrated over the CMB, by (Gubbins and Roberts 1983):

$$\begin{aligned} d_{n}(\mathbf {r}_n,t_n) = \oint _{S'} G_{r}(\mathbf {r}_n, \mathbf {r}')B_r(\mathbf {r}',t_n) dS', \end{aligned}$$
(6)

where the surface element is \(dS'=\mathrm{sin}\theta 'd\theta 'd\phi '\). The data kernel \(G_{r}(\mathbf {r}_n, \mathbf {r}')\) is the radial derivative with respect to \(\mathbf {r}\), of the Green’s functions for the exterior Neumann boundary value problem (e.g., Gubbins and Roberts 1983; Barton 1989):

$$\begin{aligned} G_{r} = \frac{1}{4\pi }\frac{h_n^2(1-h_n^2)}{f_n^3}, \end{aligned}$$
(7)

where \(h_n=r'/r_n\) and \(r'\) is the CMB radius, \(f_n=R_n/r_n\) where \(R_n=\sqrt{r_n^2+r'^2-2r_nr'\zeta _n}\) and \(\zeta _n = \mathrm{cos}\, \gamma _n=\mathrm{cos}\,\theta _n \, \mathrm{cos}\,\theta '+\mathrm{sin}\,\theta _n \mathrm{sin}\,\theta '\, \mathrm{cos(\phi _n-\phi ')}\), where \(\gamma _n\) is the angular distance between a measurement at position \((\theta _n,\phi _n)\) and a position on the CMB \((\theta ',\phi ')\). The data kernel describes how a particular measurement samples the CMB radial field; radial magnetic measurements sample the CMB radial field most strongly directly below the measurement position. Regarding the time-dependence, we use a first order Taylor expansion around a reference time \(t_0\), such that

$$\begin{aligned} d_{n}(\mathbf {r}_n,t_n) \approx \oint _{S'} G_{r}(\mathbf {r}_n,\mathbf {r}') \left[ B_r(\mathbf {r}',t_0)+\dot{B_r}(\mathbf {r}',t_0)\Delta t_n\right] dS'. \end{aligned}$$
(8)

The time difference, \(\Delta t_n = t_n-t_0\), is computed with respect to the target time, \(t_0\). Inserting Eq.(8) into Eq.(5), we obtain

$$\begin{aligned} \widehat{B}_r(\mathbf {r}_0,t_0) &= \oint _{S'} {\mathcal {K}}(\mathbf {r}_0, t_0, \mathbf {r}') \, B_r(\mathbf {r}',t_0)dS' \nonumber \\ & + \oint _{S'} \dot{{\mathcal {K}}} (\mathbf {r}_0, t_0, \mathbf {r}') \, \dot{B_r}(\mathbf {r}',t_0) dS', \end{aligned}$$
(9)

where \(\mathcal {K}(\mathbf {r}_0, {t_0}, \mathbf {r}')\) and \(\dot{{\mathcal {K}}}(\mathbf {r}_0, t_0, \mathbf {r}')\) are spatial averaging kernels for the CMB field and secular variation, respectively, constructed from the weighting coefficients and the data kernels

$$\begin{aligned} {\mathcal {K}}(\mathbf {r}_0, t_0, \mathbf {r}')&= \sum _n^N q_n(\mathbf {r}_0,t_0) \, G_r(\mathbf {r}_n,\mathbf {r}') \\ \dot{{\mathcal {K}}}(\mathbf {r}_0, t_0, \mathbf {r}') &= \sum _n^N q_n(\mathbf {r}_0,t_0) \, G_r(\mathbf {r}_n,\mathbf {r}')\Delta t_n. \end{aligned}$$
(10)

By varying the weight coefficients, \(q_n\), the shape of the averaging kernels change. Notice that time differences \(\Delta t_n\), between the measurement times and the target time, are effectively additional temporal weights applied to the kernel \({\mathcal {K}}\) in order to obtain \(\dot{{\mathcal {K}}}\).

In order to obtain estimates of the secular variation of the radial field on the CMB, at the target location and time \(\widehat{\dot{B}}_r(\mathbf {r}_0,t_0)\), we minimize the following objective function:

$$\begin{aligned} \Theta &= \oint _{S'} [\dot{{\mathcal {K}}}(\mathbf {r}_0, t_0, \mathbf {r}') -\dot{{\mathcal {T}}}(\mathbf {r}_0, \mathbf {r}')]^2dS' \nonumber \\ & \oint _{S'} [{\mathcal {K}}(\mathbf {r}_0, t_0, \mathbf {r}')]^2dS' + \lambda ^2 \mathbf {q}^T \underline{\underline{\mathbf {E}}} \mathbf {q}, \end{aligned}$$
(11)

where \(\lambda\) is a trade-off parameter (units of \([\mathrm{nT}^{-1}]\)), \(\mathbf {q}\) is vector of the weighting coefficients, \(\underline{\underline{\mathbf {E}}}\) is the data error covariance matrix which we define below and \(\dot{{\mathcal {T}}}\) is an SV target kernel that we choose to be a Fisher distribution on the sphere (Fisher 1953):

$$\begin{aligned} \dot{{\mathcal {T}}}(\mathbf {r}_0, \mathbf {r}') =\frac{\kappa }{4 \pi \mathrm{sinh}\kappa } e^{\kappa \cos \, \gamma _0}, \end{aligned}$$
(12)

where \(\gamma _0\) is the angular distance on the CMB between the target position \((\theta _0,\phi _0)\) and another position \((\theta ',\phi ')\). On the basis of tests carried out by Hammer and Finlay (2019), we initially set \(\kappa =600\) corresponding to a target kernel width of \(15^{\circ }\); this is narrower than can be achieved for \(\dot{{\mathcal {K}}}\) with the available data, but it avoids excessive ringing associated with taking a Dirac delta function as the target kernel. When computing SOLA estimates for a given time window, we select a subset (\(n=1,...,N\)) of the measurements. Using this data subset the data error covariance matrix \(\mathbf {E}\) is defined as follows. Using all available measurements (\(m=1,...,M\)), for each satellite mission within 2 degree bins of quasi-dipole (QD) latitude (Richmond 1995), we first derived robust data error variances as a function of QD latitude:

$$\begin{aligned} \sigma ^2(\theta _{QD}) = { \sum \limits _{m=1}^M w_m\left( \epsilon _m-\mu \right) ^2 \over \sum \limits _{m=1}^M w_m }, \end{aligned}$$
(13)

where \(\epsilon _m\) are residuals with respect to predictions of the CHAOS-7.2 internal field model for SH degrees \(n \in [1,13]\), \(\mu\) are robust mean residuals within the considered bin and \(w_m\) are Huber weights (e.g., Constable 1988) for the data within each bin. Here we use QD coordinates, as this is appropriate for characterizing processes related to unmodeled ionospheric currents which we consider to be a likely source of contamination, especially at high latitudes. Figure 6 presents the resulting QD-latitude-dependent error estimates \(\sigma (\theta _{QD})\) for the radial field component used in this study, comparing the values for the Ørsted, CHAMP, CryoSat-2 and Swarm datasets. When computing SOLA estimates for a specified time window of, e.g., 2 years, we select a data subset of dataset #2. Using this data subset of N measurements, a specific data error covariance matrix E is computed. Diagonal elements of this data error covariance matrix \(\mathbf {E}\) are finally defined as:

$$\begin{aligned} E_{n,n}=\sigma ^2(\theta _{QD}) / w_n, \end{aligned}$$
(14)

where \(w_n\) are additional robust (Huber) weights determined a priori for each datum (\(n=1,...,N\)), based on their residual to CHAOS-7.2, in order to account for the expected long-tailed error distribution. Off-diagonal elements of \(\mathbf {E}\) are set to zero.

In addition to minimizing Eq.(11), we simultaneously impose the following constraint:

$$\begin{aligned} \oint _{S'}{\mathcal {K}}(\mathbf {r}_0, t_0, \mathbf {r}') dS' + \oint _{S'} \dot{{\mathcal {K}}}(\mathbf {r}_0, t_0, \mathbf {r}')dS' = 1, \end{aligned}$$
(15)

where the first term is in practice very small when estimating the SV, since it is minimized in Eq. (11). This constraint ensures that a valid averaging kernel is obtained.

Discretization of integrals over the CMB was carried out using Lebedev quadrature (Lebedev and Laikov 1999) and the system of equations was solved for the coefficients, \(q_n\), using a Lagrange multiplier method, see Hammer and Finlay (2019) for further details.

Once SOLA estimates of the CMB radial field SV at the chosen target location and epoch are obtained, by minimizing Eq. (11) subject Eq. (15), we are able to easily appraise them based on (i) their averaging kernel width, which we define as the angular distance between the points at which the averaging kernel first reaches zero amplitude moving away from its maximum value, and (ii) the variance of the SOLA estimate which we computed as:

$$\begin{aligned} \hat{\sigma }^2(\mathbf {r}_0,t_0) = \mathbf {q}^T \underline{\underline{\mathbf {E}}}\mathbf {q}. \end{aligned}$$
(16)

By changing the parameter \(\lambda\) of Eq.(11), a range of solutions can be computed which describes a trade-off between having an averaging kernel width as small as possible and the variance of the estimate being as small as possible (Parker 1977). Below we discuss the effect of changing \(\lambda\) on our results.

Fig. 6
figure 6

Latitude-dependent standard deviations, \(\sigma (\theta _{QD})\), contributing to the data error budget. For the radial field component in \(2^{\circ }\) bins (Northern hemisphere having positive QD latitude) for each satellite mission

Full size image

Results: SOLA estimates of CMB SV and SA from Ørsted, CHAMP, CryoSat-2 and Swarm data

We begin by first comparing SOLA estimates for the CMB radial field SV obtained using separate data subsets from the Swarm and CryoSat-2 missions, respectively. First, Swarm and CryoSat-2 data subsets are extracted from the main dataset #2 described in Sect. 2, so that each cover the same 2-year time window from 2015.0 to 2017.0 Next, in order to obtain data subsets with suitable spatial and temporal coverage, we considered bins surrounding each point in an approximately equal-distance grid at satellite altitude of \(\approx 2.5^{\circ }\) spacing, based on the partitioning algorithm of Leopardi (2006), and randomly sampled one datapoint from each bin, resetting the bins every 2 months. Data subsets spanning the full 2-year window from 2015.0 to 2017.0 were produced by accumulating these 2-monthly globally-distributed subsets. The resulting Swarm and CryoSat-2 data subsets spanning 2015.0 to 2017.0 consisted of 62469 and 54685 radial field observations, respectively.

In Fig. 7, we compare maps collecting SOLA CMB radial field SV estimates centered on epoch 2016.0, derived using the CryoSat-2 and Swarm data subsets spanning 2015.0 to 2017.0. To ensure that we obtained SOLA estimates of comparable resolution, we first computed SOLA SV estimates using the Swarm data subset and taking \(\lambda =3 \times 10^{-3}\ \mathrm{nT}^{-1}\). This resulted in well-behaved averaging kernels with widths \(\approx 38^{\circ }\). Next, we used these averaging kernels as the target kernels in order to derive similar estimates using the Swarm and CryoSat-2 data subsets, thus effectively seeking Swarm and CryoSat-2 SV estimates with the same spatial resolution. The maps in the top row of Fig. 7 show the resulting global collections of SOLA SV estimates, obtained on a \(1^{\circ }\) grid of target locations at the CMB, based on the CryoSat-2 (left plot) and Swarm (right plot) data subsets, respectively.

Fig. 7
figure 7

SOLA CMB radial field SV estimates for epoch 2016.0, derived using data from 2015.0 to 2017.0 (top plots), associated error estimates (middle plots) and averaging kernel widths (bottom plots). Estimates are derived using CryoSat-2 (left plots) and Swarm (right plots) data, respectively

Full size image

Maps collecting the formal standard deviations for each SOLA estimate (derived using Eq.16) and their averaging kernel widths are presented in middle and bottom rows of Fig. 7. The standard deviations of the Swarm-based SOLA SV estimates are fairly homogeneous with values of \(0.3-0.4\mu \hbox {T}/\hbox {yr}\); those for CryoSat-2 SOLA are somewhat larger, being in the range \(\approx 0.6-1.8\mu \hbox {T}/\hbox {yr}\). We note that the CryoSat-2 error estimates are slightly lower at higher latitudes. The same is true for the Swarm map, but the variations in errors estimates are in that case much smaller. Kernel widths in both cases are also fairly homogeneous except at auroral latitudes where distinct behavior of the kernels are found related to the data error estimates having increased amplitude, as seen in Fig. 6.

As seen from the kernel widths in the bottom row of Fig. 7, very similar resolution has been obtained (by construction) in the Swarm and CryoSat-2 SV maps. The same field features are clearly identified in both maps. For instance, we notice high latitude SV patches in the northern hemisphere which have been associated with a high latitude jet of core flow (Livermore et al. 2017), and there is increased amplitude of SV over the hemisphere centered on the Atlantic in comparison with the Pacific hemisphere. This first test gives us confidence that the CryoSat-2 measurements can indeed be used to reliably map SV features at the CMB on a timescale of 2 years, and with a spatial resolution down to \(\approx 38^{\circ }\) degrees.

Next, we go further and investigate the second time derivative or secular acceleration (SA) of the radial field at the CMB, which is of great interest for investigating the dynamics of the core (e.g., Chulliat et al. 2010; Finlay et al. 2016; Chi-Durán et al. 2020). Again we first compare maps based on CryoSat-2 and Swarm data. To obtain SA estimates we initially use the accumulated change between SV estimates 2 years apart. In particular, in Fig. 8 we show the SA in 2017 based on the difference between SOLA SV estimates in 2016.0 and 2018.0. Note, that this is not an instantaneous secular acceleration but a centered difference in SV estimates 2 years apart, each based on 2 years of data. To study the SA, we computed SOLA SV estimates from Swarm data taking \(\lambda =1 \times 10^{-2}\ \mathrm{nT}^{-1}\), and used the resulting associated averaging kernels as the target kernels for the SOLA estimates from both Swarm and CryoSat-2 data.

Figure 8 presents global grids of SOLA CMB radial field SA estimates centered on 2017.0, again with a \(1^{\circ }\) spacing, derived from CryoSat-2 (left plot) and Swarm (right plot) data subsets. Here the map is centered on the Pacific region where there has been interesting SA activity during the past 6 years (Finlay et al. 2016, 2020). As for the SV maps, we find error estimates, computed assuming the contributing SV estimates have independent errors, to be fairly homogeneous, with values ranging between \(0.11-0.27\mu \mathrm{T}/\mathrm{yr}^2\) for the estimates derived using CryoSat-2 data and \(0.06-0.08\mu \mathrm{T}/\mathrm{yr}^2\) for the estimates derived using Swarm data. In both cases kernel widths are close to \(\approx 42^{\circ }\), except in the auroral region. Comparing the CryoSat-2 and Swarm-based SA maps, similar features can be observed. In particular, this is the case for the features seen under Asia and Indonesia. A distinctive feature reproduced in both maps is the sequence of intense patches of SA at low latitudes in a localized region below central America, having amplitudes of approximately \(1.9\pm 0.3\mu \mathrm{T}/\mathrm{yr}^2\) and \(1.7\pm 0.1\mu \mathrm{T}/\mathrm{yr}^2\) for the CryoSat-2 and Swarm maps, respectively. The location and amplitude of these features are similar in the two maps, confirming that CryoSat-2 data can be used to track such SA structures. Because the SOLA estimates are local averages, distant high latitude measurements, where ionospheric electrical currents may be prominent even during dark quiet times, will have little influence on such low latitude SV and SA estimates. Finally, we note a strong patch under the Bering Sea of amplitude approximately \(1.4\pm 0.3\mu \mathrm{T}/\mathrm{yr}^2\) and \(1.0\pm 0.1\mu \mathrm{T}/\mathrm{yr}^2\) for the CryoSat-2 and Swarm maps, respectively.

Fig. 8
figure 8

Maps of radial SA at the CMB centered on 2017.0 derived from differences of SOLA SV estimates in 2016.0 and 2018.0 using data from 2015–2017 and 2017–2019

Full size image

These initial investigations of the CMB radial field SV and SA using the SOLA technique indicate, as seen earlier in the GVO time series, that low latitude regions experience significant sub-decadal core field variations. We are therefore motivated to study the field time-dependence in the equatorial region in more detail. We do this in Fig. 9 by constructing time-longitude (TL) plots of SOLA CMB radial field SA estimates along the geographic equator from 2002 to 2019, centered on the Pacific. We again compute our SA estimates based on differences of SOLA SV estimates 2 years apart, each derived from 2-year data windows, and sliding the windows in 2-month steps. Here we use radial field data from the Ørsted, CHAMP, CryoSat-2 and Swarm satellites.

The top left plot in Fig. 9 presents estimates based on averaging kernels obtained from Swarm data using \(\lambda =3 \times 10^{-2}\ \mathrm{nT}^{-1}\). These have associated error estimates between \(0.02-0.08\mu \mathrm{T}/\mathrm{yr}^2\) and kernel widths \(\approx 50^{\circ }\). Their resolution is lower than that shown in Fig. 8 and corresponds to approximately SH degree 8. For comparison the SA predicted by the CHAOS-7.2 model for SH degrees \(n \in [1,8]\) are shown on the top right plot. Note here that the CHAOS-7.2 model makes use of uncalibrated vector magnetic data from CryoSat-2 between 2010 and 2014, and co-estimates magnetometer calibration parameters (Finlay et al. 2020). The SOLA and CHAOS TL plots in the top row of Fig. 9 show largely the same SA features, illustrating the convergence of the two techniques at long wavelengths of the SA and when constructing SOLA SA estimates from SV differences between consecutive 2-year time windows. For instance, the evolution of the SA features observed in Fig. 8 under the central Americas, can be identified ranging from longitudes \(240^{\circ }\) to \(320^{\circ }\) centered on 2017. In addition, we find strong SA patches in the CryoSat-2 data around 2013 at longitudes \(70^{\circ }\) to \(160^{\circ }\) and \(280^{\circ }\) to \(320^{\circ }\). Notice, that there seems to be a sign changing sequence occurring at longitudes \(240^{\circ }\) to \(320^{\circ }\) going from 2005 to 2019.

Next, we increase the spatial resolution by instead deriving SV estimates using \(\lambda =1 \times 10^{-2}\ \mathrm{nT}^{-1}\), which leads to slightly larger error estimates in the range \(0.1-0.5\mu \mathrm{T}/\mathrm{yr}^2\) and kernel widths \(\approx 42^{\circ }\), i.e., similar to Fig. 8. This is shown in the bottom left plot while the CHAOS-7.2 model predictions for SH degrees \(n \in [1,10]\), matching approximately the kernel width, are shown on the bottom right plot. Although the SOLA TL-plot looks somewhat noisier than the CHAOS plot, similar coherent evolving structures having higher amplitudes can clearly be identified. The noisier appearance in the interval 2010-2014 likely indicates the limitations of the CryoSat-2 data, but they clearly provide useful information during this period.

Fig. 9
figure 9

Time-longitude plots of the SA along the geographical equator at the CMB. Left column: SOLA radial field SA estimates derived from differences between consecutive SOLA SV estimates based on 2-year data windows moving in 2-month steps and using \(\lambda =3 \times 10^{-2}\ \mathrm{nT}^{-1}\) (top left), and \(\lambda =1 \times 10^{-2}\ \mathrm{nT}^{-1}\) (bottom left). Right column: SA predictions from the CHAOS-7.2 model for SH degrees up to 8 (top right) and 10 (bottom right); note that temporal regularization causes significant time-averaging of the high-degree SA in the CHAOS model

Full size image

With data from the Swarm mission, it is possible to go further and also increase the temporal resolution of the SA by taking 1-year differences of SOLA SV estimates derived from 1-year data windows and sliding in 1-month steps. The result of applying this procedure to obtain SA estimates on the geographic equator between 2015.0 and 2019.5 is shown in the left plot of Fig. 10. These SOLA estimates have associated errors of \(0.3-0.6\mu \mathrm{T}/\mathrm{yr}^2\) and kernel widths \(\approx 42^{\circ }\). The right plot shows similar CHAOS-7.2 model predictions for SH degrees \(n \in [1,10]\). Both TL-plots shows the similar large-scale features, for instance, the features under central America from longitudes \(240^{\circ }\) to \(320^{\circ }\), which are elongated compared with Fig. 9 due to the change in scale of the y-axis (time). However the SOLA results show significantly more time-dependence, revealing features that were smoothed out by the temporal regularization of CHAOS-7.2. Changes of sign in the SA within about 1 year can be observed.

Fig. 10
figure 10

Time–longitude plots along the geographical equator at the CMB, showing: 1-year differences of SOLA SV estimates derived from 1-year data windows moving in 1-month steps (left), and SA predictions from the CHAOS-7.2 model for SH degrees \(n \in [1,10]\) (right)

Full size image

Particularly interesting is the appearance in the Pacific region around 2017, at longitudes \(150^{\circ }\) to \(220^{\circ }\), of side-by-side positive and negative intense SA features, that have subsequently drifted westwards. This SA change coincides with the peak in the radial SV field observed in the Pacific region during Swarm time seen in the GVO map (Fig. 5). We note the presence of features in Fig. 10 that appear to drift rapidly both eastwards and westwards, for example from \(160^{\circ }\) East in 2015 to \(220^{\circ }\) East in 2017. Such rapidly drifting behavior of low latitude SA patches is difficult to explain in terms of simple core flow advection processes. They may instead be a signature of wave propagation close to the core surface. A range of possible candidates for fast waves in the core have recently been described, some requiring only a strong magnetic field and rotation (Aubert and Finlay 2019; Gerick et al. 2020) while others rely on the presence of a possible stratified layer at the top of the core (Buffett and Matsui 2019). Though tempting, it may be dangerous to interpret such features that are at the limit of the present spatial resolution and temporal resolution (Gillet 2019). It will be important to assess whether such features remain coherent in the future, as the resolution of the SA increases.

Conclusions

In this article, we have studied global patterns and sub-decadal changes in geomagnetic secular variation during the past 20 years. We have shown that continuous coverage of magnetic field measurements from low-Earth orbiting satellite missions is now available during this period, provided one takes advantage of calibrated platform magnetometer data from the CryoSat-2 satellite. Using vector magnetic field measurements from Ørsted, CHAMP, Cryosat-2 and Swarm we have constructed Geomagnetic Virtual Observatory (GVO) time series that track sub-decadal changes in the core field secular variation at satellite altitude, and we have used the Subtractive Optimally Localized Averages (SOLA) technique to study the secular variation of the radial field and its time changes down at the core–mantle boundary. These are local methods whereby field measurements in the vicinity of the location of interest are combined so as to estimate the field at that point; measurements far away from the site of interest, have little or no influence on the field estimates.

Using the GVO method, we derived composite time series of geomagnetic secular variation, spanning nearly 20 years, on a global grid of 300 GVOs. GVO time series derived from IGRF-13 and CHAOS-6-x9 calibrated CryoSat-2 data, show similar sub-decadal SV features, and comparable levels of scatter. We found a scatter level for radial field SV of \(3.5\,\hbox {nT}/\hbox {yr}\) for the CryoSat-2 GVOs compared with \(2.5\,\hbox {nT}/\hbox {yr}\) for CHAMP and \(1.0\,\hbox {nT}/\hbox {yr}\) for Swarm. Comparing GVOs with overlapping epochs from 2014 to 2018, derived from calibrated CryoSat-2 data and Swarm data, we find similar sub-decadal SV changes, thus confirming the possibility of using CryoSat-2 for core field studies. In our 20-year-long composite GVO records, we observe fluctuations in the radial SV field of up to 20 nT/yr in amplitude occurring at low latitudes over time periods of 5–10 years. For instance, we see a rapid change of slope in the radial field SV over Indonesia around 2014, over South America and the South Atlantic region around 2007, 2011 and 2014, and in the central Pacific around 2017. Some of these events have previously been discussed in ground observatory records (Brown et al. 2013; Torta et al. 2015). They have the distinct “\(\Lambda\)” and “V”-shapes that are often associated with geomagnetic jerks, but as indicated in earlier investigations (Olsen and Mandea 2007, 2008; Chulliat and Maus 2014) we find these events are localized rather than global features.

Using the SOLA technique, we mapped the radial field SV directly at the core–mantle boundary using satellite data, computing spatially localized averages of the SV time-averaged over specified windows. Taking differences between consecutive SV estimates we obtained estimates for the secular acceleration (SA) at the CMB. Using only CryoSat-2 measurements we are able to successfully map the SA at the CMB, down to spatial averaging widths of \(\approx 42^{\circ }\), corresponding approximately to SH degree 10 or length scales of 2500 km at the CMB. Comparing SV and also SA field maps at the CMB, derived from CryoSat-2 and Swarm measurements, the same features can be identified having similar amplitudes and latitude/longitude extent.

In time–longitude plots of SOLA-based radial field SA estimates along the geographic equator at the CMB we find strong SA features, with amplitudes of \(\pm 2.5\mu \mathrm{T}/\mathrm{yr}^2\), under Indonesia at longitudes \(70^{\circ }\) to \(160^{\circ }\) from 2011 to 2014 during CryoSat-2 time, under central America at longitudes \(240^{\circ }\) to \(320^{\circ }\) from 2015 to 2019, and in sequences with alternating signs of radial SA under South America and the South Atlantic region at longitudes \(240^{\circ }\) to \(360^{\circ }\) during 2004-2019. The imaged SA features around 2013 at longitudes \(70^{\circ }\) to \(160^{\circ }\) and \(280^{\circ }\) to \(360^{\circ }\) demonstrates the usefulness of the CryoSat-2 measurements. Our results lend further support to a sign changing sequence of SA, for length scales down to 2500 km, at longitudes \(240^{\circ }\) to \(360^{\circ }\) from 2005 to 2019; this has been noticed in previous studies (Chulliat et al. 2015; Alken et al. 2020a). We have shown it is possible to increase the temporal resolution of SA estimates during Swarm era, compared to that seen for example in the CHAOS-7 model, by computing SA estimates from the differences of consecutive SOLA SV estimates derived using 1-yearly time windows. We find similar coherent structures as seen in TL plots constructed using 2-yearly time window, but also see that changes of sign in the SA occur within 1 year. In the central Pacific region at longitudes \(150^{\circ }\) to \(220^{\circ }\) we find strong positive and negative SA features appearing side-by-side in late 2017 that drift westwards until 2020. The results presented in Figs. 9 and 10, demonstrate that estimates of core field SA different from those found in the CHAOS model can be obtained despite employing similar data selection and external field modeling schemes. This is due to the important role played by the model parameterization and regularization in the SA recovered in the CHAOS model, especially for the small length-scales and fast time changes which are towards the limits of what can be reliably resolved from the data.

The rapid fluctuations of the core magnetic field described in this study are likely caused by time variations in the motions of the liquid metal outer core. In particular, changes in secular acceleration patterns at low latitudes provide constraints on the equatorial dynamics of the outer core (e.g. Aubert and Finlay 2019; Kloss and Finlay 2019). This is a topic of active research with various possible phenomenon recently proposed including equatorially trapped MAC waves in a stratified layer close to the core surface (Buffett and Matsui 2019; Chi-Durán et al. 2020) or the equatorial focusing of hydrodynamic waves driven by turbulent convection deep within the core (Gerick et al. 2020; Aubert and Finlay 2019).

It is undeniable that much core dynamics occurs on timescales either much longer, or much shorter, than can be can be resolved using the available satellite and ground observations (Gillet 2019). However, with 20 years of low-Earth orbit satellite measurements of the vector magnetic field, and with tools similar to those presented here, it is now possible to probe and characterize core field changes with increasing spatial and temporal resolution. We have shown here that platform magnetometer data can contribute to this activity, provided they are appropriately calibrated.