# Magnetic reversals in a geodynamo model with a stably-stratified layer

Nicolás P Müller, Christophe Gissinger, François Pétrélis

► To cite this version:

Nicolás P Müller, Christophe Gissinger, François Pétrélis. Magnetic reversals in a geodynamo model with a stably-stratified layer. Physics of the Earth and Planetary Interiors, 2026, 371, pp.107502. 10.1016/j.pepi.2026.107502. hal-05490592v2

HAL Id: hal-05490592

https://hal.science/hal-05490592v2

Submitted on 11 Feb 2026

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 CC BY 4.0 - Attribution - International License
Magnetic reversals in a geodynamo model with a stably-stratified layer

Nicolás P. Müller<sup>a,*</sup>, Christophe Gissinger<sup>a,b</sup> and François Pétrélis<sup>a</sup>

<sup>a</sup>Laboratoire de Physique de l'École normale supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, F-75005, Paris, France

<sup>b</sup>Institute Universitaire de France (IUF), Paris, France

ARTICLE INFO

Keywords:

Geodynamo

Polarity reversals

Magnetohydrodynamics

ABSTRACT

We study the process of magnetic reversals in the presence of a stably-stratified layer below the core-mantle boundary using direct numerical simulations of the incompressible magnetohydrodynamics equations under the Boussinesq approximation in a spherical shell. We show that the dipolar-multipolar transition shifts to larger Rayleigh numbers in the presence of a stably-stratified layer, and that the dipolar strength of the magnetic field at the core-mantle boundary increases due to the skin effect. By imposing an heterogeneous heat flux at the outer boundary, we break the equatorial symmetry of the flow, and show that different heat flux patterns can trigger different dynamo solutions, such as hemispheric dynamos and polarity reversals. Using kinematic dynamo simulations, we show that the stably-stratified layer leads to similar growth rates of the dipole and quadrupole components of the magnetic field, playing the role of a conducting boundary layer, favouring magnetic reversals, and a dynamics predicted by low-dimensional models.

# 1. Introduction

Planetary magnetic fields are generated and maintained through the dynamo process, mechanism through which kinetic energy is converted into magnetic energy, being thus able to sustain a large scale magnetic field (Moffatt, 1978; Landeau et al., 2022; Jones, 2011). This mechanism takes place in the core of planets, where the turbulent motion of an electrically conducting fluid is mostly driven by convection and strongly affected by rotation. In the solar system, most planets produce a self-sustained magnetic field, such as the Earth, Mercury and the gaseous giants. Due to differences in the internal composition and size of these planets, the properties of their magnetic fields vary, including their intensity, tilt angle relative to the rotation axis, degree of axial symmetry, and north-south symmetry, among others. At present, the Earth has a magnetic field with a strong dipolar component with tilting angle of about 11° and an intensity at its surface that varies between 0.25 and 0.65 G.

One of the most interesting phenomena of the geomagnetic field is its ability to reverse polarity in an irregular and chaotic manner (Lowrie and Kent, 2013). Over the past 80 million years, Earth's magnetic field has reversed on average with a frequency of approximately 4 Myr⁻¹, with individual reversal events typically lasting around 20,000 years. However, both the frequency and duration of reversals vary significantly: some intervals, known as superchrons, can span tens of millions of years without any reversals, while others show more frequent and irregular polarity changes. Understanding this process has been one of the main objectives of direct numerical simulations (DNS) of the geodynamo model (Braginsky and Roberts, 1995; Glatzmaier and Roberts, 1995; Christensen et al., 2001; Kutzner and Christensen, 2002; Li et al., 2002; Wicht and Olson, 2004; Takahashi et al., 2005; Driscoll and Olson, 2009; Olson et al., 2009; Amit et al., 2010; Olson et al., 2011;

Olson and Amit, 2014; Sheyko et al., 2016; Menu et al., 2020). Parametric studies on this problem suggested that the transition from stable dynamos dominated by a strong dipolar component to reversing solutions is controlled by the relative amplitude of inertia and Coriolis force (Kutzner and Christensen, 2002; Christensen, 2006). As the flow becomes more turbulent, the reversal frequency increases. More precisely, this transition has been proposed to be controlled by the local Rossby number defined as $Ro_ℓ = Ro\bar{e}_u/π$, with $Ro = U/(LΩ)$ the Rossby number, with $U$ the characteristic velocity, $L$ the characteristic length scale, $Ω$ the rotation rate, and $\bar{e}_u$ the characteristic spherical harmonic degree of the flow. Numerical simulations first suggested that the transition takes place around a critical value $Ro_ℓ ≈ 0.1$, with smaller values of $Ro_ℓ$ corresponding to stable dipolar solutions. The reversals found for values $Ro_ℓ > 0.1$ are typically non-dipolar, meaning that the dipole component is relatively small with respect to the quadrupole and higher order modes, which make these kind of reversals not Earth-like. This regime is usually known as multipolar dynamo. There has been since then an interest to find a regime of parameters in which Earth-like magnetic reversals are observed. In particular, simulations show that Earth-like magnetic field reversals occur only within a very narrow range of local Rossby number around 0.1, a degree of fine-tuning that seems unlikely to have been maintained over hundreds of millions of years in Earth's core. Furthermore, it was later observed that the local Rossby number fails to describe the dipolar-multipolar transition for some range of parameters, and other control parameters have been proposed, such as the ratio between kinetic and magnetic energies (Menu et al., 2020; Zaire et al., 2022; Frasson et al., 2025). In particular, Tassin et al. (2021) showed that reversing dynamos occur when this ratio is larger than $O(1)$, which is much larger than the typical values estimated for the Earth's core of $10^{-3}$. In the recent work of Jones and Tsang (2025), it was shown that it is possible to obtain Earth-like reversals while maintaining realistic kinetic-to-magnetic energy ratios by increasing the

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

nicolas.muller@lpp.polytechnique.fr (N.P. Müller)

ORCID(s):
magnetic Prandtl number together with a corresponding increase in the thermal Prandtl number, preserving realistic transport levels.

The exact nature of convective motion in Earth's outer core remains uncertain. Standard models typically assume that convection occurs throughout the entire region between the inner-core boundary (ICB) and the core-mantle boundary (CMB). However, recent seismic and geomagnetic measurements have questioned this assumption, suggesting the existence of a stably-stratified layer (SSL) below the CMB (Helffrich and Kaneshima, 2010, 2013; Buffett et al., 2016; Kaneshima, 2018). The origin of this layer may be either thermal or compositional in nature, potentially resulting from the accumulation of light elements below the CMB or the high thermal conductivity of the liquid metal relative to the heat flux at the CMB, leading to a sub-adiabatic temperature gradient (Pozzo et al., 2012). The existence of such layer in Earth's core remains a subject of debate, with its estimated size ranging from $H_s = 0$ km to $H_s = 450$ km, and stratification strengths ranging from $N = 0$ to $N = 50\Omega$, where $N$ is the buoyancy frequency (Christensen, 2018; Gastine et al., 2020). The presence of this layer has also been predicted in other planetary systems, such as Mercury, Saturn or Jupiter, often used to explain different large-scale structures of the magnetic field (Christensen and Aubert, 2006; Moore et al., 2018; Takahashi et al., 2019; Yan and Stanley, 2021; Guervilly, 2022). One of the main effects of a SSL below the CMB is that high-order modes of the magnetic field are smoothed out due to the skin effect, leading to a magnetic field dominated by low-order modes. The recent work of Aubert et al. (2025) proposed a mechanism that differs from the classical approaches based on the levels of inertia of the flow. The proposed kinematic mechanism suggests that in the presence of a stable top layer, there is a competition between the subsurface upwelling (generated in the convective region) and the surface circulation (generated in the stable top layer).

In addition to the possible presence of a SSL below the CMB, seismic measurements and laboratory experiments suggest that the heat flux at the CMB is not uniform, but it rather exhibits a complex spatial distribution (Su et al., 1994; Sumita and Olson, 2002). DNS of the geodynamo problem have shown that an heterogeneous heat flux at the CMB can modify the morphology of the magnetic field, trigger polarity reversals, and control their frequency (Glatzmaier et al., 1999; Olson et al., 2010; Gissinger et al., 2012; Christensen, 2018; Yan et al., 2023; Terra-Nova and Amit, 2024; Frasson et al., 2025). In particular, the non-axisymmetric structure of the geomagnetic field at the CMB are often associated to lateral heterogeneous heat flux.

In this work, we study the influence of a stably-stratified layer below the CMB on the geodynamo process, and its effects on polarity reversals. We characterise the dipolar-multipolar transition for different stably-stratified layers, and explore the effects of an axisymmetric heterogeneous heat flux at the CMB that breaks the equatorial symmetry of

<figure><img src="image_3.png" /><figcaption>**Figure 1:** Schematic picture of the static temperature gradient used to reproduce a mixed system, with a convective region and a stably-stratified layer described by (5).</figcaption></figure>

the flow. The goal is to understand how breaking equatorial symmetry influences the large-scale magnetic field and promotes the emergence of Earth-like polarity reversals. We finally use a low-dimensional model and kinematic dynamo simulations to provide a better interpretation of our results.

## 2. Model

We model the dynamics of an electrically conducting fluid using the incompressible magnetohydrodynamics (MHD) equations under the Boussinesq approximation in a spherical shell of inner radius $r_i$ and outer radius $r_o$. The aspect ratio is defined as $\chi = r_i/r_o$ and is fixed to 0.35 as the case of Earth. We express the length in units of the gap $L = r_o - r_i$, and the time in units of $L^2/v$, with $v$ the kinematic viscosity. The magnetic field $B$ is expressed in units of $\sqrt{\rho\mu\eta\Omega}$, with $\rho$ the mean fluid density, $\mu$ the magnetic permeability, $\eta$ the magnetic diffusivity and $\Omega$ the rotation rate. The velocity field $u$ is expressed in units of $v/L$ and the temperature $T$ in units of $q_o L/k$ with $q_o$ the mean heat flux at the CMB and $k$ the thermal conductivity. These units are used to obtain the dimensionless MHD equations (Braginsky and Roberts, 1995)

$$ \partial_t u + u \cdot \nabla u = -\nabla p + \nabla^2 u - \frac{2}{E} \hat{z} \times u + \frac{\text{Ra}}{E} \frac{r}{r_o} \theta \hat{r} + \frac{1}{\text{EPm}} (\nabla \times B) \times B, \quad (1) $$

$$ \partial_t B = \nabla \times (u \times B) + \frac{1}{\text{Pm}} \nabla^2 B, \quad (2) $$

$$ \partial_t \theta + u \cdot \nabla \theta = -u_r \frac{dT_s}{dr} + \frac{1}{\text{Pr}} \nabla^2 \theta, \quad (3) $$

$$ \nabla \cdot u = \nabla \cdot B = 0, \quad (4) $$

where $p$ the total pressure, and $T = T_s + \theta$ the total temperature with $T_s$ the static temperature profile and $\theta$ the temperature fluctuations. We define the dimensionless numbers using the standard convention: the Ekman number $E = v/(\Omega L^2)$, the thermal Prandtl number $Pr = v/\kappa$, the magnetic Prandtl number $Pm = v/\eta$ and the Rayleigh number $Ra = \alpha g_o q_o L^2/(\nu k \Omega)$ with $q_o$ the mean heat flux
<figure><img src="image_4.png" /><figcaption>**Figure 2:** (a) Time averaged dipolar strength $f_{dip}$ as a function of the size of the stably-stratified layer $H_s$ for Pm = 10 and E = 10⁻³ for different Ra. As the size of the layer $H_s$ increases, the dipolar component of the magnetic field at the CMB increases due to the skin effect. Error bars correspond to the standard deviation. (Right panels) Visualisation for Ra = 15Ra$_c$ of the radial (left) velocity and (right) magnetic fields for (b) a fully convective system $H_s = 0$, and (c) $H_s = 0.24L$.</figcaption></figure>

at the CMB, $k$ the thermal conductivity, $\alpha$ the thermal expansion coefficient, and the gravity acceleration defined as $g = g_0 r/r_0 \hat{r}$ with $g_0$ the value at the outer core. We use no-slip boundary conditions for the velocity field at both the ICB and CMB, constant temperature at the ICB, fixed heat flux $q_0$ at the CMB, a conducting inner core, and an insulating mantle.

We solve this system of equations using the PARODY-JA code (Dormy et al., 1998; Aubert et al., 2008), that uses a pseudo-spectral method in the poloidal and azimuthal directions, and a second-order finite-differences scheme for the radial direction, with a refined number of collocation points near the boundaries. The code implements a MPI parallelisation in the radial direction, and it is coupled to the open-source library SHTns (Schaeffer, 2013) that implements an openMP parallelisation in the angular directions to improve the performance of the spherical harmonic transforms. The resolution in the angular directions is defined such that the ratio between the maximum and minimum values of the magnetic energy spectrum is larger than 50. The resolution

in the radial direction is chosen so that there are at least 10 points in the Ekman layer at the boundaries.

We introduce an SSL in our system of equations in an ad-hoc manner by implementing a static temperature gradient background given by a piecewise function. In the convective region, this gradient is given by the solution of $\nabla^2 T_s = 0$ for $r_i < r < r_s$, while we propose a simple linear profile for the stably-stratified region $r_s < r < r_o$, leading to (Nakagawa, 2011, 2015; Couston et al., 2017; Christensen, 2018; Gastine et al., 2020)

$$ \frac{dT_s}{dr} = \begin{cases} -\frac{r_o^2}{r^2} & \text{for } r_i < r < r_s \\ \Gamma \frac{r-r_s}{H_s} - \frac{r_o^2 r_o - r}{r^2 H_s} & \text{for } r_s < r < r_o \end{cases} \qquad (5) $$

with $H_s = r_o - r_s$ the size of the SSL, and $\Gamma$ an extra free parameter that controls the maximum stratification strength as $N_{max}/\Omega = \sqrt{RaE\Gamma}$. We choose this value so that the maximum stratification strength is between 5 and 15, corresponding to a relatively strong stratification. We vary the size of the layer from $H_s = 0$ to $H_s = 0.24L$. The shape of the static temperature gradient is shown schematically in Fig. 1. Other choices of the static temperature gradient are possible, such as a continuous profile given by an hyperbolic tangent (Nakagawa, 2011, 2015; Takahashi et al., 2019). The drawback of these choices is that they introduce an extra parameter that determines the size of the transition between the convective and stably-stratified regions. This layer can also be modelled using double-diffusive convection, considering both thermal and compositional buoyancies, that leads to different regimes such as classical or finger convection (Monville et al., 2019; Guervilly, 2022). In this work, we restrict ourselves to a linear profile for the SSL below the CMB.

On top of the SSL, we explore the effects of an heterogeneous heat flux at the CMB. For this purpose, we impose $q_{CMB}(\theta, \phi) = q_0 + F(\theta, \phi)$ with $F(\theta, \phi)$ a function proportional to the spherical harmonics $Y_l^m$ of degree $l$ and order $m$. In particular, we study the effects of the axisymmetric patterns $Y_1^0$ or $Y_3^0$, and a linear combination of them, such that $F(\theta, \phi) = \bar{F}(\theta)$. We control the amplitude of this heat flux as

$$ \delta q = \frac{q_{\max} - q_{\min}}{2q_0}. \qquad (6) $$

### 3. Outputs

The code PARODY-JA solves the incompressible MHD Eqs. (1)-(4) using the poloidal-toroidal decomposition, where the total magnetic field $B$ is expressed in terms of two scalar potentials $B_{pol}$ and $B_{tor}$ such that $B = \nabla \times \nabla \times (B_{pol}r) + \nabla \times (B_{tor}r)$. The total magnetic field $B$ can be described as a linear combination of spherical harmonics $Y_l^m$. In particular, we can extract the axisymmetric dipole and quadrupole components of the magnetic field at the
<figure><img src="image_5.png" /><figcaption>**Figure 3:** Dipolar-multipolar transition for different sizes $H_s$ of the stably-stratified layer as a function of (a) the Rayleigh number and (b) the local Rossby number for $E = 10^{-3}$, $Pr = 1$, and $Pm = 10$. Circles indicate stable solutions, while square markers indicate reversing solutions.</figcaption></figure>

CMB from the poloidal and toroidal decomposition as

$$ 
\begin{aligned}
D &= B_{\text{pol}}^{\text{CMB}}(l=1, m=0), \\
Q &= B_{\text{pol}}^{\text{CMB}}(l=2, m=0).
\end{aligned}
\qquad (7) 
$$

Note that the toroidal component of the magnetic field is zero at the CMB due to the insulating mantle.

The code is expressed in units of the viscous time $t_v = L^2/\nu$, which is related to the magnetic diffusion time $t_\eta = L^2/\eta = Pmt_v$. To study magnetic field reversals, it is more appropriate to express the time in terms of the dipole diffusion time

$$ t_{\text{dip}} = \frac{r_o^2}{\pi^2 \eta} = \frac{Pm}{\pi^2 (1 - \chi)^2} t_v, \qquad (8) $$

with $\chi = 0.35$. The dipole diffusion time is typically $t_{\text{dip}} \sim 38.000$ years for Earth's core.

We define the total dimensionless kinetic and magnetic energies as

$$ E_{\text{kin}} = \frac{1}{2V} \int_V u^2 dV, \qquad (9) $$

$$ E_{\text{mag}} = \frac{1}{2V} \int_V \left( \frac{\mathbf{B}}{\sqrt{\text{EPm}}} \right)^2 dV, \qquad (10) $$

with $V$ the volume of the spherical shell. The amplitudes of the velocity and magnetic fields are quantified by the Reynolds number $Re = UL/\nu$ and the Elsasser number $\Lambda = B^2/(\rho\mu\eta\Omega)$ with $U$ and $B$ the volume-averaged RMS amplitudes of the velocity and magnetic fields, respectively. For the adimensionalisation used in our simulations, these numbers are computed as

$$ Re = \sqrt{2E_{\text{kin}}}, \qquad (11) $$

$$ \Lambda = 2PmE_{\text{mag}}, \qquad (12) $$

from which we can compute the magnetic Reynolds number as

$$ Rm = PmRe. \qquad (13) $$

To characterize the dipolar-multipolar transition, we define the relative dipole field strength at the outer boundary as

$$ f_{\text{dip}} = \frac{\sqrt{E_{\text{mag}}^{\text{CMB}}(l=1)}}{\sqrt{\sum_{l=1}^{12} E_{\text{mag}}^{\text{CMB}}(l)}} \qquad (14) $$

with $E_{\text{mag}}^{\text{CMB}}(l) = \sum_m E_{\text{mag}}^{\text{CMB}}(l, m)$ the magnetic energy spectrum evaluated at the CMB. The denominator is truncated at $l=12$ to compare with geomagnetic observations, that lead to a typical value for the Earth $f_{\text{dip}} = 0.68$ (Christensen, 2006; Menu et al., 2020). The transition is typically controlled by the local Rossby number

$$ R_o = \frac{\bar{\ell}_u R_o}{\pi}, \qquad (15) $$

with $\bar{\ell}_u = \sum_l l E_{\text{kin}}(l) / \sum_l E_{\text{kin}}(l)$ the characteristic spherical harmonic degree of the flow, and $R_o = ReE$ the Rossby number with $Re = \sqrt{2E_{\text{kin}}}$ the Reynolds number.

To quantify the equatorial symmetry of the magnetic field, we use the mean hemisphericity at the CMB

$$ H_{\text{CMB}} = \frac{\langle B_r^N \rangle^2 - \langle B_r^S \rangle^2}{\langle B_r^N \rangle^2 + \langle B_r^S \rangle^2}, \qquad (16) $$

with $B_r^N$ and $B_r^S$ the radial magnetic field evaluated at the CMB in the northern and southern hemispheres, respectively. For an equatorially symmetric system, the mean hemisphericity is zero, while if the magnetic field is concentrated in one hemisphere, it is $\pm 1$. We quantify the symmetry breaking for the flow with respect to the equator from the antisymmetric toroidal and poloidal components of the kinetic energy as

$$ E_{\text{kin}}^{\text{asym}} = \sum_{l+m=2n} E_{\text{kin}}^{\text{tor}}(l, m) + \sum_{l+m=2n+1} E_{\text{kin}}^{\text{pol}}(l, m). \qquad (17) $$
<caption>**Table 1**</caption>

Table of simulations homogeneous boundary conditions. We report the size of the stably-stratified layer $H_s$, the normalised Rayleigh number $Ra/Ra_c$, the kinetic energy $E_{kin}$, the magnetic energy $E_{mag}$, the ratio of kinetic and magnetic energies $E_{kin}/E_{mag}$, the magnetic Reynolds number $Rm$, the antisymmetric kinetic energy $E_k^{asym}$, the dipolar strength $f_{dip}$, whether there are reversals, and the final time of the simulations $t_{run}$ expressed in diffusive time units. All simulations correspond to $E = 10^{-3}$, $Pr = 1$, $Pm = 10$, and $r_i/r_o = 0.35$.

<table><thead><tr><th>H<sub>s</sub>/L</th><th>Ra/Ra<sub>c</sub></th><th>Λ</th><th>E<sub>kin</sub></th><th>E<sub>mag</sub></th><th>E<sub>kin</sub>/E<sub>mag</sub></th><th>Rm</th><th>Ra<sub>c</sub></th><th>E<sub>k</sub><sup>asym</sup></th><th>f<sub>dip</sub></th><th>reversals?</th><th>t<sub>run</sub></th></tr></thead><tbody><tr><td>0</td><td>5</td><td>16.36</td><td>94</td><td>818</td><td>0.12</td><td>137</td><td>47.8</td><td>0.16</td><td>0.68</td><td>no</td><td>4.2</td></tr><tr><td>0</td><td>10</td><td>24.12</td><td>247</td><td>1206</td><td>0.20</td><td>222</td><td>47.8</td><td>0.19</td><td>0.56</td><td>no</td><td>2</td></tr><tr><td>0</td><td>15</td><td>23.98</td><td>428</td><td>1198</td><td>0.36</td><td>292</td><td>47.8</td><td>0.21</td><td>0.49</td><td>no</td><td>2</td></tr><tr><td>0</td><td>20</td><td>23.04</td><td>625</td><td>1152</td><td>0.54</td><td>353</td><td>47.8</td><td>0.22</td><td>0.44</td><td>no</td><td>2</td></tr><tr><td>0</td><td>22</td><td>17.08</td><td>759</td><td>853</td><td>0.89</td><td>389</td><td>47.8</td><td>0.19</td><td>0.39</td><td>no</td><td>5</td></tr><tr><td>0</td><td>23</td><td>15.29</td><td>819</td><td>764</td><td>1.07</td><td>404</td><td>47.8</td><td>0.18</td><td>0.36</td><td>yes</td><td>3.75</td></tr><tr><td>0</td><td>25</td><td>12.89</td><td>929</td><td>644</td><td>1.44</td><td>431</td><td>47.8</td><td>0.18</td><td>0.30</td><td>yes</td><td>8</td></tr><tr><td>0</td><td>30</td><td>12.63</td><td>1164</td><td>631</td><td>1.84</td><td>482</td><td>47.8</td><td>0.19</td><td>0.20</td><td>yes</td><td>2</td></tr><tr><td>0.02</td><td>20</td><td>16.90</td><td>680</td><td>845</td><td>0.81</td><td>369</td><td>47.7</td><td>0.17</td><td>0.41</td><td>no</td><td>3</td></tr><tr><td>0.02</td><td>25</td><td>11.59</td><td>962</td><td>579</td><td>1.66</td><td>438</td><td>47.7</td><td>0.17</td><td>0.30</td><td>yes</td><td>6</td></tr><tr><td>0.04</td><td>15</td><td>22.09</td><td>417</td><td>1104</td><td>0.38</td><td>289</td><td>48.0</td><td>0.17</td><td>0.59</td><td>no</td><td>4</td></tr><tr><td>0.04</td><td>20</td><td>17.25</td><td>677</td><td>862</td><td>0.79</td><td>368</td><td>48.0</td><td>0.16</td><td>0.52</td><td>no</td><td>3</td></tr><tr><td>0.04</td><td>25</td><td>11.59</td><td>995</td><td>579</td><td>1.72</td><td>446</td><td>48.0</td><td>0.15</td><td>0.41</td><td>yes</td><td>9</td></tr><tr><td>0.04</td><td>30</td><td>11.96</td><td>1228</td><td>598</td><td>2.05</td><td>495</td><td>48.0</td><td>0.16</td><td>0.37</td><td>yes</td><td>3</td></tr><tr><td>0.09</td><td>15</td><td>22.48</td><td>373</td><td>1124</td><td>0.33</td><td>273</td><td>50.22</td><td>0.18</td><td>0.80</td><td>no</td><td>4</td></tr><tr><td>0.09</td><td>20</td><td>16.72</td><td>599</td><td>836</td><td>0.72</td><td>346</td><td>50.22</td><td>0.16</td><td>0.77</td><td>no</td><td>2.5</td></tr><tr><td>0.09</td><td>25</td><td>14.05</td><td>841</td><td>702</td><td>1.20</td><td>410</td><td>50.22</td><td>0.16</td><td>0.74</td><td>no</td><td>3</td></tr><tr><td>0.09</td><td>30</td><td>9.68</td><td>1175</td><td>484</td><td>2.43</td><td>484</td><td>50.22</td><td>0.16</td><td>0.53</td><td>yes</td><td>6.75</td></tr><tr><td>0.14</td><td>15</td><td>21.35</td><td>371</td><td>1067</td><td>0.35</td><td>272</td><td>56.81</td><td>0.18</td><td>0.87</td><td>no</td><td>2.49</td></tr><tr><td>0.14</td><td>20</td><td>17.29</td><td>596</td><td>864</td><td>0.69</td><td>345</td><td>56.81</td><td>0.17</td><td>0.84</td><td>no</td><td>2.54</td></tr><tr><td>0.14</td><td>25</td><td>15.21</td><td>841</td><td>760</td><td>1.11</td><td>410</td><td>56.81</td><td>0.16</td><td>0.83</td><td>no</td><td>2.52</td></tr><tr><td>0.14</td><td>30</td><td>7.76</td><td>1257</td><td>387</td><td>3.24</td><td>501</td><td>56.81</td><td>0.14</td><td>0.56</td><td>yes</td><td>2.51</td></tr><tr><td>0.14</td><td>35</td><td>8.74</td><td>1435</td><td>437</td><td>3.28</td><td>535</td><td>56.81</td><td>0.17</td><td>0.31</td><td>yes</td><td>1.67</td></tr><tr><td>0.24</td><td>15</td><td>17.43</td><td>374</td><td>871</td><td>0.43</td><td>273</td><td>66.4</td><td>0.18</td><td>0.93</td><td>no</td><td>3</td></tr><tr><td>0.24</td><td>20</td><td>14.28</td><td>584</td><td>713</td><td>0.82</td><td>341</td><td>66.4</td><td>0.17</td><td>0.92</td><td>no</td><td>1.27</td></tr><tr><td>0.24</td><td>22</td><td>13.23</td><td>681</td><td>661</td><td>1.03</td><td>369</td><td>66.4</td><td>0.17</td><td>0.92</td><td>no</td><td>3</td></tr><tr><td>0.24</td><td>25</td><td>9.25</td><td>927</td><td>462</td><td>2.01</td><td>430</td><td>66.4</td><td>0.15</td><td>0.89</td><td>no</td><td>7.56</td></tr><tr><td>0.24</td><td>26</td><td>9.35</td><td>979</td><td>467</td><td>2.10</td><td>442</td><td>66.4</td><td>0.15</td><td>0.79</td><td>yes</td><td>5.28</td></tr><tr><td>0.24</td><td>27.5</td><td>7.32</td><td>1100</td><td>366</td><td>3.00</td><td>469</td><td>66.4</td><td>0.15</td><td>0.56</td><td>yes</td><td>3.56</td></tr><tr><td>0.24</td><td>30</td><td>7.47</td><td>1212</td><td>373</td><td>3.25</td><td>492</td><td>66.4</td><td>0.16</td><td>0.39</td><td>yes</td><td>3</td></tr><tr><td>0.24</td><td>35</td><td>9.75</td><td>1342</td><td>487</td><td>2.75</td><td>518</td><td>66.4</td><td>0.19</td><td>0.34</td><td>yes</td><td>2.25</td></tr><tr><td>0.24</td><td>40</td><td>10.72</td><td>1508</td><td>536</td><td>2.81</td><td>549</td><td>66.4</td><td>0.20</td><td>0.33</td><td>yes</td><td>1.5</td></tr></tbody></table>

# 4. Results

## 4.1. Parametric study

We perform a parametric study of the geodynamo by fixing $Pr = 1$, $E = 10^{-3}$, $Pm = 10$, the aspect ratio of the system $\chi = 0.35$, and varying $Ra$, and the size of the SSL $H_s = r_o - r_s$ keeping the stratification strength between $N_{max}/\Omega \in (5, 15)$. Some simulations are initialised with a fluid at rest, imposing a small temperature and magnetic perturbations, while others are restarted from runs with similar parameters. The relative large values of $E$ and $Pm$ are chosen to enable an extensive parametric study of the system, allowing us to identify different regimes of geodynamo behaviour. The input and output quantities of the simulations are shown in Table 1.

The effects of the SSL in the magnetic field morphology are shown in Fig. 2. The dipolar strength of the magnetic field at the CMB $f_{dip}$ defined in Eq. (14) increases with the

size of the SSL for a fixed value of $Ra/Ra_c$, where $Ra_c$ is the numerically determined critical Rayleigh number for the onset of convection. This is a consequence of the skin effect, which attenuates the contribution of the more spatially structured modes of the magnetic field. The visualisation of the radial velocity and magnetic fields in Figs. 2(b)-(c) show how they are attenuated in the presence of the SSL.

Stable dipolar dynamos are typically characterised by a dipolar strength $f_{dip} > 0.5$, whereas lower values correspond to multipolar solutions, in which the dipole component is no longer dominant. The transition between these two dynamo solutions remains a major unresolved issue in our understanding of the geodynamo. Figure 3 (a) shows the dipolar-multipolar transition as a function of the Rayleigh number. In the case of a fully convective system ($H_s = 0$), the transition takes place at around $Ra/Ra_c = 15$. The introduction of the SSL shifts this transition to larger values, around $Ra/Ra_c = 27.5$ for $H_s = 0.24L$, and
**Table 2**

Table of simulations with heterogeneous boundary conditions. We report the same parameters as in Table 1. The heat flux pattern at the CMB is determined by the spherical harmonic $Y_l^m$, and the amplitude of the heterogeneity is given by $\delta q$ defined in Eq. (6). All simulations correspond to $E = 10^{-3}$, Pr = 1, Pm = 10, Ra = 25Ra_c, and a SSL of size $H_s = 0.24L$ and $N_{max}/\Omega \approx 10$.

<table><thead><tr><th>Y<sub>l</sub><sup>m</sup></th><th>δq (%)</th><th>Λ</th><th>E<sub>kin</sub></th><th>E<sub>mag</sub></th><th>E<sub>kin</sub>/E<sub>mag</sub></th><th>Rm</th><th>E<sub>kin</sub><sup>asym</sup></th><th>f<sub>dip</sub></th><th>Reversals?</th><th>t<sub>run</sub></th></tr></thead><tbody><tr><td>Y<sub>l</sub><sup>0</sup></td><td>-3.5</td><td>10.61</td><td>888</td><td>530</td><td>1.67</td><td>420</td><td>0.17</td><td>0.76</td><td>no</td><td>3.10</td></tr><tr><td>Y<sub>l</sub><sup>0</sup></td><td>0.0</td><td>9.59</td><td>916</td><td>480</td><td>1.91</td><td>427</td><td>0.15</td><td>0.90</td><td>no</td><td>7.56</td></tr><tr><td>Y<sub>l</sub><sup>0</sup></td><td>0.17</td><td>9.36</td><td>932</td><td>468</td><td>1.99</td><td>431</td><td>0.14</td><td>0.82</td><td>no</td><td>3.25</td></tr><tr><td>Y<sub>l</sub><sup>0</sup></td><td>0.35</td><td>7.77</td><td>985</td><td>389</td><td>2.53</td><td>443</td><td>0.13</td><td>0.81</td><td>no</td><td>4.45</td></tr><tr><td>Y<sub>l</sub><sup>0</sup></td><td>1.05</td><td>9.12</td><td>939</td><td>456</td><td>2.06</td><td>432</td><td>0.14</td><td>0.82</td><td>no</td><td>3.15</td></tr><tr><td>Y<sub>l</sub><sup>0</sup></td><td>1.75</td><td>10.71</td><td>887</td><td>536</td><td>1.66</td><td>420</td><td>0.16</td><td>0.86</td><td>no</td><td>2.00</td></tr><tr><td>Y<sub>l</sub><sup>0</sup></td><td>2.1</td><td>11.18</td><td>876</td><td>559</td><td>1.57</td><td>418</td><td>0.17</td><td>0.81</td><td>no</td><td>4.00</td></tr><tr><td>Y<sub>l</sub><sup>0</sup></td><td>3.5</td><td>12.23</td><td>832</td><td>612</td><td>1.36</td><td>407</td><td>0.19</td><td>0.78</td><td>no</td><td>6.00</td></tr><tr><td>Y<sub>l</sub><sup>0</sup></td><td>8.75</td><td>12.85</td><td>828</td><td>643</td><td>1.29</td><td>406</td><td>0.24</td><td>0.56</td><td>no</td><td>1.50</td></tr><tr><td>Y<sub>l</sub><sup>0</sup></td><td>-5.3</td><td>11.20</td><td>880</td><td>560</td><td>1.57</td><td>418</td><td>0.17</td><td>0.86</td><td>yes</td><td>9.06</td></tr><tr><td>Y<sub>l</sub><sup>0</sup></td><td>0.0</td><td>9.42</td><td>922</td><td>471</td><td>1.96</td><td>428</td><td>0.15</td><td>0.89</td><td>no</td><td>7.56</td></tr><tr><td>Y<sub>l</sub><sup>0</sup></td><td>0.53</td><td>10.33</td><td>884</td><td>516</td><td>1.71</td><td>420</td><td>0.16</td><td>0.92</td><td>no</td><td>5.45</td></tr><tr><td>Y<sub>l</sub><sup>0</sup></td><td>2.65</td><td>10.72</td><td>867</td><td>536</td><td>1.64</td><td>415</td><td>0.17</td><td>0.91</td><td>no</td><td>13.44</td></tr><tr><td>Y<sub>l</sub><sup>0</sup></td><td>5.3</td><td>9.82</td><td>908</td><td>491</td><td>1.85</td><td>425</td><td>0.16</td><td>0.83</td><td>yes</td><td>29.50</td></tr><tr><td>Y<sub>l</sub><sup>0</sup> + Y<sub>l</sub><sup>0</sup></td><td>1.75</td><td>8.54</td><td>957</td><td>427</td><td>2.24</td><td>438</td><td>0.14</td><td>0.75</td><td>yes</td><td>27.65</td></tr></tbody></table>

makes the dipolar-multipolar transition more abrupt. Figure 3 (b) shows the same transition as a function of the local Rossby number defined in Eq. (15). The critical value of $R_o$ for the transition depends on the size of the layer. Zaire et al. (2022) proposed that the transition between stable and reversing dynamos is controlled by the ratio between kinetic and magnetic energies, that we report in Table 1. Their study in fully convective systems suggests that the transition takes place at around $E_{kin}/E_{mag} \approx 0.7$. In our case, as the size of the stably-stratified layer increases, this transition seems to shift to larger values, being around 2 for $H_s = 0.24L$.

## 4.2. Effects of an heterogeneous heat flux at the CMB

Having established the influence of the SSL on the dipolar-multipolar transition, we now turn to exploring the effects of an imposed axisymmetric heterogeneous heat flux at the CMB. We choose the case with a large SSL with $H_s = 0.24L$ and a value of $Ra/Ra_c = 25$ such that it is close to the transition with a stable dipolar solution. The input and output quantities of these simulations are summarised in Table 2. We start by imposing the $Y_l^0$ pattern at the CMB, and vary the amplitude of the symmetry breaking $\delta q$. We verified that when imposing a heat flux at the CMB lower than 10%, the equatorial symmetry of the flow does not change significantly. Figure 4 shows the phase diagram of the axial dipole *D* and quadrupole *Q* components of the magnetic field at the CMB, defined in Eq. (7), for different values of $\delta q$ in the presence of an SSL. In the case of no symmetry breaking (black curve), the magnetic field fluctuates around a stable state with a strong dipole component, and a quadrupole that fluctuates around zero. As the amplitude of the heterogeneous heat flux increases (red curve), the quadrupole is enhanced, and the magnetic field fluctuates around a mixed state, and does not display reversals. A movie of the simulation with $\delta q = 8.5\%$ can be found

in the Supplementary Material *sup*. Moreover, when the amplitude of the heterogeneous heat flux is negative (blue curve), the quadrupole also changes its sign. The straight line connecting to the origin is to indicate the angle between the dipole and quadrupole components. This enhancement of the quadrupole component is associated with an increase in the hemisphericity of the magnetic field, defined in Eq. (16) and shown in Fig. 4-(b). It is remarkable that breaking the heat flux symmetry by only a few percent - leading to a minimal modification of the velocity field ($E_{kin}^{asym}/E_{kin}^{tot} < 20\%$) - can result in a strongly hemispherical magnetic field, in line with the results of an $\alpha^2$ model of dynamo (Gallet and Pétrélis, 2009).

We also investigate the effects of a more complex heterogeneous heat flux pattern at the CMB given by the spherical harmonic $Y_3^0$. This choice is motivated primarily by theoretical and methodological considerations: $Y_3^0$ provides a simple axisymmetric pattern with increased latitudinal complexity compared to $Y_1^0$, allowing us to investigate the effects of higher-degree axisymmetric heterogeneity. The results of the simulations are shown in Fig. 5. As the amplitude of the heterogeneous heat flux increases, the system eventually reaches a regime in which there are magnetic reversals, switching from a $B \rightarrow -B$ solution, where both the dipole and quadrupole components change sign (red curve). Arrows indicate the typical trajectories of reversals. The time series of the dipole and quadrupole for $\delta q = 5.3\%$ is shown in Fig. 5 (b). It is found systematically that during the reversal process, the dipole reverses first, and the quadrupole follows slightly later. A movie of this simulation can be found in the Supplementary Material *sup*. Reversals follow the opposite trajectory when the amplitude of the heterogeneous heat flux is negative (blue curve). A long evolution of the dipolar angle and the dipolar strength $f_{dip}$ for $\delta q = 5.3\%$ is shown in Fig. 6. In the period of 120 dipole diffusion times, we observe several non-periodic chaotic reversals of
<figure><img src="image_6.png" /></figure>

**Figure 4:** (a) Evolution of the magnetic field at the CMB in the D-Q phase space for different amplitudes of the heterogeneous heat flux $\delta q$ for the simplest pattern $Y_1^0$. The points correspond to a portion of the time series (50$\tau_v$) taken from the statistically steady state of the dynamo solution. (b) Time averaged hemisphericity, defined as in Eq. (16), as a function of $\delta q$. All simulations correspond to $E = 10^{-3}$, $Pr = 1$, $Pm = 10$, and $Ra/Ra_c = 25$. Error bars indicate the standard deviation.

the magnetic field. During the reversal, the dipolar strength decreases from $f_{dip} \approx 0.9$ to a value $f_{dip} < 0.5$, consistent with the reduction of the dipole component during this event shown in Fig. 5(b). The high value of $f_{dip}$ between reversals can be attributed to the presence of the SSL, allowing for a strong dipolar magnetic field that reverses, in contrast with the magnetic reversals observed in the absence of a layer in the multipolar regime, dominated by high order modes. We also verified the robustness of these solutions by introducing small variations (up to 10%) in both the Rayleigh number and the amplitude of the heterogeneous heat-flux pattern $Y_3^0$. In these cases, we recover reversals with similar properties to those shown in Figs. 5 and 6.

We can also consider a linear combination of the two patterns $Y_1^0$ and $Y_3^0$ as a heterogeneous heat flux at the CMB. Figure 7 shows the evolution of the magnetic field in the D-Q phase space for a linear combination of these two patterns, both with the same amplitude and a total flux variation $\delta q = 1.75\%$. A movie of this simulation can be found in the Supplementary Material *sup*. This heterogeneous heat flux leads to polarity reversals of the magnetic field, with a trajectory that is qualitatively similar to the one described by $Y_3^0$, but here the time delay between the reversal of the

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

**Figure 5:** (a) Evolution of the magnetic field at the CMB in the D-Q phase space for different amplitudes of the heterogeneous heat flux $\delta q$ for the pattern $Y_3^0$. Arrows indicate the orientation of the typical trajectories of reversals. (b) Time series of $D$ and $Q$ for the case $\delta q = 5.3\%$. All simulations correspond to $E = 10^{-3}$, $Pr = 1$, $Pm = 10$, and $Ra/Ra_c = 25$.

dipole and quadrupole is larger, and the trajectories deviate significantly from the origin. The arrows in the figure show the direction of the reversal. The phenomenology observed here differs from the one described in the absence of a SSL, where both the dipole and quadrupole components fluctuate around the origin, as shown in the inset of Fig. 7(a). We performed also a set of simulations (not shown here) to explore if it is possible to recover reversals in the absence of a SSL and for $Ra/Ra_c = 5$. For high amplitudes of $\delta q$, we do not recover dynamo solutions as a consequence of a high symmetry breaking of the flow. For low amplitudes of $\delta q$ we do observe an enhancement of the quadrupolar component of the magnetic field. In the next section, we provide an interpretation of this phenomenology based on a low-dimensional model and the difference of the dipole and quadrupole growth rates.

# 5. Discussion

The phenomenology presented in the previous section, involving the emergence of hemispheric dynamos and polarity reversals due to an equatorial symmetry breaking, can be interpreted using low-dimensional models. These models consider that the axial dipole $D$ and axial quadrupole $Q$ are the two dominant components of the magnetic field, which
<figure><img src="image_8.png" /><figcaption>**Figure 6:** Polarity reversals as a function of dipolar diffusion time for a symmetry breaking $Y_3^0$ with an amplitude $\delta q = 5.3\%$, $E = 10^{-3}$, $Pm = 10$, $Pr = 1$ and $Ra = 25Ra_c$ in the presence of a stably-stratified layer.</figcaption></figure>

interact in a non-linear manner coupled by the hydrodynamic properties of the system. Taking into account the $B \rightarrow -B$ invariance of the MHD equations, this model is written as (Pétrélis and Fauve, 2008; Pétrélis et al., 2009)

$$ \begin{aligned} \dot{D} &= \sigma_D D + \alpha Q + C_{11} D^3 + C_{12} D^2 Q + C_{13} D Q^2 + C_{14} Q^3 \\ \dot{Q} &= \sigma_Q Q + \beta D + C_{21} D^3 + C_{22} D^2 Q + C_{23} D Q^2 + C_{24} Q^3 \end{aligned} \quad (18) $$

where all the parameters $C_{ij}$, $\alpha$, $\beta$ and $\sigma_{D,Q}$ are real numbers. In particular, the coefficients $\sigma_D$ and $\sigma_Q$ determine the dipole and quadrupole growth rates in the limit in which these two modes are decoupled and the non-linear interactions are negligible. When $\sigma_D$ and $\sigma_Q$ are similar, an equatorial symmetry breaking will result either in periodic polarity reversals or in hemispheric dynamos (Gallet and Pétrélis, 2009). In the full model, stochastic term is introduced to account for the effects of turbulence (Pétrélis and Fauve, 2008). Alternatively, Eqs (18) can be coupled to a third mode modelling the asymmetric velocity mode generated by the heterogeneous heat flux, resulting in a deterministic model for chaotic inversions (Gissinger, 2012).

<figure><img src="image_9.png" /><figcaption>**Figure 7:** Evolution of the dipole and quadrupole magnetic field at the CMB for a linear combination of the heterogeneous heat flux patterns $Y_1^0 + Y_3^0$ with an amplitude $\delta q = 1.75\%$. The arrows indicate the orientation of reversals. The inset shows a typical simulation with reversals in a fully convective domain. All simulations correspond to $E = 10^{-3}$, $Pr = 1$, $Pm = 10$, and $Ra/Ra_c = 25$.</figcaption></figure>

It has been shown that such low-dimensional models based on dipole-quadrupole interactions can reproduce fairly well the magnetic field reversals observed on Earth (Pétrélis et al., 2009, 2011) and in experimental dynamos (Berhanu et al., 2007; Gissinger et al., 2010). However, this behaviour has never been observed, so far, in numerical modelling of the geodynamo, even in the presence of a heterogeneous heat flux (Gissinger et al., 2012). Interestingly, it is therefore the presence of an SSL in the present simulations that gives rise to this long-anticipated phenomenology. This could be a consequence of the attenuation of fast modes of the magnetic field due to the skin effect in the SSL. A more appealing interpretation involves the relative growth rates of the dipole and quadrupole components and their non-linear interaction. To understand the behaviour of these quantities, we perform DNS of the kinematic dynamo model, which corresponds to the same system of equations (1)-(4) eliminating the Lorentz force from Eq. (1). We point out that the velocity field still evolves in time under the same boundary conditions (including the heterogeneous CMB heat flux) as in the self-consistent simulations, but now independently of the magnetic field. The magnetic field therefore responds passively to the temporal fluctuations of this flow. We first focus on the hydrodynamic properties of the system, in particular on the
<figure><img src="image_10.png" /><figcaption>**Figure 8:** Kinematic dynamo simulations for different sizes of the stably-stratified layer $H_s$ and different Rayleigh numbers, for fixed $E = 10^{-3}$, $Pr = 1$, and $Pm = 10$. Time-averaged (a) relative antisymmetric kinetic energy and (b) relative growth rate difference between the dipole and quadrupole. Error bars represent the standard deviation.</figcaption></figure>

equatorial symmetry of the flow as a function of the Rayleigh number. Figure 8(a) shows the relative antisymmetric kinetic energy for different values of the size of the SSL defined in Eq. (17). For low values of $Ra/Ra_c$, this quantity is zero indicating that the flow is symmetric with respect to the equator. In this regime, the dipole and quadrupole components are independent as they are decoupled. At around $Ra = 10Ra_c$, the antisymmetric kinetic energy undergoes a transition and the equatorial symmetry of the flow is broken. When this happens, $D$ and $Q$ are coupled. When the two modes have similar growth rates, this favours magnetic reversals. We compute the dipole $\sigma_D$ and quadrupole $\sigma_Q$ growth rates as

$$D(t) = D_0 \exp(\sigma_D t), \qquad (19)$$

$$Q(t) = Q_0 \exp(\sigma_Q t). \qquad (20)$$

It is interesting to note that this transition occurs simultaneously with the dipolar-multipolar transition in the fully convective system, and explains why magnetic reversals are obtained in this regime. This result is consistent with Garcia et al. (2017), where it was shown that the loss of dipolarity is associated to a breaking of the equatorial symmetry of the flow.

More interestingly, we now focus on the effect of the SSL on the kinematic growth rates of the dipole and quadrupole modes below the symmetry breaking threshold, in the

regime where the two linear growth rates remain decoupled. Figure 8(b) shows the dipole and quadrupole growth rates for $Ra = 7Ra_c$ for different values of $H_s$. When the size of the SSL increases, the relative difference between the kinematic growth rates of the dipole and quadrupolar modes decreases. This behaviour can be understood with the prediction by Proctor (1977), who showed that in $\alpha^2$-dynamos surrounded by a perfectly conducting, motionless outer layer, a degeneracy arises between the growth rates of the dipole and quadrupole components. In this context, the SSL beneath the core-mantle boundary can be interpreted as an effective conducting boundary, leading to a similar degeneracy as its thickness increases (Favier and Proctor, 2013). These conditions strongly favour the emergence of the low-dimensional dynamics described above: as the dipole and quadrupole growth rates $\sigma_D$ and $\sigma_Q$ approach one another, even a weak symmetry breaking induces a spectacularly strong coupling between the two modes, enhancing the dynamics predicted by low-dimensional models - including both chaotic reversals and the development of hemispherical magnetic fields.

## 6. Conclusions

In this work, we studied the effects of a stably-stratified layer (SSL) below the core-mantle boundary (CMB) on the morphology of the magnetic field, with particular focus on polarity reversals. Using direct numerical simulations of the incompressible magnetohydrodynamics equations in a rotating spherical shell driven by convection, we showed that the SSL enhances the dipolar strength of the magnetic field at the CMB due to the skin effect, that attenuates rapidly fluctuating fields. As a result, the transition from dipolar to multipolar dynamos occurs at higher values of the Rayleigh number when the SSL is present, and the transition becomes more abrupt.

By imposing both an SSL and a heterogeneous heat flux at the CMB to break the equatorial symmetry of the flow, we observed a variety of dynamo regimes, including hemispheric dynamos and polarity reversals. Hemispheric dynamos arise from an increase in the quadrupolar component of the magnetic field, which is due to the imposed heterogeneous heat flux and the resulting weak equatorial symmetry breaking of the flow. For the parameters studied in this work, we found that the simplest heat flux pattern, described by the spherical harmonic $Y_1^0$, controls the hemisphericity of the magnetic field. In contrast, a more complex pattern $Y_3^0$ is responsible for triggering polarity reversals. When a linear combination of these two patterns is imposed, reversals follow a clear trajectory in the dipole-quadrupole ($D - Q$) phase space, with the dipole reversing systematically before the quadrupole. This reversal behaviour differs from that observed in fully convective systems within the multipolar regime, where both the dipole and quadrupole components fluctuate without following a well-defined trajectory in the $D - Q$ plane. We also highlight that between the reversals, the dipolar strength of the magnetic field takes relatively
large values $f_{\text{dip}} \approx 0.9$, that are typically not observed in simulations.

We interpret these results and the role of the SSL by means of kinematic dynamo simulations. In fully convective dynamos, we associate the dipolar-multipolar transition with a hydrodynamical instability - controlled by the Rayleigh number - that breaks the equatorial symmetry of the flow. When the flow remains symmetric with respect to the equatorial plane, the dipole and the quadrupole become decoupled, and the growth rates of the dipole $\sigma_D$ and the quadrupole $\sigma_Q$ approach each other as the thickness of the SSL increases. This behaviour is consistent with the prediction by Proctor (1977), who showed that in $a^2$-dynamos surrounded by a perfect conductor, a degeneracy arises between the dipolar and quadrupolar growth rates.

The mechanism based on the interaction of the dipole and quadrupole components of the magnetic field provides an explanation for the mechanism of magnetic reversals observed in the geodynamo. It explains the observations of Garcia et al. (2017), where it was shown that the loss of dipolarity is associated to a breaking of the equatorial symmetry of the flow. Recent studies of Aubert et al. (2025) and Jones and Tsang (2025) provide other alternative for which the geodynamo model can reproduce Earth-like magnetic reversals. In particular, Aubert et al. (2025) provide a kinematic mechanism controlled by the relative strength of subsurface upwelling and horizontal circulation at the core surface. In Jones and Tsang (2025), the authors find a regime of parameters in which the magnetic field reverses even in the low inertia limit by fixing the ratio between Pm and Pr. Our results propose a different mechanism based on a low-dimensional model used to explain reversals in dynamo experiments Pétrélis et al. (2009); Pétrélis and Fauve (2008). The simulations showing reversing solutions in this work operate close to the dipolar-multipolar transition, but we expect this mechanism to apply also for low-inertia systems.

This interpretation of the SSL, and its critical role in reducing the difference between dipolar and quadrupolar growth rates ($\sigma_D - \sigma_Q$), provides an elegant resolution to a long-standing question regarding the robustness of magnetic field reversals. Indeed, despite considerable variability in core conditions - such as inner-core growth, fluctuations in the heat flux, and other long-term evolutionary processes - the presence of a substantial SSL acts as a stabilizing factor. By maintaining the near-degeneracy of $\sigma_D$ and $\sigma_Q$, it ensures that the low-dimensional dynamics driven by the dipole-quadrupole interaction remains robust, even in the presence of other complex small scale modes. Consequently, even small symmetry breaking effects become sufficient to sustain the observed reversal behaviour over geological timescales, making Earth's dynamo remarkably resilient to changes in its internal properties. Finally, we remark the results shown in this work were obtained for a relatively high value of the Ekman number of $E = 10^{-3}$ and the magnetic Prandtl number $Pm = 10$, which are quite far from realistic values for the Earth of $E = 10^{-15}$ and $Pm = 10^{-6}$. Further studies could

investigate the dependence of this phenomenology with the Ekman number and also the magnetic Prandtl number.

## Acknowledgments

We acknowledge support from Agence Nationale de la Recherche (Grant No. ANR-19-CE31-0019-01). This project was provided with computer and storage resources by GENCI at CINES thanks to the grant A0170515642 on the supercomputer Adastra GENOA. This work was also granted access to the HPC resources of MesoPSL funded by the Region Île-de-France.

## References

See Supplemental Material at URL-will-be-inserted-by-publisher for visualizations of the numerical simulations.

Amit, H., Leonhardt, R., Wicht, J., 2010. Polarity Reversals from Paleomagnetic Observations and Numerical Dynamo Simulations. Space Science Reviews 155, 293–335. doi:10.1007/s11214-010-9695-2.

Aubert, J., Aurnou, J., Wicht, J., 2008. The magnetic structure of convection-driven numerical dynamos. Geophysical Journal International 172, 945–956. doi:10.1111/j.1365-246X.2007.03693.x.

Aubert, J., Landeau, M., Fournier, A., Gastine, T., 2025. Core-surface kinematic control of polarity reversals in advanced geodynamo simulations. Physics of the Earth and Planetary Interiors 364, 107365. doi:10.1016/j.pepi.2025.107365.

Berhanu, M., Monchaux, R., Fauve, S., Mordant, N., Pétrélis, F., Chiffaudel, A., Daviaud, F., Dubrulle, B., Marié, L., Ravelet, F., Bourgoin, M., Odier, P., Pinton, J.F., Volk, R., 2007. Magnetic field reversals in an experimental turbulent dynamo. Europhysics Letters (EPL) 77, 59001. doi:10.1209/0295-5075/77/59001.

Braginsky, S.I., Roberts, P.H., 1995. Equations governing convection in earth's core and the geodynamo. Geophysical & Astrophysical Fluid Dynamics 79, 1–97. doi:10.1080/03091929508228992.

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. Geophysical Journal International 204, 1789–1800. doi:10.1093/gji/ggy552.

Christensen, U., 2018. Geodynamo models with a stable layer and heterogeneous heat flow at the top of the core. Geophysical Journal International 215, 1338–1351. doi:10.1093/gji/ggy352.

Christensen, U., Aubert, J., Cardin, P., Dormy, E., Gibbons, S., Glatzmaier, G., Grote, E., Honkura, Y., Jones, C., Kono, M., Matsushima, M., Sakuraba, A., Takahashi, F., Tilgner, A., Wicht, J., Zhang, K., 2001. A numerical dynamo benchmark. Physics of the Earth and Planetary Interiors 128, 25–34. doi:10.1016/S0031-9201(01)00275-8.

Christensen, U.R., 2006. A deep dynamo generating Mercury's magnetic field. Nature 444, 1056–1058. doi:10.1038/nature05342.

Christensen, U.R., Aubert, J., 2006. Scaling properties of convection-driven dynamos in rotating spherical shells and application to planetary magnetic fields. Geophysical Journal International 166, 97–114. doi:10.1111/j.1365-246X.2006.03009.x.

Couston, L.A., Lecoanet, D., Favier, B., Le Bars, M., 2017. Dynamics of mixed convective-stably-stratified fluids. Physical Review Fluids 2, 094804. doi:10.1103/PhysRevFluids.2.094804.

Dormy, E., Cardin, P., Jault, D., 1998. MHD flow in a slightly differentially rotating spherical shell, with conducting inner core, in a dipolar magnetic field. Earth and Planetary Science Letters 160, 15–30. doi:10.1016/S0012-821X(98)00078-8.

Driscoll, P., Olson, P., 2009. Effects of buoyancy and rotation on the polarity reversal frequency of gravitationally driven numerical dynamos. Geophysical Journal International 178, 1337–1350. doi:10.1111/j.1365-246X.2009.04234.x.
Favier, B., Proctor, M.R.E., 2013. Growth rate degeneracies in kinematic dynamos. Physical Review E 88, 031001. doi:10.1103/PhysRevE.88.031001.

Frasson, T., Schaeffer, N., Nataf, H.C., Labrosse, S., 2025. Geomagnetic dipole stability and zonal flows controlled by mantle heat flux heterogeneities. Geophysical Journal International 240, 1481–1504. doi:10.1093/gji/ggae457.

Gallet, B., Pétrélis, F., 2009. From reversing to hemispherical dynamos. Physical Review E 80, 035302. doi:10.1103/PhysRevE.80.035302.

Garcia, F., Oruba, L., Dormy, E., 2017. Equatorial symmetry breaking and the loss of dipolarity in rapidly rotating dynamos. Geophysical & Astrophysical Fluid Dynamics 111, 380–393. doi:10.1080/03091929.2017.1347785.

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

Gissinger, C., 2012. A new deterministic model for chaotic reversals. The European Physical Journal B 85, 137. doi:10.1140/epjb/e2012-20799-5.

Gissinger, C., Dormy, E., Fauve, S., 2010. Morphology of field reversals in turbulent dynamos. EPL (Europhysics Letters) 90, 49001. doi:10.1209/0295-5075/90/49001.

Gissinger, C., Petitdemange, L., Schrinner, M., Dormy, E., 2012. Bistability between Equatorial and Axial Dipoles during Magnetic Field Reversals. Physical Review Letters 108, 234501. doi:10.1103/PhysRevLett.108.234501.

Glatzmaier, G.A., Coe, R.S., Hongre, L., Roberts, P.H., 1999. The role of the Earth's mantle in controlling the frequency of geomagnetic reversals. Nature 401, 885–890. doi:10.1038/44776.

Glatzmaier, G.A., Roberts, P.H., 1995. A three-dimensional self-consistent computer simulation of a geomagnetic field reversal. Nature 377, 203–209. doi:10.1038/377203a0.

Guervilly, C., 2022. Fingering Convection in the Stably Stratified Layers of Planetary Cores. Journal of Geophysical Research: Planets 127, e2022JE007350. doi:10.1029/2022JE007350.

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

Helffrich, G., Kaneshima, S., 2013. Causes and consequences of outer core stratification. Physics of the Earth and Planetary Interiors 223, 2–7. doi:10.1016/j.pepi.2013.07.005.

Jones, C.A., 2011. Planetary Magnetic Fields and Fluid Dynamos. Annual Review of Fluid Mechanics 43, 583–614. doi:10.1146/annurev-fluid-122109-160727.

Jones, C.A., Tsang, Y.K., 2025. Low inertia reversing geodynamos. Physics of the Earth and Planetary Interiors 360, 107303. doi:10.1016/j.pepi.2024.107303.

Kaneshima, S., 2018. Array analyses of SmKS waves and the stratification of Earth's outermost core. Physics of the Earth and Planetary Interiors 276, 234–246. doi:10.1016/j.pepi.2017.03.006.

Kutzner, C., Christensen, U., 2002. From stable dipolar towards reversing numerical dynamos. Physics of the Earth and Planetary Interiors 131, 29–45. doi:10.1016/S0031-9201(02)00016-X.

Landeau, M., Fournier, A., Nataf, H.C., Cébron, D., Schaeffer, N., 2022. Sustaining Earth's magnetic dynamo. Nature Reviews Earth & Environment 3, 255–269. doi:10.1038/s43017-022-00264-1.

Li, J., Sato, T., Kageyama, A., 2002. Repeated and Sudden Reversals of the Dipole Field Generated by a Spherical Dynamo Action. Science 295, 1887–1890. doi:10.1126/science.1066959.

Lowrie, W., Kent, D.V., 2013. Geomagnetic Polarity Timescales and Reversal Frequency Regimes, in: Channell, J., Kent, D., Lowrie, W., Meert, J. (Eds.), Geophysical Monograph Series. American Geophysical Union, Washington, D. C., pp. 117–129. doi:10.1029/145GM09.

Menu, M.D., Petitdemange, L., Galtier, S., 2020. Magnetic effects on fields morphologies and reversals in geodynamo simulations. Physics of the Earth and Planetary Interiors 307, 106542. doi:10.1016/j.pepi.2020.106542.

Moffatt, H.K., 1978. Magnetic Field Generation in Electrically Conducting Fluids. Cambridge Monographs on Mechanics and Applied Mathematics, Cambridge University Press, Cambridge [Eng.]; New York.

Monville, R., Vidal, J., Cébron, D., Schaeffer, N., 2019. Rotating double-diffusive convection in stably stratified planetary cores. Geophysical Journal International 219, S195–S218. doi:10.1093/gji/ggz347.

Moore, K.M., Yadav, R.K., Kulowski, L., Cao, H., Bloxham, J., Connerney, J.E.P., Kotsiaros, S., Jørgensen, J.L., Merayo, J.M.G., Stevenson, D.J., Bolton, S.J., Levin, S.M., 2018. A complex dynamo inferred from the hemispheric dichotomy of Jupiter's magnetic field. Nature 561, 76–78. doi:10.1038/s41586-018-0468-5.

Nakagawa, T., 2011. Effect of a stably stratified layer near the outer boundary in numerical simulations of a magnetohydrodynamic dynamo in a rotating spherical shell and its implications for Earth's core. Physics of the Earth and Planetary Interiors 187, 342–352. doi:10.1016/j.pepi.2011.06.001.

Nakagawa, T., 2015. An implication for the origin of stratification below the core–mantle boundary region in numerical dynamo simulations in a rotating spherical shell. Physics of the Earth and Planetary Interiors 247, 94–104. doi:10.1016/j.pepi.2015.02.007.

Olson, P., Amit, H., 2014. Magnetic reversal frequency scaling in dynamos with thermochemical convection. Physics of the Earth and Planetary Interiors 229, 122–133. doi:10.1016/j.pepi.2014.01.009.

Olson, P., Driscoll, P., Amit, H., 2009. Dipole collapse and reversal precursors in a numerical dynamo. Physics of the Earth and Planetary Interiors 173, 121–140. doi:10.1016/j.pepi.2008.11.010.

Olson, P.L., Coe, R.S., Driscoll, P.E., Glatzmaier, G.A., Roberts, P.H., 2010. Geodynamo reversal frequency and heterogeneous core–mantle boundary heat flow. Physics of the Earth and Planetary Interiors 180, 66–79. doi:10.1016/j.pepi.2010.02.010.

Olson, P.L., Glatzmaier, G.A., Coe, R.S., 2011. Complex polarity reversals in a geodynamo model. Earth and Planetary Science Letters 304, 168–179. doi:10.1016/j.eps1.2011.01.031.

Pétrélis, F., Besse, J., Valet, J.P., 2011. Plate tectonics may control geomagnetic reversal frequency. Geophysical Research Letters 38, L19303. doi:10.1029/2011GL048784.

Pétrélis, F., Fauve, S., 2008. Chaotic dynamics of the magnetic field generated by dynamo action in a turbulent flow. Journal of Physics: Condensed Matter 20, 494203. doi:10.1088/0953-8984/20/49/494203.

Pétrélis, F., Fauve, S., Dormy, E., Valet, J.P., 2009. Simple Mechanism for Reversals of Earth's Magnetic Field. Physical Review Letters 102, 144503. doi:10.1103/PhysRevLett.102.144503.

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

Proctor, M.R.E., 1977. On the eigenvalues of kinematic α-effect dynamos. Astronomische Nachrichten 298, 19–25. doi:10.1002/asna.19772980104.

Schaeffer, N., 2013. Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations. Geochemistry, Geophysics, Geosystems 14, 751–758. doi:10.1002/ggge.20071.

Sheyko, A., Finlay, C.C., Jackson, A., 2016. Magnetic reversals from planetary dynamo waves. Nature 539, 551–554. doi:10.1038/nature19842.

Su, W.J., Woodward, R.L., Dziewonski, A.M., 1994. Degree 12 model of shear velocity heterogeneity in the mantle. Journal of Geophysical Research: Solid Earth 99, 6945–6980. doi:10.1029/93JB03408.

Sumita, I., Olson, P., 2002. Rotating thermal convection experiments in a hemispherical shell with heterogeneous boundary heat flux: Implications for the Earth's core. Journal of Geophysical Research: Solid Earth 107. doi:10.1029/2001JB000548.

Takahashi, F., Matsushima, M., Honkura, Y., 2005. Simulations of a QuasiTaylor State Geomagnetic Field Including Polarity Reversals on the Earth Simulator. Science 309, 459–461. doi:10.1126/science.1111831.

Takahashi, F., Shimizu, H., Tsunakawa, H., 2019. Mercury's anomalous magnetic field caused by a symmetry-breaking self-regulating dynamo. Nature Communications 10, 208. doi:10.1038/s41467-018-08213-7.

Tassin, T., Gastine, T., Fournier, A., 2021. Geomagnetic semblance and dipolar–multipolar transition in top-heavy double-diffusive geodynamo
models. Geophysical Journal International 226, 1897–1919. doi:10.1093/gji/ggab161.

Terra-Nova, F., Amit, H., 2024. Regionally-triggered geomagnetic reversals. Scientific Reports 14, 9639. doi:10.1038/s41598-024-59849-z.

Wicht, J., Olson, P., 2004. A detailed study of the polarity reversal mechanism in a numerical dynamo model. Geochemistry, Geophysics, Geosystems 5, 2003GC000602. doi:10.1029/2003GC000602.

Yan, C., Barik, A., Stanley, S., Leung, J.S.Y., Mittelholz, A., Johnson, C.L., Plesa, A.C., Rivoldini, A., 2023. An Ancient Martian Dynamo Driven by Hemispheric Heating: Effect of Thermal Boundary Conditions. The Planetary Science Journal 4, 11. doi:10.3847/PSJ/acaec93.

Yan, C., Stanley, S., 2021. Recipe for a Saturn-Like Dynamo. AGU Advances 2, e2020AV000318. doi:10.1029/2020AV000318.

Zaire, B., Jouve, L., Gastine, T., Donati, J.F., Morin, J., Landin, N., Folsom, C.P., 2022. Transition from multipolar to dipolar dynamos in stratified systems. Monthly Notices of the Royal Astronomical Society 517, 3392–3406. doi:10.1093/mnras/stac2769.