# Waves in Earth's core and geomagnetic field forecast

N. Gillet<sup>a,*</sup>, F. Dall'Asta<sup>a</sup>, P.-O. Amblard<sup>b</sup>, R. Claveau<sup>a</sup>, J. Aubert<sup>c</sup>

## ARTICLE INFO

### Keywords:

Earth's core

Geomagnetic secular variation

Inverse problem

Hydromagnetic waves

## ABSTRACT

We use advanced numerical geodynamo series to derive a reduced stochastic model of the dynamics at the surface of Earth's core. Considering order 3 autoregressive (AR-3) processes allows to replicate the simulated spatiotemporal spectrum over a broad range of time-scales, spanning millennia to a fraction of year, including the cut-off found for periods shorter than approximately 2 years and associated with magnetic dissipation. We show how to derive such a forward model from a variety of input simulation series, and present its implementation into the pygeodyn data assimilation algorithm, based on a sequential ensemble method. The updated scheme is applied to perform magnetic field hindcasts and core flow reanalyses. For all observable length-scales, the rate of change of the observed magnetic field is most of the time accounted for within the spread of the forward model trajectories. AR-3 predictions on average supersede by about 35 % linear extrapolations on short (2 yr) time-scales, reducing high-frequency spurious variations in reanalysed flow motions. This improvement is reduced to ≈ 10% for 5 yr increments, with a large variability from one epoch to the other depending on the overall curvature of the magnetic field evolution. We perform a reanalysis over the period 1880–2023 covered by observatory and satellite records. We find enhanced kinetic energy in three period ranges around 12.5, 6.5 and 3.5 years. At all three periods, fluid motions share geometrical properties compatible with quasi-geostrophic magneto-Coriolis waves: equatorial symmetry, larger amplitude near the equator, flow dominated by low azimuthal wave number and modulated in longitude, phase speed much faster than the fluid velocity and decreasing with the period. At 6.5 yr period we trace back to the mid-1990’s the patterns previously detected from satellite data. We also find in the 1960–70’s a similar wave-train, possibly in link with the 1969 geomagnetic jerk. The AR-3 model, in conjunction with early satellite records, likely helps isolate such coherent features on interannual time-scales. Similar wave-like motions also show up at 3.5 yr period around 1970 and during the past decades. At periods around 12.5 yr we detect recurrent patterns starting as far back as 1920, and modulated over decadal time-scales. Our results show growing evidence for core dynamics governed by the presence of hydro-magnetic waves over a wide range of periods. This may allow deterministic and/or empirical descriptions of the signal that may help sound deep Earth’s properties, and improve predictions of the magnetic field evolution.

## 1. Introduction

One of the major challenges when modeling the geodynamo is to embrace the wide variety of spatial and temporal scales covered by the Earth's core magneto-hydro-dynamics (Nataf and Schaeffer, 2015). Recent advances in computation facilities have enlarged the range of time-scales over which the outputs from numerical experiments can be confronted with geomagnetic observations. With early computations, static properties could be challenged (Christensen et al., 2010).

By approaching the input parameters of direct numerical simulations (DNS) to the geophysical configuration, force balances as well as transient motions are more and more accurately reproduced (Schaeffer et al., 2017; Schwaiger et al., 2019). However, with DNS one can hardly model the physics relevant for the observed changes at decadal and shorter time-scales. A parameterization of the small length-scales interactions must be adopted such as large-eddy simulations (LES, see Aubert et al., 2017) in order to reach a dynamical regime directly comparable with observatory and satellite geomagnetic records. This involves at first

<footnote><sup>a</sup> Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, Univ. Gustave Eiffel, ISTERe, Grenoble 38000, France</footnote>

<footnote><sup>b</sup> CNRS, Univ. Grenoble Alpes, Grenoble INP, GIPSA-Lab, Grenoble 38000, France</footnote>

<footnote><sup>c</sup> Univ. Paris Cité, Institut de Physique du Globe de Paris, CNRS, Paris 75005, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: nicolas.gillet@univ-grenoble-alpes.fr (N. Gillet).</footnote>
order a dominant Coriolis force, equilibrated by the pressure gradients (quasi-geostrophy, or QG, rendering fluid motions in the core almost invariant along the Earth's rotation axis), and at second order a balance between the magnetic, Archimedes and Coriolis forces (the QG-MAC balance, see Aubert, 2019).

Using LES, Earth-like conditions have recently been reached with geodynamo simulations (Aubert, 2023), within the “path” set-up designed by Aubert et al. (2017). However, it still remains computationally demanding to integrate the primitive equations over several decades to centuries (the time covered by historical and modern magnetic records, e.g. Jackson and Finlay, 2015; Hulot et al., 2015). Meanwhile, the confrontation of forward Earth’s core models to geomagnetic observations via data assimilation algorithms (Fournier et al., 2010; Sanchez et al., 2019) requires the integration of a large number of model trajectories, in order to obtain a proper measure of the uncertainties on the core state variables. In this context, some first attempts at developing reduced QG models where magnetic forces play a key role have arisen (Canet et al., 2009; Jackson and Maffei, 2020), but these are not yet operational.

Barrois et al. (2017) then Gillet et al. (2019) followed an alternative avenue of research, considering empirical reduced models of the core surface dynamics anchored to advanced numerical simulations of the geodynamo (see also Ropp and Lesur, 2023). Based on order 1 autoregressive (or “AR-1”) stochastic processes to model the core flow, they not only replicate the spatiotemporal second-order statistics for the simulated dynamo flow on decadal and longer periods. They also extend towards short periods the core surface flow power spectrum $S_U(f)$ as $f^{-2}$ (where $f$ is the frequency), in a manner that is consistent with the frequency spectrum $S_B(f)\propto f^{-4}$ observed for the magnetic field from ground-based and satellite records on periods longer than a few years (Gillet et al., 2015a; Lesur et al., 2018; Baerenzung et al., 2018). This latter property is welcome as geodynamo simulations away from Earth-like conditions present a too low magnetic to kinetic energy ratio. As a consequence they do not separate enough the magnetic time (or Alfvén time, $\tau_A \approx 2$ yr in the core) from the turn-over time ($\tau_U \approx 130$ yr). A too large value for $\tau_A$ in simulations tends to shrink the frequency spectrum at high frequencies.

Stochastic forward models have been implemented into the pygeodyn data assimilation software (Huder et al., 2019), where they are plugged into an augmented state ensemble Kalman filter (EnKF, see Evensen, 2003): some hundreds of model state trajectories are calculated in parallel, in order to approximate the model state second-order statistics (a quantity required by the EnKF algorithm). The stochastic approach has proven useful for instance in isolating wave patterns from the satellite magnetic observations (Gillet et al., 2022b; Ropp and Lesur, 2023). Interpreted as QG Magneto-Coriolis (or QG-MC) waves by Gillet et al., they are responsible for geomagnetic jerks (i.e. abrupt changes in the rate of change of the magnetic field), or equivalently pulses in the norm of the second time-derivative of the field (e.g. Chulliat and Maus, 2014; Kloss and Finlay, 2019). The interpretation in term of QG-MC waves is reinforced by the analysis of advanced numerical simulations (Aubert et al., 2022; Schwaiger et al., 2024). Nevertheless, one cannot rule out an alternative explanation to transient SV changes based on Magneto-Archimedes-Coriolis (or MAC) waves, in the presence of a stable layer below the core surface. Although these have been invoked to interpret patterns at both high and low latitudes (Chi-Durán and Buffett, 2023; Buffett and Matsui, 2019), the morphology of the observed magnetic field puts an upper bound to the degree of stratification, which shall be mild enough to allow for deep convection columns to penetrate the top layer (Gastine et al., 2020).

Reduced stochastic models come at a limited numerical cost, for two main reasons: first they only require a model state of limited dimensionality to mimic the observable dynamics (some hundreds of parameters), and second they allow to consider a relatively large time-step (typically a fraction of year) for the forward integration of the dynamical

model. This renders possible the implementation of dedicated algorithms, capable of separating as much as possible the information from the noise within the ensemble statistics (namely the Graphical LASSO, see Friedman et al., 2007), and avoids drawbacks often met with assimilation methods such as ensemble collapse (for an illustration of these with pygeodyn, see Istas et al., 2023).

However, the characteristic spectral behavior of AR-1 models is not appropriate to reproduce dissipative effects observed on periods shorter than $\tau_A$ in geodynamo simulations, in particular a cut-off in the spectra $S_U(f)$ and $S_B(f)$ for the flow and the magnetic field, in link with magnetic diffusion (Aubert and Gillet, 2021). As a consequence, AR-1 models allow a too large variance on time-scales less than $\tau_A$, and tend to generate noisy inverted core flow series. In order to address this issue, we propose in the present study to fit higher-order AR processes to dynamo simulation series. This is made possible thank to a few available geodynamo runs at Earth-like conditions (the 100p dynamo, see Aubert, 2023). Since these cover relatively short durations, the information they carry on short time-scales must be merged with that from a lower resolution simulation computed over a long time-span (the 71p dynamo by Aubert and Gillet, 2021). The order 3 (“AR-3”) model thereafter produced is integrated into the pygeodyn tool, and applied to the reanalysis of Gauss coefficient magnetic field data over the observatory and satellite eras.

The present study is organized as follows. We present in Section 2 the field models considered as magnetic data, and the dynamo series used to construct our new stochastic model. In Section 3 we present the derivation of the AR-3 forward model, and the modifications it implies for the pygeodyn assimilation tool, including the analysis step of the Kalman filter. We illustrate in Section 4 this revised algorithm with two applications: the prediction of the magnetic field, and the recovery of wave-like motions within Earth’s core. In Section 5 we discuss implications from our developments, and future ways forward.

## 2. Data

### 2.1. Geomagnetic data

In this study we consider as geomagnetic data Gauss coefficients from two field models, COV-OBS-x2 (Huder et al., 2020) and CHAOS-7 (Finlay et al., 2020), that we consider over the periods 1880–2020 and 1997–2023, respectively. Main field (MF) and secular variation (SV) Schmidt semi-normalized Gauss coefficients describing the Earth’s core magnetic field evolution are integrated into the pygeodyn tool (both considered up to spherical-harmonic degree $n = 13$), with analyses performed every $\Delta t^d = 1$ yr (see details in Section 3.3).

The construction of COV-OBS-x2 relies on a priori information coherent with the occurrence of abrupt changes in the rate of change of the field (geomagnetic jerks, see Mandea et al., 2010). This is operated by considering a priori cross-covariances from autoregressive processes of order 2 (Gillet et al., 2013), in agreement with a temporal power spectrum of MF series behaving as $S_B(f)\propto f^{-4}$ in the range $2\lesssim 1/f \lesssim 100$ years (see Lesur et al., 2018). This permits the derivation of a posteriori model error variances $\sigma_{MF,m}^2$ and $\sigma_{SV,m}^2$ on the MF and SV Gauss coefficients. However, by construction these do not account for the rapid field evolution (time-scales $\lesssim 3$ years), due to the low-pass filtering induced by the projection onto splines functions with knots every 2 yr (Pick et al., 2019). We thus add in quadrature errors of representativeness (of variances $\sigma_{MF,r}^2$ and $\sigma_{SV,r}^2$), which correspond to the unmodelled rapid changes, as prescribed by Huder et al. (2020). The total error budget for MF and SV COV-OVS-x2 data is then $\sigma_{MF}^2 = \sigma_{MF,m}^2 + \sigma_{MF,r}^2$ and $\sigma_{SV}^2 = \sigma_{SV,m}^2 + \sigma_{SV,r}^2$.

The time evolution of the CHAOS-7 core field is regularised by penalizing the third time-derivative of internal Gauss coefficients up to degree 20, plus adding a penalty on the second time-derivative at both end-points. The penalization varies with the spherical harmonic degree
and order, reducing the impact on core field changes of unmodelled electrical currents towards auroral regions. This construction does not allow the derivation of model uncertainties on the parameters that describe the MF and its SV. We thus associate to CHAOS-7 Gauss coefficients MF and SV uncertainties as estimated with the Kalmag model (Baerenzung et al., 2022). We consider these as a reasonable error estimate, as their power is comparable to that of the difference between CHAOS and Kalmag models for all considered length-scales (see the Lowes spectra in Fig. 3 of Baerenzung et al., 2020). Kalmag being based on a sequential Kalman filter approach, these are considered as instantaneous error estimates, with no need to add an error budget for unmodelled processes.

## 2.2. Geodynamo simulation data

The stochastic model plugged into the assimilation tool pygeodyn is anchored to spatiotemporal statistics learnt from various numerical simulations of the geodynamo. We base this knowledge on several members of the “path” dynamos (Aubert et al., 2017). They impose a mantle control associated with a heterogeneous mass distribution in the mantle, gravitationally coupled to a non-spherical inner core (Aubert et al., 2013), plus an electro-magnetic torque associated with a thin conducting layer at the base of the mantle. These are adjusted in order to generate an overall westward drift commensurate with that inferred from geomagnetic observations (Pichon et al., 2016; Schwaiger et al., 2024). The path dynamos are all dipole-dominated, and “Earth-like” in the sense of the static criteria put forward by Christensen et al. (2010).

The path dynamos follow a one-dimensional trajectory in the parameter space defined by the input dimensionless numbers, namely the magnetic Ekman number $E_m(\epsilon) = \eta/\Omega D^2 = \epsilon E_m^0$, the magnetic Prandtl number $P_m(\epsilon) = \nu/\eta = \sqrt{\epsilon} P_m^0$ (plus a Rayleigh number, a measure of the buoyancy force, and the thermal Prandtl number $Pr = \nu/\kappa$ set to 1), with $\nu$ the kinematic viscosity, $\eta$ the magnetic diffusivity, $\kappa$ the thermal diffusivity, $\Omega$ the Earth’s rotation rate, and $D$ the gap between the inner and outer core. The single variable $\epsilon$ defines how all the input parameters vary along the path, starting from “standard” simulations run at $E_m^0 = 1.2 \times 10^{-5}$ and $P_m^0 = 2.5$ (or $\epsilon = 1$). Aubert et al. (2017) consider $\epsilon = 10^{-7}$ (or $E_m = 2.37 \times 10^{-9}$, $P_m = 1.26 \times 10^{-3}$) as being representative of the Earth’s core conditions. Simulations along the path are named with respect to their relative position on a logarithmic scale between $\epsilon = 1$ (0 %-path or $0p$) and $\epsilon = 10^{-7}$ (100 %-path or $100p$). Such low values of the Ekman number are achievable thanks to a parameterization of the small-scale turbulence, based on a hyperdiffusion of momentum and heat at spherical harmonic degrees $n \ge 30$.

We consider in the present study several sets of simulation data, consisting of Gauss coefficients series for the magnetic field (up to spherical harmonic degree $n = 13$) and the flow (up to $n = 18$) at the surface of the core. First a 10 kyr-long series from the $71p$ dynamo (Aubert and Gillet, 2021), composed of $N^a + 1$ snapshots sampled every $\Delta t^* = 0.2$ yr. Second, 8 series ($j \in [1, 8]$) of the $100p$ dynamo (Aubert, 2023), each about 80 yr long and composed of $N^{b,j} + 1$ snapshots, also sampled every $\Delta t^*$. They are separated in two subsets (named $100p-A$ and $100p-B$), run with two different values of $\eta$, thus leading to slightly different values of the output parameters (see below). Considering these two subsets allows to consider a wider range of behaviors, given the limited duration of the $100p$ runs.

The main characteristics of the two dynamos are summarized in Table 1. Path dynamos are scaled according to the magnetic diffusion time $\tau_\eta = D^2/\eta$. This quantity is approximately preserved along the path, with $\tau_\eta \approx 135$ kyr for $71p$, 245 kyr for $100p-A$ and 185 kyr for $100p-B$. The different values correspond to various estimates of the magnetic diffusivity, all in line with geophysical estimates (one has $\tau_\eta \sim 200$ kyr for $\eta \sim 0.7$ m$^2$/s, see Davies et al., 2015). The turn-over time $\tau_U = D/U$ is also preserved along the path at values representative of Earth’s conditions (where $\tau_U \sim 150$ yr for a characteristic fluid velocity within Earth’s core $U \approx 15$ km/yr, see Finlay et al., 2010). Consequently the magnetic Reynolds number $R_m = \tau_\eta/\tau_U$ is also approximately preserved along the path, at a realistic value $R_m \approx 1100$ for $71p, \approx 1900$ for $100p-A$ and $\approx 1400$ for $100p-B$.

Advancing along the path, the magnetic energy $E_M$ becomes predominant in comparison with the kinetic energy $E_K$, as measured by the Alfvén number $Al = \tau_A/\tau_U \sim \sqrt{E_K/E_M}$. Here $\tau_A = D/U_A$ is the Alfvén time, with $U_A = B/\sqrt{\rho\mu}$ the Alfvén speed ($\rho$ the core density, $\mu$ the magnetic permeability of free-space and $B$ the magnetic field intensity within the core). At Earth’s core conditions, $Al \approx 10^{-2}$ is found for values of the magnetic field $B \approx 5$ mT (Gillet et al., 2010), corresponding to $\tau_A \sim 2$ yr. While the $100p$ dynamo runs lead to similar values $(\tau_A^{100p} \approx \tau_A)$, less advanced simulations over-estimate the Alfvén time, with $\tau_A^{71p} = 5.8 \approx 3\tau_A$ yr in $71p$. As a consequence, the flow variability on time-scales shorter than $O(\tau_A^{71p})$ will be much attenuated in comparison to the Earth’s core, rendering this dynamo non-geophysical at periods $T \lesssim \tau_A^{71p}$. This point is illustrated in Fig. 1, where we superimpose idealized temporal spectra for the flow over periods $0.1 < T < 10^4$ kyr, of the form

<figure><img src="image_1.png" /><figcaption>Fig. 1. Schematic view of the idealized flow temporal PSD $S_U(f)$ from Eq. (1), for the available path dynamos. The parameters considered are $\tau_U = 135$ yr, and $\tau_A$ as documented in Aubert (2018, 2023); Aubert and Gillet (2021). For comparison, the spectrum referenced to as “AR-1” corresponds to $S_U(f) \propto \frac{1}{1+(\tau_U f)^2}$, and is independent of $\tau_A$.</figcaption></figure>

**Table 1**

Main input and output dimensionless numbers and time-scales for the numerical simulations considered in this study. Times are given in years.

<table><thead><tr><th>dynamo</th><th>E<sub>m</sub></th><th>P<sub>m</sub></th><th>R<sub>m</sub></th><th>Al</th><th>τ<sub>A</sub></th><th>τ<sub>U</sub></th><th>τ<sub>η</sub></th></tr></thead><tbody><tr><td>71p</td><td>3.8 × 10<sup>-8</sup></td><td>7.9 × 10<sup>-3</sup></td><td>1100</td><td>4.6 × 10<sup>-2</sup></td><td>5.8</td><td>125</td><td>135,000</td></tr><tr><td>100p-A</td><td>2.4 × 10<sup>-9</sup></td><td>1.3 × 10<sup>-3</sup></td><td>1900</td><td>1.3 × 10<sup>-2</sup></td><td>1.7</td><td>130</td><td>247,000</td></tr><tr><td>100p-B</td><td>2.4 × 10<sup>-9</sup></td><td>1.3 × 10<sup>-3</sup></td><td>1400</td><td>1.2 × 10<sup>-2</sup></td><td>1.5</td><td>130</td><td>184,000</td></tr><tr><td>Earth</td><td>O(10<sup>-9</sup>)</td><td>O(10<sup>-6</sup>)</td><td>O(10<sup>3</sup>)</td><td>O(10<sup>-2</sup>)</td><td>≈ 2</td><td>≈ 130</td><td>≈ 200,000</td></tr></tbody></table>
$$S_U(f) \propto \frac{1}{1 + (\tau_{uf})^2 + Al^{-2}(\tau_A f)^6}, \quad (1)$$

and for several values of $\epsilon$ spanning $[10^0, 10^{-7}]$ – i.e. from 0 % to 100 %-path. Their trend changes from $S_U(f) \propto f^0$ on long periods, to $S_U(f) \propto f^{-6}$ on short time-scales, with $S_U(f) \propto f^{-2}$ in the intermediate range. The choice for such a spectral shape is motivated by the following. The cut-off between the two ranges on long and intermediate time-scales is governed by the turn-over time $\tau_U$, which is common to all dynamos along the path (due to the preservation of $\tau_U$). This is coherent with spectra of the secular acceleration $\partial_n^2 B_r$ that behave like $S_B(f) \propto f^2$ and are independent of the position along the path towards low frequencies (Aubert and Gillet, 2021). The second cut-off in (1) between the ranges at intermediate and short time-scales is controlled by the Alfvén time $\tau_A$, and moves towards high frequencies while advancing along the path and lowering the Alfvén number $Al$. This is coherent with the analysis of path dynamo data by Aubert and Gillet, which shows a cut-off at high frequencies involving magnetic diffusion and governed by the Lundquist number $S = \tau_\eta/\tau_A$. Since the diffusion time $\tau_\eta$ is comparable for all values of $\epsilon$, this cut-off thus evolves with $\tau_A$. Finally the intermediate range between $\tau_A$ and $\tau_U$ is coherent with the temporal spectrum of geomagnetic series (Gillet et al., 2015a; Lesur et al., 2018).

The above idealized spectra present a steeper decay on short time-scales than those discussed by Aubert and Gillet (2021). We further motivate this choice in Section 3.1, with our strategy to model the core surface dynamics using AR-3 stochastic equations. In their analysis of the 71p data, Aubert and Gillet indeed propose $S_B(f) \propto f^{-2}$ for periods down to about one year, which would correspond instead to $S_U(f) \propto f^{-4}$. However, the slope observed in the flow and field power spectra steepens towards short time-scales. Our choice with Eq. (1) of a -6 slope is motivated by our need to model the core dynamics time-scales as short as ≈ 0.1 yr. Considering even shorter periods would require to target even steeper spectral slopes.

Finally, we highlight that $R_m$ and $\tau_U$ are not the only quantities that are preserved along the path. It is also the case of the Elsasser number $\Lambda = \tau_\Omega \tau_U / \tau_A^2 = O(1)$, with $\tau_\Omega = 1/\Omega$ the time-scale associated with global rotation. It happens that the period of MC waves, for given wave-numbers, scales as $\tau_{MC} \propto \tau_A^2 / \tau_\Omega = \Lambda \tau_U$ (Gillet et al., 2022b). This time-scale is by deduction also preserved when varying $\epsilon$. Aubert and Gillet (2021) had shown that the convective dynamics on periods longer than $\tau_U$ is invariant along the path, and for periods between $\tau_U$ and $\tau_A$ the dynamics at play is in a QG-MAC balance (Aubert, 2019). Based on the above arguments, the dynamics on periods longer than the $O(\tau_A^{71p})$ is likely similar between the 100p and the 71p dynamos. We thus conclude that it is legitimate to consider simultaneously data from the 71p and 100p dynamos, provided that we account for the filtering of the short period dynamics in 71p. Covering 10 kyr, this simulation will thus constrain the reduced stochastic model from decadal to millennial time-scales. Conversely, it would be computationally exhausting to integrate over millennial time-scales a simulation at the parameters of 100p; this dynamo will thus essentially control decadal and shorter periods. We describe in the next section how we fit a linear autoregressive model of the core surface dynamics simultaneously to all these dynamo data. The data from the series are assumed independent the one from the other. Since the path dynamos are computed using a stress-free boundary condition at the core-mantle boundary (CMB), the flow data considered here are those from the top of the simulated core.

### 3. Method

The pygeodyn data assimilation tool (Huder et al., 2019) is based upon an EnKF (Evensen, 2003), made of a succession of forecast and analysis steps for an ensemble of realization of the model state. Forecasts consist in the forward integration of a dynamical system of equations for

the model state, here reduced to large-scale projections, at the core surface, of the magnetic field, the fluid flow, plus an error term accounting for unmodelled processes (the so-called “augmented state”). The forward model is embedding multivariate stochastic equations, based upon spatiotemporal statistics extracted from geodynamo simulations. The analysis step comes down to a regression of the model state, knowing observational data together with cross-covariances within the ensemble of realizations of the forecast state. We present in this section how we improve the assimilation tool pygeodyn from previous versions (Gillet et al., 2019; Istas et al., 2023). The main modification consists in extending the forward stochastic model to AR-3 processes, instead of AR-1 previously, in an attempt at better representing the rapid evolution of the magnetic field and flow (see Section 3.2). This is made possible thanks to the availability of new numerical simulations of the geodynamo run at closer to Earth-like parameters (Aubert, 2023). Such changes in the forward model impact the analysis step, described in Section 3.3. In particular, we now resort to a simultaneous analysis of data at three successive epochs, in order to access the first and second time-derivatives of the model state. It is conceptually equivalent to implementing a fixed-lag ensemble Kalman smoother (e.g. Evensen and Van Leeuwen, 2000; Cosme et al., 2012).

#### 3.1. The augmented model state

We use below the spherical coordinate system $(r, \theta, \phi)$. Considering an insulating mantle, the magnetic field vector $\mathbf{B}$ derives from a potential $V$ (i.e. $\mathbf{B} = -\nabla V$), and its evolution above the core surface is entirely prescribed by that of the radial field $B_r$ at the CMB, of radius $r = c$. Only the largest length-scales (spherical harmonic degrees $n \le L_B = 13$) are accessible, due to the steeper decay with radius of the core field power towards short length-scales, and to the presence of a dominant lithospheric field for $n \gtrsim L_B$ (Langel and Estes, 1982). We describe the time evolution of $B_r$ through the radial induction equation projected onto large observable length-scales (truncating the spherical harmonic expansion to degree $L_B$, as denoted by overlines):

$$\partial_t \bar{B}_r = - \overline{\nabla_h \cdot (\mathbf{u}_h \bar{B}_r)} + e_r. \quad (2)$$

Here $\mathbf{u}_h$ is the large-scale component of the horizontal velocity field at the core surface, while the errors of representativeness $e_r$ incorporate all contributions from unresolved processes, whether stemming from subgrid processes (that overall dominate this error budget) or from magnetic diffusion (for which one would need $\partial_r V_{r=c}$, a quantity that is not directly accessible from surface observations).

All fields are represented using Schmidt semi-normalized spherical harmonics of degree $n$ and order $m$. Gauss coefficients describing the MF $(g_n^m, h_n^m)$, its SV $(\partial_n g_n^m, \partial_n h_n^m)$, and the errors or representativeness are stored into vectors $\mathbf{b}$, $\partial_n \mathbf{b}$ and $\mathbf{e}$ respectively. All three vectors are of size $N_B = L_B(L_B + 2) = 195$. Spherical harmonic toroidal $(t_n^{m,c}, t_n^{m,s})$ and poloidal $(s_n^{m,c}, s_n^{m,s})$ coefficients describing the large scale part of the core surface flow are stored in a vector $\mathbf{u}$. For this field the expansion is restricted to degree $L_U = 18$, resulting in a flow vector of size $N_U = 2L_U(L_U + 2) = 720$. This choice for $L_U$ is a compromise chosen in order to limit aliasing effects on the inverted flow, while keeping a tractable model state dimension and mitigating the role of $e_r$. We refer to Gillet et al. (2019) for further details. $\mathbf{u}$ and $\mathbf{e}$ are described as perturbations with respect to a background state $(\mathbf{u}_0$ and $\mathbf{e}_0$, from the time-average of the 71p dynamo simulation).

The derivation of our empirical AR-3 model for the flow and errors of representativeness requires obtaining pairs of eigenvalues and eigenvectors from a finite number of dynamo samples (see below §3.2). Ideally, having access to many more independent samples than the model dimension would lead to a converged inference. This condition cannot be achieved here for the flow without reducing the dimension of the flow model state. This is less of an issue for the error term (recall that $N_U \gg N_B$). The flow perturbation is then projected onto a subspace, using
a principal component analysis (PCA). It is performed from the 71p dynamo samples (spanning 10 kyr every $\Delta t^* = 0.2$ yr), because within our set-up, this simulation is the one constraining long periods, which dominate the flow kinetic energy. We then have

$$ \begin{cases} \mathbf{u}(t) = \mathbf{u}_0 + \mathbf{P}\mathbf{v}(t) \\ \mathbf{e}(t) = \mathbf{e}_0 + \varepsilon(t) \end{cases}, \qquad (3) $$

with P a projection matrix of dimension $N_U \times N_V$, where $N_V$ is the number of retained principal components. Using $N_V = 100$, we keep $\approx 90\%$ of the surface kinetic energy for degrees $n \le 18$ in the 71p dynamo. In matrix form Eq. (2) then reads

$$ \partial_t \mathbf{b} - A(\mathbf{b})\mathbf{u}_0 - \mathbf{e}_0 = \hat{A}(\mathbf{b})\mathbf{v} + \varepsilon, \text{ with } \hat{A}(\mathbf{b}) = A(\mathbf{b})P. \qquad (4) $$

The projected flow and the error vectors are concatenated to build the augmented model state vector $\mathbf{z}$, of size $N_z = N_V + N_B$, such that

$$ \mathbf{z}^T = [\mathbf{v}^T, \varepsilon^T]. \qquad (5) $$

Eq. (2) is coupled to a stochastic model for $\mathbf{z}(t)$, as described in the next section.

## 3.2. An AR-3 forward model of the core surface dynamics

We recall below some general notions relative to stochastic differential equations, and needed for our purpose. For a general textbook on this wide topic, we refer for instance to Van Kampen (1992). Realizations of autoregressive processes of order $p$ (or AR-p) have power spectral densities (PSD) $S(f) \propto f^{-2p}$ towards high frequencies. We show in Fig. 2 examples of PSD for two flow coefficients obtained from dynamo simulations. The PSD behave approximately as $S(f) \propto f^{-6}$ for periods $T = 1/f \gtrsim \tau_A$ the Alfvén time. This observation motivates the choice of a multivariate AR-3 process in order to represent the behavior of $\mathbf{v}(t)$, and consequently $\varepsilon(t)$ for which similar slopes in the PSD have been observed. For the sake of simplicity, we follow Gillet et al. (2019) and consider for the forecast the models $\mathbf{v}$ and $\varepsilon$ as being uncorrelated.

We present below the derivation of the AR-3 model for $\mathbf{v}$. The same approach shall be followed for the error model $\varepsilon$. A continuous, linear AR-3 process takes the form

$$ d\mathbf{v}''(t) + \tilde{D}_2d\mathbf{v}'(t) + \tilde{D}_1d\mathbf{v}(t) + \tilde{D}_0\mathbf{v}(t)dt = \tilde{B}\mathbf{d}\zeta(t), \qquad (6) $$

with the notation $\mathbf{v}' = \partial_t \mathbf{v}$, $\mathbf{v}'' = \partial_t^2 \mathbf{v}$, and $\zeta$ a multivariate Wiener process. Operators $\tilde{D}_{i, \varepsilon} \in \{0, 1, 2\}$ and $\tilde{B}$ are of dimension $N_V \times N_V$. The continuous AR-3 process (6) is first transformed into an equivalent discrete AR-1 process of the form

$$ \Delta\mathbf{x}_{t_i} = \mathbf{x}_{t_{i+1}} - \mathbf{x}_{t_i} = -\Delta t \mathbf{D}\mathbf{x}_{t_i} + \sqrt{\Delta t} \mathbf{w}_{t_i}, \qquad (7) $$

for the (time-discrete) modified model state $\mathbf{x}_{t_i} = [\mathbf{v}_{t_i}^T, \mathbf{v}_{t_i}^T, \mathbf{v}_{t_i}^T]^T$, as detailed in Appendix A. $\mathbf{w}_{t_i}$ is a multivariate white noise with covariance matrix S. Here our goal is to identify matrices D and S from several dynamo series.

If the 100p dynamo can be considered as representative of the geophysical configuration, we should look at 71p as a low pass version of the dynamo state $\mathbf{x}(t)$, that is

$$ \tilde{\mathbf{x}}(t) = \int_{-\tau_w/2}^{\tau_w/2} w(\tau)\mathbf{x}(t+\tau)d\tau, \qquad (8) $$

where $w(\tau)$ is a weighting function with compact support of width $\tau_w$. In practice, we use a Blackman window of width $\tau_w = 3$ year, cutting periods shorter than $\approx \tau_A^{71p}$ for the 71p dynamo's dynamics (see Fig. 2). This latter is characterized by the presence of inertial waves, responsible for the peaks observed at high frequencies for non-zonal flow coefficients, and whose periods are too long once the dynamo is scaled into

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. PSD (in $log_{10}$ units) for the toroidal flow coefficient $\varepsilon_1^0$ (top) and $\varepsilon_3^{2c}$ (bottom), for a 10 kyr realization of the AR3 model (red). Superimposed are the median PSD of the 100p runs (black), and the PSD for the 71p simulations before (cyan) and after (grey) filtering with a 3-yr Blackman window (see text for details). The two lines correspond to -2 (dash-dotted red) and -6 (dotted red) slopes, characteristic of the PSD for respectively AR-1 and AR-3 stochastic processes towards high frequencies. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

geophysical units (Aubert and Gillet, 2021). We show in Appendix B how we consider simulation data as time-filtered representations of the geodynamo process. We next present in Appendix C how to get a least-squares estimate of the stochastic model parameters from 71p and 100p dynamo series sampled at $\Delta t^* = 0.2$ yr, under Gaussian and independence hypotheses. This provides the entries of the matrices D and S needed for the time integration of the discrete process (7) for any forecast time-step $\Delta t$, provided $\Delta t \ll 2\Delta t^*$ the Nyquist period. We also illustrate in Appendix C the main characteristics of the obtained empirical process. In particular, eigenvalues of D are complex conjugates (the model is real), and all eigenvalues have a positive real part (the model is stationary). Their imaginary parts suggest oscillators presenting period ranging from a few years to millennia, with relatively low quality factors (see Fig. C.13).

Overall, the AR-3 stochastic model obtained this way can be seen as a compromise between the information contained into the 100p and the 71p simulations. This is illustrated in Fig. 2, where we compare the PSD witnessed for the two employed dynamos and for our AR-3 model, for two large length-scale flow coefficients, one zonal ($t_0^z$) and one non zonal ($t_3^{2c}$). A qualitatively similar behavior is found for higher harmonic
degree coefficients. To further describe the behavior of the obtained AR-3 forward model, we show in Fig. 3 the temporal correlation functions obtained for some examples of core surface flow coefficients, after integration of a 10 kyr long free-run of the AR-3 model. These are compared to empirical correlation functions found from the 71p dynamo, and for the multivariate AR-1 model derived from 71p by Gillet et al. (2019) and previously used in pygeodyn. In comparison with the AR-1 model, the AR-3 model provides cross-correlations closer to those observed for the dynamo, showing a less smooth behavior when focusing on decadal lags. It also replicates for non-zonal coefficients the centennial and longer oscillations observed in 71p, with much less spurious oscillations seen in the correlation function in comparison to the AR-1 process. Finally, the cross-correlation matrix for the flow coefficients, presented in Fig. 4, proves that the spatial cross-covariances seen in the dynamo simulations are also recovered by the AR-3 model, as they were for the simpler AR-1 model (Fig. A1 in Gillet et al., 2019). The off-diagonal lines are reminiscent of the quasi-geostrophic nature of the flow in the dynamo, which couples for a given order *m* spherical harmonic coefficients of degrees *n*, *n* + 2, *n* + 4... (Pais and Jault, 2008; Gillet et al., 2015b). The above confrontation of the spatiotemporal characteristics found for the dynamo simulation and for the stochastic process derived from it validates the use of our AR-3 empirical model for mimicking the Earth's core dynamics, and its implementation into the pygeodyn data assimilation algorithm. We present in the next Section 3.3 how these developments lead to changes in the analysis step of the EnKF.

### 3.3. Analysis step for an AR-3 forward model

Forecasting an AR-3 process requires to know, at each initialization of the forward integration, the model state at 3 successive epochs. When using such a model in an EnKF (succession of forecast and analysis

<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Lower-left matrix: cross-correlation matrix for the toroidal and poloidal flow coefficients from the 71p dynamo. The coefficients are stored as follows: first the toroidal ones ($t_1^0, t_1^{1c}, t_1^{1s}, t_2^0, t_2^{1c}, t_2^{1s}, t_2^{2c}, t_2^{2s}, t_3^0, t_3^{...}$), and then the poloidal ones, stored the same way. Upper-left matrix: difference between the cross-correlation matrix from a 10 kyr free run of the AR-3 process and that from the 71p dynamo.</figcaption></figure>

steps), this request will hold for the initialization step, as well as for each time we restart the forecast from an analysis. At each analysis epoch $t_a$, we thus need to know $z_{t_a}$, $z_{t_a-\Delta t}$ and $z_{t_a-2\Delta t}$ in order to start again the forecast, with $\Delta t$ the forecast time-step. To this aim, we invert for $x_{t_a}$ (the augmented state and its successive time derivatives), and start the forecast using a Taylor expansion:

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Temporal cross-correlation for several examples of toroidal flow coefficient: $t_1^0$ (top left), $t_3^0$ (top right), $t_3^{2c}$ (bottom left), $t_4^{3s}$ (top right), for the 71p simulation (black), and 10 kyr long free runs from AR-1 (orange dashed) and AR-3 (blue dashed) core surface flow models learnt from dynamo data (see text for details). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
$$ \left\{ \begin{array}{l} \mathbf{z}_{t_a - \Delta t} = \mathbf{z}_{t_a} - \Delta t \mathbf{z}'_{t_a} + \frac{1}{2} \Delta t^2 \mathbf{z}''_{t_a} \\ \mathbf{z}_{t_a - 2\Delta t} = \mathbf{z}_{t_a} - 2\Delta t \mathbf{z}'_{t_a} + 2\Delta t^2 \mathbf{z}''_{t_a} \end{array} \right. \qquad (9) $$

Then, in order to invert the model state and its derivatives at an epoch $t = t_a$, we consider observations from three successive epochs: $t_a^-$ = $t_a - \Delta t_a$, $t_a^+$ and $t_a^+$ = $t_a + \Delta t_a$. Over the interval $[t_a^-, t_a^+]$, we forecast using the coupled system of Eqs. (4, A.4) an ensemble of $N^e$ realizations of the model state $(\mathbf{z}, \mathbf{b})$. The analysis of $\mathbf{x}_{t_a}$ is then performed in two steps: first for $\mathbf{b}$ from MF data, and second for the augmented state $\mathbf{z}$ from SV data, as in previous versions of the pygeodyn assimilation tool (e.g., Gillet et al., 2019). The details of the adaptation of the analysis step to the AR-3 forward model are provided in Appendix D, in particular regarding the modification of the matrices entering the Kalman gain matrix such as the observation operator or the forecast covariance matrix. Following Istas et al. (2023), this latter is empirically estimated from the ensemble of forecast trajectories using the Graphical LASSO method (Friedman et al., 2007).

## 4. Results

The reanalyses presented below have been obtained using an ensemble of dimension $N^e = 400$, an analysis step $\Delta t_a = 1$ yr and a forecast time-step $\Delta t = 1$ month $\simeq 0.0833$ yr. This value is to be compared with the Nyquist period $2\Delta t^* = 0.4$ yr. The above assimilation scheme has been applied to Gauss coefficients data from the COV-OBS-x2 geomagnetic field model (Huder et al., 2020) and from the CHAOS-7 model (Finlay et al., 2020). For the analysis of MF and SV data (respectively Eqs. (D.2 and D.6)), we consider an ensemble of noisy observations, as it was the case for previous versions of pygeodyn (see Gillet et al., 2019). For COV-OBS-x2 these correspond to an ensemble of $N^e$ realization of the model (obtained following Gillet et al., 2013). When considering CHAOS-7, we generate an ensemble of $N^e$ noisy versions of this model, as follows: for each realization, and at all considered analysis epochs, we add to CHAOS-7 the differences between one Kalmag realization and the ensemble average Kalmag model. We first present in Section 4.1 an evaluation of the ability of the above scheme to predict the observed SV. Next we show in Section 4.2 some examples of wave recovery in several period ranges, from the reanalysis of Gauss coefficients data.

### 4.1. Prediction of the secular variation

In order to validate the behavior of the scheme described in Section 3, we set-up the following test. The AR-3 algorithm is used to perform a reanalysis from COV-OBS-x2 data over 1880–2020. The initial condition is a random draw within the ensemble of accessible dynamo states. Starting from this reanalysis, 10 yr long hindcasts experiments are operated, beginning at $t = t_0$ every two years for $t_0 \in [1960, 2020]$. For each of the hindcasts, the reanalysis at $t_0$ is re-operated with the left-hand scheme involving Eqs. (D.7–D.9), so that the ‘predicted’ trajectory does not account for observations at $t > t_0$. The predicted trajectories are then compared with the observed one from COV-OBS-x2. The goal of this exercise is to measure whether the magnetic field predicted this way behaves in a geophysical manner. Because of the stochastic nature of the forward model, we do not expect the ensemble average SV prediction $\langle \partial_t g_n^m \rangle(t)$ for the Gauss coefficient of degree $n$ and order $m$ to closely match the observed one, denoted $\partial_t g_n^{mobs}(t)$ – here the mean COV-OBS-x2 coefficients. However, this latter should (at least in most of the cases) be well represented by the spread within the ensemble of SV predictions. To this purpose, we introduce for each Gauss coefficient a normalized measure of the SV prediction misfit,

$$ \epsilon_n^m(t_0 \tau) = \frac{\partial_t g_n^{mobs}(t_0 + \tau) - \langle g_n^m \rangle(t_0 + \tau)}{\sigma_n^m(t_0 + \tau)}, \qquad (10) $$

$$ \text{where } \sigma_n^m(t)^2 = \frac{1}{N^e} \sum_{k=1}^{N^e} \left( \langle \dot{g}_n^m \rangle(t) - \dot{g}_{n,k}^m(t) \right)^2 $$

and $k$ indicates the ensemble member.

The normalized misfit (10) is evaluated at $t - t_0 = 5$ yr for the four experiments starting in 1970, 1980, 1990 and 2000. Fig. 5 indicates no systematic outliers. Neither does it show any obvious misfit patterns as a function of degree and/or order (e.g. towards zonal or sectorial coefficients, or towards short length-scales). Over all degrees and orders, and over all four considered epochs, the $\epsilon_n^m$ (5 yr) present a standard deviation of about 0.40. This points to a behavior of the ensemble of model trajectories consistent with the observed field changes. This is further evidenced with some examples of Gauss coefficients series, as illustrated in Fig. 6. Despite the occurrence of some rare off-set trajectories (for instance $dh_2^1/dt$), we visually check that the observed SV is most of the time accounted for within one standard deviation, and for all considered prediction era. This is the case on all length-scales, including the largest ones (the dipole field), whose large rapid variations witnessed over the past decades are not easily captured naturally by dynamical models (Aubert, 2023). The miss-match found for $dh_2^1/dt$ may be related to the intense values observed for this particular coefficient, which may not fall within the usual range of values for the AR-3 model. It may also be in link with the eccentric dipole evolution. Indeed, the energy associated with changes in $g_2^1$ and $h_2^1$ is by construction almost stationary in the path set-up, in relation with an azimuthal westward drift of the eccentric gyre on millennial time-scales (Aubert et al., 2013). Conversely, geophysical observations suggest significant variations of this quantity over the past decades (see Fig. 6) or millennia (Nilsson et al., 2022). This may highlight a discrepancy between the dynamo model considered here and the geophysical configuration, as discussed in §5.

We further evaluate the prediction ability of our AR-3 model, by estimating the r.m.s. radial field prediction error over the Earth’s surface $\Sigma$, of reference radius $a = 6371.2$ km, for a forecast started in $t_0$ over a time window $\tau$. It is directly calculated from Gauss coefficients as (Lowes, 1966)

$$ \begin{aligned} \Delta B_r^{AR3}(t_0 \tau) &= \left[ \frac{1}{\Sigma} \int_{\Sigma} \left[ B_r^{obs}(t_0 + \tau) - \langle B_r^{AR3} \rangle(t_0 + \tau) \right]^2 d\Sigma \right]^{1/2} \\ &= \left[ \sum_{n=1}^{L_8} \frac{(n+1)^2}{(2n+1)} \sum_{m=0}^{n} (\Delta g_n^{mAR3}(t_0 + \tau))^2 + (\Delta h_n^{mAR3}(t_0 + \tau))^2 \right]^{1/2} \end{aligned} \qquad (11) $$

where $\Delta g_n^{mAR3}(t) = g_n^{mobs}(t) - g_n^{mAR3}(t)$, and the upperscripts “obs” and “AR3” stand respectively for the observed field and the AR-3 prediction. We only consider here the radial component, as this contains all the information for a potential magnetic field vector.

The median value of the above diagnostic for time increments $\tau = 2$, 5 and 10 years are given in Table 2, where it is compared to the similar diagnostic $\Delta B_r^{LE}$ (resp. $\Delta B_r^{AR1}$) obtained by replacing the AR-3 prediction in (11) by the linear extrapolation $B_r^{LE}(t_0 + \tau) = B_r^{obs}(t_0) + \tau \partial_t B_r^{obs}(t_0)$ at the starting epoch (resp. the prediction using the AR-1 forward scheme as in Gillet et al., 2019). With a median around 7 nT, AR-3 MF prediction errors after 2 years are lower by $\approx 35\%$ (resp. 45 %) in comparison with those based on a linear extrapolation (resp. on AR-1 forecasts). After 5 yr the improvement w.r.t. the linear extrapolation or the AR-1 forecasts is only about 10%. It furthermore presents a large variability, with the best scores as low as $\approx 30$ nT and the worst up to 80 nT, over the considered period (against a range typically between 35 and 70 nT for both the AR-1 and the linear extrapolation). The quality of the AR-3 prediction depends indeed on how representative of the up-coming era is the
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Normalized SV misfit (10) for 5 yr predictions starting in 1970 (top left), 1980 (top right), 1990 (bottom left) and 2000 (bottom right). In green values away from ±2σ. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

curvature at the last analysis epoch. The occurrence of successive jerks (or recurrent SA pulses), as witnessed over the recent satellite era, render 5 yr AR-3 forecasts particularly sensitive. After 10 yr lags, the AR-3 forecasts bring no improvement, meaning that the information from the curvature of $B_r(t = t_0)$ has been lost. The fact that all three predictions present similar scores for lags $\tau > T_A$ is expected. First, our AR-1 and AR-3 models are both imprinted by the same dynamo data, so they will share strong similarities on timescales $T > T_A$. Second, for a signal of time-scale $T \gg \tau > T_A$, the average SV trajectory in our AR processes will evolve approximately as $\exp(-\tau/T)$. The SV signal is dominated by long periods (for which $\tau/T \ll 1$), so that the average SV forecast from our stochastic models is necessarily close to a linear extrapolation. The main quality of our AR-3 model is to provide a realistic estimate of the trajectory spread on a broad range of time-scales. Improving the average prediction would require to introduce further deterministic information and/or some knowledge about the current geodynamo state, as discussed in §5.

The above 5 yr MF prediction errors have to be compared with a similar diagnostic operated with alternative tools. Fournier et al. (2021) compared forecasts using geodynamo simulations under various configurations and a linear extrapolation of the MF. For 5 yr predictions starting in 2005, 2009 and 2015, they found (see their Table 1) MF prediction errors typically within 60 to 70 nT/yr for linear MF extrapolations, with comparable values (within typically ±3 nT) for dynamo based predictions. Similarly Sanchez et al. (2020) found prediction errors for geodynamo-based hindcasts slightly larger than when assuming a linear extrapolation. Alternatively, Beggan and Whaler (2018) obtained $\Delta B_r \approx 70$ nT when performing 5 yr predictions starting in 2010 based on inverted core-flows. Despite a large variability from epoch to epoch, the above comparison is overall favorable to our AR-3 tool. The significant improvement for short lags (2 yr) is even encouraging for reducing the bias in the flow trajectories for reanalyses purposes.

## 4.2. Decadal and interannual motions within Earth's core

We now turn to the second focus, which treats of transient (inter-annual and decadal) core surface motions. The algorithm described in Section 3 is first applied to perform a reanalysis of COV-OBS-x2 data over the interval 1880–2020. The ensemble of solutions in 1999 is then used as an initial condition for a reanalysis of CHAOS-7 data over 1999–2023. The two series are merged to provide a continuous series over 1880–2023, using a sine taper of half-period $\tau_{tap} = 2020 - 1999 = 21$ yr (the time-span over which both models overlap). The use of such a taper allows (i) to analyse the longest available era (given the two models at hand), thus refining the spectral density when performing signal analysis, and (ii) to prevent ad hoc high frequency content around 1999 that would show up while merging abruptly the two series, which would consequently pollute our interpretations. We illustrate the results with band-pass filtered core surface flow motions in the three period ranges $[T_{min}, T_{max}] = [10, 16]$ yr, $[5, 8]$ yr and $[2.75, 4.4]$ yr, all presenting the same relative period width $(T_{max} - T_{min})/\sqrt{T_{max}T_{min}} \approx 0.47$. We use a causal Butterworth filter of order 2. The choice of such windows is motivated by Fig. 8, where we witness bumps around approximately 12.5, 6.5 and 3.5 yr periods in the PSD of the azimuthal velocity $u_\phi$ at the core surface equator, where motions are the largest (e.g. Gillet et al., 2015b, 2022b; Kloss and Finlay, 2019). At similar periods enhanced signals are also witnessed in ground-based magnetic time series recorded at low latitudes (Lesur et al., 2022, see their Fig. 3 for the Y component at Hermanus, South-Africa). We focus below on the azimuthal velocity $u_\phi$ for the ensemble mean, with time-longitude and time-latitude diagrams. In each of the considered period range, we show (see Figs. 9, 10 and 11) a time-longitude diagram at the equator, and time-latitude diagrams at two longitudes.

We first discuss the slowest time-scales around 12.5 yr period. Westward propagating patterns appear starting already from ~1920, with an amplitude of the order of 3 km/yr. An azimuthal propagation shows up under the Pacific hemisphere at longitudes from ~60°W to
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** 10 yr SV predictions (in nT/yr) for several examples of Gauss coefficients, starting every 4 yr from 1960 to 2020, superimposed with COV-OBS-x2 (dashed black). Thick lines: ensemble mean. Shaded area: ±1σ (spread)</figcaption></figure>

~ 100' E, with a speed of the order of 1000 km/yr. It appears clearer than on similar diagrams obtained using an AR-1 forward model (Istas et al., 2023, see their Fig. 12). The patterns are modulated over decadal time-scales, with a lower magnitude around 1960, as already noticed by

Istas et al.. For the last 15 yr the signature appears faint, consequence of the bandpass filtering towards endpoints. The decadal modulation is obvious on time-latitudes diagrams, as is the latitudinal propagation. The latter corresponds to an average cylindrical radial speed of ~ 100
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Comparisons of the r.m.s. prediction MF errors $\Delta B_r^{AR2}$, $\Delta B_r^{LE}$ and $\Delta B_r^{AR1}$ for lags $\tau = 2, 5$ and $10$ yr (in nT), for hincasts starting every 2 yr from 1960 until $2020 - \tau$.</figcaption></figure>

**Table 2**

For various values of the lag $\tau$, and from trajectories obtained using the AR-3 model, the AR-1 model or a linear extrapolation (LE), the median values (over hincasts from Fig. 7) and the standard deviation for $\Delta B_r$ (see Eq. (11)) the r.m.s. MF prediction error at the Earth's surface (in nT).

<table><thead><tr><th>forecast</th><th>LE</th><th>AR-1</th><th>AR-3</th></tr></thead><tbody><tr><td>τ = 2 yr</td><td>10.2 ± 1.8</td><td>12.1 ± 2.4</td><td>6.8 ± 1.5</td></tr><tr><td>τ = 5 yr</td><td>47.3 ± 7.1</td><td>47.8 ± 8.4</td><td>42.0 ± 13.6</td></tr><tr><td>τ = 10 yr</td><td>147.9 ± 20.9</td><td>148.2 ± 21.0</td><td>148.9 ± 36.1</td></tr></tbody></table>

km/yr, for motions predominantly equatorially symmetric. The wave-forms vary with longitude, as it is the case for QG-MC eigenmodes travelling on a non-axisymmetric background magnetic field (Gerick et al., 2021; Gillet et al., 2022b). All these characteristics recall wave-trains extracted from the 71p numerical simulation (Schwaiger et al., 2024).

We recover at periods around 6.5 yr the westward wave-like motions first highlighted by Gillet et al. (2022b), and then confirmed by Ropp and Lesur (2023) and Li et al. (2024). If their intensity is lower prior to the modern satellite era, our new analysis with the AR-3 model allows to

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** PSD (in $log_{10}$ units) of the azimuthal flow acceleration $\partial_t u_\phi$ at the core equator, averaged over longitudes. The three shaded bands correspond to the three bandwith documented in Section 4.2. In order to base our estimate on a period covered by enough observation sites, the PSD has been obtained starting from 1958 (the International Geophysical Year), until 2023.</figcaption></figure>

trace the actual wave train as far back as the 1990's (see Fig. 10). Furthermore, it reveals a similar coherent wave-train occurring during the 1960's-70's, which was not caught by our previous reanalyses using AR-1 forecasts (Istas et al., 2023). This period coincides with the global coverage offered by the POGO satellite mission (1965–1970, Cain and Sweeney, 1973). Despite less accurate (and intensity only) measurements, these data might allow the detection of planetary-scale waves in the core. As for the 12.5 yr signal, the 6.5 yr waves are modulated in longitude. The amplitude of the wave-train witnessed over the recent era is particularly intense in the central Pacific, where it reaches velocities as high as ~4 km/yr (to be compared with an amplitude about twice lower over most of the equatorial belt). We do not recover such a localized extremum for the signal detected in the 1960's, although this might be related to a weaker observational constraint. The time-latitude diagrams show latitudinal propagation, which would correspond to a cylindrical radial average speed in the range ≈ 150 – 250 km/yr, about twice larger than for the 12.5 yr signal, and possibly weaker while approaching the equator.

We finally focus on the signal with periods around 3.5 yr, where motions are of the order of 0.5 km/yr. This is about twice less in comparison with motions recovered on the same time-scales from GVO data by Istas et al. (2023), possibly in link with (i) the smoother effect from the AR-3 scheme, and (ii) a larger time variability contained into GVO data (in comparison with magnetic Gauss coefficients series). Apart from a localized signal around 1970 near 120 E, those are in general more intense over the satellite era, when the taper used to merge the COV-OBS-x2 and CHAOS-7 reanalyses favors the latter (see Fig. 11). This makes sense since by construction COV-OBS-x2 field changes are progressively damped on periods shorter than ≈ 3 yr (due to the projection onto splines with 2 yr know spacing) while large-scale SV patterns from CHAOS-7 are less constrained on short time-scales, thus allowing a stronger power. No obvious westward propagation is found. One may see eastward moving patterns (starting from ~ 120° E), although this observation remains elusive. These would correspond to a rapid azimuthal propagation (~ 4000 km/yr at the core equator), although it may as well be (i) the consequence of the band-pass filtering, and (ii) alternatively interpreted as standing waves. As already highlighted by Istas et al. (2023), transient motions in this period range are primarily equatorially symmetric, and present obvious latitudinal propagation at a speed of the order of ~ 300 – 400 km/yr, i.e. faster than observed for patterns isolated on 6.5-yr time-scales. A similar coherent propagation is also recovered in the eastern hemisphere around 1970. Localized in longitude, it recalls the behavior of QG Alfvén waves detected in numerical simulations (Aubert and Finlay, 2019; Aubert et al., 2022). This
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Band-pass filtered azimuthal core surface flow (in km/yr) in the period range between 10 and 16 years. Time-longitude diagram at the core equator (top). Time-latitude diagram at 120° E (middle) and 150° W (bottom).</figcaption></figure>

latter feature was not apparent in the reanalysis of the same data by Istas et al. while considering an AR-1 forward model (see their Fig. 11). This suggest, as for the 6.5 yr signal detected in the 1960’s–70’s, an improvement from our new implementation based on the AR-3 scheme.

One common observation to all three investigated period ranges is the more complex structure of the azimuthal flow patterns at longitudes around 90° E. All signals seem to travel faster in the azimuthal direction than in the cylindrical direction. The recovered propagation speed decreases with the considered time-scales. In all cases phase speed of the inverted patterns is much larger than the fluid velocity. All these characteristics are shared with the analysis performed by Schwaiger et al. (2024) from the 71p dynamo, reinforcing the geophysical relevance of the path dynamos (Aubert et al., 2017) for this aspect of the dynamics.

# 5. Summary and discussion

## 5.1. Geomagnetic data assimilation based on an Earth-like stochastic model

We have presented a new scheme for building stochastic dynamical models of the dynamics taking place at the surface of Earth’s core. It allows to benefit from the information contained on short time-scales in geodynamo simulations run in Earth-like configurations (Aubert, 2023). Such a reduced model, based on multivariate autoregressive processes of order 3, not only replicates the spatiotemporal statistics on decadal and longer periods, as offered by previous lower-order models (Barrois et al., 2017; Gillet et al., 2019; Baerenzung et al., 2020), but it also mimics the cut-off in the frequency spectrum of the core surface magnetic field and flow observed in dynamo simulations for periods shorter than ≈ τ<sub>A</sub>, and related to the magnetic dissipation (Aubert and Gillet, 2021). Integrated
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Band-pass filtered azimuthal core surface flow (in km/yr) in the period range between 5 and 8 years. Time-longitude diagram at the core equator (top). Time-latitude diagrams at 170° E (middle) and 80° W (bottom).</figcaption></figure>

to the pygeodyn geomagnetic assimilation tool (Huder et al., 2019) as a kind of stochastic ensemble Kalman smoother, we illustrate its use for two purposes: the prediction of the geomagnetic field, and the recovery of transient wave-like motions.

By construction, the tool does not offer deterministic forecasts. Meanwhile, we show that it provides realistic range of SV forecasts, and this on all observable length-scales. It furthermore supersedes AR-1 or linear extrapolations on short periods (2 yr), thus reducing spurious short-term oscillations in reanalysed core flow series. Applied to 5 yr predictions, our tool leads in average to a ≈ 10% decrease in the r.m.s. radial magnetic field forecast error, in comparison with a linear extrapolation. Despite a large variability found in the 5 yr prediction scores from one epoch to the other, we propose with our AR-3 scheme an alternative to the forecasts based on geodynamo-based reanalyses (Fournier et al., 2021) or involving inverted core flows as operated by Beggan and Whaler (2018). An immediate application of our method would be the production of 5 yr predictions in the context of the

international geomagnetic reference field model (Alken et al., 2021). It could be performed either from the reanalysis of Gauss coefficients series, as illustrated here, or of virtual observatory data (Hammer et al., 2021), since both kinds of observations can be considered when using the pygeodyn software (Istas et al., 2023).

## 5.2. Wave motions on various time-scales

We performed a reanalysis of Gauss coefficients data from the COV-OBS-x2 and CHAOS-7 field models. We observe an enhanced flow intensity near the equatorial belt on periods around 12.5, 6.5 and 3.5 yr. These are reminiscent of similar spectral lines found in ground-based series recorded at low latitudes (Lesur et al., 2022). The signal previously isolated around 6.5 yr periods over the past two decades where continuous satellite coverage is available (Gillet et al., 2022b; Ropp and Lesur, 2023; Li et al., 2024) is confirmed by our reanalysis. But our AR-3 scheme helped extract a new evidence for interannual wave-like
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Band-pass filtered azimuthal core surface flow (in km/yr) in the period range between 2.75 and 4.5 years. Time-longitude diagram at the core equator (top). Time-latitude diagrams at 160° E (middle) and 120° E (bottom).</figcaption></figure>

motions prior to the recent satellite era, in particular during the 1960's-70's around both 3.5 and 6.5 yr periods. This discovery may be related to the presence of continuous intensity records from space during the second half of the 1960's, but also to the occurrence of particularly intense SV changes in conjunction with the jerk around 1969 (e.g. Mandea et al., 2010; Jackson and Finlay, 2015). Indeed, the flow recovery is sensitive to the data quality and coverage, but as well to the magnitude of the flow changes, as characterized with a synthetic study by Schwaiger et al. (2024). This is particularly true for the dynamics on short periods, for which motions are relatively less intense. Around the period of 6.5 yr, coherent flow patterns are triggered westward from below central America starting from about 1960. These also present equatorward propagation from mid-latitudes below the Americas. Such an evolution is to be related to the strongest SV changes at the CMB seen under Northern to Central America around the 1969 jerk (Blangsbøll et al., 2022). At the same epoch, flow changes at periods around 3.5 yr, although present over all longitudes, are the strongest in a more narrow

band of longitudes under eastern Asia around 1970. These patterns may suggest the presence of two distinct events, as discussed by Blangsbøll et al. (2022).

On longer periods around 12.5 yr, we extracted modulated flow oscillations as far back as 1920. These were present, although with much less spatial coherence, in the previous reanalysis by Istas et al. (2023), also involving pygeodyn but considering for the forecast step an order 1 stochastic model. Again the smoothing effect coming with the AR-3 scheme helps isolate transient events. This decadal oscillation might have to do with the ~10 yr mode extracted by Li et al. (2024) using a dynamic mode decomposition (the discrepancy between 12.5 and 10 yr periods may come from the limited time-span (22 yr) considered by Li et al.). Performing a PCA on the band-pass filtered ensemble average flow series, we isolate a single component responsible for about 40 % (resp. 36 %) of the core surface kinetic energy in the period range [5, 8] yr (resp. [10, 16] yr). These are shown in Fig. 12. The associated 6.5 yr oscillation intensifies over the past 30 yr, where it amounts to 55 % of
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 12.** Mollweide projection of the azimuthal velocity $u_\phi$ (colormap, normalized to $\max(|u_\phi|)$) and streamlines, for the first PC obtained from band-pass filtered core surface flow series, in the period ranges [10, 16] yr (top), [5, 8] yr (middle) and [2.75, 4.4] yr (bottom). The PCA has been performed over the period 1920–2023 for the longer periods, and 1958–2023 (i.e. starting from the International Geophysical Year) for the latter two. The first PC corresponds respectively to 36 %, 40 % and 18 % of the core surface kinetic energy around respectively 12.5, 6.5 and 3.5 yr periods. Focusing on period 1990–2023, it accounts for 55 % (resp. 27 %) of the flow variance at 6.5 yr (resp. 3.5 yr).</figcaption></figure>

the flow variance, while the 12.5 yr oscillation dominates with an almost constant amplitude over the past century or so. This difference in behavior could arise from a natural variability in the core dynamics. It may also be partly related to flow motions being more accurately recovered at decadal than interannual time-scales, as for a fixed flow magnitude, changes on longer periods are associated with a larger magnetic signal (Schwaiger et al., 2024). The two most energetic PC shown in Fig. 12 are predominantly equatorially symmetric, with a stronger azimuthal than orthRADial velocity. More intense near the equator, they are dominated by low harmonic orders and present a relatively short latitudinal extent. All these properties are shared by QG-MC modes (Gerick et al., 2021; Gillet et al., 2022b). The large fraction of kinetic energy carried by a single PC is another evidence for the presence of QG-MC waves at various periods, as observed in numerical

simulations (Gillet et al., 2022a; Schwaiger et al., 2024). The situation is less clear for 3.5 yr signals, where the flow energy is more spread within the various PC. At these short periods the most energetic PC is still mainly equatorially symmetric and more intense towards the equator, although we notice relatively larger orthoradial flows. One could speculate on the presence of higher order QG-MC waves, although in this period range closer to the Alfvén time, the dynamics might be inprinted by localized QG Alfvén waves (Aubert and Finlay, 2019; Gillet et al., 2022b; Aubert et al., 2022).

We have seen similarities between inverted oscillatory patterns and waves isolated in the 71p dynamo simulation. One may wonder how much our choice of dynamo prior impacts the reconstruction. Regarding waves of period ≈ 6 yr, we recall that similar patterns have been imaged over the satellite era when using different prior information and inversion set-up by Ropp and Lesur (2023) and Li et al. (2024). For earlier epochs where the observational constraint is less tight, the dynamo prior will necessarily play a larger role, as it reduces non-uniqueness when and where data constraints are sparser. Meanwhile, the phase and amplitude of the recovered decadal flow changes are imposed by the SV data (otherwise the average model would tend to decay towards zero, with the ensemble average patterns presenting weak amplitudes). Schwaiger et al. (2024) have shown that it is possible to recover decadal motions prior to the modern satellite era, provided the signal is intense enough. Nevertheless, definitely determining the nature of transient flow motions calls for assessing the sensitivity of inverted core flows to the choice of dynamo prior.

In any case, the growing evidence for a variety of hydro-magnetic waves in Earth's core opens a window towards a more deterministic description of the geomagnetic signal, and with it the possibility of improved predictions of the field evolution. One may envision two strategies in this direction. The first is based on an empirical description of the waves signature at the core surface, as initiated here or in the work for instance by Li et al. (2024) or Wardinski et al. (2020). One would next need realistic models of the forcing, which may arise with stochastic descriptions as in the present work (see also Buffett and Knezek, 2018, for another family of waves). Alternatively, one could rely on deterministic models of the whole fluid cavity, considering either eigenmodes as a basis for projecting the magnetic signal in space and time (Gerick et al., 2021; Luo et al., 2022), or the forward integration of primitive equations at extreme conditions (Aubert, 2023), possibly in a linearized set-up (Barrois and Aubert, 2024). However, this avenue may require the knowledge of a minimal complexity for the base magnetic field. This issue requests dedicated sensitivity studies, as initiated by Gerick and Livermore (2024). Furthermore, calculating eigenmodes or solving the dynamo equations on complex base states will come at an enhanced numerical cost.

### 5.3. Assessing geodynamical issues

Finally, our AR-3 modeling strategy may help address geodynamical issues, by confronting geomagnetic data to various dynamo priors (via pygeodyn) for hypothesis testing. We give below two examples. First, we have seen evidence for possible systematic bias between the observed and predicted SV for the coefficient $h_2^1$. This likely highlights some discrepancy between the considered geodynamo model and the Earth. Indeed, the path dynamo set-up does not favor for this coefficient changes as large as 30 nT/yr on decadal time-scales, as seen in Fig. 6. One may be tempted to associate this difference in behavior to the actual presence of the South Atlantic Anomaly. There is some evidence, from paleomagnetic data, for a relationship between low latitude intensity anomalies and minima of the dipole field (Nilsson et al., 2022). There is furthermore some connection between the actual dipole decay and the eccentric westward gyre (Finlay et al., 2016), while this latter appears related to the SAA (Aubert, 2015). It thus seems reasonable to connect the mismatch in the $h_2^1$ SV coefficient to the presence of a low latitude,
low intensity anomaly presenting a behavior seen as unusual when considered within the frame of the simulated dynamo.

Second, the analysis of length-of-day (LOD) records has highlighted two interannual spectral lines around 6 and 8.5 yr (Duan and Huang, 2020). The origin of these latter is found to be primarily associated with axisymmetric, axially invariant motions within Earth's core (Istas et al., 2023; Rosat and Gillet, 2023). However, the underlying physical mechanism remains debated. A natural candidate comes with the propagation of purely axisymmetric torsional Alfvén waves (Gillet et al., 2010). Alternatively QG-MC waves, although mostly non-axisymmetric, may also carry angular momentum as they couple all spherical harmonic orders when propagating over a non-zonal background field (potentially carrying a $m=0$ by-product contribution, as discussed by Istas et al., 2023). The analysis presented here shows evidence for enhanced wave-like patterns, mainly non-zonal, around 6.5 yr but not at 8.5 yr (see also Li et al., 2024). This observation does not favor a scenario where QG-MC waves would be the source of the 8.5 yr LOD signal. Given the currently available frequency resolution (prescribed by the time-span covered by accurate magnetic data), we do not exclude QG-MC waves for the 6 yr LOD signal – although both observed LOD oscillations could be associated with torsional waves. The excitation of these waves may involve the azimuthally averaged axial magnetic torque in the bulk of the core (Teed et al., 2014; Gillet et al., 2017), and/or a mantle-inner core gravitational coupling (MICG, see Mound and Buffett, 2003, 2005). To further disentangle the several possibilities, one would need high-resolution geodynamo simulations in various geodynamical conditions (e.g. with/without MICG).

The updated AR-3 scheme implemented here into the pygeodyn tool, by extracting statistical information on period as short as a fraction of $\tau_A$, offers a way to confront the simulated rapid dynamics to geophysical (magnetic and LOD) observables. Geodynamo simulations have put forward natural temporal and longitudinal modulations of wave-like motions within Earth's core (Aubert et al., 2022; Schwaiger et al., 2024). So did we with the reconstruction of wave-trains from both recent and early satellite and ground-based magnetic data. These modulations may bring extra constraints on the geodynamo processes, such as dissipation within the fluid core cavity or at its boundary (via the quality factor of the QG-MC or torsional oscillators), thus possibly involving deep mantle properties (electrical conductivity, density anomalies). Long-lived global coverage from space is crucial to best document the natural variability within Earth's core, as offered by a

Swarm mission extended as much as possible, and the Macau Science Satellites (Zhang, 2022).

## Authors statement

The AR3 algorithm was first designed by NG, with the help of POA. NG designed a python prototype, which has been implemented into the pygeodyn reanalysis tool by FD, then refined by RC. Reanalyses have been run by FD and RC. NG performed the post-processing. Numerical dynamo data have been provided by JA. The manuscript has been primarily written by NG, and proof-read by all co-authors.

## CRediT authorship contribution statement

**N. Gillet:** Writing – original draft, Visualization, Validation, Super-vision, Software, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization.
**F. Dall'Asta:** Writing – review & editing, Validation, Software, Methodology, Investigation, Data curation. **P.-O. Amblard:** Writing – review & editing, Methodology, Formal analysis. **R. Claveau:** Writing – review & editing, Validation, Software, Methodology, Investigation, Data curation. **J. Aubert:** Writing – review & editing, Validation, Data curation.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

We thank the two anonymous referees for their constructive comments. This work has been funded by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (GRACEFUL Synergy Grant agreement No 855677) and by ESA in the framework of EO Science for Society, through contract 4000127193/19/NL/IA (Swarm +4D Deep Earth: Core). This project was supported by CNES as an application of the Swarm mission of ESA. NG, FD and RC are part of Labex OSUG@2020 (ANR10 LABX56).

## Appendix A. AR-1 discrete formulation of an AR-3 process

The continuous AR-3 process (6) can be equivalently written in the form of the AR-1 model

$$d\mathbf{x}(t) = - \begin{bmatrix} 0 & -I & 0 \\ 0 & 0 & -I \\ \bar{D}_0 & \bar{D}_1 & \bar{D}_2 \end{bmatrix} \mathbf{x}(t)dt + \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \\ d\zeta(t) \end{bmatrix}, \text{ with } \mathbf{x}(t) = \begin{bmatrix} \mathbf{v}(t) \\ \mathbf{v}^{\ddot{}}(t) \\ \mathbf{v}^{\dddot{}}(t) \end{bmatrix}. \quad (A.1)$$

A 1st order differentiation gives

$$\mathbf{x}_{t_i} = \begin{bmatrix} \frac{\mathbf{v}_{t_{i+1}}}{\Delta t} \\ \frac{\mathbf{v}_{t_{i+1}} - \mathbf{v}_{t_i}}{\Delta t} \\ \frac{\mathbf{v}_{t_{i+2}} - 2\mathbf{v}_{t_{i+1}} + \mathbf{v}_{t_i}}{\Delta t^2} \end{bmatrix} \quad \text{and} \quad \Delta \mathbf{x}_{t_i} = \begin{bmatrix} \frac{\mathbf{v}_{t_{i+1}} - \mathbf{v}_{t_i}}{\Delta t} \\ \frac{\mathbf{v}_{t_{i+2}} - 2\mathbf{v}_{t_{i+1}} + \mathbf{v}_{t_i}}{\Delta t} \\ \frac{\mathbf{v}_{t_{i+3}} - 3\mathbf{v}_{t_{i+2}} + 3\mathbf{v}_{t_{i+1}} - \mathbf{v}_{t_i}}{\Delta t^2} \end{bmatrix}. \quad (A.2)$$

Discretized using an Euler-Maruyama scheme, the sampled process then becomes

$$\Delta \mathbf{x}_{t_i} = -\Delta t \begin{bmatrix} 0 & -I & 0 \\ 0 & 0 & -I \\ \bar{D}_0 & \bar{D}_1 & \bar{D}_2 \end{bmatrix} \mathbf{x}_{t_i} + \sqrt{\Delta t} \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \\ \eta_{t_i} \end{bmatrix}, \quad (A.3)$$
with $\Delta\mathbf{x}_{ti} = \mathbf{x}_{ti+1} - \mathbf{x}_{ti}$ and $\eta_{ti}$ a white multivariate random noise vector, centered and whose covariance matrix is the identity matrix. This discrete process (A.3) may equivalently be written as a function of the $\mathbf{v}_{ti}$ at 4 successive epochs:

$$ \mathbf{v}_{ti+3} + D_2 \mathbf{v}_{ti+2} + D_1 \mathbf{v}_{ti+1} + D_0 \mathbf{v}_{ti} = B \eta_{ti} \quad (A.4) $$

The link between matrices entering (A.3) and (A.4) is

$$ \tilde{D}_2 = \frac{D_2 + 3I}{\Delta t}, \quad \tilde{D}_1 = \frac{2D_2 + D_1 + 3I}{\Delta t^2}, \quad \tilde{D}_0 = \frac{D_2 + D_1 + D_0 + I}{\Delta t^3}, \quad \tilde{B} = \frac{B}{\Delta t^{5/2}} \quad (A.5) $$

In practice we will time-step (A.4). To this purpose we invert (A.5) to get

$$ \left\{ \begin{array}{l} D_2 = \Delta t \tilde{D}_2 - 3I \\ D_1 = \Delta t^2 \tilde{D}_1 - 2\Delta t \tilde{D}_2 + 3I \\ D_0 = \Delta t^3 \tilde{D}_0 - \Delta t^2 \tilde{D}_1 + \Delta t \tilde{D}_2 - I \\ B = \Delta t^{5/2} \tilde{B} \end{array} \right. \quad (A.6) $$

The next step requires to fit matrices $\tilde{D}_2, \tilde{D}_1, \tilde{D}_0, \tilde{B}$ using geodynamo samples, before integrating with an Euler scheme the radial induction Eq. (4), coupled to the evolution Eq. (A.4) for $\mathbf{v}(t)$, and the equation for $\varepsilon(t)$ obtained in a similar manner. To this purpose, we shall rewrite Eq. (A.3) as

$$ \Delta \mathbf{x}_{ti} = - \Delta t \mathbf{D} \mathbf{x}_{ti} + \sqrt{\Delta t} \mathbf{w}_{ti}, \quad (A.7) $$

with $\mathbf{w}_{ti} = B\eta_{ti}$ a multivariate white noise of covariance matrix $S = BB^T$. For the process to be stationary, all eigenvalues of D should have their real part positive.

## Appendix B. Considering simulation data as time-filtered representations of the geodynamo process

Simulation data are sampled every $\Delta t^* = 0.2$ yr. The time-filtered model state (8) discretized at this rate is written

$$ \tilde{\mathbf{x}}_{ti} = \sum_{j=-N_w}^{N_w} w_j \mathbf{x}_{ti,j}, \quad \text{with } N_w = (\tau_w / \Delta t^* + 1) / 2, \quad (B.1) $$

where the $w_j$ are the weights of the considered window once discretized. For $\Delta t^* = 0.2$ yr, the discretized filter that enters (B.1) is based on $2N_w + 1 = 17$ successive samples. The discrete process (A.7), once filtered in time and sampled every $\Delta t^*$, becomes

$$ \Delta \tilde{\mathbf{x}}_{ti} = -\Delta t^* D \tilde{\mathbf{x}}_{ti} + \sqrt{\Delta t^*} \tilde{\mathbf{w}}_{ti}, \quad \text{with } \tilde{\mathbf{w}}_{ti} = \sum_{j=-N_w}^{N_w} w_j \mathbf{w}_{ti,j}. \quad (B.2) $$

$\mathbf{w}_t$ being a white noise, the covariance matrix for $\tilde{\mathbf{w}}_t$ is simply

$$ \tilde{S} = \mathbb{E} \left( \tilde{\mathbf{w}}_t \tilde{\mathbf{w}}_t^T \right) = \mu_w S, \quad \text{with} \quad \mu_w = \sum_{j=-N_w}^{N_w} w_j^2. \quad (B.3) $$

We write below $Q = S^{-1}$, so that

$$ \tilde{Q} = \tilde{S}^{-1} = Q / \mu_w. \quad (B.4) $$

For the 71p dynamo data, we then introduce the following notations:

$$ \begin{gathered} \tilde{\mathbf{X}}^a = [\tilde{\mathbf{x}}_{t_{Nw}}^a \dots \tilde{\mathbf{x}}_{t_{Nw}^a}^a], \quad \Delta \tilde{\mathbf{X}}^a = [\Delta \tilde{\mathbf{x}}_{t_{Nw}}^a \dots \Delta \tilde{\mathbf{x}}_{t_{Nw}^a}^a] \quad \text{and} \\ \tilde{\mathbf{W}}^a = (\Delta \tilde{\mathbf{X}}^a + \Delta t^* \Delta \tilde{\mathbf{X}}^a) / \sqrt{\Delta t^*}. \end{gathered} \quad (B.5) $$

All three matrices are of dimension $P \times (N^a - 2N_w)$. Regarding the 8 realizations of the 100p dynamo, we have

$$ \mathbf{X}^{b,j} = [\mathbf{x}_{t_1}^{b,j} \dots \mathbf{x}_{t_{Nw}^{b,j}}^{b,j}] \quad \text{and} \quad \Delta \mathbf{X}^{b,j} = [\Delta \mathbf{x}_{t_1}^{b,j} \dots \Delta \mathbf{x}_{t_{Nw}^{b,j}}^{b,j}], \quad \text{for } j \in [18]. \quad (B.6) $$

Those are concatenated as

$$ \begin{gathered} \mathbf{X}^b = [\mathbf{X}^{b,1} \dots \mathbf{X}^{b,8}], \quad \Delta \mathbf{X}^b = [\Delta \mathbf{X}^{b,1} \dots \Delta \mathbf{X}^{b,8}], \quad \text{so that} \\ \mathbf{W}^b = (\Delta \mathbf{X}^b + \Delta t^* \Delta \mathbf{X}^b) / \sqrt{\Delta t^*}. \end{gathered} \quad (B.7) $$

We write $N^b = \sum_{j=1}^8 N^{b,j}$ the total number of flow increments available from the 8 realizations of the 100p simulations. Note that one extreme event, for
which the 0.2 yr time resolution does not allow to properly estimate the second time-derivative, has been removed from the 100p series.

## Appendix C. Least squares estimation of the model parameters

In order to have at hand independent time-averaged noise samples $\tilde{w}_{ti}$ from 71p, let first consider epochs equally spaced every $N_w$, i.e. the $\tilde{N}^a = \text{int}(\{N^a - 2N_w\}/(2N_w + 1))$ time-filtered estimates stored into

$$ \tilde{\mathbf{X}}_s^a = [\tilde{\mathbf{x}}_{s,1}^a \dots \tilde{\mathbf{x}}_{s,N^a}^a], \quad \text{with} \quad \tilde{\mathbf{x}}_{s,k}^a = \tilde{\mathbf{x}}_{k(2N_w+1)-N_w}^a, \quad \text{for } k \in [1, \tilde{N}^a]. \qquad \text{(C.1)} $$

Similarly we define $\Delta\tilde{\mathbf{X}}_s^a = [\Delta\tilde{\mathbf{x}}_{s,1}^a \dots \Delta\tilde{\mathbf{x}}_{s,N^a}^a]$ (and)

$$ \tilde{\mathbf{W}}_s^a = (\Delta\tilde{\mathbf{X}}_s^a + \Delta t^* \mathrm{D}\tilde{\mathbf{X}}_s^a) / \sqrt{\Delta t^*}. \qquad \text{(C.2)} $$

The noise samples contained into $\tilde{\mathbf{W}}_s^a$ being uncorrelated, the log-likelihood (under Gaussian assumption) of the flow state defined from both the filtered 71p samples and the instantaneous samples from 100p (denoted logp) can be formulated using (B.4) as

$$ -2\log p = -\tilde{N}^a \log|Q| / \mu_w - N^b \log|Q| + \sum_{i=1}^{N^b} (\Delta\tilde{\mathbf{x}}_{ti}^b + \Delta t^* \mathrm{D}\tilde{\mathbf{x}}_{ti}^b)^T \frac{Q}{\Delta t^*} (\Delta\tilde{\mathbf{x}}_{ti}^b + \Delta t^* \mathrm{D}\tilde{\mathbf{x}}_{ti}^b) + \sum_{k=1}^{\tilde{N}^a} (\Delta\tilde{\mathbf{x}}_{s,k}^a + \Delta t^* \mathrm{D}\tilde{\mathbf{x}}_{s,k}^a)^T \frac{Q}{\mu_w \Delta t^*} (\Delta\tilde{\mathbf{x}}_{s,k}^a + \Delta t^* \mathrm{D}\tilde{\mathbf{x}}_{s,k}^a). \qquad \text{(C.3)} $$

Using textbook formulas (e.g. Lütkepohl, 2005, p. 87), we rewrite (C.3) as

$$ -2\log p = \tilde{N}^a \log(\mu_w) - (\tilde{N}^a + N^b) \log|Q| + \text{Tr} \left[ (\Delta\tilde{\mathbf{X}}^b + \Delta t^* \mathrm{D}\tilde{\mathbf{X}}^b) (\Delta\tilde{\mathbf{X}}^b + \Delta t^* \mathrm{D}\tilde{\mathbf{X}}^b)^T \frac{Q}{\Delta t^*} \right] + \text{Tr} \left[ (\Delta\tilde{\mathbf{X}}_s^a + \Delta t^* \mathrm{D}\tilde{\mathbf{X}}_s^a) (\Delta\tilde{\mathbf{X}}_s^a + \Delta t^* \mathrm{D}\tilde{\mathbf{X}}_s^a)^T \frac{Q}{\mu_w \Delta t^*} \right]. \qquad \text{(C.4)} $$

where Tr(M) means the trace of the matrix M. Deriving (C.4) w.r.t. D and Q, we get

$$ \left\{ \begin{aligned} \partial_D(-2\log p) &= \frac{Q}{\Delta t^*} \left[ \frac{\Delta\tilde{\mathbf{X}}_s^a \tilde{\mathbf{X}}_s^{aT}}{\mu_w} + \frac{\Delta t^* \mathrm{D}\tilde{\mathbf{X}}_s^a \tilde{\mathbf{X}}_s^{aT}}{\mu_w} + \Delta\mathbf{X}^b \mathbf{X}^b T + \Delta t^* \mathrm{D}\mathbf{X}^b \mathbf{X}^b T \right] \\ \partial_Q(-2\log p) &= - (\tilde{N}^a + N^b) Q^{-1} + \frac{\tilde{\mathbf{W}}_s^a \tilde{\mathbf{W}}_s^{aT}}{\mu_w} + \mathbf{W}^b \mathbf{W}^b T \end{aligned} \right. . \qquad \text{(C.5)} $$

Maximizing the log-likelihood w.r.t. D and Q then results in

$$ \left\{ \begin{aligned} D^* &= -\frac{1}{\Delta t^*} \left[ \frac{\Delta\tilde{\mathbf{X}}_s^a \tilde{\mathbf{X}}_s^{aT}}{\mu_w} + \Delta\mathbf{X}^b \mathbf{X}^b T \right] \left[ \frac{\tilde{\mathbf{X}}_s^a \tilde{\mathbf{X}}_s^{aT}}{\mu_w} + \mathbf{X}^b \mathbf{X}^b T \right]^{-1} \\ S^* &= \frac{1}{\tilde{N}^a + N^b} \left[ \frac{\tilde{\mathbf{W}}_s^a \tilde{\mathbf{W}}_s^{aT}}{\mu_w} + \mathbf{W}^b \mathbf{W}^b T \right] \end{aligned} \right. . \qquad \text{(C.6)} $$

In practice, we obtain better conditioned matrices by stacking the information contained in all the $N^a$ (instead of $\tilde{N}^a$) snapshot, considering that we have approximately one independent noise estimate per packets of $N_w$ data, i.e. $\approx \tilde{N}^a$ independent data from the 71p dynamo. This leads to approximate the log-likelihood with

$$ -2\log p = -\tilde{N}^a \log|Q| / \mu_w - N^b \log|Q| + \sum_{i=1}^{N^b} (\Delta\tilde{\mathbf{x}}_{ti}^b + \Delta t^* \mathrm{D}\tilde{\mathbf{x}}_{ti}^b)^T \frac{Q}{\Delta t^*} (\Delta\tilde{\mathbf{x}}_{ti}^b + \Delta t^* \mathrm{D}\tilde{\mathbf{x}}_{ti}^b) + \frac{1}{N_w} \sum_{i=N_w}^{N^a-N_w} (\Delta\tilde{\mathbf{x}}_{ti}^a + \Delta t^* \mathrm{D}\tilde{\mathbf{x}}_{ti}^a)^T \frac{Q}{\mu_w \Delta t^*} (\Delta\tilde{\mathbf{x}}_{ti}^a + \Delta t^* \mathrm{D}\tilde{\mathbf{x}}_{ti}^a). \qquad \text{(C.7)} $$

The optimal estimates for D and S then become
$$
\begin{equation}
\begin{aligned}
D^* &= -\frac{1}{\Delta t^*} \left[ \frac{\Delta \tilde{\mathbf{X}}^a \tilde{\mathbf{X}}^{aT}}{\mu_w N_w} + \Delta \mathbf{X}^b \mathbf{X}^{bT} \right] \left[ \frac{\tilde{\mathbf{X}}^a \tilde{\mathbf{X}}^{aT}}{\mu_w N_w} + \mathbf{X}^b \mathbf{X}^{bT} \right]^{-1} \\
S^* &= \frac{1}{\tilde{N}^a + N^b} \left[ \frac{\tilde{\mathbf{W}}^a \tilde{\mathbf{W}}^{aT}}{\mu_w N_w} + \mathbf{W}^b \mathbf{W}^{bT} \right]
\end{aligned}
\tag{C.8}
\end{equation}
$$

By extracting blocks from $D^*$ and $S^*$ in (C.8), we get the estimated matrices $\tilde{D}_{t \in \{0,1,2\}}^*$ and $\tilde{B}^*$. Then using (A.6) we finally obtain the operators $D_{0,1,2}$ and $B$ required for the time integration of $v_t$ with the forward model (A.4), for any arbitrary forecast time-step $\Delta t$, provided $\Delta t \ll 1/\max(|\lambda|)$, where the vector $\lambda$ contains the eigenvalues of the drift matrix $D$. Typically $1/\max(|\lambda|)$ is close to the Nyquist period $2\Delta t^*$. We show in Fig. C.13 the range of decay times $1/\text{Re}(\lambda)$ and oscillatory periods $2\pi/\text{Re}(\lambda)$ for $D^*$ obtained using (C.8). These span a few years to millennia. The associated oscillators appear significantly damped, with decay times of magnitude comparable to the period.

<figure><img src="image_17.png" /><figcaption>Fig. C.13. Period $(2\pi/\text{Im}(\lambda))$ as a function of decay time $(1/\text{Re}(\lambda))$ for the eigenvalues $\lambda$ of the operator $D^*$ for the AR-3 process. The dashed line corresponds to $y = x$.</figcaption></figure>

## Appendix D. Adaptation of the analysis steps to the AR-3 model

### D.1. Analysis of the magnetic field

We consider below a reanalysis performed every $\Delta t_a$ over the time-span $[t_a^{start}, t_a^{end}]$. We note $H_b$ the observation operator relating MF observations to the field model coefficient (stored in a vector $b_t$) at a given epoch $t$. $H_b$ is the identity matrix in the case of Gauss coefficients data, and would be replaced by data kernels for GVO or ground-based data, as in Istas et al. (2023). We write $\mathbf{y}_t^{MF}$ a MF observation vector at epoch $t$, attached of errors $\mathbf{e}_t^{MF}$ characterized by a cross-covariance matrix $R_t^{MF} = \mathbb{E}(\mathbf{e}_t^{MF} \mathbf{e}_t^{MF^T})$. For all analysis epochs $t_a^{start} < t_a < t_a^{end}$, we relate MF observations at the three nearby analysis epochs to the model state $\mathbf{B}_{t_a}^T = [\mathbf{b}_{t_a}^T, \mathbf{b}_{t_a}^{T}, \mathbf{b}_{t_a}^{T}]$ using a 2nd order Taylor expansion, synthesized as

$$ \mathbf{Y}_{t_a}^{MF} = H_B \mathbf{B}_{t_a} + \mathbf{E}_{t_a}^{MF}, \quad \text{where} \quad \mathbf{Y}_{t_a}^{MF} = \begin{bmatrix} \mathbf{y}_{t_a}^{MF} \\ \mathbf{y}_{t_a}^{MF} \\ \mathbf{y}_{t_a}^{MF} \end{bmatrix}, \quad \mathbf{E}_{t_a}^{MF} = \begin{bmatrix} \mathbf{e}_{t_a}^{MF} \\ \mathbf{e}_{t_a}^{MF} \\ \mathbf{e}_{t_a}^{MF} \end{bmatrix} \quad \text{and} \quad \mathbf{H}_B = \begin{bmatrix} H_b & -\Delta t_a H_b & \frac{1}{2} \Delta t_a^2 H_b \\ H_b & 0 & 0 \\ H_b & +\Delta t_a H_b & \frac{1}{2} \Delta t_a^2 H_b \end{bmatrix}. \tag{D.1} $$

From the ensemble of forecast realizations $\mathbf{b}_t^f$ over the interval $[t_a^-, t_a^+]$, we derive an ensemble of $\mathbf{B}_{t_a}^f = [\mathbf{b}_{t_a}^T, \mathbf{b}_{t_a}^{T}, \mathbf{b}_{t_a}^{T}]^T$ using a cubic polynomial fit to each of the MF forecast trajectories. From this latter we compute the empirical forecast covariance matrix $\mathbf{P}_{BB}^f(t_a) = \mathbb{E}(\delta\mathbf{B}_{t_a}^f \delta\mathbf{B}_{t_a}^f)$, with the notation $\delta\mathbf{B} = \mathbf{B} - \langle \mathbf{B} \rangle$ and $\langle \mathbf{B} \rangle$ the ensemble average. We apply on $\mathbf{P}_{BB}^f$ the Graphical LASSO method (e.g. Friedman et al., 2007), with the free parameter $\lambda = 0.1$ chosen so as to reduce spurious cross-correlations while keeping significant cross-covariances in $\mathbf{P}_{BB}^f$ and maintaining it positive definite (see Istas et al., 2023). From the resulting $\mathbf{P}_{BB}^f$ and $H_B$ we then obtain for each ensemble member an analysed MF model as

$$ \mathbf{B}_{t_a}^a = \mathbf{B}_{t_a}^f + \mathbf{P}_{BB}^f(t_a) \mathbf{H}_B^T [\mathbf{H}_B \mathbf{P}_{BB}^f(t_a) \mathbf{H}_B^T + \mathbf{R}^{MF}(t_a)]^{-1} (\mathbf{Y}_{t_a}^{MF} - \mathbf{H}_B \mathbf{B}_{t_a}^f). \tag{D.2} $$

The MF data error covariance matrix $\mathbf{R}^{MF}(t_a)$ is here made of three diagonal blocks, each composed of the MF observation error matrices $\mathbf{R}_{t_a, t_a, t_a}^{MF}$ at the three successive analysis epochs. The analysed MF at the three epochs can then be deduced from $\mathbf{B}_{t_a}$ using a Taylor expansion:
$$ \begin{bmatrix} \mathbf{b}_{t_a}^a \\ \mathbf{b}_{t_a^c}^a \\ \mathbf{b}_{t_a^c}^a \end{bmatrix} = \begin{bmatrix} I_b & -\Delta t_a I_b & \frac{1}{2} \Delta t_a^2 I_b \\ I_b & 0 & 0 \\ I_b & +\Delta t_a I_b & \frac{1}{2} \Delta t_a^2 I_b \end{bmatrix} \mathbf{B}_{t_a}. \qquad (D.3) $$

$I_b$ is the identity matrix of the size of the vector **b**.

## D.2. Analysis of the augmented state **z**

A similar strategy is considered for the subsequent analysis of **z**. Let write $\mathbf{Z}_{t_a} = [\mathbf{z}_{t_a}^T, \mathbf{z}_{t_a}^{f^T}, \mathbf{z}_{t_a}^{f^T}]^T$. SV observations at the three nearby analysis epochs are related to the model state $\mathbf{Z}_{t_a}$ using a 2nd order Taylor expansion, assuming that we know the MF state at the three epochs of interest from the MF analysis described above. This is expressed as

$$ \Delta \mathbf{Y}_{t_a}^{SV} = H_Z(\mathbf{B}_{t_a}) \mathbf{Z}_{t_a} + \mathbf{E}_{t_a}^{SV}, \text{ where } \Delta \mathbf{Y}_{t_a}^{SV} = \begin{bmatrix} \Delta \mathbf{y}_{t_a}^{SV} \\ \Delta \mathbf{y}_{t_a}^{SV} \\ \Delta \mathbf{y}_{t_a}^{SV} \end{bmatrix}, \mathbf{E}_{t_a}^{SV} = \begin{bmatrix} \mathbf{e}_{t_a}^{SV} \\ \mathbf{e}_{t_a}^{SV} \\ \mathbf{e}_{t_a}^{SV} \end{bmatrix} \text{ and} \\ H_Z(\mathbf{B}_{t_a}) = \begin{bmatrix} H_z(\mathbf{b}_{t_a}^a) & -\Delta t_a H_z(\mathbf{b}_{t_a}^a) & \frac{1}{2} \Delta t_a^2 H_z(\mathbf{b}_{t_a}^a) \\ H_z(\mathbf{b}_{t_a}^a) & 0 & 0 \\ H_z(\mathbf{b}_{t_a}^a) & +\Delta t_a H_z(\mathbf{b}_{t_a}^a) & \frac{1}{2} \Delta t_a^2 H_z(\mathbf{b}_{t_a}^a) \end{bmatrix}. \qquad (D.4) $$

$\Delta \mathbf{y}_{t_a}^{SV}$ is a SV observation vector at epoch $t$, once corrected from the SV prediction from the average dynamo state – see Eq. (4). $\mathbf{e}_{t_a}^{SV}$ is the associated observation error vector, of covariance matrix $R_t^{SV} = \mathbb{E}(\mathbf{e}_{t_a}^{SV} \mathbf{e}_{t_a}^{SV^T})$. The observation operator

$$ H_z(\mathbf{b}_t^a) = [H_b \hat{A}(\mathbf{b}_t^a) H_b] \qquad (D.5) $$

depends on the previously analysed MF, with $\hat{A}$ coding for the elements of the matrix that enters the electro-motive force from the large-scale field and flow, once projected into the PCA world (see Eq. (4)).

From the ensemble of forecast realizations $\mathbf{z}_t^f$ over the interval $[t_a, t_a^+]$, we derive an ensemble of $\mathbf{Z}_{t_a} = [\mathbf{z}_{t_a}^{f^T}, \mathbf{z}_{t_a}^{f^T}, \mathbf{z}_{t_a}^{f^T}]^T$, using a cubic polynomial fit. From the latter we compute the empirical forecast covariance matrix $\mathbf{P}_{ZZ}^f(t_a) = \mathbb{E}(\delta\mathbf{Z}_{t_a}^f \delta\mathbf{Z}_{t_a}^f)$, on which we apply the Graphical LASSO method (with the free parameter $\lambda = 0.1$) in order to reduce the impact of spurious cross-correlations. From $\mathbf{P}_{ZZ}^f$ and $H_Z$ we obtain for each ensemble member the analysed augmented state as

$$ \mathbf{Z}_{t_a}^a = \mathbf{Z}_{t_a}^f + \mathbf{P}_{ZZ}^f(t_a) \mathbf{H}_Z^T [\mathbf{H}_Z \mathbf{P}_{ZZ}^f(t_a) \mathbf{H}_Z^T + \mathbf{R}^{SV}(t_a)]^{-1} (\Delta \mathbf{Y}_{t_a}^{SV} - \mathbf{H}_Z \mathbf{Z}_{t_a}^f). \qquad (D.6) $$

The SV data error covariance matrix $\mathbf{R}^{SV}$ is made of three diagonal blocks, each composed of the SV observation error matrices $\mathbf{R}_{t_a, t_a, t_a}^{SV}$ at the three successive analysis epochs.

## D.3. Specific case of the analysis at the last observation epoch

The above analysis scheme is slightly modified for the last analysis at $t_a = t_a^{end}$. As there is no available observation afterward, one has to adapt from the centered Taylor expansion used in (D.3-D.4). The modification detailed below is also considered for the last analysis epoch when operating hindcast experiments in Section 4.1. When no data is available after $t_a$, we consider a 2nd order Taylor expansion on the left, thus replacing Eq. (D.1) by

$$ H_B = \begin{bmatrix} H_b & -2\Delta t_a H_b & 2\Delta t_a^2 H_b \\ H_b & -\Delta t_a H_b & \frac{1}{2}\Delta t_a^2 H_b \\ H_b & 0 & 0 \end{bmatrix}, \text{ with } \mathbf{Y}_{t_a}^{MF} = \begin{bmatrix} \mathbf{y}_{t_a}^{MF} \\ \mathbf{y}_{t_a}^{MF} \\ \mathbf{y}_{t_a}^{MF} \end{bmatrix} \text{ and } \mathbf{E}_{t_a}^{MF} = \begin{bmatrix} \mathbf{e}_{t_a}^{MF} \\ \mathbf{e}_{t_a}^{MF} \\ \mathbf{e}_{t_a}^{MF} \end{bmatrix}. \qquad (D.7) $$

Here we introduce the notation $t_a^{2-} = t_a - 2\Delta t_a$. The analysed MF at the last three analysis epochs is then obtained as

$$ \begin{bmatrix} \mathbf{b}_{t_a}^a \\ \mathbf{b}_{t_a}^a \\ \mathbf{b}_{t_a}^a \end{bmatrix} = \begin{bmatrix} I_b & -2\Delta t_a I_b & 2\Delta t_a^2 I_b \\ I_b & -\Delta t_a I_b & \frac{1}{2}\Delta t_a^2 I_b \\ I_b & 0 & 0 \end{bmatrix} \mathbf{B}_{t_a}, \qquad (D.8) $$

instead of Eq. (D.3). From it we generate the observation operator for the augmented state analysis as
$$H_Z(\mathbf{B}_{ta}) = \begin{bmatrix} H_z(\mathbf{b}_{ta}^a) & -2\Delta t_a H_z(\mathbf{b}_{ta}^a) & 2\Delta t_a^2 H_z(\mathbf{b}_{ta}^a) \\ H_z(\mathbf{b}_{ta}^{a-}) & -\Delta t_a H_z(\mathbf{b}_{ta}^a) & \frac{1}{2} \Delta t_a^2 H_z(\mathbf{b}_{ta}^a) \\ H_z(\mathbf{b}_{ta}^b) & 0 & 0 \end{bmatrix}, \qquad (D.9)$$

to replace Eq. (D.4).

## Data availability

Data will be made available on request.

## References

Alken, P., Thébault, E., Beggan, C.D., Amit, H., Aubert, J., Baerenzung, J., Bondar, T., Brown, W., Calififf, S., Chambodut, A., et al., 2021. International geomagnetic reference field: the thirteenth generation. Earth Planets Space 73, 1–25.

Aubert, J., 2015. Geomagnetic forecasts driven by thermal wind dynamics in the earth's core. Geophys. J. Int. 203, 1738–1751.

Aubert, J., 2018. Geomagnetic acceleration and rapid hydromagnetic wave dynamics in advanced numerical simulations of the geodynamo. Geophys. J. Int. 214, 531–547.

Aubert, J., 2019. Approaching earth's core conditions in high-resolution geodynamo simulations. Geophys. J. Int. 219, S137–S151.

Aubert, J., 2023. State and evolution of the geodynamo from numerical models reaching the physical conditions of earth's core. Geophys. J. Int. 235, 468–487.

Aubert, J., Finlay, C.C., 2019. Geomagnetic jerks and rapid hydromagnetic waves focusing at Earth's core surface. Nat. Geosci. 12, 393–398.

Aubert, J., Gillet, N., 2021. The interplay of fast waves and slow convection in geodynamo simulations nearing Earth's core conditions. Geophys. J. Int. 225, 1854–1873.

Aubert, J., Finlay, C.C., Fournier, A., 2013. Bottom-up control of geomagnetic secular variation by the earth's inner core. Nature 502, 219–223.

Aubert, J., Gastine, T., Fournier, A., 2017. Spherical convective dynamos in the rapidly rotating asymptotic regime. J. Fluid Mech. 813, 558–593.

Aubert, J., Livermore, P.W., Finlay, C.C., Fournier, A., Gillet, N., 2022. A taxonomy of simulated geomagnetic jerks. Geophys. J. Int. 231, 650–672.

Baerenzung, J., Holschneider, M., Wicht, J., Sanchez, S., Lesur, V., 2018. Modeling and predicting the short term evolution of the geomagnetic field. J. Geophys. Res. Solid Earth 123, 4539–4560.

Baerenzung, J., Holschneider, M., Wicht, J., Lesur, V., Sanchez, S., 2020. The kalmag model as a candidate for igrf-13. Earth Planets Space 72, 1–13.

Baerenzung, J., Holschneider, M., Saynisch-Wagner, J., Maik, T., 2022. Kalmag: a high spatio-temporal model of the geomagnetic field. Earth Planets Space 74, 139.

Barrois, O., Aubert, J., 2024. Characterisation of hydromagnetic waves propagating over a steady, non-axisymmetric background magnetic field. Proc. Roy. Soc. A 480, 20240067.

Barrois, O., Gillet, N., Aubert, J., 2017. Contributions to the geomagnetic secular variation from a reanalysis of core surface dynamics. Geophys. J. Int. 211, 50–68.

Beggan, C.D., Whaler, K.A., 2018. Ensemble kalman filter analysis of magnetic field models during the champ-swarm gap. Phys. Earth Planet. Inter. 281, 103–110.

Blangsball, R.M., Finlay, C.C., Kloss, C., 2022. Localized origin at the core-mantle boundary of the 1969 geomagnetic pulse. Geophys. Res. Lett. 49, e2022GL101070.

Buffett, B., Knezek, N., 2018. Stochastic generation of mac waves and implications for convection in earth's core. Geophys. J. Int. 212, 1523–1535.

Buffett, B., Matsui, H., 2019. Equatorially trapped waves in earth's core. Geophys. J. Int. 218, 1210–1225.

Cain, J., Sweeney, R., 1973. The pogo data. J. Atmos. Terr. Phys. 35, 1231–1247.

Canet, E., Fournier, A., Jault, D., 2009. Forward and adjoint quasi-geostrophic models of the geomagnetic secular variation. J. Geophys. Res. Solid Earth 114.

Chi-Durán, R., Buffett, B.A., 2023. Extracting spatial–temporal coherent patterns in geomagnetic secular variation using dynamic mode decomposition. Geophys. Res. Lett. 50, e2022GL101288.

Christensen, U.R., Aubert, J., Hulot, G., 2010. Conditions for earth-like geodynamo models. Earth Planet. Sci. Lett. 296, 487–496.

Chulliat, A., Maus, S., 2014. Geomagnetic secular acceleration, jerks, and a localized standing wave at the core surface from 2000 to 2010. J. Geophys. Res. Solid Earth 119, 1531–1543.

Cosme, E., Verron, J., Brasseur, P., Blum, J., Auroux, D., 2012. Smoothing problems in a bayesian framework and their linear gaussian solutions. Mon. Weather Rev. 140, 683–695.

Davies, C., Pozzo, M., Gubbins, D., Alfe, D., 2015. Constraints from material properties on the dynamics and evolution of earth's core. Nat. Geosci. 8, 678–685.

Duan, P., Huang, C., 2020. Intradecalad variations in length of day and their correspondence with geomagnetic jerks. Nat. Commun. 11, 1–8.

Evensen, G., 2003. The ensemble Kalman filter: theoretical formulation and practical implementation. Ocean Dyn. 53, 343–367.

Evensen, G., Van Leeuwen, P.J., 2000. An ensemble kalman smoother for nonlinear dynamics. Mon. Weather Rev. 128, 1852–1867.

Finlay, C., Dumberry, M., Chulliat, A., Pais, M., 2010. Short timescale core dynamics: theory and observations. Space Sci. Rev. 155, 177–218.

Finlay, C.C., Aubert, J., Gillet, N., 2016. Gyre-driven decay of the earth's magnetic dipole. Nat. Commun. 7, 10422.

Finlay, C., Kloss, C., Olsen, N., Hammer, M., Toffner-Clausen, L., Grayver, A., Kuvshinov, A., 2020. The CHAOS-7 geomagnetic field model and observed changes in the South Atlantic anomaly. Earth Planets Space 72, 1–31.

Fournier, A., Hulot, G., Jault, D., Kuang, W., Tangborn, A., Gillet, N., Canet, E., Aubert, J., Lhuillier, F., 2010. An introduction to data assimilation and predictability in geomagnetism. Space Sci. Rev. 155, 247–291.

Fournier, A., Aubert, J., Lesur, V., Thébault, E., 2021. Physics-based secular variation candidate models for the igrf. Earth Planets Space 73, 1–18.

Friedman, J., Hastie, T., Tibshirani, R., 2007. Sparse inverse covariance estimation with the graphical lasso. Biostatistics 9, 432–441.

Gastine, T., Aubert, J., Fournier, A., 2020. Dynamo-based limit to the extent of a stable layer atop earth's core. Geophys. J. Int. 222, 1433–1448.

Gerick, F., Livermore, P., 2024. Interannual magneto-Coriolis modes and their sensitivity on the magnetic field within the Earth's core. Proc. R. Soc. A. https://doi.org/10.1098/rspa.2024.0184, 48020240184.

Gerick, F., Jault, D., Noir, J., 2021. Fast quasi-geostrophic magneto-Coriolis modes in the Earth's Core. Geophys. Res. Lett. 48, e2020GL090803.

Gillet, N., Jault, D., Canet, E., Fournier, A., 2010. Fast torsional waves and strong magnetic field within the Earth's core. Nature 465, 74–77.

Gillet, N., Jault, D., Finlay, C., Olsen, N., 2013. Stochastic modeling of the earth's magnetic field: inversion for covariances over the observatory era. Geochem. Geophys. Geosyst. 14, 766–786.

Gillet, N., Barrois, O., Finlay, C.C., 2015a. Stochastic forecasting of the geomagnetic field from the cov-obs. x1 geomagnetic field model, and candidate models for igrf-12. Earth Planets Space 67, 1–14.

Gillet, N., Jault, D., Finlay, C.C., 2015b. Planetary gyre, time-dependent eddies, torsional waves, and equatorial jets at the Earth's core surface. J. Geophys. Res. Solid Earth 120, 3991–4013.

Gillet, N., Jault, D., Canet, E., 2017. Excitation of travelling torsional normal modes in an earth's core model. Geophys. J. Int. 210, 1503–1516.

Gillet, N., Huder, L., Aubert, J., 2019. A reduced stochastic model of core surface dynamics based on geodynamo simulations. Geophys. J. Int. 219, 522–539.

Gillet, N., Gerick, F., Angappan, R., Jault, D., 2022a. A dynamical prospective on interannual geomagnetic field changes. Surv. Geophys. 43, 71–105.

Gillet, N., Gerick, F., Jault, D., Schwaiger, T., Aubert, J., Istas, M., 2022b. Satellite magnetic data reveal interannual waves in Earth's core. PNAS 119, e2115258119.

Hammer, M.D., Finlay, C.C., Olsen, N., 2021. Applications for cryosat-2 satellite magnetic data in studies of Earth's core field variations. Earth Planets Space 73, 1–22.

Huder, L., Gillet, N., Thollard, F., 2019. Pygeodyn 1.1: 0: a python package for geomagnetic data assimilation. Geosci. Model Dev. 12, 3795–3803.

Huder, L., Gillet, N., Finlay, C.C., Hammer, M.D., Tchoungui, H., 2020. COV-OBS.x2: 180 years of geomagnetic field evolution from ground-based and satellite observations. Earth Planets Space 72, 160.

Hulot, G., Sabaka, T.J., Olsen, N., Fournier, A., 2015. The present and future geomagnetic field. In: Schubert, G. (Ed.), Treatise on Geophysics, Second edition 5. Elsevier, Oxford, pp. 33–78. chapter 5.02.

Istas, M., Gillet, N., Finlay, C., Hammer, M., Huder, L., 2023. Transient core surface dynamics from ground and satellite geomagnetic data. Geophys. J. Int. 233, 1890–1915.

Jackson, N., Finlay, C.C., 2015. Geomagnetic secular variation and its applications to the core. In: Schubert, G. (Ed.), Treatise on Geophysics, Second edition 8. Elsevier, pp. 137–184. chapter 5.05.

Jackson, A., Maffei, S., 2020. Plesio-geostrophy for earth's core: I. Basic equations, inertial modes and induction. Proceed. Royal Soc. A 476, 20200513.

Kloss, C., Finlay, C.C., 2019. Time-dependent low-latitude core flow and geomagnetic field acceleration pulses. Geophys. J. Int. 217, 140–168.

Langel, R., Estes, R., 1982. A geomagnetic field spectrum. Geophys. Res. Lett. 9, 250–253.

Lesur, V., Wardinski, I., Baerenzung, J., Holschneider, M., 2018. On the frequency spectra of the core magnetic field gauss coefficients. Phys. Earth Planet. Inter. 276, 145–158.

Lesur, V., Gillet, N., Hammer, M.D., Mandea, M., 2022. Rapid variations of earth's core magnetic field. Surv. Geophys. 43, 41–69.

Li, J., Lin, Y., Zhang, K., 2024. Dynamic mode decomposition of the core surface flow inverted from geomagnetic field models. Geophys. Res. Lett. 51 e2023GL106362.

Lowes, F.J., 1966. Mean-square values on sphere of spherical harmonic vector fields. J. Geophys. Res. 71 (8), 2179–2179.

Luo, J., Marti, P., Jackson, A., 2022. Waves in the Earth's core. II. Magneto-Coriolis modes. Proceed. Royal Soc. A 478, 20220108.
Lütkepohl, H., 2005. *New Introduction to Multiple Time Series Analysis*. Springer Science & Business Media.

Mandea, M., Holme, R., Pais, A., Pinheiro, K., Jackson, A., Verbanac, G., 2010. Geomagnetic jerks: rapid core field variations and core dynamics. *Space Sci. Rev.* 155, 147–175.

Mound, J., Buffett, B., 2003. Interannual oscillations in length of day: implications for the structure of the mantle and core. *J. Geophys. Res. Solid Earth* 108.

Mound, J., Buffett, B., 2005. Mechanisms of core-mantle angular momentum exchange and the observed spectral properties of torsional oscillations. *J. Geophys. Res. Solid Earth* 110.

Nataf, H.C., Schaeffer, N., 2015. Turbulence in the core. In: Schubert, G. (Ed.), *Treatise on Geophysics, Second edition* 8. Elsevier, pp. 161–181. chapter 8.06.

Nilsson, A., Suttie, N., Stoner, J.S., Muscheler, R., 2022. Recurrent ancient geomagnetic field anomalies shed light on future evolution of the South Atlantic anomaly. *Proc. Natl. Acad. Sci.* 119, e2200749119.

Pais, M.A., Jault, D., 2008. *Quasi-geostrophic flows responsible for the secular variation of the Earth's magnetic field*. *Geophys. J. Int.* 173, 421–443. https://doi.org/10.1111/j.1365-246X.2008.03741.x

Pichon, G., Aubert, J., Fournier, A., 2016. Coupled dynamics of earth's geomagnetic westward drift and inner core super-rotation. *Earth Planet. Sci. Lett.* 437, 114–126.

Pick, L., Korte, M., Thomas, Y., Krivova, N., Wu, C.J., 2019. Evolution of large-scale magnetic fields from near-earth space during the last 11 solar cycles. *J. Geophys. Res. Space Physics* 124, 2527–2540.

Ropp, G., Lesur, V., 2023. Mid-latitude and equatorial core surface flow variations derived from observatory and satellite magnetic data. *Geophys. J. Int.* 234, 1191–1204.

Rosat, S., Gillet, N., 2023. Intradecadal variations in length of day: coherence with models of the earth's core dynamics. *Phys. Earth Planet. Inter.* 341, 107053.

Sanchez, S., Wicht, J., Bärenzung, J., Holschneider, M., 2019. Sequential assimilation of geomagnetic observations: perspectives for the reconstruction and prediction of core dynamics. *Geophys. J. Int.* 217, 1434–1450.

Sanchez, S., Wicht, J., Bärenzung, J., 2020. Predictions of the geomagnetic secular variation based on the ensemble sequential assimilation of geomagnetic field models by dynamo simulations. *Earth Planets Space* 72, 1–20.

Schaeffer, N., Jault, D., Nataf, H.C., Fournier, A., 2017. Turbulent geodynamo simulations: a leap towards Earth's core. *Geophys. J. Int.* 211, 1–29.

Schwaiger, T., Gastine, T., Aubert, J., 2019. Force balance in numerical geodynamo simulations: a systematic study. *Geophys. J. Int.* 219, S101–S114.

Schwaiger, T., Gillet, N., Jault, D., Istas, M., Mandea, M., 2024. Wave-like motions and torques in earth's core as inferred from geomagnetic data: a synthetic study. *Phys. Earth Planet. Inter.* 346, 107104.

Teed, R.J., Jones, C.A., Tobias, S.M., 2014. The dynamics and excitation of torsional waves in geodynamo simulations. *Geophys. J. Int.* 196, 724–735.

Van Kampen, N.G., 1992. *Stochastic Processes in Physics and Chemistry*, 1. Elsevier.

Wardinski, I., Saturnino, D., Amit, H., Chambodut, A., Langlais, B., Mandea, M., Thebault, E., 2020. Geomagnetic core field models and secular variation forecasts for the 13th international geomagnetic reference field (igrf-13). *Earth Planets Space* 72, 1–22.

Zhang, K., 2022. A novel geomagnetic satellite constellation: science and applications. *Earth Planet. Phys.* 7, 4–21.