Assimilation of ground and satellite magnetic measurements: inference of core surface magnetic and velocity field changes–References

Assimilation of ground and satellite magnetic measurements: inference of core surface magnetic and velocity field changes

We jointly invert for magnetic and velocity fields at the core surface over the period 1997–2017, directly using ground-based observatory time series and measurements from the CHAMP and Swarm satellites. Satellite data are reduced to the form of virtual observatory time series distributed on a regular grid in space. Such a sequential storage helps incorporate voluminous modern magnetic data into a stochastic Kalman filter, whereby spatial constraints are incorporated based on a norm derived from statistics of a numerical geodynamo model. Our algorithm produces consistent solutions both in terms of the misfit to the data and the estimated posterior model uncertainties. We retrieve core flow features previously documented from the analysis of spherical harmonic field models, such as the eccentric anti-cyclonic gyre. We find enhanced diffusion patterns under both Indonesia and Africa. In contrast to a steady flow that is strong under the Atlantic hemisphere but very weak below the Pacific, interannual motions appear evenly distributed over the two hemispheres. Recovered interannual to decadal flow changes are predominantly symmetrical with respect to the equator outside the tangent cylinder. In contrast, under the Northern Pacific we find an intensification of a high latitude jet, but see no evidence for a corresponding feature in the Southern hemisphere. The largest flow accelerations that we isolate over the studied era are associated with meanders, attached to the equatorward meridional branch of the planetary gyre in the Eastern hemisphere, that is linked to the appearance of an eastward equatorial jet below the Western Pacific.

1 Introduction

Inferring information on the motions of the liquid outer core of the Earth requires properly separating the numerous sources of observed magnetic fields (geodynamo, crustal magnetisation, ionospheric and magnetospheric currents and their Earth induced counterparts). To circumvent some of the leakage issues, magnetic field models are often built using regularizations, to ensure spectral convergence of the core field and its time variations. This prevents a proper assessment of a posteriori errors on model coefficients. When these are used as data in reconstructions of the core dynamics, it can lead to biased estimates. Furthermore, by proceeding in successive steps (to a field model and then on to the core flow), one loses information.

From the early 1990’s alternative avenues of research arose, through which field models were built under topological constraints derived from physical insights. Constable et al.(1993)Constable, Parker, & Stark and O’Brien et al.(1997)O’Brien, Constable, & Parker proposed algorithms to apply, on single epoch pairs of models, magnetic flux conservation conditions at the core-mantle boundary (CMB) that are appropriate assuming that magnetic diffusion is negligible. Along the same lines, Jackson et al.(2007)Jackson, Constable, Walker, & Parker added a constraint on the radial vorticity. They showed that it was possible for a magnetic model to satisfy both these topological conditions, and the constraint from magnetic observations, from the late 19th century onwards.

Conversely, Chulliat & Olsen(2010) tested the validity of the frozen flux hypothesis using data from Magsat, Oersted and CHAMP satellite missions. They found an increase of the data misfit in some areas, potentially suggesting local failures of the constraint. Such studies motivated the co-estimation, from magnetic observations, of both the field and the flow, imposing with a weak formalism the frozen flux radial induction equation at the CMB (Lesur et al.(2010)Lesur, Wardinski, Asari, Minchev, & Mandea; Wardinski & Lesur(2012)). They concluded that the frozen flux constraint remained compatible with ground-based and satellite magnetic records. Pursuing an alternative approach, Beggan & Whaler(2009) and Whaler & Beggan(2015) obtained piecewise constant or linear flow models directly from magnetic data (see also Whaler et al.(2016)Whaler, Olsen, & Finlay).

One limitation though of such approaches is related to the uncertainties associated with the large scale induction equation itself (and associated null-flux curves), assuming models truncated at spherical harmonic degree $n\simeq 13$ (Gillet et al.(2009)Gillet, Pais, & Jault). Subgrid-scale effects arising due to the non-linear induction process (e.g. Eymin & Hulot(2005); Pais & Jault(2008); Gillet et al.(2009)Gillet, Pais, & Jault; Baerenzung et al.(2016)Baerenzung, Holschneider, & Lesur) turn out to be the main source of uncertainty in the recovery of core surface flows from modern geomagnetic records. Barrois et al.(2017)Barrois, Gillet, & Aubert – hereafter referred as BGA17 – illustrate how ignoring their impact leads to severely biased flow models (see also Baerenzung et al.(2017)Baerenzung, Holschneider, Wicht, Sanchez, & Lesur, on the reliability of core flow reconstructions).

BGA17 furthermore show from the analysis of geodynamo simulations that magnetic diffusion at the core surface, enslaved to poloidal flow below the CMB, affects the recorded field changes at all time-scales including rapid changes. This may seem at odds with the often used assumption of negligible magnetic diffusion that follows the argument of a high magnetic Reynolds number for large-scale motions in the core (see Holme(2015)).

In the present work we invert, from magnetic field observations collected at and above the Earth’s surface, for both the magnetic and velocity fields at the core surface, taking into account both magnetic diffusion and subgrid induction. We merge spatial information provided by numerical simulations, specifically from the Coupled Earth dynamo (CED) model (Aubert et al.(2013)Aubert, Finlay, & Fournier) and temporal constraints coming from a restriction of the field evolution to a chosen class of stochastic process. The sequential algorithm of BGA17, which considers as input data time series of spherical harmonic coefficients of the main field, is extended to account for both virtual observatory (Mandea & Olsen(2006)) and ground observatory time series that cover the period 1997–2017. Our approach has similarities with the previous works of Gillet et al.(2015a)Gillet, Barrois, & Finlay and Baerenzung et al.(2016)Baerenzung, Holschneider, & Lesur, which favoured flat flow spatial spectra at the CMB, since the spatial dynamo norm employed here departs from the norms often employed to ensure spectral convergence. In addition, our stochastic framework allows us to discuss posterior model errors for both the flow and the magnetic field.

The paper is organised as follows. In section §2 we describe the ground-based observatory data and satellite-based virtual observatory data, and the methodology we follow to recover magnetic and velocity fields at the CMB. In section §3.1, we present our resulting geomagnetic model and its associated uncertainties, before we analyse in §3.2 our core flow solutions. Finally, implications for our understanding of the core dynamics and possible further improvements for the algorithm are given in section §4.

2 Methodology

2.1 Ground-based and virtual observatory data

2.1.1 Ground observatory data

We use magnetic measurements made at 186 ground observatories (GO) covering the period 1997–2017. Hourly mean values are taken from the BGS database, version 0111, using Intermagnet and WDC Edinburgh data as available in May 2017. The data have been checked and corrected for known baseline jumps (Macmillan & Olsen(2013)). ‘Revised monthly means’ were then derived from these hourly means, following the procedure described by Olsen et al.(2014)Olsen, Lühr, Finlay, Sabaka, Michaelis, Rauberg, & Tøffner-Clausen. Briefly, predictions of the large-scale magnetospheric field (and the associated induced field) from the CHAOS-6 field model, as well as predictions for the ionospheric Sq field (and the associated induced field) from the CM4 model (Sabaka et al.(2004)Sabaka, Olsen, & Purucker) are subtracted from the hourly mean values, and then robust (Huber-weighted) monthly mean values are computed using an iterative-reweighting procedure. Annual differences of such revised monthly means are routinely used in deriving the CHAOS series of field models and in order to study high resolution secular variation since, compared with simple monthly means, they are less contaminated by external field effects. Here, since we also wish to use the field itself for model construction, the median difference between each series and CHAOS-6 predictions was removed, in order to account in a simple way for the bias due to unmodelled crustal fields. In order to obtain the same sampling rate as that adopted for the virtual observatory series described below, the revised monthly mean series were finally averaged over 4 months windows to obtain the GO series used in our data assimilation scheme.

2.1.2 Virtual observatory data

In addition to GO data, we make use of satellite measurements from the CHAMP and Swarm missions covering respectively 2000–2010 and 2014–2017, through so-called virtual observatory (VO) data (Mandea & Olsen(2006); Olsen & Mandea(2007)). These provide a regular spatial and temporal sampling of the global field, convenient for our Kalman filter algorithm (detailed in §2.2) and involve estimates from an easily manageable number of locations, which has computational advantages.

VO data were computed using measurements collected by the CHAMP vector field magnetometer between July 2000 and September 2010 and from the Swarm vector field magnetometers, onboard all three satellites (Alpha, Bravo, Charlie), between January 2014 and April 2017. Starting from the CHAMP MAG-L3 and Swarm Level 1b MAG-L, version 0501, data products, we sub-sampled at 15s intervals the data in the vector field magnetometer (VFM) frame. Using the Euler rotation angles as given by the CHAOS-6-x3 model (which was based on Swarm and ground observation data up until April 2017), we rotated the VFM data into an Earth-Centered Earth-Fixed (ECEF) coordinate frame.

Measurements from known problematic days were removed, for instance where satellite manoeuvres happened. Furthermore, gross data outliers with deviations more than 500 nT from CHAOS-6-x3 field model predictions were rejected. Based on previous studies of VO estimates (e.g., Beggan et al.(2009)Beggan, Whaler, & Macmillan), we then employed data selection criteria retaining only data for which:

1. the sun was at maximum $10^{\circ}$ above horizon;

2. geomagnetic activity index $\mathrm{K_{p}}<3^{\mathrm{o}}$;

3. the $RC$ disturbance index (Olsen et al.(2014)Olsen, Lühr, Finlay, Sabaka, Michaelis, Rauberg, & Tøffner-Clausen) had $|d\mathrm{RC}/dt|<3$ nT/h;

4. merging electric field at the magnetopause $\mathrm{E_{m}}\leq 0.8$ mV/m, with $E_{m}=0.33v^{4/3}B_{t}^{2/3}\mathrm{sin}(|\Theta|/2)$. $v$ is the solar wind speed, $\Theta=\mathrm{arctan}(B_{y}/B_{z})$ and $B_{t}=\sqrt{B_{y}^{2}+B_{z}^{2}}$. $B_{y}$ and $B_{z}$ are components of the interplanetary magnetic field (IMF) in the geocentric solar magnetospheric (GSM) coordinate system, calculated using 2 hourly means of 1-min values of the IMF and solar wind extracted from the OMNI database;

5. IMF $B_{z}>0$ nT and IMF $|B_{y}|<10$ nT, again based on 2 hourly mean of 1 minute values.

Following this data selection, estimates of the fields due to various unmodelled sources were next removed from the data:

1. the magnetospheric and its induced fields as given by the CHAOS-6-x3 model;

2. the ionospheric and its induced fields as given by the CM4 model (Sabaka et al.(2004)Sabaka, Olsen, & Purucker);

3. the static internal field for spherical harmonic degrees $n>20$ given by the CHAOS-6-x3 model.

Although imperfect, in our opinion it is more consistent to remove such estimates rather to ignore known field sources.

Based on this data we then carried robust inversions for time-averaged point estimates (i.e. VOs) using data windows of 4 months width (60 days each side of an epoch $t_{j}$). In order to aid the robust inversion procedure in identifying and downweighting outliers, following Olsen & Mandea(2007) as a pre-processing step, we also removed a time-dependent internal field, here taken from the CHAOS-6-x3 model (Finlay et al.(2016b)Finlay, Olsen, Kotsiaros, Gillet, & Tøffner-Clausen), for spherical harmonic degrees 1 to 20, within each four month window. The CHAOS-6x-3 prediction at the target point and time was then added back at the end of the analysis. Note that this does not prevent our 4-monthly VO series, and the derived SV series from departing from CHAOS-6x-3; information about the time-dependence within each 4 month window is however lost.

We assume that the residual field $\tilde{\bf B}$, after the removal of the time-dependent internal field from the CHAOS-6-x3, can be represented as the gradient of a scalar potential $V$, i.e.

$$\tilde{\bf B}=-\nabla V\,.$$ (1)

The residual field and associated positions are transformed into a local Cartesian coordinate system with origin at the VO points of interest, with $x$ pointing towards geographic South, $y$ pointing towards East and $z$ pointing upwards. We use an expansion of the local potential up to cubic terms. Because the geomagnetic field is irrotational ($\nabla\times\tilde{\bf B}=0$) and solenoidal ($\nabla\cdot\tilde{\bf B}=0$), this local potential is entirely determined by 15 independent parameters:

$$V(x,y,z)=v_{x}x+v_{y}y+v_{z}z+v_{xx}x^{2}+v_{yy}y^{2}-(v_{xx}+v_{yy})z^{2}$$ (2)

$$+2v_{xy}xy+2v_{xz}xz+2v_{yz}yz-(v_{xyy}+v_{xzz})x^{3}$$

$$+3v_{xxy}x^{2}y+3v_{xxz}x^{2}z+3v_{xyy}xy^{2}+3v_{xzz}xz^{2}+6v_{xyz}xyz$$

$$-(v_{xxy}-v_{yzz})y^{3}+3v_{yyz}y^{2}z+3v_{yzz}yz^{2}-(v_{xxz}+v_{yyz})z^{3}\,.$$

For each VO position vector ${\bf r}_{k}=(\theta_{k},\phi_{k},r_{k})$ and at epoch $t_{j}$, all data positioned within a cylinder of radius 850km ($\approx 7.5^{\circ}$) of the VO target ${\bf r}_{k}$, and within 60 days either side of $t_{j}$ were used to build a local data vector ${\bf d}^{k,j}$. These data are then related to the 15 parameters defining the VO potential model ${\bf m}_{vo}^{k,j}$ at that site and epoch via ${\bf d}^{k,j}=\underline{\underline{{\bf g}}}^{k,j}{\bf m}_{vo}^{k,j}$, where the elements of the matrix $\underline{\underline{{\bf g}}}^{k,j}$ are determined from (1) and (2).

Rather than working directly with ${\bf d}^{k,j}$ in deriving ${\bf m}_{vo}^{k,j}$ we make use of along-track and East-West (using Swarm Alpha and Charlie only) sums and differences of the magnetic field components. An advantage of using field differences is that these have a reduced sensitivity to large-scale external signals, although data sums also need to be included in order to ensure sufficient information on the longer wavelengths core field. Using sums and differences has been found advantageous in a number of other field modelling efforts (Sabaka et al.(2015)Sabaka, Olsen, Tyler, & Kuvshinov; Olsen et al.(2015)Olsen, Hulot, Lesur, Finlay, Beggan, Chulliat, Sabaka, Floberghagen, Friis-Christensen, Haagmans, et al.). We calculate along-track (AT) sums ($\Sigma$) and differences ($\Delta$) as

$$\left\{\begin{array}[]{rl}\Sigma d_{i}^{AT}=&[\tilde{B}_{i}(\mathbf{r},t)+\tilde{B}_{i}(\mathbf{r}+\delta\mathbf{r},t+15\mathrm{s})]/2\\ \Delta d_{i}^{\mathrm{AT}}=&[\tilde{B}_{i}(\mathbf{r},t)-\tilde{B}_{i}(\mathbf{r}+\delta\mathbf{r},t+15\mathrm{s})]\end{array}\right.\,.$$ (5)

$\tilde{B}_{i}={\mathbf{1}_{i}}\cdot\mathbf{\tilde{B}}(\mathbf{r})$ are the residual magnetic field components in spherical polar coordinates (where $i=r,\theta$ or $\phi$, and ${\mathbf{1}_{i}}$ are unit vectors). The East-West cross-track (CT) sums and differences between are calculated as

$$\left\{\begin{array}[]{rl}\Sigma d_{i}^{\mathrm{CT}}=&[\tilde{B}_{i}^{\mathrm{Alpha}}(\mathbf{r}_{1},t_{1})+\tilde{B}_{i}^{\mathrm{Charlie}}(\mathbf{r}_{2},t_{2})]/2\\ \Delta d_{i}^{\mathrm{CT}}=&[\tilde{B}_{i}^{\mathrm{Alpha}}(\mathbf{r}_{1},t_{1})-\tilde{B}_{i}^{\mathrm{Charlie}}(\mathbf{r}_{2},t_{2})]\end{array}\right.\,.$$ (8)

Here, for a given orbit of Alpha we select the corresponding Charlie data to be the one closest in colatitude such that $|\delta t|=|t_{1}-t_{2}|<50s$. Crucially, in order to relate these sums and differences to the VO model parameters, we also take sums and differences of the elements of the design matrices $\underline{\underline{{\bf g}}}^{k,j}$ associated with the predictions of the VO model for the field components at the individual data locations. This results in a design matrix

$$\underline{\underline{{\bf G}}}^{k,j}=\begin{bmatrix}\underline{\underline{\Sigma{\bf g}}}^{k,j}\\ \underline{\underline{\Delta{\bf g}}}^{k,j}\end{bmatrix}$$ (9)

associated with the data vector ${\bf D}^{k,j}=\left[\Sigma{\bf d}^{k,j}\,\,\Delta{\bf d}^{k,j}\right]^{T}$. In this way we fully account for the change in the unit vectors associated with the two locations contributing to the sums and differences when deriving the parameters ${\bf m}_{vo}^{k,j}$. The inversion for each ${\bf m}_{vo}^{k,j}$ is carried out via a robust Huber weighted least-square fit

$${\bf m}_{vo}^{k,j}=\left[(\underline{\underline{{\bf G}}}^{k,j})^{T}{\bf W}\underline{\underline{{\bf G}}}^{k,j}\right]^{-1}(\underline{\underline{{\bf G}}}^{k,j})^{T}{\bf D}^{k,j}$$ (10)

where ${\bf W}$ is a diagonal vector of Huber weights that ensures a robust solution (Olsen(2002); Sabaka et al.(2004)Sabaka, Olsen, & Purucker) and are iteratively updated until convergence. Once ${\bf m}_{vo}^{k,j}$ is determined, the three field components at the site and epoch of interest, $\tilde{\bf B}_{k}({\bf r}_{k},t_{j})=-\nabla V_{k}({\bf r}_{k},t_{j})$, are computed and added back on to the CHAOS-6-x3 prediction for the internal field (for degrees 1-14 only, to avoid as far as possible the lithospheric field) at the target location.

We constructed VO estimates at $P_{\mathrm{VO}}=200$ locations, with a spacing of about 1600 km ($\approx 14^{\circ}$, see dots in Figure 1), located in an approximately equal area grid based on the spherical surface partition algorithm of Leopardi(2006). The altitude of the VOs are 300km and 500km during the CHAMP and Swarm periods, respectively. Using predictions of the three components ($B_{r}$, $B_{\theta}$, $B_{\phi}$) of the magnetic field at $P_{\mathrm{VO}}$ locations, we finally obtain $3P_{\mathrm{VO}}$ time series (i.e. one point every 4 months during CHAMP and Swarm times, 48 epochs in all), stored in a vector ${\bf y}_{\mathrm{VO}}(t)$. The SV was computed as annual differences of the 4 month time series.

2.1.3 Uncertainty estimates for the GO and VO series

In order to obtain as much information as possible from the GO and VO data, while at the same time seeking to avoid overfitting them, it is important that appropriate uncertainty estimates are specified for each time series. We define ${\bf C}_{\mathrm{GO}}$ and ${\bf C}_{\mathrm{VO}}$ to be the measurement error cross-covariance matrices for GO and VO data at each epoch, of sizes respectively $3P_{GO}\times 3P_{GO}$ and $3P_{VO}\times 3P_{VO}$. Data errors are supposed independent of time. Different data uncertainties are assigned for the VO’s derived from CHAMP and Swarm respectively.

Regarding the GO time series described above, we follow a similar approach to that used in CHAOS field model series (Olsen et al.(2014)Olsen, Lühr, Finlay, Sabaka, Michaelis, Rauberg, & Tøffner-Clausen; Finlay et al.(2016b)Finlay, Olsen, Kotsiaros, Gillet, & Tøffner-Clausen) and derive uncertainty estimates as follows. A three-by three covariance matrix was computed for each observatory location from the time-series of the three components, after removing the predictions of the CHAOS-6 field model and de-trending. The square root of the diagonal elements of these covariance matrices were taken to be the uncertainty estimates for each component at each observatory. The same procedure was applied to both the MF and SV series.

For consistency, a very similar procedure was also applied to the VO series in order to obtain their uncertainty estimates. For each VO location, covariances were calculated between the time series of the three components (after removing from each series the predictions of the CHAOS-6 model and de-trending), in order to obtain a three-by three covariance matrix. A robust procedure for calculating the covariances (using the Minimum Covariance Determinant estimator, Verboven & Hubert(2005)) was employed. However, only the square-root of the diagonal elements of the covariance matrices were taken to be the uncertainty estimates for each series, with similar procedures applied to both MF and SV series. To illustrate the range of the adopted uncertainty estimates, we show in Figure 1 the r.m.s. SV uncertainty estimates for all locations where data (GO or VO) are used in this study.

Figure 1: SV observation error estimates (colorscale in $nT/y$) at all location where GOs (hexagons) and VOs (circles) are used in this study, for the 3 components $\dot{B}_{r}$, $\dot{B}_{\theta}$ and $\dot{B}_{\phi}$ (from top to bottom). The size of the markers is proportional to the magnitude of the a-priori error estimates.

Note that by using only the diagonal elements of ${\bf C}_{\mathrm{GO}}$ and ${\bf C}_{\mathrm{VO}}$ we effectively consider the errors on each GO and VO series to be uncorrelated with the errors on other series. In reality errors between components and between series will be correlated. This can be taken into account using full (i.e. dense) covariance matrices. It is however challenging to estimate cross-covariances for matrices of size larger than the length of the contributing times series (consisting of one sample every four months). We therefore postpone this step to future studies. Instead, by restricting to only 200 VO locations and ensuring that there was little overlap between the VO search radii we reduce as far as possible the correlations between distinct VO series.

Finally, we concatenate the above GO and VO main field data vectors for each epoch into ${\bf y}^{o}(t)=[{\bf y}_{\mathrm{VO}}^{T}\,{\bf y}_{\mathrm{GO}}^{T}]^{T}$. The associated observation errors covariance matrix ${\bf R}_{yy}$, of rank $P=3P_{\mathrm{VO}}+3P_{\mathrm{GO}}$, is thus derived from the diagonals of ${\bf C}_{\mathrm{VO}}$ and ${\bf C}_{\mathrm{GO}}$. In the next section we will consider both main field and secular variation data. SV data $\dot{\bf y}^{o}(t)$ are computed as annual differences of the four monthly (GO or VO) series. We follow the same approach as above to estimate the SV data errors variances (shown in Figure 1) that are stored in a diagonal matrix ${\bf R}_{\dot{y}\dot{y}}$ of rank $P$.

2.2 Re-analysis of GO and VO data ground and satellite magnetic observations

The assimilation algorithm used in the present study is essentially the one derived by BGA17 (see their table 2 for a summary). It is a sequential tool, consisting of a succession of forecast and analysis steps. The main modifications concern the direct integration of observations at and above the Earth’s surface, while BGA17 considered data in the form of MF and SV spherical harmonic coefficients. We begin by recalling the main points of our stochastic forecast model, before we go on to describe the changes implemented in the present study regarding the analysis step. These essentially concern the observation operator linking the state variables to the observations.

2.2.1 Stochastic forecast model

We forecast the evolution of the radial magnetic field, $B_{r}$, at the CMB using the radial component of the induction equation, written as

$$\frac{\partial\overline{B}_{r}}{\partial t}=-\overline{\nabla_{h}\cdot\left({\bm{u}}_{H}\overline{B}_{r}\right)}+e_{r}+d_{r}({\bm{u}}_{H},\overline{B}_{r})\,,$$ (11)

where overlines mean the projection onto large length-scales. $e_{r}$ stands for the subgrid induction processes arising due to the unresolved magnetic field at small length-scales, ${\bm{u}}_{H}$ is the horizontal flow, and $d_{r}$, enslaved to $\overline{B}_{r}$ and ${\bm{u}}_{H}$, approximates the radial component of the diffusion operator (see below). The evolutions of $e_{r}$ and ${\bm{u}}_{H}$ are governed by order one auto-regressive stochastic processes,

$$\frac{\mathrm{d}e_{r}}{\mathrm{d}t}+\frac{e_{r}}{\tau_{e}}=\zeta_{e}\,,$$ (12)

$$\frac{\mathrm{d}{\bm{u}}_{H}}{\mathrm{d}t}+\frac{({\bm{u}}_{H}-\hat{\bm{u}}_{H})}{\tau_{u}}=\zeta_{u}\,,$$ (13)

with $\zeta_{e}$ and $\zeta_{u}$ white noise processes, and $\hat{\bm{u}}_{H}$ the background flow model (obtained as the time-averaged flow from the CED model). These processes come from the same family of process as employed by Baerenzung et al.(2017)Baerenzung, Holschneider, Wicht, Sanchez, & Lesur. For each process, an effective restoring force is implemented via single time scales that we respectively fix as $\tau_{e}=10$ yrs and $\tau_{u}=30$ yrs. Spatial cross-covariances of the two above fields are derived from statistics of a free run of the CED (Aubert et al.(2013)Aubert, Finlay, & Fournier).

The advected fields $e_{r}$, ${\bm{u}}_{H}$, $\overline{B}_{r}$ and $d_{r}$ are represented through spherical harmonics, whose coefficients are stored in vectors ${\bf e}(t)$, ${\bf u}(t)$, ${\bf b}(t)$ and ${\bf d}(t)$, respectively. Diffusion in equation (11), and its dependence on $e_{r}$ and ${\bf u}_{H}$, is also an expression of cross-covariances extracted from the CED (involving the radial magnetic field below the CMB). The projection onto large length-scales is processed in the spectral domain, restricting the induction equation (and thus the expansion of the fields $e_{r}$, $\overline{B}_{r}$ and $d_{r}$) to spherical harmonic degrees $n\leq n_{b}=14$, while the velocity field is truncated at $n_{u}=18$. We write as $\dot{\bf b}(t)$ the vector of SV spherical harmonic coefficients.

2.2.2 Integrating ground and satellite data in the assimilation tool

We write as ${\bf M}$ the operator that links the vector ${\bf b}(t)$ to the three components main field observations ${\bf y}(t)$ in the spatial domain (e.g. Olsen et al.(2010)Olsen, Glassmeier, & Jia):

$${\bf y}(t)={\bf M}{\bf b}(t)\,.$$ (14)

At each epoch it is of size $n_{o}\times n_{b}(n_{b}+2)$, with $n_{o}=3(P_{\mathrm{VO}}+P_{\mathrm{GO}})$ the size of the observation vector. The matrix ${\bf M}$ is composed of sub-matrices ${\bf M}_{r}$, ${\bf M}_{\theta}$ and ${\bf M}_{\phi}$, depending on the considered component of the magnetic field. In practice, elements of the matrix are, for a column $j$ corresponding to a coefficient $g_{n_{j}}^{m_{j}}$, and a line $i$ to an observation at a coordinate ${\bf r}_{i}=(r_{i},\theta_{i},\phi_{i})$,

$${\bf M}_{ri,j} = (n_{j}+1)\left(\dfrac{a^{\oplus}}{r_{i}}\right)^{n_{j}+2}\mathscr{P}_{n}^{m}(\theta_{i})\cos(m_{j}\phi_{i})\,,$$ (15)

$${\bf M}_{\theta i,j} = \left(\dfrac{a^{\oplus}}{r_{i}}\right)^{n_{j}+2}\dfrac{d\mathscr{P}_{n}^{m}(\theta_{i})}{d\theta}\cos(m_{j}\phi_{i})\,,$$ (16)

$${\bf M}_{\phi i,j} = \left(\dfrac{a^{\oplus}}{r_{i}}\right)^{n_{j}+2}\dfrac{m_{j}\mathscr{P}_{n}^{m}(\theta_{i})}{\sin(\theta_{i})}(-1)\sin(m_{j}\phi_{i})\,.$$ (17)

For a line $j$ corresponding to a coefficient $h_{n_{j}}^{m_{j}}$, the function $\sin$ replaces $\cos$ in (15) and (16), and $\cos$ replaces $(-1)\sin$ in (17). $a^{\oplus}=6371.2$ km is the Earth’s spherical reference radius and $\mathscr{P}_{n}^{m}$ are the Legendre polynomials.

The analysis in the Kalman filter algorithm employed by BGA17 consists of two steps: first an analysis of the vector ${\bf b}$ containing MF spherical harmonic coefficients from MF spherical harmonic coefficients data, and second an analysis of the vector ${\bf z}$ (that concatenates ${\bf u}$ and ${\bf e}$) from SV spherical harmonic coefficients data. Writing as ${\bf P}_{bb}^{f}$ the forecast model covariance matrix for ${\bf b}$, the first analysis (equation (19) of BGA17) is replaced here by

$$\forall k\in[1,N_{m}],\;$$ (18)

$${\bf b}^{ka}(t_{a})={\bf b}^{kf}(t_{a})+{\bf P}_{bb}^{f}{\bf M}^{T}\left[{\bf M}{\bf P}_{bb}^{f}{\bf M}^{T}+{\bf R}_{yy}\right]^{-1}$$

$$\cdot\left({\bf y}^{ko}(t_{a})-{\bf M}{\bf b}^{kf}(t_{a})\right)\,,$$

with $t_{a}$ the analysis epoch and the superscript $k$ referring to the $k^{th}$ realization within an ensemble chosen to be of size $N_{m}=50$. Writing as ${\bf P}_{zz}^{f}$ the forecast model covariance matrix for ${\bf z}$, the second analysis (equation (20) of BGA17) is replaced here by

$$\forall k\in[1,N_{m}],\;$$ (19)

$${\bf z}^{ka}(t_{a})={\bf z}^{kf}(t_{a})+{\bf P}_{zz}^{f}{\bf G}^{kT}\left[{\bf G}^{k}{\bf P}_{zz}^{f}{\bf G}^{kT}+{\bf R}_{\dot{y}\dot{y}}\right]^{-1}$$

$$\cdot\left(\delta\dot{\bf y}^{ko}(t_{a})-{\bf G}^{k}{\bf z}^{kf}(t_{a})\right)\,,$$

where the new observation operator is ${\bf G}^{k}={\bf M}{\bf H}({\bf b}^{ka})$, with ${\bf H}$ as defined in BGA17. Here $\delta\dot{\bf y}^{ko}(t_{a})=\dot{\bf y}^{ko}(t_{a})-{\bf M}{\bf d}^{kf}(t_{a})$ are the direct SV observations corrected by the forecast contribution from diffusion to the radial induction equation. This latter is sought iteratively at each analysis step, as in BGA17. Note that we consider an ensemble of observations ${\bf y}^{o}$ and $\dot{\bf y}^{o}$, which are perturbed by random noise according to respectively ${\bf R}_{yy}$ and ${\bf R}_{\dot{y}\dot{y}}$. We recall that we consider in equations (18–19) forecast covariance matrices ${\bf P}_{zz}$ and ${\bf P}_{bb}^{f}$ that are frozen throughout the re-analysis period. These are derived directly from the CED cross-covariances on ${\bf b}$, ${\bf u}$ and ${\bf e}$ spherical harmonic coefficients, involving scaling prefactors obtained analytically from the stochastic model presented in §2.2.1 (see BGA17 for details). For comparison, Baerenzung et al.(2017)Baerenzung, Holschneider, Wicht, Sanchez, & Lesur employ a full implementation of the Ensemble Kalman filter (Evensen(2003)), i.e. they update the cross-covariances at each analysis step, requiring many more realizations to obtain well-conditioned matrices.

Finally, an extra complexity arises because the number of observation sites changes over time. Indeed, because of the selection criteria, the number of satellite data available may not always be sufficient to make a reliable VO estimate. Under these conditions the VO data point is considered to be absent: the associated elements of the data vector ${\bf y}^{o}(t)$ at a given time $t$ are removed, together with the corresponding lines and columns of ${\bf R}_{yy}$, and the corresponding lines of the matrix ${\bf M}$ (and thus ${\bf G}$). This procedure is performed during each analysis. Thus, the size $P$ of the data vector changes through time, reflecting the changing number of available satellite observations through time (see Figure 2).

Figure 2: Time evolution of the number of SV data points (VOs in red, GOs in blue).

To summarise, in this study we work with predictions made by spherical harmonic coefficients that are projected in physical space, where they are adjusted during the analysis step according to the observations and the covariance matrices. As such, our algorithm is still based almost entirely on the spectral domain; only the analysis steps are performed in physical space, in order to match the observed magnetic field data. Notice that we corrected for two mistakes in the implementation of the algorithm by BGA17: a sign error in the background flow $\hat{\bf u}$, and off-diagonal elements of the covariance matrix for ${\bf e}$ were non intentionally ignored. Performing comparisons between re-analyses before and after correction, we found two consequences: a reduction of the dispersion within the ensemble of realizations, and an (almost stationary) shift in the analysed diffusion for some coefficients (including the axial dipole, see §3.1.2). This latter is almost entirely compensated by a shift in the analysed $e_{r}$, with minor impact on the recovered flow. Otherwise, the qualitative conclusions of BGA17 remain unaltered.

2.3 Posterior diagnostics

We now define several diagnostics used to evaluate the quality and the consistency of our results. We shall compare a quantity ${\bm{x}}$ (MF, SV, subgrid error, diffusion… in the spatial or spectral domain) with observations ${\bm{x}}^{o}$ (when available), or with the same quantities ${\bm{x}}^{c}$ from the CHAOS-6 geomagnetic model (Finlay et al.(2016b)Finlay, Olsen, Kotsiaros, Gillet, & Tøffner-Clausen). We define its time average

$$\hat{\bm{x}}=\frac{1}{t_{f}-t_{i}}\int_{t_{i}}^{t_{f}}{\bm{x}}(t)dt\,,$$ (20)

with $t_{i}$ and $t_{f}$ the initial and final epochs, its ensemble mean

$$\left<{\bm{x}}(t)\right>=\dfrac{1}{N_{m}}\sum_{k=1}^{N_{m}}{\bm{x}}^{k}(t)\,,$$ (21)

the dispersion within the ensemble

$${\bm{\sigma}}_{x}(t)=\displaystyle\sqrt{\dfrac{1}{N_{m}-1}\sum_{k=1}^{N_{m}}\left({\bm{x}}^{k}(t)-\left<{\bm{x}}(t)\right>\right)^{2}}\,,$$ (22)

and finally the bias between our ensemble mean model and the reference ${\bm{x}}^{c}$,

$${\bm{\delta}}_{x}(t)=\displaystyle{\bm{x}}^{c}-\left<{\bm{x}}(t)\right>\,.$$ (23)

We also define spatial power spectra of any magnetic trajectory ${\bf b}(t)$ as

$${\cal R}_{b}(n,t)=(n+1)\left(\frac{a^{\oplus}}{c}\right)^{2n+4}\sum_{m=0}^{n}\left[{g_{n}^{m}(t)}^{2}+{h_{n}^{m}(t)}^{2}\right]\,,$$ (24)

with similar notations for $\dot{\bf b}(t)$, ${\bf d}(t)$ and ${\bf e}(t)$. $c=3485$ km is the Earth’s core radius, and $g_{n}^{m}$ and $h_{n}^{m}$ are Schmidt semi-normalised spherical harmonic coefficients for the magnetic field at the Earth’s surface. Finally, the spatial power spectrum for core flow trajectories ${\bf u}$ writes

$${\cal S}(n,t)=\frac{n(n+1)}{2n+1}\sum_{m=0}^{n}\left[{{t_{c}}_{n}^{m}(t)}^{2}+{{t_{s}}_{n}^{m}(t)}^{2}+{{s_{c}}_{n}^{m}(t)}^{2}+{{s_{s}}_{n}^{m}(t)}^{2}\right]\,,$$ (25)

with ${{t_{c,s}}_{n}^{m}}$ and ${{s_{c,s}}_{n}^{m}}$ Schmidt semi-normalised spherical harmonic coefficients for the toroidal and poloidal components of the flow. We also define the flow norm

$${\cal N}=\sum_{n=1}^{n_{u}}\frac{n(n+1)}{2n+1}\sum_{m=0}^{n}\left[{{t_{c}}_{n}^{m}}^{2}+{{t_{s}}_{n}^{m}}^{2}+{{s_{c}}_{n}^{m}}^{2}+{{s_{s}}_{n}^{m}}^{2}\right]\,.$$ (26)

The above power spectra can be considered for the ensemble mean or the dispersion within the ensemble, in which case they are respectively noted ${\cal R}_{<x>}(n,t)$ and ${\cal R}_{\delta x}(n,t)$. Additionally, all those quantities may be averaged in time and/or computed only at analysis periods. For example, the time-averaged spatial power spectrum of the dispersion of magnetic field solutions at analysis epochs is $\hat{\cal R}_{\delta b}^{a}(n)$. The same convention as above holds for core flow spectra.

3 Results

We apply our algorithm to VO and GO magnetic field observations over a period spanning from $t_{i}=1996.92$ to $t_{f}=2016.92$. We recall that since we use satellite measurements from CHAMP and Swarm missions, VOs are available only over the period 2000–2010 and 2014–2017, whereas GOs are available over the whole time-span. Analysis are performed every $\Delta t^{a}=4$ months. The sequences of analyses and forecasts between 1997 and 2001 are used to warm-up the filter (see Figure 7 in BGA17), avoiding an increase in the ensemble spread over the first years of the targeted satellite era. This warm-up period is not considered below when interpreting the ensemble of inverted magnetic field and flow. We first describe predictions from our re-analysis for observations in the physical domain (§3.1.1), before we present the resulting magnetic model (§3.1.2), and insights on core flows over various time-scales (§3.2).

3.1 Geomagnetic field models

3.1.1 Predictions for GO and VO series

We compare in Figure 3 our series of SV forecasts and analysis with two examples of observation series (one VO and one GO), and with the predictions from CHAOS-6. The large spread of the SV forecasts is to be expected given the large uncertainties associated with subgrid errors and the large-scale flow (see BGA17). At both sites, the dispersion within the ensemble of SV trajectories encompasses most of the time the observations. Moreover the predictions from CHAOS-6 and from our ensemble of SV models are generally consistent. Our algorithm thus seems able to provide a coherent estimate of the SV probability density function (PDF) at the Earth’s surface and at satellite altitude. In addition, we highlight that even during the period 2010-2014 where no VO data are available, the trajectory of SV model, controlled by the stochastic prior and GO data only, remains reasonable, with a slight increase in the ensemble spread that always contains CHAOS-6. Note that our algorithm tends to drive the system toward low SV values (see the saw-tooth patterns in Figure 3). This feature is to be expected given our choice of the stochastic models for ${\bm{u}}_{H}$ and $e_{r}$, which control the evolution of the SV. In the absence of data constraints, the process will drift back the ensemble average trajectories for ${\bm{u}}_{H}$ and $e_{r}$ towards the average dynamo state, which by construction is responsible for a weak SV. This is not a major drawback as soon as we analyse frequently enough, though it does limit the prediction capabilities of our tool (as discussed in BGA17).

Figure 3: SV time series for the three components $(dB_{r}/dt,dB_{\theta}/dt,dB_{\phi}/dt)$, at one VO location $\left\{r=6671\;km,\theta=90^{\circ},\phi=88,8^{\circ}\right\}$ (top), and at Chambon-la-forêt $\left\{r=6366\;km,\theta=42^{\circ},\phi=2^{\circ}\right\}$ (bottom). SV observations are shown in black, CHAOS-6 predictions in green, predictions from our analysis in red. The shaded area correspond to $\pm\sigma_{\dot{b}}$, see equation (22).

We check in Figure 4 the accuracy with which our model fits MF and SV observations, with the histograms of the prediction errors (over all analyses) normalised to the observation errors, for the three components of the magnetic field. Concerning the MF, prediction errors are only weakly biased, excepted for $B_{\theta}$ (normalised biases on the three components are $\mu_{r}=-0.02$, $\mu_{\theta}=-0.23$ and $\mu_{\phi}=0.0$). The histograms of prediction errors are reasonably close to Gaussian for the three components with observation errors that appear to be under-estimated on average, in particular on $B_{r}$ (normalised r.m.s. errors on the three components are $\sigma_{r}=2.18$, $\sigma_{\theta}=1.55$ and $\sigma_{\phi}=1.63$). The SV predictions errors are remarkably consistent with the a priori errors with small biases and standard deviation close to unity for the three components ($\mu_{r}=-0.06$, $\sigma_{r}=1.01$; $\mu_{\theta}=-0.09$, $\sigma_{\theta}=1.11$ and $\mu_{\phi}=0.03$, $\sigma_{\phi}=1.14$), even though the distributions appears more peaked than a Gaussian. The Kalman filter employed here implicitly assumes Gaussian distributed data errors. However, the above remark suggests that alternative treatments of data residuals may be worth considering in future studies (e.g. L1 or Huber norms, see Constable(1988); Farquharson & Oldenburg(1998)).

Figure 4: Top: Histograms of MF prediction errors $\delta_{MF}$ (eq. 23), accumulated over all analysis epochs, normalised to the observation errors, for the components $B_{r}$ (left), $B_{\theta}$ (middle) and $B_{\phi}$ (right). Superimposed in black are the Gaussian distribution fits obtained with the mean $\mu$ and the variance $\sigma^{2}$ for each of the three distributions. Bottom: same histograms for the SV prediction errors $\delta_{SV}$.

3.1.2 Field models, and contributions to the SV

We now describe in more detail our MF and SV models. We present in Figure 5 MF and SV maps for our ensemble average model at the CMB truncated at spherical harmonic degree $n=14$. Comparing it to a more traditional field model CHAOS-6, which is temporally regularised, the overall agreement is very good, indicating that our tool is indeed capable of producing reasonable field models. MF discrepancies to CHAOS-6 are relatively small, with peak to peak values less than $10\%$ of the total amplitude for a field truncated at degree 14. They are dominated by isotropically distributed, small length-scale patterns. As well as being dominated by small length-scales, the disagreements are larger for the SV, with peak to peak differences about $30\%$ of the total amplitude, which is to be expected given the blue SV spectrum at the CMB, meaning that small length scales dominate. Interestingly, the largest differences are localised under South America and the Indian Ocean, where the planetary gyre respectively detaches from and joins the equatorial belt (Pais & Jault(2008)) and where rapid time-dependence is observed (Finlay et al.(2016a)Finlay, Aubert, & Gillet).

Figure 5: Top: CMB maps of the ensemble average radial magnetic field $\left<{\bf b}\right>$ (eq. 21) in 2017 (left: MF in mT; right: SV in $\mu$T/yr), as estimated with our algorithm. Bottom: MF (left) and SV (right) maps of the difference of our ensemble average field with CHAOS-6 (truncated at degree 14) at the CMB (with the same colorscales).

In Figure 6 we show the various contributions in our model to two SV spherical harmonic coefficient series. The dispersion within the ensemble of models is large enough to include time changes as estimated by CHAOS-6, with some exceptions during the high solar activity era, e.g. in 2002 for $h_{6}^{6}$, and at the very end of the CHAOS-6 era (this latter possibly in link with the damping of SA towards end-points in the regularised field model). We notice a larger spread of the analysis for the axial dipole than for non-zonal coefficients of intermediate length-scale such as $h_{6}^{6}$. This may be a consequence of the weaker constraint on zonal coefficients from surface observations (e.g. Kotsiaros & Olsen(2012)), although we only notice such behaviour for $g_{1}^{0}$. An enhancement of the dispersion is noticeable between 2010 and 2014, displaying in the spectral domain the impact of the decreasing number of data during this era when no vector satellite data were available. Over 2001–2006, the ensemble average $h_{6}^{6}$ trajectory shows distinctive square shaped variations, probably partly related to variations in the number of data satisfying selection criteria during this interval of enhanced solar activity when only CHAMP data were available.

Figure 6: SV spherical harmonic coefficient time series for $\dot{g}_{1}^{0}$ (left) and $\dot{h}_{6}^{6}$ (right). Predictions from our ensemble average model are shown in dark red ($\pm 2\sigma_{\dot{b}}$ in red), and CHAOS-6 in black. Contributions from subgrid errors and diffusion extracted from our ensemble of realizations are superimposed in respectively blue and yellow (with dispersions $\pm 1\sigma_{diff}$ and $\pm 1\sigma_{er}$ in the corresponding colors).

Spatial spectra shown in Figure 7 summarise the characteristics of our model in the spectral domain. We find excellent agreement with CHAOS-6 for the main field and its secular variation, except at the small length scales of the SV ($n>10$), which are more likely to be affected by the different data set chosen and by the different temporal kernel used (short time windows in our case against whole time-span for CHAOS-6). The ensemble spread gives a good approximation of the characteristic distance between our model and CHAOS-6. Diffusion and subgrid errors in the SV have approximately the same amplitude except for the dipole. The power stored in these two SV sources represents about $10$ to $20\%$ of the total SV energy at all scales.

Even though the dispersion within the model predictions is large enough to encompass most of the MF and SV observations, the dispersion within the ensemble of realizations is lower, by a factor about 2.5, than the distance between the ensemble average model and CHAOS-6 for both the MF (at all length-scales) and the SV (towards small length-scales only). A complete account of SV errors from all subgrid interactions (see Baerenzung et al.(2017)Baerenzung, Holschneider, Wicht, Sanchez, & Lesur) may help reduce the above under-estimation. Our current estimate is nevertheless larger than that obtained for the COV-OBS.x1 model Gillet et al.(2015a)Gillet, Barrois, & Finlay. We suspect that the accumulation of data (assumed independent) during the construction of this latter field model involved too strong a decrease of the posterior error within the COV-OBS framework. The more consistent approach to error propagation developed here and presented in Figure 7 favours larger uncertainties on spherical harmonic coefficients during the satellite era.

Figure 7: Top: time averaged spatial power spectra at the Earth surface of the magnetic field of CHAOS-6 ($\hat{\cal R}_{b^{c}}^{a}$, eq. (24), in green), our estimate ($\hat{\cal R}_{\left<b\right>}^{a}$, in red), the difference between the two ($\hat{\cal R}_{\delta b}^{a}$, red thin line) and the dispersion within our ensemble of analyses ($\hat{\cal R}_{\sigma_{b}}^{a}$, dotted line). Bottom: idem for the SV, superimposed with the spectra of the contributions from subgrid errors (blue) and from diffusion (yellow).

Overall, we are generally able to retrieve earlier well-established results. For instance the contribution from advection dominates (over diffusion) the axial dipole decay (Finlay et al.(2016b)Finlay, Olsen, Kotsiaros, Gillet, & Tøffner-Clausen; Barrois et al.(2017)Barrois, Gillet, & Aubert) and its fluctuations – even though our estimate for the contribution from diffusion to $dg_{1}^{0}/dt$, shifted upward by a couple of nT/yr in comparison with the results of BGA17 (see §2.2.2), amounts to a relatively larger fraction over the latest years where the dipole decay tends to be weaker. The ensemble average SV originating from diffusion is presented in Figure 8 for 2017: the most significant contributions appear below Africa and Indonesia. The strongest diffusion appears linked to intense patches of up/down-wellings in the equatorial belt at the CMB (see Figure 8) and/or where strong gradient of ${\bf B}$ occur. This is a direct consequence of our estimation of diffusion through cross-covariances involving core surface velocity and magnetic fields (see BGA17 and Amit & Christensen(2008)). In the framework of our modelling, such diffusion patterns seem to be required by magnetic observations rather by the imposed prior cross-covariances (or if it is the case, it does not show up in the background state).

Figure 8: Magnetic diffusion at the CMB (top, color-scale in $\mu$T.y${}^{-1}$), and horizontal divergence $\nabla_{h}\cdot{\bf u}_{h}$ (bottom, color-scale in $10^{-3}$ yrs${}^{-1}$) superimposed with passive tracers trajectories (black, tracer size scale in km/yr), for the ensemble average model in 2017. Core flow visualisations are performed using the tools provided at https://geodyn.univ-grenoble-alpes.fr/. The size of the tracers is proportional to the velocity field (see the legend). The initial positions of the tracers is random; each trajectory is advected by the velocity field for a fixed time; along each trajectory, the late (early) positions are darker (lighter).

3.2 Core flows solutions

Next, we study with more details the temporal information contained in our core flow solutions. The idea is to extract an average signal and a linear acceleration, together with the flow at different periods, to check if we witness any preferential frequency, or if the characteristics of the flow change with the period. To do so, we apply a least-squares regression to our core flow solution with a function of the form

$${\bf u}(t)=\hat{\bf A}+{\bf A}_{L}(t-t_{0})+\displaystyle\sum_{k=1}^{11}\left[{\bf A}_{k}^{s}\sin\left(2\pi(t-t_{0})\dfrac{k}{T}\right)\right.$$ (27)

$$\left.+{\bf A}_{k}^{c}\cos\left(2\pi(t-t_{0})\dfrac{k}{T}\right)\right]\,,$$

with $t_{0}=(t_{i}+t_{f})/2=2008.92$ and $T=t_{f}-t_{i}=16$ yrs. Vectors $\hat{\bf A}$, ${\bf A}_{L}$, ${\bf A}_{k}^{c}$ and ${\bf A}_{k}^{s}$ store respectively the spherical harmonic coefficients of the time average velocity, time average flow acceleration, and cosines and sines from periods $16$ yrs (for $k=1)$ to $1.45$ yrs (for $k=11$) – of course the longer periods are not well constrained given the short time-span considered here.

We show in Figure 9 the norm (26) of all flow constituents for the ensemble average solution. The flow is dominated by long periods, translating onto core surface motions the red SV temporal spectrum (see Gillet et al.(2015a)Gillet, Barrois, & Finlay; Lesur et al.(2017)Lesur, Wardinski, Baerenzung, & Holschneider). In comparison with a r.m.s. time average flow of 11.1 km/yr, the linear acceleration ${\bf A}_{L}$ corresponds, integrated over 16 yrs, to a r.m.s. flow increment of 6.6 km/yr.

Figure 9: Core flow norm ${\cal N}$ for all flow constituents that enter equation (27). The norm $\mathscr{N}$ for the linear flow acceleration is obtained by integrating the linear trend over the 16 yrs.

3.2.1 Stationary motions, and flow model uncertainties

We show in Figure 10 core surface maps of the flow intensity and tracers trajectories for the ensemble average flow constituents $\hat{\bf A}$. We retrieve on the map for the time average flow classical features, such as the westward gyre offset towards the Atlantic Ocean found in many studies (e.g. Pais & Jault(2008); Gillet et al.(2015b)Gillet, Jault, & Finlay; Aubert(2014); Baerenzung et al.(2017)Baerenzung, Holschneider, Wicht, Sanchez, & Lesur), with a Pacific hemisphere that is on average much less energetic. The most energetic flow features are associated with (i) azimuthal motions in the equatorial belt below Africa, (ii) high latitudes azimuthal jets in the Pacific hemisphere and (iii) meridional circulations, poleward (resp. equatorward) around 90${}^{\circ}$W (resp. 90${}^{\circ}$E).

Figure 10: Intensity maps at the CMB of the flow constituents $\hat{\bf A}$ (in km/yr) for the ensemble average flow solution, superimposed with passive tracers trajectories (black). Top: Aitoff projection. Middle: North (right) and South (left) polar projections. Bottom: Aitoff projection for equatorially symmetric (left) and antisymmetric (right) components. The colorscale and tracer size scale are the same for all sub-figures.

Our solution is dominated by equatorially symmetric features (see Figure 10, bottom), as expected outside the tangent cylinder (or TC, the cylinder tangent to the inner core, whose axis coincides with to the rotation axis) when rotation forces dominates the momentum balance (e.g. Pais & Jault(2008)). Nevertheless, the symmetry may be locally broken. The most striking examples of this are anti-cyclonic circulations within the TC, retrieved in both the Northern and Southern hemispheres (Figure 10, middle). In contrast with polar vortices previously inferred from geomagnetic observations (Olson & Aurnou(1999); Amit & Olson(2006)), features we isolate here are off-set to one side of the polar caps (i.e. they contain an important $m=1$ contribution). This is a common configuration for polar vortices found in the most up to date numerical simulations (Schaeffer et al.(2017)Schaeffer, Jault, Nataf, & Fournier), which show much variability through epochs.

We show in Figure 11 the time-average spatial power spectra for the ensemble average solution and for the dispersion within the ensemble of models. The former is comparable with the spectrum of the prior CED. The latter indicates that uncertainties, as measured by the ensemble spread, constitute a large fraction of the flow magnitude for degrees $n\geq 10$. The oscillation in the power seen between odd and even degrees might be magnified by possibly too low subgrid error budget (see §3.1.2).

Figure 11: Time averaged spatial power spectra for ensemble average core flows ($\hat{\cal S}_{\left<u\right>}^{a}$, red thick line) and the spectra for the ensemble average of each realization ($\left<\hat{\cal S}_{uk}^{a}\right>$, yellow thick dotted line), eq. 25): obtained from the re-analysis of VO and GO data. Spectra for the corresponding dispersion within the ensembles of models are displayed in dotted lines. In green is shown the averaged spectrum for the prior CED.

3.2.2 On transient core surface motions

We now explore transient flow motions. We particularly focus on the amount of equatorial symmetry of our solutions inside and outside the TC, in order to detect if our model is sensitive to the specific geometry of the Earth’s core (does it hold a signature of the TC?). As for the time-average flow, the linear acceleration over the past 16 yrs is primarily symmetric with respect to the equator (see Figure 12). The largest contributions consist of accelerating circulations around the meridional, Eastern branch of the gyre. Associated with these time-changing eddies around the equatorward branch of the planetary gyre, an Eastward equatorial jet intensifies under the Western Pacific. This suggests an underlying dynamics more complex than a simple longitudinal shift of the planetary gyre.

Figure 12: Same as Figure 10 for the flow constituent ${\bf A}_{L}$ (in km/yr${}^{2}$).

Interestingly, our average solution does not show a major intensification of equatorially symmetric azimuthal jets at high latitudes in the Pacific hemisphere, as inferred by Livermore et al.(2017)Livermore, Hollerbach, & Finlay. Indeed, we see an increase of the Northern jet only, by about $67\%$ in average (the one $\sigma$ dispersion within the ensemble of flow realizations allowing for an increase up to $100\%$). Although still an appreciable acceleration, it is significantly less than the factor of 3 found by Livermore et al.(2017)Livermore, Hollerbach, & Finlay. The disagreement is likely due to our global inversion (in opposition to their local model). The difference seems to be related with anti-symmetric circulations within the TC. One should keep in mind that in these high and low latitude areas, gradients of $B_{r}$ are much larger in the Northern hemisphere, meaning that the signature of any motions near the TC below the Southern Pacific are significantly weaker. As for the stationary constituent, the equatorial symmetry is not perfectly respected, and we retrieve the largest anti-symmetrical features within the TC, associated with polar jets.

We give in Figure 13 an example of one interannual flow constituent at the CMB for a period of 5.3 yrs. In this case, the most energetic flows are concentrated into non-axisymmetric azimuthal jets near the equator (already highlighted by Gillet et al.(2015b)Gillet, Jault, & Finlay; Finlay et al.(2016b)Finlay, Olsen, Kotsiaros, Gillet, & Tøffner-Clausen), and into localised circulations at mid and high latitudes. These are not confined to the Atlantic hemisphere: despite being less energetic on average, the Pacific hemisphere shows interesting interannual flow variations. At these sub-decadal periods, we have not detected any obvious propagation of non-zonal flow patterns, which might be interpreted as the signature of azimuthally propagating waves (as advocated for by Chulliat & Maus(2014); Chulliat et al.(2015)Chulliat, Alken, & Maus). The other periods display globally the same kind of features and no particular behaviour is found at any period. At these time-scales also show up less intense anti-symmetric features; the most significant shows up in the equatorial area (for instance under the Atlantic ocean and the Western Pacific), and towards high latitudes on the edge of the TC.

Figure 13: Same as Figure 10 for the flow constituent ${\bf A}_{3}^{c}$ (in km/yr).

Figure 14 summarises the amount of equatorial symmetry found in regions inside and outside the TC, for our core flow solutions at all periods. It appears almost independent of the considered period: outside the TC, it is within $90$ to $95\%$ of the surface energy for all flow constituents of equation (27). The partition of energy between symmetric and anti-symmetric flow components is more balanced inside the TC where, depending on the considered time-scale, $\approx 55\pm 15\%$ of the energy is contained in equatorially symmetric flows. This latter observation could be expected because the presence of the inner core is intended to partially break the equatorial symmetry However, it is remarkable that the algorithm appears accurate enough to detect a specific behaviour within the tiny areas covered by polar caps. Moreover, although our ensemble average model and the CED show very similar amounts of equatorial symmetry outside the TC (the value for the CED model is $95\%$ of symmetrical flows inside and outside TC), they differ significantly inside the TC (it is much less in the inverted flows). As a consequence, the larger proportion of equatorial antisymmetry inside the TC is driven by observations (against the prior information).

Figure 14: Fraction of energy contained into the equatorial symmetric part of the flow, inside (blue line) and outside (yellow line) of the tangent cylinder (TC), for each of the flow constituent that enter equation (27). The total symmetric part of the flow is also displayed in green. The value for the CED dynamo used as a prior, is $0.95$ both inside and outside the TC.

4 Summary and discussion

Following earlier strategies for geomagnetic field model reconstruction (e.g. Jackson et al.(2007)Jackson, Constable, Walker, & Parker; Lesur et al.(2010)Lesur, Wardinski, Asari, Minchev, & Mandea), and moving towards geomagnetic data assimilation (Aubert(2015); Gillet et al.(2015a)Gillet, Barrois, & Finlay; Baerenzung et al.(2017)Baerenzung, Holschneider, Wicht, Sanchez, & Lesur), we continue the work initiated in BGA17. We retain their idea of combining spatial information from numerical simulations of the geodynamo with temporal information implemented through stochastic equations, chosen to replicate the frequency spectrum of ground-based geomagnetic series. However, instead of considering spherical harmonic coefficients of the main field as data, here we have inverted observations (GOs and VOs) directly, at and above the Earth’s surface. In this respect we follow the studies by Beggan & Whaler(2009) and Whaler & Beggan(2015), although we account for subgrid processes (of great importance, as shown by BGA17 or Baerenzung et al.(2016)Baerenzung, Holschneider, & Lesur) and for surface magnetic diffusion. This avenue allows us to propose PDFs for the main field and its secular variation, as well as for the recovered core motions.

4.1 Geophysical insights

The MF models presented here are consistent both with observations and with the imposed dynamical prior. The model uncertainties, as suggested by the ensemble spread, are slightly less than the distance of the average model to CHAOS-6. We recover in our core flow solutions a westward gyre that circulates around the TC at high latitudes in the Pacific hemisphere, and flows closer to the equator in the Atlantic hemisphere. The largest contributions from magnetic diffusion are associated with up/down-wellings where the gyre meets the equatorial region (under Indonesia) and in the equatorial region below Africa. At all time-scales, the flow is predominantly symmetric with respect to the equator, except inside the TC where the situation is more balanced (contrary to our dynamo prior that is mostly symmetric everywhere).

The most intense time-average flow acceleration over the past 16 years is linked with evolving meanders around the equatorward branch of the gyre in the Eastern hemisphere, also associated with the appearance of an Eastward equatorial jet under the Western Pacific. We do find a decadal intensification of jets near the TC, although the magnitude of the acceleration we infer is lower than that estimated by Livermore et al.(2017)Livermore, Hollerbach, & Finlay with their reduced model. In our study, it is furthermore confined to the Northern hemisphere. This equatorial asymmetry may be interpreted as the signature of an ageostrophic acceleration, keeping in mind that main field gradients are weak in the Southern Pacific, implying a weaker constraint on flow motions there (see Figure 7 in Baerenzung et al.(2016)Baerenzung, Holschneider, & Lesur). However, because our prior does not show any particular bias in those areas, it is likely that those features are mostly driven by the data. On interannual periods, we find relatively energetic flow changes in both the Atlantic and the Pacific hemispheres, with both non-zonal equatorial jets and time-dependent mid-to-high latitudes eddies evident.

4.2 Future work

We currently lack a physical understanding for the features described above, whether it be through quasi-geostrophic flows (e.g. Labbé et al.(2015)Labbé, Jault, & Gillet), motions within a stratified layer (e.g. Buffett & Knezek(2017)), or any other interpretation through a reduced model. We also lack suitable long coverage by high quality satellite records to perform spectral analyses with a refined sampling in the frequency domain, which would allow us to isolate possible waves at interannual periods. Development of such reduced models, and their coupling with stochastic processes for modelling unresolved processes, will be an important next step in our ability to understand and predict geomagnetic field changes.

Meanwhile, our stochastic model itself could be improved; in particular it is desirable to avoid driving back the average trajectory towards an average dynamo simulation. This is indeed an unlikely state for the current era (say over decadal to centennial time-scales), which might be better represented by a re-analysis of for instance centennial motions from historical records (Jonkers et al.(2003)Jonkers, Jackson, & Murray). Furthermore, because of the short time-span covered today by satellite data, we found it challenging to derive well-conditioned matrices for VO uncertainties. This is a key-point for such data assimilation studies, which calls for further developments, e.g. through projections onto reduced basis in the data space. Alternatively, we may wish to co-estimate, together with the core state, time-dependent external fields. Although possible, this calls for a severe re-encoding of both the forecast and analysis steps, in order to integrate satellite measurements along the tracks.

The general philosophy of our work is to retrieve information on the state of the Earth’s core, and to provide realistic uncertainties on all state variables in a simple way. The encouraging magnetic models obtained with this approach render our algorithm suitable for deriving candidates to the International Geomagnetic Reference Field (Thébault et al.(2015)Thébault, Finlay, Beggan, Alken, Aubert, Barrois, Bertrand, Bondar, Boness, Brocco, et al.). Remaining in a stochastic framework, modifications of the forward model parametrisation – such as accounting for a background state closer to the flow responsible for the magnetic field over the past decades – may extend the prediction capability of our algorithm. However, targeting accurate field predictions one will have to resort to deterministic (i.e. dynamically based) equations for the core state.

5 Acknowledgements

We thank Julien Aubert for providing the Coupled-Earth dynamo series used to build the core state statistics, and Loic Huder for spotting two errors in the code at the origin of the results of BGA17. We also thank Julien Baerenzung and an anonymous referee for their useful comments that helped improve the quality of our manuscript. We would like to thank as well GFZ German Research Centre for Geoscience for providing access to the CHAMP MAG-L3 data and to ESA for providing access to the Swarm L1b MAG-L data. We also like to thank the staff of the geomagnetic observatories and INTERMAGNET for supplying high-quality observatory data. NG and OB were partially supported by the French Centre National d’Etudes Spatiales (CNES) for the study of Earth’s core dynamics in the context of the Swarm mission of ESA. ISTerre is part of Labex OSUG@2020 (ANR10 LABX56), which with the CNES also finance the phd grant of OB. Numerical computations were performed at the Froggy platform of the CIMENT infrastructure (https://ciment.ujf-grenoble.fr) supported by the Rhône-Alpes region (GRANT CPER07 13 CIRA), the OSUG@2020 Labex (reference ANR10 LABX56) and the Equip@Meso project (referenceANR-10-EQPX-29-01). MH and CF were supported by the Danish Council for Independent Research - Natural Sciences, Grant DFF-4002-00366.