# Rapid geomagnetic variations and stable stratification at the top of Earth's core

Julien Aubert

► To cite this version:

Julien Aubert. Rapid geomagnetic variations and stable stratification at the top of Earth's core. Physics of the Earth and Planetary Interiors, 2025, 362, pp.107335. 10.1016/j.pepi.2025.107335. hal-04997794

HAL Id: hal-04997794

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

Submitted on 19 Mar 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Rapid geomagnetic variations and stable stratification at the top of Earth's core

Julien Aubert*

Université Paris Cité, Institut de physique du globe de Paris, CNRS, F-75005 Paris, France

## ARTICLE INFO

### Keywords:

- Earth's core

- Geomagnetism

- Geodynamo

- Magnetohydrodynamics

## ABSTRACT

Probing the possible presence and physical properties of a stably stratified layer atop Earth's core is crucial to better determine the past history and heat budget of the planet. This has previously been done by ascribing a variety of interannual to decadal geomagnetic variations to hydromagnetic waves internal to the layer. This study presents the first self-consistent simulation of the stratified layer dynamics in interplay with the underlying core convection, in physical conditions matching those of Earth's core. Magneto-Archimedes-Coriolis waves of decadal periods appear in stratified layers deeper than a few tens of kilometers and with Brunt-Väisälä frequency matching the rotation rate of the planet. However, the level at which core convection excites these waves is generally insufficient to account for observed geomagnetic variations in this period range. Strong stratification is furthermore deleterious to a number of observed features that unstratified models are successful at reproducing. Fluid flow at the core surface decouples from the interior and becomes strongly dissimilar to geomagnetic inferences. Magnetic jerks and their corresponding near-equatorial, rapidly alternating magnetic acceleration patterns also disappear, because the supporting interannual magneto-Coriolis waves are impeded by the stratified layer. This negative impact on the reproduction of the observed rapid geomagnetic variations limits the possible extent of a stable top layer to the first few tens of kilometers beneath Earth's core surface.

## 1. Introduction

One of the fundamental goals of geodynamics is to build a consistent cooling history for Earth from the observations of its past and present internal motion. In this respect, it is crucial to assess whether the top of Earth's core is presently stable against convective perturbations, because this has important implications on the rate at which the core presently cools and reconstruction scenarios of the planet thermal history (e.g. Labrosse, 2015; Nimmo, 2015). A stably stratified top layer can indeed arise because the core surface heat flow is below the value that can be conducted adiabatically. This hypothesis has gained traction following an upward re-evaluation of the electrical and thermal conductivities of Earth's core based on first-principle calculations (Pozzo et al., 2012; De Koker et al., 2012) and experiments (Gomi et al., 2013; Ohta et al., 2016; Hasegawa et al., 2024). This led to estimates of the adiabatic heat flow escaping the core as high as 15 TW, matching or exceeding the estimates for the total core heat flow (see e.g. Davies et al., 2015; Frost et al., 2022). The experimental evidence for high core conductivity is however not completely consensual, as there are also

reports leading to values about three times smaller (Konõpková et al., 2016). Stable top stratification of the core can also arise through chemical effects. Light elements coming from below (Gubbins and Davies, 2013; Bouffard et al., 2019), from chemical core-mantle exchanges (Buffett and Seagle, 2010; Nimmo, 2015; Russell et al., 2024), or from primordial impacts (Landeau et al., 2016) can indeed accumulate near the outer core surface. Density anomalies arising from compositional stratification are in principle larger than thermal anomalies and more likely to be seismically detectable. Strongly stable layers several hundred kilometers deep have been inferred from seismic anomalies in the outermost core (e.g. Hellfrich and Kaneshima, 2010). The interpretation of the seismological evidence however crucially depends on the base model chosen as reference (Irving et al., 2018), and there exists an ambiguity on the exact location (deep mantle or top core) to which the anomalies are ascribed (Zhang and Wen, 2023). As a result, there is no current consensus on the seismic detection of a stable layer at the top of the core.

The Earth's core hosts the geodynamo, through which the planetary magnetic field is regenerated, and the geomagnetic signal is also a

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: aubert@ipgp.fr.</footnote>
potentially important probe for top core stratification. The first-order effect of a stable top layer is to impede radial motion and to damp geomagnetic variations because of an electromagnetic skin effect (e.g. Christensen, 2006). The resulting smoothing of the core surface magnetic field makes it very dissimilar to observations of the past centuries, which strongly limits the depth and strength of stable layers and favours the absence of top stratification (Gastine et al., 2020). Because the lower mantle structure is presumably thermally heterogeneous (e.g. Dannberg et al., 2024), the top core could still be partly unstable, while remaining compliant to some extent with the observed geomagnetic structure (Christensen, 2018; Mound et al., 2019; Mound and Davies, 2023). Also, the presence of a restoring buoyancy force may enrich the dynamics within the layer by allowing for hydromagnetic waves not seen in unstratified environments. Buffett (2014) and Chi-Durán and Buffett (2023) ascribed decadal geomagnetic variations to axisymmetric azimuthal Magneto-Archimedes-Coriolis (zonal MAC) waves within a layer about 100 km thick. Likewise, the details of interannual to decadal geomagnetic signals observed at low and high latitudes were also accounted for with MAC waves in a similar layer (Buffett and Matsui, 2019; Chi-Durán et al., 2021; Nicolas and Buffett, 2023). It remains however difficult to explain the core-originated variations in the length of the day through motion restricted to an external layer (Buffett et al., 2016). Chulliat et al. (2015) have also proposed that interannual geomagnetic variations may be explained by equatorial magnetic Rossby waves propagating within a stratified layer similar to that of Buffett (2014).

This approach is potentially powerful in constraining the physical properties of stable layers using direct observations. However, it has so far rested on a modal analysis. A family of wave functions is first obtained by solving a linear eigenvalue problem, after which the observed signal is fitted with a number of these functions. When attempting to determine whether the proposed waves are actually excited in a realistic context, and the amplitude level at which they settle, modal analysis furthermore requires a parameterised description of the convective source (e.g. Jaupart and Buffett, 2017; Nicolas and Buffett, 2023) and an idealised description of the background magnetic field in which the waves propagate. This is legitimate, given that a self-consistent simulation of core dynamics that includes both waves and their convective excitation source requires numerically challenging conditions of rapid rotation, strong power input, low dissipation and strong stratification. Direct numerical simulations nevertheless remain needed in order to ensure a fully self-consistent description of the waves, the convection and the background dynamo, and to discriminate between the different possible types of waves. This last point is important because there exist competing interpretations to MAC or magnetic Rossby waves, that do not require stable top core stratification. Rapid magneto-Coriolis (MC) waves have for instance been shown to account well for interannual geomagnetic variations observed by satellites (Gerick et al., 2021; Gillet et al., 2022), a result confirmed in self-consistent simulations (Aubert et al., 2022; Aubert, 2023).

Within the wide spectrum of geomagnetic variations, the interannual to decadal range therefore probably offers the largest dynamical richness and sensitivity to the background convective stability at the top of the core. Also, the temporal coverage of rapid geomagnetic variations offered by satellites (e.g. Finlay et al., 2023) is more detailed than the historical coverage of slowly-evolving (secular and longer) variations. Finally, the recent development of numerical geodynamo simulations able to reach the physical conditions of Earth's core for the simulation of rapid magnetohydrodynamics (Aubert et al., 2017; Aubert, 2023) now offers the possibility to simulate the entirety of this hydromagnetic wave dynamics, realistically and self-consistently coupled with the slowly-evolving underlying core convection and dynamo action. The purpose of this study is therefore to investigate this coupled dynamics and reexamine its signatures. In order to obtain better constraints on top core stratification using the rapid geomagnetic signal, we wish to assess whether the waves that have been proposed in a stably stratified context

exist at the level at which they have been inferred, and conversely, whether the waves that do not need stratification can survive if a stable layer is present.

This manuscript is laid out as follows: Section 2 introduces the numerical geodynamo model and categorises the rapid hydromagnetic waves of the system. The results are presented in Section 3, and discussed in the light of available data in Section 4.

## 2. Models and methods

### 2.1. Model description

We consider a spherical shell of thickness $D = r_o - r_i$, aspect ratio $r_i/r_o = 0.35$, representing the Earth's outer core, filled with a fluid of density $\rho$, viscosity $\nu$, thermal diffusivity $\kappa$ and electrical conductivity $\sigma$. Within the fluid domain, we numerically solve for Boussinesq convection, density anomaly transport and magnetic induction in the magnetohydrodynamic approximation. The model equations may be found in Aubert et al. (2017). The unknown fields are the velocity field $u$, magnetic field $B$ and density anomaly field $C$. The fluid domain is electromagnetically coupled with a solid inner core of radius $r_i$ and electrical conductivity $\sigma$ equal to that of the outer core, and a solid mantle between radii $r_o$ and $a = 1.83r_o$ (Earth's surface), featuring a thin conducting layer at its base of thickness $\Delta$ and conductivity $\sigma_m$ such that $\Delta\sigma_m/D\sigma = 10^{-4}$. The inner core and mantle are linked by a gravitational remote torque $\Gamma_G = -\Gamma\tau(\Omega_{IC} - \Omega_M)$ (expressed here from the point of view of the inner core) depending on their axial differential rotation $\Omega_{IC}$, $\Omega_M$ and introducing a coupling constant $\Gamma\tau$. The moments of inertia assigned to the inner core, outer core and mantle respect Earth-like proportions, and the ensemble rotates at a constant angular velocity $\Omega = \Omega_e\epsilon$ with respect to a Galilean reference frame.

The fluid shell boundaries at $r = r_i$ and $r = r_o$ respect mechanically stress-free and electromagnetically conducting boundary conditions. A homogeneous mass anomaly flux $F$ is imposed at radius $r = r_i$ to drive convection. The reference model imposes a vanishing mass anomaly flux at $r = r_o$, and mass conservation is ensured by a uniform volumetric buoyancy sink term in the density anomaly transport equation (see Aubert et al., 2009). This describes a situation where buoyancy is released at the bottom of the shell from inner core freezing, and where the core surface heat flow is exactly adiabatic. The neutrally buoyant top core of this reference situation can be stably stratified by adding an adverse radial background density anomaly profile within a layer of thickness $H$:

$$ \frac{\partial C_0}{\partial r} (r \ge r_o - H) = - \frac{N^2 \rho (r - r_o + H)}{g_o H} \quad (1) $$

where $N$ is the maximal Brunt-Väisälä frequency reached at the surface of the core. In this new situation, the mass anomaly flux at $r = r_o$ is now non-vanishing, as for instance is the case if the core surface heat flow is subadiabatic. To ensure mass conservation, a corresponding correction $-\kappa\nabla^2 C_0$ is added to the uniform sink term of the density anomaly transport equation.

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

$$ Ra_F = \frac{g_o F}{4\pi\rho\Omega^3 D^4}, \quad E_\eta = \frac{\eta}{\Omega D^2}, \quad Pr = \frac{\nu}{\kappa}, \quad Pm = \frac{\nu}{\eta} \quad (2) $$

Here $g_o$ is gravity at the core surface, and the magnetic diffusivity is $\eta = 1/\mu\sigma$ with $\mu$ the magnetic permeability. Adding to these main parameters are the secondary parameters $H/D$, $N/\Omega$ and $\xi = \Gamma\tau/\rho D^5\Omega$. All input parameters for our models are summarised in supplementary Table 1. The choice of parameters is guided by the path theory (Aubert et al., 2017; Aubert, 2023), where these relate to a single master parameter $\epsilon$ through $Ra_F(\epsilon) = \epsilon Ra_F(0)$, $E_\eta(\epsilon) = \sqrt{\epsilon}E_\eta(0)$, $Pr(\epsilon) = 1$, $Pm(\epsilon) = \sqrt{\epsilon}Pm(0)$ and $\xi = \sqrt{\epsilon}\xi(0)$. Here we choose $Ra_F(0) = 2.61 \times$
$10^{-5}$, $E_{\eta}(0) = 6 \times 10^{-6}$, $Pm(0) = 5$ and $\xi(0) = 2.1$. Some of our model cases explore the physical conditions of Earth's core, which are reached for $\epsilon = 10^{-7}$ (100% of the path). We also report on solutions obtained for $\epsilon = 10^{-5}$ (71% of the path) that are less numerically costly and therefore enable to further explore the $(H/D, N/\Omega)$ parameter space. Solutions obtained at 71% and 100% in principle share the same dynamics over time scales longer than the convective overturn time of the system, a defining property of the parameter space path along which they are computed (Aubert et al., 2017). Moving along this path towards Earth's core condition, the solution however enriches in rapid hydro-magnetic waves, the nature of which is recalled in Section 2.4. The 71% and 100% path positions are sufficiently close to each other for a number of wave properties to be converged and representative of Earth's core (for instance the frequency and amplitude of magneto-Coriolis waves, Aubert et al., 2022; Aubert, 2023).

We choose a spherical coordinate system $(r, \theta, \varphi)$ with unit vectors $e_r$, $e_{\theta}$, $e_{\varphi}$, and $e_z$ the axial vector. Within this system, the numerical implementation proceeds with a spectral spherical harmonic decomposition of $\mathbf{u}$, $\mathbf{B}$ and $\mathbf{C}$ up to degree and order $\ell_{\max} = 133$, a discretisation in the radial direction on a second-order finite-differencing scheme with NR grid points, and a second-order, semi-implicit time stepping scheme. The numerical implementation uses the spherical harmonics transform library SHTns (Schaeffer, 2013) freely available at https://bitbucket.org/nschaeff/shtns. All the solutions reported here are computed at extreme physical conditions and therefore resort to a hyperdiffusive treatment applied to the velocity and density anomaly fields at small scales, but not to the magnetic field which remains fully resolved. For harmonic degrees $\ell$ larger than a cut-off $l_h = 30$, the molecular viscous and thermochemical diffusivities $\nu, \kappa$ are replaced by effective values $(\nu_{\text{eff}}, \kappa_{\text{eff}}) = (\nu, \kappa) q_h^{l_h-k_h}$. This scheme is physically justified and validated in Aubert et al. (2017); Aubert (2019).

To mitigate the numerical cost associated with obtaining solutions approaching or matching the conditions of Earth's core, the following procedure is adopted. To obtain well-equilibrated solutions at 71 and 100% of the path in the neutrally stratified case, we cast less advanced solutions that were obtained previously by applying the path scaling laws (similar to Aubert and Gillet, 2021) and apply a dynamical balancing step to reduce the transients, following the procedure described in Aubert (2023). A stably stratified layer is then introduced in the system, after which the numerical integration is continued for a duration $\tau_{\text{sim}}$ of a few decades to centuries of physical time. A statistically steady state is typically obtained for the wave dynamics during the first 20 years of this integration time, after which the properties of the system are diagnosed. The background magnetic field over which wave propagation is observed then remains similar to that of the neutrally stratified case, and is reasonably Earth-like from the standpoint of quantitative semblance criteria ($\chi^2 < 2$ as defined in Christensen et al., 2010). The potentially deleterious long-term effects of top stratification of the magnetic field morphology, as for instance previously documented by Gastine et al. (2020), are therefore left outside the scope of this study, and the focus is set on the rapid (interannual to decadal) dynamics only.

## 2.2. Outputs

Numerical models are first characterised by their root-mean squared velocity $\mathbf{U}$ and magnetic field amplitude $\mathbf{B}$ in the shell. From this the overturn time $\tau_U = D/U$ and Alfvén time $\tau_A = D\sqrt{\rho\mu}/B$ are introduced. We also introduce two other important time scales, the rotational time $\tau_{\Omega} = 1/\Omega$ and magnetic diffusion time $\tau_{\eta} = D^2/\eta$. In supplementary Table 1, the relative values of these four times are set by the magnetic Ekman number $E_{\eta} = \tau_{\Omega}/\tau_{\eta}$, the Lehnert number $\lambda = \tau_A/\tau_{\eta}$ and the magnetic Reynolds number $Rm = \tau_{\eta}/\tau_U$. Models computed at the end of the parameter space path respect the values of $E_{\eta}$, $\lambda$ and $Rm$ relevant to Earth's core.

The magnetic field vector $\mathbf{B}$ is monitored at $r = r_o$, the model surface, up to spherical harmonic degree 13, and can be upward continued to the Earth surface by assuming the mantle to be an electrical insulator. From this, the magnetic variation and acceleration vectors $\ddot{\mathbf{B}} = \partial\mathbf{B}/\partial t$ and $\ddot{\mathbf{B}} = \partial^2\mathbf{B}/\partial t^2$ are obtained. At an arbitrary altitude above the model surface $r = r_o$ we can compute the time-averaged spectra $\langle B^2(\ell) \rangle$ and $\langle \ddot{B}^2(\ell) \rangle$ of the magnetic field and its temporal variations as functions of the spherical harmonic degree $\ell$. The classical secular variation time scale $\tau_{\text{SV}}^{\text{1}}$ is then obtained by adjusting by least-squares fitting a functional form $\tau_{\text{SV}}^{\ell}/\ell$ to $\tau_{\text{SV}}(\ell) = \sqrt{\langle B^2(\ell) \rangle / \langle \ddot{B}^2(\ell) \rangle}$ between degrees $\ell = 2$ and $\ell = 13$.

The magnetic acceleration is studied in detail at two key locations: the equator of the model surface, where we report on the root-mean-squared value (with the square brackets standing for time averaging in the following):

$$SA_{\text{eq}} = \left( \frac{1}{2\pi} \int_{0}^{2\pi} \left[ \ddot{\mathbf{B}}(r_o, \pi/2, \varphi) \cdot \mathbf{e}_r \right]^2 d\varphi \right)^{1/2} \quad (3)$$

and a higher latitude, where we report on the quantity

$$SA_{\text{TC}} = \left( \frac{1}{2\pi} \int_{0}^{2\pi} \frac{1}{2} \left[ \left( \ddot{\mathbf{B}}(r_o, \theta_{\text{TC}}, \varphi) \cdot \mathbf{e}_r \right)^2 + \left( \ddot{\mathbf{B}}(r_o, \pi - \theta_{\text{TC}}, \varphi) \cdot \mathbf{e}_r \right)^2 \right] d\varphi \right)^{1/2} \quad (4)$$

This averages the accelerations at latitudes $\theta_{\text{TC}}$ and $\pi - \theta_{\text{TC}}$ where the axial cylinder tangent to the inner core intersects the model surface. We also wish to track abrupt changes in $\ddot{\mathbf{B}}$ known as magnetic jerks. To this end, we use the jerk energy metric introduced in a previous study (Aubert and Finlay, 2019):

$$E_j(t) = \frac{1}{S_E} \int_{S_E} \left( \left[ \ddot{\mathbf{B}} \right]_t^{t+\Delta t} - \left[ \ddot{\mathbf{B}} \right]_{t-\Delta t}^t \right)^2 dS \quad (5)$$

where $S_E$ is the surface of the Earth, and time averaging is performed with running averages over windows of $\Delta t = 3$ yr, an interval useful for studying global variations in the secular acceleration (e.g. Finlay et al., 2020; Aubert et al., 2022).

We also monitor the velocity field $\mathbf{u}$ at the model surface up to spherical harmonic degree 30. From this we extract the rapid part $\mathbf{u}^{\text{fast}}$ by highpass filtering $\mathbf{u}$ with a fourth order Butterworth filter below periods of 10 years of physical time. We also define the scalar equatorial diagnostic

$$u_{\varphi}^{\text{fast}} = \left( \frac{1}{2\pi} \int_{0}^{2\pi} \left[ \mathbf{u}^{\text{fast}}(r_o, \pi/2, \varphi) \cdot \mathbf{e}_{\varphi} \right]^2 d\varphi \right)^{1/2} \quad (6)$$

We finally report on the convective power integrated over the core volume $V$:

$$P = g_o \int_V (\mathbf{u} \cdot \mathbf{e}_r) dV \quad (7)$$

## 2.3. Dimensioning

To cast dimensionless inputs and outputs to dimensional values, we recall here the procedure previously detailed in Aubert (2023), which fundamentally relies on the invariance along the parameter space path, of a quasi-geostrophic, magneto-Archimedes-Coriolis force balance that governs the dynamics (Aubert et al., 2017). Dimensionless lengths expressed relative to the shell gap $D$ are dimensioned by multiplying with $D = 2260$ km. The magnetic diffusion time is obtained from the end-of-path magnetic Ekman number $E_{\eta}(\text{EOP}) = \sqrt{10^{-7}E_{\eta}(0)}$ through $\tau_{\eta} = (\Omega_E E_{\eta}(\text{EOP}))^{-1}$, with the Earth rotation rate $\Omega_E = 2\pi/(1 \text{ day})$.
Dimensionless times expressed relative to the magnetic diffusion time are then dimensioned by multiplying with the resulting value $\tau_{\eta} = 229800$ yr (corresponding to $\eta = 0.70\text{m}^2/\text{s}$). Dimensionless velocities expressed in units of $\eta/D$ (magnetic Reynolds number units) are dimensioned by multiplying with $\eta/D = 9.83 \times 10^{-3}\text{km}\cdot\text{yr}^{-1}$. Dimensionless magnetic field amplitudes expressed in units of $\sqrt{\rho\mu\eta\Omega}$ (Elsser units) are dimensioned by multiplying with $\sqrt{\rho\mu\eta\Omega_E} = 0.84\text{mT}$, using $\rho = 11000\text{kg/m}^3$ and $\mu = 4\pi10^{-7}\text{H}\cdot\text{m}^{-1}$. Dimensionless density anomalies expressed in units of $\rho\Omega\eta/g_0D$ are dimensioned by multiplying with $\rho\Omega_E\eta/g_0D = 2.49 \times 10^{-8}\text{kg}\cdot\text{m}^{-3}$, using $g_0 = 10\text{m}\cdot\text{s}^{-2}$.

Supplementary Table 2 reports the dimensional values of key inputs and outputs. Our choice of parameters (particulary $Ra_F$ and $Pm$) ensures that at 100 % of the path, the Alfvén, rotational, overturn, and magnetic diffusion times $\tau_A \approx 2$ yr, $2\pi\tau\Omega = 1$ day, $\tau_U \approx 120$ yr, $\tau_{\eta} \approx 2.3 \times 10^5$ yr are all consistent with their actual or estimated values in Earth's core (Finlay et al., 2023). Concerning the magnetic diffusion time in particular, we note that the corresponding core conductivity $\sigma = 1/\mu\eta = 1.13 \times 10^6 \text{S}\cdot\text{m}^{-1}$ corresponds to the high values determined using first-principle calculations (e.g Pozzo et al., 2012). We also note that the convective power $P \approx 3$ TW of our two cases with neutrally buoyant top (without stratification, cases 71pn and 100pn) corresponds to a core heat flow $Q_{CMB}$ equal to the adiabatic value $Q_{ad} = 15$ TW consistent with such high core electrical and thermal conductivities (Aubert, 2023). A value $Q_{CMB} = 15$ TW is also reasonable in the light of recent multidisciplinary determinations (Frost et al., 2022). Starting from this reference situation, cases with increasing levels of top stratification represent situations of increasingly subadiabatic core heat flow $Q_{CMB} < Q_{ad}$, with the convective power decreasing accordingly.

## 2.4. Summary of expected hydromagnetic waves

As in our previous studies (Aubert et al., 2022; Aubert, 2023), the two cases without top stratification (71pn and 100pn) produce two types of rapid hydromagnetic waves, both previously identified in flows inferred from geomagnetic variations: torsional magneto-inertial Alfvén waves (Gillet et al., 2010) and quasi-geostrophic, magneto-Coriolis (QGMC) waves (Gillet et al., 2022). Torsional Alfvén waves are axisymmetric transverse waves propagating in the cylindrical radial direction. Their interannual period is commensurate with the Alfvén time $\tau_A$. These waves are non-dispersive. Rapid QGMC waves are non-axisymmetric, and propagate in the cylindrical radial direction as well as westwards in the azimuthal direction at the surface of the core. Their period is also interannual when they are excited by Alfvén waves (Aubert et al., 2022), and these waves are dispersive. Both types of waves feature axially invariant (columnar) wavefronts, and propagate at speeds commensurate with the Alfvén velocity $D/\tau_A \approx 1000\text{km}\cdot\text{yr}^{-1}$, largely exceeding the typical speed $D/\tau_U \approx 20\text{km}\cdot\text{yr}^{-1}$ of material advection in the core. We also expect our cases with stable top stratification to feature highly anisotropic and dispersive magneto-Archimedes-Coriolis (MAC) waves. With reference to geomagnetic variations, we are more specifically interested in 'slow' MAC waves (Buffett, 2014; Jaupart and Buffett, 2017; Buffett and Matsui, 2019; Nicolas and Buffett, 2023) selecting interannual to decadal periods much slower than the planetary rotation period, and wavevectors nearly orthogonal to the rotation and gravity vectors, in order to balance Coriolis, magnetic and buoyancy disturbances.

## 2.5. Penetration depth for rapid hydromagnetic waves

Here we finally introduce a relevant penetration depth $d_p$ for motion in the stratified layer. Studying the penetration of convection near its linear onset, Takehiro and Lister (2001) theoretically predicted that $d_p = (N/\Omega)^{-1}\delta_R$, where $\delta_R$ is the length scale of the onset convective cells. Noting that the space scale $\delta_R$ and time scale $\tau_R$ of thermal Rossby waves

at the onset of rotating convection are linked through $\delta_R/D = (\Omega\tau_R)^{-1}$ (e.g. Busse, 1970), this can be written

$$d_p = (N\tau_R)^{-1}D \qquad (8)$$

This simply states the general principle that the penetration depth is governed by the time scale ratio between disturbances at the stratification front and gravitational return to rest. Gastine et al. (2020) have shown that this theory is also valid in a situation of fully developed convection and self-sustained dynamo action. In their set-up, buoyancy sources are present immediately beneath the stratification front, such that the fluid there is slightly supercritical and thermal Rossby waves are indeed the only remaining source of motion at this location. This explains the somewhat counter-intuitive nature of the result (8), which states that penetration is ultimately governed by motion at a time scale rather different from (and much shorter than) that of the fully developed bulk convection.

Our set-up differs from Gastine et al. (2020) in the sense that there are no buoyancy sources near the stratification front, such that thermal Rossby waves are absent. At the physical conditions explored by our models, hydromagnetic waves with fundamental time scale $\tau_A$ are the most rapid source of motion there. We shall therefore explore the physics of the system in relationship with the associated penetration depth

$$d_p = (N\tau_A)^{-1}D \qquad (9)$$

and the ratio

$$H/d_p = N\tau_A H/D.$$

## 3. Results

We first study the effect of top stratification on the total core flow comprising both rapid (interannual waves) and slow (secular convection) components. In the absence of stratification (Fig. 1a), the core surface flow comprises a longitudinally hemispheric signature (see black latitudinal contours) denoting the presence of an eccentric planetary-scale gyre in the bulk (Finlay et al., 2023), high-latitude westwards (blue) polar vortices, and rapidly alternating signals at low latitude denoting the presence of QGMC waves. The QGMC waves are the fastest evolving component, with westward phase velocities of about 8 degrees per year that are much faster than the convective westward drift of about 0.2 degrees per year. With top stratification being present and $H/d_p \approx 20$ (Fig. 1b), the slowly evolving features (eccentric gyre and polar vortices) of the bulk circulation are already no longer reflected in the surface flow. The rapidly alternating, westward-drifting QGMC waves persist near the equator but with significantly larger azimuthal length scale than those obtained without stratification. This increasingly axisymmetric structure results from the increasing impediment of radial flow by stratification. At mid- to high latitudes, new flow structures are observed, with westward drifting velocities of about 4 degrees per year also much faster than the convective westward drift, and a slower equatorward-drifting component. Propagation at speeds exceeding material advection, in a direction nearly orthogonal to the rotation and gravity vectors, together with a force balance analysis (see Fig. 2 below) indicate that these correspond to slow MAC waves of decadal periods, propagating in the stably stratified layer. When $H/d_p \approx 200$, the layer dynamics is almost completely decoupled from the bulk (Fig. 1c). The only waves from the bulk that can pass through are torsional waves, because they are axisymmetric and therefore do not involve any impeded radial component. The layer dynamics otherwise only comprises azimuthally-elongated MAC waves at medium to high-latitudes. The equatorial region generally remains devoid of azimuthal MAC wave signals in Figs. 1b,c and 2, suggesting that equatorial trapping (Buffett and Matsui, 2019) is not enforced in this situation.

To ascertain the nature of the observed waves, the rapid part (time
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Snapshots of azimuthal flow **u·e**<sub>φ</sub> at the core surface in Hammer projection (top) and of axisymmetric azimuthal flow in meridional cut (bottom, together with a close-up, with the base of the stratified layer marked as a dashed line), from cases 100pn (a), 4 (b) and 100phi (c). These cases feature the same stratified layer thickness *H* = 139 km and increasing levels of top stratification. Flows are represented at the native spatial resolution of the model.</figcaption></figure>

scales below 10 yr) of their underlying force balance (more precisely the axial components of curled forces) is analysed in Fig. 2 near the top of the core. In the case where this region is neutrally buoyant, a clear magneto-Coriolis equilibrium is observed near the equator, with buoyancy being insignificant as expected. This confirms the QGMC nature of rapid waves observed there. When the top of the core is stably stratified, significant rapidly evolving buoyancy disturbances are observed at midlatitudes, in equilibrium with the Coriolis and Lorentz disturbances. This demonstrates the MAC nature of the rapid waves observed there. Some near-equatorial structures still present a magneto-Coriolis equilibrium indicating that QGMC waves also persist in this case, consistent with the observations made in Fig. 1b.

The effects of top stratification on the magnetic variation time scale τ<sup>1</sup><sub>SV</sub> are well characterised by the parameter *H*/*d*<sub>p</sub> (Fig. 3). Weakly stratified models with *H*/*d*<sub>p</sub> < 10 do not depart from the value τ<sup>1</sup><sub>SV</sub> ≈ 400 yr featured by non-stratified models (Fig. 1a), and also typical of Earth's core in the past two decades (e.g. Lhuillier et al., 2011; Christensen et al., 2012). For *H*/*d*<sub>p</sub> > 10, the changes in core surface flow

observed in Figs. 1b,c are responsible for a significant departure of τ<sup>1</sup><sub>SV</sub> above the Earth value.

Analysing the rapid (time scales below 10 years) part of core surface flow in the equatorial region (Fig. 4a) together with the magnetic signals that it creates (Fig. 5a) further reveals the typical signature of QGMC waves in the absence of top core stratification. This signature consists in latitudinally localised, rapidly alternating flows and magnetic accelerations, with foci that present a specific long-range longitudinal synchronisation (Figs. 1a, 5d, see also Aubert et al., 2022). The westward drift of QGMC wave patterns at speeds of about 8 degrees per year is also illustrated here (Fig. 4a). Moderate levels of stratification cause the selection of larger azimuthal wavelengths for the QGMC waves (Fig. 4b), following the hindrance of radial flow at the nodal points of these waves. At high stratification levels, because the remaining torsional waves have a low relative amplitude respectively to QGMC waves (Aubert and Gillet, 2021; Aubert, 2023), the patterns obtained for flow (Fig. 4c) and magnetic acceleration (Fig. 5b,e) show significantly decreased amplitudes as well as a complete longitudinal synchronisation of temporal

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Hammer projections of the curls of the buoyancy, Coriolis and Lorentz forces obtained 2 km beneath the core surface. Represented are snapshots of the *e*<sub>z</sub> components expanded up to spherical harmonic degree 30, obtained after highpass filtering the time series of forces below periods of 10 yr with a fourth order Butterworth filter. The top row presents case 71pn with a non-stratified top core and the bottom row presents case 5, with a stably stratified layer of thickness *H* = 139 km and *H*/*d*<sub>p</sub> = 72 (*N*/*Ω* = 0.89). Both cases are computed at 71 % of the parameter space path. QGMC and MAC wave signatures are identified at selected locations. Amplitudes are normalised respectively to the maximum Coriolis force in the non-stratified case.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Evolution of the time scale $\tau_{\text{SV}}^1$ for magnetic variations as a function of the ratio $H/d_p = N\tau_A H/D$ between the stratified layer thickness $H$ and the penetration distance $d_p$ at the Alfvén frequency. Half-circles outside the graph denote the values reached in non-stratified cases ($N=0$). Cases are color-coded according to their position along the parameter space path and thickness of the stable layer (see legend). The horizontal dashed line marks the value obtained from geomagnetic field model CHAOS-8.1 (Kloss et al., 2024). Grey arrows locate the models presented in Fig. 1.</figcaption></figure>

alternations that is indicative of a fully axisymmetric structure. Comparing the simulated magnetic accelerations over a few decades to the observed equatorial geomagnetic acceleration signal (Fig. 5c,f), the best match in both amplitude and pattern is obtained with Fig. 5a,d i.e. at low levels of stratification ($H/d_p < 10$).

The evolution of the rapid core surface flow amplitude $u_{\varphi}^{\text{fast}}$ with the strength of stratification can again be analysed with $H/d_p$ as the main control parameter (Fig. 6a). In particular, this reveals an intermediate range near $H/d_p \approx 10$ where $u_{\varphi}^{\text{fast}}$ can actually be stronger than in the absence of top stratification (though with some scatter). This is another consequence of the selection of larger azimuthal length scales seen in Figs. 1b and 4b. At strong levels of stratification, however, we recover the typically weak amplitudes of torsional waves. In the absence of stratification, the corresponding magnetic acceleration amplitude $SA_{\text{eq}}$

is about invariant between 71 % and 100 % of the path (Fig. 6b), similar to the flow amplitude (Fig. 6a). This is expected because the properties (amplitude and period) of the QGMC waves are already converged at 71 % of the path (Aubert, 2023). At high levels of stratification, however, $u_{\varphi}^{\text{fast}}$ remains about the same between 71 % and 100 % while $SA_{\text{eq}}$ is significantly larger at 100 % of the path. This is expected because the fundamental period $\tau_A$ of torsional waves becomes about 3 times shorter at the end of the path than at 71 % (supplementary Table 2), causing larger acceleration signals despite the flow amplitude being the same. Up to the highest levels of stratification that we investigated, the presence of a 10 km-thick layer did not appear to affect $u_{\varphi}^{\text{fast}}$ and $SA_{\text{eq}}$ (Figs. 6a,b). In this case the radial length scale of QGMC waves indeed becomes comparable to the layer thickness and hence insensitive to its presence. Fig. 6b finally confirms that non- or weakly stratified situations with $H/d_p < 10$ best reproduce the amplitude of the near-equatorial magnetic acceleration signal seen in satellite geomagnetic data.

Similar results are obtained when analysing the amplitude $SA_{\text{TC}}$ of high-latitude magnetic acceleration (Fig. 7). In the absence of top stratification, $SA_{\text{TC}}$ is about similar at 71 % and 100 % of the path, and things change at high stratification levels, with the value at 100 % of the path being significantly stronger. At such levels, the high-latitude magnetic acceleration signal reveals the signature of a MAC wave regularly drifting westwards at about 5 degrees per year (Fig. 8b), not seen in the absence of stratification (Fig. 8a). The geomagnetic observations (Fig. 8c) show some indications of a similar westward drift, as previously noted by Chi-Durán et al. (2020). However, neither the amplitude (Fig. 7) nor the irregular drift of the observed pattern are in line with the strongly stratified model of Fig. 8b, while the non-stratified case in Fig. 8a again provides a better, but this time slightly underpowered match for both. This suggests that the forcing of the reference unstratified model needs modification in order to capture a slightly higher flow power near and inside the tangent cylinder. Finally, similar to the equatorial acceleration signal, the presence of a 10-km thick layer does not appear to affect the high-latitude magnetic acceleration (Fig. 7).

The frequency-domain power spectral density of magnetic acceleration (Fig. 9a) confirms that compared to the observations, the output from strongly stratified models is underpowered at frequencies lower

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Time-longitude diagrams of the rapid equatorial flow $\mathbf{u}^{\text{fast}} \cdot \mathbf{e}_{\varphi}$ obtained in cases 100pn (a), 4 (b) and 100phi (c), similar to Fig. 1 and featuring a stratified layer of thickness $H = 139$ km. The green line in (a) is indicative of a westward drift at 8 degrees per year. Flows are represented up to spherical harmonic degree and order 30. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 5. Time-longitude diagrams of the equatorial radial magnetic acceleration $\ddot{\mathbf{B}} \cdot \mathbf{e}$, at the core surface from cases 100pn (a), 100phi (b) computed at 100 % of the path, and (c) for Earth from geomagnetic field model CHAOS-8.1 (Kloss et al., 2024). The two model cases are computed without (a) and with (b) a stratified layer of thickness $H = 139$ km. The rapid flow maps corresponding to a,b are shown in Fig. 4a,c. d-f: Hammer projections of the core surface radial magnetic acceleration at times located by horizontal dashed lines in a-c. Magnetic accelerations are represented up to spherical harmonic degree and order 13.</figcaption></figure>

than the Alfvén frequency $f_A = 1/\tau_A \approx 0.5$ yr$^{-1}$ (time scales longer than the Alfvén time $\tau_A \approx 2$ yr). At shorter time scales (higher frequencies), the spectral falloff of models at 100 % of the path is less steep than that of the geomagnetic field model CHAOS-8.1 (Kloss et al., 2024), reflecting the resolution limit of the geomagnetic field model (Aubert, 2023). The corresponding jerk energy time series (Fig. 9b) shows that the non-stratified case at 100 % of the path accounts remarkably well for the amplitude of observed geomagnetic jerks in the past two decades, as well as the 3-year recurrence of jerk energy peaks. Jerks are severely damped in the two stratified cases, which are unable to match the geomagnetic observations. Because QGMC waves (Fig. 4a) are responsible for the rapidly alternating magnetic accelerations (Fig. 5a) that create magnetic jerk events in the model, it is indeed expected that these can only occur if the core surface flow can sustain those waves.

Zonal MAC waves have been proposed as the source of a multi-decadal oscillation in temporal variations of the geomagnetic dipole (Buffett, 2014; Buffett et al., 2016; Chi-Durán and Buffett, 2023). When a stratified layer is present, the axisymmetrically averaged surface flow indeed contains the signature of a zonal MAC wave (Fig. 10b) not seen in the absence of stratification (Fig. 10a). Although the axisymmetric latitudinal (poloidal) component of the wave can in principle influence the temporal variations of the magnetic dipole, its amplitude is too weak to yield an observable oscillation on this signal (Fig. 10c). Because of the surface flow alterations brought by top stratification (see Fig. 1b,c), dipole variations of advective origin are strongly damped in Fig. 10c, revealing the smooth variations of diffusive origin.

Core-originated variations in the length of day $\Delta$LOD can be recovered using the core surface flows inferred from geomagnetic variations (Jault et al., 1988). Denoting as $t'_l$ the surface flow toroidal harmonic coefficient of degree $l$ and order $m$, the following approximation holds to a very good degree of accuracy (Jault and Finlay, 2015):

$$ \Delta\text{LOD} \approx \frac{T^2}{2\pi r_o} \frac{I_c}{I_m} \left( t'_1 + \frac{12}{7} t'_3 \right) \quad (10) $$

<figure><img src="image_7.png" /><figcaption>Fig. 6. Evolution of the root-mean-squared, rapid equatorial core surface flow $u_{\varphi}^{\text{fast}}$ (a) and magnetic acceleration $SA_{\text{eq}}$ (b) with $H/d_p$. Half-circles outside the graph denote the values reached in non-stratified cases ($H = 0$). Color coding of cases is the same as in Fig. 3. The horizontal dashed line in (b) denotes the value obtained in the geomagnetic field model CHAOS-8.1 (Kloss et al., 2024) between epochs 1997 and 2024.5. Grey arrows locate the models presented in Figs. 4 and 5.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Evolution of the root-mean-squared, high-latitude radial magnetic acceleration $SA_{TC}$ at the core surface with $H/d_p$. Half-circles outside the graph denote the values reached in non-stratified cases ($H = 0$). Color coding of cases is the same as in Fig. 3. The horizontal dashed line denotes the value obtained in the geomagnetic field model CHAOS-8.1 (Kloss et al., 2024) between epochs 1997 and 2024.5. Grey arrows locate the models presented in Fig. 8.</figcaption></figure>

where $T = 86400$ s is the reference length of the day, and $I_c, I_m$ are the respective moments of inertia of the core and mantle. In Fig. 11 we test the recovery from core surface flows against the true evolution of the length-of-day, directly computed in the simulation from the variations in the angular momentum of the core. Regardless of the strength of stratification, the agreement is found to be excellent at rapid (interannual) time scales (Fig. 11a-c), as previously reported in inferences of core surface flows from high-quality satellite geomagnetic data (Istas et al., 2023). At these time scales indeed, the magnetic coupling between the outer core and the stratified layer is efficient and torsional waves can always carry their angular momentum across the layer (Figs. 1c, 4c). At decadal and longer time scales, however, magnetic forces tend to get minimised and the magnetic coupling cannot be as efficient. The excellent recovery obtained in simulations without a stratified top layer (Fig. 11a and Schwaiger et al., 2024) is degraded with a moderate level of stratification (Fig. 11b), and lost with strong stratification levels (Fig. 11c).

## 4. Discussion

In general, the self-consistent simulation results confirm the existence of a number of waves previously proposed by modal analyses, but place tight constraints on their amplitudes. When excited by the underlying core convection, a buoyantly stable top layer hosts MAC waves at medium to high latitudes (Figs. 1,2,8,10), with interannual to multidecadal time scales. This matches earlier conclusions (Jaupart and Buffett, 2017; Nicolas and Buffett, 2023) obtained by coupling the analysis of non-stratified convective dynamos and a linear model for wave generation in the stratified layer. The dynamo simulations reported here offer the advantage of self-consistently coupling the layer with the core convection, allowing for saturation of the wave amplitude, and producing a prediction for the inductive magnetic signature associated with the waves. Comparing this latter signature with geomagnetic observations, we are however unable to confirm a number of points that were previously made in support for a geomagnetic evidence of core stratification with $N \approx \Omega_E$ and $H > 100$ km (e.g. Buffett, 2014; Chi-Durán et al., 2021). In such conditions, the amplitude of the zonal MAC wave is indeed too low in models to account for the decadal oscillations of the dipole observed in the geomagnetic record, while models without top core stratification are able to produce such oscillations at the correct level (Fig. 10). The same holds when examining the high-latitude magnetic acceleration signal (Fig. 8). High-latitude MAC waves in a stratified environment cannot match the observed amplitudes (Fig. 7), and the conservative hypothesis of an absence of stratification is here again sufficient to account for most of the signal. The models are, as always, admittedly subject to some a-priori choice such as the particular configuration of the imposed buoyancy profile (Eq. (1)) and a formulation merging thermal and chemical density anomalies. In their defense, however, they are not subject to extrapolation uncertainty as the relevant physical parameters for the simulation of rapid core dynamics are reached in some of our cases computed at 100 % of the parameter space path.

Though surface flows are more indirect observations of core dynamics than the geomagnetic variations they derive from, it is instructive to examine how much of the current understanding of the general circulation of the core is compatible with the presence of a sizeable top stratification (again, $N \approx \Omega_E$ and $H > 100$ km). At decadal time scales and longer, stable top stratification obscures the signature of the eccentric gyre and polar vortices (Fig. 1), two robustly observed features

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Time-longitude diagrams of the high-latitude radial magnetic acceleration $\ddot{B}_e$ at the northern intersection of the core surface with the inner core tangent cylinder, from cases 100pn (a), 100phi (b) computed at 100 % of the path (see Figs. 5d,e for corresponding Hammer projections) and (c) for Earth from the geomagnetic field model CHAOS-8.1 (Kloss et al., 2024). The two model cases are computed without (a) and with (b) a stratified layer of thickness $H = 139$ km. The green line in (a) is indicative of a westward drift at 5 degrees per year. Magnetic accelerations are expanded up to spherical harmonic degree and order 13. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** a: Spectral density in the frequency domain of the squared magnetic acceleration $\ddot{B}^2$ averaged over Earth's surface, represented for the geomagnetic field model CHAOS-8.1 (Kloss et al., 2024) and model cases 100pn, 100pmed and 100phi computed at 100 % of the path, with increasing levels of top stratification within a layer of thickness $H = 139$ km. These cases have a common value for the Alfvén frequency $f_A = 1/\tau_A$, which is reported as a vertical dashed line. Power spectral densities are obtained with a Thomson multitaper method of concentration half-bandwidth $4/\tau_{sim}$. The grey line marks a $f^{-6}$ spectral falloff slope. b: time series of the jerk energy $E_j$ for CHAOS-8.1 between 2000 and 2022, and the same model cases as in (a).</figcaption></figure>

of the core surface flow (see Finlay et al., 2023), and degrades the qualitative agreement between recovered and true variations in the length of the day (Fig. 11). At interannual time scales, the loss of quasi-geostrophic magneto-Coriolis waves in the surface flow (Fig. 4) is also responsible for the disappearance of rapidly alternating equatorial magnetic accelerations (Fig. 5). Weakening the equatorial rapid flow and associated magnetic acceleration (Fig. 6) causes the magnetic jerks to also disappear (Fig. 9). Equatorial magnetic Rossby waves (Chulliat et al., 2015) and equatorially trapped MAC waves (Buffett and Matsui, 2019), two alternatives that have been proposed to explain the same signal within a stably stratified environment, were not observed here. Again, a negative result is not the demonstration of the impossibility of a positive result, and geomagnetic jerks could for instance arise from another, yet uncovered physical mechanism that would be immune to top core stratification. Still, a principle of parsimony favours the absence

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** a,b: Time-latitude maps of the axisymmetrically-averaged value of the latitudinal flow $\mathbf{u} \cdot \mathbf{e}_\phi$, in cases 71pn (a) and 6 (b), computed at 71 % of the parameter space path, without (a) and with (b) a stable top layer of thickness $H = 139$ km. c: Time series of the first time derivative of the Gauss coefficient for the axial dipole (as defined e.g. in Alken et al., 2021) for the same cases.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** True (black) and recovered (green) length-of-day variations in models 100pn (a), 100pmed (b) and 100phi (c), computed at 100 % of the parameter space path, without (a) and with (b,c) a stable top layer of thickness $H = 139$ km. The recovered time series is obtained from the surface flow and Eq. (10). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

of strong stratification in any case, since the gyre, the polar vortices, the QGMC waves, their magnetic acceleration patterns and the jerks are naturally produced by the non-stratified models, at amplitude levels that match the geomagnetic inferences.

In earlier studies (e.g Christensen, 2018; Gastine et al., 2020), the role of the electromagnetic skin effect was put forward when describing the effect of the stratified layer on the core surface magnetic field structure. This effect acts on the magnetic field, and does therefore not
explain the disappearance of flow structures such as the eccentric gyre and polar vortices (Fig. 1), which are primarily sustained by a non-magnetic thermal wind balance (e.g Aubert, 2020). Our results show that it is rather the dynamical decoupling between the outer core and the top stratified layer that is responsible here, with the layer for instance sustaining lateral density gradients different from the rest of the outer core. Likewise, we emphasize that the disappearance of magnetic jerk events (Fig. 9) does not result from a skin effect either, but from the fact that the stratified layer can no longer support the associated non-axisymmetric QGMC wave flow (Fig. 4), with the layer instead supporting waves in MAC balance (Fig. 2).

We next consider how much top stratification can be allowed without degrading the agreement between the model output and geomagnetic observations. Layers with $H = 60$ km and thicker typically do not strongly alter the magnetic signal coming from below until $H/d_p \approx 10$ (Figs. 3–8) i.e. $d_p = D/N\tau_A \ge 6$ km. With $\tau_A \approx 2$ yr in the models (supplementary Table 1) as well as in Earth’s core (Gillet et al., 2010), we obtain $N \lesssim 6 \times 10^{-6}$ s$^{-1}$, or $N/\Omega_E \lesssim 0.1$. This is an order of magnitude weaker than the typical earlier geomagnetic inferences (e.g. Buffett, 2014; Chulliat et al., 2015, and subsequent references), and more than two orders of magnitude weaker than the seismic inferences (e.g. Hellfrich and Kaneshima, 2010). We therefore conclude that the evidence of rapid geomagnetic variations favours the hypothesis of vanishing stratification at the top of Earth’s core. Considering how the surface magnetic field morphology is affected on longer time scales (secular overturn and longer), Gastine et al. (2020) previously reached the same conclusion. Our study is possibly somewhat less restrictive, though, as thin layers may offer an alternative here. With $H = 10$ km, significantly stronger levels of stratification ($N \approx 3\Omega_E$) can be achieved without degrading the agreement (Figs. 6,7). In this case the radial wavelength of QGMC waves becomes larger than the layer thickness, and the waves are comparatively much less affected by stratification.

This geomagnetic constraint has direct implications on Earth’s heat budget. In the case of a thermally stable layer, we illustrate this by relating the physical properties of the stable layer to the deviation of the core surface heat flow $Q_{CMB}$ below the adiabatic value $Q_{ad}$. At the core surface, the density gradient $\partial C/\partial r$ in Eq. (1) then corresponds to a thermal anomaly gradient $\partial T/\partial r(r_o) = -N^2/\alpha g_o$, leading to

$$N^2 = \frac{\alpha g_o}{4\pi r_o^2 k} (Q_{ad} - Q_{CMB}) \quad (11)$$

where $k$ is thermal conductivity at the core surface. When $Q_{CMB}$ is near $Q_{ad}$, the depth of the stable layer can be approximated through (e.g. Lister and Buffett, 1998):

$$H = \beta r_o (1 - Q_{CMB}/Q_{ad}) \quad (12)$$

where $\beta \approx 0.3$ can be estimated from the turning point in the background density profile of the core. For illustrative purposes, here we use $k = 100 \text{ W.m}^{-1}.\text{K}^{-1}$ (Davies et al., 2015), $\alpha = 10^{-5} \text{ 1/K}$, $g_o = 10 \text{ m.s}^{-2}$, $r_o = 3485 \text{ km}$, $Q_{ad} = 15 \text{ TW}$ and a subadiabatic deviation $Q_{ad} - Q_{CMB} = 0.5 \text{ TW}$. This leads to $H = 35 \text{ km}$ and $N = 6 \times 10^{-5} \text{ s}^{-1}$, near or at the limits mentioned above. This means that only minute subadiabatic departures of the core heat flow are allowed. While this is a strong constraint, it is quite plausible that $Q_{CMB}$ could presently lie near or above 15 TW, as estimated from the mantle side (Frost et al., 2022) and from the power needed in dynamo models to reproduce the secular geomagnetic variations (Aubert, 2023). The main problem is in fact the early thermal history that this implies, and particularly the hot initial core temperature (Labrosse, 2015). We note that stronger departures of $Q_{CMB}$ below the adiabat can be allowed if the stable thermal layer is weakened by a destabilising flux of chemical species across the core-mantle boundary (e.g. an exsolution flux, Badro et al., 2018). Large lateral variations of the core heat flow relative to its mean surface average (e.g. Olson et al., 2015; Mound et al., 2019; Dannberg et al., 2024) are however unlikely to affect this constraint. When a stratified

layer of sufficient strength is present, these lateral variations indeed tend to be buffered and confined to an upper layer a few hundred meters thick (Christensen, 2018).

The geomagnetic constraint finally also severely restricts the possibility of a chemically stable top layer. This should also be limited to the first tens of kilometers beneath the core surface, because $N/\Omega_E$ values in excess of 10 are likely as soon as light components can accumulate at the top of the core (Buffett and Seagle, 2010; Gubbins and Davies, 2013). As previously noted by Gastine et al. (2020), the geomagnetic evidence is here again incompatible with the seismic evidence of a thick ($H = 300$ km) and strong ($N > 7\Omega_E$) stable layer at the top of the core (Hellfrich and Kaneshima, 2010), and favours an explanation of the corresponding anomalies by a refined equation of state (e.g. Irving et al., 2018).

## Code and data availability

The numerical code is available from the corresponding author upon reasonable request. Data supporting the findings of this study are available together with analysis codes within the supplementary Tables, at doi: 10.18715/IPGP.2024.m4gvpp59, and from the corresponding author upon request.

## CRediT authorship contribution statement

Julien Aubert: Writing – original draft, Investigation, Conceptualization.

## Author statement

Julien Aubert: conception; formal analysis; investigation; writing – original draft, writing – review and editing.

## Declaration of competing interest

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

## Acknowledgements

JA thanks two anonymous referees for useful review comments. 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 A0140402122 and A0160402122). The results presented in this work rely on data collected at magnetic observatories. The national institutes that support them and INTER-MAGNET are thanked for promoting high standards of magnetic observatory practice (www.intermagnet.org).

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pepi.2025.107335.

## Data availability

Data will be made available on request.

## References

Alken, P., Thébault, E., Beggan, C.D., Amit, H., Aubert, J., Baerenzung, J., Bondar, T.N., Brown, W.J., Calif, S., Chambodut, A., Chulliat, A., Cox, G.A., Finlay, C.C., Fournier, A., Gillet, N., Grayver, A., Hammer, M.D., Holschneider, M., Huder, L., Hulot, G., Jager, T., Kloss, C., Korte, M., Kuang, W., Kuvshinov, A., Langlais, B.,
Léger, J.M., Lesur, V., Livermore, P.W., Lowes, F.J., Macmillan, S., Magnes, W., Mandea, M., Marsal, S., Matzka, J., Metman, M.C., Minami, A., Morschhauser, A., Mound, J.E., Nair, M., Nakano, S., Olsen, N., Pavón-Carrasco, F.J., Petrov, V.G., Ropp, G., Rother, M., Sabaka, T.J., Sanchez, S., Saturnino, D., Schnepf, N.R., Shen, X., Stolle, C., Tangborn, A., Toffner-Clausen, L., Toh, H., Torta, J.M., Varner, J., Vervelidou, F., Vigneron, P., Wardinski, I., Wicht, J., Woods, A., Yang, Y., Zeren, Z., Zhou, B., 2021. International geomagnetic reference field: the thirteenth generation. Earth Planets Space 73, 49. https://doi.org/10.1186/s40623-020-01288-x.

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

Aubert, J., 2020. Recent geomagnetic variations and the force balance in Earth's core. Geophys. J. Int. 221, 378-393. https://doi.org/10.1093/gji/ggaa007.

Aubert, J., 2023. State and evolution of the geodynamo from numerical models reaching the physical conditions of Earth's core. Geophys. J. Int. 235, 468-487. https://doi.org/10.1093/gji/ggzad229.

Aubert, J., Finlay, C.C., 2019. Geomagnetic jerks and rapid hydromagnetic waves focusing at Earth's core surface. Nat. Geosci. 12, 393-398. https://doi.org/10.1038/s41561-019-0355-1.

Aubert, J., Gillet, N., 2021. The interplay of fast waves and slow convection in geodynamo simulations nearing Earth's core conditions. Geophys. J. Int. 225, 1854-1873. https://doi.org/10.1093/gji/ggab054.

Aubert, J., Labrosse, S., Poitou, C., 2009. Modelling the palaeo-evolution of the geodynamo. Geophys. J. Int. 179, 1414-1428. https://doi.org/10.1111/j.1365-246X.2009.04361.x.

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

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

Badro, J., Aubert, J., Hirose, K., Nomura, R., Blanchard, I., Borensztajn, S., Siebert, J., 2018. Magnesium partitioning between earth's mantle and core and its potential to drive an early exsulution geodynamo. Geophys. Res. Lett. 45, 13,240-13,248. doi: https://doi.org/10.1029/2018GL080405.

Bouffard, M., Choblet, G., Labrosse, S., Wicht, J., 2019. Chemical convection and stratification in the earth's outer core. Front. Earth Sci. 7. https://doi.org/10.3389/feart.2019.00099.

Buffett, B.A., 2014. Geomagnetic fluctuations reveal stable stratification at the top of the Earth's core. Nature 507, 484-487. https://doi.org/10.1038/nature13122.

Buffett, B., Matsui, H., 2019. Equatorially trapped waves in Earth's core. Geophys. J. Int. 218, 1210-1225. https://doi.org/10.1093/gji/ggz233.

Buffett, B.A., Seagle, C.T., 2010. Stratification of the top of the core due to chemical interactions with the mantle. J. Geophys. Res. 115. https://doi.org/10.1029/2009JB006751.

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

Busse, F.H., 1970. Thermal instabilities in rapidly rotating systems. J. Fluid Mech. 44, 441-460.

Chi-Durán, R., Buffett, B.A., 2023. Extracting spatial-temporal coherent patterns in geomagnetic secular variation using dynamic mode decomposition. Geophys. Res. Lett. 50, e2022GL101288. https://doi.org/10.1029/2022GL101288.

Chi-Durán, R., Avery, M.S., Knezek, N., Buffett, B.A., 2020. Decomposition of geomagnetic secular acceleration into traveling waves using complex empirical orthogonal functions. Geophys. Res. Lett. 47, e2020GL087940. https://doi.org/10.1029/2020GL087940.

Chi-Durán, R., Avery, M.S., Buffett, B.A., 2021. Signatures of high-latitude waves in observations of geomagnetic acceleration. Geophys. Res. Lett. 48, e2021GL094692. https://doi.org/10.1029/2021GL094692.

Christensen, U.R., 2006. A deep dynamo generating Mercury's magnetic field. Nature 444, 1056-1058.

Christensen, U., 2018. Geodynamo models with a stable layer and heterogeneous heat flow at the top of the core. Geophys. J. Int. 215, 1338-1351. https://doi.org/10.1093/gji/ggy352.

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

Christensen, U.R., Wardinski, I., Lesur, V., 2012. Timescales of geomagnetic secular acceleration in satellite field models and geodynamo models. Geophys. J. Int. 190, 243-254. https://doi.org/10.1111/j.1365-246X.2012.05508.x.

Chulliat, A., Alken, P., Maus, S., 2015. Fast equatorial waves propagating at the top of the Earth's core. Geophys. Res. Lett. 42, 3321-3329. https://doi.org/10.1002/2015GL064067.

Dannberg, J., Gassmoeller, R., Thallner, D., LaCombe, F., Sprain, C., 2024. Changes in core-mantle boundary heat flux patterns throughout the supercontinent cycle. Geophys. J. Int. https://doi.org/10.1093/gji/ggae075.

Davies, C., Pozzo, M., Gubbins, D., Alfe, D., 2015. Constraints from material properties on the dynamics and evolution of Earth's core. Nat. Geosci. 8, 678-685. https://doi.org/10.1038/nge02492.

De Koker, N., Steinle-Neumann, G., Vleck, V., 2012. Electrical resistivity and thermal conductivity of liquid Fe alloys at high P and T, and heat flux in Earth's core. PNAS 109, 4070-4073. https://doi.org/10.1073/pnas.1111841109.

Finlay, C.C., Kloss, C., Olsen, N., Hammer, M.D., Toffner-Clausen, L., Grayver, A., Kuvshinov, A., 2020. The CHAOS-7 geomagnetic field model and observed changes

in the South Atlantic anomaly. Earth Planets Space 72, 156. https://doi.org/10.1186/s40623-020-01252-9.

Finlay, C.C., Gillet, N., Aubert, J., Livermore, P.W., Jault, D., 2023. Gyres, jets and waves in the Earth's core. Nat. Rev. Earth Environ. 4, 377-392. https://doi.org/10.1038/s43017-023-00425-w.

Frost, D.A., Avery, M.S., Buffett, B.A., Chidester, B.A., Deng, J., Dorfman, S.M., Li, Z., Liu, L., Lv, M., Martin, J.F., 2022. Multidisciplinary constraints on the thermal-chemical boundary between Earth's core and mantle. Geochem. Geophys. Geosyst. 23, e2021GC009764. https://doi.org/10.1029/2021GC009764.

Gastine, T., Aubert, J., Fournier, A., 2020. Dynamo-based limit to the extent of a stable layer atop Earth's core. Geophys. J. Int. 222, 1433-1448. https://doi.org/10.1093/gji/ggaa250.

Gerick, F., Jault, D., Noir, J., 2021. Fast quasi-geostrophic magneto-coriolis modes in the Earth's Core. Geophys. Res. Lett. 48, e2020GL090803. https://doi.org/10.1029/2020GL090803.

Gillet, N., Jault, D., Canet, E., Fournier, A., 2010. Fast torsional waves and strong magnetic field within the Earth's core. Nature 465, 74-77. https://doi.org/10.1038/nature09010.

Gillet, N., Gerick, F., Jault, D., Schwaiger, T., Aubert, J., Istas, M., 2022. Satellite magnetic wave reveal interannual waves in Earth's core. PNAS 119, e2115258119. https://doi.org/10.1073/pnas.2115258119.

Gomi, H., Ohta, K., Hirose, K., Labrosse, S., Caracas, R., Verstraete, M.J., Hernlund, J.W., 2013. The high conductivity of iron and thermal evolution of the Earth's core. Phys. Earth Planet. Inter. 224, 88-103. https://doi.org/10.1016/j.pepi.2013.07.010.

Gubbins, D., Davies, C., 2013. The stratified layer at the core-mantle boundary caused by barodiffusion of oxygen, sulphur and silicon. Phys. Earth Planet. Inter. 215, 21-28. https://doi.org/10.1016/j.pepi.2012.11.001.

Hasegawa, A., Ohta, K., Yagi, T., Hirose, K., Yamashita, Y., 2024. Inversion of the temperature dependence of thermal conductivity of hcp iron under high pressure. Sci. Rep. 14, 23582. https://doi.org/10.1038/s41598-024-74110-3.

Helffrich, G., Kaneshima, S., 2010. Outer-core compositional stratification from observed core wave speed profiles. Nature 468, 807-810. https://doi.org/10.1038/nature09636.

Irving, J.C.E., Cottaar, S., Lekić, V., 2018. Seismically determined elastic parameters for Earth's outer core. Sci. Adv. 4, eaar2538. https://doi.org/10.1126/sciadv.aar2538.

Istas, M., Gillet, N., Finlay, C.C., Hammer, M.D., Huder, L., 2023. Transient core surface dynamics from ground and satellite geomagnetic data. Geophys. J. Int. 233, 1890-1915. https://doi.org/10.1093/gji/ggzad039.

Jault, D., Finlay, C., 2015. 8.09 - waves in the core and mechanical core-mantle interactions. In: Schubert, G. (Ed.), Treatise on Geophysics, (Second Edition). Elsevier, Oxford, pp. 225-244. https://doi.org/10.1016/B978-0-444-53802-4.00150-0.

Jault, D., Gire, C., Le Mouël, J.L., 1988. Westward drift, core motions and exchanges of angular momentum between core and mantle. Nature 333, 353-356.

Jaupart, E., Buffett, B., 2017. Generation of mac waves by convection in earth's core. Geophys. J. Int. 209, 1326-1336. https://doi.org/10.1093/gji/ggx088.

Kloss, C., Finlay, C.C., Olsen, N., Toffner-Clausen, L., Gillet, N., Grayver, A., 2024. CHAOS-8 Geomagnetic Field Model. https://doi.org/10.5281/zenodo.13950012.

Konopkova, Z., McWilliams, R.S., Gómez-Pérez, N., Goncharov, A.F., 2016. Direct measurement of thermal conductivity in solid iron at planetary core conditions. Nature 534, 99-101.

Labrosse, S., 2015. Thermal evolution of the core with a high thermal conductivity. Phys. Earth Planet. Inter. 247, 36-55. https://doi.org/10.1016/j.pepi.2015.02.002.

Landeau, M., Olson, P., Deguen, R., Hirsh, B.H., 2016. Core merging and stratification following giant impact. Nat. Geosci. 9, 786-789. https://doi.org/10.1038/ngeo2808.

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

Lister, J.R., Buffett, B.A., 1998. Stratification of the outer core at the core-mantle boundary. Phys. Earth Planet. Inter. 105, 5-19. https://doi.org/10.1016/S0031-9201(97)00082-4.

Mound, J.E., Davies, C.J., 2023. Longitudinal structure of Earth's magnetic field controlled by lower mantle heat flow. Nat. Geosci. 16, 380-385. https://doi.org/10.1038/s41561-023-01148-9.

Mound, J., Davies, C., Rost, S., Burnou, J., 2019. Regional stratification at the top of Earth's core due to core-mantle boundary heat flux variations. Nat. Geosci. 12, 575-580. https://doi.org/10.1038/s41561-019-0381-z.

Nicolas, Q., Buffett, B., 2023. Excitation of high-latitude mac waves in earth's core. Geophys. J. Int. 233, 1961-1973. https://doi.org/10.1093/gji/ggzad047.

Nimmo, F., 2015. 9.08 - thermal and compositional evolution of the core. In: Schubert, G. (Ed.), Treatise on Geophysics, (Second Edition). Elsevier, Oxford, pp. 201-219. https://doi.org/10.1016/B978-0-444-53802-4.00160-3.

Ohta, K., Kuwayama, Y., Hirose, K., Shimizu, K., Ohishi, Y., 2016. Experimental determination of the electrical resistivity of iron at Earth's core conditions. Nature 534, 95-98.

Olson, P., Deguen, R., Rudolph, M.L., Zhong, S., 2015. Core evolution driven by mantle global circulation. Phys. Earth Planet. Inter. 243, 44-55. https://doi.org/10.1016/j.pepi.2015.03.002.

Pozzo, M., Davies, C.J., Gubbins, D., Alfe, D., 2012. Thermal and electrical conductivity of iron at Earth's core conditions. Nature 485, 355-358. https://doi.org/10.1038/nature11031.

Russell, S., Irving, J.C.E., Myhill, R., Cottaar, S., 2024. The emerging picture of a complex core-mantle boundary. Nat. Commun. 15, 4569. https://doi.org/10.1038/s41467-024-48939-1.
Schaeffer, N., 2013. Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations. Geophys. Geochem. Geosystems. 14, 751–758. https://doi.org/10.1002/ggge.20071.

Schwaiger, T., Gillet, N., Jault, D., Istas, M., Mandea, M., 2024. Wave-like motions and torques in Earth’s core as inferred from geomagnetic data: a synthetic study. Phys. Earth Planet. Inter. 346, 107104. https://doi.org/10.1016/j.pepi.2023.107104.

Takehiro, S., Lister, J., 2001. Penetration of columnar convection into an outer stably stratified layer in rapidly rotating spherical fluid shells. Earth. Plant. Sci. Lett. 187, 357–366. https://doi.org/10.1016/S0012-821X(01)00283-7.

Zhang, X., Wen, L., 2023. Prem-like velocity structure in the outermost core from global sks and scs waveform modeling. Phys. Earth Planet. Inter. 344, 107091. https://doi.org/10.1016/j.pepi.2023.107091.