# Effects of geodynamo priors and geomagnetic data on inverted core surface flows

Hannah F Rogers, Nicolas Gillet, Julien Aubert, Paolo Personnettaz, Mioara Mandea

► To cite this version:

Hannah F Rogers, Nicolas Gillet, Julien Aubert, Paolo Personnettaz, Mioara Mandea. Effects of geodynamo priors and geomagnetic data on inverted core surface flows. Physics of the Earth and Planetary Interiors, 2025, 364, pp.107323. 10.1016/j.pepi.2025.107323. hal-05346261

HAL Id: hal-05346261

https://hal.science/hal-05346261v1

Submitted on 4 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Effects of geodynamo priors and geomagnetic data on inverted core surface flows

Hannah F. Rogers<sup>a,b,*</sup>, Nicolas Gillet<sup>a</sup>, Julien Aubert<sup>c</sup>, Paolo Personnettaz<sup>a,b</sup>, Mioara Mandea<sup>b</sup>

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

## ARTICLE INFO

### Keywords:

Core flow inversion

Geomagnetic secular variation

A priori information

Geodynamo simulations

## ABSTRACT

High quality continuous satellite magnetic data gives unprecedented insights into core surface motions but the inversion process remains under-determined, thus requiring additional information. We make use of the pygeodyn inversion tool, in which the core surface dynamics is forecasted with stochastic equations anchored to prior information from a geodynamo free run. There is no clear overview of the relative importance of the different ‘ingredients’ entering in the core flow inversion process, or of the robustness of the inferred flow features. We document the impact of three factors (or ‘ingredients’): (1) choice of geodynamo prior, (2) choice of geomagnetic field model, and (3) choice of maximum spherical harmonic degree for secular variation. We compare three field models (CHAOS-7, COV-OBS.x2 and Kalmag) and three truncation degree for secular variation. The six dynamo priors allow us to analyse the influence of changing the input parameters towards Earth-like values, changing the dynamo set-up (isotropic boundary conditions or not, stratified or not), and the torques at play.

We quantitatively show that changing the prior has a greater impact on the core surface flow than the choice of magnetic field model used as observations. This difference becomes less pronounced in the past due to the data quality decreasing and the field models showing more divergence. Furthermore, we find that the slowly evolving ‘background’ component of flow is more variable among models than the transient component. Transient flow motions appear very coherent for all considered dynamo priors and field models. The recovery from satellite measurements of subdecadal wave dynamics and of a growing eastward motion under the Pacific are robust and appear data driven. Finally, we quantify the efficiency of core flows for generating magnetic variations as a function of frequency. Normalized to their magnitude, rapid changes generate relatively more magnetic signal than slow variations. This is interpreted as the signature of a slow core dynamics approaching a state where induction is mitigated because of dynamical adjustment. As a consequence, a larger part of core surface motions enters the null-space of the kinematic core flow inverse problem towards long time-scales.

## 1. Introduction

The magnetic field generated by motion of liquid iron in the outer core of the Earth (also called the main field, or MF) is the largest contribution to the magnetic field measured on the surface of the Earth. Assuming the mantle to be insulating, the MF vector **B** above the core surface (of radius $c = 3485$ km) is entirely defined by its radial component $B_r$, whose evolution is governed by the radial induction equation at the top of the core:

$$ \frac{\partial B_r}{\partial t} = - \nabla_h \cdot (\mathbf{u}_h B_r) + \frac{\eta}{r} \nabla^2 (r B_r), \quad (1) $$

where $\eta$ is the magnetic diffusivity and $\mathbf{u}_h$ is the horizontal core motion (e.g. Holme, 2015). $\nabla_h$ is the horizontal gradient operator. Over the years, various inversion techniques have been designed for investigating core motions from measurements of the main field.

Traditionally, the high conductivity in the core allows for the magnetic diffusion in Eq. (1) to be considered negligible on periods of years to decades for the largest scales ($\approx 1000$ km) (Roberts and Scott, 1965). This happens to be also the case over a wide range of time-scales in

<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> Centre National d'Etudes Spatiales, 2 place Maurice Quentin, Paris Cedex 01 75039, France</footnote>

<footnote><sup>c</sup> Institut de Physique du Globe de Paris, CNRS, Université Paris Cite, Paris 75005, France</footnote>

<footnote>* Corresponding author at: School of Earth and Environment, University of Leeds, Woodhouse, Leeds LS2 9JT, UK.
E-mail address: h.f.rogers@leeds.ac.uk (H.F. Rogers).</footnote>
advanced numerical simulations of the geodynamo (Aubert, 2018; Aubert and Gillet, 2021). Neglecting diffusion comes down to the ‘frozen-flux’ assumption where observed large scale temporal changes in the magnetic field over the period of years to centuries can be used as a ‘tracer’ for the fluid motion at the core-mantle boundary (CMB) – e.g. Backus and Mouël (1986); Bloxham and Jackson (1991). Roberts and Scott (1965) recognised that flow determination is a non-unique problem, as later quantified by Backus (1968). On top of the frozen-flux hypothesis, Eq. (1) requires additional assumptions about the structure of the flow in order to reduce the ambiguity generated from inverting for two unknown parameters (the northward and eastward flow components) from one observational constraint (the radial magnetic field and its time derivative – known as Secular Variation, SV). For further discussion of the associated issues in the case of helical and quasi-geostrophic flows, we refer respectively to Amit and Olson (2004) and Schwaiger et al. (2023).

Flow models can be constructed under different assumptions (for instance toroidal, steady, tangentially-geostrophic, quasi-geostrophic etc., see Finlay et al., 2010a), which reduces ambiguity but alter the flow structure determined. With the advances in numerical simulations, it has become possible to build this required prior information based on second-order statistics from free runs of the primitive geodynamo equations (e.g. Fournier et al., 2011; Barrois et al., 2017; Sanchez et al., 2020). To build converged statistics, these should be run for long enough so as to cover O(100) vortex turn-over times $\tau_U$ (≈ 130 yr in Earth's core), i.e. some 10 kyr. Generating such simulations at Earth-like conditions still constitutes a challenge, and prior dynamo norms are possibly sensitive to the input parameters in the dynamo set-up (in particular the considered boundary conditions and core-mantle couplings). In this study we attempt to quantify the effect of the three most important factors (or ‘ingredients’), with the pygeodyn geomagnetic data assimilation tool (Huder et al., 2019). We are interested in how the choice of geodynamo prior, the choice of geomagnetic field (and SV) model, and choice of maximum spherical harmonic degree (therefore spatial resolution) for secular variation will affect the resulting flow model.

The paper is structured by first giving a detailed description of the data and methodology used. In Section 2 we provide a summary of the field models we consider in this work, before we describe the geodynamo models used to build the prior information in the inversion process. We then present the pygeodyn inversion methodology in Section 3.1, recalling the particular points where changing the ‘input ingredients’ (i.e. dynamo prior, SV model, and its truncation) could affect the flow model outcome. In Section 3.2 we introduce the diagnostics used to compare the resulting models from the flow inversion. Our results are split into two main types: firstly the effect of the input ingredients on flow estimates at a given epoch (Section 4), then on transient motions (Section 5). Our paper ends with a discussion (Section 6) and a conclusion (Section 7).

## 2. Magnetic data and geodynamo priors

The comparison proposed in this study is based on three magnetic field and secular variation models, and six geodynamo prior models, which are input into the pygeodyn tool for core flow reconstruction. One parameter to choose when utilising secular variation models is the maximum spherical harmonic truncation degree, which we discuss

further in Section 2.1. The data and priors have mostly been published in other studies, which we reference thoroughly throughout. To help the reader, mathematical notations are summarized in Appendix A.

### 2.1. Magnetic field and secular variation models

For investigating the effect of field model on the core surface flow inversion, we use three different magnetic field and SV models: Kalmag (Baerenzung et al., 2020), CHAOS-7 (Finlay et al., 2020) and COV-OBS. x2 (Huder et al., 2020). CHAOS-7 is not constructed prior to the continuous satellite era so we have three models for comparison starting from 1999, and two models over the pre-satellite era. We recall below their main characteristics and differing constructions, also summarized in Table 1. All describe the time evolution of the radial field at the core surface, assuming a potential field outside Earth's core. These models are represented by means of spherical harmonics as:

$$B_r(r, \theta, \phi, t) = \sum_{l=1}^{L} (l+1) \left(\frac{a}{r}\right)^{l+2} \sum_{m=0}^{l} (g_l^m(t) \cos(m\phi) + h_l^m(t) \sin(m\phi)) P_l^m(\cos\theta), \quad (2)$$

where L is the maximum spherical harmonic degree, a = 6371.2 km is the Earth's reference radius, $P_l^m$ the associated Legendre polynomials of degree l and order m, and $(g_l^m, h_l^m)$ the Schmidt semi-normalized Gauss coefficients. SV models are represented in a similar way with SV Gauss coefficients $g_l^m(t)$ and $h_l^m(t)$.

Traditionally, in the expansion of Eq. (2) the maximum spherical harmonic degree for SV ($L_{SV}$) is chosen to be the same as for MF ($L_B$). This latter is fixed to $L_B = 13$ due to the difficulty in separating the main field from the crustal magnetic field (Langel and Estes, 1982). However, core field time-changes are much larger than magnetic field variations associated with the induced magnetization (Thébault et al., 2009; Hulot et al., 2009), and operate on shorter time-scales. Therefore, there is the potential to gain additional information by considering the SV to a higher maximum spherical harmonic degree than the main field, at least when continuous satellite data are available. Here we consider $L_{SV}$ to be 13, 16 and 18 to see the effect of increasing SV spatial resolution with CHAOS-7 and Kalmag (for COV-OBS.x2, proposed up to spherical harmonic degree 14, only $L_{SV} = 13$ is considered).

From this spectral representation we define the main field spatial spectrum at the Earth's surface (Lowes, 1966),

$$S_B(l) = (l+1) \sum_{m=0}^{l} \left( (g_l^m)^2 + (h_l^m)^2 \right), \quad (3)$$

a similar spectrum $S_{SV}(l)$ for the SV, and two norms for the field at the core surface, (e.g. Shure et al., 1982)

$$\begin{cases} \|B_r\|^2 = \sum_{l=1}^{L_B} \frac{l+1}{2l+1} \left(\frac{a}{c}\right)^{2l+4} S_B(l) \\ \|\nabla_h B_r\|^2 = \frac{1}{c^2} \sum_{l=1}^{L_B} \frac{l(l+1)^2}{2l+1} \left(\frac{a}{c}\right)^{2l+4} S_B(l) \end{cases} \quad (4)$$

where $\|X\|^2 = \frac{1}{\zeta} \int_{\Sigma} X^2 d\Sigma$ and $\|X_h\|^2 = \frac{1}{\zeta} \int_{\Sigma} (X_{\theta}^2 + X_{\phi}^2) d\Sigma$, with $\Sigma$ the core surface.

<table><caption>**Table 1**

Summary of the main characteristics considered for the three magnetic field models used in this study.</caption><thead><tr><th></th><th>L<sub>B</sub></th><th>L<sub>SV</sub></th><th>Time-span</th><th>A priori</th><th>Time Modelling</th></tr></thead><tbody><tr><td>CHAOS-7</td><td>13</td><td>13, 16, 18</td><td>1999–2023</td><td>Regularized</td><td>order 6 splines every 6 months</td></tr><tr><td>COV-OBS.x2</td><td>13</td><td>13</td><td>1880–2020</td><td>AR-2</td><td>order 4 splines every 2 years</td></tr><tr><td>Kalmag</td><td>13</td><td>13, 16, 18</td><td>1900–2023</td><td>AR-2</td><td>sequence of snapshots</td></tr></tbody></table>
## 2.1.1. CHAOS-7

The CHAOS suite of models describes the time-dependent near-Earth geomagnetic field, starting from 1999 (Olsen et al., 2006, 2014; Finlay et al., 2016, 2020). We use here the version CHAOS-7.16 that ends in 2023. It is based on measurements from the Oersted, CHAMP, Cryosat-2 and Swarm satellite missions, on top of ground-based magnetic observations. CHAOS-7 consists of a time-dependent internal field expanded up to spherical harmonic degree 20, a static internal field that merges to the LCS-1 lithospheric field model above degree 25 (Olsen et al., 2017), and a model of the magnetospheric field and its induced counterpart. Strict geomagnetically quiet-time criteria are applied to magnetic data and an iteratively reweighted regularized least-squares procedure is applied to calculate the model parameters. Gauss coefficients describing the core field are projected in time using splines of order 6, with knots every 6 months. Their inversion relies on a regularization procedure, which involves for the core field a penalty on $\partial^3 B_r/\partial t^3$ for low spherical harmonic degrees, plus a penalty on $\partial^2 B_r/\partial t^2$ at the start and end points. Note that the damping function of internal field variations varies with the degree and order. In particular a stronger constraint is imposed to zonal coefficients, as they are particularly sensitive to unmodelled magnetic sources towards auroral latitudes. CHAOS-7 is regularly updated using the latest satellite and ground data, and improving regularization sensitivities. Since Gauss coefficients from CHAOS-7 do not come with associated uncertainties, we consider when inverting CHAOS-7 data uncertainties from the Kalmag model (see below).

## 2.1.2. COV-OBS.x2

The COV-OBS.x2 model (Huder et al., 2020) is the latest update to a long-term field model covering the observatory and satellite eras (1840–2020). Note that we choose 1880 as the start of our investigations from this model, in order to mitigate the effect of a degraded spatiotemporal resolution at the early epochs, due to the weaker observational constraint. This timing corresponds to the increasing number of geomagnetic observatories associated with the first International Polar Year 1882–83 (e.g. Hulot et al., 2015). COV-OBS.x2 is produced from historical, surveys and ground-based observations, plus data from the POGO, Magsat, CHAMP and Swarm satellite missions. The latter two are integrated by means of 4-monthly virtual observatories (Hammer et al., 2021). Ground observatories data are integrated as annual differences of annual means prior to 1997, and afterwards annual differences of 4-monthly revised means. The model describes the evolution of the core field (up to degree 14), plus the magnetospheric axial dipole field in dipole coordinates. Gauss coefficients are projected onto splines of order 4, with knots every 2 years. A consequence of this compromise is that variations on periods shorter than 3–4 yr are filtered out (Pick et al., 2019). The a priori constraint on core field parameters needed to reduce non-uniqueness is based on cross-covariance functions associated with second-order auto-regressive (or “AR-2”) processes (Gillet et al., 2013). These are coherent with the occurrence of geomagnetic jerks (discontinuities in the second time-derivative of the field, see Mandea et al., 2010), or a MF power spectral density $S_B(f) \propto f^{-4}$, with f the frequency, on time-scales in the approximate range 2–100 yr (Lesur et al., 2018). Using realistic prior information allows to propose a posteriori uncertainties on the field coefficients.

## 2.1.3. Kalmag

The Kalmag model (Baerenzung et al., 2020, 2022) covers the time-span starting from 1900. It is constrained by ground observatory, satellite (POGO, MagSat, CHAMP, Oersted and Swarm), as well as land, marine and airborne surveys. It co-estimates the internal core and lithospheric fields together with a variety of external sources (near and far magnetospheric fields, field aligned currents...). The core field and secular variation is expanded up to spherical harmonic degree 20. It is constructed using a sequential modelling, with a Kalman filter approach (as initiated by Ropp et al., 2020). This parameterization allows for a

description of shorter times-scales of the core field in comparison with COV-OBS.x2. Based on AR-2 stochastic equations similar to the ones introduced by Gillet et al. (2013), the time integration of an ensemble of internal models allows for the descriptions of uncertainties on the Gauss coefficients. As for all the models, the time resolution decreases for shorter length-scales (see Fig. 5 in Baerenzung et al., 2022). We consider the latest Kalmag update from https://ionocovar.agnld.uni-potsdam.de/Kalmag/ (accessed in December 2023) so that the reanalysis will end in 2023.

## 2.2. Geodynamo prior models

Numerical dynamo modelling is a way of simulating core dynamics from the governing primitive equations that control the physical processes of the evolution of the flow, magnetic field, and co-density in the Earth’s core (Braginsky and Roberts, 1995; Wicht et al., 2010). In general, the parameters of dynamo models are not representative of Earth’s values (e.g. Christensen et al., 2010; Schaeffer et al., 2017). The main dimensionless parameters considered here are the Ekman number $Ek$ (rotation over viscous diffusion times), magnetic Prandtl number $P_m$ (magnetic over viscous diffusion times), magnetic Reynolds number $R_m$ (magnetic diffusion over advection times), Alfvén number $Al$ (Alfvén over advection times), Prandtl number $Pr$ (codensity diffusion time over viscous time), and flux Rayleigh number $Ra_F$:

$$Ek = \frac{\nu}{\Omega D^2}, \quad P_m = \frac{\nu}{\eta}, \quad R_m = \frac{UD}{\eta}, \quad Al = \frac{U}{U_A}, \quad Pr = \frac{\nu}{\kappa}, \quad Ra_F = \frac{g_0 F}{4\pi \rho \Omega^3 D^4} \quad (5)$$

$\Omega$ is the rotation rate of the planet, $D = c - b$ the gap between the outer and inner core boundaries (of radius $b = 1220$ km), $\nu$ the kinematic viscosity, $\eta$ the magnetic diffusivity, $\kappa$ the co-density diffusivity, $g_0$ is the gravitational acceleration at the CMB, $F$ the mass anomaly flux, $U$ and $U_A$ the typical advection and Alfvén velocities ($U_A = B/\sqrt{\rho\mu}$, where $B$ is the typical magnetic field intensity in the bulk of the core, $\rho$ the core density, $\mu$ the magnetic permeability of free-space). Reasonable values of these parameters for the considered simulations and the Earth are given in Table 2 – but keep in mind that we have a limited knowledge of the exact value of some of these dimensionless numbers (Davies et al., 2021). Note that all simulations considered in this study use stress-free mechanical boundary conditions. This is for two reasons. Firstly, this simplifies the determination of the flow at the top of the free-stream, which is in this case the top of the core (Firsov, 2024). Secondly, we need to consider geophysically relevant torques between the core and the mantle to study angular momentum changes, but the viscous torque would be much too high for the large values of Ekman number used in our study (Hollerbach, 1994).

Despite being in the wrong parameter regime, numerical dynamos can produce Earth-like features such as a predominantly dipolar field (Christensen et al., 2010), possibly presenting reversals (e.g. Aubert et al., 2008). They may include an impact of polar wander (Frasson et al., 2024) and mantle control on the geodynamo (Kutzner and Christensen, 2002; Gubbins et al., 2007; Olson et al., 2010), and can be confronted to further geophysical observables such as latitudinal variation of the magnetic dipole direction and strength over paleomagnetic time-scales (Sprain et al., 2019). Criteria for comparing geodynamo simulations to the Earth have been developed based in particular on the field morphology at the CMB, including the fraction of energy stored in the axial dipole component, the degree of equatorial symmetry, the ratio between zonal and non-zonal non-dipole components, and the amount of magnetic flux concentration (Christensen et al., 2010). Diagnostics based on force balances can also be used to assess the Earth-likeness of simulated dynamos (Aubert, 2019; Schwaiger et al., 2020).

Calculating prior information from dynamo simulations requires long series, in order to describe enough independent configurations of
<caption>**Table 2**</caption>

Relevant parameter values for each of the dynamo models and Earth. Earth values have been calculated based on $\nu = 10^{-6}$ m$^2$/s (de Wijs et al., 1998), $\eta = 0.6$ m$^2$/s (Davies et al., 2015), $U = 5 \times 10^{-4}$ m/s (Finlay et al., 2010b), and $B = 4$ mT (Gillet et al., 2010). Time-span is re-scaled to be in years, $N'$ is the number of sampled timesteps, $\Gamma$ is whether the coupling is present (Yes or No, Electromagnetic or Gravitational). The geometry of forcing can be either 'iso' i.e. with homogeneous or isotropic codensity flux at both boundaries or 'het' i.e. with spatially heterogeneous codensity flux boundary conditions at both boundaries. All simulations have $Pr = 1$. If 'LES' is 'Y' then hyperdiffusion has been applied to the velocity and density anomaly fields.

<table><thead><tr><th></th><th>Ek</th><th>P<sub>m</sub></th><th>Ra<sub>F</sub></th><th>R<sub>m</sub></th><th>Al</th><th>Time-span</th><th>N'</th><th>Γ<sub>EM</sub></th><th>Γ<sub>G</sub></th><th>Forcing</th><th>LES?</th></tr></thead><tbody><tr><td>71p</td><td>3 × 10<sup>-10</sup></td><td>7.9 × 10<sup>-3</sup></td><td>2.7 × 10<sup>-10</sup></td><td>1150</td><td>0.049</td><td>10,000</td><td>50,001</td><td>Y</td><td>Y</td><td>het</td><td>Y</td></tr><tr><td>50p</td><td>1 × 10<sup>-8</sup></td><td>4.5 × 10<sup>-2</sup></td><td>8.54 × 10<sup>-9</sup></td><td>1100</td><td>0.114</td><td>20,207.4</td><td>101,037</td><td>Y</td><td>Y</td><td>het</td><td>Y</td></tr><tr><td>0p</td><td>3 × 10<sup>-5</sup></td><td>2.5</td><td>2.7 × 10<sup>-5</sup></td><td>950</td><td>0.7</td><td>135,360</td><td>1,505</td><td>Y</td><td>Y</td><td>het</td><td>N</td></tr><tr><td>Neutral_top1</td><td>3 × 10<sup>-5</sup></td><td>4</td><td>2.7 × 10<sup>-5</sup></td><td>1450</td><td>0.68</td><td>90,530</td><td>9,054</td><td>Y</td><td>Y</td><td>iso</td><td>N</td></tr><tr><td>Stable_top1</td><td>3 × 10<sup>-5</sup></td><td>4</td><td>2.7 × 10<sup>-5</sup></td><td>1500</td><td>0.81</td><td>32,200</td><td>3,221</td><td>Y</td><td>Y</td><td>iso</td><td>N</td></tr><tr><td>S1<sup>γ</sup></td><td>1 × 10<sup>-6</sup></td><td>0.2</td><td>3 × 10<sup>-7</sup></td><td>540</td><td>1.2</td><td>9,418</td><td>3,803</td><td>Y</td><td>N</td><td>iso</td><td>Y</td></tr><tr><td>Earth</td><td>~ 10<sup>-15</sup></td><td>~ 10<sup>-6</sup></td><td>~ 10<sup>-12</sup></td><td>~ 1500</td><td>~ 10<sup>-2</sup></td><td></td><td></td><td>Y</td><td>Y</td><td>het</td><td></td></tr></tbody></table>

the core flow. Considering a turn-over time in the core of the order of 100 yr and basing our statistics on ~100 independent states requires a series spanning O(10) kyr. We restrict our dynamo prior to use maximum spherical harmonic degree for fluid motions modelled in the core flow inversion ($L_U = 18$), of the MF and the SV extracted from the magnetic model ($L_B = 13$ and $L_{SV} = 13, 16$, or 18). The truncated core surface flow is represented via toroidal and poloidal scalars $T$ and $S$ as in Holme (2015).

$$ \mathbf{u}_h = \nabla \wedge (\mathbf{T}\mathbf{r}) + \nabla_H(\mathbf{r}\mathbf{S}), \quad (6) $$

expanded with spherical harmonics as

$$ \left\{ \begin{aligned} T(\theta, \phi, t) &= \sum_{l=1}^{L_U} \sum_{m=0}^{l} [c_l^m(t) \cos(m\phi) + s_l^m(t) \sin(m\phi)] P_l^m(\cos\theta) \\ S(\theta, \phi, t) &= \sum_{l=1}^{L_U} \sum_{m=0}^{l} [c_l^m(t) \cos(m\phi) + s_l^m(t) \sin(m\phi)] P_l^m(\cos\theta). \end{aligned} \right. \quad (7) $$

where $c_l^m$ and $s_l^m$ are the Schmidt semi-normalized toroidal and poloidal flow coefficients in km/yr. From those we define the spatial power spectrum of the core surface kinetic energy

$$ S_U(n) = \frac{l(l+1)}{2l+1} \sum_{m=0}^{l} [c_l^{m2} + s_l^{m2} + c_s^{m2} + s_s^{m2}] \quad (8) $$

and the two norms (Shure et al., 1982; Holme, 2015)

$$ \|\mathbf{u}_h\|^2 = \sum_{l=1}^{L_U} S_U(l) \quad \text{and} \quad \|\nabla_h \cdot \mathbf{u}_h\|^2 = \frac{1}{c^2} \sum_{l=1}^{L_U} \frac{l^2(l+1)^2}{2l+1} \sum_{m=0}^{l} (c_l^{m2} + s_l^{m2}). \quad (9) $$

Numerical dynamos are run in non-Earth-like parameters and need to be re-scaled to allow a direct comparison of the magnetic field spatial spectra from the simulations to geomagnetic field models of the modern era (Lhuillier et al., 2011; Aubert and Gillet, 2021). This allows for a direct incorporation of the simulated Gauss coefficients when calculating the prior information needed for the *pygeodyn* data assimilation tool. There are four main timescales we can use for this: (i) magnetic diffusion timescale ($\tau_\eta = D^2/\eta$), (ii) overturn timescale ($\tau_U = D/U$), (iii) Alfvén timescale ($\tau_A = \sqrt{\rho\bar{\mu}D}/B$), and (iv) rotational timescale ($\tau_\Omega = 2\pi/\Omega$). We scale all dynamos but $S1'$ using the invariants of the parameter space path as $\tau_U/\tau_\eta = 1/R_m$ is preserved along the path (Aubert et al., 2017). The velocity field ($\mathbf{u}$) and the Alfvén velocity ($\mathbf{B}/\sqrt{\rho\bar{\mu}}$) are expressed in units of $\eta/D$ and then rescaled by Earth-like values ($B/\sqrt{\rho\bar{\mu}\Omega\eta} \approx 4.5$). Because $R_m$ for $S1'$ is significantly less than the value expected for the Earth, the rescaling of $S1'$ from numerical to geophysical units has been performed so that a fit of $\tau_{SV}(l) = \sqrt{S_B/S_{SV}}$ (see Eqs. (3) and (8)), of the form $3\tau_U/l$ for $2 \le l \le 13$ (see Bouligand et al., 2016), gives a value of $\tau_U$ comparable to the observed one ($\tau_U = 135$ yr). The Gauss coefficients describing the $S1'$ magnetic field above the core surface are scaled so as to reproduce the spatial power

spectrum $S_B(l)$ found from the CHAOS-7 field model for $2 \le l \le 13$.

Below we describe the two families of models we are considering, along the dynamo path designed by Aubert et al. (2017), hereafter referred to as the AGF17 path, and across this dynamo path. Many dynamo paths exist so the distinction between the two families will show the difference between (i) decreasing Ek and $P_m$ number while keeping a similar set-up (boundary conditions and forcing), and (ii) changing the set-up while keeping comparable values of the Ekman number. Each prior has been run over a time-span long enough to cover at least about 100 turn-over times. The overall periods covered by each considered dynamo are summarized in Table 2, together with the number of snapshots considered to calculate the first- and second-order statistics (mean and cross-covariances) and the key parameters for dynamo generation. Dynamos run at less extreme input parameters (i.e. for larger values of Ek and $P_m$) tend to have longer timespans due to lower computation requirements.

### 2.2.1. Dynamos along the path designed by Aubert et al. (2017)

Aubert et al. (2017) proposed a numerical path towards the geodynamo conditions, suggesting that the MAC (Magneto-Archimedes-Coriolis) force balance in the Earth's core can be used to obtain some dimensionless power-driven viscosity-free scalings. Dynamos along this path use mechanically stress-free, electromagnetically conducting and fixed codensity flux boundary conditions at both boundaries. The fluid shell is electromagnetically coupled to the inner core and to an outer shell representing the mantle. The base of the mantle is made electrically conducting with conductance $\sigma_m\Delta_m$, and a conductance ratio $\Sigma = (\sigma_m\Delta_m)/(\sigma_cD) = 10^{-4}$, where $\sigma_c$ is the electrical conductivity of the fluid core. Both the inner core and mantle feature a time-dependent axial differential rotation with respect to the outer core, and are coupled together through a gravitational restoring torque, following the formulation of Aubert et al. (2013) which involves the differences of angular rotation rate between the inner core and mantle. The surface-averaged mass anomaly flux $F$ emerging from the inner boundary serves to define the flux Rayleigh number $Ra_F$, while $F$ vanishes at the outer boundary. At the inner-core boundary, a longitudinal hemispheric perturbation is superimposed to the homogeneous mass-anomaly flux. At the outer boundary, an heterogeneous perturbation reflecting a thermal control from the deep mantle is also superimposed to the vanishing mass anomaly flux. The magnitudes and patterns of both perturbations are fully described in Aubert et al. (2013, 2017). A volumetric codensity source term is present, to ensure mass conservation. This setup generates, on average, an eccentric gyre that is entirely bottom-driven (see Pichon et al., 2016). In our study, we consider three dynamo simulations along this path (here from the higher to lower values of Ek and $P_m$): the $0p$ dynamo, also known as 'Coupled-Earth' (Aubert et al., 2013), the $50p$ dynamo (Aubert et al., 2017), and the $71p$ dynamo (Aubert and Gillet, 2021). Their main parameters are summarized in Table 2.

By advancing along the AGF17 dynamo path, the value of the mag-
netic Reynolds number $R_m$ is almost preserved, setting a similar relative importance of advection w.r.t. diffusion of the magnetic field for all these simulations. Meanwhile, the various time-scales characteristic of the core dynamics tend to be more and more separated when lowering the input parameters *Ek* and *Pm*. This is measured by an Alfvén number *Al* close to unity for *Op* (i.e. order unity magnetic to kinetic energy ratio), while *Al*≪1 is obtained for *71p* (i.e. a dominant magnetic energy). Reaching such low values of *Al* is rendered possible thanks to the use of hyperdiffusion for the small length-scales of the velocity and codensity fields (large eddy simulation, or “LES”). To rescale dimensionless simulated fields into geophysical units for the AGF17 path models, following Aubert and Gillet (2021) we use $D = 2260$ km as core thickness, the magnetic diffusion time $\tau_\eta = D^2/\eta = 135000$ years as time scale and the Elsasser unit $\sqrt{\rho\mu\eta\Omega} = 0.9$ mT as magnetic field unit.

### 2.2.2. Dynamos across the path designed by Aubert et al. (2017)

The path described above is conditioned to a starting position in the parameter space, and by the specific parameterization of the boundary conditions, torques and forcings. Various such paths may be followed and in this study we also consider prior information built from simulations set on positions across the AGF17 path. As an alternative we use three simulations, run with Ekman numbers similar to the ‘Op’ geodynamo.

Our first variation across the AGF17 path is to consider a dynamo similar to *Op* but run with homogeneous codensity flux boundary conditions (i.e. stripping out the heterogeneous parts imposed at the inner and outer boundaries of the AGF17 models), while the conductance ratio $\Sigma$, EM and gravitational torques coupling the various spherical enve- lopes being kept the same as above. It is named *Neutral_top1*. Note that $P_m$ has been increased (from 2.5 to 4) in order to obtain a larger value of $R_m$ (from $\simeq 950$ to $\simeq 1450$) consistent with estimates of core electrical conductivity (Davies et al., 2015). Secondly, the *Stable_top1* dynamo has been derived from the *Neutral_top1* dynamo, where the only change consists in modifications to the codensity profile so as to generate a stable layer at the top of the fluid core (following Eq. (9) of Gastine et al., 2020). The resulting fixed stratified layer has a thickness of 140 km and contains a stabilising density anomaly corresponding to a Brünt-Väisälä frequency $N = \sqrt{\frac{g_0}{\rho} \frac{d\rho'}{dr}} \simeq 0.89\Omega$ at $r = c$ ($g_0$ the gravity acceleration at the CMB, $\rho'$ the profile of density anomaly w.r.t. the adiabat). This corresponds to a relatively mild stratification in comparison with the configuration investigated for instance by Buffett et al. (2016). To redimensionalize the output of simulations *Neutral_top1* and *Stable_top1*, we use $D = 2260$ km as core thickness, $\tau_\eta = D^2/\eta = 183852$ yr as magnetic diffusion time, and the Elsasser magnetic field unit $\sqrt{\rho\mu\eta\Omega} = 0.84$ mT. This accounts for the difference of magnetic Prandtl number

between these models and the AGF17 models (see Table 2).

Finally, we also consider the $S1^\gamma$ simulation which uses a setup derived from the $S1$ simulation of Schaeffer et al. (2017). $S1^\gamma$ uses fixed-codensity flux boundary condition. The mass anomaly fluxes *F* crossing the inner and outer boundaries are homogeneous, and define the flux-based Rayleigh number $Ra_F$ as above. No volumetric codensity source is present. The inner core boundary is considered electrically insulating. Different from $S1$, $S1^\gamma$ uses stress-free boundary conditions, and includes a solid and electrically conducting layer at the bottom of the mantle of thickness $\Delta_m = 100$km and conductance ratio $\Sigma = (\sigma_m\Delta_m)/(\sigma_c D) = 0.018$. Due to the presence of this layer an electromagnetic (EM) torque is exerted by the conducting mantle on the outer core. In this simulation, it is the only toque to balance the time variation of the angular momentum. Due to the stress free boundary conditions, no viscous torque is present. Furthermore, the insulating inner core does not sustain any EM torque at ICB. The different set-up compared to the original $S1$ does not affect so much the output value of $R_m$ ($\simeq 540$). Note that we did not consider the alternative stress-free $S1^\gamma$ model (run with the same value of *Ek* but with $P_m = 0.8$, leading to $R_m \simeq 1800$, and no core-mantle coupling, see Schwaiger et al., 2023), as for this latter the angular momentum torque was imposed to be constant, leading to very specific prior configuration in particular for the $c^1t_1^0$ and $c^1t_3^0$ coefficients responsible for the length-of-day (LOD) changes. The computation of $S1^\gamma$ relies on hyperdiffusion of the codensity and velocity fields for harmonic degrees larger than 63, using the formulation proposed in Nataf and Schaeffer (2015), with a multiplication factor of 30.

### 2.2.3. Comparison of the dynamo priors

Before entering the inversion process with *pygeodyn*, which involves the magnetic models as well as the above dynamo simulations, we present a brief comparison of the several priors. In Fig. 1 we show the spatial power spectra given by Eq. (8) for the temporal mean $\langle u \rangle$ and for the deviation to the mean $\Delta u = u - \langle u \rangle$ for the dynamo flows. By far the largest time average coefficient is $t_1^0$ for all dynamos gravitationally coupled to the mantle (all but $S1^\gamma$) and thus presenting a large stationary westward drift (Pichon et al., 2016). This is clearly visible on the median flow map for the *Neutral_top1* prior shown in Fig. 2 (b). This westward flow is responsible for the large $S_U(1)$ values. The odd-even oscillation present in spectra for the time-average flow is reminiscent of a predominantly equatorially symmetric, and mostly zonal flow. This feature is particularly pronounced for $S1^\gamma$, related to the presence of intense polar vortices shown in Fig. 2(a). Such features were already documented from the $S1$ dynamo by Schaeffer et al. (2017). They are related with a larger contribution to the surface kinetic energy for the odd degrees zonal motions between $l = 5$ and $l = 13$. The median flow for the $S1^\gamma$ prior is strongly westward within the cylinder tangent to the inner

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Spatial power spectra (Eq. (8)) for the time average flow (a) and the departure from this average (b) of each prior considered in this study.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Comparison of the median flow of the $S1^\gamma$ (a) and *Neutral_top1* (b) priors. (c): profiles of the median azimuthal velocity as a function of latitude for those two priors, superimposed with profiles for the first and last 10 % distribution (dotted lines). Please note that the y-axis limits for (c) are different for the two priors. Continents are shown for reference and can be rotated to any longitude due to the isotropic forcing.</figcaption></figure>

core and aligned with the rotation axis (approximately 20° from the poles), with a median speed up to ≈ 150 km/yr at 80° latitudes, as illustrated with the zonal flow profiles in Fig. 2(c). Note that remaining North-South and longitudinal asymmetries in the median values for those prior with an isotropic forcing are related to the need of longer runs to properly converge the statistics (see also Appendix B). Running the computations for ~100$\tau_U$ may not be enough as geodynamo simulations show natural fluctuations at multi-millenial time-scales (e.g. Schaeffer et al., 2017). We also witness eastward zonal motion at the rim of the tangent cylinder. The signature of westward polar vortices is also visible for *Neutral_top1*, but with a much weaker amplitude (typically some 10 km/yr). For this dynamo the zonal motion, generally westward, is in average close to -5 km/yr near the equator. For both dynamos the zonal flow presents a large time variability, with some eastward equatorial flow sometimes present in *Neutral_top1*, and alternating zonal bands occurring in $S1^\gamma$ between the tangent cylinder and the equator (e.g. light red background at 60° north and south on the orthographic plots in Fig. 2(a) indicating slow eastward motions). We also note that more flow crosses the tangent cylinder within the $S1^\gamma$ dynamo than the *Neutral_top1*, in which the flow moves parallel to the tangent cylinder.

In Fig. 1(b) we show the spectra for the departure to the temporal mean. Except for degrees 1 and 2 (and for the odd degree zonal coefficients for $S1^\gamma$ discussed above), the spatial spectrum for the deviation to the mean is larger than that of the mean for most coefficients and dynamo priors, denoting a large variability of the flow solutions around the median state. Towards short length-scales, for all priors the spectra for $\Delta u$ are rather flat. They are comparable from one dynamo to the other when considering the several priors with gravitational torques. It is slightly weaker for dynamos along the AGF17 path, partly due to a scaling of time based on various magnetic diffusion times. The $S1^\gamma$ prior has more energy across all spherical harmonic degrees than all other

considered dynamos. This is probably due to the difference in the chosen rescaling procedure between $S1^\gamma$ and other dynamo priors (see above). Other factors might also be at the origin of the larger velocities, such as the absence of torque transfer at the ICB in $S1^\gamma$, due to the mechanical and electrical boundary conditions (no gravitational or viscous torque and an insulating inner core). At low degrees, the spectra $S_U$ for $\Delta u$ present a wider range within the several dynamos, the stratified prior *Stable_top1* presenting significantly lower variance for degree 1. We also notice that dynamos along the AGF17 path do not show a marked peak at degree 2, contrary to *Neutral_top1* and *Stable_top1* that do not consider an anisotropic forcing.

The probability distribution function (pdf) of some zonal and non-zonal flow coefficients for the various considered priors are further described in Appendix B. In particular we illustrate a relatively weak change of the distributions while moving along the AGF17 path (a bit more spread while decreasing *Ek* and $P_m$). The use of isotropic boundary forcing for *Neutral_top1* or *Stable_top1* (as opposed to mantle control generating preferential longitudes as in the AGF17 set-up) tends to center the pdf of non-zonal coefficients. The spread of the coefficients pdf is slightly reduced for the stably stratified prior (see Fig. B.1). In the case of $S1^\gamma$ the spread is enhanced, related to the choice for the unit of time when dimensionalizing dynamo outputs into geophysical dimensions.

We also document in Fig. B.2 the power spectra and pdf from the errors of representativeness arising from unmodelled processes (see Section 3.1). These are comparable for one dynamo to the other, although slightly larger on large-scales for *Neutral_top1* and $S1^\gamma$, slightly weaker for *Stable_top1* (but in all cases by a factor less than 2 in r.m.s.). It indicates that any differences seen when changing the dynamo prior will be mostly due to the difference in the flow coefficients priors, and not the error of representativeness associated with them.
# 3. Inversion methodology and diagnostics

## 3.1. Geomagnetic data assimilation with pygeodyn

The pygeodyn inversion methodology has been covered in previous papers (e.g. Huder et al., 2019; Gillet et al., 2019; Istas et al., 2023). Here we only recall the main steps of the employed algorithm, introducing the required notations. We also highlight the locations where our choice of prior dynamo and of SV observation (including the maximum spherical harmonic degree) are influencing the resulting core surface flow output. The pygeodyn inversion tool is based on a sequential data assimilation algorithm (an ensemble Kalman filter, see Evensen, 2003), i.e. constructed from a succession of forecasts for the model state and analysis steps. Analyses are conducted at regular intervals of a multiple of the forecast timestep $\Delta t^f$.

Our forward dynamical model describes the time-evolution of the large scale potential magnetic field $\vec{B}_r$ at $r = c$, by projecting onto large length-scales (as denoted by overlines) the radial induction Eq. (1) to get

$$ \frac{\partial \vec{B}_r}{\partial t} = - \nabla_h' (\mathbf{u}_h \cdot \vec{B}_r) + e_r, \quad (10) $$

where $e_r$ stands for error of representativeness that contains the contributions from the subgrid induction and diffusion processes (see Gillet et al., 2019). All fields are expanded in the spectral domain using spherical harmonics, and Eq. (10) in matrix form becomes

$$ \dot{\mathbf{b}} = A(\mathbf{b})\mathbf{u} + \mathbf{e}, \quad (11) $$

where $\dot{\mathbf{b}}, \mathbf{b}, \mathbf{u}$ and $\mathbf{e}$ are vectors that store the spherical harmonic coefficients for the radial SV, the main field, the core surface flow and the errors of representativeness respectively. The Gaunt-Elsasser matrix $A(\mathbf{b})$ relates the SV to the flow (Whaler, 1986). Below we shall invert simultaneously for the two unknowns $\mathbf{u}$ and $\mathbf{e}$, in an augmented state approach. The MF is truncated to degree $L_B = 13$, the limit over which lithospheric magnetic sources dominate over the core signal (Langel and Estes, 1982). The SV is truncated to different spherical harmonic degrees depending on the model, while $e_r$ is by definition truncated to the same degree as the SV. The flow model is truncated at degree $L_U = 18$, following Gillet et al. (2019). This choice, somewhat arbitrary, is chosen as a trade-off between a risk of aliasing (for values of $L_U$ too close to $L_B$) and the wish to limit the dimension of the unknown vector (flow coefficients of degrees much larger than $L_B$ being unresolved, see Baerenzung et al., 2016; Barrois et al., 2017).

Eq. (10) is coupled to two order 1 auto-regressive (or AR-1) multivariate equations for the time evolution of $\mathbf{e}$ and $\mathbf{u}$. Using an Euler-Maruyama scheme, their time integration takes the form

$$ \left\{ \begin{aligned} \mathbf{u}(t_{k+1}) &= \mathbf{u}(t_k) - \Delta t^f D_u(\mathbf{u}(t_k) - \langle \mathbf{u} \rangle) + \sqrt{\Delta t^f} B_u \mathbf{w}_u(t_k), \\ \mathbf{e}(t_{k+1}) &= \mathbf{e}(t_k) - \Delta t^f D_e(\mathbf{e}(t_k) - \langle \mathbf{e} \rangle) + \sqrt{\Delta t^f} B_e \mathbf{w}_e(t_k) \end{aligned} \right. , \quad (12) $$

where $D_{u,e}$ and $B_{u,e}$ are the drift and diffusion matrices. Vectors $\mathbf{w}_{u,e}$ are built from centered unit variance normal random draws. We want to compare the use of all priors at hand with the same algorithm. For this reason we are unable to use the dense methodology of pygeodyn proposed by Gillet et al. (2019) due to not all priors being run to low enough values of $E_k$ and $P_m$. Simulations such as $Op$ or $S1^r$, run at moderate $E_k$ values, once scaled to geophysical units, indeed allow to model decadal and longer changes of the flow, but present too large Alfvén numbers to replicate shorter time-scales (see Fig. 1 in Gillet et al., 2024), preventing the use of the dense algorithm. We thus employ the diagonal algorithm from pygeodyn (Huder et al., 2019) where the drift matrices simply take the form $D_{u,e} = I_{u,e}/\tau_{u,e}^{AR}$, with the parameters $\tau_{u,e}^{AR} = 30$ yr and $\tau_{e}^{AR} = 10$ yr and $I_{u,e}$ are identity matrices of the same dimension as $\mathbf{u}$ and $\mathbf{e}$. This means all dynamo timesteps are considered independently when constructing the prior. Diffusion matrices $B_{u,e}$ are chosen such that a free run

of Eq. (12) generates vectors $\mathbf{u}$ and $\mathbf{e}$ whose spatial cross-covariances converge towards those extracted from the considered dynamo prior (denoted by stars):

$$ \left\{ \begin{aligned} P_{uu} &= \mathbb{E}(\Delta \mathbf{u} \Delta \mathbf{u}^T) = \frac{1}{N^* - 1} \sum_{j=1}^{N^*} \Delta \mathbf{u}^*(t_j) \Delta \mathbf{u}^*(t_j)^T \\ P_{ee} &= \mathbb{E}(\Delta \mathbf{e} \Delta \mathbf{e}^T) = \frac{1}{N^* - 1} \sum_{j=1}^{N^*} \Delta \mathbf{e}^*(t_j) \Delta \mathbf{e}^*(t_j)^T \end{aligned} \right. , \quad (13) $$

Here the $t_j$ are the epochs where dynamo samples are provided, and $N^*$ is the number of available samples from the dynamo prior (see Table 2). To obtain the second-order statistics $P_{ee}$ from the geodynamo simulations, we calculate $\mathbf{e}$ at each epoch as the difference between the simulated SV and the induction term from the truncated field and flow (it thus contains by construction all unmodelled processes: subgrid induction and diffusion). The dimension of $P_{ee}$ is governed by the truncation degree $L_{SV}$. It is important to note that we rely only on the geodynamo prior information to complete the forecasting steps.

The forecast is operated for an ensemble of $N^e = 200$ realisations of the core state, using $\Delta t^f = 1$ month ($\sim 0.0833$ yr). Each ensemble member is initialized as a random draw from the dynamo statistics. Analyses are performed every $\Delta t^a = 0.5$ yr (6 forecasting steps) over the continuous satellite era (starting 1 yr before 2000), in two steps. Firstly an analysis of the ensemble of realisations of $\mathbf{b}$ is performed from magnetic field observations (Gauss coefficients from a field model, up to degree $L_B = 13$). Next an analysis of the ensemble of realisations of the augmented state $\mathbf{z}$ is performed, where $\mathbf{z}^T = [\mathbf{u}^T \mathbf{e}^T]^T$, using the Gauss coefficient observations of the SV. These are considered for various SV truncation levels $L_{SV}$, depending on what is accessible for each field model.

The analysis step relies on the estimate of several matrices. We only give here the specific details relevant to our present applications. The observation error covariance matrix, which enters the Kalman gain matrix in our sequential assimilation, is built from the dispersion within the ensemble of main field and SV models, when provided, as it is the case for Kalmag and COV-OBS.x2. It is considered as diagonal, and vary with time. For CHAOS-7 data, for which no data error matrix are given, we consider uncertainties from Kalmag. At each analysis epoch, the forecast covariance matrix $P_{zz}^f = \mathbb{E}(\Delta \mathbf{z}^f \Delta \mathbf{z}^f)^T$ is obtained by applying the Graphical-LASSO (Friedman et al., 2008) to the empirical estimate $P_{zz}^f(t^a) = \frac{1}{N^e-1} \sum_{k=1}^{N^e} \Delta \mathbf{z}_k^f(t^a) \Delta \mathbf{z}_k^f(t^a)^T$, where $k$ is the index within the ensemble, the upper-script $f$ stands for the forecast state, and $\Delta \mathbf{z}$ means here the distance to the ensemble mean. As prescribed by Istas et al. (2023), for $N^e = 200$ we consider $\lambda = 0.15$ for the free parameter entering the G-LASSO estimate. In comparison with previous pygeodyn applications, the SV observation operator, which involves the matrix $A(\mathbf{b})$, is expanded to higher degrees when considering $L_{SV} > L_B$.

## 3.2. Diagnostics

In order to compare the influence of the choice of prior dynamo on the flow solution, we estimate the normalized core surface kinetic energy of the difference to a reference model,

$$ \mathcal{Q}(\%) = \frac{\| \mathbf{u} - \mathbf{u}_{ref} \|^2}{\| \mathbf{u}_{ref} \|^2} \times 100. \quad (14) $$

We choose the model calculated from the Neutral top1 prior and using as observations Kalmag up to $L_{SV} = 13$ as a reference model. This choice is motivated by the fact that (i) Neutral top1 is the one off-path dynamo the closest to the path set-up, (ii) Kalmag is the only model that is updated to the most recent epoch (for comparison with the full time-span covered by CHAOS-7) and extends to the pre-satellite era (for comparison with COV-OBS.x2), and (iii) the $L_{SV}$ is often chosen to be the
same as $L_B = 13$.

We furthermore introduce, for a snapshot estimate of $\mathbf{u}_h$ and a given field model $B_r$ at any epoch $t$, a normalized measure of the amount of SV induced by the field and flow at the core surface,

$$ \chi(t) = \frac{\|\partial_t B_r\|}{\|\mathbf{u}_h\| \|\nabla_h B_r\| + \|\tilde{B}_r\| \|\nabla_h \cdot \mathbf{u}_h\| + \|e_r\|} \quad (15) $$

calculated from the norms defined in Eqs. (4) and (9). The errors of representativeness $e_r$ are treated in the same way as the SV in terms of truncation and calculation of the norms. Here the norms are calculated by truncating the fields at degrees $L_B^\dagger = 13$, $L_{SV}^\dagger = 8$ and various values of $L_U^\dagger = 8, 13$ and $18$, in order to mitigate the impact of unresolved flow and SV towards short length-scales and short time-scales. The choice of a reduced $L_{SV}^\dagger$ in the definition of $\chi$ is motivated by limited spatial resolution of the SV further back in time (and by the limited access to SV changes on short time-scales, see below the diagnostic $\tilde{\chi}(f)$). The denominator gives a measure of the potential amount of SV that can be produced for given magnitudes of the flow, the magnetic field and its horizontal gradient, in absence of cancellation (from both terms in $\nabla_h \cdot (\mathbf{u}_h B_r) = \mathbf{u}_h \cdot \nabla_h B_r + B_r \nabla_h \cdot \mathbf{u}_h$). A value of $\chi$ close to zero indicates that the flow, for a given configuration of $B_r$ at the CMB, is very inefficient in inducing SV, i.e., in the null space of the operator $A(\mathbf{b})$. Conversely, high values of $\chi$ suggest an efficient induction for a given flow and field. This measure is an alternative to the angle $\gamma$ between the flow and the horizontal gradient of $B_r$, considered for instance by Finlay and Amit (2011) under the toroidal flow assumption. A value of $\cos\gamma$ near 1 maximizes induction, under the assumption of a negligible stretching term in (15). Conversely, a flow following isocontours of $B_r$ leads to $\cos\gamma = 0$. Finlay and Amit (2011) found weak (resp. high) values of $\cos\gamma$ for high (resp. low) latitude flux patches.

The term $A(\mathbf{b})\mathbf{u}$ is the projection along magnetic observations (the radial SV at the core surface) of the electromotive term $\nabla \times (\mathbf{u} \times \mathbf{B})$ present everywhere in the core. We thus consider the null space of $A(\mathbf{b})$ (documented starting from Roberts and Scott, 1965; Backus, 1968), as the kinematic signature in the observation world of a general tendency for the magnetic and velocity fields to self-organize so as to minimize induction (for subdominant diffusion), or $|\partial_t \mathbf{B}| = |\nabla \times (\mathbf{u} \times \mathbf{B})| \ll \mathbf{UB}/\ell$, where $\mathbf{U}, \mathbf{B}$ and $\ell$ are characteristic values of the flow, the field, and the length-scale. This situation, a consequence of a dynamical equilibration, was observed for the steady flow in geodynamo simulations (Aubert, 2005). It somewhat generalizes Ferraro's (1937) law of isorotation to all three directions of the flow (one way of mitigating induction being the dynamic alignment documented in MHD turbulence, e.g. Tobias et al., 2012; Yokoi (2013)).

An equivalent diagnostic may also be estimated for the dynamics within a band of periods around an epoch of interest ($T_k \in [T_{inf}, T_{up}]$), in which case we separate between the band-pass filtered flow and field perturbations ($\tilde{\mathbf{u}}_h(t)$ and $\tilde{B}_r(t)$) and the quasi-static flow and field ($\langle \mathbf{u}_h \rangle$ and $\langle B_r \rangle$, brackets meaning here the time-average over a considered time-span $\tau = [t_s, t_e]$). Neglecting the quadratic term involving both $\tilde{\mathbf{u}}_h$ and $\tilde{B}_r$, the band-pass filtered SV takes its origin from 5 terms:

$$ \partial_t \tilde{B}_r \approx \tilde{\mathbf{u}}_h \cdot \nabla_h \langle B_r \rangle + \langle B_r \rangle \nabla_h \tilde{\mathbf{u}}_h + \langle \mathbf{u}_h \rangle \cdot \nabla_h \tilde{B}_r + \tilde{B}_r \nabla_h \cdot \langle \mathbf{u}_h \rangle + \tilde{e}_r. \quad (16) $$

and we define

$$ \tilde{\chi}(T_k) = \|\partial_t \tilde{B}_r\|_{\tau} \left( \left\| \tilde{\mathbf{u}}_h \right\|_{\tau} \|\nabla_h \langle B_r \rangle\|_{\tau} + \| \langle B_r \rangle\|_{\tau} \left\| \nabla_h \cdot \tilde{\mathbf{u}}_h \right\|_{\tau} \right. \\ \left. + \| \langle \mathbf{u}_h \rangle \|_{\tau} \left\| \nabla_h \tilde{B}_r \right\|_{\tau} + \left\| \tilde{B}_r \right\|_{\tau} \left\| \nabla_h \cdot \langle \mathbf{u}_h \rangle \right\|_{\tau} + \left\| \tilde{e}_r \right\|_{\tau} \right)^{-1}, \quad (17) $$

where $\|X\|_{\tau} = \frac{1}{t_s-t_e} \int_{t_e}^{t_s} \|X(t)\| dt$ is the time average of the considered norm over the investigated time-span. On interannual time-scales, we anticipate that time-variations of the SV should be primarily due to the

flow changes (and less to field changes, e.g. Lesur et al., 2010; Aubert, 2018), although this might depend on the kind of dynamical event responsible for the observed abrupt SV changes (Aubert et al., 2022). On longer time-scales, magnetic field changes are expected to play a relatively larger role. Errors of representativeness are expected to play a role on all time-scales (see Fig. 2 in Gillet et al., 2019).

We have split our results below into two sections: the impact of the 'ingredients' on the recovered flow at a given epoch (Section 4), and on the time-changes in the flow, with a focus on the best documented satellite era (Section 5).

## 4. Influence of ingredients on the recovery of surface flows at a given epoch

For all our investigations, we choose Kalmag, Neutral top1 dynamo prior for $L_{SV} = 13$ to be our reference and change one factor at a time (choice of prior, field model or $L_{SV}$) to investigate its effect. Only flow models obtained by varying those factors relative to the chosen reference are shown in the figures, in order to focus on the main messages. We investigate first in Section 4.1 the effect of the prior, before we focus in Section 4.2 on the influence of the field model and $L_{SV}$. The comparisons for priors are first carried out for January 2017 (during the Swarm era) to provide the best observational constraint while having all three SV models available. We also examine snapshots in 1970 and 1920 to see the effect of field model data quality, before and after the International Geophysical Year in 1958 that led to much improve the quality and quantity of ground-based magnetic records (Hulot et al., 2015). We invert for core surface flow using pygeodyn for each of the considered set-up using 200 ensemble realisations to estimate the forecast covariance matrix. We produce flow at 1 yearly intervals between 1900 and 1999, and 6 month intervals between 1999 and 2023. Results presented below are based on the ensemble mean model.

### 4.1. Influence of the choice of prior dynamo

For this section we test the different priors with the Kalmag SV field model truncated to $L_{SV} = 13$. We extract from the inverted flow series the snapshot for January 2017. All models closely fit the data, and to a similar level (see SV spectra in Fig. C.1). Furthermore, the inverted errors of representativeness $e_r$ amount to an overall similar contribution to the observed SV for all considered priors.

Flow spectra shown in Fig. 3(a) exhibit significant variability. The $S1^\gamma$ prior presents the largest distance to the reference Neutral top1 for harmonic degrees $l \ge 4$. Flows inverted while considering for the prior Stable top1 (resp. dynamos along the AGF17 path) are slightly closer to the reference on long (resp. short) length-scales. Indeed using a prior from Stable top1 tends to slightly enhance the power towards high harmonic degrees. The power for the difference between the models is relatively large (in comparison with the power for the reference flow) already from spherical harmonic degree $l = 3$. It is comparable with the uncertainty level of inverted flow models, consequence of the large contribution of the errors of representativeness to the observed SV (Barrois et al., 2017; Baerenzung et al., 2016), and to the overall tendency in under-estimating the flow magnitude (Schwaiger et al., 2023). However, this does not prevent from possibly recovering transient motions, as demonstrated by Schwaiger et al. (2024), and discussed below.

These differences are further illustrated in the flow maps for the distance to the reference model shown in Fig. D.1. The difference to the reference for all the flow models based on priors constructed from dynamos along the AGF17 path present a similar magnitude, of the order of a few km/yr (a bit larger for $\mathcal{O}p$). Overall these present relative kinetic energy difference to the reference $\mathcal{O}\lesssim 8\%$ r.m.s., as summarized in Table 3. Within this subset of dynamos, the largest differences are seen within the tangent cylinder and towards the equator. Despite rather different overall set-up (e.g. homogeneous forcing, inclusion of stratified layer), comparable low difference values are found when using
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Spatial spectral energy for each inverted flow model (left), and for the difference to the reference inverted flow model (right) in 2017. (Top): influence of the choice of dynamo prior. (Bottom): influence of the choice of SV field model as data and to the choice of truncation level $L_{SV}$.</figcaption></figure>

**Stable_top1.** For all these dynamos the prior involves both gravitational and EM torques. The largest difference is by far found for flow inverted using the $S1^\gamma$ prior (with $\varrho > 30\%$), the only one considering just an EM torque. Differences there amount up to 10 km/yr locally, with structures much elongated in the azimuthal direction, as illustrated in Fig. D.1. It is reminiscent of the median flow for this dynamo prior (see Section 2 and Fig. 2(a)). This is further supported by the larger power seen for odd degrees on the core surface flow spectra for the $S1^\gamma$ prior shown in Fig. 3 (a). However, there is a lower difference between model coefficients inverted using $S1^\gamma$ and Neutral_top1 as prior, than within those two prior means themselves, which is indicative of the impact of the observational constraints w.r.t. the prior in the core flow inverse problem. An Atlantic-Pacific hemisphere dichotomy is present in all inverted flows (see Appendix D), suggesting that this feature is most likely data-driven and not imposed by the dynamo prior employed.

We summarize in Table 3 the influence of these various ingredients by calculating the kinetic energy $\varrho$ for the difference in flow models, as a percentage of the reference (see Eq. (14)). Regardless of the considered epoch, using dynamo priors across the AGF17 path leads to more difference from the reference than along this path. Flow models obtained using $S1^\gamma$ as the prior always present the largest kinetic energy of flow difference.

## 4.2. Influence of the magnetic field model and their truncation level

We now turn to the influence of the SV model and to their truncation

level $L_{SV}$. For the remainder of this section we test the different field models with the Neutral_top1 prior. We continue to use for the reference the flow model inverted from the Kalmag magnetic model with $L_{SV} = 13$. Field models are less constrained before the satellite era due to the lower precision of the sensors and the uneven global coverage of ground observatories (e.g. Hulot et al., 2015). We investigate how much the larger uncertainties within the field models prior to the satellite era influence the flows. We consider 3 snapshots, 2017, 1970 and 1920, which should cover the major jumps in data quality, in particular the international geophysical year in 1958 (with a net increase in the number of ground observatories and the popularization of proton magnetometers early after), and the continuous satellite monitoring from 1999.

First we present the results for January 2017 in order to sample a Swarm data epoch, while avoiding edge effects that may occur in some field models (in particular for COV-OBS.x2 that ends in 2020). We check that spatial spectra for the predicted SV are extremely similar, and those for the errors of representativeness overlap over all degrees (see Fig. C.1). Overall, all recovered flow models closely fit the SV models considered as observations in the inversion process. No significant difference in the fit to SV coefficients of degree $l \le 13$ is found while changing $L_{SV}$ from 13 to 18, or between CHAOS-7 and Kalmag. A slightly larger prediction misfit (although still small) is observed when inverting COV-OVS-x2 data, which may come from CHAOS-7 and Kalmag using the same SV data error.

Flow spectral energy shown in Fig. 3(c) overlap for all considered
<caption>**Table 3**</caption>

Relative kinetic energy difference $\rho$ from the reference for different choices of field model, dynamo prior, $L_{SV}$ and epochs. The reference model for each time period is the Kalmag, $L_{SV} = 13$, Neutral_top1 prior. Unless explicitly stated on the left column, $L_{SV} = 13$. All boxes containing 'N/A' indicate the combinations that we cannot calculate the kinetic energy for.

<table><thead><tr><th>January 2017</th><th>CHAOS-7</th><th>Kalmag</th><th>COV-OBS.x2</th></tr></thead><tbody><tr><td>71p</td><td>7.5 %</td><td>5.7 %</td><td>6.7 %</td></tr><tr><td>50p</td><td>8.0 %</td><td>6.1 %</td><td>7.4 %</td></tr><tr><td>0p</td><td>7.7 %</td><td>7.8 %</td><td>8.7 %</td></tr><tr><td>Neutral_top1</td><td>2.8 %</td><td>Reference</td><td>1.6 %</td></tr><tr><td>Stable_top1</td><td>11.9 %</td><td>7.3 %</td><td>8.9 %</td></tr><tr><td>S1<sup>γ</sup></td><td>32.7 %</td><td>34.1 %</td><td>36.6 %</td></tr><tr><td>Neutral_top1, L<sub>SV</sub> = 16</td><td>10.0 %</td><td>7.7 %</td><td>N/A</td></tr><tr><td>Neutral_top1, L<sub>SV</sub> = 18</td><td>13.4 %</td><td>8.2 %</td><td>N/A</td></tr></tbody><thead><tr><th>January 1970</th><th>CHAOS-7</th><th>Kalmag</th><th>COV-OBS.x2</th></tr></thead><tbody><tr><td>71p</td><td>N/A</td><td>6.3 %</td><td>10.5 %</td></tr><tr><td>50p</td><td>N/A</td><td>4.8 %</td><td>12.0 %</td></tr><tr><td>0p</td><td>N/A</td><td>8.9 %</td><td>21.7 %</td></tr><tr><td>Neutral_top1</td><td>N/A</td><td>Reference</td><td>7.2 %</td></tr><tr><td>Stable_top1</td><td>N/A</td><td>8.5 %</td><td>17.1 %</td></tr><tr><td>S1<sup>γ</sup></td><td>N/A</td><td>31.0 %</td><td>43.8 %</td></tr><tr><td>Neutral_top1, L<sub>SV</sub> = 16</td><td>N/A</td><td>2.2 %</td><td>N/A</td></tr><tr><td>Neutral_top1, L<sub>SV</sub> = 18</td><td>N/A</td><td>2.7 %</td><td>N/A</td></tr></tbody><thead><tr><th>January 1920</th><th>CHAOS-7</th><th>Kalmag</th><th>COV-OBS.x2</th></tr></thead><tbody><tr><td>71p</td><td>N/A</td><td>6.8 %</td><td>20.1 %</td></tr><tr><td>50p</td><td>N/A</td><td>6.7 %</td><td>20.2 %</td></tr><tr><td>0p</td><td>N/A</td><td>12.9 %</td><td>31.1 %</td></tr><tr><td>Neutral_top1</td><td>N/A</td><td>Reference</td><td>17.4 %</td></tr><tr><td>Stable_top1</td><td>N/A</td><td>10.6 %</td><td>31.0 %</td></tr><tr><td>S1<sup>γ</sup></td><td>N/A</td><td>32.3 %</td><td>57.6 %</td></tr><tr><td>Neutral_top1, L<sub>SV</sub> = 16</td><td>N/A</td><td>2.2 %</td><td>N/A</td></tr><tr><td>Neutral_top1, L<sub>SV</sub> = 18</td><td>N/A</td><td>2.6 %</td><td>N/A</td></tr></tbody></table>

models (with slightly more dispersion for harmonic degrees above 10). Spectra for the flow difference to the reference indicate that increasing $L_{SV}$ from 13 to 16 or 18 has a greater impact on the flow model than considering a different SV field model. The spectral energy of the difference in flow models for various $L_{SV}$ is as large as that for the reference models for $l \ge 10$, suggesting that the flow for larger degrees are less well constrained. Overall the effect on flow spectra of changing the field model (for high quality satellite data) is lower than the effect of changing the prior, which is comparable to changing $L_{SV}$.

The above results are summarized in Table 3 where we collect the relative CMB kinetic energy $\rho$ for the various cases and in Fig. 3(d). Less than 2% difference is found while changing the SV field model, while it amounts up to about 8% when changing prior or increasing $L_{SV}$. This is further described in the spatial domain with maps of the difference (see Fig. D.1). Changing the field model to CHAOS-7 or COV-OBS.x2 produce the smallest difference from the reference ($\le 2$ km/yr). There is not a specific large-scale pattern related to the change of field model but the difference is very symmetrical about the equator. When increasing $L_{SV}$ from the reference, the main flow difference is in the Atlantic hemisphere at small lengthscale structures. The differences in magnitude for $L_{SV} = 18$ are a bit larger than those for $L_{SV} = 16$ (not shown) but show the same spatial patterns.

After considering the effect during the satellite era, we then investigate the effect of changing the field model during periods of less accurate data collection. Appendix D contains flow maps at two additional epochs (1920 and 1970), which show that the difference from the reference increases with increasing time from the present. The weaker effect of $L_{SV}$ during the pre-satellite is related to badly constrained high degree SV coefficients at these epochs: the average SV model tends to be smoothed on short length-scales, and thus the ensemble average of the inverted flows models is not so much sensitive to these small-scale patterns. With lower data quality, the data constraint on the field models is weaker and there is likely to be more influence from the dynamo prior in the core flow inversion.

When accurate satellite data are available, changing the prior has the greatest impact, while the choice of SV observation (field model) has the smallest effect. Meanwhile, increasing the $L_{SV}$ changes the kinetic energy difference to an equivalent amount as changing the prior along the dynamo path. There is thus potential information on core flow models from SV at degrees $l \ge 13$.

The effect of changing the prior along the AGF17 path results in a smaller difference than moving across this path, with the $S1^\gamma$ prior being noticeably offset from all other solutions. This is further illustrated in Fig. 4(a) and (b), where we show the time average inverted flow over the satellite era for both the Neutral_top and $S1^\gamma$. Several differences are visible, such as:

- the presence on an equatorially asymmetric and westward flow below the eastern equatorial Pacific when using $S1^\gamma$ as a prior, whereas it is mainly equatorially symmetric and eastward for Neutral_top1;

- alternating eastward and westward azimuthal motions around the tangent cylinder when using the $S1^\gamma$ prior (shown by light red background at 60° north and south on the orthographic plots in Fig. 4 (a)), which recalls the signature from the median dynamo prior flow (Fig. 2(a));

- stronger flow motions perpendicular to the tangent cylinder with the $S1^\gamma$ prior, while it is closer to a jet-like feature along the tangent cylinder with Neutral_top1.

- much more intense westward polar motions near the tangent cylinder when using the $S1^\gamma$ prior;

This latter point is clearly visible with the latitudinal profiles of the zonal velocity in Fig. 4(c), as we did with the prior flow (Fig. 2(c)). If both models show westward polar vortices, larger fluctuations of this velocity component along $\theta$ are observed when using $S1^\gamma$, together with some more equatorial asymmetry.

Fig. 4(d) shows a comparison of the pdf of the zonal motions at 80° latitudes for the prior dynamos and for the posterior (after inverting SV data) which illustrates this equatorial asymmetry for $S1^\gamma$ well. It also highlights the bias between the two flow solutions, in relation with the prior spread. Significantly stronger westward polar vortices inside the northern tangent cylinder may be related to stronger convective activity (Olson and Aurnou, 1999; Livermore et al., 2017; Lézin et al., 2023) but there is a strong disagreement between our models about the maximum a priori amplitude of this jet. The posterior pdf over the satellite era covers a range of the parameter space in the vicinity of the prior maximum probability when using the Neutral_top1. It is not the case with $S1^\gamma$, for which inverted polar vortex zonal motions (resp. about 17 and 35 km/yr in the southern and northern hemispheres) are much weaker than the most likely values from the prior (resp. $\approx 150$ and 80 km/yr), shown in Fig. 4(d) and Fig. 2(c).

With increasing time from the present, the difference when changing the field model increases and becomes more important for the inversion than the effect of changing the prior. In 1970 the effect of changing field model (COV-OBS.x2 versus Kalmag) is comparable to changing prior from the reference to for instance Stable_top1 or 71p. Conversely, in 1920 the effect of the field model is second largest after changing to the $S1^\gamma$ prior. The one exception to the increasing relative kinetic energy difference with time is the effect of $L_{SV}$. This is related to weakly constrained high degree spherical harmonic coefficients backward in time (the ensemble average flow model is little sensitive to almost unconstrained SV data), leading to a small effect of increasing $L_{SV}$.

## 5. Influence of ingredients on the recovery of transient motions

In addition to the effects on the total ensemble mean flow motions inverted with the pygeodyn algorithm, we document below the effects of choice of the geodynamo prior, of secular variation model, and of $L_{SV}$ on transient flow motions. Firstly, we investigate in Section 5.1 the effect
<figure><img src="image_6.png" /><figcaption>Fig. 4. Maps (with polar views) of median flow over the satellite era (2000–2023) for Kalmag $L_{SV} = 13$ with the $S1^γ$ (a) and the reference Neutral_top1 (b) priors. (c) latitudinal profiles of the zonal velocity in both cases, with solid line for the median over the satellite era, and dashed lines for the 1st and 4th quartiles. Over the same era. (d) distribution of the zonal velocity at 80° latitudes (solid lines: North, dashed lines: South) for the two prior dynamos, and for the inverted flow models. The pdf values for the $S1^γ$ prior is < 0.008.</figcaption></figure>

on a selection of core flow Gauss coefficients and on the LOD over the entire study period (1900–2023). Then we focus in Section 5.2 on the satellite period (1999–2023), where the observational constraint is higher, illustrating the sensitivity of the overall flow change and of interannual wave-like motions. We present below a selection of models relative to the reference that are representative when choosing a different field model, prior or $L_{SV}$. However, we have verified that this relationship holds true when changing the reference and/or multiple ingredients.

## 5.1. Flow coefficients time changes and length-of-day variations

In Fig. 5(a), we show the time evolution of some of the main zonal and non-zonal toroidal coefficients (the same that are also highlighted with the pdfs from the prior dynamos in Fig. B.1, plus two higher degree coefficients). We focus on coefficients representing equatorially symmetric flows, that overall dominate core motions. The overall trends are similar for all considered priors/SV model/$L_{SV}$, although some of the patterns in $c^t_1$ variations are slightly attenuated (e.g. Neutral_top1 with COV-OBS.x2 data) or enhanced (see for instance 71p with Kalmag data), depending on the inversion set-up. Features on top of the trend can also be seen across all analysed models (e.g. the extremum in the $s^t_2$ series around 1980, or decadal oscillation in $c^t_2$). Some offset covering the entire studied era can be detected between two reanalyses using different dynamo priors, for instance on $c^t_3$ from Kalmag data with the 71p and Neutral_top1 priors. The dispersion within the various coefficients series is usually largest for the earliest epochs, reflecting that the improved data quality reduces the spread resulting from the various choices of parameters. The choice of field model does affect the shape of the trend of the coefficient variations, e.g. using the COV-OBS.x2 SV

model a negative peak appears in the late 1930s for the $c^t_3$ coefficient but not until the 1950s for a model inverted from Kalmag. Also, some of the shorter time-scale deviations seen in the Kalmag model are not seen in the COV-OBS.x2 model, likely due to the temporal parameterization used for the construction of this latter (Huder et al., 2020; Pick et al., 2019). The choice of $L_{SV}$ (from 13 to 18) induces little differences on the large-scale flow coefficients, although we witness an off-set for $s^t_2$ over 1990–2015. Importantly, the $S1^γ$ prior does not produce major differences in time changes of these large-scale coefficients, in comparison with the dynamo priors considering a gravitational torque, and this despite displaying a significantly different behaviour in the statistics for the recovered snapshots (see previous Section).

In Table 4 we provide, for each coefficient series shown in Fig. 5, the correlation coefficient between each set-up and the reference (Kalmag, Neutral_top1, $L_{SV} = 13$). The obtained correlations are extremely high for all coefficients and models. Generally, the correlations are slightly less towards larger degrees (e.g. $c^t_5$ and $c^t_7$ to $c^t_1$ and $c^t_3$). The lower correlations are found for the model inverted from COV-OBS.x2. This is the consequence of the worsening data quality backward in time, implying a greater variation between field models, and consequently between the flow coefficients, pre-1970. Indeed, we do not observe the same loss of correlation over the period covered by CHAOS-7, due to the high data quality over the satellite era.

We now discuss how each inversion set-up affects the prediction of geodetic LOD series. Predictions, calculated following Jault and Finlay (2015) as

$$ \delta LOD = 1.232(\delta^c t_1^0 + 1.776 \delta^c t_3^0 + 0.080 \delta^c t_5^0 + 0.002 \delta^c t_7^0) \quad (18) $$

(in ms, where the flow coefficients are measured in km/yr) are shown in
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Effect of choice of parameters (dynamo prior, SV data, *L*<sub>SV</sub>) on some of the main zonal and non-zonal flow coefficients (a) and on the predicted LOD variation (b), over 1900–2023. The geodetic LOD series shown for comparison (black) has been corrected for solid tides, atmospheric contributions when available, and long term trend (for details, see Istas et al., 2023).</figcaption></figure>

Fig. 5(b). For all predicted curves the temporal mean has been removed (as the reference for LOD variations is undetermined). Overall the several models present a very decent fit to the observed decadal LOD changes. As with the flow coefficients documented above, changing the prior appears to affect the vertical scaling of the LOD. For instance 71p shows scaled-up values compared to *Stable_top1* or *Neutral_top1*

(consequence of reduced variations seen in $e^{c_1 t_1^0}$ or $e^{c_3 t_3^0}$). LOD predictions for flows inverted using $S1^γ$ are not significantly different from most other priors, despite its singular signature on spatial spectra. The effect of increasing $L_{SV}$ has almost no impact on the fit to LOD as increasing $L_{SV}$ only has a little effect on the lower degree zonal coefficients that enter the LOD prediction. Curiously, despite the increase in data quality since
<caption>**Table 4**</caption>

Correlation coefficient, for the Gauss coefficients shown in Fig. 5, between each series and the reference from Kalmag, Neutral_top1, L<sub>SV</sub> = 13. For the LOD, the correlation is calculated between the predicted series and the geodetic one.

<table><thead><tr><th>Correlations for series in Fig. 5</th><th>c<sub>1</sub>t<sub>1</sub><sup>0</sup></th><th>c<sub>3</sub>t<sub>3</sub><sup>0</sup></th><th>c<sub>1</sub>t<sub>2</sub><sup>1</sup></th><th>c<sub>3</sub>t<sub>2</sub><sup>1</sup></th><th>c<sub>1</sub>t<sub>5</sub><sup>2</sup></th><th>c<sub>3</sub>t<sub>7</sub><sup>0</sup></th><th>LOD</th></tr></thead><tbody><tr><td>Kalmag, 71p, L<sub>SV</sub> = 13</td><td>0.99</td><td>0.92</td><td>0.99</td><td>0.92</td><td>0.91</td><td>0.91</td><td>0.88</td></tr><tr><td>Kalmag, 50p, L<sub>SV</sub> = 13</td><td>1.00</td><td>0.97</td><td>1.00</td><td>0.94</td><td>0.95</td><td>0.92</td><td>0.88</td></tr><tr><td>Kalmag, Neutral_top1, L<sub>SV</sub> = 13</td><td>Reference</td><td>Reference</td><td>Reference</td><td>Reference</td><td>Reference</td><td>Reference</td><td>0.84</td></tr><tr><td>Kalmag, Stable_top1, L<sub>SV</sub> = 13</td><td>1.00</td><td>0.98</td><td>1.00</td><td>0.98</td><td>0.98</td><td>0.94</td><td>0.85</td></tr><tr><td>Kalmag, S1<sup>Y</sup>, L<sub>SV</sub> = 13</td><td>0.99</td><td>0.94</td><td>0.96</td><td>0.94</td><td>0.97</td><td>0.92</td><td>0.85</td></tr><tr><td>CHAOS-7, Neutral_top1, L<sub>SV</sub> = 13</td><td>0.92</td><td>0.99</td><td>0.99</td><td>0.98</td><td>0.99</td><td>0.98</td><td>0.43</td></tr><tr><td>COV-OBS.x2, Neutral_top1, L<sub>SV</sub> = 13</td><td>0.97</td><td>0.87</td><td>0.95</td><td>0.71</td><td>0.76</td><td>0.80</td><td>0.84</td></tr><tr><td>Kalmag, Neutral_top1, L<sub>SV</sub> = 18</td><td>1.00</td><td>0.99</td><td>1.0</td><td>0.93</td><td>0.99</td><td>0.98</td><td>0.85</td></tr></tbody></table>

the start of the satellite era, the best fit to the LOD series does not occur since 2000: features like the peak in LOD in the late 2000s are missing, and none of the models predict the significant decrease measured before 2020.

We show in Table 4 the correlation between the observed and predicted LOD. Correlations obtained comparing LOD predictions and observations are lower than when comparing flow coefficients obtained with different set-ups, due to the flow coefficients being calculated with the same methodology whereas the LOD dataset is independently observed. A correlation of around 0.85 is found for all models starting in 1900. The correlation involving CHAOS-7 is significantly lower as it is only calculated over the past 24 years, but this is true for all models. Further work is needed to investigate why the LOD fit is not improved over the satellite era where the observational constraint is best. One

possible reason may simply be associated with the rather large uncertainties inherent to the core flow inverse problem, since a large fraction of the SV is due to unresolved processes (Gillet et al., 2019). Large errorbars on decadal LOD predictions are indeed proposed in many studies, without being significantly reduced over the satellite era (see for instance Bärenzung et al., 2018; Finlay et al., 2023).

## 5.2. Flow change and interannual waves over the satellite era

Now we focus on changes in flow motions over the satellite era (2000–2023), that benefit from the best observational constraint. We concentrate on the azimuthal component $u_\phi$ at the equator, where some of the strongest flow variations have been previously documented (Finlay et al., 2023). Indeed this area is of particular interest as it is at

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Time-longitude diagrams at the equator for the azimuthal flow anomaly w.r.t. the time-average over the satellite period 2000–2023, $\Delta u_\phi = u_\phi - \langle u_\phi \rangle$, for a selection of model parameter choices. The reference model, shown in the top left grey box, is inverted from Kalmag using L<sub>SV</sub> = 13 and the Neutral_top1 prior. The effect of changing L<sub>SV</sub> is shown in the second left-hand row. The effect of changing the field model is shown in the bottom left. The effect of changing the dynamo prior is shown in the right-hand column.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Time-dependent azimuthal motions at the equator, bandpassed for periods between 4.5 and 9.5 years, over the satellite period (2000–2023), for a selection of inversion set-ups. The reference model (shown in the top left grey box) is inverted from Kalmag using $L_{SV} = 13$ and the Neutral_top1 prior. The effect of changing $L_{SV}$ is shown in the second left-hand row. The effect of changing the field model is shown in the bottom left. The effect of changing the dynamo prior is shown in the right hand column.</figcaption></figure>

<caption>**Table 6**</caption>

Correlation coefficient for the flow minus its temporal mean (top) or for band-passed flow between 4.5 and 9.5 years (bottom) for the equatorial time-longitude plots (Fig. 6 and 7 respectively) over the satellite period (2000–2023), w.r.t. the reference flow model for different model parameters (dynamo prior, SV data, $L_{SV}$). Unless explicitly stated on the left column, $L_{SV} = 13$ is used.

<table><thead><tr><th>Equatorial T-L Correlation<br>Coefficient u<sub>φ</sub> - &lt;u<sub>φ</sub>&gt;</th><th>CHAOS-7</th><th>Kalmag</th><th>COV-OBS.x2</th></tr></thead><tbody><tr><td>71p</td><td>0.97</td><td>0.98</td><td>0.92</td></tr><tr><td>50p</td><td>0.94</td><td>0.98</td><td>0.94</td></tr><tr><td>0p</td><td>0.98</td><td>0.98</td><td>0.95</td></tr><tr><td>Neutral_top1</td><td>0.98</td><td>Reference</td><td>0.94</td></tr><tr><td>Stable_top1</td><td>0.97</td><td>0.99</td><td>0.96</td></tr><tr><td>S1<sup>γ</sup></td><td>0.94</td><td>0.95</td><td>0.97</td></tr><tr><td>Neutral_top1, L<sub>SV</sub> = 16</td><td>0.97</td><td>0.99</td><td>N/A</td></tr><tr><td>Neutral_top1, L<sub>SV</sub> = 18</td><td>0.97</td><td>0.99</td><td>N/A</td></tr><tr><td>Same as above when u&#771;<sub>φ</sub><br>for T &isin; [4.5, 9.5] yr</td><td>CHAOS-7</td><td>Kalmag</td><td>COV-OBS.x2</td></tr><tr><td>71p</td><td>0.69</td><td>0.95</td><td>0.72</td></tr><tr><td>50p</td><td>0.66</td><td>0.97</td><td>0.82</td></tr><tr><td>0p</td><td>0.83</td><td>0.99</td><td>0.87</td></tr><tr><td>Neutral_top1</td><td>0.83</td><td>Reference</td><td>0.87</td></tr><tr><td>Stable_top1</td><td>0.8</td><td>0.98</td><td>0.83</td></tr><tr><td>S1<sup>γ</sup></td><td>0.73</td><td>0.97</td><td>0.83</td></tr><tr><td>Neutral_top1, L<sub>SV</sub> = 16</td><td>0.79</td><td>0.99</td><td>N/A</td></tr><tr><td>Neutral_top1, L<sub>SV</sub> = 18</td><td>0.83</td><td>0.99</td><td>N/A</td></tr></tbody></table>

inverted using a priori statistics from various numerical geodynamo simulation. It suggests a flow model sensitivity to observational constraints stronger towards short time-scales, despite a decreasing amplitude of core motions with frequency. Snapshot flows are indeed dominated by long periods, because of the red temporal spectrum of core motions: more power towards long periods, consequence of a similar property for the SV (see Gillet et al., 2015; Lesur et al., 2018). Said differently, the slowly evolving ‘background’ flow would tend to fall in the null space of the observation operator in Eq. (10), contrary to rapid transient motions.

In order to test this hypothesis, we show in Fig. 8 the diagnostic $\tilde{\chi}(T_k)$ for 6 bandpass windows in $T_k \in [5, 20]$ yr (see Eq. (17)) for flow models over the satellite era (here from the Kalmag field model using 71p and $L_{SV} = 13$). Order-2 Butterworth bandpass filter windows are computed for each period $T_k$ over the bandwidth $[2T_k/3, 3T_k/2]$, with

$$ \log(T_k) = \log(T_{inf}) + (\log(T_{up}) - \log(T_{inf})) / k, \quad (19) $$

so that each band covers a similar spectral range on a log scale, with $T_{up} = 20$ yr and $T_{inf} = 5$ yr (or 10 and 2.5 $\tau_A$, considering $\tau_A = 2$ yr in Earth’s core). $\tilde{\chi}(T_k)$ measures the efficiency of SV generation from a given flow and configuration of $B_r$ at the CMB for bandpassed motions over the satellite era. We observe a clear decrease of $\tilde{\chi}$ with the period, which we interpret as the observable signature of a dynamo field that tends, for long time-scales, towards a state where $\nabla \times (\mathbf{u} \times \mathbf{B})$ is much less in amplitude than a potential e.m.f. amplitude $UB/\ell$ based on a typical length-scale and magnitudes for the field and flow (see § 3.2). This condition is never exactly achieved in practice (see the value $\chi \sim$
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Top: efficiency of the flow to generate SV for a given field model $B_r$, as measured by the diagnostic $\chi$ (Eq. (15)) for snapshot flows and $\tilde{\chi}$ (Eq. (17)) for bandpassed motions from the inverted model obtained using Kalmag with $L_{SV} = 13$ and the 71p prior, over the satellite era. See text for details. Bottom: relative fraction for the 5 contributions to the SV norm in the denominator of Eq. (17).</figcaption></figure>

0.3 found for snapshot motions), likely because the geodynamo is perpetually disturbed on a broad range of periods. Our results resonate with the findings by Rau et al. (2000), who put forward that most of the snapshot flow lies within the null space of the core flow inverse problem in dynamo simulations. Regardless of truncation degree $L_U^\dagger$ used for the flow in Eqs. (15) and (17), the values of $\tilde{\chi}$ are larger than for $\chi$. If wave-like motions clearly live on a quasi-steady magnetic state with high values of $\tilde{\chi}$ on periods of a few $\tau_A$, it is not the case for slower motions for which the fields have time to adapt in order to reduce induction, thus reducing $\tilde{\chi}$, consequence of a dynamical adjustment. We check in Fig. 8 (top) that the above trend is found for various flow truncation degrees, mitigating the effect of the finite resolution in flow models. The diagnostic $\tilde{\chi}$ naturally decreases when increasing $L_U^\dagger$, because this inflates the denominator in Eq. (17). Induction associated with the large scale field and flow is dominated by the contribution from the time-changing flow for periods from 5 to 20 yr. Its amplitude is comparable to that of $e_r$, as already shown in Gillet et al. (2019). The part associated with the transient field gradually increases towards lower frequencies (see the green bars in Fig. 8, bottom).

A discussion along the same line is presented by Aubert et al. (2022), where wave-like and convective motions are clearly distinguished. The analysis proposed here puts forward a smooth transition between the two. This is illustrated in Fig. 9 where we show $\tilde{\chi}(T_k)$ and $\chi$ as estimated from the 71p dynamo. We consider periods $T_k$ ranging from $2\tau_A$ to $200\tau_A$. We observe that $\tilde{\chi}$ continuously decreases (whatever the considered truncation levels) from periods of a few $\tau_A$ onward, to reach

values as low as 0.1. Values on short timescales are somewhat weaker than the same diagnostic from geophysical field and flow models. This may be explained by a lesser amount of high energy jerk events in the 71p dynamo, generating an overall r.m.s. wave motions weaker in this simulation than observed in the Earth, as discussed in Aubert and Gillet (2021). The contribution of errors of representativeness to the denominator in Eq. (17) is relatively little, in comparison with what is observed for the inverted solutions, shown in Fig. 8. It may be that the augmented state approach from pygeodyn tends to magnify the role of $e_r$. It may also be that with the (imperfectly) inverted field and flow models, the norms associated with the nonlinear terms are downgraded. Over all considered periods, the main driver of the transient SV is associated with the terms from $u_h \cdot \nabla_h B_r$. The dominant role of transient motions on short time-scales is smoothly replaced by that of the time-varying field on periods longer than $\sim 30\tau_A$ (see the relative contributions to the denominator of Eq. (17) from green and red bars in Fig. 9, bottom).

We can postulate core dynamics in a scenario where damped hydro-magnetic waves travel above the magnetic field over a broad range of time-scales. These tend to adjust the field and flow to a situation where the two would generate as little induction as possible. The situation for inverted surface motions is a bit different, as the flow in the null space (i. e. generating no radial SV) helps enforce the a priori constraint. However, this cannot explain the decrease of the diagnostic $\tilde{\chi}$ towards short time-scales, because both slow and rapid flow changes are essentially associated with the same observation operator. Our present analysis, anchored to observational constraints, is limited to the core surface. It
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Efficiency of the flow to generate SV for a given field model $B_r$, as measured by the diagnostic $\chi$ (Eq. (15)) for snapshot flows and $\tilde{\chi}$ (Eq. (17)) for bandpassed motions from the 71p dynamo. Bandpass windows are computed from Eq. (19), with $T_{up} = 200\tau_A$ and $T_{inf} = 2\tau_A$. Diagnostics have been calculated over the time-span referred to as 5000–10,000 years in Aubert and Gillet (2021). The snapshot estimate is obtained as the mean over 30 independent epochs. Bottom: relative fraction for the 5 contributions to the SV norm in the denominator of Eq. (17).</figcaption></figure>

calls for dedicated numerical studies of the induction process in the bulk of the core as a function of time-scale.

## 6.2. Testing geodynamical hypotheses with magnetic data

We have presented here a confrontation of magnetic data with dynamo simulations, via the core flow inverse problem. Our discussion involved various computations run in rather distinct geodynamical configurations (different torques set-ups, axial symmetry of the forcing, or stratification levels at the core surface). Even though all considered priors allow to fit magnetic data equally well (including the outlier $S1'$), we have illustrated how the posterior flow distribution can in some locations be in agreement or not with the dynamo prior (recall Fig. 4). With such kind of analysis one would tend to consider the $S1'$ prior as unlikely, calling for alternative hypotheses – e.g. re-instating an EM coupling with the inner core, possibly rotating and/or gravitationally coupled to the mantle. This calls for further hypothesis testing regarding the Earth-likeness of geodynamo models. To this purpose, it would be insightful to revisit the various coupling set-ups explored for instance by Pichon et al. (2016), but at more extreme parameters.

Indeed, the input parameters used for most of the dynamos (Ek, $P_m$) reduce the frequency range over which these simulations can be compared to the Earth (Aubert and Gillet, 2021). It makes impossible to

use their temporal statistics on short time-scales while generating the reduced stochastic models in *pygeodyn*, as operated with 71p by Gillet et al. (2019) and Istas et al. (2023), because the turn-over and Alfvén time-scales are not separated enough. We stress that future comparisons of prior set-ups would need dynamo models run to conditions where spatio-temporal statistics can be constructed, for extracting information on the transient dynamics from the priors. We believe our study stresses the importance and need to construct a pool of accessible coefficients from geodynamo priors for hypothesis testing. This is one target of the IAGA Fair-Mod Working Group, as of IUGG 2023 meeting (IUGG AR23, 2024).

## 7. Conclusions

We proposed an attempt at quantifying the impact of prior dynamo information and the magnetic data constraint when inverting for core motions using the *pygeodyn* algorithm. Our results indicate that regardless of the considered dynamo prior, the recovery of interannual wave dynamics from satellite measurements (Gillet et al., 2022; Ropp and Lesur, 2023; Li et al., 2024) appears robust. So is the overall decadal trend over the satellite era, with in particular the emergence of an Eastward jet under the equatorial Pacific (Ropp and Lesur, 2023; Finlay et al., 2023; Firsov et al., 2023). The quasi-static component of flow
motions is more strongly influenced by our choice of prior, especially to the choice of core-mantle interaction (e.g. the existence of gravitational torque involving the inner core, and the presence or not of a persistent westward drift, see Pichon et al., 2016). This is a consequence of slow motions generating a weak induction relative to their potential, consequence of a dynamical adjustment. This adjustment is operated by transient motions, likely quasi-geostrophic Alfvén and Magneto-Coriolis waves at the conditions of Earth's core (Gillet et al., 2022, 2024; Aubert et al., 2022; Gerick and Livermore, 2024). These tend to bring back the perturbed field and flow towards the background state. Fast waves, almost equipartitioned in magnetic and kinetic energy for periods close to $τ_A$, travel over the background field to propagate the information of a perturbation. Longer period waves, more magnetic than kinetic, may also alter the field on long time-scales.

Prior to the satellite era, the effect of changing the field model is comparable to changing prior along the dynamo path proposed by Aubert et al. (2017). The effect of moving across this path seems to be a dominant factor regardless of the considered snapshot. We also show that some information on the core dynamics can potentially be gained by considering the secular variation at spherical harmonic degrees above 13. This part of the SV signal is only accessible over the period covered by satellite data, although restricted to a time-filtered picture towards shorter length-scales (e.g. Baerenzung et al., 2022), because of the ill-posed nature of the geomagnetic inverse problem where the magnetic signal steeply decreases with the harmonic degree. Accessing such short length-scales on long time-scales will require to extend the era covered by highly accurate satellite data. Further work on building the optimal geomagnetic field model and numerical dynamo prior is required. Nevertheless, our results show consistent fit to the geomagnetic readings and many similar flow features. Finally, it is important to remember that these flow inversions are dependent on accurate field measurements from satellite and ground observatories, and further support to ensure continuous high quality satellite measurements is required over the Swarm mission and beyond.

### CRediT authorship contribution statement

**Hannah F. Rogers:** Writing – review & editing, Writing – original draft, Methodology, Investigation, Formal analysis. **Nicolas Gillet:** Writing – review & editing, Writing – original draft, Supervision, Methodology, Formal analysis, Conceptualization. **Julien Aubert:** Writing – review & editing, Software, Data curation. **Paolo Person-nettaz:** Writing – review & editing, Software, Data curation. **Mioara Mandea:** Writing – review & editing, Supervision, Conceptualization.

### Declaration of competing interest

None.

### Acknowledgments

We thank Hagay Amit and an anonymous referee for their thorough reviews and constructive comments about our manuscript. We thank Nathanaël Schaeffer for the access to the XSHELLS GPU code (Schaeffer et al., 2017; Schaeffer, 2013), to run the $S1^γ$ dynamo, and for the fruitful discussions. The $S1^γ$ dynamo computations were performed using HPC resources of IDRIS under allocation A0140407382 attributed by GENCI (Grand Equipement National de Calcul Intensif). All other dynamo performed at S- CAPAD, IPGP and using HPC resources from GENCI-TGCC and, GENCI-IDRIS and GENCI-CINES (grant number A0140402122). 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 financially supported by CNES as an application of the Core_Swarm mission. HR, NG and PP are part of Labex OSUG@2020 (ANR10 LABX56).

## Appendix A. Summary of mathematical notations

We split our mathematical notation into two tables: (i) Table A.1 contains all functions that can be applied to an example variable $x$, and (ii) Table A.2 summarizes the mathematical symbols used.

<caption>**Table A.1**</caption>

Summary of the functions that can be applied to an example variable $x$.

<table><thead><tr><th>Function</th><th>Description</th></tr></thead><tbody><tr><td>$\dot{x}$</td><td>Time derivative of x (sometimes $\dot{\phi}_x$)</td></tr><tr><td>$||x||$</td><td>Norm of x (see Eqs. (4), (9))</td></tr><tr><td>$\bar{x}$</td><td>Large length-scale projection of x</td></tr><tr><td>$\Delta x$</td><td>Distance to the mean of x</td></tr><tr><td>$x^*$</td><td>x extracted from the dynamo prior</td></tr><tr><td>$\langle x \rangle$</td><td>Temporal mean of x over a time-span $[t_s, t_e]$</td></tr><tr><td>$\tilde{x}(T)$</td><td>x band-pass filtered around a period T</td></tr></tbody></table>
<table><caption>**Table A.2**
Summary of the mathematical notations.</caption><thead><tr><th>Symbol</th><th>Description</th><th>Symbol</th><th>Description</th></tr></thead><tbody><tr><td>a</td><td>Earth's reference radius, 6371.2 km</td><td>S<sub>B</sub>(l)</td><td>MF spatial spectrum at Earth's surface</td></tr><tr><td>A(b)</td><td>Gaunt-Elsasser matrix</td><td>S<sub>SV</sub>(l)</td><td>SV spatial spectrum at Earth's surface</td></tr><tr><td>b</td><td>Inner core radius, 1220 km</td><td>S<sub>U</sub>(l)</td><td>Spatial spectrum of the core surface kinetic energy</td></tr><tr><td>b</td><td>Vector of MF spherical harmonic coefficients</td><td>c<sub>s</sub>g<sup>m</sup><sub>l</sub></td><td>Gauss coefficients for the poloidal scalar</td></tr><tr><td>b̂</td><td>Vector of SV spherical harmonic coefficients</td><td>c<sub>s</sub>t<sup>m</sup><sub>l</sub></td><td>Gauss coefficients for the toroidal scalar</td></tr><tr><td>B</td><td>Typical magnetic field intensity in the core</td><td>T(θ, φ)</td><td>surface flow toroidal scalars</td></tr><tr><td>B<sub>u,e</sub></td><td>Diffusion matrix for time evolution of u or e</td><td>T<sub>k</sub></td><td>Center period for bandpass filter</td></tr><tr><td>B(r, θ, φ)</td><td>magnetic field vector</td><td>U</td><td>Typical flow velocity</td></tr><tr><td>B<sub>r</sub>(r, θ, φ)</td><td>Radial component of the magnetic field</td><td>U<sub>A</sub></td><td>Typical Alfvén speed</td></tr><tr><td>c</td><td>Outer core radius, 3483 km</td><td>u</td><td>Vector of core surface flow coefficients</td></tr><tr><td>D</td><td>gap between outer and inner core boundaries</td><td>u<sub>h</sub>(θ, φ)</td><td>Horizontal core motion at the CMB</td></tr><tr><td>D<sub>u,e</sub></td><td>Drift matrix in AR models for u or e</td><td>w<sub>u,e</sub></td><td>random forcing vector in the AR equations for u and e</td></tr><tr><td>e</td><td>Vector of spherical harmonic coefficients for e<sub>r</sub></td><td>z</td><td>Augmented state for pygeodyn, z = |u, e|</td></tr><tr><td>e<sub>r</sub>(θ, φ)</td><td>Errors of representativeness</td><td>Δm</td><td>Thickness of the electrically conducting mantle layer</td></tr><tr><td>F</td><td>Mass anomaly flux</td><td>Δt<sup>a</sup></td><td>Analysis timestep for pygeodyn</td></tr><tr><td>g<sub>0</sub></td><td>Gravitational acceleration at the CMB</td><td>Δt<sup>f</sup></td><td>Forecast timestep for pygeodyn</td></tr><tr><td>g<sup>m</sup><sub>S</sub>, h<sup>m</sup><sub>l</sub></td><td>Schmidt semi-normalized MF Gauss coefficients</td><td>η</td><td>Magnetic Diffusivity</td></tr><tr><td>g<sup>m</sup><sub>S</sub>, h<sup>m</sup><sub>l</sub></td><td>Schmidt semi-normalized SV Gauss coefficients</td><td>κ</td><td>Co-density diffusivity</td></tr><tr><td>l<sub>u,e</sub></td><td>Identity matrix of the same dimensions as u or e</td><td>λ</td><td>G-LASSO free-parameter</td></tr><tr><td>l</td><td>Spherical harmonic degree</td><td>μ</td><td>Magnetic permeability of free-space</td></tr><tr><td>ℓ</td><td>Characteristic length-scale for the flow</td><td>ν</td><td>Kinematic viscosity</td></tr><tr><td>L<sub>B</sub></td><td>MF truncation degree in pygeodyn</td><td>ρ</td><td>Core density</td></tr><tr><td>L<sup>l</sup><sub>B</sub></td><td>MF truncation degree in χ and χ̃</td><td>ρ'</td><td>radial density anomaly profile</td></tr><tr><td>L<sub>SV</sub></td><td>SV truncation degree in pygeodyn</td><td>Q</td><td>Normalized core surface Kinetic energy difference, Eq. (14)</td></tr><tr><td>L<sup>l</sup><sub>SV</sub></td><td>SV truncation degree in χ and χ̃</td><td>Σ</td><td>Ratio of mantle to core conductances</td></tr><tr><td>L<sub>U</sub></td><td>Core flow truncation degree in pygeodyn</td><td>σ<sub>c</sub></td><td>Electrical conductivity of the core</td></tr><tr><td>L<sup>l</sup><sub>U</sub></td><td>Core flow truncation degree in χ and χ̃</td><td>σ<sub>m</sub></td><td>Electrical conductivity of the mantle</td></tr><tr><td>m</td><td>Spherical harmonic order</td><td>τ</td><td>timespan of reanalysis</td></tr><tr><td>N(r)</td><td>Brünt-Väisala frequency</td><td>τ<sub>A</sub></td><td>Alfvén time</td></tr><tr><td>N<sup>a</sup></td><td>Ensemble size in pygeodyn</td><td>τ<sub>SV</sub>(l)</td><td>SV timescale</td></tr><tr><td>N<sup>f</sup></td><td>Number of samples from dynamo prior</td><td>τ<sub>U</sub></td><td>Vortex overturn time</td></tr><tr><td>P<sup>m</sup><sub>l</sub>(cosθ)</td><td>Legendre polynomials</td><td>τ<sub>u,e</sub><sup>AR</sup></td><td>restoring time of the drift matrix for u and e</td></tr><tr><td>P<sub>ee</sub></td><td>cross-covariance matrix of e from dynamo prior</td><td>τ<sub>η</sub></td><td>Magnetic diffusion time</td></tr><tr><td>P<sub>uu</sub></td><td>cross-covariance matrix of u from dynamo prior</td><td>τ<sub>Ω</sub></td><td>Rotation timescale</td></tr><tr><td>P<sup>f</sup><sub>zz</sub></td><td>Forecast covariance matrix for z</td><td>χ</td><td>Normalized measure of the SV, Eq. (15)</td></tr><tr><td>R<sub>Ar</sub></td><td>Flux Rayleigh Number</td><td>χ̃(T<sub>k</sub>)</td><td>Normalized measure of the time filtered SV, Eq. (17)</td></tr><tr><td>S(θ, φ)</td><td>Poloidal scalars for the core surface flow</td><td>Ω</td><td>Rotation rate of the planet</td></tr></tbody></table>

## Appendix B. Comparison of the dynamo priors

Fig. B.1 presents the probability density function (pdf) of the largest toroidal zonal ($c_{l_1}^0$ against $c_{l_3}^0$) and non-zonal ($c_{l_2}^1$ against $c_{l_2}^1$) components. These latter two coefficients capture a significant fraction of the large eccentric westward gyre imaged from satellite magnetic observations (Pais and Jault, 2008; Gillet et al., 2015; Bärenzung et al., 2018). Their joint distribution highlights some aspects of the forcing imposed to each dynamo. All dynamo along the AGF17 path are centered on the same values, and present a comparable spread, only slightly increased between *Op* and *50p*. PDFs are clearly off-set w.r.t. 0, characteristic of the mean westward drift (consequence of the gravitational and EM couplings, see Pichon et al., 2016). The off-set $c_{l_2}^1$ is reminiscent of the eccentric gyre predominant in the Atlantic hemisphere for the path dynamos. Note the plateau in the distribution of the $c_{s}t_{2}^{1}$, clear for 71p, characteristic of an azimuthal drift of the eccentric gyre in this dynamo.
<figure><img src="image_12.png" /><figcaption>**Fig. B.1.** Top: probability density function of the largest (equatorially symmetric) zonal (left) and non-zonal (right) coefficients for dynamos along the AGF17 path and *Neutral_top1*. Bottom: the same for dynamos across the AGF17 path and *Op*. The same color refers to the same dynamo in all subplots. The solid lines and dashed lines show one and two standard deviations from the mean respectively.</figcaption></figure>

The same is found for the *Stable_top1* dynamo, except that the pdf of $c_2^s t_2^1$ is centered on 0, because the forcing for this dynamo is isotropic instead. Dynamo off the AGF17 path and coupled to the mantle show a $c_1^t_1$ slightly less negative in comparison with the path dynamo. As witnessed on the power spectra at low degrees, the stratified dynamo is characterized by a reduced spread of the pdfs. But the most striking difference between all these priors is associated with the distribution observed for $S1^\gamma$. The $S1^\gamma$ pdfs reflect the mean and variance spectra above by showing a larger spread compared to the dynamos with gravitational torques as well as having a larger negative $c_3^t$ coefficient. The distribution of the $t_2^1$ coefficients are centered on 0, reflecting the isotropic forcing and boundary conditions, but the lack of circular structures in the $S1^\gamma$ pdf may reflect a too short time series for obtaining converged statistics – especially since it presents millennial scale variations, as found with the $S1$ simulation by Schaeffer et al. (2017).

We repeated such analysis for the errors of representativeness from the several dynamo prior (Fig. B.2). The spectra for $e_r$ show similar order of magnitude within the various dynamos, although it is a bit weaker for *Stable_top1* over all length-scales, and a little larger for *Neutral_top1* and $S1^\gamma$ on large length-scales. The spatial spectrum for the spread is larger than that of the mean for all degrees. The pdfs broadly overlap for all of the coefficients and show that for this aspect, the priors are rather similar over all of the spherical harmonic degrees and for the $S1^\gamma$ prior the pdfs are not larger than those with gravitational torques.
<figure><img src="image_13.png" /><figcaption>**Fig. B.2.** Top: spatial power spectra (Eq. (3)) for the time average errors of representativeness (left) and the departure from this average (right) of each prior considered in this study. Middle: probability density function of the largest zonal (left) and some sectorial (right) coefficients for the AGF17 path dynamo and *Neutral_top1*. Bottom: the same for dynamos across the AGF17 path and *Op*. The same color refers to the same dynamo in all subplots. The solid lines and dashed lines show one and two standard deviations from the mean respectively.</figcaption></figure>

## Appendix C. Comparison of SV spectra

We provide the SV spectra at the Earth's surface in Fig. C.1 for the predicted SV, the errors of representativeness and the SV prediction error for various dynamo priors and choices for the SV data and its truncation level $L_{SV}$.
<figure><img src="image_14.png" /><figcaption>(b) SV and ER for Changing SV Model in January 2017</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. C.1.** Spatial spectra of the predicted SV (solid circled lines), the errors of representativeness (dashed circled lines), and the model prediction errors (dotted circled lines), for epoch 2017. Left: influence of the choice of dynamo prior. Right: influence of the choice of SV model considered as data and of truncation level $L_{SV}$.</figcaption></figure>

## Appendix D. Flow maps for the difference to the reference inverted flow

Here we present flow maps at the CMB from the reference inverted flow model. Additional to the relative kinetic energy statistics presented in Table 3, we also display maps for the differences to the reference case while varying the dynamo prior, the SV data and its truncation level $L_{SV}$. These are provided for three epochs representative of the satellite era (2017, Fig. D.1), an epoch after 1958 the international geomagnetic year (1970, Fig. D.2), and one prior to the advent of modern instruments (1920, Fig. D.3).

<figure><img src="image_16.png" /><figcaption>**Fig. D.1.** Influence of SV data and choice of dynamo prior for flow models inverted in 2017. Grey box: flow maps of the reference model (Kalmag with the *Neutral_top1* prior and $L_{SV}$ = 13). Top: difference to the reference model while changing the dynamo prior along the AGF17 path. Middle: difference to the reference model while changing the dynamo prior across the AGF17 path. Bottom: difference to the reference model while changing the SV field model and $L_{SV}$. The colours show the magnitude of $u_\phi$ and the streamlines show the flow speed and direction. Note that the scale for maps of the reference model $u_\phi$ is twice larger than for maps of the difference to the reference (but the same magnitude is used for streamlines).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. D.2. Same as Fig. D.1 in 1970.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. D.3. Same as Fig. D.1 in 1920.</figcaption></figure>
## Data availability

The pygeodyn data assimilation code is available from https://geodyn.univ-grenoble-alpes.fr/.

The dynamo priors are accessible from https://gricad-gitlab.univ-grenoble-alpes.fr/Geodynamo/pygeodyn_data

The core flow inversions will be available at a gricad-gitlab repository when the paper is accepted.

Also the zenodo link to the model repositories and plotting files are: https://doi.org/10.5281/zenodo.15020035

## References

Amit, H., Olson, P., 2004. Helical core flow from geomagnetic secular variation. Phys. Earth Planet. Inter. 147, 1–25.

Aubert, J., 2005. Steady zonal flows in spherical shell dynamos. J. Fluid Mech. 542, 53–67.

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

Aubert, J., 2019. Approaching Earth's core conditions in high-resolution geodynamo simulations. Geophys. J. Int. 219 (Supplement 1), S137–S151. https://doi.org/10.1093/gji/ggz232.

Aubert, J., Gillet, N., 2021. The interplay of fast waves and slow convection in geodynamo simulations near Earth's core conditions. Geophys. J. Int. 225 (3), 1854–1873. ISSN 0956-540X. https://doi.org/10.1093/gji/ggab054.

Aubert, J., Aurnou, J., Wicht, J., 2008. The magnetic structure of convection-driven numerical dynamos. Geophys. J. Int. 172 (3), 945–956. ISSN 0956540X. https://doi.org/10.1111/j.1365-246X.2007.03693.x.

Aubert, J., Finlay, C.C., Fournier, A., 2013. Bottom-up control of geomagnetic secular variation by the Earth's inner core. Nature 502 (7470), 219–223. https://doi.org/10.1038/nature12574. URL: http://www.nature.com/doifinder/10.1038/nature12574.

Aubert, J., Gastine, T., Fournier, A., 2017. Spherical convective dynamics in the rapidly rotating asymptotic regime. J. Fluid Mech. 813, 558–593. ISSN 0022-1120. https://doi.org/10.1017/jfm.2016.789.

Aubert, J., Livermore, P.W., Finlay, C.C., Fournier, A., Gillet, N., 2022. A taxonomy of simulated geomagnetic jerks. Geophys. J. Int. 231 (1), 650–672. ISSN 0956-540X. https://doi.org/10.1093/gji/ggac212.

Backus, G., 1968. Kinematics of geomagnetic secular variation in a perfectly conducting core, philosophical transactions of the Royal Society of London a: mathematical. Phys. Eng. Sci. 263 (1141), 239–266. ISSN 0080-4614. https://doi.org/10.1098/rsta.1968.0014. URL: http://rsta.royalsocietypublishing.org/content/263/1141/239.

Backus, G.E., Mouël, J.-L.L., 1986. The region on the core—mantle boundary where a geostrophic velocity field can be determined from frozen-flux magnetic data. Geophys. J. Int. 85 (3), 617–628.

Baerenzung, J., Holschneider, M., Lesur, V., 2016. The flow at the Earth's core-mantle boundary under weak prior constraints. J. Geophys. Res. Solid Earth 121 (3), 1343–1364. ISSN 2169-9356. https://doi.org/10.1002/2015JB012464.

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), 163. ISSN 1880-5981. https://doi.org/10.1186/s40623-020-01295-y.

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

Bärenzung, 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 (6), 4539–4560.

Barrois, O., Gillet, N., Aubert, J., 2017. Contributions to the geomagnetic secular variation from a reanalysis of core surface dynamics. Geophys. J. Int. 50–68. ISSN 0956-540X. URL: https://academic.oup.com/gji/article-lookup/doi/10.1093/gji/ggx280.

Bloxham, J., Jackson, A., 1991. Fluid flow near the surface of Earth's outer core. Rev. Geophys. 29 (1), 97–120.

Bouligand, C., Gillet, N., Jault, D., Schaeffer, N., Fournier, A., Aubert, J., 2016. Frequency spectrum of the geomagnetic field harmonic coefficients from dynamo simulations. Geophys. J. Int. 207 (2), 1142–1157.

Braginsky, S.I., Roberts, P.H., 1995. Equations governing convection in earth's core and the geodynamo. Geophys. Astrophys. Fluid Dyn. 79 (1–4), 1–97. ISSN 0309-1929. https://doi.org/10.1080/03091929508228992.

Buffett, B., Knezek, N., Holme, R., 2016. Evidence for MAC waves at the top of Earth's core and implications for length of day. Geophys. J. Int. 204 (3), 1789–1800. https://doi.org/10.1093/gji/ggv552.

Christensen, U.R., Aubert, J., Hulot, G., 2010. Conditions for earth-like geodynamo models. Earth Planet. Sci. Lett. 296 (3–4), 487–496. https://doi.org/10.1016/j.epsl.2010.06.009.

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 (3), 1531–1543. https://doi.org/10.1002/2013JB010604.

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 (9), 678–685. https://doi.org/10.1038/ngeo2492.

Davies, C.J., Bono, R.K., Meduri, D.G., Aubert, J., Greenwood, S., Biggin, A.J., 2021. Dynamo constraints on the long-term evolution of Earth's magnetic field strength. Geophys. J. Int. 228 (1), 316–336. https://doi.org/10.1093/gji/ggab342.

de Wijns, G.A., Kresse, G., Vocadlo, L., Dobson, D., Alfe, D., Gillan, M.J., Price, G.D., 1998. The viscosity of liquid iron at the physical conditions of the Earth's core. Nature 392 (6678), 805–807.

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

Ferraro, V.C., 1937. The non-uniform rotation of the sun and its magnetic field. In: Monthly Notices of the Royal Astronomical Society, vol. 97, pp. 458–497.

Finlay, C.C., Amit, H., 2011. On flow magnitude and field-flow alignment at Earth's core surface. Geophys. J. Int. 186 (1), 175–192.

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

Finlay, C.C., Kloss, C., Olsen, N., Hammer, M.D., Toffner-Clausen, L., Grayver, A., Kuvshinov, A., 2010b. The CHAOS-7 geomagnetic field model and observed changes in the South Atlantic Anomaly. Earth Planets Space 72. https://doi.org/10.1186/s40623-020-01252-9.

Finlay, C.C., Kloss, C., Olsen, N., Hammer, M.D., 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). https://doi.org/10.1186/s40623-020-01252-9.

Finlay, C.C., Olsen, N., Kotsiaros, S., Gillet, N., Toffner-Clausen, L., 2016. Recent geomagnetic secular variation from swarm and ground observatories as estimated in the CHAOS-6 geomagnetic field model. Earth Planets Space 68 (1), 112. ISSN 1880-5981. https://doi.org/10.1186/s40623-016-0486-1.

Finlay, C.C., Gillet, N., Aubert, J., Livermore, P.W., Jault, D., 2023. Gyres, jets and waves in the Earth's core. Nat. Rev. Earth Environ. 1–16.

Firsov, I., 2024. Geomagnetic View on the Flow and Radial Shear at the Top of the Earth's Core. Ph.D. thesis. Université Grenoble Alpes [2020–2024].

Firsov, I., Jault, D., Gillet, N., Aubert, J., Mandea, M., 2023. Radial shear in the flow at the Earth's core surface. Geophys. J. Int. 235 (3), 2524–2539.

Fournier, A., Aubert, J., Thébault, E., 2011. Inference on core surface flow from observations and 3-D dynamo modelling. Geophys. J. Int. 186 (1), 118–136.

Frasson, T., Schaeffer, N., Nataf, H.-C., Labrosse, S., 2024. Geomagnetic dipole stability and zonal flow changes controlled by mantle heat flux heterogeneities. Geophys. J. Int. 240 (3), 1481-1504. https://doi.org/10.1093/gji/ggae457.

Friedman, J., Hastie, T., Tibshirani, R., 2008. Sparse inverse covariance estimation with the graphical lasso. Biostatistics 9 (3), 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 (2), 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, 48020240184. https://doi.org/10.1098/rspa.2024.0184.

Gillet, N., Lesur, V., Olsen, N., 2010. Geomagnetic core field secular variation models. Space Sci. Rev. 155 (1–4), 129–145. ISSN 00386308. https://doi.org/10.1007/s11214-009-9586-6.

Gillet, N., Jault, D., Finlay, C.C., Olsen, N., 2013. Stochastic modeling of the Earth's magnetic field: inversion for covariances over the observatory era. Geochem. Geophys. Geosyst. 14 (4), 766–786. ISSN 1525-2027. https://doi.org/10.1002/ggg.e.20041.

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

Gillet, N., Huder, L., Aubert, J., 2019. A reduced stochastic model of core surface dynamics based on geodynamo simulations. Geophys. J. Int. 219 (1), 522–539. ISSN 0956-540X. https://doi.org/10.1093/gji/ggz313.

Gillet, N., Gerick, F., Jault, D., Schwaiger, T., Aubert, J., Istas, M., 2022. Satellite magnetic data reveal interannual waves in Earth's core. Proc. Natl. Acad. Sci. 119 (13), e2115258119.

Gillet, N., Dall'Asta, F., Amblard, P.-O., Claveau, J., Aubert, J., 2024. Waves in Earth's core and geomagnetic field forecast. Phys. Earth Planet. Inter. 357, 107284.

Gubbins, D., Willis, A.P., Sreenivasan, B., 2007. Correlation of Earth's magnetic field with lower mantle thermal and seismic structure. Phys. Earth Planet. Inter. 162 (3–4), 256–260. https://doi.org/10.1016/j.pepi.2007.04.014.

Hammer, M.D., Cox, G.A., Brown, W.J., Beggan, C.D., Finlay, C.C., 2021. Geomagnetic Virtual Observatories: monitoring geomagnetic secular variation with the Swarm satellites. Earth Planets Space 73 (1). https://doi.org/10.1186/s40623-021-01357-9.

Hollerbach, R., 1994. Magnetohydrodynamic Ekman and Stewartson layers in a rotating spherical shell. Proc. R. Soc. Lond. A 444 (1921), 333–346.

Holme, R., 2015. Large-Scale Flow in the Core, vol. 8. Elsevier B.V. https://doi.org/10.1016/B978-044452748-6.00127-9. ISBN 9780444527486.

Huder, L., Gillet, N., Thollard, F., 2019. pygeodyn 1.1.0: a Python package for geomagnetic data assimilation. Geosci. Model Dev. 12 (8), 3795–3803. ISSN 1991-9603. https://doi.org/10.5194/gmd-12-3795-2019.

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 (1), 160. https://doi.org/10.1186/s40623-020-01194-2.

Hulot, G., Olsen, N., Thébault, E., Hemant, K., 2009. Crustal concealing of small-scale core-field secular variation. Geophys. J. Int. 177 (2), 361–366. https://doi.org/10.1111/j.1365-246X.2009.04119.x.

Hulot, G., Sabaka, T., Olsen, N., Fournier, A., 2015. 5.02 - The present and future geomagnetic field. In: Schubert, G. (Ed.), Treatise on Geophysics, Second edition. Elsevier, Oxford, pp. 33–78. https://doi.org/10.1016/B978-0-444-53802-4.00096-8. ISBN 978-0-444-53803-1. URL: https://www.sciencedirect.com/science/article/pii/B9780444538024000968.
Istas, M., Gillet, N., Finlay, C.C., Hammer, M.D., Huder, L., 2023. Transient core surface dynamics from ground and satellite geomagnetic data. Geophys. J. Int. 233 (3), 1890-1915. ISSN 0956-540X. https://doi.org/10.1093/gji/ggad039.

IUGG AR23, 2024. Annual Report 2023, International Union of Geodesy and Geophysics, Potsdam, 135 p.

Jault, D., Finlay, C., 2015. Waves in the Core and Mechanical Core-Mantle Interactions. chap. 8.09., Elsevier, pp. 225-244. https://doi.org/10.1016/B978-0-444-53802-4.00150-0.

Kloss, C., Finlay, C.C., 2019. Time-dependent low-latitude core flow and geomagnetic field acceleration pulses. Geophys. J. Int. 217 (1), 140-168. ISSN 0956-540X. https://doi.org/10.1093/gji/ggy545.

Kutzner, C., Christensen, U.R., 2002. From stable dipolar towards reversing numerical dynamos. Phys. Earth Planet. Inter. 131 (1), 29-45.

Langel, R., Estes, R., 1982. A geomagnetic field spectrum. Geophys. Res. Lett. 9 (4), 250-253.

Lesur, V., Wardinski, I., Asari, S., Minchev, B., Mandea, M., 2010. Modelling the Earth's core magnetic field under flow constraints. Earth Planets Space 62 (6), 503-516. https://doi.org/10.5047/eps.2010.02.010.

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. https://doi.org/10.1007/s10712-021-09662-4.

Lézin, M., Amit, H., Terra-Nova, F., Wardinski, I., 2023. Mantle-driven north-south dichotomy in geomagnetic polar minima. Phys. Earth Planet. Inter. 337. https://doi.org/10.1016/j.pepi.2023.107000.

Lhuillier, F., Fournier, A., Hulot, G., Aubert, J., 2011. The geomagnetic secular variation timescale in observations and numerical dynamo models. Geophys. Res. Lett. 38 (9), 1-5. ISSN 00948276., https://doi.org/10.1029/2011GL047356.

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

Livermore, P., Hollerbach, R., Finlay, C., 2017. An accelerating high-latitude jet in Earth's core. Nat. Geosci. 10 (1), 62-68. https://doi.org/10.1038/ngeo2859. URL http://www.nature.com/doifinder/10.1038/ngeo2859.

Lowes, F.J., 1966. Mean-square values on sphere of spherical harmonic vector fields. J. Geophys. Res. 71 (8), 2179. https://doi.org/10.1029/JZ071i008p02179.

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 (1-4), 147-175. https://doi.org/10.1007/s11214-010-9663-x.

Nataf, H.-C., Schaeffer, N., 2015. 8.06 - Turbulence in the core. In: Schubert, G. (Ed.), Treatise on Geophysics, Second edition. Elsevier, Oxford, pp. 161-181. https://doi.org/10.1016/B978-0-444-53802-4.00142-1. ISBN 978-0-444-53803-1. URL https://www.sciencedirect.com/science/article/pii/B9780444538024001421.

Olsen, N., Lühr, H., Sabaka, T.J., Mandea, M., Rother, M., Töffner-Clausen, L., Choi, S., 2006. CHAOS-a model of the Earth's magnetic field derived from CHAMP, Ørsted, and SAC-C magnetic satellite data. Geophys. J. Int. 166 (1), 67-75. https://doi.org/10.1111/j.1365-246X.2006.02959.x.

Olsen, N., Lühr, H., Finlay, C.C., Sabaka, T.J., Michaelis, I., Rauberg, J., Töffner-Clausen, L., 2014. The CHAOS-4 geomagnetic field model. Geophys. J. Int. 197 (2), 815-827. https://doi.org/10.1093/gji/ggu033.

Olsen, N., Ravat, D., Finlay, C.C., Kother, L.K., 2017. LCS-1: a high-resolution global model of the lithospheric magnetic field derived from CHAMP and swarm satellite observations. Geophys. J. Int. 211 (3), 1461-1477. https://doi.org/10.1093/gji/ggx381.

Olson, P., Aurnou, J., 1999. A polar vortex in the Earth's core. Nature 402 (6758), 170-173. ISSN 1476-4687. https://doi.org/10.1038/46017.

Olson, P.L., Coe, R.S., Driscoll, P.E., Glatzmaier, G.A., Roberts, P.H., 2010. Geodynamo reversal frequency and heterogeneous core-mantle boundary heat flow. Phys. Earth Planet. Inter. 180 (1-2), 66-79.

Pais, M.A., Jault, D., 2008. Quasi-geostrophic flows responsible for the secular variation of the (E)arth's magnetic field. Geophys. J. Int. 173 (2), 421-443. ISSN 0956-540X. 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 (4), 2527-2540.

Rau, S., Christensen, U., Jackson, A., Wicht, J., 2000. Core flow inversion tested with numerical dynamo models. Geophys. J. Int. 141 (2), 485-497.

Roberts, P.H., Scott, S., 1965. On analysis of the secular variation, I. A hydromagnetic constraint: theory. J. Geomagn. Geoelectr. 17 (2), 137-151.

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

Ropp, G., Lesur, V., Baerenzung, J., Holschneider, M., 2020. Sequential modelling of the Earth's core magnetic field. Earth Planets Space 72 (1), 153. https://doi.org/10.1186/s40623-020-01230-1.

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), 157. https://doi.org/10.1186/s40623-020-01279-y.

Schaeffer, N., 2013. Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations. Geochem. Geophys. Geosyst. 14 (3), 751-758.

Schaeffer, N., Jault, D., Nataf, H.-C., Fournier, A., 2017. Turbulent geodynamo simulations: a leap towards Earth's core. Geophys. J. Int. 211 (1), 1-29. ISSN 0956-540X. https://doi.org/10.1093/gji/ggx265.

Schwaiger, T., Gastine, T., Aubert, J., 2020. Relating force balances and flow length scales in geodynamo simulations. Geophys. J. Int. 224 (3), 1890-1904. https://doi.org/10.1093/gji/ggaa545.

Schwaiger, T., Jault, D., Gillet, N., Schaeffer, N., Mandea, M., 2023. Local estimation of quasi-geostrophic flows in Earth's core. Geophys. J. Int. 234 (1), 494-511. ISSN 0956-540X. https://doi.org/10.1093/gji/ggad089.

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 (1), 107104. ISSN 00319201. https://doi.org/10.1016/j.pepi.2023.107104.

Shure, L., Parker, R.L., Backus, G.E., 1982. Harmonic splines for geomagnetic modelling. Phys. Earth Planet. Inter. 28 (3), 215-229. ISSN 00319201. https://doi.org/10.1016/0031-9201(82)90003-6.

Sprain, C.J., Biggin, A.J., Davies, C.J., Bono, R.K., Meduri, D.G., 2019. An assessment of long duration geodynamo simulations using new paleomagnetic modeling criteria (QPM). Earth Planet. Sci. Lett. 526, 115758. https://doi.org/10.1016/j.epsl.2019.115758.

Thébault, E., Hemant, K., Hulot, G., Olsen, N., 2009. On the geographical distribution of induced time-varying crustal magnetic fields. Geophys. Res. Lett. 36 (1).

Tobias, S.M., Cattaneo, F., Boldyrev, S., 2012. MHD Dynamos and Turbulence. In: Davidson, P.A., Kaneda, Y., Sreenivasan, K.R. (Eds.), *Ten Chapters in Turbulence*. Cambridge University Press, pp. 351-404.

Whaler, K.A., 1986. Geomagnetic evidence for fluid upwelling at the core-mantle boundary. Geophys. J. R. Astron. Soc. 86 (2), 563-588. ISSN 1365-246X. https://doi.org/10.1111/j.1365-246X.1986.tb03844.x.

Whaler, K., Hammer, M., Finlay, C., Olsen, N., 2022. Core surface flow changes associated with the 2017 Pacific geomagnetic jerk. Geophys. Res. Lett. 49 (15), e2022GL098616.

Wicht, J., Stellmach, S., Harder, H., 2010. Numerical dynamo simulations: from basic concepts to realistic models. In: Freeden, W., Nashed, M.Z., Sonar, T. (Eds.), Handbook of Geomathematics. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 459-502. https://doi.org/10.1007/978-3-642-01546-5 16.

Yokoi, N., 2013. Cross helicity and related dynamo. Geophys. Astrophys. Fluid Dyn. 107 (1-2), 114-184.