State and evolution of the geodynamo from numerical models reaching the physical conditions of Earth’s core

SUMMARY

Our understanding of the geodynamo has recently progressed thanks to geomagnetic data of improved quality, and analyses resting on numerical simulations of increasing realism. Here, these two advances are combined in order to diagnose the state and present dynamics of Earth’s core in physically realistic conditions. A sequential, ensemble-based framework assimilates the output of geomagnetic field models covering the past 180 yr into a numerical geodynamo simulation, the physical realism of which is also advanced as data is assimilated. The internal dynamical structure estimated for the geodynamo at present reproduces previously widely documented features such as a planetary-scale, eccentric westwards gyre and localization of buoyancy release beneath the Eastern (0°E−180°E) hemisphere. Relating the typical magnetic variation timescale of the assimilated states to the power at which they operate, the present convective power of the geodynamo is estimated at 2.95 ± 0.2 TW, corresponding to an adiabatic heat flow out of the core of 14.8 ± 1 TW if the top of the core is convectively neutrally stratified at present. For the first time, morphologically and dynamically relevant trajectories are obtained by integrating the estimated states forward for a few decades of physical time using a model reaching the physical conditions of Earth’s core. Such simulations accurately account for the spatio-temporal content of high-resolution satellite geomagnetic field models and confirm earlier interpretations in terms of rapid core dynamics. The enforcement of a realistic force balance approaching a Taylor state allows for propagation of weak (velocity perturbation of about 0.6 |$\mathrm{km\,yr^{-1}}$|⁠) axisymmetric torsional waves with period about 5 yr, supported by a magnetic field of root-mean-squared amplitude of 5.6 mT inside the core. Quasi-geostrophic magneto-Coriolis waves of interannual periods and significantly stronger velocity perturbation (about |$7 \mathrm{km\,yr^{-1}}$|⁠) are also reproduced, with properties that converge towards those recently retrieved from the analysis of geomagnetic variations before fully achieving Earth’s core conditions. The power spectral density of magnetic variations falls off rapidly at frequencies exceeding the inverse Alfvén time (about |$0.6\mathrm{yr^{-1}}$|⁠), which indicates that the excitation of hydromagnetic waves occurs preferentially at large spatial scales. The possibility to account for geomagnetic variations from years to centuries in physically realistic models opens the perspective of better constraining properties of the deep Earth through geomagnetic data assimilation.

Dynamo: theories and simulations, Rapid time variations, Satellite magnetics

1 INTRODUCTION

The geomagnetic signal documents the dynamics of Earth’s convecting fluid core, where it is generated through a self-sustained dynamo mechanism. It also contains important geophysical constraints on the structure and past evolution of the deep Earth (e.g. Landeau et al. 2022). In comparison with other geophysical signals of internal origin, the exploitation of geomagnetic data has traditionally been more complex because of the difficulty to elaborate a theory and numerical simulations able to reproduce the generation of this signal in physically realistic conditions. Yet, the rapid pace at which ground observatory and satellite geomagnetic data improve the core sounding for rapid dynamics (Lesur et al. 2022), thereby revealing new and previously unsuspected physical phenomena (Gerick et al. 2021; Gillet et al. 2022), holds the promise to learn significantly more about deep Earth and the geodynamo from the standpoint of observations, despite their intrinsic spatial and temporal limitations.

The progress of geomagnetic observations has largely fuelled the development of numerical geodynamo simulations (see e.g. Wicht & Sanchez 2019). Early in the development of this discipline, the output of such simulations has presented a large degree of morphological similarity with the geomagnetic field (Glatzmaier & Roberts 1995; Kuang & Bloxham 1997; Christensen et al. 1998). The success of the early models has raised the question of whether they were obtaining the right answers for the wrong physical reasons, because the conditions at which these models operated were very remote from those of the Earth’s core. The development of a standardized setup for simulations (Christensen et al. 2001) and systematic surveys have led to a more precise mapping of the accessible parameter space, delineating the regions where the morphology of the magnetic field output could be successfully compared to observations on a quantitative basis, and over a variety of timescales (Christensen et al. 2010; Davies & Constable 2014; Mound et al. 2015; Sprain et al. 2019). Over secular timescales typical of core convection, it appeared that morphologically compliant models could be obtained at two basic conditions: (i) the convective overturn time τ_U = D/U (where D is the core thickness and U a typical flow velocity) and the magnetic diffusion time τ_η = D²/η (where η is the magnetic diffusivity in the core) should be set in an Earth-like ratio Rm = τ_η/τ_U ≈ 1000 (the magnetic Reynolds number); (ii) the planetary rotation timescale τ_Ω = 1/Ω obtained from the rotation rate Ω should be sufficiently shorter than the core overturn time τ_U, that is, |$Ro=\tau _\Omega /\tau _{U}\le 10^{-2}$| (the Rossby number), such that the convective motions that sustain and shape the magnetic field should be placed under the rotational constraint set by the influence of the Coriolis force.

On this basis, it appeared plausible that a continuum of models could connect the physical conditions at which most classical models operate to those of the Earth’s core. This perspective later materialized (Aubert et al. 2017) with the formulation of a theoretical ‘path’ in parameter space that achieves this connection by following a line of constant Rm and decreasing Ro. An approximated numerical method was also supplied that enabled the exploration of this path up to conditions closely approaching those of Earth’s core (Aubert 2019; Aubert & Gillet 2021; Aubert et al. 2022). Along this path, it was shown that the large-scale (hundreds of kilometres) and long-period (secular) output of models remains invariant, thereby contributing to settle the debate on the realism of early models. At shorter (interannual to decadal timescales), however, the increasing magnetic and rotational constraints present in models advancing along the path give rise to rapid hydromagnetic waves that are absent from classical models, and that superimpose to the invariant and slower core convection (Aubert & Gillet 2021). The new quasi-geostrophic (QG) magneto-inertial (Alfvén) and magneto-Coriolis waves exhibited that way could successfully explain the occurrence of geomagnetic jerks, or abrupt changes in the acceleration of Earth’s magnetic field (Aubert & Finlay 2019; Aubert et al. 2022). A tight interconnection between numerical and observational research also enabled the first theoretical characterization (Gerick et al. 2021) and observation (Gillet et al. 2022; Istas et al. 2023) of rapid, QG magneto-Coriolis (QGMC) waves in core flows inferred from geomagnetic variations.

With the stumbling block of physical modelling realism being now largely removed, it is possible to envision the next steps for the geophysical exploitation of ground observatory and satellite geomagnetic data. This appears all the more promising that data assimilation, which has come to the fore in geomagnetism more than a decade ago (Kuang et al. 2009; Fournier et al. 2010), has reached some degree of maturity through the formulation of steadily improving frameworks (e.g. Aubert 2015, 2020) able to estimate the internal dynamical state of the core from geomagnetic surface data and statistical covariances linking the surface to the interior. Data assimilation has led to physics-based predictions of the future evolution of the geodynamo at the range of a few years (Kuang et al. 2010; Fournier et al. 2015, 2021a). The ability of such predictions to exceed traditional mathematical extrapolations however crucially rests on the possibility to correctly describe the rapid core dynamics and its magnetic perturbations that include geomagnetic jerks (Fournier et al. 2021b).

This study attempts to lay the groundwork for the assimilation of geomagnetic data into physically realistic numerical geodynamo simulations. This is a difficult task, because it imposes the requirement to simultaneously handle morphology and dynamics, that is, to reconcile an essentially statistical (and hence possibly out of dynamical equilibrium), data-based inference of the core internal state with the delicate dynamics that is enforced close to Earth’s core conditions. Our first goal is therefore to obtain an ensemble of states that satisfy this double constraint. Using such states as starting points, forecasts a few decades long in physical time are numerically achievable in a simulation operating (for the first time) at the physical conditions of Earth’s core regarding the relative position of all key timescales, which forms the second goal of this study. This opens the way towards our third goal, performing a detailed comparison between the signatures of rapid hydromagnetic waves and the observations in a context that now dispenses with the needs of scaling laws and extrapolation. In laying down a unified numerical model that accounts for both slow (secular) and rapid (interannual) geomagnetic variations in physically and morphologically realistic conditions, the next goal is to constrain deep Earth properties relevant to global geodynamics from the standpoint of various timescales. Finally, we wish to inform the future design of geomagnetic field models, as well as that of future satellite missions (e.g. Hulot et al. 2018), by constraining the spatial and temporal properties of rapid geomagnetic dynamics.

This manuscript is laid out as follows: Section 2 introduces the numerical geodynamo model, the data assimilation framework as well as the sources of data. The results are presented in Section 3, and discussed in Section 4.

2 MODELS AND METHODS

2.1 Numerical model

The dynamical model used here solves for Boussinesq convection, thermochemical density anomaly transport and magnetic induction in the magnetohydrodynamic approximation. The set of equations and numerical method can be found in Aubert et al. (2013, 2017). We search for the velocity field |$\boldsymbol{\mathrm{u}}$|⁠, magnetic field |$\boldsymbol{\mathrm{B}}$| and density anomaly field C, within an electrically conducting and rotating spherical shell of thickness D = r_o − r_i and aspect ratio r_i/r_o = 0.35 representing the Earth’s outer core. This shell is electromagnetically coupled to a solid inner core of radius r_i and electrical conductivity σ equal to that of the outer core, and a solid mantle between radii r_o and a = 1.83r_o (the surface of the Earth), featuring a thin conducting layer at its base of thickness Δ and conductivity σ_m such that Δσ_m/Dσ = 10⁻⁴. The inner core and mantle are furthermore gravitationally coupled, and can present axial differential rotations Ω_IC, Ω_M with respect to the planetary rotation frame rotating at a constant angular velocity |$\boldsymbol{\mathrm{\Omega }}=\Omega \boldsymbol{\mathrm{e}}_{z}$|⁠. The moments of inertia assigned to the inner core, outer core and mantle respect Earth-like proportions. Mechanically stress-free and electromagnetically conducting boundary conditions are imposed at r = r_i and r = r_o.

Convection is driven by a homogeneous mass anomaly flux F imposed at radius r = r_i, a volumetric buoyancy sink term in the volume to ensure mass conservation, and a vanishing mass anomaly flux imposed at r = r_o. Throughout this work, it is therefore assumed that the top of the core is neutrally stratified, that is, the heat flow out of the core is exactly adiabatic, a reasonable median hypothesis given the existing geophysical constraints (Davies et al. 2015; Frost et al. 2022). The modelling setup differs from Aubert et al. (2013, 2017) as we do not superimpose spatial heterogeneities on top of the homogeneous inner and outer mass anomaly fluxes. Here, we will rather let the possible heterogeneities arise from the constraints set by adherence to the geomagnetic data (Aubert 2014).

Within a spherical coordinate system (r, θ, φ) with unit vectors |$\boldsymbol{\mathrm{e}}_{r},\boldsymbol{\mathrm{e}}_{\theta },\boldsymbol{\mathrm{e}}_{\varphi }$|⁠, the numerical implementation involves a spectral spherical harmonic decomposition of |$\boldsymbol{\mathrm{u}}$|⁠, |$\boldsymbol{\mathrm{B}}$| and C up to degree and order ℓ_max = 133 and a discretization in the radial direction on a second-order finite-differencing scheme with NR gridpoints. Alternatively, we will also use a cylindrical system (s, z, φ) of unit vectors |$\boldsymbol{\mathrm{e}}_{s},\boldsymbol{\mathrm{e}}_{z},\boldsymbol{\mathrm{e}}_{\varphi }$|⁠. We use the spherical harmonics transform library SHTns (Schaeffer 2013) freely available at https://bitbucket.org/nschaeff/shtns. We also use a second-order, semi-implicit time stepping scheme. Solutions at extreme conditions are approximated by using a hyperdiffusive treatment applied to small scales of the velocity and density anomaly fields, but not the magnetic field which remains fully resolved. The functional form of hyperdiffusivity involves a cut-off ℓ_h = 30 below which hyperdiffusivity is not applied, and a parameter q_h such that the effective viscous and thermo-chemical diffusivities ν_eff, κ_eff smoothly increase over their molecular values ν, κ as |$(\nu _\text{eff},\kappa _\text{eff})=(\nu ,\kappa )\, q_{h}^{\ell -\ell _{h}}$| for spherical harmonic degrees ℓ ≥ ℓ_h. Full details on the physical justification of this approximation may be found in Aubert et al. (2017). Values of NR and q_h used in our models are reported in Table 1. All models produce a self-sustained, dipole-dominated magnetic field without polarity reversals during their integration. At the outer boundary of the model, the morphology of this field is deemed Earth-like, according to the standard rating χ² < 2 of free runs unconstrained by the data (see values in Table 1 and definition in Christensen et al. 2010).

Table 1.

Main parameters of the numerical models. All models have Pr = 1 and ℓ_max = 133. Earth’s core estimates are obtained from Pichon et al. (2016) and Aubert et al. (2017).

Label	Per cent of path	ϵ	RaF	Eη	Pm	ξ	χ2	NR	qh
Prior PA	0	1	2.7 × 10−5	7.5 ×  10−6	4	0.75	1.1	160
A1–A42	14	10−1	2.7 × 10−6	2.37 × 10−6	1.26	0.75		200	1.05
A1–A42	29	10−2	2.7 × 10−7	7.5 × 10−7	0.4	0.75		320	1.09
A1–A42	43	10−3	2.7 × 10−8	2.37 × 10−7	0.126	0.75		504	1.12
A1–A4	71	10−5	2.7 × 10−10	2.37 × 10−8	1.26 × 10−2	0.75		1280	1.19
A1–A4	100	10−7	2.7 × 10−12	2.37 × 10−9	1.26 × 10−3	0.75		2496	1.25
Prior PB	0	1	1.68 ×  10−5	7.5 × 10−6	4	2.1	0.72	160
B1–B42	14	10−1	1.68 × 10−6	2.37 × 10−6	1.26	0.66		200	1.05
B1–B42	29	10−2	1.68 × 10−7	7.5 × 10−7	0.4	0.21		320	1.09
B1–B42	43	10−3	1.68 × 10−8	2.37 × 10−7	0.126	6.6 × 10−2		504	1.12
B1–B4	71	10−5	1.68 × 10−10	2.37 × 10−8	1.26 × 10−2	6.6 × 10−3		1280	1.19
B1–B4	100	10−7	1.68 × 10−12	2.37 × 10−9	1.26 × 10−3	6.6 × 10−4		2496	1.25
Earth			|$\mathcal {O}(10^{-12})$|	|$\mathcal {O}(10^{-9})$|	|$\mathcal {O}(10^{-6})$|	|$\mathcal {O}(10^{-4})$|
A1–A4	114	10−8	2.7 × 10−13	7.5 × 10−10	4 × 10−4	0.75		3192	1.40
B1–B4	114	10−8	1.68 × 10−13	7.5 × 10−10	4 × 10−4	2.1 × 10−4		3192	1.40

2.2 Model dimensionless parameters and path theory

The four main dimensionless parameters are the flux-based Rayleigh, magnetic Ekman, Prandtl and magnetic Prandtl numbers:

$$\begin{eqnarray} Ra_{F}=\dfrac{g_{o}F}{4\pi \rho \Omega ^{3}D^{4}},~ E_{\eta }=\dfrac{\eta }{\Omega D^{2}},~ Pr=\dfrac{\nu }{\kappa },~ Pm=\dfrac{\nu }{\eta }. \end{eqnarray}$$

(1)

Here g_o is gravity at the core surface, ρ is the outer core density, the magnetic diffusivity η = 1/μσ relates to the conductivity σ through the magnetic permeability μ, and we recall that ν, κ are the molecular viscous and thermo-chemical diffusivities. The gravitational torque exerted on the inner core by the mantle writes Γ_G = −Γτ(Ω_IC − Ω_M), with Γ being the gravitational coupling strength and τ the viscous relaxation time of the solid inner core (Pichon et al. 2016). This therefore involves a fifth dimensionless parameter ξ = Γτ/ρD⁵Ω.

This study uses the path theory (Aubert et al. 2017) that connects the parameter space region where classical models computed at moderate conditions are found to the region with extreme parameters where Earth’s core resides. Starting from a parent model at the beginning of a path (hereafter labelled as 0 per cent of this path), a sequence of child models can be built that regularly advances along this path. The child models preserve the large-scale (spherical harmonic degree below 30) and long-timescale (secular and longer) output of the parent model, while enriching their shorter-timescale content with the rapid hydromagnetic wave dynamics that is expected to take place in Earth’s core (Aubert & Gillet 2021; Aubert et al. 2022). Labelling the parent model parameters with ‘(0)’, any position along the path may be represented by a path parameter ϵ, or alternatively a percentage of path −100 log₁₀(ϵ)/7, with the main parameters of the child model being:

$$\begin{eqnarray} Ra_{F}(\epsilon )=\epsilon Ra_{F} (0),~ E_{\eta }(\epsilon )=\sqrt{\epsilon } E_{\eta } (0),\\ Pr(\epsilon )=1,~ Pm(\epsilon )=\sqrt{\epsilon } Pm (0). \end{eqnarray}$$

(2)

The start of path corresponds to ϵ = 1, and Aubert et al. (2017) have shown that Earth’s core conditions are adequately described by the end of path (100 per cent) obtained at ϵ = 10⁻⁷. Table 1 lists the path positions and parameters used in the parent and child models of this study and confirms that Ra_F and E_η indeed reach Earth-like values at the end of the path. Despite being the smallest achieved to date in a a numerical dynamo simulation, the magnetic Prandtl number Pm however does not reach the value Pm ≈ 10⁻⁶ that is usually quoted for Earth’s core. As in Aubert et al. (2017), we state again that the end of path nevertheless ensures that the magnetic diffusion time τ_η is far shorter than the viscous and thermochemical diffusion times τ_ν = D²/ν and τ_κ = D²/κ, such that the two latter associated diffusivities are no longer relevant to the dynamics. Considering the fifth dimensionless number ξ describing the strength of gravitational coupling, Aubert et al. (2017) and subsequent studies have adopted a constant ξ = ξ(0) along the path, a choice that we will also adopt in one sequence of parent and child models (sequence A). Table 1 however shows that this choice overshoots the expected gravitational coupling at Earth’s core conditions, and we shall see that it is also not optimal for long-timescale invariance of the solution along the path. Along the path, it appears more relevant (Aubert et al. 2022) to preserve the ratio of the core overturn time τ_U = D/U with the gravitational restoring time of the mantle I_M/Γτ, with I_M being the moment of inertia of the mantle. This implies that ξ = Γτ/ρD⁵Ω should scale like the Rossby number Ro = τ_Ω/τ_U of the simulation. Using the path scaling law |$Ro\sim \sqrt{\epsilon }$| given in Aubert et al. (2017), this choice then prompts |$\xi =\sqrt{\epsilon }~\xi (0)$|⁠, and ensures an Earth-like value for ξ at the end of the path (Table 1). Our second sequence B therefore implements this dependency. Aside from this, the main difference between sequences A and B is the flux Rayleigh number Ra_F at which they operate, which is lower in sequence B.

2.3 Dimensioning of model quantities

The numerical dynamo model solves for dimensionless quantities, which need to be dimensioned in order to compare these with the geodynamo and to perform data assimilation. For the first time, we introduce models situated at the end of the parameter space path, where this operation is greatly facilitated. Indeed, once standard values |$\rho =11000\mathrm{kg\,m^{-3}}$|⁠, |$D=2260\mathrm{km}$| and |$g_{o}=10\mathrm{m\,s^{-2}}$| are adopted for the density, thickness of the core and for the gravity at the core surface, the sole specification of the magnetic diffusion time |$\tau_\eta=D^2/\eta$| then enables a unique dimensional rescaling of all quantities. The dimensionless time expressed in magnetic diffusion time units of D²/η can be dimensioned straightforwardly. The same goes for the dimensionless velocity |$|\boldsymbol{\mathrm{u}}|$| once it is expressed in magnetic Reynolds units of η/D. The rotation rate Ω is obtained from the magnetic Ekman number through Ω = (E_ητ_η)⁻¹. The dimensionless magnetic field |$|\boldsymbol{\mathrm{B}}|$| expressed in Elsasser units of |$\sqrt{\rho \mu \eta \Omega }$| is then also straightforwardly dimensioned. Expressing the dimensionless density anomaly C in units of ρΩη/g_oD finally also leads to a straightforward dimensioning.

For a model situated before the end of the path, the same dimensioning scheme can be applied. Though this yields geophysically consistent values for distances, times and velocities, this will however severely underestimate the magnetic and density anomaly field amplitudes, because the rotation rate of these models is slower than that of Earth. This poses a problem to data assimilation, as the real and synthetic magnetic field amplitudes then cannot be directly compared. A way to handle this problem (Aubert 2020) is to make use of the invariance of the leading-order QG, magneto-Archimedes-Coriolis (MAC) force balance (Aubert et al. 2017; Schwaiger et al. 2019) between pressure, Coriolis, magnetic and buoyancy forces along the path, and to propose a correction to the dimensional rescaling of |$|\boldsymbol{\mathrm{B}}|$| and C that preserves this force balance once the quantities are cast into the dimensional world. This amounts to expressing their respective dimensionless values in units of |$\sqrt{\rho \mu \eta \Omega }$| and ρΩη/g_oD (as above), but then to multiplying the results with |$\sqrt{\rho \mu \eta \Omega _\mathrm{EOP}}$| and ρΩ_EOPη/g_oD, with the rotation rate Ω_EOP at the end of the path being used here in place of the actual rotation rate Ω. Stating this in terms of the path parameter ϵ (with ϵ = 10⁻⁷ representing the end of path), we recall that Ω = (E_ητ_η)⁻¹ with constant τ_η, and invoke the path rule |$E_{\eta }\sim \sqrt{\epsilon }$| (eq. 2). The dimensional values of |$|\boldsymbol{\mathrm{B}}|$| and C obtained through the standard procedure should then be corrected by multiplying them with |$\sqrt{\Omega _\mathrm{EOP}/\Omega }=(10^{7}\epsilon )^{1/4}$|⁠, and Ω_EOP/Ω = (10⁷ϵ)^1/2, respectively. Because the same dimensional force balance is then consistently respected, the dimensional amplitudes |$|\boldsymbol{\mathrm{u}}|$|⁠, |$|\boldsymbol{\mathrm{B}}|$| and C are in principle invariant along the path, and therefore match their geophysically relevant end-of-path values. In the following, and unless otherwise stated, we will generally use this correction that enables data assimilation at any path position.

2.4 Geomagnetic data assimilation framework

The basis of our data assimilation framework is the ensemble-based scheme introduced in Aubert (2015). Here, this scheme is simply modified in order to allow for sequential analyses at multiple epochs, as well as for advancing the position of the underlying dynamical model along the path as the assimilation proceeds. The numerical implementation uses a representation of |$\boldsymbol{\mathrm{u}}$| and |$\boldsymbol{\mathrm{B}}$| with toroidal and poloidal scalars:

$$\begin{eqnarray} \boldsymbol{\mathrm{u}}=\nabla \times (T \boldsymbol{\mathrm{r}})+\nabla \times \nabla \times (P \boldsymbol{\mathrm{r}}), \end{eqnarray}$$

(3)

$$\begin{eqnarray} \boldsymbol{\mathrm{B}}=\nabla \times (Z \boldsymbol{\mathrm{r}})+\nabla \times \nabla \times (W \boldsymbol{\mathrm{r}}). \end{eqnarray}$$

(4)

Here |$\boldsymbol{\mathrm{r}}=r\cdot \boldsymbol{\mathrm{e}}_{r}$| is the radius vector. The scalar fields P, T, W, Z, C representing the velocity, magnetic and density anomaly fields are expanded on a discrete radial grid r_j, j = 1,..., NR, and in a spherical harmonics over a fully normalized basis of complex functions |$Y_{\ell }^{m} (\theta ,\varphi )$| of degree ℓ and order m according to :

$$\begin{eqnarray} P = \sum _{\ell =0}^{\ell _\mathrm{max}} \sum _{m=-\ell }^{\ell } P_{\ell }^{m} (r_j)~Y_{\ell }^{m} (\theta ,\varphi ), \end{eqnarray}$$

(5)

in the case of, for example, the poloidal velocity scalar P (see Table 1 for values of NR and the maximal degree ℓ_max).

2.4.1 Geomagnetic data

We use as input data the output of geomagnetic field model COV-OBS (Gillet et al. 2013) between epochs 1840 and 1990, and CHAOS-7.9 (Finlay et al. 2020) at epoch 2000. This output is used up to spherical harmonic degree ℓ_obs = 13, which can safely be assumed to represent the field of core origin. Geomagnetic field models are sampled at discrete epochs for the poloidal magnetic field at the core surface |$W_{\ell }^{m}(r_{o})$| (the main field or MF data), and its first time derivative evaluated on a spline basis |$\dot{W}_{\ell }^m(r_{o})$| (the secular variation or SV data), which relate to the supplied classical Gauss coefficients |$g_{\ell }^{m},h_{\ell }^{m}$| through :

$$\begin{eqnarray} W_{\ell }^m(r_{o})=\dfrac{r_{o}}{\ell }\left(\dfrac{a}{r_{o}}\right)^{\ell +2}\dfrac{\left(g_{\ell }^{m}-i h_{\ell }^{m}\right)}{\sqrt{2\ell +1}}. \end{eqnarray}$$

(6)

Here, we recall that |$a=6371.2 \mathrm{km}$| is Earth’s radius, and a similar expression holds for |$\dot{W}_{\ell }^{m}(r_{o})$|⁠. The MF and SV data vectors

$$\begin{eqnarray} \boldsymbol{\mathrm{d}}=\left[W_{\ell }^m(r_{o})\right]^\intercal , 1\le \ell \le \ell _\mathrm{obs}, \end{eqnarray}$$

(7)

$$\begin{eqnarray} \dot{\boldsymbol{\mathrm{d}}}=\left[\dot{W}_{\ell }^m(r_{o})\right]^\intercal , 1\le \ell \le \ell _\mathrm{obs}, \end{eqnarray}$$

(8)

where superscript ^⊺ denotes the transpose, are introduced in the assimilation framework at epochs of analysis. Also needed are the data error covariance matrices |$\boldsymbol{\mathrm{R}}$| and |$\dot{\boldsymbol{\mathrm{R}}}$| statistically representing the deviation |$\boldsymbol{\mathrm{e}},\dot{\boldsymbol{\mathrm{e}}}$| of |$\boldsymbol{\mathrm{d}},\dot{\boldsymbol{\mathrm{d}}}$| to their true value in Earth’s core:

$$\begin{eqnarray} \boldsymbol{\mathrm{R}}=E(\boldsymbol{\mathrm{e}}\cdot \boldsymbol{\mathrm{e}^{\prime }})~~,~~ \dot{\boldsymbol{\mathrm{R}}}=E(\dot{\boldsymbol{\mathrm{e}}}\cdot \dot{\boldsymbol{\mathrm{e}}}^{\prime }), \end{eqnarray}$$

(9)

where E is the expected value and the prime denotes the transpose complex conjugate. For the MF data error covariance matrix |$\boldsymbol{\mathrm{R}}$|⁠, we use the posterior matrix supplied with COV-OBS between epochs 1840 and 1990. Such a matrix is however not supplied with CHAOS-7.9. At epoch 2000 we therefore specify a diagonal matrix with the following coefficients:

$$\begin{eqnarray} R_{\ell }^m=\dfrac{r_{o}^{2}}{\ell ^{2}(\ell +1)^{2}(2\ell +1)}\left(\dfrac{a}{r_{o}}\right)^{2\ell +4} \dfrac{e_{B}^{2}}{\ell _\mathrm{obs}}. \end{eqnarray}$$

(10)

Here |$e_{B}=1\mathrm{nT}$| is the error level, the energy of which is evenly distributed at the Earth surface among spherical harmonic degrees 1–13.

The SV data error covariance matrix |$\dot{\boldsymbol{\mathrm{R}}}$| is also obtained from the posterior matrix supplied with COV-OBS for the first analyses that use an underlying dynamical model situated at the start of the path. Subsequent analyses use a dynamical model that advances along the path and progressively enforce a tighter adherence to the QG-MAC force balance (Aubert et al. 2017). Similarly to Aubert (2020), using covariance matrices as tight as those supplied by COV-OBS while enforcing these constraints results in overfitting, that is, producing statistical estimates with overenergetic subsequent dynamics. For this reason, we prescribe a diagonal data error covariance matrix with similar form as (10), that is,

$$\begin{eqnarray} \dot{R}_{\ell }^m=\dfrac{r_{o}^{2}}{\ell ^{2}(\ell +1)^{2}(2\ell +1)}\left(\dfrac{a}{r_{o}}\right)^{2\ell +4} \dfrac{\dot{e}_{B}^{2}}{\ell _\mathrm{obs}}. \end{eqnarray}$$

(11)

Here, we choose |$\dot{e}_{B}=20 \mathrm{nT\,yr^{-1}}$|⁠, an error level that relaxes the fit to the data in order to reach an acceptable trade-off between this fit and the progressively tighter adherence to the QG-MAC balance.

2.4.2 Prior model and free run

Embedded into the data assimilation framework is an inverse modelling scheme that estimates the internal fields of the model from the surface data. This requires statistical knowledge of the correlations between the surface and the depth, that are encoded into model covariance matrices. Their computation is performed through a free run (unconstrained by the data) of our prior models PA and PB, situated at the start of the parameter space path (Table 1). During these runs, N = 1900 samples of the model state are collected, over the entire radial grid r_j, 1 ≤ j ≤ NR, and up to spherical harmonic degree ℓ_asm = 60. The time spacing being two samples is roughly an overturn time τ_U, such that samples may be considered statistically decorrelated. The magnetic and hydrodynamic state vectors are defined as

$$\begin{eqnarray} \boldsymbol{\mathrm{b}}=[W_\ell ^m(r_{j}),Z_\ell ^m(r_{j})]^\intercal , 1\le j\le NR~,~1\le \ell \le \ell _{\mathrm{asm}}, \end{eqnarray}$$

(12)

$$\begin{eqnarray} \boldsymbol{\mathrm{x}}=[S_\ell ^m(r_{j}), T_\ell ^m(r_{j}), C_\ell ^m(r_{j})]^\intercal , 1\le j\le NR~,~0\le \ell \le \ell _{\mathrm{asm}}. \end{eqnarray}$$

(13)

For practical reasons of numerical implementation, the spheroidal flow potential |$S_\ell ^m=\mathrm{d}(r P_\ell ^m)/r\mathrm{d} r$| is preferred over the equivalent poloidal flow potential |$P_\ell ^m$|⁠. The average over the ensemble of samples 1 ≤ i ≤ N is defined as, for example, |$\overline{\boldsymbol{\mathrm{x}}}=\sum _{i=1}^{N} \boldsymbol{\mathrm{x}}_{i}/N$|⁠, and the covariance matrices |$\boldsymbol{\mathrm{P}},\boldsymbol{\mathrm{Q}}$| for |$\boldsymbol{\mathrm{b}},\boldsymbol{\mathrm{x}}$| are approximated by

$$\begin{eqnarray} \boldsymbol{\mathrm{P}}\approx \dfrac{1}{N-1}\sum _{i=1}^{N} \left((\boldsymbol{\mathrm{b}}_{i}-\overline{\boldsymbol{\mathrm{b}}})(\boldsymbol{\mathrm{b}}_{i}-\overline{\boldsymbol{\mathrm{b}}})^{\prime }\right), \end{eqnarray}$$

(14)

$$\begin{eqnarray} \boldsymbol{\mathrm{Q}}\approx \dfrac{1}{N-1}\sum _{i=1}^{N} \left((\boldsymbol{\mathrm{x}}_{i}-\overline{\boldsymbol{\mathrm{x}}})(\boldsymbol{\mathrm{x}}_{i}-\overline{\boldsymbol{\mathrm{x}}})^{\prime }\right). \end{eqnarray}$$

(15)

2.4.3 Algorithm

The assimilation proceeds in an alternation of forecasts and analyses. Starting from a previously analysed state |$\boldsymbol{\mathrm{b}}_{i}^{a}$|⁠, |$\boldsymbol{\mathrm{x}}_{i}^{a}$|⁠, the forecast step advances an ensemble of states to obtain the forecast vectors |$\boldsymbol{\mathrm{b}}_{i}^{f}$|⁠, |$\boldsymbol{\mathrm{x}}_{i}^{f}$|⁠. At discrete epochs where data |$\boldsymbol{\mathrm{d}},\dot{\boldsymbol{\mathrm{d}}}$| are available, the analysis step corrects the states with an estimation derived from the data to obtain the newly analysed states |$\boldsymbol{\mathrm{b}}_{i}^{a}$|⁠, |$\boldsymbol{\mathrm{x}}_{i}^{a}$|⁠.

The analysis of the magnetic state vector |$\boldsymbol{\mathrm{b}}$| with the MF data |$\boldsymbol{\mathrm{d}}$| is purely statistical. We introduce the observation operator |$\boldsymbol{\mathrm{H}}$|⁠, a matrix with a simple form consisting of ones in entries corresponding to the observed poloidal magnetic field at the core surface, and zeros otherwise. The direct problem then writes |$\boldsymbol{\mathrm{H}}\boldsymbol{\mathrm{b}}=\boldsymbol{\mathrm{d}}$|⁠. As in Aubert (2015, 2020), the ensemble of analysed (corrected) solutions |$\boldsymbol{\mathrm{b}}^{a}_{i}$| is obtained from the forecasts |$\boldsymbol{\mathrm{b}}^{a}_{f}$| through classical Bayesian inference:

$$\begin{eqnarray} \boldsymbol{\mathrm{b}}^{a}_{i}=\boldsymbol{\mathrm{b}}_{i}^{f}+\dfrac{\boldsymbol{\mathrm{P}}\boldsymbol{\mathrm{H}^{\prime }}}{\boldsymbol{\mathrm{H}}\boldsymbol{\mathrm{P}}\boldsymbol{\mathrm{H}^{\prime }}+\boldsymbol{\mathrm{R}}} \left(\boldsymbol{\mathrm{d}}-\boldsymbol{\mathrm{H}}\boldsymbol{\mathrm{b}}_{i}^{f}\right). \end{eqnarray}$$

(16)

The ensemble |$\boldsymbol{\mathrm{b}}^{a}_{i}$| estimates the variability of the magnetic field hidden inside the core (including the unobservable toroidal field and small-scale components), given the surface poloidal observations |$\boldsymbol{\mathrm{d}}$|⁠, their error level |$\boldsymbol{\mathrm{R}}$|⁠, the previous model state |$\boldsymbol{\mathrm{b}}_{i}^{f}$| and the covariance properties |$\boldsymbol{\mathrm{P}}$| of the underlying dynamical model. We then compute an ensemble of estimates |$\boldsymbol{\mathrm{\delta }}_{i}^{a}$| for magnetic diffusion at the core surface through

$$\begin{eqnarray} \boldsymbol{\mathrm{\delta }}_{i}^{a}=\eta \left[(\nabla ^2 W^{a}_{i})_\ell ^{m}\right]^\intercal , r=r_{o}, 1\le \ell \le \ell _\mathrm{asm}, \end{eqnarray}$$

(17)

where |$W_{i}^a$| is the poloidal magnetic potential embedded in |$\boldsymbol{\mathrm{b}}^{a}_{i}$|⁠.

As in Aubert (2015, 2020) again, the analysis of the hydrodynamic state vector with the SV data |$\dot{\boldsymbol{\mathrm{d}}}$| combines statistical estimates with an inverse resolution of the magnetic induction equation at the core surface:

$$\begin{eqnarray} \displaystyle \dfrac{\displaystyle \partial \boldsymbol{\mathrm{B}}}{\displaystyle \partial t}=\nabla \times \left(\boldsymbol{\mathrm{u}}\times \boldsymbol{\mathrm{B}}\right)+\eta \nabla ^2\boldsymbol{\mathrm{B}}. \end{eqnarray}$$

(18)

We introduce the matrices |$\boldsymbol{\mathrm{M}}_{i}^\mathrm{a}$| that encode the spectral representations of the magnetic induction terms |$\nabla \times (\boldsymbol{\mathrm{u}}\times \boldsymbol{\mathrm{B}}_{i}^{a})$| from eq. (18), taken at the core surface (r = r_o). Their coefficients involve the analysed magnetic field |$\boldsymbol{\mathrm{b}}^{a}_{i}(r_{o})$| up to degree ℓ_asm, together with Elsasser and Adams–Gaunt coupling integrals (see Aubert 2013). The direct problem corresponding to eq. (18) then writes |$\boldsymbol{\mathrm{M}}_{i} \boldsymbol{\mathrm{x}}=\dot{\boldsymbol{\mathrm{d}}} - \boldsymbol{\mathrm{\delta }}_{i}^{a}$|⁠, and the analysed solution again obtained by Bayesian inference then writes

$$\begin{eqnarray} \boldsymbol{\mathrm{x}}^{a}_{i}=\left\langle \boldsymbol{\mathrm{x}}_{i}^{f}\right\rangle +\dfrac{\boldsymbol{\mathrm{Q}}\boldsymbol{\mathrm{M}}^{a^{\prime }}_{i}}{\boldsymbol{\mathrm{M}}_{i}^{a}\boldsymbol{\mathrm{Q}}\boldsymbol{\mathrm{M}}^{a^{\prime }}_{i} +\dot{\boldsymbol{\mathrm{R}}}} \left(\dot{\boldsymbol{\mathrm{d}}} - \boldsymbol{\mathrm{\delta }}_{i}^{a} -\boldsymbol{\mathrm{M}}_{i}^{a}\left\langle \boldsymbol{\mathrm{x}}_{i}^{f}\right\rangle \right). \end{eqnarray}$$

(19)

Here, the analysis differs from a classical sequential procedure such as that in eq. (16) in the sense that it is based on a time average (denoted by the angle brackets) of the hydrodynamic state vector over the duration of the previous forecast run. This modification is justified in Section 2.4.4.

To monitor the progression of the assimilation, after each analysis we report on the classical diagnostic quantities previously introduced in Aubert (2020). The normalized prior misfits to the MF and SV data are defined as

$$\begin{eqnarray} \Delta _\mathrm{MF}(i) = \sqrt{\dfrac{1}{N_{13}} \left(\boldsymbol{\mathrm{d}}-\boldsymbol{\mathrm{H}}\boldsymbol{\mathrm{b}}_{i}^{a}\right)^{\prime } \boldsymbol{\mathrm{R}}^{-1}\left(\boldsymbol{\mathrm{d}}-\boldsymbol{\mathrm{H}}\boldsymbol{\mathrm{b}}_{i}^{a}\right)}, \end{eqnarray}$$

(20)

$$\begin{eqnarray} \Delta _\mathrm{SV}(i) = \sqrt{\dfrac{1}{N_{13}} \left(\dot{\boldsymbol{\mathrm{d}}}-\boldsymbol{\mathrm{\delta }}_{i}^{a}-\boldsymbol{\mathrm{M}}_{i}^{a}\boldsymbol{\mathrm{x}}_{i}^{a}\right)^{\prime } \dot{\boldsymbol{\mathrm{R}}}^{-1}\left(\dot{\boldsymbol{\mathrm{d}}}-\boldsymbol{\mathrm{\delta }}_{i}^{a}-\boldsymbol{\mathrm{M}}_{i}^{a}\boldsymbol{\mathrm{x}}_{i}^{a}\right)}. \end{eqnarray}$$

(21)

Here, N₁₃ = 105 is the number of complex spherical harmonic coefficients involved in a decomposition up to degree and order ℓ_obs = 13. The normalized posterior increments brought by the analysis to the core surface magnetic and velocity fields are defined by

$$\begin{eqnarray} \Delta _{B}(i)=\sqrt{\dfrac{1}{N_{60}} \left(\boldsymbol{\mathrm{b}}_{i}^{a}-\boldsymbol{\mathrm{b}}_{i}^{f}\right)_{B(r_{o})}^{\prime } \boldsymbol{\mathrm{P}}_{B(r_{o})}^{-1}\left(\boldsymbol{\mathrm{b}}_{i}^{a}-\boldsymbol{\mathrm{b}}_{i}^{f}\right)_{B(r_{o})}}, \end{eqnarray}$$

(22)

$$\begin{eqnarray} \Delta _{u}(i)=\sqrt{\dfrac{1}{2 N_{60}} \left(\boldsymbol{\mathrm{x}}_{i}^{a}-\left\langle \boldsymbol{\mathrm{x}}_{i}^{f}\right\rangle \right)_{u(r_{o})}^{\prime } \boldsymbol{\mathrm{Q}}^{-1}_{u(r_{o})}\left(\boldsymbol{\mathrm{x}}_{i}^{a}-\left\langle \boldsymbol{\mathrm{x}}_{i}^{f}\right\rangle \right)_{u(r_{o})}}. \end{eqnarray}$$

(23)

Here, N₆₀ = 1891 is the number of complex spherical harmonic coefficients involved in a decomposition up to degree and order ℓ_asm = 60, and the subscripts B(r_o), u(r_o) indicate that the computation is restricted to the core surface state vector and covariance matrix components.

2.4.4 Changes of path position during assimilation

As the assimilation proceeds, the main idea of this study is to also advance the physical realism, or position of the underlying dynamical model along the path, at specific points in time. The theoretical scaling laws (Aubert et al. 2017) that should be applied to dimensionless fields when changing path positions are embedded by design into our dimensioning scheme (Section 2.3) that preserves the QG-MAC force balance. Formulating the algorithm in the dimensional world therefore dispenses with the need to apply these laws. By virtue of the invariance of the large-scale structure of models along the path (Aubert et al. 2017), the same prior covariance matrices |$\boldsymbol{\mathrm{P}}$| and |$\boldsymbol{\mathrm{Q}}$| can also be used at all path positions. Models at various path positions however involve different radial resolutions NR (Table 1). We therefore proceed by interpolating the forecast to the radial resolution at the start of the path, computing the correction at this resolution and interpolating this correction back to the actual radial resolution of the current path position. Spherical harmonic degrees ℓ_asm + 1 to ℓ_max are not analysed and therefore carried over through analysis without correction.

After a change of path position, preserving the QG-MAC force balance also implies that convective transients are minimized. This beneficial property has been used in our previous studies to initiate the computation of advanced models (e.g. Aubert 2019; Aubert et al. 2022). As the underlying model advances along the path, inertial forces coming at the next order also decrease relatively to the QG-MAC forces (Aubert et al. 2017). This gradually enforces a tighter adherence to the QG-MAC force balance in the full volume, a progress over Aubert (2020) where this balance could only be enforced at the surface. Still, with decreasing levels of inertial forces, a precise dynamical equilibration becomes increasingly at odds with the statistical corrections brought by analyses. The system mostly responds by exciting rapid inertial wave transients of timescales commensurate with τ_Ω, which are otherwise subdominant in properly equilibrated dynamics (Aubert & Gillet 2021). In order to mitigate this effect, in eq. (19) we base our analyses on the time average |$\left\langle \boldsymbol{\mathrm{x}_{i}^{f}}\right\rangle $| of the hydrodynamic state vector over the duration of the previous forecast rather than its actual value. Because such waves typically have smaller magnetic than kinetic energy (e.g. Gerick et al. 2021), a similar correction is not needed for the magnetic state vector |$\boldsymbol{\mathrm{b}}$|⁠, the analysis of which is based on the actual value at the forecast end in eq. (16).

2.5 Assimilation sequences and dynamically balanced final forecast runs

Each of our sequences A and B consist in an ensemble of trajectories Ai, Bi, 1 ≤ i ≤ 42 where geomagnetic data are assimilated in the epoch range 1840–2000. The assimilation starts with an analysis at epoch 1840, where individual members of the free run ensemble are used as forecast vectors to initiate the algorithm (16) and (19). The states are advanced using an underlying model situated at the start of the path in the epoch range 1840–1900, with analyses being performed every 20 yr in this range. The underlying model is then advanced to 14 per cent in the range 1900–1930, 29 per cent in the range 1940–1960 and 43 per cent in the range 1970–2000. Analyses are performed every 10 yr in the range 1900–2000.

Once 42 final analysed states are obtained at epoch 2000 for each sequence A, B, we terminate the trajectories by computing 80-yr long forecasts until epoch 2080 at 43 per cent of the path, with no further analysis. A subset of 8 final states (A1–A4 and B1–B4) at epoch 2000 are also advanced until 2002 using an underlying model jumping up to 114 per cent of the path, that is, beyond its end (Table 1). Using the 2000–2002 time-averaged fields as starting points in 2001 removes a large part of the imbalances triggering inertial transients, at the minor cost of losing a bit of adherence to the geomagnetic data. These states, which we hereafter refer to as ‘dynamically balanced’, are used to initiate forecasts until 2081 carried out at 43, 71 and 100 per cent of the path.

2.6 Dimensional inputs and outputs

Table 2 lists a number of dimensional quantities for free prior runs and final forecast runs at 43 and 100 per cent of the path. Aside from the magnetic diffusion time τ_η set as a fixed parameter for each sequence, the corresponding core conductivity σ = 1/μη, planetary rotation period 2πτ_Ω = 2π/Ω and the core overturn time τ_U = D/U (where U are the rms velocity inside the shell) are reported. The Alfvén timescale |$\tau _{A}=\sqrt{\rho \mu }D/B$| (where B is the rms magnetic field inside the shell) is presented without applying the dimensioning correction discussed in Section 2.3, in order to keep in mind that the typical velocity D/τ_A of rapid hydromagnetic waves is variable along the path and reaches a geophysically relevant value only at its end. With the dimensioning correction this time being applied, Table 2 also reports on the total convective power |$p=\int _{V} (\boldsymbol{\mathrm{u}}\cdot \boldsymbol{\mathrm{e}}_{r})C g_{o} r /r_{o}\, \mathrm{d}V$|⁠, where V is the volume of the outer core, the rms field intensity B inside the shell and its intensity B_o contained in spherical harmonic degrees below 12 at core surface.

Table 2.

Dimensional inputs and outputs for the prior free runs and final forecast runs, in the two sequences A, B (see the main text for definitions). The reported quantities are averages over the reported time span and ensemble members. Earth’s core estimates are obtained from Gillet et al. (2010), Lhuillier et al. (2011), Davies et al. (2015) and the CHAOS-7.9 geomagnetic field model (Finlay et al. 2020).

Sequence A	Prior free run	Forecast	Forecast	Earth’s core
Per cent of path	0	43	100
Time span	366 300 yr	2000–2080	2001–2081
Ensemble		A1–A42	A1–A4
τη (yr)	247 000	247 000	247 000	≈2 × 105
|$\sigma \, (10^{6} \mathrm{S\,m^{-1})}$|	1.21	1.21	1.21	≈1.1
|$\tau _\mathrm{SV}^{1}$| (yr)	550	412	393	420 ± 50
τU (yr)	169	139	131	≈130
τA (yr)	116	17.5	1.74	≈2
2πτΩ (d)	4248	134	1.34	1
p (TW)	1.27	2.19	2.28	≈3
B (mT)	4.1	4.8	4.8	≈4
|$B_\mathrm{o}^\mathrm{12}$| (mT)	0.60	0.40	0.40	0.40
Sequence B	Prior free run	Forecast	Forecast	Earth’s core
Per cent of path	0	43	100
Time span	308 000 yr	2000–2080	2001–2081
Ensemble		B1–B42	B1–B4
τη (yr)	183 852	183 852	183 852	≈2 × 105
|$\sigma \, (10^{6} \mathrm{S\,m^{-1})}$|	0.9	0.9	0.9	≈1.1
|$\tau _\mathrm{SV}^{1}$| (yr)	497	343	376	420 ± 50
τU (yr)	158	124	129	≈130
τA (yr)	89	15	1.5	≈2
2πτΩ (d)	3162	100	1	1
p (TW)	1.88	3.31	3.19	≈3
B (mT)	5.29	5.62	5.60	≈4
|$B_\mathrm{o}^\mathrm{12}$| (mT)	0.69	0.41	0.41	0.40

Our choices for τ_η in Table 2 assume typically high electrical conductivities σ in the outer core. This hypothesis has recently gained traction following progress in ab-initio computations of material properties at core conditions (see a review in Davies et al. 2015). This implies that our models operate in a magnetic Reynolds number range Rm = τ_η/τ_U = 1200–1900 higher than Rm ≈ 1000 used in Aubert et al. (2017) and following studies. The scaling theory presented in these studies predicts that Rm remains close to invariant along the path, such that τ_U should be approximately constant at constant τ_η. Values of τ_U reported in Table 2 emphasize deviations to this theory that are well below an order of magnitude, with a drift towards shorter overturn times being observed at advanced path positions. Likewise, the Alfvén timescale decreases a bit more than the theoretical prediction τ_A = ϵ^1/4τ_A(0). Similar deviations have been previously identified in free runs (Aubert 2019; Aubert & Gillet 2021), but they are made slightly stronger here because of the energy injected as the states are corrected by data assimilation. To handle these, the power of the start-of-path prior models is therefore reduced relatively to the value that would be appropriate if the scaling theory was perfectly respected. As also largely anticipated from this theory, our two sequences A and B terminate the path with nearly Earth-like values for all properties listed in Table 2. The remaining differences between models and Earth reflect slightly different design choices for sequences A and B. Sequence A does not fully target Earth’s rotation rate and expected convective power at the end of the path. It does however match the expected amplitude B of the rms magnetic field inside the core more closely than sequence B, which is constructed to exactly reach |$2\pi \tau _{\Omega }=1 \mathrm{d}$|⁠, at the expense of a stronger convective power and internal magnetic field. Both free runs of prior models PA and PB overestimate by about 50 per cent the large-scale magnetic field |$B_\mathrm{o}^\mathrm{12}$| at the core surface, which is then brought to the observed value by assimilation.

At or above the core surface, we also examine the variation (SV) |$\dot{\boldsymbol{\mathrm{B}}}$| and acceleration (SA) |$\ddot{\boldsymbol{\mathrm{B}}}$| of the magnetic field. As in Aubert (2018), we compute the spectrum of typical SV and SA timescales

$$\begin{eqnarray} \tau _{\mathrm{SV}} (\ell ) = \displaystyle \sqrt{\left\langle E_{\mathrm{B}}(\ell ,t)\right\rangle / \left\langle E_{\mathrm{\dot{B}}}(\ell ,t) \right\rangle }, \end{eqnarray}$$

(24)

$$\begin{eqnarray} \tau _{\mathrm{SA}} (\ell ) = \displaystyle \sqrt{\left\langle E_{\mathrm{\dot{B}}}(\ell ,t)\right\rangle / \left\langle E_{\mathrm{\ddot{B}}}(\ell ,t) \right\rangle }, \end{eqnarray}$$

(25)

where t is time, the angle brackets denote a time averaging and E_B(ℓ, t), |$E_{\mathrm{\dot{B}}}(\ell ,t)$| and |$E_{\mathrm{\ddot{B}}}(\ell ,t)$| are the energies respectively contained into the magnetic field |$\boldsymbol{\mathrm{B}}$|⁠, its SV and SA any height above the core surface. A least-squares regression |$\tau _{\mathrm{SV}} (\ell )=\tau _\mathrm{SV}^{1}/\ell$| is performed in the range ℓ = 2–13 on the SV timescales to obtain the master secular variation timescale |$\tau _\mathrm{SV}^{1}$|⁠. Values reported in Table 2 confirm the relationship |$\tau _\mathrm{SV}^{1}\approx 3\tau _{U}$| obtained by Christensen et al. (2012) in a systematic survey, and are Earth-like at the end of the path.

Similarly to earlier studies (Aubert & Finlay 2019; Aubert et al. 2022), at the surface of the Earth S_E we define the time-series E_SA(t) of the rms energy contained into the SA as

$$\begin{eqnarray} E_\mathrm{SA}(t)= \dfrac{1}{S} \int _{S} \ddot{\boldsymbol{\mathrm{B}}}^{2}(a,\theta ,\varphi ,t)~\mathrm{d}S, \end{eqnarray}$$

(26)

and the jerk energy E_J as a sliding energy difference of Earth surface magnetic acceleration averaged over consecutive 3-yr windows:

$$\begin{eqnarray} E_\mathrm{J} (t) = \dfrac{1}{S_\mathrm{E}} \int _{S_{\mathrm{E}}} \left(\left\langle \ddot{\boldsymbol{\mathrm{B}}}\right\rangle _{t}^{t+3 \mathrm{yr}}-\left\langle \ddot{\boldsymbol{\mathrm{B}}}\right\rangle _{t-3 \mathrm{yr}}^{t}\right)^{2}(a,\theta ,\varphi ,t)~\mathrm{d}S. \end{eqnarray}$$

(27)

Finally, we report on the level of enforcement of the Taylor constraint (Taylor 1963) in our models, by computing the Taylor integral |$\mathcal {T}(s)$| defined as in Wicht & Christensen (2010), Teed et al. (2014), Aubert et al. (2017) and Aubert & Gillet (2021):

$$\begin{eqnarray} \mathcal {T}(s)=\dfrac{ \displaystyle \int _{z_{-}}^{z_{+}} \boldsymbol{\mathrm{e}}_{\varphi }\cdot \left\langle \left(\nabla \times \boldsymbol{\mathrm{B}}\right)\times \boldsymbol{\mathrm{B}}\right\rangle _{\varphi }\, \mathrm{d}z}{\displaystyle \int _{z_{-}}^{z_{+}} \left|\, \boldsymbol{\mathrm{e}}_{\varphi }\cdot \left\langle \left(\nabla \times \boldsymbol{\mathrm{B}}\right)\times \boldsymbol{\mathrm{B}}\right\rangle _{\varphi }\, \right|\, \mathrm{d}z}. \end{eqnarray}$$

(28)

Here, the brackets |$\left<\right>_\varphi$| refer to an average in the azimuthal direction, and the integrals are evaluated between the lower and upper heights z_{−, +} at which an axial cylinder of cylindrical radius s intersects the spherical shell boundaries. We report on the average |$\mathcal {T}$| over cylindrical radii of |$|\mathcal {T}(s)|$|⁠.

3 RESULTS

3.1 Dynamically and morphologically relevant trajectories

To first assess the dynamical relevance of trajectories obtained in sequence A, we represent in Fig. 1 the root-mean-squared amplitude of each force as a function of the spherical harmonic degree ℓ (see computation details in Aubert et al. 2017). The results largely comply with the expected theoretical behaviour previously described in Aubert et al. (2017), Aubert (2019) and Aubert & Gillet (2021). Throughout the sequence, a leading-order QG balance is observed between the pressure and Coriolis force. The MAC balance follows at next order between the buoyancy, Lorentz and the part of Coriolis force not balanced by pressure. This QG-MAC balance remains remarkably stable across epochs and path positions. Likewise, the harmonic degree at which the Lorentz and buoyancy forces are in balance remains close to ℓ_⊥ ≈ 10 throughout the sequence, indicating a weakly variable dominant length scale of convection d_⊥ = πD/ℓ_⊥ in the fluid shell (Schwaiger et al. 2021). Deviations from the QG-MAC balance are balanced at the next order by inertia, with viscosity coming last in the hierarchy and being essentially irrelevant to the large-scale dynamics of the system. The relative level of inertia decreases as the path position is advanced during the progression of the sequence, accordingly with the theoretical evolution of the Rossby number Ro = U/ΩD ∼ ϵ^1/2 that compares the amplitudes of the inertial and Coriolis forces (Aubert et al. 2017). The level at which the QG-MAC balance is satisfied (i.e. the amplitude ratio of the MAC forces over inertia which forms the residue of the balance) therefore increases as we advance along the path.

Figure 1.

Root-mean-squared amplitude of forces obtained within the fluid shell for ensemble member A3, during the course of sequence A. The forces are presented as functions of the spherical harmonic degree ℓ, and normalized respectively to the maximum reached by the pressure gradient (see methodological details in Aubert et al. 2017). Panels (a)–(c) are obtained during the assimilation sequence at the end of forecasts, and panels (d) and (e) during the final forecast part of the sequence started from dynamically balanced states.

Dynamical relevance can be further checked by examining in Fig. 2 the enforcement of the Taylor constraint (eq. 28) during the assimilation. The quantity |$\mathcal {T}$| evaluates the axisymmetric azimuthal Lorentz force acting on axial cylinders in the fluid shell. Because the Coriolis, pressure and buoyancy forces vanish on these cylinders, only the next-order forces (inertia and viscosity) can balance this component, and |$\mathcal {T}$| is therefore expected to decrease accordingly on approach of Earth’s core conditions. The Taylor constraint is restored against deviations by the propagation of torsional Alfvén waves of typical timescales τ_A. In the assimilation part of the sequence (Fig. 2), |$\mathcal {T}$| is computed immediately after analyses, in order to weigh this dynamical compensation against the deviations introduced by the statistical corrections brought by these analyses. Fig. 2 shows that sequential assimilation is successful in absorbing these deviations and gradually decreasing |$\mathcal {T}$|⁠, though at levels that do not fully reach those typical of free runs (Aubert & Gillet 2021). This is expected as the time span between analyses (20 or 10 yr) tends to not leave enough time for a sufficient number of torsional wave periods to occur (see values of τ_A in Table 2). This problem is partly solved in final forecasts, by further dynamically balancing the last analysis through a 2-yr run performed at 114 per cent of the path (Section 2.5), the time average of which averages out a significant number of torsional wave overtones of fundamental timescale |$\tau _{A}\approx 1\mathrm{yr}$|⁠. As seen in Fig. 2, the time integration of a balanced state at 43, 71 and 100 per cent of the path then produces levels of |$\mathcal {T}$| in line with the values expected from free runs, except at 100 per cent of the path where |$\mathcal {T}$| remains slightly above an extrapolation of the results of Aubert & Gillet (2021). This is presumably the effect of the moderate resolution and high viscous hyperdiffusivity at small scale (Fig. 1), which have been adopted here as a trade-off for computational tractability.

Figure 2.

Average level |$\mathcal {T}$| of the absolute Taylor integral |$|\mathcal {T}(s)|$| (eq. 28) within the fluid shell, presented as a function of epoch during the evolution of sequence A. In the assimilation part of the sequence (1840–2000), |$\mathcal {T}$| is reported as an ensemble average and ±1 std. dev. of A1–A42 values computed immediately after analyses. In the forecast part of the sequence (2001–2080), |$\mathcal {T}$| is reported every 8 yr from 2009 onwards during the numerical integration of the dynamically balanced ensemble member A3 at 43, 71 and 100 per cent of the path. The horizontal coloured lines mark the expected values of |$\mathcal {T}$| interpolated and extrapolated from previous free runs at path positions up to 71 per cent (fig. 1f in Aubert & Gillet 2021).

To assess the morphological relevance of data assimilation trajectories, in Fig. 3 we next examine the temporal evolution of the classical assimilation diagnostics (eqs 20–23). The normalized prior misfit to MF data Δ_MF (Fig. 3a) remains less than unity throughout the sequence, meaning that between 1840 and 1990, the analyses adhere to COV-OBS within the supplied error bounds, and in 2000 to CHAOS-7.9 within less than |$e_{B}=1\mathrm{nT}$| at Earth’s surface. The normalized posterior increment Δ_B also remains less than unity, meaning that the corrections brought to the core surface magnetic state vector are within the variability of the prior model. The largest correction logically needs to be applied to the initial model state that is unconstrained by the data, after which subsequent corrections are more minor. The normalized prior misfit to SV data Δ_SV (Fig. 3b) decreases during the sequence and remains in the vicinity of unity, meaning that the inverse core flow problem framework is able to fit the SV within the error bounds supplied for COV-OBS for the first analyses at the start of the path (1840, 1860 and 1880), and within |$e_{\dot{B}}=20\mathrm{nT\,yr^{-1}}$| at Earth’s surface for all subsequent analyses using an underlying model that advances along the path. The normalized posterior increment to core surface flow Δ_u is kept in the vicinity of 0.25, meaning that the correction brought by each analysis to the core flow is well below the natural variability of the prior. As anticipated in Section 2.4.1, obtaining reasonable levels for the inertial response in Fig. 1 is tied to the specification of a loose margin of error for |$e_{\dot{B}}$|⁠.

Figure 3.

Assimilation diagnostics (eqs 20–23, ensemble averages over A1–A42 and ±1 std. dev.) obtained at the successive analyses of sequence A. (a) Normalized prior misfit Δ_MF to MF data, and posterior increment Δ_B brought on core surface MF. (b) Normalized prior misfit Δ_SV to SV data, and posterior increment Δ_u brought on core surface flow. The coloured lines recalls the path position of the underlying dynamical model.

Fig. 4 presents the temporal evolution of two large-scale Gauss coefficients during sequence A versus that of the input geomagnetic field models. As the SV is assimilated together with the MF through the resolution of the core flow inverse problem at each analysis step (eq. 19), it is instructive to diagnose both the evolution of MF (Figs 4a and c) and SV (Figs 4b and d) coefficients, with the SV enhancing deviations between truth and assimilation that appear otherwise minimal on MF trajectories. As previously discussed with single-epoch forecasts (Aubert 2015), the observed rates of axial dipole decay |$\dot{g}^{1}_{0}$| in the epoch range 1840–2010 are frequently larger than those typically obtained in free runs of the parent prior model PA and its child versions along the path, such that forecasts runs tend to favour a return to more moderate values of dipole decay (Fig. 4b). While still present in the axial dipole decay after 1990, this tendency is more in line with the observed geomagnetic evolution in this epoch range. Coefficients other than the axial dipole (e.g. Fig. 4d) do not appear to suffer from a similar forecast bias. Forecasting the ensemble from 2000 until 2080 at 43 per cent of the path (light blue shaded regions in Fig. 4) predicts an evolution essentially in line with that of the geomagnetic field until 2020. As in Aubert (2015) and Finlay et al. (2016), the ensemble average subsequently predicts that dipole decay appears set to continue in the next decade, after which the ensemble spread renders the prediction less certain. Ensemble averages should indeed be interpreted with caution, as individual trajectories such as that of ensemble member A3 in Figs 4(b) and (d) can show a rich, non-monotonous behaviour. Despite the different parameters used for A3 forecasts at 43, 71 and 100 per cent of the path, a significant overlap is seen in the first 50 yr of their decadal evolution, after which the trajectories diverge owing to the typical e-folding time of such systems (Hulot et al. 2010). This is expected because a similar QG-MAC balance is enforced in the three runs and governs the slow evolution of the system. The less obvious evolution of the frequency content associated to interannual dynamics will be examined in greater detail in Section 3.3.

Figure 4.

(a) and (c) Temporal evolution of Gauss coefficients |$g_{1}^{0}(t)$| and |$h_{1}^{1}(t)$| and (b) and (d) their time derivatives |$\dot{g}_{1}^{0}(t)$| and |$\dot{h}_{1}^{1}(t)$| during sequence A. The geomagnetic field models COV-OBS and CHAOS-7 serving as input to data assimilation are shown in black, together with the A1–A42 ensemble averages and ±1 std. dev. of assimilation trajectories. The sequence is terminated by 80-yr forecasts of ensemble members A1–A42 at 43 per cent of the path, and three individual forecasts at 43, 71 and 100 per cent started from the dynamically balanced final analysis of ensemble member A3.

As seen in Figs 4(b) and (d), fast inertial transients of sub-annual period commensurate with τ_Ω are of low amplitude at the beginning of forecasts at 100 per cent of the path, and vanish after about a decade. This is a benefit of the dynamical balancing procedure introduced in Section 2.5. Fig. 5(a) shows that about 15 yr after their initiation at epoch 2001, the trajectories feature a core surface MF, Earth surface SV and SA that can relevantly be compared to the observations extracted from the geomagnetic field model CHAOS-7.9 used as input. The amplitude and structure of the MF and SV forecasts are in broad agreement with the geomagnetic observations, as are also those of the SA, a quantity that is not assimilated but entirely an output of the underlying dynamical model. Likewise (Fig. 5b), the typical SV timescales τ_SV(ℓ) of forecast runs are close to those derived from the MF and SV of CHAOS-7.9, as also attested by their adherence to the power law |$\tau _\mathrm{SV}(\ell )=\tau _\mathrm{SV}^{1}/\ell$| for ℓ = 2–13, with master timescales |$\tau _\mathrm{SV}^1$| that are close to each other (⁠|$\tau _\mathrm{SV}^1=410 \mathrm{yr}$| for CHAOS-7.9 versus 379 yr for the average of forecast runs). At large scale ℓ = 1–5, the underlying dynamical model at 100 per cent of the path also produces typical SA timescales τ_SA(ℓ) in line with those derived from the SV and SA of CHAOS-7.9. At larger degrees, the divergence between the SA timescale of forecast runs and that of CHAOS-7.9 is presumably informative on the regularization and source separation hypotheses adopted when constructing the geomagnetic field model (Finlay et al. 2020). The SA timescale is, for instance, highly sensitive to the level of regularization (Lesur et al. 2022), with typical values of τ_SA(ℓ) below 10 yr being observed in weakly regularized models.

Figure 5.

(a) and (b) From top to bottom, Hammer projections of the radial magnetic field at the core surface, of the radial SV and SA at Earth’s surface, in (a) the geomagnetic field model CHAOS-7.9 in 2016.5, and (b) the forecast run of ensemble member A4 using an underlying model at 100 per cent of the path, at the same epoch. All fields are presented up to spherical harmonic degree and order 13. (c) Magnetic variation and acceleration timescales τ_SV(ℓ) and τ_SA(ℓ) (eqs 24 and 25), presented as functions of the spherical harmonic degree ℓ, for geomagnetic field model CHAOS-7.9 (black) averaged over 1997–2022, and for the average of forecast runs A1–A4, B1–B4 over 2011–2080 (with the latter curves being smoother because of averaging over a longer time and over ensemble members). Also presented are best fits of τ_SV(ℓ) to the functional form |$\tau _\mathrm{SV}^{1}/\ell$| to the data and model results, with the corresponding values of |$\tau _\mathrm{SV}^{1}$| being also reported.

To summarize, trajectories obtained at 100 per cent of the path can be considered morphologically relevant as they are anchored to the past evolution of the geomagnetic signal and forecast its present and future evolution reasonably well (Figs 3–5). They can also be considered dynamically relevant as they respect the leading-order QG-MAC force balance expected in Earth’s core, with inertial deviations to this balance essentially in line with theoretical expectations (Figs 1 and 2). That these trajectories continue the system evolution into the future in a manner that respects the characteristics of past temporal variations (Figs 4 and 5) is another manifestation of their combined morphological and dynamical relevance.

3.2 Present geodynamical state of Earth’s core

Fig. 6 documents the ensemble-averaged estimated state of Earth’s core at epoch 2000, immediately after the last analysis of sequence B. The assimilation essentially preserves the results obtained using earlier inversion, assimilation schemes and dynamo models (Aubert 2013, 2014, 2015, 2020). The core surface magnetic field at this epoch (Fig. 6a) comprises the large-scale assimilated structure, together with some small-scale information brought by the integration of the dynamical model. The structure of the azimuthal magnetic field below the core surface (Figs 6b and c) comprises low-latitude flux concentrations, as part of the general organization of magnetic field lines that are pushed towards the outer boundary by convective plumes (see e.g. fig. 8 in Aubert 2019), together with concentrations in the vicinity of the axial cylinder tangent to Earth’s inner core. The general flow circulation in the core (Figs 6d, e and g–i) confirms previous assimilation-based determinations as well as earlier studies using different methods (Pais & Jault 2008; Gillet et al. 2013, 2015; Pais et al. 2015; Barrois et al. 2017; Livermore et al. 2017; Bärenzung et al. 2018; Gillet et al. 2019). Its structure comprises a planetary-scale, roughly axially columnar westward eccentric gyre (best seen in Fig. 6h) touching the core surface in the Atlantic hemisphere, and the inner core surface in the Pacific hemisphere, associated with upwellings beneath Asia and downwelling beneath America (Fig. 6g). This gyre accounts for significant near-equatorial westward core surface flow in the Atlantic (Fig. 6d), as well as a high-latitude jet beneath the Bering Straits (Fig. 6e). Deviations from purely axially columnar flow have been previously documented in, for example, Livermore et al. (2017) and Aubert (2020), but appear somewhat stronger with the present sequential assimilation scheme, mostly as an effect of the rather loose fit to the SV prescribed to preserve dynamical consistency at advanced path positions. As also discussed in Aubert (2015, 2020), the gyre is in thermal wind balance with a buoyancy distribution (Fig. 6f) comprising excess release in the Eastern (0°E-180°E) hemisphere. Here, we recall that the prior models used in this study do not impose any longitudinal buoyancy heterogeneity, such that the localization of this structure entirely emerges from the prescription of the geomagnetic data, and most notably the hemispherical pattern of the SV.

Figure 6.

State of Earth’s core immediately after the last analysis at epoch 2000, as estimated from the B1–B42 ensemble average of sequence B at this epoch. (a) and (b) Hammer projections of the radial magnetic field at the core surface and azimuthal field below the core surface. (c) Meridional cut of the axisymmetric azimuthal magnetic field. (d) and (e) Hammer and North polar view of the core surface flow (arrows, grey shading indicates flow velocity). (f)–(h) Equatorial cuts of the density anomaly, radial and azimuthal velocities. (i) Meridional cut of the axisymmetric azimuthal (zonal) flow. All fields are presented at native model resolution.

The assimilated states match the observed geomagnetic variations during the assimilation phase (Fig. 4), but as they are no longer analysed after 2000, the forecasts subsequently evolve freely and can present different magnetic variation levels in relationship to the different levels of convective power that they have reached. Examining these is a way to check the compatibility between forced and free evolution, and also to quantitatively estimate a few quantities of crucial geodynamical interest. Fig. 7 represents the convective power p against the the master SV timescale |$\tau _\mathrm{SV}^{1}$| and the internal magnetic field amplitude B during 80-yr runs of ensemble members computed at 43 and 100 per cent of the path. By virtue of the invariance of long-term dynamics and using the dimensioning scheme outlined in Section 2.3, similar dimensional values of p, |$\tau _\mathrm{SV}^{1}$| and B are obtained for given ensemble members irrespectively of the path position of their underlying model. The largest deviations between 43 and 100 per cent are observed for |$\tau _\mathrm{SV}^{1}$| in sequence A, mostly because a suboptimal specification of gravitational coupling scaling along the path (see Section 2.2) introduces additional SV from solid-body rotations. Using the correct scaling in sequence B largely alleviates these deviations.

Figure 7.

Convective power p as a function of (a) master SV timescale |$\tau _\mathrm{SV}^{1}$| and (b) internal magnetic field amplitude B in the 84 final 80-yr forecasts of ensemble members A1–A42, B1–B42 performed at 43 per cent of the path, as well as the 8 forecasts A1–A4, B1–B4 carried out at 100 per cent of the path and started from dynamically balanced states. Solid lines connect same ensemble members at 43 and 100 per cent of the path. The presented quantities are time averages over the duration of forecasts. Dimensional results of sequence A are extrapolated slightly beyond the end of path in order to bring its rotation rate exactly to an Earth day (see the text). Shaded grey regions locate the models compatible with an estimate |$\tau _\mathrm{SV}^{1}=370-470\mathrm{yr}$| for the geodynamo at present (Lhuillier et al. 2011), and dashed lines locate the average estimates of |$\tau _\mathrm{SV}^{1}$|⁠, p and B.

The results obtained here with models situated before the end of the path are therefore relevant to Earth’s core conditions, but as can be anticipated from Section 2.3, a precise quantitative application to these conditions requires that the end-of-path planetary rotation rate Ω_EOP of the model should be exactly that of the Earth. This is the case in sequence B but with the results of sequence A need to be further extrapolated slightly beyond the end of path, because 2π/Ω_EOP = 1.34 d for this sequence (Table 2). Since the path scaling rules are embedded into the dimensioning procedure of Section 2.3, this is simply achieved in Fig. 7 by using the exact rotation rate of Earth in place of Ω_EOP in this procedure. This brings the results of sequence A in line with those of B, and in Fig. 7(a), this reveals a correlation between convective power and master SV timescale. Retaining models in the range |$\tau _\mathrm{SV}^{1}=370-470\mathrm{yr}$| inferred for the geodynamo at present (Lhuillier et al. 2011) yields a power range 2.5–3.4 TW for compliant ensemble members, among which a mean and std. dev. |$p=2.95\pm {0.2}\mathrm{TW}$| can be estimated for the current power of Earth’s core convection. Under our assumption that the core–mantle boundary heat flow is exactly adiabatic, convection is entirely driven by inner core crystallization, in which case the thermodynamic efficiency factor of the geodynamo is estimated at f = 0.2 (Lister 2003; Aubert et al. 2009). Our estimate for p therefore translates into a core–mantle boundary heat flow |$Q=p/f=14.75 \pm 1\mathrm{TW}$|⁠. Retaining models in the power range 2.5–3.4 TW, a mean and std. dev. |$B=5.58 \pm 0.15\mathrm{mT}$| is also estimated among compliant ensemble members for the internal strength of Earth’s magnetic field (Fig. 7b). This is similar to, if somewhat above the estimate |$B\approx 4\mathrm{mT}$| previously retrieved from the observation of torsional Alfvén waves by Gillet et al. (2010), and |$B\approx 2.5\mathrm{mT}$| estimated from the dissipation of nutations by Buffett (2010).

3.3 Rapid hydromagnetic wave dynamics in forecast runs

As already seen in Fig. 4, the slow, (multidecadal to secular) evolution of forecasts initialized with the same dynamically balanced state is broadly similar across path positions. This is confirmed by time-series of SV at a given Earth surface location (Fig. 8a) as well as in frequency-domain power spectra of the Earth surface SA (Fig. 8b). Fig. 8 however shows that the subdecadal frequency content of forecast runs is variable across path positions. At 43 per cent of the path, transient inertial waves of subannual periods commensurate with the rotational time τ_Ω are still visible for a few decades in Fig. 8(a), but their frequency increases, while their duration and amplitudes drastically decrease on approach of 100 per cent. Likewise, the SA energy density spectral tails of Fig. 8(b) move to higher frequencies and lower levels as the model path position is advanced. This behaviour is consistent with the increase of τ_Ω along the path (Table 2) as well as the decreasing kinetic over magnetic energy ratio of inertial waves as their frequency increases relatively to the Alfvén frequency (e.g. Gerick et al. 2021). At 100 per cent of the path, the SV and SA therefore no longer contain an appreciable contribution from inertial waves within the interannual (0.1–1 yr⁻¹) frequency range of interest.

Figure 8.

(a) East (Y) component of SV at Earth surface location 0°N, 45°E, during forecast runs of ensemble member A2 at 43, 71 and 100 per cent of the path, plotted together with the prediction of CHAOS-7.9 at this location. The three numerical models are initiated from the same dynamically balanced state at epoch 2001. (b) Frequency-domain power spectral density of Earth surface SA, averaged over forecast runs of ensemble members A1–4, B1,2,4 and obtained from their magnetic field outputs in the range 2010–2080. Run B3 is studied separately as it produces a strong magnetic jerk event during its integration (see Fig. 11). Also reported is the SA power spectral density obtained from CHAOS-7.9. Power spectral densities are obtained with a Thomson multitaper method of concentration half-bandwidth 4/Δt, with a time-series duration |$\Delta t=70\mathrm{yr}$| for forecast runs and |$\Delta t=25\mathrm{yr}$| for CHAOS-7.9 between 1997 and 2022. The grey lines mark |$f^0$|⁠, |$f^{-2}$| and |$f^{-6}$| spectral fall-off slopes, with a transition frequency between spectral indices −2 and −6 being reported here in the vicinity of the Alfvén frequency |$f_{A}=1/\tau _{A}\approx 0.6 \mathrm{yr^{-1}}$| at 100 per cent of the path.

Advancing along the path also, the |$0.1-1\mathrm{{yr}^{-1}}$| frequency range is enriched in SA power (Fig. 8b), such that the transition between a |$f^0$| plateau and a |$f^{-2}$| fall-off occurs at a gradually increasing frequency. The increase of transition frequency and of power at given frequency are broadly similar to those observed in free runs along the path by Aubert & Gillet (2021), where they have been related to the increase of the Lundquist number S = τ_η/τ_A describing dissipative mechanisms acting on hydromagnetic waves. In this study also, a transition frequency at about 0.1|$\mathrm{yr^{-1}}$| was extrapolated for Earth’s core conditions, which is confirmed by the present results obtained at 100 per cent of the path.

At frequencies in excess of |$1/\tau _{A} \approx 0.6\mathrm{yr^{-1}}$|⁠, a sharp fall-off steeper than |$f^{-6}$| is observed in the SA power spectral density at 100 per cent of the path. A similarly steep fall-off (also with similar Alfvén-time related corner frequency) was previously observed in the flow acceleration spectra of axisymmetric torsional waves (see fig. 7 in Aubert & Gillet 2021). This hints at a wave excitation mechanism that is only efficient at large spatial scales because of the axially columnar pattern of QG waves, with little energy being observed in overtones past a few times the Alfvén frequency. Fig. 8(b) finally shows that forecast runs performed at 100 per cent of the path very well account for the SA power spectral density contained in CHAOS-7.9 at interannual frequencies. Both observational and numerical spectra also feature similarly steep power decays on approach of annual frequencies. There, the power deficit seen in CHAOS respectively to the numerical prediction could again be indicative of the regularization and source separation hypotheses adopted for this model (Finlay et al. 2020). The increase of this power deficit at frequencies in excess of about |$1/3\mathrm{yr^{-1}}$| provides a posterior justification for the use of a 3-yr running average window (e.g. Aubert & Finlay 2019; Finlay et al. 2020) to compare numerical and observational signatures of rapid core dynamics (for instance through the jerk energy time-series, eq. 27).

The insets in Fig. 8(a) illustrate in the time domain the SA power increase seen at interannual timescales in the spectral space (Fig. 8b), highlighting the gradual appearance of SV oscillations with typical periods of a few years to amplitudes of about |$10\mathrm{nT\,yr^{-1}}$|⁠. Signals of similar amplitudes and periods of about 7 yr have also been identified at geomagnetic observatories (Gillet et al. 2022), and theoretically interpreted as QGMC waves propagating at the core surface (Gerick et al. 2021). Such waves have also been reproduced in a free numerical simulation run performed at 71 per cent of the path (Aubert et al. 2022), where they are constantly excited by deep core convection. In the simulation, they have also been shown to account for interannual magnetic variations that include regularly recurring magnetic jerk events.

Forecast runs carried out at 71 and 100 per cent of the path (Figs 9a and b) confirm the presence of expected characteristics of such wave flows. Particularly prominent are wedge-shaped structures of small latitudinal (or cylindrically radial), large azimuthal scale, with patterns converging towards equator in time-latitude plots (Figs 9c and d) while featuring a westward phase drift of a few degrees per year in time–longitude representations (Figs 9e and f). The highly equator-symmetric patterns of rapid wave flows (Figs 9a–d) also logically reflects an underlying axially columnar structure that is strongly imposed by the Taylor–Proudman constraint, in sharp contrast with deviations from equator symmetry seen in total flows (e.g. Fig. 6d). Indeed, because the force balance of rapid flows is a low-level deviation from the higher-level QG-MAC balance enforced for convective flows (Fig. 1, see also Aubert & Gillet 2021), the level of quasi-geostrophy of the former is therefore much higher than that of the latter. The amplitude of waves increases from 71 to 100 per cent of the path (compare the two columns of Fig. 9) but their patterns and temporal characteristics remain remarkably stable, demonstrating a convergence of wave characteristics before the end of the path is reached. According to the path scaling rules (Aubert et al. 2017), the fundamental dimensioning timescale involved in the QGMC wave dispersion relationship of Gillet et al. (2022) is |$\tau _{A}^{2}/\tau _{\Omega }\sim (\epsilon ^{0.25})^{2}/\epsilon ^{0.5}\sim \epsilon ^{0}$| and indeed remains constant between the two runs. From 71 per cent of the path onwards, the results are in agreement with typical periods of 7 yr and amplitudes of about |$7\mathrm{km\,yr^{-1}}$| retrieved from recent geomagnetic variations by Gillet et al. (2022) and Ropp & Lesur (2023).

Figure 9.

Rapid core surface azimuthal flow during forecast runs of ensemble member A4 performed at (a), (c) and (e) 71 per cent and (b), (d) and (f) 100 per cent of the path. Flows are presented up to spherical harmonic degree 30 and high-pass filtered below periods of 10 yr. (a) and (b) Hammer projection of flows at epoch 2028. Green arrows delineate the phase propagation direction of wave patterns. (c) and (d) Time–latitude diagram of flows at longitude 52.5°W marked by a black curve in (a) and (b). (e) and (f) Time–longitude diagram of flows at equator. The slanted green lines mark phase propagation at speed |$9^o \mathrm{W}.\mathrm{yr}^{-1}$|⁠. Note the residual inertial wave transients seen in panel (e) as streaks of extremely rapid eastward propagation in the beginning of forecasts at 71 per cent of the path, and at even higher frequency and lower amplitude in panel (f).

Similarly to the behaviour seen in free model runs (Aubert et al. 2022), the temporally alternating QGMC wave flows of Fig. 9 are responsible for recurring pulses in the SA energy E_SA and magnetic jerk events with corresponding peaks in the jerk energy E_J that are out of phase (Fig. 10a). Consistently with the increase of wave amplitude, the energy of these jerks increases along the path (compare Figs 10a and b), but the timing of jerks remains consistent between 71 and 100 per cent of the path, again witnessing the convergence of wave characteristics seen in Fig. 9. At 100 per cent of the path (Fig. 10a), the SA and jerk energy time-series are highly similar in amplitude and pattern to those that can be retrieved from recent geomagnetic field models (see e.g. fig. 1c in Aubert & Finlay 2019). This therefore further strengthens an interpretation of these observations in terms of QGMC waves (Fig. 9) constantly excited by deep core convection.

Figure 10.

(a) and (c) Time-series of the SA and jerk energies E_SA and E_J in forecast runs of ensemble members B2 and B3 at 100 per cent of the path. (b) and (d) Time-series of E_J for the same ensemble members, at 43 and 71 per cent of the path. (e) Hammer projection of rapid core surface azimuthal flows (same spatial and temporal truncation is Figs 9a and b) at epoch 2048 in forecast runs of ensemble member B3 at 43 and 100 per cent of the path. Green arrows delineate the phase propagation direction of wave patterns. Orange arrows locate the origin of the perturbation at the origin of the jerk event.

In addition to regularly recurring jerks of weaker amplitude, one of the forecast runs (run B3, Figs 10c and d) produced a strong event consistently and simultaneously seen at path positions 43, 71 and 100 per cent. As categorized by Aubert et al. (2022), this ‘shallow wave’ event was triggered near the core surface by the arrival of a buoyancy plume and expulsion of magnetic flux. The occurrence of such an event therefore requires a specific morphological configuration that appears infrequent within the ensemble of states produced by the assimilation. As part of the long-term convective dynamics that is invariant along the path, the trigger remains similar along the path (Fig. 10e), but the nature and amplitude of the subsequent waves causing the jerk signature changes. At 100 per cent of the path, the jerk perturbation emits non-axisymmetric QGMC waves together with axisymmetric torsional Alfvén waves that propagate inwards in the cylindrical radial direction, while none of these waves are seen at 43 per cent of the path. This difference in wave content accounts for a significant increase in jerk energy (Figs   10c and d) along the path. Compared to other forecasts carried out at 100 per cent of the path that show a more standard behaviour, the SA frequency-domain power spectrum of ensemble member B3 shows a significantly increased power in the |$0.1-1\mathrm{yr^{-1}}$| range as well as a transition from a |$f^{-2}$| to a |$f^{-6}$| fall-off also occurring at a higher frequency (Fig. 11). Both are consistent with an excitation source (Fig. 10e) of smaller spatial scale and larger amplitude than the deep convection that routinely excites the QGMC waves at the origin of standard rapid geomagnetic variations.

Figure 11.

Frequency-domain power spectral density of Earth surface SA for the 2010–2080 forecast run of ensemble member B3 at 100 per cent of the path (orange). Also reported from Fig. 8 are the average spectra of A1–4, B1,2,4 (green) and the spectrum obtained from CHAOS-7.9 (black, see the legend of Fig. 8 for details).

While Fig. 10(e) shows the appearance of strong axisymmetric torsional Alfvén waves in the vicinity of an extreme event, such waves can also routinely be observed in forecast runs (Fig. 12), although at amplitudes about an order of magnitude below those of QGMC waves. Torsional Alfvén waves propagating in and outwards along magnetically coupled axial fluid cylinders are present at and beyond 43 per cent of the path, with propagation speeds evolving accordingly with the evolution of the Alfvén timescale τ_A of these runs (see Table 2 and reported velocities in Fig. 12). As previously illustrated in Aubert & Gillet (2021), torsional waves are of weak amplitude at 43 per cent of the path (Fig. 12a) and become more prominent at more advanced path positions. At 71 and 100 per cent of the path (Figs 12b and c), the waves are strongest during the first 10 yr of the forecast, as a result of the equilibration of residual force balance deviations. The fundamental mode at period |$\sqrt{3}r_{o}\tau _{A}/D=4.7$| yr is clearly seen at 100 per cent of the path (Fig. 12c). Its typical period, flow amplitude (about |$0.6\mathrm{km\,yr^{-1}}$|⁠), and pattern are in line with those retrieved from geomagnetic variations by Gillet et al. (2010, 2015). Note that these studies however favour a slightly longer Alfvén time |$\tau _{A}\approx 2\mathrm{yr}$|⁠, fundamental torsional mode period |$\sqrt{3}r_{o}\tau _{A}/D\approx 5-6\mathrm{yr}$|⁠, associated with an internal field amplitude |$B=4\mathrm{mT}$| slightly smaller than that obtained in the forecast runs. At less advanced path positions of 71 and 43 per cent, Figs 12(a) and (b) rather show torsional wave overtones, because fundamental modes of respective periods 14.7 and 47 yr are filtered out by the high-pass filter used in Fig. 12 to reveal rapid flow components. Positive torsional wave reflections at the equatorial core surface are best seen at 71 per cent of the path (Fig. 12b). Reflections are however significantly less obvious at 100 per cent (Fig. 12c), and appear negative in a few instances. The quality factor associated to the conducting layer present above the model surface increases as we advance along the path (Aubert & Gillet 2021), which indeed causes a decrease of the reflection coefficient as well as a reversal of its sign (Schaeffer & Jault 2016).

Figure 12.

Equator-symmetric part of the core surface axisymmetric azimuthal flow retained up to spherical harmonic degree 30, high-pass filtered below periods of 10 yr, and represented as a function of cylindrical radius (as an equivalent of latitude). This quantity accurately represents the columnar flows supporting the torsional waves inside the core (Aubert & Gillet 2021). Shown are forecast runs of ensemble member A2 at (a) 43 per cent, (b) 71 per cent and (c) 100 per cent of the path. Green lines follow the expected propagation of torsional Alfvén waves, with typical 1-D propagation speeds |$D/\sqrt{3}\tau _{A}$| reported on each panel.

4 DISCUSSION

4.1 Accounting for geomagnetic variations over years to centuries in physically realistic numerical models

Here, we have used physical models with conditions reaching those of Earth’s core to assimilate historical geomagnetic data over the past 180 yr. The success of this approach resides in trajectories that adhere to the force balance expected in Earth’s core (Figs 1 and 2) while at the same time accounting for most observed features of geomagnetic variations at periods of centuries down to years (Figs 4, 5, 8 and 10). In the assimilation part of the sequence, the inverse modelling framework embedded in the assimilation scheme has found states within the natural variability of the prior model that account reasonably well for the past geomagnetic field and SV configurations (Fig. 3), and also reproduce (Fig. 6) most of the flow, density anomaly and magnetic features previously identified in Earth’s core present state (most recently, Livermore et al. 2017; Bärenzung et al. 2018; Aubert 2020). The real compatibility check between the physical model and the geomagnetic evolution however takes place once the assimilation phase is finished, during the course of forecasts that are no longer constrained by the data. There, the evolution at secular (Figs 4 and 7) down to interannual (Figs 5, 8 and 10) timescales was shown to remain in line with past and present geomagnetic variations.

Forecasts carried out over several decades feature a rich interannual magnetic variation content, comprising the signature of interannual QG magneto-Coriolis waves with amplitudes of a few kilometres per year (Fig. 9) and significantly weaker torsional Alfvén waves (Fig. 12) with amplitudes of a few tenths of kilometres per year. This gives further support to earlier observations of these waves in Earth’s core (Gillet et al. 2010, 2015, 2022; Istas et al. 2023), with amplitude and period characteristics that are very similar to those obtained in the present numerical simulations. The ubiquity of QGMC waves at core conditions also gives support to an explanation of routinely recurring recent geomagnetic jerks in terms of a sea of waves consistently excited by deep core convection and propagating to and over the core surface (Aubert et al. 2022). Despite having been observed during only one out of eight forecast runs, a strong jerk event triggered by the arrival at the core surface of a buoyancy plume and expulsion of magnetic flux (i.e. a ‘shallow wave’ event in the terminology of Aubert et al. 2022) could be produced in morphologically and physically realistic conditions (Fig. 10). This is an important result because it strengthens the recent interpretation (Blangsbøll et al. 2022) of the 1969 historical geomagnetic jerk in terms of a similar mechanism. Because such events have remained so far elusive in the satellite geomagnetic record, demonstrating the possibility of their existence in a realistic context also encourages efforts towards maintaining low Earth orbiting, dedicated geomagnetic satellites (see Hulot et al. 2018), which in the near-future could document new occurrences of these events.

The information brought by realistic models on the frequency range of rapid geomagnetic variations (Fig. 8) will also be useful to the design and preparation of future satellite missions, and to the elaboration of next-generation global geomagnetic field models such as those of the CHAOS series (Finlay et al. 2020). Of significant importance is the existence of a rapid fall-off of power beyond the Alfvén frequency |$f_{A}=1/\tau _{A}\sim 1 \mathrm{yr^{-1}}$|⁠, which presumably indicates that the excitation of wave-driven rapid variations intrinsically occurs at large scales. This property was found to be robust across levels of convective forcing (for instance between sequences A and B), provided that these lead to similar values of the Alfvén time. The agreement of power spectra from realistic forecasts with that of the geomagnetic field model CHAOS-7.9 suggests that this model presently resolves reasonably well the fastest magnetic variations carried by hydromagnetic waves, and could do even better if spectral shapes from numerical models would be used as prior constraints. This could presumably offer new, physically informed temporal regularization hypotheses and help to better separate the signal from internal origin from external field contamination in such representations.

With the availability of numerical models that sustain the observed morphology and variations of the geomagnetic field from years to centuries, while reaching for the first time the physical conditions of Earth’s core regarding the separation of key timescales (Table 2), the remaining debates regarding the relevance of numerical simulations to describe the geodynamo should be further settled. Applying the path theory (Aubert et al. 2017) in the practical context of geomagnetic data assimilation also informs on how realistic the models should be, depending on the timescale range of geomagnetic variations that one aims at reproducing. Starting forecast runs from the same initial state and integrating them with physical models at various positions along the parameter space path, we have seen that the long-term (multidecadal to secular) output of the simulation is already representative of end-of-path (or Earth’s core) conditions at 43 per cent of this path (Figs 4, 7 and 8). This is an interesting property as it means that forecasting the long-term evolution of the geodynamo does not involve the need to simulate the system at fully realistic physical conditions. The main requirement here is that the QG-MAC force balance (Fig. 1) that describes the system evolution in this timescale range (Aubert 2020; Aubert & Gillet 2021) should be enforced with decently correct morphological conditions for the internal velocity, density anomaly and magnetic fields (Fig. 6). Reproducing the interannual range of geomagnetic variations however requires to further advance along the path, but here again, the convergence of QGMC wave properties is achieved from 71 per cent onwards (Fig. 9) and also does not require to reach fully realistic conditions. Reaching the end of path is finally only necessary if one wishes to realistically reproduce the torsional Alfvén wave dynamics (Fig. 12).

4.2 Current limitations

Our current implementation of geomagnetic data assimilation at realistic core conditions is not without weaknesses. As was announced in Introduction, the main difficulty here is to obtain states that are at the same time statistically consistent with the data and dynamically equilibrated. To address this difficulty, we have adopted a strategy consisting in letting the sequential scheme gradually enforce the relevant dynamical balance as data is assimilated. This has proved rather successful (Figs 1 and 2), but this was achieved at some cost. First, we have opted for sub-correction of the system state at analysis time (Section 2.4.1), by specifying rather loose levels of adherence to the SV data. Second, we have introduced in eq. (19) an additional time averaging of the hydrodynamic state vector over the course of assimilation forecasts, in order to further mitigate the inertial transients caused by the dynamical imbalance of statistically estimated states. Finally, we have also introduced a dynamical balancing step (Section 2.5) prior to running the final forecasts at advanced path positions, at the cost of losing some more adherence to the geomagnetic data. For all these reasons, our report has elided the traditional study of rms forecast errors against the true geomagnetic evolution (e.g. Fournier et al. 2021b), to rather focus on the overall morphological and energetic agreement between simulated and observed geomagnetic variations (Figs 5 and 8), as well as on the underlying physical mechanisms at work in the simulations (Figs 9 and 12). Despite this careful handling of transients, the key timescales of models (Table 2) tend to drift a bit more than (but remain well within the order of magnitude of) scaling predictions from the path theory (Aubert et al. 2017). Further improvements in data assimilation schemes aiming at a better simultaneous exploitation of combined statistical and dynamical constraints are needed here, and could also be useful to further overcome the inherently sparse observational sampling of the geodynamo system.

Some well-known statistical peculiarities of the historical geomagnetic field, such as large, apparently periodic decadal swings in axial dipole variations, are still not correctly rendered during the assimilation phase, though subsequent multidecadal forecasts can produce these (Fig. 4). This long-standing problem of geomagnetic data assimilation frameworks (Aubert 2015, 2020) stems from typically insufficient power in temporal variations of the start-of-path prior models, which is not entirely solved by advancing the underlying physical model along the path. Further work should therefore focus on the implementation of first-order mechanisms that are possibly still missing in start-of-path parent models and their child versions along the path. Among these, it will be for instance interesting to investigate at which levels new classes of waves hosted within a convectively stable layer at the top of the core created by density stratification could act a a source for additional power in dipole variations, as was proposed by Buffett (2014).

4.3 Strengthening the geomagnetically derived geodynamical constraints

The possibility to run models at Earth’s core conditions and in a morphologically relevant state over about a century of physical time (i.e. an overturn time of the core) has enabled (Fig. 7) the formulation of constraints linking the level of geomagnetic variations to the power of Earth’s core convection. The resulting geomagnetically derived estimates |$p=2.95 \pm 0.2\mathrm{TW}$| for convective power, and |$Q=14.8 \pm 1\mathrm{TW}$| for the core–mantle boundary heat flow (assuming that this boundary is neutrally stratified with respect to convection), are within the range of earlier determinations obtained from other various disciplines (Frost et al. 2022). The suite of physical models underlying the assimilation has been formulated consistently with such rather high expected values of core–mantle boundary heat flow, by implementing values of electrical conductivity (Table 2) in the upper range of current estimates (Davies et al. 2015). However, because magnetic diffusion plays only a secondary role on the generation of magnetic variations (e.g. Christensen et al. 2012), the correlation presented in Fig. 7(a) is essentially a relationship between convective power and flow amplitude, which itself is also largely diffusivity-invariant according to theoretical and numerical scaling studies (Davidson 2013; Aubert et al. 2017). For this reason, the above power and heat flow estimates are fairly robust against the value of electrical conductivity assumed for the underlying models, provided that the magnetic Reynolds number Rm = τ_η/τ_U reaches values in excess of a thousand.

The availability of morphologically and physically realistic models of the geodynamo now appears in position to place such geomagnetically derived geodynamical constraints at similar levels of reliability as those traditionally obtained from other geophysical signals. In particular, it should be possible to use data assimilation in realistic models in order to carry further systematic investigations of physical properties of the deep Earth that are not well constrained at present. Among these, the electrical conductivity of the core could for instance be possibly investigated through the level at which it affects the damping of rapid magneto-Coriolis waves, and that of the deep mantle through the reflexion properties of torsional waves at the core surface (see also Schaeffer & Jault 2016). As stated above also, the presence of density stratification near the core surface could also create new rapid geomagnetic signatures corresponding to the propagation of MAC waves within a stable layer (e.g. Buffett 2014; Chi-Durán et al. 2021). Including these ingredients and systematically studying their effects within the present assimilation framework should provide new ways to constrain the corresponding physical properties from the standpoint of geomagnetic variations.

More generally, one can think of a suite of models computed along the parameter space path as describing different timescale ranges of the same parent physical system. An emerging challenge for numerical geodynamo modelling is therefore to design a holistic parent system, the child versions of which along the path should be able to account for the whole spectrum of geomagnetic variations, from years to billion years. The prospect of using this whole spectrum of data to constrain the physical properties of a unique system, and thereby those of Earth’s core, also certainly appears geophysically desirable. First-order stumbling blocks however still need to be removed. On millenial timescales, recent palaeomagnetic results (Nilsson et al. 2022) suggest that longitudinally hemispheric geomagnetic anomalies have been persistent but not necessarily favoured at a specific longitude, in apparent contrast with mechanisms relying on control from buoyancy heterogeneities imposed at the inner core boundary (Aubert et al. 2013; Deguen et al. 2014) or core–mantle boundary (Aubert et al. 2008; Gubbins et al. 2011; Mound & Davies 2023). As this suggests that the eccentric gyre imaged at present (Fig. 6) is an intrinsic property of core dynamics rather than a consequence of boundary control, further studies of the dynamics giving rise to this structure and its temporal evolution seem warranted. Likewise, on timescales of million years, current geodynamo models are still unable to reproduce geomagnetic reversals in asymptotic conditions of low inertia (Tassin et al. 2021). Here also, new mechanisms need to be sought, that are able to produce these key features with appropriate morphological properties as measured by relevant compliance criteria (Sprain et al. 2019; Meduri et al. 2021). By presenting here a first step focused on the reproduction of timescales from years to centuries, this study is meant to stimulate efforts towards achieving this goal.

Acknowledgement

JA thanks two anonymous referees for reviews and Christopher C. Finlay for discussions and comments on the manuscript. This project has been funded by ESA in the framework of EO Science for Society, through contract 4000127193/19/NL/IA (SWARM + 4D Deep Earth: Core). Numerical computations were performed at S-CAPAD, IPGP and using HPC resources from GENCI-TGCC and, GENCI-IDRIS and GENCI-CINES (grant numbers A0120402122 and A0140402122). The results presented in this work rely on data collected at magnetic observatories. The national institutes that support them and INTERMAGNET are thanked for promoting high standards of magnetic observatory practice (www.intermagnet.org).

DATA AVAILABILITY

The numerical code and simulation data analysed in this study are available from the corresponding author upon reasonable request. Data from the numerical model are also available online at the URLs https://4d-earth-swarm.univ-grenoble-alpes.fr/datahttps://www.ipgp.fr/~aubert/4dearth.