Paleomagnetic signatures of core-mantle interactions inferred from top-heavy thermochemical geodynamo simulations

We employ a numerical geodynamo model of a Boussinesq fluid [Willis et al., 2007, Davies et al., 2011, Matsui et al., 2016]. A spherical coordinate system $(r,\theta,\phi)$ is used to represent the domain bounded by the inner and outer boundaries, with radii $r_{i}$ and $r_{o}$, respectively, with the shell gap $h=r_{o}-r_{i}$. The whole system rotates with a constant angular velocity $\boldsymbol{\Omega}=\Omega\hat{\boldsymbol{z}}$ about the vertical $\hat{\boldsymbol{z}}$ axis, and gravity varies linearly with radius as $\boldsymbol{g}=-g_{o}(r/r_{o})\hat{\boldsymbol{r}}$ where $g_{o}$ is the gravitational acceleration at the outer boundary. The fluid is a mixture of light components dissolved in a comparatively heavy liquid (e.g., oxygen mixed in liquid iron in Earth’s outer core). The relevant physical properties of the mixture are the kinematic viscosity, $\nu$, the thermal and chemical diffusivities, $\kappa_{T}$ and $\kappa_{\xi}$, and the coefficients of thermal and chemical expansion, $\alpha_{T}$ and $\alpha_{\xi}$. The thermal diffusivity is defined as $\kappa_{T}=k_{T}/\rho_{o}c_{p}$, where $k_{T}$ is the thermal conductivity, $\rho_{o}$ is the reference density, and $c_{p}$ is the specific heat of the mixture. The fluid is electrically conducting, with magnetic permeability $\mu_{o}$ and magnetic diffusivity $\eta$.

The boundaries are assumed to be no-slip and electrically insulating at $r_{i}$ and $r_{o}$. Fixed-flux thermal boundary conditions are imposed at the inner and outer boundaries such that the total radial heat flow is equal through the inner and outer surfaces ($\boldsymbol{Q}_{T,i}=\boldsymbol{Q}_{T,o}$). The temperature gradient at the boundaries is expressed as $\boldsymbol{\nabla}T_{c}=-(\beta_{T}/r^{2})\hat{\boldsymbol{r}}$ and related to the heat flux at the ICB through Fourier’s Law as, $\boldsymbol{Q}_{T,i}=4\pi r_{i}^{2}(-k_{T}\boldsymbol{\nabla}T_{c})=4\pi k_{T}\beta_{T}\hat{\boldsymbol{r}}$. A standard setup [e.g., Kutzner and Christensen, 2002] is used to model the release of light elements from the inner core boundary, with fixed flux conditions imposed at $r_{i}$ such that $\boldsymbol{\nabla}\xi_{c,i}=-(\beta_{\xi}/r_{i}^{2})\hat{\boldsymbol{r}}$. To ensure stationary solutions, we assume that the chemical flux from the inner core is balanced by a spatially homogeneous sink ($S_{\xi}$) that maintains the global balance of lighter elements [see e.g., Kono and Roberts, 2001].

The nondimensional governing equations determining the velocity field ($\boldsymbol{u}$), the magnetic field ($\boldsymbol{B}$), temperature ($T$) and chemical composition ($\xi$) can be found in equations A15-A19 of the Supplementary Information. The non-dimensional numbers appearing in these equations are the Ekman number ($E$), the thermal and chemical flux Rayleigh numbers ($Ra_{T}$ and $Ra_{\xi}$), and the thermal, chemical, and magnetic Prandtl numbers ($Pr_{T}$, $Pr_{\xi}$ and $Pm$) defined as

The modified flux Rayleigh number used in the subsequent sections relates to the flux Rayleigh number as $\widetilde{Ra}_{T}=Ra_{T}E$ (and similarly $\widetilde{Ra}_{\xi}=Ra_{\xi}E$). To isolate the effect of CMB heat flux heterogeneity, we run simulation pairs with homogeneous thermal boundary conditions and laterally heterogeneous thermal flux at the CMB, following Mound and Davies [2017]. The parameter space now also includes the pattern and amplitude of lateral variation in the CMB heat flux. In this study, the pattern is fixed using an inferred heat flux pattern from seismic tomography following Masters et al. [1996]. The amplitude of heterogeneity is characterised as $q^{*}=(q^{T}_{max}-q^{T}_{min})/q^{T}_{avg}$ where $q^{T}_{max}$, $q^{T}_{min}$, and $q^{T}_{avg}$ are the maximum, minimum, and horizontally averaged heat flux through the CMB, respectively. The amplitude is also fixed at $q^{*}=5$, where $q^{*}>2$ has been reported to produce RILs for pure thermal convection [Mound et al., 2019]. The homogeneous models correspond to $q^{*}=0$. The average heat flux $q^{T}_{avg}$ is the same for homogeneous and heterogeneous model pairs.

We have performed $14$ new geodynamo simulations by varying $\widetilde{Ra}_{T}$, $\widetilde{Ra}_{\xi}$ and $q^{*}$ while keeping the other parameters fixed at $r_{i}/r_{o}=0.35$, $E=10^{-5}$, $Pr_{T}=1$, $Pr_{\xi}=10$ and $Pm=5$. The aspect ratio $r_{i}/r_{o}=0.35$ is representative of Earth’s present-day core configuration. The Ekman number is chosen small enough to ensure rotationally constrained dynamics, while allowing sufficient runtime to assess the typical magnetic field behaviour. We also adopt the common practice of keeping $Pr_{T}=1$, although this is considerably higher than the estimated value for the core ($Pr_{T}\sim 0.05$, see Pozzo et al. [2013]). Conversely, our choice of $Pr_{\xi}=10$ is considerably lower than the core estimate ($Pr_{\xi}\sim 100$, Pozzo et al. [2013]), and is selected solely for computational benefit. The Prandtl number ratio $Pr_{\xi}/Pr_{T}$, which controls the scale separation between the two scalar fields, is therefore much lower in our simulations compared to Earth’s core.

Of the $14$ geodynamo simulations, $7$ are run with homogeneous CMB heat flux conditions ($q^{*}=0$). Four of these simulations with $\widetilde{Ra}_{T}\in\{90,550,1200,4000\}$ are run at $\widetilde{Ra}_{\xi}=10^{4}$ and three with $\widetilde{Ra}_{T}\in\{90,550,1200\}$ are run at $\widetilde{Ra}_{\xi}=10^{5}$. At $E=10^{-5}$, onset of pure thermal convection with $Pr_{T}=1$ and pure chemical convection with $Pr_{\xi}=10$ occurs at $\widetilde{Ra}_{T}^{c}=24.7$ and $\widetilde{Ra}_{\xi}^{c}=34.1$, respectively [Naskar et al., 2025]. The range of Rayleigh numbers are chosen from our previous parameter survey of non-magnetic double-diffusive convection simulations at $E=10^{-5}$ [Naskar et al., 2026], which identifies the dynamical regime where vigorous turbulence coexists with strong rotational constraint. Along with $7$ homogeneous cases, $7$ more simulations are run for the same ($\widetilde{Ra}_{T}$, $\widetilde{Ra}_{\xi}$) pairs, with a tomographic pattern of heat flux at the CMB with $q^{*}=5$. In the subsequent sections, the simulations with homogeneous CMB heat flux conditions (i.e. $q^{*}=0$), and those with a tomographic pattern of CMB heat flux with $q^{*}=5$ are simply referred to as homogeneous and heterogeneous simulations, respectively. The simulation parameters are listed in B.

Hyperdiffusion has been used to enable long simulation runs [see e.g., Aubert et al., 2017], where enhanced diffusivity is applied to scales smaller than a certain ”cut-off” length-scale, represented by a spherical harmonic degree $\ell_{h}=128$. The diffusion operator applied to the momentum, composition, and magnetic field remains as $\nabla^{2}$ for $\ell<\ell_{h}$, and is replaced by $q_{h}^{\ell-\ell_{h}}\nabla^{2}$ for $\ell\geq\ell_{h}$. The value of the parameter $q_{h}=1.0325$ is chosen to ensure a smooth increase of hyperdiffusion with wavenumber. After some initial transient, all simulations have been run for at least $0.2$ magnetic diffusion times, which has been reported to be long enough to determine the signature of CMB heterogeneity in the time-average morphology of the magnetic field [Davies and Constable, 2014, Mound et al., 2015].

We use the following dimensionless output diagnostics to analyse the simulations. The global root mean square flow velocity is represented by the magnetic Reynolds number

where $U$ is the characteristic velocity, $V_{s}$ is the volume of the outer core, and the kinetic energy is defined as

The magnetic Reynolds number is related to the Reynolds number through $Re=Uh/\nu=Rm/Pm$. Also, the total buoyant power ($P_{tot}$) is defined as,

where $u_{r}$ is the radial velocity component. Additionally, the Elsasser number $\Lambda$ measures the dimensionless magnetic energy $E_{\textrm{mag}}$ and is given by

where

Finally, dipolarity of the magnetic field at the outer boundary is calculated by

Earth’s core is expected to exhibit a primary quasi-geostrophic balance between Coriolis and Pressure forces, with the residual Coriolis force balanced by Lorentz and buoyancy forces in the secondary balance, constituting a QG-MAC balance of forces governing the dynamics, while inertial and viscous effects remain small [Aurnou and King, 2017, Calkins, 2018, Davidson, 2013, Yadav et al., 2016, Aubert et al., 2017, Aubert, 2019, Schwaiger et al., 2019, 2020]. We estimate the force balance in our simulations from the scale-dependent forces along with the volume-averaged ratio of various force terms in the fluctuating momentum equation (equation A.2), following Naskar et al. [2025]. In Naskar et al. [2025], the forces and curled forces are partitioned into mean (i.e. azimuthal average) and corresponding fluctuating parts (see also e.g.Calkins et al. [2021], Nicoski et al. [2024]), and we only consider the fluctuating part of the forces here (indicated by a prime). All the forces are integrated over the bulk fluid, excluding regions within a radial thickness corresponding to $10$ times the Ekman-layer depths adjacent to the upper and lower boundaries (defined using the linear intersection method). The metrics to verify the dynamical balances are based on the ratio of the magnitudes of the various forces, and their curls [Naskar et al., 2025, Clarke et al., 2026]. The inertial contribution to the primary balance is quantified by the ratio of inertia to Coriolis forces $F_{I/C}$, and the ratios of curled Lorentz (and Archimedean) to Coriolis forces $CF_{M/C}$ (and $CF_{A/C}$) have been used to check the secondary balance. Here, the curled force representation ensures reliable estimates of the secondary forces by eliminating the dynamically irrelevant gradient parts of the forces that are balanced by the pressure gradient[Teed and Dormy, 2023, 2025]. The metrics are defined as

where the total buoyancy is $\boldsymbol{A^{\prime}}=\boldsymbol{A^{\prime}_{T}}+\boldsymbol{A^{\prime}_{\xi}}$. Additionally, we use another metric to estimate the relative importance of two buoyancy forces,

We compute existing criteria to quantify the spatial structure and temporal variability of the simulated magnetic fields and compare them to geomagnetic and paleomagnetic observations across various time scales. Our aim is to assess differences between homogeneous and heterogeneous simulations, and between thermally dominated and chemically dominated simulations. We further seek to evaluate whether such differences are reflected in observed features of Earth’s magnetic field. We consider the gufm1 field model [Jackson et al., 2000] spanning the past 400 years, the CALS10k.2 model [Constable et al., 2016] of the last 10,000 years, and the GGF100k model [Panovska et al., 2018] spanning the last 100,000 years. This range of timescales is both numerically achievable with the present set of simulation parameters and sufficient to investigate the signature of CMB heterogeneity in the non-dipole time-averaged field (TAF) [Davies and Constable, 2014]. To compare with the field models, the magnetic diffusive timescale in the simulations are rescaled assuming a magnetic diffusivity of $\eta=1.6\ \text{m}^{2}/\text{s}$, which corresponds to a magnetic diffusion time of roughly 100,000 years. We compare simulated and observed fields based on:

• 1.

  Temporal field variability on decadal timescales using the master secular variation time-scale, $\tau_{SV}$, following Lhuillier et al. [2011];

• 2.

  Morphological semblance on the timescale of centuries using the four compliance criteria proposed by Christensen et al. [2010];

• 3.

  $10-100$ kyr variability using the paleo-secular variation index ($P_{i}$), with $\ell_{max}=10$ [Panovska and Constable, 2017];

• 4.

  Statistical properties of field morphology and temporal variability over $>100$ kyr timescales using the five $Q_{PM}$ criteria proposed by Sprain et al. [2019].

Precise definitions of these criteria are given in section A.3 of the SI, though they are identical to those used in the original publications.

Thermal RILs that form under LLVPs or light elements that accumulate preferentially near the poles may imprint regional signatures on the magnetic field. If sufficiently strong and persistent, such signatures would be expressed as departures from a geocentric axial dipole field in the long-term behaviour of the Earth. Therefore, we look for measures of such behaviour that are both potentially observable and quantitatively distinguishable between homogeneous and heterogeneous (or thermally/chemically dominated) cases. In particular, we consider the following:

• 1.

  The dipole-dominance and longitudinal variation of the TAF at Earth’s surface quantified using the ratios $(AD/NAD)_{TAF}$ and $(ZND/ND)_{TAF}$, calculated from the Lowes power spectra for the magnetic energy [Lowes, 1974], following Biggin et al. [2026];

• 2.

  The polar minima in the TAF quantified by the time-averaged radial field magnitude at the normalised by its maximum value in the northern hemisphere;

• 3.

  Deviation from GAD measured by the ratios of specific Gauss coefficients relative to the axial dipole (e.g., $g^{1}_{1}/g^{0}_{1}$ and $g^{0}_{3}/g^{0}_{1}$).

• 4.

  Eccentricity of the dipole, calculated as an offset from the centre of Earth in the equatorial plane ($r_{dip}$) [James and Winch, 1967, Olson and Deguen, 2012].

• 5.

  Longitudinal structures near the equator [Mound and Davies, 2023] and their variability, quantified using the median ($S_{|30|median}$) and inter-quartile range ($S_{|30|iqr}$) of dispersion of virtual geomagnetic poles (VGPs) at low latitudes following Biggin et al. [2026].

• 6.

  Longitudinal variation of the time- and latitudinally-averaged paleo-secular variation index $P_{i}$ following [Mason et al., 2024].

• 7.

  The magnitude and spatial variability of the inclination anomaly ($\Delta I$), on Earth’s surface. The longitudinal variability of inclination anomaly is calculated from the time-averaged $\Delta I$ map at the surface. At each latitude, we define the range of $\Delta I$ as the difference between the maximum and minimum $\Delta I$ in the longitudinal distribution.

All of the metrics outlined in this section will not be discussed in detail below; however, they are reported in full in B to assist further comparison with previous work.

We begin by examining several global diagnostic quantities to characterise the dynamo runs. Generally, the simulations produce dipole-dominated fields with dipole fractions $f_{\textrm{dip}}\geq 0.48$ (Figure 1), except for one heterogeneous case ($q^{*}=5$, $\widetilde{Ra}_{T}=90$ and $\widetilde{Ra}_{\xi}=10^{5}$) with $f_{\textrm{dip}}=0.27$ that exhibits multipolar reversing behaviour. Most of our simulations access Earth-like values of the magnetic Reynolds number $Rm$, which is estimated to vary between $\sim 500$ and several thousand in Earth’s core depending on assumed values of the core’s electrical conductivity and the root mean square flow speed [Davies et al., 2015]. In our suite, the heterogeneous simulations (square symbols, Figure 1) systematically yield smaller $f_{\textrm{dip}}$ values compared to their homogeneous counterpart, despite having comparable $Rm$. The magnetic-to-kinetic energy ratio, $M_{ratio}=E_{\textrm{mag}}/E_{\textrm{kin}}$, which is large in Earth’s core, ranges between $2.9\text{---}39.4$ in our simulations (see B) with no clear trend as the control parameters are varied. With the exception of three simulations at $\widetilde{Ra}_{T}=90$, the suite of simulations sits within a range of $f_{\textrm{dip}}$ and $Rm$, that has been found to produce magnetic fields with good compliance with the observed characteristics of Earth’s magnetic field [Nakagawa and Davies, 2022].

The dynamical regime of our simulations can be further assessed using scale-dependent force spectra as well as volume-averaged force ratios [Naskar et al., 2025]. In the heterogeneous simulations (e.g., Figure A1b), thermal buoyancy is systematically enhanced at spherical harmonic degree $\ell=2$, the dominant scale of the imposed thermal heterogeneity, relative to their homogeneous counterparts (e.g., Figure A1a). Furthermore, simulations with high $\widetilde{Ra}_{\xi}$ (e.g., Figure A1c) show stronger chemical buoyancy with the crossover between dominant buoyancy and Lorentz force occurring at a smaller spatial scale (i.e. higher $\ell$) than corresponding low $\widetilde{Ra}_{\xi}$ cases (e.g., Figure A1b). Nevertheless, across all cases, the simulations exhibit the characteristic QG-MAC balance [Aubert et al., 2017, Schwaiger et al., 2019, Nakagawa and Davies, 2022], regardless of the heat flux boundary condition or the dominant buoyancy source.

The dependence of the global force balance on the control parameters can be quantified with volume-averaged force measures [Naskar et al., 2025, Teed and Dormy, 2025, Clarke et al., 2026]. The integrated ratios of curled Lorentz to curled Coriolis forces remain near unity ($CF_{M/C}\sim 1$), indicating strong field solutions for all cases [see also Teed and Dormy, 2025]. Furthermore, the curled total buoyancy to curled Coriolis forces ($CF_{A/C}$) fall within $\mathcal{O}(0.1-1)$ across the parameter range examined (B), consistent with a secondary MAC balance. The ratio of inertia to Coriolis forces ($F_{I/C}$) increases with convective driving in the simulations but remains $\mathcal{O}(0.01-0.1)$, confirming that inertia plays a subdominant dynamical role. The ratio of thermal to chemical buoyancy force ($A_{T/\xi}$) varies by nearly three orders of magnitude over the parameter space explored here (B), spanning regimes from thermally dominated convection (e.g., $A_{T/\xi}=41.5$ at $\widetilde{Ra}_{T}=4000$, $\widetilde{Ra}_{\xi}=10^{4}$ and $q^{*}=0$) to chemically dominated convection (e.g., $A_{T/\xi}=0.09$, at $\widetilde{Ra}_{T}=90$, $\widetilde{Ra}_{\xi}=10^{5}$ and $q^{*}=0$).

We now consider metrics describing the general spatial morphology and time variability of the simulated magnetic fields, with particular attention to whether systematic differences arise between homogeneous and heterogeneous simulations, or between thermally and chemically dominated regimes. These metrics were developed taking into account constraints from different observational models and thus provide complementary expectations for field behaviour over a range of timescales (recall Section 2.5).

The SV timescale, $\tau_{SV}$, applies to relatively short-term temporal variability and decreases with increasing thermal or chemical forcing in our simulations. The homogeneous simulations have $\tau_{SV}$ values that are somewhat (approximately $5-15$%) smaller than their heterogeneous counterparts; however, all simulations follow an approximately inverse scaling relation between $\tau_{SV}$ and the magnetic Reynolds number (Figure A2). An empirical fit to the data gives coefficients closely matching the scaling predictions of Christensen et al. [2012] (see their equation $18$).

Comparison to the spatial morphology of Earth’s recent magnetic field is covered by the $\chi^{2}_{total}$ compliance computed following Christensen et al. [2010]. Time-averaged values for the homogeneous simulations generally have a lower average $\chi^{2}_{total}$ than their heterogeneous counterparts (B). This contrast is mainly due to the propensity of heterogeneous simulations to have patches of concentrated flux that are stronger than both the homogeneous simulations and those inferred for modern Earth. However, since the flux concentration factor (see section A.3) is not well constrained for the geomagnetic field on millennial and longer timescales [Panovska et al., 2019, Wardinski et al., 2025], this does not necessarily imply a shortcoming of the heterogeneous simulations. Indeed, most (12 of 14) simulations have periods of time for which they are in good or excellent morphological agreement with the modern geomagnetic field (i.e., $\chi^{2}_{total,min}<4$). Other than the formation of strong flux patches in heterogeneous simulations, the $\tau_{SV}$ and $\chi^{2}_{total}$ metrics, derived based on models of Earth’s modern geomagnetic field, indicate relatively small differences between homogeneous and heterogeneous simulations. Furthermore, the simulations can achieve excellent compliance ($\chi^{2}_{total,min}\leq 2$) regardless of whether the driving mechanism is thermal or compositional, similar to the findings of Tassin et al. [2021].

The long-term variability of the magnetic field can be characterised by the median of paleo-secular variation index [Panovska and Constable, 2017], which decreases with increasing dipolarity of the simulated fields (Figure 2b). Conversely, the $\Delta Q_{PM}$ misfits (Figure 2a) generally decrease with decreasing $f_{\textrm{dip}}$ reaching values of $\Delta Q_{PM}\sim 6-8$, around $f_{\textrm{dip}}\sim 0.5-0.6$. The multipolar reversing case departs from this trend, exhibiting a much lower $f_{\textrm{dip}}$. A similar dependence of misfit with $f_{\textrm{dip}}$ was identified by Meduri et al. [2020], who reported a minimum $\Delta Q_{PM}\sim 5$ near $f_{\textrm{dip}}\sim 0.5-0.6$. However, given that our simulations are considerably shorter ($\sim 20\text{---}100$ kyr) than the paleofield observations on which the $Q_{PM}$ criteria are based, we don’t place emphasis on absolute $\Delta Q_{PM}$ values and instead focus on inter-simulation comparisons. Since $f_{\textrm{dip}}$ for the homogeneous simulations is systematically higher than their heterogeneous counterparts, the homogeneous simulations have systematically lower values for their median $P_{i}$ and systematically larger values of $\Delta Q_{PM}$. No systematic dependence of these metrics on the relative thermal versus chemical forcing is evident.

The criteria so far considered are based on global measures of the magnetic field’s spatial morphology or temporal variability, and thus are not well suited to detect the potential regional signatures of CMB heat flux heterogeneity or preferential LEA in polar regions. To illustrate the typical differences between our simulations, we focus on three simulations: one homogeneous ($q^{*}=0,\ \widetilde{Ra}_{\xi}=10^{4}$) and two heterogeneous cases ($q^{*}=5,\ \widetilde{Ra}_{\xi}=10^{4}$ and $\widetilde{Ra}_{\xi}=10^{5}$), all at fixed $\widetilde{Ra}_{T}=550$. The time-averaged field of the $\widetilde{Ra}_{\xi}=10^{4}$ heterogeneous model (Figure 3b) exhibits significant longitudinal variation compared to the corresponding homogeneous model (Figure 3a). This propensity for heterogeneous boundary forcing to induce longitudinal field structure holds across all our simulations and agrees with previous results for purely thermal simulations [Olson and Christensen, 2002, Olson et al., 2018a, Mound and Davies, 2023, Biggin et al., 2026].

The most prominent signature of the heterogeneous simulations is the presence of high-latitude magnetic flux patch(es) in the TAF that are absent in the homogeneous simulations; although the heterogeneous simulations also exhibit longitudinal structures near the equator, as found by Mound and Davies [2023]. The precise structure and position of the longitudinal features in the TAF vary with the simulation control parameters; however, their importance can be estimated by the ratio of axial dipole to non-axial dipole power in the time-averaged field, $(AD/NAD)_{TAF}$. If we discount the case with $f_{\textrm{dip}}=0.27$, the homogeneous simulations have consistently larger values of this ratio than their heterogeneous counterparts (by a factor of $\sim 5-20$). Additionally, all of our simulations exhibit a polar minimum in the radial magnetic field ($B_{r}$), which tends to deepen with increasing $\widetilde{Ra}_{\xi}$ and simulations dominated by chemical buoyancy tend to have a more pronounced polar minimum than thermally dominated simulations at comparable total buoyant power (Figure A3). Also, heterogeneous simulations have a stronger polar minimum than their homogeneous counterparts.

Long-term observational models (e.g., GGF100k, CALS10k.2) are constrained for low spatial resolutions; we therefore compare the simulations and field models in terms of the low-degree ($\ell\leq 4$) time-averaged Gauss coefficients (normalised by $g^{0}_{1}$). The observational models are consistent with a zero mean for the non-dipole coefficients considered, although GGF100k displays more variability (Figure 4). Heterogeneous simulations tend to have non‑zero mean values, whereas the values for homogeneous cases approach zero (e.g., Figure 4a). That is, CMB heterogeneity tends to induce higher dipole tilt than the homogeneous cases, although a zero‑tilt scenario always lies within one standard deviation of the simulation average. The axial octupole component ($g^{0}_{3}/g^{0}_{1}$) has a larger magnitude in heterogeneous simulations than its homogeneous counterparts. Simulations with higher chemical buoyancy forcing (i.e., lower $A_{T/\xi}$) tend to produce smaller values of $g^{0}_{3}/g^{0}_{1}$ than simulations with comparable buoyant power (Figure 4b), and thus are more likely to overlap with the observational ranges.

Given the tendency of heterogeneous simulations to induce dipole tilt, we also examine whether the dipole remains concentric with the rotation axis. The dipole eccentricity ($r_{dip}$) is generally higher for the heterogeneous dynamos than in their homogeneous counterpart (e.g., Figure 5) and tends to be offset towards the Pacific hemisphere. However, the mean eccentricity of the dipole across all simulations remains within the observational range, defined by the standard deviation of $r_{dip}$ in the field models (Figure A4), and no systematic relationship is found between the dipole eccentricity and the thermal‑to‑chemical buoyancy ratio.

Among the quantities we examined for potential observable signatures of longitudinal field variation, the inclination anomaly ($\Delta I$) provides clear and systematic behaviour. The time-averaged inclination anomaly maps of our three example cases (Figure 3) are compared with the corresponding map from the GGF100k field model (Figure 6). Uncertainties in field models arising from heterogeneous spatio-temporal sampling of the geomagnetic field together with limitations in the simulation’s capacity to access the physical conditions of Earth’s core mean that we do not expect simulated fields to reproduce the specific morphology of GGF100k. However, all homogeneous simulations (e.g., Figure 6a) fail to reproduce both the magnitude and the longitudinal variability evident in the observational model (Figure 6d). In contrast, heterogeneous simulations (e.g., 6b,c) exhibit substantial longitudinal variability and a magnitude of $\Delta I$ values that more closely resemble GGF100k. Quantifying the variability by the equatorial range of inclination anomaly (Figure 7), most homogeneous cases falls below the observational ranges, while heterogeneous cases readily reproduce anomaly ranges inferred from observations.