# Necessary Core-Mantle Boundary Heat Flux Patterns for Recovering the Latitude of the South Atlantic Anomaly

Filipe Terra-nova, Hagay Amit

► To cite this version:

Filipe Terra-nova, Hagay Amit. Necessary Core-Mantle Boundary Heat Flux Patterns for Recovering the Latitude of the South Atlantic Anomaly. Journal of Geophysical Research: Solid Earth, 2025, 130 (11), pp.e2025JB031193. 10.1029/2025JB031193. hal-05368906

HAL Id: hal-05368906

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

Submitted on 17 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# JGR Solid Earth

## RESEARCH ARTICLE

10.1029/2025JB031193

**Key Points:**

* Equatorially anti-symmetric CMB heat flux is required to produce large latitudes of surface magnetic intensity minima in numerical dynamos

* Dynamo models with large latitudes of surface intensity minima are often not dipole dominated

* Larger amplitude of anti-symmetric CMB heat flux heterogeneity combined with faster rotating dynamo models may reconcile these two observations

**Correspondence to:**

F. Terra-Nova,
terranova@ipgp.fr

**Citation:**

Terra-Nova, F., & Amit, H. (2025). Necessary core-mantle boundary heat flux patterns for recovering the latitude of the South Atlantic Anomaly. *Journal of Geophysical Research: Solid Earth*, 130, e2025JB031193. https://doi.org/10.1029/2025JB031193

Received 17 JAN 2025

Accepted 3 NOV 2025

**Author Contributions:**

**Conceptualization:** Filipe Terra-Nova, Hagay Amit

**Data curation:** Filipe Terra-Nova

**Formal analysis:** Filipe Terra-Nova

**Funding acquisition:** Filipe Terra-Nova, Hagay Amit

**Investigation:** Filipe Terra-Nova

**Methodology:** Filipe Terra-Nova, Hagay Amit

**Project administration:** Hagay Amit

**Software:** Filipe Terra-Nova

**Validation:** Filipe Terra-Nova, Hagay Amit

**Visualization:** Filipe Terra-Nova

**Writing – original draft:** Filipe Terra-Nova, Hagay Amit

**Writing – review & editing:** Filipe Terra-Nova, Hagay Amit

© 2025. The Author(s).
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

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

## Necessary Core-Mantle Boundary Heat Flux Patterns for Recovering the Latitude of the South Atlantic Anomaly

Filipe Terra-Nova¹,² Ⓓ and Hagay Amit¹ Ⓓ

¹Nantes Université, University Angers, Le Mans Université, CNRS, Laboratoire de Planétologie et Géosciences, LPG UMR 6112, Nantes, France, ²Université Paris Cité, Institut de physique du globe de Paris, CNRS, F-75005, Paris, France

**Abstract** The present-day geomagnetic field at Earth's surface is characterized by an anomalously weak region known as the South Atlantic Anomaly (SAA). Numerical dynamo simulations with an imposed outer boundary heat flux pattern inferred from lowermost mantle seismic anomalies reproduce well the longitude of the SAA, but the latitudes of the modeled surface intensity minima are significantly lower than the observed (Terra-Nova et al., 2019, https://doi.org/10.1093/gji/ggy519). Here we explore numerical dynamos with a variety of large-scale synthetic patterns of outer boundary heat flux in order to gain fundamental insight into what patterns may yield large latitudes of surface intensity minima. We find that equatorially anti-symmetric heat flux patterns are required in order to produce large latitudes of surface intensity minima. However, these dynamo models often reverse and their dipole dominance is lost. This trade-off between dipole-dominated models with low latitudes of surface intensity minima as opposed to multipolar models with large latitudes of surface intensity minima is also evident upon increase in the strength of the convection vigor. Snapshots of dynamo models that satisfy the geomagnetic observations in terms of both large latitudes of surface intensity minima as well as large dipolarity exist, but are rare. We speculate that reconciling these two fundamental geomagnetic constraints may be achieved with a larger amplitude of an equatorially anti-symmetric outer boundary heat flux combined with a more rapidly rotating dynamo model.

**Plain Language Summary** The present-day geomagnetic field at Earth's surface is characterized by an anomalously weak region known as the South Atlantic Anomaly (SAA). In numerical models of core convection such features appear very close to the equator, in contrast to the center of the SAA that resides nowadays in Patagonia. Here we explore numerical dynamo simulations with a variety of large-scale synthetic patterns of core-mantle boundary (CMB) heat flux in order to understand what patterns may yield large latitudes of surface intensity minima. We find that equatorially anti-symmetric heat flux patterns result in large latitudes of surface intensity minima. However, in these models the field morphology is very different from that of the geomagnetic field. We speculate that reconciling these two fundamental geomagnetic constraints, field morphology on the CMB and large latitude of surface intensity minima, may be achieved with a larger amplitude of an equatorially anti-symmetric outer boundary heat flux heterogeneity combined with a more rapidly rotating dynamo model.

## 1. Introduction

Some of the most outstanding features of the geomagnetic field are surface intensity minima (Domingos et al., 2017; He et al., 2021; Tarduno et al., 2015; Terra-Nova et al., 2017). Monitoring geomagnetic surface intensity minima and unraveling their dynamical origin is crucial for understanding the dynamics at Earth's liquid outer core and the interactions with its surrounding solid envelopes. Studying surface intensity minima also has a major societal impact due to the effects of the weak field on the Van Allen radiation belts where more energetic particles from space penetrate Earth's atmosphere (Heirtzler, 2002; Olson & Amit, 2006). Satellites and other spacecrafts passing through regions of low field intensity are bombarded by solar energetic particles that degrade electronics, optics, solar panels and other critical systems by breaking chemical bonds and disrupting crystalline structure (e.g. Lean, 2005).

Since the advent of intensity measurements, global field models constructed from direct observations showed that reversed geomagnetic flux patches on the core-mantle boundary (CMB) have been expanding and intensifying together with the decreasing axial dipole (e.g., Gubbins, 1987; Olson & Amit, 2006). The historical dipole intensity decrease has been linked to the time-dependent area of reversed flux on the CMB (Metman et al., 2018), to
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a, b): Radial geomagnetic field at the CMB based on the CHAOS7.10 (Finlay et al., 2020) model expanded until maximum spherical harmonic degree and order $ℓ_{max} = 14$. (c, d): Surface intensity with black and red contours corresponding to $1.2F_{min}$ and $1.08F_{min}$, respectively, where $F_{min}$ is the minimum surface intensity. The 1.2 and 1.08 factors were chosen to highlight the emergence of the second minimum in 2020. Black dots denote local minima. (e, f): Zoom-ins to the region of the surface intensity minimum or minima. Note the different scales between (e, f) versus (c, d).</figcaption></figure>

the absence of a poleward drifting intense high-latitude normal polarity flux patch below Patagonia (Finlay et al., 2016) as well as to a complex influence of small-scale magnetic field features in the southern hemisphere (Finlay et al., 2012). At Earth's surface, a region of anomalously weak field intensity termed the South Atlantic Anomaly (SAA) is observed (Heirtzler, 2002). At high altitudes the SAA is approximated by an eccentric dipole (Domingos et al., 2017). Detailed analysis of global geomagnetic field models shows that the SAA is indeed related to prominent flux patches on the CMB (Bloxham et al., 1989) albeit in a non-trivial way (Terra-Nova et al., 2017).

Although the satellite era is too short to accommodate significant changes in the geomagnetic field, within a period of 20 years a transition from one surface intensity minimum to two minima is observed (Figures 1c–1f). The radial magnetic field at the CMB for a model CHAOS-7.10 of the satellite era (Finlay et al., 2020) in 2000 and 2020 is shown in Figures 1a and 1b. Intense normal flux patches (NFPs) are found at high-latitudes close to the intersection of the inner core Tangent Cylinder (TC) with the CMB as well as close to the equator (e.g., Jackson, 2003). Prominent reversed flux patches (RFPs) are found below the South Atlantic region (e.g., Bloxham et al., 1989). From 2000 to 2020, the normal flux below mid latitudes of the Indian ocean became patchier. At the same time the RFP below South Africa moved westward and the RFP below in between South Africa and Antarctica intensified. These changes in the radial geomagnetic field on the CMB between 2000 and 2020 culminated in corresponding surface intensity changes for the same years (Figures 1c–1f) from one minimum in 2000 to two local minima in 2020, with the stronger anomaly (i.e., weaker geomagnetic field intensity) below South America and another less pronounced minimum at a higher latitude close to Africa (see Alken et al., 2021; Finlay et al., 2020; Rother et al., 2021; Terra-Nova et al., 2019).
A field forecasting study by data assimilation schemes predicts that the African minimum will become weaker than the South American (Sanchez et al., 2020). Importantly, the African minimum is centered at a higher latitude than the South American. According to refined measures that better characterize the SAA area and center by accounting for surface intensity outside the SAA region and shape anisotropy, the evolution of the SAA area and center during the historical era has not been monotonous, including episodes of steady area, eastward drift and rapid southward drift (Amit et al., 2021). Much attention has been given to the absolute minimum of the SAA in past studies, whilst considerably less to other identified minima (Terra-Nova et al., 2019). Nevertheless, the low intensity West Pacific Anomaly found in historical and archeological global geomagnetic field models (Tema et al., 2023; Yue et al., 2024) has been linked to an independent database of equatorial aurorae recorded in Korea (He et al., 2021). These studies highlight the fact that surface intensity minima may emerge at multiple locations, providing a new constraint on models of core dynamics and suggesting that the societal impact of surface intensity minima might not be restricted to the South Atlantic in the future.

Based on analysis of archeomagnetic field models, it has been proposed that the location of the SAA is recurrent over millennial timescales (Brown et al., 2018; Campuzano et al., 2019; Hartmann et al., 2019; Hare et al., 2018; Hellio & Gillet, 2018; Panovska et al., 2019; Trindade et al., 2018). The persistence of the SAA is possibly related to the impact of the lower mantle heterogeneity on the geodynamo. Thermal core-mantle interactions affect the convection pattern in the outer core and thus the morphology of the generated geomagnetic field. Correlations between longitudes of lowermost mantle seismic anomalies and prominent geomagnetic flux patches (normal and reversed) provide evidence for mantle control on the geodynamo (Gubbins, 2003; Terra-Nova et al., 2016). Archeo- and paleomagnetic studies based on materials from three specific regions, South America (Engbers et al., 2022; Hartmann et al., 2010, 2011, 2019; Poletti et al., 2016; Shah et al., 2016; Trindade et al., 2018), Central Africa (Hare et al., 2018; Tarduno et al., 2015) and West Pacific (Cai et al., 2021; Li et al., 2023; Luo et al., 2025; Tema et al., 2023), registered anomalous weak surface field that may be associated to SAA-like features passing through the measurements sites. Alternatively, the current SAA location could be chaotically variable. Based on a data assimilation scheme, the SAA will drift in the near future to the Pacific (Aubert, 2015), possibly indicating a transient geodynamo feature. In contrast, the forecast of Sanchez et al. (2020) suggests that the African minimum will overcome the South American in the near future, that is, the global minimum will remain below the South Atlantic region. These different forecasts likely arise due to the different dynamo models that they rely on. In the coupled Earth setup (Aubert et al., 2013) used by Aubert (2015) the gravitational coupling between the mantle and inner core leads to westward drift, whereas in Sanchez et al. (2020) the prior information is based on a dynamo model without such gravitational coupling.

Terra-Nova et al. (2019) designed a quantitative algorithm to identify objectively the locations of surface intensity minima. They applied their algorithm to a set of geomagnetic field models as well as to numerical dynamo simulations with heterogeneous outer boundary heat flux inferred from lower mantle seismic anomalies. Their results demonstrated that the SAA longitude is consistent with mantle control on the geodynamo. However, the latitude of the SAA at present, around 26.5°, is significantly larger than the <6° latitude of surface intensity minima in some archeomagnetic field models (Korte et al., 2009; Korte & Constable, 2011; Licht et al., 2013; Nilsson et al., 2014) as well as in dynamo models with a tomographic outer boundary heat flux pattern (Terra-Nova et al., 2019). It therefore remains a challenge to explain the large latitude of the SAA.

Most studies of numerical dynamos that incorporate CMB thermal heterogeneity impose an outer boundary heat flux pattern proportional to lowermost mantle seismic anomalies (e.g., Aubert et al., 2008; Gubbins et al., 2007; Mound & Davies, 2023; Olson & Christensen, 2002; Sahoo & Sreenivasan, 2020). However, non-thermal contributions to the seismic heterogeneity (Choblet et al., 2023; Nakagawa & Tackley, 2008) may affect the recovery of persistent geodynamo features (Amit & Choblet, 2009, 2012; Amit, Deschamps & Choblet, 2015). Identifying the required CMB heat flux ingredient for recovering the SAA latitude may constrain CMB heat flux models.

In this study we aim to understand what type of CMB heat flux pattern may lead to instantaneous large latitude of surface intensity minima. For this purpose, we explore numerical dynamo simulations with synthetic single harmonic outer boundary heat flux patterns, an approach that has been often used to study various geodynamo features (Amit et al., 2015), for example field morphology (Aubert et al., 2007; Olson & Christensen, 2002), field intensity (Takahashi et al., 2008) and paleomagnetic reversal frequency (Frasson et al., 2025; Glatzmaier et al., 1999; Olson et al., 2010). We set the target value for the latitude of the surface intensity anomaly based on geomagnetic field models for the satellite (e.g., Figures 1c–1f) as well as the observatory eras.
In parallel to the SAA latitude, we also require that the dynamo models should exhibit a fundamental Earth-like morphology, that is, its high level of dipolarity (Christensen et al., 2010). This property is commonly quantified by the ratio between the power of the dipole and the power of the total field $f_{dip}$. Based on numerical dynamo simulations, the lower bound of $f_{dip}$ values for dipole-dominated models has been proposed to be between 0.35 – 0.50 (e.g., Christensen & Aubert, 2006; Frasson et al., 2025; Oruba & Dormy, 2014; Tassin et al., 2021). Geomagnetic field models may set an upper bound for Earth-like $f_{dip}$. These observation based models are characterized by time-average dipolarities of $\langle f_{dip} \rangle = 0.714 \pm 0.021$ and $0.684 \pm 0.002$ for the observatory (Huder et al., 2020) and satellite (Finlay et al., 2020) eras (where $\langle \rangle$ denotes time averaging and $\pm$ denotes the standard deviation which represents time dependence), respectively. Accordingly, various ranges for Earth-like $f_{dip}$ values have been proposed, for example, 0.35 – 0.70 (Aubert et al., 2009), 0.45 – 0.75 (Christensen, 2010) or 0.35 – 0.75 (Davies et al., 2022). Because we focus on the SAA latitude at present day, we select a large lower bound of $f_{dip} = 0.5$ as a requirement for an Earth-like modeled field morphology.

The paper is outlined as follows. In Section 2.1 we describe the numerical dynamo models including the imposed outer boundary heat flux patterns. In Section 2.2 we recall the algorithm of Terra-Nova et al. (2019) for the determination of the latitude of the surface intensity minima. The results are presented in Section 3 and geophysical implications are discussed in Section 4.

## 2. Methods

### 2.1. Numerical Dynamos

Our numerical dynamo simulations have the same generic setup of thermo-chemical convection as in Terra-Nova et al. (2019). Here we recall the governing equations, the control parameters and the boundary conditions. We then introduce the single harmonic patterns used as imposed outer boundary heat flux in this study.

The numerical dynamo models solve the set of magnetohydrodynamic (MHD) equations, which correspond to conservation of momentum, electromagnetic induction, conservation of heat (or co-density), incompressibility and non-existent magnetic monopoles. The numerical technique relies on a pseudo-spectral approach in which non-linear terms are solved on a grid whereas diffusive terms are treated implicitly using spherical harmonic expansions (e.g., Olson et al., 1999). These equations can be expressed in non-dimensional form (e.g., Wicht, 2002) as follows:

$$E\left(\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} - \nabla^2 \mathbf{u}\right) + 2\hat{z} \times \mathbf{u} + \nabla P = Ra^* \frac{\mathbf{r}}{r_o} C + \frac{1}{Pm} (\nabla \times \mathbf{B}) \times \mathbf{B}, \quad (1)$$

$$\frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{u} \times \mathbf{B}) + \frac{1}{Pm} \nabla^2 \mathbf{B}, \quad (2)$$

$$\frac{\partial C}{\partial t} + \mathbf{u} \cdot \nabla C = \frac{1}{Pr} \nabla^2 C + \epsilon, \quad (3)$$

$$\nabla \cdot \mathbf{u} = 0, \quad (4)$$

$$\nabla \cdot \mathbf{B} = 0, \quad (5)$$

where $\mathbf{u}$ is the fluid velocity, $\hat{z}$ the direction of the axis of rotation, $P$ the modified pressure that includes the centrifugal force, $\mathbf{r}$ the radial position vector, $r_o$ the core radius, $C$ the co-density, $\mathbf{B}$ the magnetic field and $\epsilon$ the buoyancy source or sink. The co-density is given by $C = \alpha T + \beta \xi$ where $T$ is temperature, $\xi$ is light elements concentration, and $\alpha$ and $\beta$ are their respective expansivities. Equations 1–3 are governed by four (internal) control parameters. The Ekman number represents the ratio of viscous to Coriolis forces:

$$E = \frac{\nu}{\Omega D^2}. \quad (6)$$
The heat flux based Rayleigh number represents the convection vigor versus retarding forces:

$$Ra = \frac{\alpha g_0 q_0 D^4}{\kappa \nu k} \qquad (7)$$

The Prandtl number and the magnetic Prandtl number are ratios of diffusivities:

$$Pr = \frac{\nu}{\kappa}, \qquad (8)$$

$$Pm = \frac{\nu}{\eta}. \qquad (9)$$

In Equations 6–9 Ω is the rotation rate, ν the kinematic viscosity, D the shell thickness, g₀ the gravitational acceleration at the outer boundary, q₀ the mean outer boundary heat flux, k the thermal conductivity, κ the thermal diffusivity and η the magnetic diffusivity. The modified Rayleigh number Ra* in Equation 1 is related to the classical Rayleigh number Ra (Equation 7) by $Ra* = RaE/Pr$.

The choice of values for the control parameters was guided by the results of Terra-Nova et al. (2019). They found that the latitude of the surface intensity minima is augmented by more vigorous convection, whereas faster rotation favors lower latitudes, while the dependence on Pm is secondary (see their Figure 18). Therefore, we chose combinations of control parameters that yield dynamos in a parametric regime not far from the limit between the non-reversing and reversing regimes, in order to allow for a significant sensitivity of the latitude of surface intensity minima to the convective conditions (Terra-Nova et al., 2019). For low Ekman numbers this regime becomes computationally costly, especially if long runs are required for meaningful statistical analysis of boundary heterogeneity effects. We therefore chose a moderate E value of 1 × 10⁻⁴ for most simulations. Additional runs with a lower E value of 3 × 10⁻⁵ were examined to test the possible effect of lower E. The Ra values were set to yield dipole-dominated non-reversing models for a reference homogeneous outer boundary heat flux case (e.g., Christensen & Aubert, 2006; Kutzner & Christensen, 2002) which resembles Earth's field morphology during a chron. The Pm value is relatively low in order to obtain long runs. For all simulations we used Pr = 1 and Pm = 2 with the exception of two models with Pm = 4 and 6 to explore the effects of changing Pm. We mostly explored the dependence of the results on Ra and the outer boundary heat flux amplitude q* and patterns (Table 1).

All dynamo simulations were carried out with rigid and electrically insulating conditions at both boundaries. The spherical shell has an inner to outer core radii ratio of 0.35, corresponding to the present-day Earth's core geometry. The models differ by their imposed outer boundary heat flux. Different heterogeneous single harmonic patterns, or combinations of two single harmonics, were tested. Because our focus is the latitude of surface intensity minima, most of these synthetic patterns are axial (Figure 2), with the exception of $Y_2$ which is the dominant large-scale pattern in seismic tomography models of the lowermost mantle (Masters et al., 2000). The amplitude of the imposed heat flux heterogeneity is quantified by Olson and Christensen (2002):

$$q^* = \frac{q_{\max} - q_{\min}}{2q_0}, \qquad (10)$$

where $q_{\max}$ and $q_{\min}$ are the maximum and minimum heat flux, respectively. We chose moderate amplitudes of outer boundary heterogeneity and we explored the dependence of the results on q*. For the inner boundary fixed co-density is imposed and no sink or source term was assigned at the volume (i.e., $\epsilon = 0$ in Equation 3), corresponding to generic thermo-chemical convection (e.g., Aubert et al., 2008; Olson et al., 2017; Terra-Nova et al., 2019). All dynamo simulations were run for more than 300 advection times, where the advection time is $\tau_u = U/D$ with U being a typical velocity scale. Following Terra-Nova et al. (2019), for each dynamo model we analyzed a long timeseries of radial magnetic field snapshots on the outer boundary corresponding to run times that range ≈ 44–262 thousand years. Typical time between analyzed snapshots ranges ≈ 0.04–0.23 advection times, which corresponds to 6–35 years in Earth's core, considering $\tau_u = 140$ years (e.g., Terra-Nova & Amit, 2020). Note that this value does not correspond to the time stepping of the simulations which is 10⁻³ smaller. This dense temporal sampling allows for practically a continuous tracking of the evolution of the radial magnetic field on the outer boundary and the intensity at Earth's surface (Terra-Nova et al., 2019). All simulations
<table><caption>**Table 1**

Numerical Dynamo Models Setup and Main Results</caption><thead><tr><th>Case</th><th>Pattern</th><th>Ra</th><th>E</th><th>Pm</th><th>q*</th><th>⟨Rm⟩</th><th>⟨f<sub>dip</sub>⟩</th><th>⟨M⟩</th><th>⟨P<sub>λ</sub>⟩</th><th>σ</th><th>⟨P<sub>λ</sub><sup>g</sup>⟩</th><th>⟨N<sub>Fmin</sub>⟩</th><th>Rev?</th><th>Duration</th></tr></thead><tbody><tr><td>U</td><td>Y<sub>0</sub><sup>0</sup></td><td>4 × 10<sup>7</sup></td><td>1 × 10<sup>-4</sup></td><td>2</td><td>1.0</td><td>393</td><td>0.72</td><td>1.42</td><td>0.05</td><td>3.09</td><td>-0.16</td><td>2.77</td><td>No</td><td>1,310</td></tr><tr><td>H10a</td><td>-Y<sub>1</sub><sup>0</sup></td><td>2 × 10<sup>7</sup></td><td>1 × 10<sup>-4</sup></td><td>2</td><td>2.0</td><td>422</td><td>0.20</td><td>0.14</td><td>21.67</td><td>26.83</td><td>21.74</td><td>2.40</td><td>Yes</td><td>1,705</td></tr><tr><td>H10b</td><td>-Y<sub>1</sub><sup>0</sup></td><td>4 × 10<sup>7</sup></td><td>1 × 10<sup>-4</sup></td><td>2</td><td>1.0</td><td>413</td><td>0.51</td><td>1.97</td><td>7.30</td><td>10.78</td><td>9.60</td><td>2.67</td><td>No</td><td>763</td></tr><tr><td>H1020</td><td>-Y<sub>1</sub><sup>0</sup> + Y<sub>2</sub><sup>0</sup></td><td>4 × 10<sup>7</sup></td><td>1 × 10<sup>-4</sup></td><td>2</td><td>1.0</td><td>372</td><td>0.54</td><td>2.78</td><td>0.60</td><td>4.88</td><td>0.41</td><td>2.53</td><td>No</td><td>503</td></tr><tr><td>H20a</td><td>Y<sub>2</sub><sup>0</sup></td><td>4 × 10<sup>7</sup></td><td>1 × 10<sup>-4</sup></td><td>2</td><td>1.0</td><td>358</td><td>0.50</td><td>3.02</td><td>0.00</td><td>5.73</td><td>0.13</td><td>2.47</td><td>No</td><td>472</td></tr><tr><td>H20b</td><td>-Y<sub>2</sub><sup>0</sup></td><td>4 × 10<sup>7</sup></td><td>1 × 10<sup>-4</sup></td><td>2</td><td>0.7</td><td>527</td><td>0.17</td><td>0.24</td><td>8.13</td><td>33.66</td><td>9.75</td><td>4.16</td><td>Yes</td><td>1,870</td></tr><tr><td>H22</td><td>Y<sub>2</sub><sup>2</sup></td><td>4 × 10<sup>7</sup></td><td>1 × 10<sup>-4</sup></td><td>2</td><td>1.0</td><td>423</td><td>0.73</td><td>1.13</td><td>0.03</td><td>4.02</td><td>0.03</td><td>2.51</td><td>No</td><td>1,470</td></tr><tr><td>H30a</td><td>-Y<sub>3</sub><sup>0</sup></td><td>2 × 10<sup>7</sup></td><td>1 × 10<sup>-4</sup></td><td>2</td><td>1.0</td><td>277</td><td>0.68</td><td>2.38</td><td>-0.82</td><td>5.35</td><td>-0.47</td><td>2.68</td><td>No</td><td>1,064</td></tr><tr><td>H30b</td><td>-Y<sub>3</sub><sup>0</sup></td><td>2 × 10<sup>7</sup></td><td>1 × 10<sup>-4</sup></td><td>2</td><td>2.0</td><td>338</td><td>0.27</td><td>0.46</td><td>15.43</td><td>34.45</td><td>17.58</td><td>3.20</td><td>Yes</td><td>1,441</td></tr><tr><td>H30c</td><td>-Y<sub>3</sub><sup>0</sup></td><td>4 × 10<sup>7</sup></td><td>1 × 10<sup>-4</sup></td><td>2</td><td>1.0</td><td>489</td><td>0.19</td><td>0.23</td><td>5.06</td><td>35.89</td><td>3.37</td><td>3.96</td><td>Yes</td><td>1,677</td></tr><tr><td>H30d</td><td>-Y<sub>3</sub><sup>0</sup></td><td>4 × 10<sup>7</sup></td><td>1 × 10<sup>-4</sup></td><td>2</td><td>2.0</td><td>492</td><td>0.23</td><td>0.26</td><td>11.52</td><td>38.11</td><td>14.84</td><td>3.68</td><td>Yes</td><td>1,648</td></tr><tr><td>H30e</td><td>-Y<sub>3</sub><sup>0</sup></td><td>2 × 10<sup>8</sup></td><td>3 × 10<sup>-5</sup></td><td>2</td><td>2.0</td><td>545</td><td>0.62</td><td>4.22</td><td>0.06</td><td>10.76</td><td>0.18</td><td>3.01</td><td>No</td><td>312</td></tr><tr><td>H30f</td><td>-Y<sub>3</sub><sup>0</sup></td><td>2 × 10<sup>8</sup></td><td>3 × 10<sup>-5</sup></td><td>4</td><td>2.0</td><td>1115</td><td>0.42</td><td>3.50</td><td>9.48</td><td>23.30</td><td>14.41</td><td>3.40</td><td>No</td><td>400</td></tr><tr><td>H30g</td><td>-Y<sub>3</sub><sup>0</sup></td><td>2 × 10<sup>8</sup></td><td>3 × 10<sup>-5</sup></td><td>6</td><td>2.0</td><td>1668</td><td>0.34</td><td>3.86</td><td>12.32</td><td>25.51</td><td>18.92</td><td>3.60</td><td>No</td><td>371</td></tr><tr><td>H3020a</td><td>-Y<sub>3</sub><sup>0</sup> + Y<sub>2</sub><sup>0</sup></td><td>4 × 10<sup>7</sup></td><td>1 × 10<sup>-4</sup></td><td>2</td><td>1.0</td><td>374</td><td>0.60</td><td>2.03</td><td>-1.02</td><td>4.74</td><td>-1.26</td><td>2.70</td><td>No</td><td>773</td></tr><tr><td>H3020b</td><td>-Y<sub>3</sub><sup>0</sup> + Y<sub>2</sub><sup>0</sup></td><td>4 × 10<sup>8</sup></td><td>3 × 10<sup>-5</sup></td><td>2</td><td>2.0</td><td>922</td><td>0.30</td><td>0.35</td><td>13.54</td><td>31.86</td><td>20.03</td><td>4.11</td><td>Yes</td><td>350</td></tr><tr><td>H3022a</td><td>-Y<sub>3</sub><sup>0</sup> + Y<sub>2</sub><sup>2</sup></td><td>2 × 10<sup>7</sup></td><td>1 × 10<sup>-4</sup></td><td>2</td><td>1.0</td><td>284</td><td>0.70</td><td>1.97</td><td>-0.96</td><td>4.52</td><td>-0.97</td><td>2.54</td><td>No</td><td>1,312</td></tr><tr><td>H3022b</td><td>-Y<sub>3</sub><sup>0</sup> + Y<sub>2</sub><sup>2</sup></td><td>4 × 10<sup>7</sup></td><td>1 × 10<sup>-4</sup></td><td>2</td><td>1.0</td><td>494</td><td>0.18</td><td>0.21</td><td>3.52</td><td>32.96</td><td>2.01</td><td>4.14</td><td>Yes</td><td>1,709</td></tr></tbody></table>

Note. U and H denote uniform and heterogeneous outer boundary heat flux, respectively, case numbers denote the imposed single harmonic heat flux pattern degree and order. Control parameters *E*, *Ra*, *Pr*, *Pm* and *q*<sup>*</sup> are defined in Equations 6–10. In all models *Pr* = 1. The number of radial grid points in the shell *N*<sub>r</sub> is 65 and 81 for models with *E* = 1 × 10<sup>-4</sup> and *E* = 3 × 10<sup>-5</sup>, respectively. The maximum spherical harmonic degree and order is *ℓ*<sub>max</sub> = 133 for models with *Pm* = 2 and *ℓ*<sub>max</sub> = 170 for models with *Pm* = 4 and 6. ⟨⟩ denotes time averaging. The magnetic Reynolds number *Rm* is calculated from the time-average kinetic energy in the shell. The outer boundary dipolarity *f*<sub>dip</sub> is defined in Equation 12. ⟨*P*<sub>λ</sub>⟩ is the weighted latitude of local surface intensity minima (Equation 18) and *σ* its standard deviation. ⟨*P*<sub>λ</sub><sup>g</sup>⟩ is the weighted latitude of global surface intensity minimum. *N*<sub>F</sub> is the mean number of local surface intensity minima per snapshot. Rev denotes if the models reversed or not. The duration of the simulations is given in terms of the number of advection times *τ*<sub>u</sub>.

were run using the pseudo spectral MagIC code (Wicht, 2002, https://github.com/magic-sph/magic) for fluid dynamics in a spherical shell coupled with efficient spherical harmonic transforms (Schaeffer, 2013).

For each dynamo model we report several output parameters (Table 1). Three main conventional output parameters were monitored. First, the magnetic Reynolds number

$$Rm = Pm\sqrt{2E_{kin}} \qquad (11)$$

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Non-dimensional normalized outer boundary heat flux anomalies versus latitude.</figcaption></figure>
is based on the time-average kinetic energy in the shell $E_{kin}$. Second, the relative dipole field strength on the outer boundary (Christensen & Aubert, 2006) is given by

$$f_{dip} = \frac{2}{\sqrt{3}} \left( \frac{((g_1^0)^2 + (g_1^1)^2 + (h_1^1)^2)}{\sum_{\ell=1}^{12} \left( \left( \frac{r_a}{r_0} \right)^{2\ell-2} \frac{(\ell+1)^2}{(2\ell+1)} \sum_{m=0}^{\ell} \left( (g_\ell^m)^2 + (h_\ell^m)^2 \right) \right)} \right)^{\frac{1}{2}} \quad (12)$$

where $\ell$ and $m$ are spherical harmonic degree and order, respectively, $g_\ell^m$ and $h_\ell^m$ are the field's internal Gauss coefficients and $r_a$ the mean radius of Earth's surface. Third, the ratio of magnetic to kinetic energies inside the shell $M$ is given by

$$M = \frac{E_{mag}}{E_{kin}} \quad (13)$$

where the non-dimensional magnetic and kinetic energies in the spherical shell are calculated by

$$E_{mag} = \frac{1}{2} \langle \mathbf{B} \cdot \mathbf{B} \rangle \quad (14)$$

and

$$E_{kin} = \frac{1}{2} \langle \mathbf{u} \cdot \mathbf{u} \rangle. \quad (15)$$

In addition we give in Table 1 the persistent latitudes of surface intensity minima which are described below.

## 2.2. Determination of the Latitude of Surface Intensity Minima

Following Terra-Nova et al. (2019), at each time step the magnetic field snapshots were calculated using the Gauss coefficients $g_\ell^m$ and $h_\ell^m$ truncated at spherical harmonic degree and order 14, corresponding to the maximal resolution of the geomagnetic core field. Smaller scales of the geomagnetic core field are not resolved because it is non-trivial to remove crustal contributions (e.g., Merrill et al., 1998). For archeomagnetic field models the data scarcity leads to truncations at lower degrees (e.g., Wardinski et al., 2025). Statistics of intensity minima latitudes at Earth's surface were calculated using the following algorithm:

1. For every snapshot, we identified the local surface intensity minima $F_{min}$ (e.g., see Figures 3, 7 and 12 of Terra-Nova et al., 2019) and extracted their latitudes $\lambda_{Fmin}$.

2. Each $\lambda_{Fmin}$ was then interpolated to an off-grid location by a second degree polynomial fit that includes the original grid point and its neighboring North and South grid points (Terra-Nova et al., 2017).

3. In order to focus on large-scale features, if in a snapshot a pair of $F_{min}$ are separated by a great circle distance smaller than 10° (which occurs seldom in the dynamo models), these two minima were linearly interpolated to a single minimum.

4. For each identified $F_{min}$, a weight was assigned based on its relative intensity

$$W_{F_{min}} = \frac{\bar{F} - F_{min}}{\bar{F}}, \quad (16)$$

where $\bar{F}$ denotes the surface intensity averaged over Earth's surface.

5. For each dynamo model, we constructed histograms representing the distributions of surface intensity minima by summing their latitudes $\lambda_{Fmin}$ in space and time over latitude bins of 2°.

6. For a given snapshot, a weighted average was obtained using the weight defined in Equation 16:
$$P_{\lambda} = \frac{\sum_{i} (W_{F_{\min}} \lambda_{F_{\min}})_{i}}{\sum_{i} (W_{F_{\min}})_{i}}, \qquad (17)$$

where i denotes summation over multiple surface intensity minima at the considered snapshot.

7. The time-average latitude was obtained by a similar weighted average, this time summing all snapshots (i.e., averaging in time):

$$\langle P_{\lambda} \rangle = \frac{\sum_{j} (W_{F_{\min}} \lambda_{F_{\min}})_{j}}{\sum_{j} (W_{F_{\min}})_{j}}, \qquad (18)$$

where j denotes summation over multiple surface intensity minima at all snapshots and ⟨⟩ denotes time averaging. In addition, we calculated the standard deviation σ of the latitude of surface intensity minima.

The main difference between the above algorithm and that of Terra-Nova et al. (2019) concerns the off-grid procedure of surface intensity minima (Terra-Nova et al., 2017) and the implementation of Equation 16. We found in some rare snapshots negative values of $W_{F_{\min}}$ corresponding to local surface intensity minima with rather large (above average) values. Because these local minima do not represent weak field, and because the weight defined in Equation 16 must be positive in order to be implemented sensibly in Equations 17 and 18, we excluded these minima from the analysis. As in Terra-Nova et al. (2019), we also report the average number of surface intensity minima ⟨$N_{F_{\min}}$⟩ found at a given snapshot in our dynamo models (Table 1).

Some of our numerical dynamo simulations have bimodal distributions of latitude of surface intensity minima. When this was the case, we implemented a Gaussian mixture model (McLachlan & Peel, 2000) to best fit a sum of two Gaussians. The output includes the latitudes of the two Gaussians ⟨$P_{\lambda_1}$⟩ and ⟨$P_{\lambda_2}$⟩, their standard deviations $\sigma_1$ and $\sigma_2$ and their relative weighted frequency amplitudes obtained by the ratio of peak heights of the histograms $W(\langle P_{\lambda_2} \rangle)/W(\langle P_{\lambda_1} \rangle)$.

In addition, we also calculated the latitude of global surface intensity minimum $\lambda_{F_{\min}}^g$ and its weighted time average

$$\langle P_{\lambda}^{g} \rangle = \frac{\sum_{k} (W_{F_{\min}} \lambda_{F_{\min}}^{g})_{k}}{\sum_{k} (W_{F_{\min}})_{k}}, \qquad (19)$$

where k denotes summation over all snapshots. Comparing this quantity with ⟨$P_{\lambda}$⟩ allows to determine whether low latitude surface intensity minima are consequences of cancellations of several minima at either opposite hemispheres or in proximity to the equator or simply due to low latitude locations of all minima.

## 2.3. Earth-Like Latitude of Surface Intensity Minima and Dipolarity

For a snapshot of a dynamo model, Earth-like latitude of surface intensity minima is considered when either $P_{\lambda}$ or $\lambda_{F_{\min}}^g$ exceed the corresponding time-average values based on the geomagnetic field model over either the satellite (Finlay et al., 2020) or the observatory (Huder et al., 2020) era. In addition, for a snapshot of a dynamo model, Earth-like dipolarity is considered when its instantaneous $f_{dip}$ value exceeds 0.5 (see Section 1).

Here we consider a dynamo model snapshot to be Earth-like when it instantaneously satisfies the two conditions mentioned above. We term $\mathcal{T}(\langle P_{\lambda}, f_{dip} \rangle)$ as the fraction of time a dynamo model satisfies the two criteria based on $P_{\lambda}$. Similarly, $\mathcal{T}(\langle \lambda_{F_{\min}}^g, f_{dip} \rangle)$ is the corresponding time fraction based on $\lambda_{F_{\min}}^g$. The two time fractions were evaluated based on the time averages of the satellite and observatory geomagnetic field models (Finlay et al., 2020; Huder et al., 2020).

## 2.4. Magnetic Flux Hemisphericities

Deviations of surface intensity minima from equatorial locations can be explained by magnetic flux hemispherical asymmetries at the CMB. Following Terra-Nova et al. (2019), we examined the hemispherical asymmetries of
reversed and normal magnetic flux. We quantified the contributions to the axial dipole moment of reversed and normal flux of the radial field at the CMB. The axial component of the magnetic dipole moment $m_z$ can be defined as follows (e.g., Moffatt, 1978):

$$m_z = \frac{3r_o}{2\mu_0} \int_S B_r \cos \theta dS \qquad (20)$$

where $\mu_0$ is the permeability of free space, S denotes the CMB surface and the spherical surface increment is $dS = r_o^2 \sin \theta d\phi d\theta$ where $\theta$ is co-latitude and $\phi$ is longitude. Because some of our dynamo models reverse, the definition of normal and reversed flux depends on the time-dependent sign of the axial dipole. To account for this polarity, we quantified for all snapshots of all models the integrated contributions to the axial dipole of each magnetic flux type, normal or reversed (denoted by N and R subscripts), at each hemisphere, northern or southern (denoted by N and S superscripts). For example, $H_N^N$ corresponds to normal flux (i.e., $B_r \cos \theta < 0$ for present-day polarity) at the northern hemisphere (i.e., $\theta < \frac{\pi}{2}$). All magnetic flux hemisphericities for the two possible global polarities are given by

$$\text{if } g_1^0 < 0: \left\{ \begin{array}{l} H_N^N \\ H_R^N \\ H_N^S \\ H_R^S \end{array} \right. = \frac{3r_o}{2\mu_0} \int_S B_r \cos \theta dS \left\{ \begin{array}{l} \text{for } \theta < \frac{\pi}{2} \text{ and } B_r \cos \theta < 0, \\ \text{for } \theta < \frac{\pi}{2} \text{ and } B_r \cos \theta > 0, \\ \text{for } \theta > \frac{\pi}{2} \text{ and } B_r \cos \theta < 0, \\ \text{for } \theta > \frac{\pi}{2} \text{ and } B_r \cos \theta > 0 \end{array} \right. \qquad (21)$$

or

$$\text{if } g_1^0 > 0: \left\{ \begin{array}{l} H_N^N \\ H_R^N \\ H_N^S \\ H_R^S \end{array} \right. = \frac{3r_o}{2\mu_0} \int_S B_r \cos \theta dS \left\{ \begin{array}{l} \text{for } \theta < \frac{\pi}{2} \text{ and } B_r \cos \theta > 0, \\ \text{for } \theta < \frac{\pi}{2} \text{ and } B_r \cos \theta < 0, \\ \text{for } \theta > \frac{\pi}{2} \text{ and } B_r \cos \theta > 0, \\ \text{for } \theta > \frac{\pi}{2} \text{ and } B_r \cos \theta < 0. \end{array} \right. \qquad (22)$$

We then quantified the flux hemisphericities by the ratios:

$$H_N^{S/N} = \frac{H_N^S - H_N^N}{H_N^S + H_N^N} \qquad (23)$$

for the normal flux and

$$H_R^{N/S} = \frac{H_R^N - H_R^S}{H_R^N + H_R^S} \qquad (24)$$

for the reversed flux. To give intuition to the way these quantities work, $H_N^{S/N} > 0$ means more normal flux in the southern hemisphere, while $H_R^{N/S} > 0$ means more reversed flux in the northern hemisphere, both favoring surface intensity minima in the northern hemisphere in this case (Terra-Nova et al., 2017, 2019).
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Histograms of the latitudinal distribution of geomagnetic surface intensity minima in panel (a) the satellite era using the CHAOS7.10 model (Finlay et al., 2020) over the period 1998–2021 with a time-step of 1 yr, (b) the observatory era using the COV-OBS. x2 model (Huder et al., 2020) over the period 1840–2020 with a time-step of 2 years, (c) the archeomagnetism era using the COV-LAKE model (Hellio & Gillet, 2018) over the period −1,000 to 1,800 with a time-step of 100 years and (d) the archeomagnetism era using the pfm9k.2 model (Nilsson et al., 2022) over the period −7,000 to 1,800 with a time-step of 50 years. The satellite and observatory era models were truncated at $ℓ_{max}$ = 14 while the archeomagnetism era models were truncated at $ℓ_{max}$ = 6. Vertical lines denote models' corresponding $\langle P_{\lambda} \rangle$ values.</figcaption></figure>

# 3. Results

## 3.1. Geomagnetic Constraints

Figure 3 shows the latitudinal distribution of surface intensity minima for a satellite era geomagnetic field model (Finlay et al., 2020) in (a), an observatory era field model (Huder et al., 2020) in (b) and two archeomagnetic field models COV-LAKE by Hellio and Gillet (2018) in (c) and pfm9k.2 by Nilsson et al. (2022) in (d). For CHAOS-7.10, two minima are located at mid latitudes of the Southern Hemisphere reflected by a $\langle P_{\lambda} \rangle \approx -32^{\circ}$. In the observatory era the majority of minima are also located in the Southern Hemisphere, but some minima are found in the Northern Hemisphere as well, which is reflected in the reduced $\langle P_{\lambda} \rangle$ value of $\approx -15^{\circ}$. In the archeomagnetic era based on the pfm9k.2 model (Nilsson et al., 2022) the latitude of surface intensity minima distribution is close to the equator (Figure 3d), in agreement with the behavior of various archeomagnetic field models (Terra-Nova et al., 2019). However, due to large data uncertainties, the results for the Holocene may vary depending on the archeomagnetic field model. For example, in the COV-LAKE model (Hellio & Gillet, 2018) there are times, especially between 600 and 1400 CE, in which surface intensity minima appear at large latitudes (Figure 3c). However, the level of dipolarity was significantly lower when the latitude of surface intensity minima was large. For the entire period of the model the time-average relative dipolarity is $f_{dip} = 0.75$, as opposed to $f_{dip} = 0.43$ in the time interval 600–1400 CE. This is in contrast to present day conditions in which large values characterize both the SAA latitude (Figures 3a and 3b) as well as the relative dipolarity (Table 1).

These values serve as the observational constraints for our study. Our goals are two-fold. First, to obtain fundamental understanding for the dependence of the persistent latitude of surface intensity minima on the pattern of the CMB heat flux. Second, to identify the pattern ingredient that may reproduce the observed latitudes of the SAA. For these purposes we analyze dynamo models with synthetic large-scale outer boundary heat flux patterns.

## 3.2. Persistent Latitude of Surface Intensity Minima in the Dynamo Models

In eight out of the 18 numerical dynamo simulations explored in this study, the histograms of the latitudinal distribution of surface intensity minima are characterized by a single peak centered at about 1° latitude from the equator or less. Out of the other ten, seven simulations present geomagnetic reversals (Table 1). Figure 4 shows histograms of the latitudinal distribution of surface intensity minima for eight of our numerical dynamos simulations, with three of those presenting geomagnetic reversals.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Histograms of the latitudinal distribution of surface intensity minima in some of the dynamos models. Black dots denote the $\langle P_\lambda \rangle$ values. Dashed blue lines in panels (f, g) are fits to two Gaussians and the dashed black lines are their sums. Red and green vertical lines denote $\langle P_\lambda \rangle$ values of the geomagnetic field models CHAOS7.10 (Finlay et al., 2020) and COV-OBS. x2 (Huder et al., 2020), see Figure 3. Point-dashed lines are merely mirrored values of $\langle P_\lambda \rangle$ to the northern hemisphere. R marks reversing dynamo simulations.</figcaption></figure>

Case U with a uniform outer boundary heat flux is our reference model. The distribution resembles a Gaussian, with a high level of equatorial symmetry. Its mean latitude of surface intensity minima $\langle P_\lambda \rangle \sim 0.05^\circ$ (Table 1 and Figure 4a) is very close to zero. The deviation from zero is merely due to the finite simulation time, so this value may serve as a practical error bar for this quantity. Its standard deviation $\sigma \approx 3^\circ$ (Figure 4a) is also low, demonstrating proximity of the surface intensity minima to the equator essentially at all times. Therefore, despite careful choices of E, Ra and Pm values (Terra-Nova et al., 2019), without boundary heterogeneity it seems practically impossible to reproduce the present-day latitude of the SAA.

Equatorially anti-symmetric patterns, that is, axial with odd degrees, are natural CMB heat flux candidates for producing significant deviation of any geomagnetic feature from the equator. Starting from the largest scale anti-symmetric single harmonic, case H10b with an imposed $-Y_1^0$ outer boundary heat flux has surface intensity minima distribution displaced from the equator by a considerable amount $\langle P_\lambda \rangle \approx 7^\circ$ (Figure 4d). The preferred latitude is at the Northern hemisphere where the outer boundary heat flux is negative (Figure 2). Moreover, case H10b has more than three times larger $\sigma$ values than case U, that is, it allows a significant likelihood to find surface intensity minima inside a relatively high latitude interval $\langle P_\lambda \rangle \pm \sigma = [-2.61; 18.09]$. Case H10a has a much larger $\langle P_\lambda \rangle \approx 22^\circ$ and $\sigma \approx 27^\circ$ values (Figure 4c) that cover with a much larger probability the observed latitude of the SAA. In fact, the tail of the histogram in this case essentially covers the entire globe, that is, this
<table><caption>**Table 2**

Two Gaussians Fit to Models With Bimodal Distribution of Surface Intensity Minima</caption><thead><tr><th>CASE</th><th>⟨P<sub>λ₁</sub>⟩</th><th>⟨P<sub>λ₂</sub>⟩</th><th>W⟨P<sub>λ₁</sub>⟩/W⟨P<sub>λ₂</sub>⟩</th></tr></thead><tbody><tr><td>H30b</td><td>-53.56 ± 11.09</td><td>18.98 ± 29.58</td><td>4.36</td></tr><tr><td>H30c</td><td>-44.96 ± 16.49</td><td>15.20 ± 30.22</td><td>2.52</td></tr><tr><td>H30d</td><td>-49.84 ± 16.04</td><td>21.85 ± 28.42</td><td>2.69</td></tr><tr><td>H3020b</td><td>24.74 ± 22.26</td><td>-33.60 ± 19.63</td><td>2.57</td></tr></tbody></table>

Note. $⟨P_{λ₁}⟩$ and $⟨P_{λ₂}⟩$ are the persistent latitudes of the two peaks surface intensity minima histograms and $W⟨P_{λ₁}⟩/W⟨P_{λ₂}⟩$ is the ratio between their weighted frequencies.

model may reproduce a surface intensity minimum at any latitude (Figure 4c). However, in this case the dipolarity drops to much below Earth-like values (Table 1).

Hemisphericity may also be obtained by more complex anti-symmetric patterns. Cases H30b and H30c correspond to imposed outer boundary heat flux patterns of the higher odd degree 3, which is the dominant single harmonic of the axial pattern of observed lower mantle seismic anomalies (Terra-Nova et al., 2019). As with case H10b, cases H30b and H30c also have much larger $⟨P_λ⟩$ ≈ 15° and ≈ 5° values, respectively (Figures 4f and 4g) than case U. In addition, cases H30b and H30c present the most spread distributions of surface intensity minima latitudes of the models in Figure 4, with large $σ$ values of ≈ 34° and ≈ 36°, respectively. Hence cases H30b and H30c have high probability (≈ 34%) to accommodate surface intensity minima at large latitudes of $⟨P_λ⟩ + σ ≈ 15°–50° and $⟨P_λ⟩ + σ ≈ 5°–41°$, respectively. A

second characteristic in the histograms of these order 3 boundary heterogeneity models is the emergence of a second peak (Figures 4f and 4g). We fit the two peaks of the latitudinal distribution of surface intensity minima using a two components Gaussian mixture model (see Section 2.2). For both cases H30b and H30c, the peak associated with the larger amplitude $⟨P_λ⟩$ Gaussian is the one closer to $⟨P_λ⟩$ whilst the secondary peak $⟨P_{λ₂}⟩$ is at much higher latitudes ($⟨P_{λ₂}⟩ ≤ -40°$) and at the opposite hemisphere to $⟨P_λ⟩$ (see Figures 4f and 4g, Table 2). In both cases the secondary peaks are at far larger latitudes than that of the SAA. However, cases H30b and H30c are reversing models with low dipolarity (Table 1) hence their field morphologies are far from Earth-like.

Motivated by the dominance of the $Y_2^2$ harmonic in tomographic models of lower mantle seismic anomalies (e.g., Masters et al., 2000), we examined case H22 with an imposed $Y_2^2$ harmonic as an outer boundary heat flux pattern. Not surprisingly, this non-axial equatorially-symmetric pattern yields a very small $⟨P_λ⟩$ value (Figure 4b), as in case U. Likewise, case H20a with a polar cooling CMB heat flux pattern has a practically zero $⟨P_λ⟩$, as expected (Figure 4). The models in Figure 4 that present $⟨P_λ⟩$ smaller than unity have also very low $σ$ values (below ≈ 5.7), hence there is a small probability (≈ 16%) to find surface intensity minima beyond latitude ≈ 5.7° or even less for cases U, H20a and H22. Case H20b with a $-Y_2^0$ outer boundary heat flux pattern corresponding to equatorial cooling has the strongest tendency to reverse (e.g., Frasson et al., 2025; Olson et al., 2010) hence its large $⟨P_λ⟩$ is probably due to lack of temporal convergence in the calculation of the latitude of surface intensity minima despite the long run (Table 1).

We also examined sums of two axial single harmonics. Adding $Y_2^0$ to either $-Y_1^0$ or $-Y_3^0$ leads to latitudinal profiles that somewhat resemble the tomographic profile (Figure 2). Comparison of case H3020a with a sum of $Y_2^0$ and $-Y_3^0$ as an imposed CMB heat flux pattern to case H30c with a $-Y_3^0$ boundary heterogeneity with all other parameters identical shows that the addition of $Y_2^0$ stabilizes the dipole (as expected from polar cooling) but the large $⟨P_λ⟩$ is lost (Figure 4h, Table 1). In contrast, when combining $-Y_3^0$ and $Y_2^2$ for the CMB heat flux, the resulting magnetic field is not significantly altered. Cases H30a and H3022a are both dipole dominated with similarly large $f_{dip}$ and surface intensity minima very close to the equator, whereas both cases H30c and H3022b have large latitudes of surface intensity minima but are multipolar (Table 1). Lastly, combining $Y_2^0$ and $-Y_1^0$ in case H1020 compared to case H10b causes the loss of the large $⟨P_λ⟩$ in the latter model (Table 1). We conclude that symmetric heat flux patterns such as $Y_2^0$ and $Y_2^2$ lead to a strong equatorial tendency of surface intensity minima, even when combined with an anti-symmetric pattern.

To reveal the impact of the choice of the convection vigor $Ra$ in producing large surface intensity minima, we compare two pairs of models, cases H30a versus H30c and cases H30b versus H30d (Table 1). Increasing $Ra$ from the dipole-dominated model H30a leads to a substantial increase in $⟨P_λ⟩$ and $σ$ for case H30c, but a loss of the Earth-like dipolarity. For a larger amplitude of boundary heterogeneity of $q^* = 2$, when the lower $Ra$ case H30b is already characterized by low dipolarity and large $⟨P_λ⟩$ and $σ$, increasing $Ra$ results in moderate changes in the already large $⟨P_λ⟩$ and $σ$ case H30d and a further decrease in $⟨f_{dip}⟩$. In summary, increasing $Ra$ while keeping all other parameters fixed may increase the latitude of surface intensity minima from too close to the equator to the large observed SAA latitude, but often at the cost of losing Earth-like dipolarity.
Next we explore the role of the amplitude of the outer boundary heat flux heterogeneity $q^*$ in producing large latitude of surface intensity minima. Again we rely on the models with a $-Y_3^0$ CMB heat flux pattern. We compare two pairs of models, this time cases H30a versus H30b with a lower Ra and cases H30c versus H30d with a larger Ra. Within each pair $q^*$ is increased. Increasing $q^*$ promotes increases in both $\langle P_\lambda \rangle$ and its standard deviation for both pairs of models. However, the increase in $q^*$ might also lead to reversals and diminishing $\langle f_{dip} \rangle$ (see cases H30a and H30b in Table 1). Interestingly, from case H30c to H30d the increase in $q^*$ leads to an increase in $\langle f_{dip} \rangle$, albeit insufficiently to reach Earth-like dipolarity. We conclude that increasing $q^*$ with an equatorially anti-symmetric heat flux pattern effectively increases the latitude of surface intensity minima as well as its temporal variability. However, it might also destabilize the dipole.

Cases H10a and H10b with a $-Y_1^0$ CMB heat flux differ in both Ra and $q^*$. Case H10a with lower Ra and larger $q^*$ has very large latitude of surface intensity minima that easily explains the location of the SAA but its dipolarity is too low. In contrast, both requirements are satisfactory for the larger Ra and lower $q^*$ case H10b, with fairly large latitude of surface intensity minima and acceptable dipolarity. Overall, the effects of increasing Ra and $q^*$ seem to be similar, that is, stronger convection and stronger boundary heterogeneity both allow for larger latitudes of surface intensity minima but often at the expense of less Earth-like field morphologies. Hence a fine tuned combination of these two parameters might be required to satisfy both the observed SAA latitude and the dipolar morphology of the geomagnetic field.

Our results indicate that obtaining large latitude of surface intensity minima requires field hemisphericity. It is well established that large dipolarity requires strong rotational effects (e.g., Christensen & Aubert, 2006). Therefore, this motivates exploring dynamo models with large amplitude of anti-symmetric outer boundary heat flux heterogeneity and a smaller Ekman number. Nevertheless, those faster rotating models still need to be located close to the limit between the non-reversing and reversing dynamos (Terra-Nova et al., 2019), so the convection vigor (Ra) should also be increased. Due to the increasing computational cost and the need for sufficiently long runs for meaningful statistical analysis, here we first present only two cases with a lower Ekman number of $E = 3 \times 10^{-5}$ and a relatively large amplitude of outer boundary heat flux heterogeneity of $q^* = 2$ (Table 1). In case H30e we imposed a $-Y_3^0$ heat flux pattern, while in case H3020b we combined a sum of $-Y_3^0$ and $Y_2^0$, in both cases with a larger Ra number compared to their corresponding $E = 1 \times 10^{-4}$ cases H30d and H3020a, respectively. The resulting histograms of surface intensity minima latitudes are shown in Figure 5. The histogram of case H30e exhibits a peak of $P_\lambda$ near the equator. However, its northern hemisphere tail extends beyond the average value based on the historical field (Figure 5a), as reflected by the large standard deviation of the latitude of surface intensity minima in this case (Table 1). In case H3020b surface intensity minima can be found in any latitude. The histogram is characterized by a major peak at mid latitudes of the northern hemisphere and a secondary peak at mid latitudes of the southern hemisphere practically where the SAA latitude resides (Figure 5b). However, this dynamo model is reversing with a low dipolarity (Table 1).

Next we investigate the dependence of the obtained latitude of surface intensity minima on $Pm$. Departing from model H30e with large amplitude of anti-symmetric outer boundary heat flux heterogeneity and a relatively small Ekman number, we explore two dynamo models with increased $Pm$ values (H30f and H30g). Somewhat counterintuitively, increasing $Pm$ leads to smaller magnetic to kinetic energy ratios $M$ (Table 1), due to non-trivial effects of the imposed large-scale equatorially anti-symmetric CMB heat flux pattern. Figure 6 shows that the surface intensity minima latitudinal distributions in the larger $Pm$ models are characterized by northern hemisphere tails which extend beyond the average value based on the historical field. A faint hint for such a tail can already be observed for the corresponding lower $Pm$ model H30e (Figure 5a). However, the frequency of snapshots with large latitude of surface intensity minima in the larger $Pm$ dynamo models is far larger (Figure 6). Indeed, these models have by far the largest standard deviations of the latitude of surface intensity minima among the non-reversing simulations (Table 1).

It is possible that the difficulties in reproducing the large latitude of surface intensity minima are due to the fact that our dynamo models are too far from Earth-like conditions. To address this issue, we further analyzed a timeseries of magnetic field output corresponding to ≈ 14 kyr from the numerical dynamo simulation 71% toward the path for Earth-likeness of Aubert and Gillet (2021). This dynamo model satisfies the Earth-like morphological criteria of Christensen et al. (2010) while approaching Earth-like control parameters, in particular a significantly lower Ekman number. However, the results in Figure 7 show that the surface intensity minima are very close to
<figure><img src="image_7.png" /><figcaption>**Figure 5.** As in Figure 4 for the dynamos models with $E = 3 \times 10^{-5}$ and $Pm = 2$ for two different heat flux patterns.</figcaption></figure>

the equator with $\langle P_{\lambda} \rangle < 1^{\circ}$, and their temporal variability is also low. Indeed, this equatorial location of surface intensity minima in a dynamo model with a very low Ekman number is in agreement with the power law of Terra-Nova et al. (2019). The 71% path dynamo therefore fails to reproduce the present latitude of the SAA.

We also explored eight forecast predictions introduced by Aubert (2023) based on the 100% percent of path dynamo model. Although these assimilated fields were initialized by the observed geomagnetic field in 2000 AD which has $P_{\lambda} \approx 25^{\circ}$S, the models predict significantly lower $P_{\lambda}$ values of $\sim 10^{\circ}$S for the next $\approx 50$ years. This compromise result between observations and modeling information (Rogers et al., 2025) further illustrates that the path dynamo models (Aubert et al., 2017) are characterized by low latitudes of surface intensity minima.

### 3.3. Instantaneous Latitude of Surface Intensity Minima and Dipolarity

Table 1 shows that the dynamo models with large time-average latitude of surface intensity minima are usually characterized by low time-average dipolarity, whereas models with Earth-like $\langle f_{dip} \rangle$ often fail in reproducing the observed $\langle P_{\lambda} \rangle$. However, this difficulty to reconcile dipole dominance and large latitude of surface intensity minima might appear in the long-term time-average properties of the dynamo models, while it is possible that some snapshots may satisfy both observational constraints.

To examine whether this is indeed the case, we explore the relation between $P_{\lambda}$ and $f_{dip}$ at long, heavily sampled timeseries of snapshots. We present results from three dynamo models (synthesis of the results for the entire timeseries of all dynamo models are summarized in Table 3). For visual clarity, in Figure 8 we present relatively short time windows of 20 advection times, corresponding to about 2.8 kyr for the geodynamo. The red and green horizontal lines denote the time-average latitudes of the observed geomagnetic surface intensity minima based on the satellite and observatory eras fields, respectively. A dynamo model snapshot is considered Earth-like if its $P_{\lambda}$ value exceeds these horizontal lines. Furthermore, the purple horizontal lines denote the critical relative dipolarity. To be considered Earth-like, a snapshot should exhibit both large enough $P_{\lambda}$ (exceeding the red or green horizontal lines) and $f_{dip}$ values above 0.5 (exceeding the purple horizontal lines). These successful periods are marked by gray boxes in Figure 8.

The dynamo model with a uniform CMB heat flux (case U, Figure 8a) is characterized by a large time-average dipolarity of $\langle f_{dip} \rangle = 0.72$, very close to the satellite and observatory era geomagnetic field values, with weak temporal variation, that is, at most snapshots the field is dipole dominated. However, in this model the surface intensity minima are strictly confined to the equatorial region. No snapshot of case U has $P_{\lambda}$ larger than $15^{\circ}$.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** As in Figure 4 for dynamos models with $E = 3 \times 10^{-5}$ and increased $Pm$ values of 4 (a) and 6 (b).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** As in Figure 4 based on a long timeseries of output from the 71% path dynamo model of Aubert and Gillet (2021).</figcaption></figure>

In the reversing dynamo model case H30c, a wide distribution of $P_λ$ values is found, including a fair fraction of snapshots with surface intensity minima at higher latitudes than those that characterize the geomagnetic field (Figure 8c). However, $\langle f_{dip} \rangle$ in case H30c is very small (see Table 1) and snapshots with $f_{dip} > 0.5$ are very rare.

Case H10b contains large enough values for both $P_λ$ and $f_{dip}$ (Figure 8e). The $P_λ$ values exhibit a clear tendency toward the northern hemisphere where the imposed outer boundary heat flux anomaly is negative. The geomagnetic threshold latitude of surface intensity minima is reached at several snapshots. The dipolarity is strongly time dependent, often exceeding 0.5. Some snapshots have both large enough $P_λ$ and $f_{dip}$ (Figure 8e).

The difficulty in reconciling large latitude of surface intensity minima with large Earth-like dipolarity might be related to the way $P_λ$ is calculated. As already speculated by Terra-Nova et al. (2019), at a given snapshot a large latitude of a surface intensity minimum might be canceled by another minimum at the other hemisphere or at least reduced by an equatorial minimum when evaluated by the weighted average in $P_λ$. We therefore repeated the analysis of Figures 8a, 8c and 8e, but instead of $P_λ$ we consider simply the single latitude of the global minimum $\lambda_{F_{min}}^g$ of each snapshot of our dynamo models (Figures 8b, 8d and 8f).

The first thing to note when comparing Figures 8a, 8c and 8e to Figures 8b, 8d and 8f concerns the observed geomagnetic field. The quantity $P_λ$ is a weighted average of several local surface intensity minima, hence it is typically lower than $\lambda_{F_{min}}^g$. The observatory era geomagnetic $\langle P_λ \rangle$ value is lower than the satellite era value because at the beginning of the observatory era some local minima appeared in the northern hemisphere (Figure 3b). In contrast, the satellite and observatory eras $\langle \lambda_{F_{min}}^g \rangle$ capture the SAA latitude and are therefore practically identical. In the dynamo models, $\lambda_{F_{min}}^g$ may substantially exceed $P_λ$ at some snapshots (compare Figure 8d vs. Figure 8c and Figure 8f vs. Figure 8e). However, the difficulty in reproducing large $\lambda_{F_{min}}^g$ for snapshots with $f_{dip} > 0.5$ persists.

Table 3 summarizes the results. Note that out of the 18 dynamo models explored in this study (Table 1), in 12 models some snapshots with both acceptable latitude of surface intensity minima and acceptable dipolarity appear (models with no such snapshots are simply not reported in Table 3). In general, such snapshots are rare—the fraction of time in which both conditions are met is about ~4.5% at most. The occurrence is more often when considering $\lambda_{F_{min}}^g$ rather than $P_λ$ when comparing with the satellite era geomagnetic field. In contrast, the occurrence is more often when considering $P_λ$ rather than $\lambda_{F_{min}}^g$ when comparing with the historical geomagnetic field. Finally, a lower dipolarity threshold of 0.45 obviously allows for more dynamo models with acceptable snapshots and larger fractions of time for each dynamo model. In case H10b these fractions of time roughly double with $f_{dip} > 0.45$ compared to with $f_{dip} > 0.5$, reaching ~8.6% based on $\lambda_{F_{min}}^g$ and the historical field model (Table 3).

Finally, we note that gradually increasing $Pm$ while keeping all other control parameters and outer boundary heat flux unchanged from $Pm = 2$ (case H30e) to $Pm = 4$ (case H30e) to $Pm = 6$ (case H30g) gives the largest latitudes of surface intensity minima (see Table 1, Figures 5a and 6). However, the relative dipolarity decreases for larger $Pm$ dynamo models, resulting in a small fraction of time when the field is both sufficiently dipolar and the latitude of surface intensity minima is sufficiently large when requiring

<table><caption>**Table 3**
Percentage of Time That Dynamo Models Mimic the Large Latitude of Surface Intensity Minima</caption><thead><tr><th rowspan="2">CASE</th><th colspan="2">CHAOS7.10</th><th colspan="2">COV-OBS.x2</th></tr><tr><th>T(P<sub>λ</sub>,f<sub>dip</sub>)</th><th>T(λ<sub>F<sub>min</sub></sub><sup>g</sup>,f<sub>dip</sub>)</th><th>T(P<sub>λ</sub>,f<sub>dip</sub>)</th><th>T(λ<sub>F<sub>min</sub></sub><sup>g</sup>,f<sub>dip</sub>)</th></tr></thead><tbody><tr><td></td><td colspan="4" style="text-align: center;">0.75 &gt; f<sub>dip</sub> &gt; 0.50</td></tr><tr><td>H10b</td><td>0.11</td><td>1.52</td><td>4.57</td><td>2.16</td></tr><tr><td>H1020</td><td>0.00</td><td>0.00</td><td>0.29</td><td>0.00</td></tr><tr><td>H22</td><td>0.00</td><td>0.01</td><td>0.07</td><td>0.03</td></tr><tr><td>H30a</td><td>0.00</td><td>0.00</td><td>0.47</td><td>0.00</td></tr><tr><td>H30b</td><td>0.03</td><td>0.04</td><td>0.56</td><td>0.04</td></tr><tr><td>H30c</td><td>0.00</td><td>0.00</td><td>0.01</td><td>0.00</td></tr><tr><td>H30d</td><td>0.00</td><td>0.10</td><td>0.15</td><td>0.13</td></tr><tr><td>H30e</td><td>0.05</td><td>1.07</td><td>2.11</td><td>1.10</td></tr><tr><td>H30f</td><td>0.00</td><td>0.15</td><td>0.15</td><td>0.20</td></tr><tr><td></td><td colspan="4" style="text-align: center;">0.75 &gt; f<sub>dip</sub> &gt; 0.45</td></tr><tr><td>H10b</td><td>0.22</td><td>2.97</td><td>8.61</td><td>4.04</td></tr><tr><td>H1020</td><td>0.00</td><td>0.00</td><td>0.29</td><td>0.00</td></tr><tr><td>H20a</td><td>0.00</td><td>0.00</td><td>0.06</td><td>0.00</td></tr><tr><td>H22</td><td>0.00</td><td>0.01</td><td>0.07</td><td>0.03</td></tr><tr><td>H30a</td><td>0.00</td><td>0.00</td><td>0.47</td><td>0.00</td></tr><tr><td>H30b</td><td>0.51</td><td>0.99</td><td>1.97</td><td>1.02</td></tr><tr><td>H30c</td><td>0.00</td><td>0.03</td><td>0.15</td><td>0.03</td></tr><tr><td>H30d</td><td>0.07</td><td>0.29</td><td>0.66</td><td>0.34</td></tr><tr><td>H30e</td><td>0.05</td><td>1.07</td><td>2.45</td><td>1.10</td></tr><tr><td>H30f</td><td>0.09</td><td>2.52</td><td>2.37</td><td>2.60</td></tr><tr><td>H30g</td><td>0.01</td><td>0.53</td><td>0.92</td><td>0.54</td></tr><tr><td>H3022b</td><td>0.00</td><td>0.01</td><td>0.01</td><td>0.01</td></tr></tbody></table>

Note. $T$ denotes the percentage of snapshots that have both large enough surface intensity minima based on the satellite or observatory era geomagnetic field models and large enough dipolarity based on the threshold value of 0.5 (top) or 0.45 (bottom). Dynamo models with all $T$ values equal to zero are not shown.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Dynamo models snapshots that satisfy instantaneously both large latitude of surface intensity minima and large relative dipolarity. Black curves (left Y-axis) denote $P_{\lambda}$ (a, c, e) and $\lambda_{F_{min}}^g$ (b, d, f), purple curves (right Y-axis) denote $f_{dip}$, all as functions of time in units of advection time for dynamo models U (a, b), H30c (c, d) and H10b (e, f). Red and green horizontal lines denote time-average latitudes of surface intensity minima based on the geomagnetic field models CHAOS7.10 (Finlay et al., 2020) and COV-OBS. x2 (Huder et al., 2020), respectively. Point-dashed red and green lines are merely mirrored values to the northern hemisphere. Purple horizontal lines denote $f_{dip} = 0.5$ (the critical Earth-like value considered here) as well as $P_{\lambda} = 0$ (i.e., equatorial location of surface intensity minima). Gray shaded boxes denote successful snapshots.</figcaption></figure>

$f_{dip} > 0.5$ (Table 3). In contrast, when allowing $f_{dip} > 0.45$ the model with $Pm = 4$ has a larger fraction of time satisfying both observations than the model with lower $Pm$ (Table 3).

## 3.4. Dynamical Origin

Persistent large latitude of any spatially-averaged dynamo quantity, including surface intensity minima, likely requires equatorial symmetry breaking. Figure 9 presents the time-average zonal profiles of temperature and flow for dynamo models with various outer boundary heat flux patterns. Due to strong rapid rotation effects, our reference dynamo model with uniform outer boundary heat flux case U (Figure 9a) is practically perfectly symmetric about the equator. Polar upwellings of warm fluid and downwellings of cold fluid at the edge of the TC are consistent with thermal wind (Lézin et al., 2023). In dynamo models with equatorially symmetric outer boundary heat flux patterns the equatorial symmetry prevails as expected (see e.g. Figure 9b). In the polar cooling
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Time-averaged non-dimensional zonal temperature anomaly $\delta T$ and azimuthal velocity $u_\phi$ (red/blue denote eastward/westward, respectively) with superimposed meridional circulation streamlines (solid/dashed contours denote anticlockwise/clockwise, respectively) for dynamo models with different outer boundary heat flux patterns.</figcaption></figure>

case H20a warmer fluid and westward flow at the equator (Figure 9e) is driven by the negative heat flux anomaly there (Figure 2).

Particularly strong symmetry breaking is obtained for the models with the equatorially anti-symmetric $-Y_1^0$ patterns (Figures 9c and 9d). In cases H10a and H10b, negative/positive heat flux anomalies in the northern/southern hemisphere (Figure 2) produce warmer/colder fluid and westward/eastward zonal flow, respectively (Figures 9c and 9d). Dynamo models with an imposed $-Y_3^0$ heat flux patterns (Figures 9f and 9g) have opposite temperature and flow patterns outside the TC compared to the cases with $-Y_1^0$ patterns (Figures 9c and 9d), simply due to the opposite signs of the heat flux anomalies of these two patterns outside the TC (Figure 2). Inside the TC the two patterns have the same signs of heat flux anomalies (see Figure 2) and indeed the thermal anomalies inside the TC at the top of the shell of case H30b (Figure 9f) resemble those of cases H10a and H10b (Figures 9c and 9d) with warm/cold northern/southern polar caps, respectively.

The thermal and velocity distributions described above (Figure 9) are clearly manifested in the morphologies of the time-average magnetic fields (Figure 10). For the calculation of the time-average magnetic field, we avoid cancellations due to reversals by accounting for the change in global polarity when the field reverses (as in e.g., Aubert et al., 2008; Kelly & Gubbins, 1997):

$$ \langle B_r \rangle = -\frac{1}{N_t} \sum_{i=1}^{N_t} \frac{g_1^0(i)}{|g_1^0(i)|} B_r(i) \quad (25) $$

where $N_t$ is the number of snapshots and the summation in Equation 25 is over all times $i$. Note the very good temporal convergence of all simulations, as evident for example, in the practically zonal time-average field patterns obtained for the imposed homogeneous and zonal outer boundary heat flux patterns (Figures 10a–10g). This adequate temporal convergence is not trivial for the small-scale low dipolarity reversing dynamo models.

The reference homogeneous outer boundary heat flux case U (Figure 10a) is characterized by a well-known time-average morphology of intense magnetic flux at high latitudes near the intersection of the TC with the outer boundary as well as polar minima (e.g., Olson & Christensen, 2002). For the cases with outer boundary heat flux
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Time-averaged non-dimensional radial component of the magnetic field at the CMB $B_r$ for dynamo models with different outer boundary heat flux patterns.</figcaption></figure>

heterogeneity, the imposed lateral outer boundary heat flux variability is well reflected in the time-average field morphologies. This is reminiscent in the persistent two intense patches at high latitudes of each hemisphere with a $Y_2^2$ pattern (Figure 10b) as was found by numerous studies (e.g., Aubert et al., 2007; Olson & Christensen, 2002). Largest scale hemispherical $-Y_1^0$ patterns concentrate the time-average field at polar regions of the large heat flux (in this case southern) hemisphere (Figures 10c and 10d), in agreement with numerical dynamo studies performed in the context of the past Martian dynamo (Amit et al., 2011). The $Y_2^0$ pattern reinforces the normal flux patches at high latitudes of both hemispheres (compare Figures 10a and 10e), as expected from polar cooling. Less attention was given in the geodynamo literature for a degree 3 zonal boundary pattern. We recover the expected pattern of intense normal flux at high latitudes of the southern hemisphere as well as reversed flux at low latitudes of the same hemisphere (Figures 10f and 10g), in concordance with peak heat flux at the south pole and local heat flux minimum at low latitudes of the southern hemisphere (Figure 2). Finally, the combined $-Y_3^0$ and $Y_2^0$ patterns, which resembles the zonal tomographic profile (Figure 2), exhibits aspects of both components. In case H3020, the $Y_2^0$ component induces a large level of equatorial symmetry, while the $-Y_3^0$ component is responsible for the more prominent polar minimum in the northern hemisphere (Figure 10h).
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Dynamo models' selected snapshots of non-dimensional magnetic surface intensity (left) and radial field at the CMB (right) with best reproduced values of the geomagnetic field model COV-OBS. x2 in 2020. (a) Best reproduced $f_{dip}$, (b) best reproduced $P_\lambda$ and (c) best reproduced pair of $f_{dip}$ and $P_\lambda$. Black circles on the left denote locations of surface intensity minima. Black and red contours on the left correspond to $1.2F_{min}$ and $1.08F_{min}$, respectively. In some cases these contours are absent if the local minima is relatively strong or practically invisible if the slope of the intensity is very steep. All magnetic fields are expanded until $\ell_{max} = 14$.</figcaption></figure>

### 3.5. Instantaneous Field Morphology on the Outer Boundary and Surface Intensity Minima

We are now in a position to relate the different dynamo quantities: imposed outer boundary heat flux, the long-term time-average temperature at the top of the shell, fluid flow, magnetic field on the outer boundary and eventually persistent locations of surface intensity minima. In case H10a, the northern hemisphere is characterized by negative heat flux (Figure 2), leading to positive temperature anomalies (Figure 9c), fluid upwelling and surface divergence, weak field on the outer boundary of the shell (Figure 10c) and significantly higher latitude surface intensity minima (Table 1). The surface intensity minima in this case appears to reside in the northern hemisphere away from the intense normal flux patches of high latitudes of the southern hemisphere rather than near the reversed flux patches of low latitudes of the same (southern) hemisphere (Figure 10c). Case H30c is less trivial because the southern hemisphere reversed flux patches are stronger (Figure 10g) due to the local minimum in the heat flux pattern (Figure 2). Nevertheless, again the effect of the normal flux patches overcomes that of the reversed flux patches, and the surface intensity minima appear in the northern hemisphere (Table 1).

Next we examine maps of the instantaneous geomagnetic field intensity at Earth's surface and the radial field at the CMB to further illustrate the relation between the two distributions in snapshots, focusing on the latitude of surface intensity minima. We consider selected snapshots that sample distinctive pairs of $P_\lambda$ and $f_{dip}$ values. Figure 11a shows a typical snapshot from case U. Here the $f_{dip}$ value is closest to that of the present geomagnetic field based on the time average of CHAOS7.10 among all snapshots of all our dynamo models. Indeed, the radial field on the CMB in this snapshot captures well the main morphological aspects of the present-day geomagnetic field, including high-latitude intense normal flux patches that result in the large dipolarity, reversed flux patches
including polar minima, and overall strong deviations from zonality. However, the surface intensity contains three local minima, all located very close to the equator (Figure 11a).

Figure 11b shows the snapshot with the closest $P_λ$ value to that of COV-OBS. x2 for the year 2020. Two local surface intensity minima are found, one significantly more anomalous (i.e., weaker intensity) at mid latitudes of the southern hemisphere and another much less weak at low latitudes of the northern hemisphere, hence the $P_λ$ value is mostly determined by the significantly weaker minimum in the southern hemisphere. In contrast to the recovery of the observed $P_λ$, the $f_{dip}$ value is far too low. The radial field at the CMB contains a strikingly large and intense reversed flux region at low- and mid-latitudes of the southern hemisphere (Figure 11b), which is responsible for both the low global dipolarity as well as for the location of the surface intensity minimum, in particular the large latitude of the latter.

Figure 11c shows the snapshot with closest pair of $P_λ$ and $f_{dip}$ values to those of COV-OBS. x2 for the year 2020 among all snapshots of all our dynamo models. We chose this snapshot by minimizing the following difference between the dynamo model snapshot and the observed:

$$ \left( \frac{P_{\lambda}^{\text{mod}} - P_{\lambda}^{\text{obs}}}{\max\langle P_{\lambda}^{\text{mod}} \rangle} \right)^2 + \left( \frac{f_{\text{dip}}^{\text{mod}} - f_{\text{dip}}^{\text{obs}}}{\max\langle f_{\text{dip}}^{\text{mod}} \rangle} \right)^2 \quad (26) $$

where max denotes maximum value and the superscripts mod and obs denote dynamo models snapshots and the observed geomagnetic field in 2020, respectively. Two local surface intensity minima are present in this snapshot, one with moderate intensity located at the equator and the weakest at mid latitudes of the Northern hemisphere. Despite the anti-symmetric outer boundary heat flux pattern imposed on this dynamo model, case H10b has a rather large time-average $f_{dip}$ = 0.51 (Table 1). In this snapshot the dipolarity $f_{dip}$ = 0.66 (Figure 11c) is very close to that of the present-day geomagnetic field. The imposed $-Y_1^0$ outer boundary heat flux, that is, southern cooling, leads to stronger concentration of intense magnetic flux at high latitudes of the colder southern hemisphere. Therefore, the weakest surface intensity minimum is found in the opposite (northern) hemisphere. The longitude of the surface intensity minimum appears at a longitude devoid of normal flux in the northern Hemisphere (Figure 11c).

Note that in both snapshots shown in Figures 11b and 11c the $P_λ$ values are reduced due to the presence of multiple surface intensity minima at distinctive latitudes. In fact, in Figure 11b the weakest minimum appears at a significantly larger latitude than that of the SAA but the other minimum is at mid latitudes of the (other) northern hemisphere. In Figure 11c the less weak minima is located at the equator, again reducing the weighted average $P_λ$. This motivates comparing the geomagnetic field and the dynamo models in terms of an alternative quantity, the latitude of the global surface intensity minimum $\lambda_{F_{min}}^g$. In Figure 12 we show maps of two snapshots with $\lambda_{F_{min}}^g$ values close to those of the geomagnetic field in the year 2020 (Figure 1), -24.5° and -26.4° for COV-OBS. x2 and CHAOS7.10, respectively. In Figure 12a the dipole is relatively weak, $f_{dip}$ = 0.39. The weakest surface intensity minimum appears above a broad region of reversed flux on the outer boundary of the shell. In Figure 12b both $f_{dip}$ and $\lambda_{F_{min}}^g$ have very close values to those of the present-day geomagnetic field. The weakest surface intensity minimum at mid latitudes of the northern hemisphere is related to a very strong RFP that extends from the equator to the edge of the TC.

## 3.6. Hemispherical Ratios of Magnetic Flux on the Core-Mantle Boundary and the Latitude of Surface Intensity Minima

In the previous section we described qualitatively the relation between distributions of reversed and normal magnetic flux on the outer boundary and the latitude of surface intensity minima. Here we attempt to quantify this relation. In Figure 13 we plot $P_λ$ versus the hemisphericities of normal (Equation 23) and reversed flux (Equation 24) for all snapshots of two dynamo models. The positive slopes (Figure 13 right) validate that more normal flux in the southern hemisphere or more reversed flux in the northern hemisphere lead to surface intensity minima in the northern hemisphere, as expected. However, when considering only one type of polarity (normal or reversed), linear fits forced to pass through the origin are very poor (see low or negative $R^2$ values in Figure 13 left and middle) because both flux types affect the location of surface intensity minima (Terra-Nova et al., 2017,
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Dynamo models' selected snapshots of non-dimensional magnetic surface intensity (left) and radial field at the CMB (right) with best reproduced values of the geomagnetic field model COV-OBS. x2 in 2020. (a) Best reproduced $\lambda_F^g$ and (b) best reproduced pair of $f_{dip}$ and $\lambda_{Fmin}^g$. Black circles on the left denote locations of surface intensity minima. Black and red contours on the left correspond to 1.2$F_{min}$ and 1.08$F_{min}$, respectively. In some cases these contours are absent if the local minima is relatively strong or practically invisible if the slope of the intensity is very steep. All magnetic fields are expanded until $\ell_{max} = 14$.</figcaption></figure>

2019). We also plot $P_\lambda$ versus their linear combinations by searching for best fits to $P_\lambda = \alpha(H_N^{S/N} + \beta H_R^{N/S})$. Here we obtain significantly improved fits (see larger positive $R^2$ values in Figure 13 right).

In the linear combinations, smaller than unity coefficients for the contributions of reversed flux hemisphericity indicate that normal flux affects the latitude of the surface intensity minima more than reversed flux, in agreement with results from dynamo models with a tomographic heat flux pattern (Terra-Nova et al., 2019). Note that the coefficient factoring the reversed flux contribution was neither constrained to be positive nor to be smaller than unity (see $\beta$ values in Table 4 for all dynamo models). However, the clouds of points demonstrate that such a simplified assessment is limited (see also low goodness of fits in Figure 13) and the specific configuration of magnetic flux might introduce some exotic relations between the radial field on the outer boundary and the latitude of the surface intensity minima. Note that the dispersion is larger for the smaller scale imposed outer boundary heat flux case H30c where the range of latitudes of surface intensity minima is far larger (compare Figures 4d and 4g). Nevertheless, the trends of all dynamo models are positive (Figure 13 right), with the exception of case H3020b which has the largest Ra rendering it the most turbulent case (Table 4). Finally, comparison between the $\langle P_\lambda \rangle$ and $\sigma$ values in Table 1 and the $\beta$ values in Table 4 suggests that reversed flux hemisphericity is relatively more/less influential in prescribing the latitude of surface intensity minima when it is large/small, respectively.

**4. Discussion**

Several measures may be used for quantifying the latitude of geomagnetic surface intensity minima. The spatio-temporal mean $\langle P_\lambda \rangle$, obtained by averaging over all minima at each snapshot, seems an obvious choice since it reflects the probability to find surface intensity minima at a certain latitude (Terra-Nova et al., 2019). A second candidate is a large standard deviation $\sigma$, which defines the limits of probable latitudes to find surface intensity within the interval $(\langle P_\lambda \rangle - \sigma; \langle P_\lambda \rangle + \sigma)$, hence may encompass large latitudes even if the mean value is located close to the geographic equator. In the dipole-dominated non-reversing dynamo models of Terra-Nova et al. (2019) with the imposed multi-harmonic tomographic outer boundary heat flux, all latitude histograms were characterized by a single peak. Here we showed that in some cases the histograms exhibit two peaks. In these cases we fitted a sum of two Gaussians in order to extract the latitude of the secondary peak, its relative amplitude and its standard deviation. Lastly, snapshots from the historical geomagnetic field and dynamo models often
<figure><img src="image_15.png" /><figcaption>**Figure 13.** 2D histograms of $P_{\lambda}$ (Equation 17) versus the hemispherical ratios of magnetic flux at the CMB $H_{N}^{S/N}$ (Equation 23) and $H_{R}^{S/N}$ (Equation 24) for all snapshots of two dynamos models. The color bar denotes the frequency of surface minima per 2D bin. Dashed curves are linear fits forced to pass through the origin.</figcaption></figure>

<table><caption>**Table 4**

Hemispherical Fits</caption><thead><tr><th>CASE</th><th>α</th><th>β</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td>U</td><td>30.53</td><td>0.17</td><td>0.28</td></tr><tr><td>H10a</td><td>77.52</td><td>0.55</td><td>0.14</td></tr><tr><td>H10b</td><td>70.86</td><td>0.34</td><td>0.50</td></tr><tr><td>H1020</td><td>49.68</td><td>0.30</td><td>0.34</td></tr><tr><td>H20a</td><td>50.53</td><td>0.24</td><td>0.30</td></tr><tr><td>H20b</td><td>55.61</td><td>0.36</td><td>0.43</td></tr><tr><td>H22</td><td>31.27</td><td>0.18</td><td>0.22</td></tr><tr><td>H30a</td><td>37.19</td><td>0.09</td><td>0.16</td></tr><tr><td>H30b</td><td>59.30</td><td>0.41</td><td>0.43</td></tr><tr><td>H30c</td><td>51.40</td><td>0.65</td><td>0.38</td></tr><tr><td>H30d</td><td>71.51</td><td>0.58</td><td>0.36</td></tr><tr><td>H30e</td><td>42.58</td><td>0.08</td><td>0.21</td></tr><tr><td>H30f</td><td>93.11</td><td>0.43</td><td>0.54</td></tr><tr><td>H30g</td><td>69.52</td><td>0.56</td><td>0.37</td></tr><tr><td>H3020a</td><td>35.57</td><td>0.17</td><td>0.15</td></tr><tr><td>H3020b</td><td>72.54</td><td>0.35</td><td>-0.23</td></tr><tr><td>H3022a</td><td>38.12</td><td>0.10</td><td>0.20</td></tr><tr><td>H3022b</td><td>55.13</td><td>0.53</td><td>0.39</td></tr></tbody></table>

Note. α and β are fitting parameters of the equation $P_{\lambda} = \alpha(H_{N}^{S/N} + \beta H_{R}^{N/S})$. R² denotes the goodness of the fit.

exhibit more than one surface intensity minimum (see $N_{F_{min}}$ in Table 1 of Terra-Nova et al., 2019, as well as in Table 1 of this study). Hence, multiple minima might lead to cancellations and reduction in the spatial mean of the latitude. Therefore, in addition we considered here the latitude of the weakest surface intensity $\lambda_{F_{min}}^g$. We found small differences between $\langle P_{\lambda}^g \rangle$ and $\langle P_{\lambda} \rangle$ values (Table 1), which provides confidence in the robustness of the quantification of the latitude of surface intensity minima. However, when comparing snapshots of the dynamo models with the satellite era geomagnetic field, considering $\lambda_{F_{min}}^g$ provides systematically larger fractions of acceptable pairs with $f_{dip}$ (Table 3 and Figure 8). It is therefore possible that the large southern latitude of $\langle P_{\lambda} \rangle \approx -31.86^{\circ}$ found for the present-day SAA is in fact due to the absence of a surface intensity minimum at the northern hemisphere or at the equator to balance it.

Reversing dynamo models may exhibit large latitudes of surface intensity minima (e.g., Figure 11b). The single peak of the histogram of case H10a is ~22° (Figure 4c and Table 1) while the secondary peaks of cases with a $-Y_3^0$ outer boundary heat flux may reach latitudes as high as ~50° (Figures 4f and 4g and Table 2). However, these reversing dynamo models are often characterized by low dipolarities (e.g., Christensen & Aubert, 2006). Conversely, in dipole-dominated non-reversing dynamo models, low latitudes of surface intensity minima (e.g., Figure 11a) are often found (Terra-Nova et al., 2019). Simultaneously satisfying both requirements, that is, large Earth-like dipolarity and large SAA-like latitude of surface intensity minima, is challenging. Within the non-reversing dynamo models, proximity to the transition to reversing dynamos favors large latitude of surface intensity minima (Terra-Nova et al., 2019). In addition, heterogeneous boundary conditions with
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Absolute value of time-average latitude of surface intensity minima |⟨$P_λ$⟩| (Equation 18) versus time-average relative dipolarity ⟨$f_{dip}$⟩ (Equation 12) for all dynamo models. Same color is used for a group of dynamo models with the same imposed outer boundary heat flux pattern, black edge-color symbols indicate dynamo models with $q^* = 2.0$ (see legend). Earth-like ranges are based on the geomagnetic field models CHAOS7.10 (Finlay et al., 2020) (point-dashed red horizontal line) and COV-OBS. x2 (point-dashed green horizontal line) for |⟨$P_λ$⟩| and on the numerical dynamo literature (e.g., Davies et al., 2022) for ⟨$f_{dip}$⟩. The gray shaded box represents Earth-like |⟨$P_λ$⟩| and ⟨$f_{dip}$⟩ values. Error bars denote standard deviations that represent time variability.</figcaption></figure>

marked north-south hemispherical anti-symmetry is key for reproducing large latitude of surface intensity minima. Therefore, a careful choice of main control parameters which localize a dynamo model close to the transition to reversals and in addition an equatorially anti-symmetric boundary condition seem to be the proper ingredients to obtain a dipole-dominated dynamo model with a large latitude of surface intensity minima, as is the case for the present-day geomagnetic field.

A synthesis of the time-average results is given in Figure 14. We find a trade-off between latitude of surface intensity minima and relative dipolarity. Some dynamo models persistently reproduce Earth-like large latitudes of surface intensity minima but are often multipolar. This is the case when large amplitude equatorially anti-symmetric outer boundary heat flux heterogeneity patterns are imposed (see black edge-color symbols in Figure 14). In contrast, other dynamo models recover large time-average Earth-like dipolarity, but the surface intensity minima in these models are often too close to the equator.

In addition to time-average properties, we also accounted for the temporal variability. In this study, we reproduced satisfactorily the historical field ⟨$P_λ$⟩ value (Figure 3b) when considering the interval (⟨$P_λ$⟩ − σ; ⟨$P_λ$⟩ + σ) while keeping reasonable values of averaged field dipolarity ⟨$f_{dip}$⟩ > 0.5 using a −$Y_1^0$ outer boundary heat flux pattern (case H10b in Figure 4d and Table 1). However, a careful examination of a long and heavily sampled timeseries of snapshots of this dynamo model reveals that instantaneous simultaneous recovery of the historical field dipolarity $f_{dip}$ and $P_λ$ is only

seldom reproduced. Moreover, when comparing with the $P_λ$ values of the satellite era field snapshots acceptable pairs of $P_λ$ and $f_{dip}$ are even more rare (Figures 8e, 8f and 11c and Table 3).

We also explored a combined outer boundary heat flux pattern of −$Y_3^0$ and $Y_2^0$, with an attempt to obtain large latitude of surface intensity minima by the former and large dipolarity by the latter. Comparison of cases H30c and H3020a with identical setup except for the addition of the polar cooling component in the latter shows that indeed the dipole was stabilized but the latitude of the surface intensity minima substantially decreased (Table 1). Adding a $Y_2^0$ pattern has little impact on the results (compare cases H3022b and H30c in Table 1), that is, this equatorially symmetric pattern produces no hemisphericity and has little impact on the dipolarity.

The trade-off between large latitude of surface intensity minima and large dipolarity is evident when considering the dependence of the results on the main convection vigor represented by Ra. In dynamo models with an anti-symmetric outer boundary heat flux pattern, increasing Ra (with all other parameters unchanged) yields larger latitudes of surface intensity minima, but on the expense of loss of dipole dominance (Table 1). Increasing amplitude of an anti-symmetric outer boundary heat flux heterogeneity $q^*$ also leads to larger latitudes of surface intensity minima. However, when the lower $q^*$ dynamo model was characterized by low dipolarity, increasing $q^*$ led to an increase in $f_{dip}$ (compare cases H30c and H30d in Table 1). This points to a non-monotonous behavior of the dynamo for changes in $q^*$. Indeed, Terra-Nova and Amit (2024) found for dynamo models with a tomographic CMB heat flux pattern an increase in reversibility with increasing $q^*$ for $q^* < 1$, but a reduction in reversibility with increasing $q^*$ for $q^* > 1$. The key to reconciling large latitude of surface intensity minima and large dipolarity may therefore be a large enough amplitude of imposed outer boundary heat flux heterogeneity.

Another possibility is that an anti-symmetric outer boundary heat flux combined with a lower Ekman number *E* which may stabilize the dipole would lead to numerical dynamo simulations that reproduce both the current large latitude of the SAA and Earth's dominant axial dipole. Due to the high computational cost, we only explored four models with a lower *E*. Although the non-reversing case H30e with a −$Y_3^0$ outer boundary heat flux pattern fails to produce a significant ⟨$P_λ$⟩, its standard deviation is about double that of the non-reversing case H30a with a larger *E* (Table 1). Moreover, the histogram of case H30e is asymmetric with a tail that extends beyond the average value based on the historical field (Figure 5a). These results hint that low *E* dynamos with a large amplitude of
equatorially anti-symmetric outer boundary heat flux heterogeneity may be promising in terms of reconciling large latitude of the surface intensity minima and large dipolarity.

Guided by the results of Terra-Nova et al. (2019), we placed our dynamo models close to the dipolar to reversing transition which is favorable for obtaining large latitude of surface intensity minima. To further augment the latitude of surface intensity minima, we considered equatorially anti-symmetric outer boundary heat flux patterns. The proximity to the transition renders the magnetic to kinetic energy ratio close to unity (Table 1). Increasing *Pm* with homogeneous boundary conditions tends to increase this ratio toward a strong field dynamo (Dormy, 2016; Schwaiger et al., 2019; Menu et al., 2020). However, in our models with an anti-symmetric CMB heat flux pattern increasing *Pm* tends to slightly decrease the energy ratio. It is possible that far from the dipolar to reversing transition increasing *Pm* would increase the field strongness. However, far from the transition it is likely that the surface intensity minima will remain very close to the equator (Terra-Nova et al., 2019).

The dynamical origin of surface intensity minima is unraveled by the relations between temperature, zonal flow, radial magnetic field at the CMB and these surface features. The coupling between temperature and azimuthal flow is described by the thermal wind equation (e.g., Aubert, 2005; Pedlosky, 1987) that governs dynamo models at large scales (Schwaiger et al., 2019). Positive/negative heat flux anomaly translates into colder/warmer fluid below the outer boundary (e.g., Gubbins, 2003) which prompts eastward/westward flow, respectively (Figure 9). The time-averaged magnetic field on the outer boundary is strong at regions of boundary-driven downwellings and reversed where boundary-driven upwellings prevail (Figure 10), in agreement with previous studies (e.g., Olson & Christensen, 2002).

The present-day latitude of the SAA is related to the vast reversed flux region below the South Atlantic (Figure 1). In general, snapshots from our dynamo models also exhibit surface intensity minima above regions where reversed flux prevails on the outer boundary (e.g., Figures 11c and 12b). However, the distribution of normal flux on the CMB also affects the locations of surface intensity minima (Terra-Nova et al., 2017). In fact, analysis of heavily sampled timeseries of snapshots from the dynamo models confirm that the surface intensity minima reside at the hemisphere with more reversed flux and with less normal flux, with the latter being more influential. However, the clouds of points in Figure 13 are evidence for the non-trivial relation between the magnetic flux on the CMB and the locations of surface intensity minima.

In summary, we showed that the latitude of the SAA is evidence for the strong impact of lower mantle thermal heterogeneity on the geodynamo. In order for a surface intensity minimum to reside in a latitude as high as that of the SAA, the CMB heat flux must have a substantial equatorially anti-symmetric component, which may contradict an equatorially symmetric core flow (e.g., Barrois et al., 2017; Finlay et al., 2023; Gillet et al., 2019). While many studies of geodynamo simulations often focus on the dominant $Y_2^2$ harmonic in mantle tomography models (e.g., Aubert et al., 2007; Olson & Christensen, 2002; Takahashi et al., 2008), the zonal part of these mantle tomography models is dominated by an equatorially anti-symmetric $-Y_3^0$ component (Figure 2). Our dynamo models with an imposed $-Y_3^0$ outer boundary heat flux pattern indeed produce very broad distributions of latitudes of surface intensity minima that easily cover the latitude of the present-day SAA. Simultaneously reproducing at the same dynamo model snapshot both the SAA latitude and the dipolarity of the present-day geomagnetic field is possible, though somewhat rare. Further exploration of dynamo models with large amplitudes of anti-symmetric outer boundary heat flux heterogeneity combined with lower Ekman numbers may improve the recovery of these two fundamental observational geomagnetic constraints.

## Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

## Data Availability Statement

The numerical code used to perform our magnetohydrodynamics spherical shell simulations is an open-source software called MagIC (Schaeffer, 2013; Wicht, 2002). It is developed on GitHub (URL: https://magic-sph.github.io/) and can be used, modified and redistributed under the terms of the GNU GPL v3 license. MagIC is available in this in-text data citation reference: Gastine et al. (2022). Data sets for this research are available in this in-text data citation reference: Terra-Nova and Amit (2025).
Acknowledgments

F.T.-N. was supported by the Centre
National d'Études Spatiales (CNES),
project MIMES. We acknowledge
financial support from the French Agence
Nationale de Recherche, project DYRE-
COMB (Grant ANR-22-CE49-0016-01).
We thank Guy Moebes for his help in
setting up the MagIC dynamo code
(https://github.com/magic-sph/magic) on
the GLICID (le Groupement Ligérien pour
le Calcul Intensif Distribué) computational
center where the numerical simulations for
this study were performed.

References

Alken, P., Thébault, E., Beggan, C., Amit, H., Aubert, J., Baerenzung, J., et al. (2021). International geomagnetic reference field: The thirteenth generation. *Earth Planets and Space*, 73(49), 49. https://doi.org/10.1186/s40623-020-01288-x

Amit, H., & Choblet, G. (2009). Mantle-driven geodynamo features—Effects of post-perovskite phase transition. *Earth Planets and Space*, *61*(11), 1255–1268. https://doi.org/10.1186/bf03352978

Amit, H., & Choblet, G. (2012). Mantle-driven geodynamo features—Effects of compositional and narrow D'' anomalies. *Earth Planet. Inter*, **190-191**, 34–43. https://doi.org/10.1016/j.pepi.2011.10.005

Amit, H., Choblet, G., Olson, P., Monteux, J., Deschamps, F., Langlais, B., & Tobie, G. (2015). Towards more realistic core-mantle boundary heat flux patterns: A source of diversity in planetary dynamos. *Progress in Earth and Planetary Science*, 2(1), 26. https://doi.org/10.1186/s40645-015-0056-3

Amit, H., Christensen, U. R., & Langlais, B. (2011). The influence of degree-1 mantle heterogeneity on the past dynamo of Mars. *Physics of the Earth and Planetary Interiors*, 189(1-2), 63–79. https://doi.org/10.1016/j.pepi.2011.07.008

Amit, H., Deschamps, F., & Choblet, G. (2015). Numerical dynamos with outer boundary heat flux inferred from probabilistic tomography—Consequences for latitudinal distribution of magnetic flux. *Geophysical Journal International*, 203(2), 840–855. https://doi.org/10.1093/gji/ggv332

Amit, H., Terra-Nova, F., Lézin, M., & Trindade, R. (2021). Non-monotonic growth and motion of the South Atlantic Anomaly. *Earth Planets and Space*, 73(38), 38. https://doi.org/10.1186/s40623-021-01356-w

Aubert, J. (2005). Steady zonal flows in spherical shell fluid dynamos. *Journal of Fluid Mechanics*, 542, 53–67. https://doi.org/10.1017/s0022112005006129

Aubert, J. (2015). Geomagnetic forecasts driven by thermal wind dynamics in the Earth's core. *Geophysical Journal International*, 203(3), 1738–1751. https://doi.org/10.1093/gji/ggv394

Aubert, J. (2023). State and evolution of the geodynamo from numerical models reaching the physical conditions of Earth's core. *Geophysical Journal International*, 235(1), 468–487. https://doi.org/10.1093/gji/ggad229

Aubert, J., Amit, H., & Hulot, G. (2007). Detecting thermal boundary control in surface flows from numerical dynamos. *Physics of the Earth and Planetary Interiors*, 160(2), 143–156. https://doi.org/10.1016/j.pepi.2006.11.003

Aubert, J., Amit, H., Hulot, G., & Olson, P. (2008). Thermo-chemical wind flows couple Earth's inner core growth to mantle heterogeneity. *Nature*, 454(7205), 758–761. https://doi.org/10.1038/nature07109

Aubert, J., Finlay, C. C., & Fournier, F. (2013). Bottom-up control of geomagnetic secular variation by the Earth's inner core. *Nature*, 502(7470), 219–223. https://doi.org/10.1038/nature12574

Aubert, J., Gastine, T., & Fournier, A. (2017). Spherical convective dynamos in the rapidly rotating asymptotic regime. *Journal of Fluid Mechanics*, 813, 558–593. https://doi.org/10.1017/jfm.2016.789

Aubert, J., & Gillet, N. (2021). The interplay of fast waves and slow convection in geodynamo simulations nearing Earth's core conditions. *Geophysical Journal International*, 214(1), 531–547.

Aubert, J., Labrosse, S., & Poitou, C. (2009). Modelling the paleo-evolution of the geodynamo. *Geophysical Journal International*, 179(3), 1414–1428. https://doi.org/10.1111/j.1365-246x.2009.04361.x

Barrois, O., Gillet, N., & Aubert, J. (2017). Contributions to the geomagnetic secular variation from a reanalysis of core surface dynamics. *Geophysical Journal International*, 211(1), 50–68. https://doi.org/10.1093/gji/ggx280

Bloxham, J., Gubbins, D., & Jackson, A. (1989). Geomagnetic secular variation. *Philosophical Transactions of the Royal Society of London*, A329, 415–502.

Brown, M., Korte, M., Holme, R., Wardinski, I., & Gunnarson, S. (2018). Earth's magnetic field is (probably) not reversing. *Proceedings of the national academy of sciences* (Vol. 115(20), pp. 5111–5116). https://doi.org/10.1073/pnas.1722110115

Cai, S., Doctor, R., Tauxe, L., Hendrickson, M., Hua, Q., Leroy, S., & Phon, K. (2021). Archaeomagnetic results from Cambodia in Southeast Asia: Evidence for possible low-latitude flux expulsion. *Proceedings of the national academy of sciences* (Vol. 118(11), p. e2022490118). https://doi.org/10.1073/pnas.2022490118

Campuzano, S. A., Gómez-Paccard, M., Pavón-Carrasco, F. J., & Osete, M. L. (2019). Emergence and evolution of the South Atlantic Anomaly revealed by the new paleomagnetic reconstruction SHAWQ2k. *Earth and Planetary Science Letters*, 512, 17–26. https://doi.org/10.1016/j.epsl. 2019.01.050

Choblet, G., Deschamps, F., Amit, H., & Lasbleis, M. (2023). Inferring the relationship between core-mantle heat flux and seismic tomography from mantle convection simulations. *Physics of the Earth and Planetary Interiors*, 342, 107072. https://doi.org/10.1016/j.pepi.2023.107072

Christensen, U. (2010). Dynamo scaling laws and applications to planets. *Space Science Reviews*, 152(1-4), 565–590. https://doi.org/10.1007/s11214-009-9553-2

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(1), 97–114. https://doi.org/10.1111/j.1365-246x.2006.03009.x

Christensen, U. R., Aubert, J., & Hulot, G. (2010). Conditions for Earth-like geodynamo models. *Earth and Planetary Science Letters*, 296(3-4), 487–496. https://doi.org/10.1016/j.epsl.2010.06.009

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

Domingos, J., Jault, D., Pais, M. A., & Mandea, M. (2017). The South Atlantic Anomaly throughout the solar cycle. *Earth and Planetary Science Letters*, 473, 154–163. https://doi.org/10.1016/j.epsl.2017.06.004

Dormy, E. (2016). Strong-field spherical dynamos. *Journal of Fluid Mechanics*, 789, 500–513. https://doi.org/10.1017/jfm.2015.747

Engbers, Y. A., Grappone, J. M., Mark, D. F., & Biggin, A. J. (2022). Low paleointensities and Ar/Ar ages from Saint Helena provide evidence for recurring magnetic field weaknesses in the South Atlantic. *Journal of Geophysical Research: Solid Earth*, 127(3), e2025JB023358. https://doi.org/10.1029/2021JB023358

Finlay, C., Kloss, C., Olsen, N., Hammer, M. D., Töffner-Clausen, L., Grayver, A., & Kuvshinov, A. (2020). The CHAOS-7 geomagnetic field model and observed changes in the South Atlantic Anomaly. *Earth Planets and Space*, 72(1), 156. https://doi.org/10.1186/s40623-020-01252-9

Finlay, C. C., Aubert, J., & Gillet, N. (2016). Gyre-driven decay of the Earth's magnetic dipole. *Nature communications*, 7(1), 10422. https://doi.org/10.1038/ncomms10422

Finlay, C. C., Gillet, N., Aubert, J., Livermore, P. W., & Jault, D. (2023). Gyres, jets and waves in the Earth's core. *Nature Reviews Earth & Environment*, 4(6), 377–392. https://doi.org/10.1038/s43017-023-00425-w

Finlay, C. C., Jackson, A., Gillet, N., & Olsen, N. (2012). Core surface magnetic field evolution 2000–2010. *Geophysical Journal International*, 189(2), 761–781. https://doi.org/10.1111/j.1365-246x.2012.05395.x
Frasson, T., Schaeffer, N., Nataf, H.-C., & Labrosse, S. (2025). Geomagnetic dipole stability and zonal flow changes controlled by mantle heat flux heterogeneities. *Geophysical Journal International*, 240(3), 1481–1504. https://doi.org/10.1093/gji/ggac457

Gastine, T., Barik, A., Rraynaud, Schwaiger, T., Putigny, B., Thtassin, J., et al. (2022). magic-sph/magic: Release magic 6.2 (v6.2). (Vol. 211). https://doi.org/10.5281/zenodo.7419346

Gillet, N., Huder, L., & Aubert, J. (2019). A reduced stochastic model of core surface dynamics based on geodynamo simulations. *Geophysical Journal International*, 219(1), 522–539. https://doi.org/10.1093/gji/ggz313

Glatzmaier, G., Coe, R., Hongre, L., & Roberts, P. (1999). The role of the Earth's mantle in controlling the frequency of geomagnetic reversals. *Nature*, 401(6756), 885–890. https://doi.org/10.1038/44776

Gubbins, D. (1987). Mechanism for geomagnetic polarity reversals. *Nature*, 326(6109), 167–169. https://doi.org/10.1038/326167a0

Gubbins, D. (2003). Thermal core-mantle interactions: Theory and observations. In V. Dehant, K. Creager, S. Karato, & S. Zatman (Eds.), *Earth's core: Dynamics, structure and rotation*. AGU Geodynamics Series—American Geophysical Union.

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

Hare, V. J., Tarduno, J. A., Huffman, T., Watkeys, M., Thebe, P. C., Manyanga, M., et al. (2018). New archeomagnetic directional records from Iron Age southern Africa (ca. 425–1550 CE) and implications for the South Atlantic Anomaly. *Geophysical Research Letters*, 45(3), 1361–1369. https://doi.org/10.1002/2017gl076007

Hartmann, G. A., Genevey, A., Gallet, Y., Trindade, R. I. F., Etchevarne, C., Le Goff, M., & Afonso, M. C. (2010). Archeointensity in Northeast Brazil over the past five centuries. *Earth and Planetary Science Letters*, 296(3–4), 340–352. https://doi.org/10.1016/j.epsl.2010.05.016

Hartmann, G. A., Genevey, A., Gallet, Y., Trindade, R. I. F., Le Goff, M., Najjar, R., et al. (2011). New historical archeointensity data from Brazil: Evidence for a large regional non-dipole field contribution over the past few centuries. *Earth and Planetary Science Letters*, 306(1–2), 66–76. https://doi.org/10.1016/j.epsl.2011.03.030

Hartmann, G. A., Poletti, W., Trindade, R. I., Ferreira, L. M., & Sanches, P. L. (2019). New archeointensity data from south Brazil and the influence of the South Atlantic Anomaly in South America. *Earth and Planetary Science Letters*, 512, 124–133. https://doi.org/10.1016/j.epsl.2019.01.049

He, F., Wei, Y., Maffei, S., Livermore, P. W., Davies, C. J., Mound, J., et al. (2021). Equatorial auroral records reveal dynamics of the paleo-West Pacific geomagnetic anomaly. *Proceedings of the national academy of sciences* (Vol. 118(20), p. e2026080118). https://doi.org/10.1073/pnas.2026080118

Heitzler, J. (2002). The future of the South Atlantic Anomaly and implications for radiation damage in space. *Journal of Atmospheric and Solar-Terrestrial Physics*, 64(16), 1701–1708. https://doi.org/10.1016/s1364-6826(02)00120-7

Hellio, G., & Gillet, N. (2018). Time-correlation-based regression of the geomagnetic field from archeological and sediment records. *Geophysical Journal International*, 214(3), 1585–1607. https://doi.org/10.1093/gji/ggy214

Huder, L., Gillet, N., Finlay, C., Hammer, M. D., & Tchoungui, H. (2020). COV-OBS.x2: 180 years of geomagnetic field evolution from ground-based and satellite observations. *Earth Planets and Space*, 72(1), 160. https://doi.org/10.1186/s40623-020-01194-2

Jackson, A. (2003). Intense equatorial flux spots on the surface of the Earth's core. *Nature*, 424(6950), 760–763. https://doi.org/10.1038/ nature01879

Kelly, P., & Gubbins, D. (1997). The geomagnetic field over the past 5 million years. *Geophysical Journal International*, 128(2), 315–330. https:// doi.org/10.1111/j.1365-246x.1997.tb01557.x

Korte, M., & Constable, C. G. (2011). Improving geomagnetic field reconstructions for 0-3 ka. *Physics of the Earth and Planetary Interiors*, 188(3–4), 247–259. https://doi.org/10.1016/j.pepi.2011.06.017

Korte, M., Donadini, F., & Constable, C. G. (2009). The geomagnetic field for 0–3 ka: 2. A new series of time-varying global models. *Journal of Geophysical Research*, 10, Q06008. https://doi.org/10.1029/2008GC002297

Kutzner, C., & Christensen, U. R. (2002). From stable dipolar towards reversing numerical dynamos. *Physics of the Earth and Planetary Interiors*, 131(1), 29–45. https://doi.org/10.1016/s0031-9201(02)00001-x

Lean, J. (2005). Living with a variable sun. *Physics Today*, 58(6), 32–36. https://doi.org/10.1063/1.1996472

Lézin, M., Amit, H., Terra-Nova, F., & Wardinski, I. (2023). Mantle-driven north-south dichotomy in geomagnetic polar minima. *Physics of the Earth and Planetary Interiors*, 337, 107000.

Li, H., Tang, J., Gai, C., Zhang, W., Humbert, F., Ni, Y., et al. (2023). Persistent influence of non-dipole geomagnetic field on East Asia over the past 4,000 years. *Journal of Geophysical Research: Solid Earth*, 128(12), e2023JB027359. https://doi.org/10.1029/2023jb027359

Licht, A., Hulot, G., Gallet, Y., & Thebault, E. (2013). Ensembles of low degree archeomagnetic field models for the past three millennia. *Physics of the Earth and Planetary Interiors*, 224, 38–67. https://doi.org/10.1016/j.pepi.2013.08.007

Luo, X., Liu, J., Wang, M., amd Schanner, H., Zhang, Y., Xie, P., et al. (2025). Low geomagnetic field intensity in southern China 6,000 years ago. *Geophysical Research Letters*, 52(9), e2024GL113552. https://doi.org/10.1029/2024gl113552

Masters, G., Laske, G., Bolton, H., & Dziewonski, A. (2000). The relative behavior of shear velocity, bulk sound velocity, and compressional velocity in the mantle: Implications for chemical and thermal structure. In S. Karato, A. Forte, R. Liebermann, G. Masters, & L. Stixrude (Eds.), *Earth's deep interior: Mineral physics and tomography from the atomic to the global Scale* (Vol. 117). AGU monograph.

McLachlan, G., & Peel, D. (2000). *Finite mixture models*. John Wiley and Sons.

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. https://doi.org/10.1016/j.pepi.2020.106542

Merrill, R., McElhinny, M., & McFadden, P. (1998). *The magnetic field of the Earth: Paleomagnetism, the core, and the deep mantle*. Academic Press.

Metman, M. C., Livermore, P. W., & Mound, J. E. (2018). The reversed and normal flux contributions to axial dipole decay for 1880–2015. *Physics of the Earth and Planetary Interiors*, 276, 106–117. https://doi.org/10.1016/j.pepi.2017.06.007

Moffatt, H. (1978). *Magnetic field generation in electrically conducting fluids*. Cambridge University Press.

Mound, J., & Davies, C. (2023). Longitudinal structure of Earth's magnetic field controlled by lower mantle heat flow. *Nature Geoscience*, 16(4), 380–385. https://doi.org/10.1038/s41561-023-01148-9

Nakagawa, T., & Tackley, P. J. (2008). Lateral variations in CMB heat flux and deep mantle seismic velocity caused by a thermal-chemical-phase boundary layer in 3D spherical convection. *Earth and Planetary Science Letters*, 271(1–4), 348–358. https://doi.org/10.1016/j.epsl.2008. 04.013

Nilsson, A., Holme, R., Korte, M., Suttie, N., & Hill, M. (2014). Reconstructing Holocene geomagnetic field variation: New methods, models and implications. *Geophysical Journal International*, 198(1), 229–248. https://doi.org/10.1093/gji/ggu120
Nilsson, A., Suttie, N., Stoner, J. S., & Muscheler, R. (2022). Recurrent ancient geomagnetic field anomalies shed light on future evolution of the south Atlantic anomaly. *Proceedings of the national academy of sciences* (Vol. 119(24), p. e2200749119). https://doi.org/10.1073/pnas.2200749119

Olson, P., & Amit, H. (2006). Changes in earth's dipole. *Naturwissenschaften*, 93(11), 519–542. https://doi.org/10.1007/s00114-006-0138-6

Olson, P., & Christensen, U. R. (2002). The time averaged magnetic field in numerical dynamos with nonuniform boundary heat flow. *Geophysical Journal International*, 151(3), 809–823. https://doi.org/10.1046/j.1365-246x.2002.01818.x

Olson, P., Christensen, U. R., & Glatzmaier, G. A. (1999). Numerical modeling of the geodynamo: Mechanisms of field generation and equilibration. *Journal of Geophysical Research*, 104(10), 10383–10404. https://doi.org/10.1029/1999jb900013

Olson, P., Landeau, M., & Reynolds, E. (2017). Dynamo tests for stratification below the core-mantle boundary. *Physics of the Earth and Planetary Interiors*, 271, 1–18. https://doi.org/10.1016/j.pepi.2017.07.003

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(1), 66–79. https://doi.org/10.1016/j.pepi.2010.02.010

Oruba, L., & Dormy, E. (2014). Predictive scaling laws for spherical rotating dynamos. *Geophysical Journal International*, 198(2), 828–847. https://doi.org/10.1093/gji/ggu159

Panovska, S., Korte, M., & Constable, C. G. (2019). One hundred thousand years of geomagnetic field evolution. *Reviews of Geophysics*, 57(4), 1289–1337. https://doi.org/10.1029/2019rg000656

Pedlosky, J. (1987). *Geophysical fluid dynamics*. Springer-Verlag.

Poletti, W., Trindade, R. I. F., Hartmann, G. A., Damiani, N., & Rech, R. M. (2016). Archeomagnetism of Jesuit missions in South Brazil (1657–1706 AD) and assessment of the South American database. *Earth and Planetary Science Letters*, 445, 36–47. https://doi.org/10.1016/j.epsl.2016.04.006

Rogers, H. F., Gillet, N., Aubert, J., Personnettaz, P., & Mandea, M. (2025). Effects of geodynamo priors and geomagnetic data on inverted core surface flows. *Physics of the Earth and Planetary Interiors*, 364, 107323. https://doi.org/10.1016/j.pepi.2025.107323

Rother, M., Korte, M., Morschhauser, A., Vervelidou, F., Matzka, J., & Stolle, C. (2021). The mag.num core field model as a parent for IGRF-13, and the recent evolution of the South Atlantic Anomaly. *Earth Planets and Space*, 73(50), 50. https://doi.org/10.1186/s40623-020-01277-0

Sahoo, S., & Sreenivasan, B. (2020). Response of Earth's magnetic field to large lower mantle heterogeneity. *Earth and Planetary Science Letters*, 549, 116507. https://doi.org/10.1016/j.epsl.2020.116507

Sanchez, S., Wicht, J., & Bärenzung, J. (2020). Predictions of the geomagnetic secular variation based on the ensemble sequential assimilation of geomagnetic field models by dynamo simulations. *Earth Planets and Space*, 72(1), 157. https://doi.org/10.1186/s40623-020-01279-y

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

Schwaiger, T., Gastine, T., & Aubert, J. (2019). Force balance in numerical geodynamo simulations: A systematic study. *Geophysical Journal International*, 219(Supplement_1), S101–S114. https://doi.org/10.1093/gji/ggz192

Shah, J., Koppers, A. A. P., Leitner, M., Leonhardt, R., Muxworthy, A. R., Heunemann, C., et al. (2016). Paleomagnetic evidence for persistence of recurrence of geomagnetic main field anomalies in the South Atlantic. *Earth and Planetary Science Letters*, 441, 113–124. https://doi.org/10.1016/j.epsl.2016.02.039

Takahashi, F., Tsunakawa, H., Matsushima, M., Mochizuki, N., & Honkura, Y. (2008). Effects of thermally heterogeneous structure in the lowermost mantle on the geomagnetic field strength. *Earth and Planetary Science Letters*, 272(3–4), 738–746. https://doi.org/10.1016/j.epsl.2008.06.017

Tarduno, J. A., Watkeys, M. K., Huffman, T. N., Cottrell, D. C., Blackman, E. G., Wendt, A., et al. (2015). Antiquity of the South Atlantic Anomaly and evidence for top-down control on the geodynamo. *Nature Communications*, 6(1), 7865. https://doi.org/10.1038/ncomms8865

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

Tema, E., Santos, Y., Trindade, R., Hartmann, G. A., Hatakeyama, T., Terra-Nova, F., et al. (2023). Archaeointensity record of weak field recurrence in Japan: New data from Late Yayoi and Kofun ceramic artefacts. *Geophysical Journal International*, 233(2), 950–963.

Terra-Nova, F., & Amit, H. (2020). Magnetic boundary layers in numerical dynamos with heterogeneous outer boundary heat flux. *Physics of the Earth and Planetary Interiors*, 309, 106589. https://doi.org/10.1016/j.pepi.2020.106589

Terra-Nova, F., & Amit, H. (2024). Regionally-triggered geomagnetic reversals. *Scientific Reports*, 14(1), 9639. https://doi.org/10.1038/s41598-024-59849-z

Terra-Nova, F., & Amit, H. (2025). Necessary core-mantle boundary heat flux patterns for large latitude of geomagnetic surface intensity minima [Dataset]. Zenodo. https://doi.org/10.5281/zenodo.14674112

Terra-Nova, F., Amit, H., & Choblet, G. (2019). Preferred locations of weak surface field in numerical dynamos with heterogeneous core-mantle boundary heat flux: Consequences for the South Atlantic Anomaly. *Geophysical Journal International*, 217(2), 1179–1199. https://doi.org/10.1093/gji/ggy519

Terra-Nova, F., Amit, H., Hartmann, G. A., & Trindade, R. I. F. (2016). Using archaeomagnetic field models to constrain the physics of the core: Robustness and preferred locations of reversed flux patches. *Geophysical Journal International*, 206(3), 1890–1913. https://doi.org/10.1093/gji/ggw248

Terra-Nova, F., Amit, H., Hartmann, G. A., Trindade, R. I. F., & Pinheiro, K. J. (2017). Relating the South Atlantic Anomaly and geomagnetic flux patches. *Physics of the Earth and Planetary Interiors*, 266, 39–53. https://doi.org/10.1016/j.pepi.2017.03.002

Trindade, R., Jaqueto, P., Terra-Nova, F., et al. (2018). Speleothem record of geomagnetic South Atlantic Anomaly recurrence. *Proceedings of the national academy of sciences* (Vol. 155(52), pp. 13198–13203).

Wardinski, I., Terra-Nova, F., & Thébault, E. (2025). Evaluation of archeo- and paleomagnetic field models and their common features. *Physics of the Earth and Planetary Interiors*, 366, 107399. https://doi.org/10.1016/j.pepi.2025.107399

Wicht, J. (2002). Inner-core conductivity in numerical dynamo simulations. *Physics of the Earth and Planetary Interiors*, 132(4), 281–302. https://doi.org/10.1016/s0031-9201(02)00078-x

Yue, Y., et al. (2024). Evolution and disappearance of the paleo-West Pacific Anomaly: Implications to the future of South Atlantic Anomaly. *Physics of the Earth and Planetary Interiors*, 353, 107–214. https://doi.org/10.1016/j.pepi.2024.107214