# Dynamo constraints on the long-term evolution of Earth's magnetic field strength

Christopher J. Davies, Richard K. Bono, Domenico G. Meduri, Julien Aubert, Samuel Greenwood, Andrew J. Biggin

► To cite this version:

Christopher J. Davies, Richard K. Bono, Domenico G. Meduri, Julien Aubert, Samuel Greenwood, et al. Dynamo constraints on the long-term evolution of Earth's magnetic field strength. Geophysical Journal International, 2021, pp.2355-2377. 10.1093/gji/gab342. insu-03590027

HAL Id: insu-03590027

https://insu.hal.science/insu-03590027

Submitted on 26 Feb 2022

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
<figure><img src="image_2.png" /></figure>

# Dynamo constraints on the long-term evolution of Earth's magnetic field strength

Christopher J. Davies<fnref n="1" />,<fnref n="1" /> Richard K. Bono,<fnref n="2" /> Domenico G. Meduri,<fnref n="2" /> Julien Aubert,<fnref n="2" />,<fnref n="3" />
Samuel Greenwood<fnref n="1" /> and Andrew J. Biggin<fnref n="2" />

Accepted 2021 August 17. Received 2021 August 5; in original form 2021 May 12

## SUMMARY

Elucidating the processes in the liquid core that have produced observed palaeointensity changes over the last 3.5 Gyr is crucial for understanding the dynamics and long-term evolution of Earth's deep interior. We combine numerical geodynamo simulations with theoretical scaling laws to investigate the variation of Earth's magnetic field strength over geological time. Our approach follows the study of Aubert *et al.*, adapted to include recent advances in numerical simulations, mineral physics and palaeomagnetism. We first compare the field strength within the dynamo region and on the core–mantle boundary (CMB) between a suite of 314 dynamo simulations and two power-based theoretical scaling laws. The scaling laws are both based on a Quasi-Geostropic (QG) force balance at leading order and a Magnetic, Archimedian, and Coriolis (MAC) balance at first order and differ in treating the characteristic length scale of the convection as fixed (QG-MAC-fixed) or determined as part of the solution (QG-MAC-free). When the data set is filtered to retain only simulations with magnetic to kinetic energy ratios greater than at least two we find that the internal field together with the root-mean-square and dipole CMB fields exhibit power-law behaviour that is compatible with both scalings within uncertainties arising from different heating modes and boundary conditions. However, while the extrapolated intensity based on the QG-MAC-free scaling matches Earth's modern CMB field, the QG-MAC-fixed prediction shoots too high and also significantly overestimates palaeointensities over the last 3.5 Gyr. We combine the QG-MAC-free scaling with outputs from 275 realizations of core–mantle thermal evolution to construct synthetic true dipole moment (TDM) curves spanning the last 3.5 Gyr. Best-fitting TDMs reproduce binned PINT data during the Bruhnes and before inner core nucleation (ICN) within observational uncertainties, but PINT does not contain the predicted strong increase and subsequent high TDMs during the early stages of inner core growth. The best-fitting models are obtained for a present-day CMB heat flow of 11–16 TW, increasing to 17–22 TW at 4 GPa, and predict a minimum TDM at ICN.

**Key words:** Dynamo: theories and simulations; Magnetic field variations through time; Palaeointensity; Palaeomagnetic secular.

## 1 INTRODUCTION

Earth has sustained a global magnetic field over most of its history. Databases of palaeointensity estimates indicate no hiatuses in the geodynamo back to 3.55 GPa (Biggin et al. 2008; Tarduno et al. 2010; Tauxe & Yamazaki 2015; Biggin et al. 2015; Bono et al. 2019), while records of a field extending back to 4.2 GPa (Tarduno et al. 2015) are currently under debate (Tang et al. 2019; Tarduno et al. 2020). These observations provide a unique probe of otherwise unobservable processes in the liquid iron core where the field is generated by a hydromagnetic dynamo. The dynamo draws its power from slow cooling due to heat extraction by the overlying mantle and so palaeointensity determinations also provide information on the nature and evolution of mantle convection (e.g. Nimmo et al. 2004; Driscoll & Bercovici

<footnote n="1">¹School of Earth and Environment, University of Leeds, Leeds LS29JT, UK. E-mail: c.davies@leeds.ac.uk</footnote>

<footnote n="2">²Geomagnetism Laboratory, Department of Earth, Ocean and Ecological Sciences, University of Liverpool, Liverpool L697ZE, UK</footnote>

<footnote n="3">³Institut de Physique du Globe de Paris, Sorbonne Paris Cité, Université Paris-Diderot, CNRS, 1 rue Jussieu, F-75005 Paris, France</footnote>
2014; O'Rourke et al. 2017). Cooling of the liquid core leads to freezing at Earth's centre and the growth of the solid inner core, which provides additional power to the dynamo through release of latent heat and gravitational energy (e.g. Gubbins et al. 2004; Nimmo 2015). By linking changes in the available power, which clearly identify inner core formation (Nimmo 2015; Davies 2015; Labrosse 2015), to variations in the observable field recent studies have attempted to date inner core formation using the palaeomagnetic record (Biggin et al. 2015; Bono et al. 2019). However, this task is hampered due to uncertainties regarding the observable expression of inner core formation (Driscoll 2016; Landeau et al. 2017). In this paper, we consider the relationship between palaeointensities and core dynamics using numerical dynamo simulations.

Detailed knowledge of geomagnetic field strength variations over geological time is hampered by the uneven spatial and temporal sampling. Spatial variations are usually treated by considering the virtual dipole moment (VDM), which normalizes the expected variation of Earth's field strength that would be produced from a dipole field. Temporal sampling is hindered because ideal magnetic recorders are rare and the laboratory efforts to recover them often end in failure, so developing a global VDM database comprising entries of approximately homogeneous fidelity is a significant challenge. The PINT database (Biggin et al. 2009, 2015) represents a community effort to develop a data set of palaeointensity observations spanning 50 ka to 3.5 Ga, compiling studies over the last 70 yr. Here, we will use an extension of the PINT database (described below) with field strength estimates extending back to ~4 Ga.

Linking palaeointensity observations to the dynamo process requires numerical simulations. These simulations produce dipole-dominated fields and spontaneous reversals and have captured large-scale features of the historical geomagnetic field (Christensen et al. 2010) and the pattern of recent secular variation (e.g. Aubert et al. 2013; Mound et al. 2015). Simulations have also reproduced some features of the Holocene field (Davies & Constable 2014); however, semblance to the palaeomagnetic field over the last 10 Myr appears harder to achieve (Sprain et al. 2019) and is sensitive to the dipole-dominance of the field and the driving mode of convection (Meduri et al. 2021). Simulations typically only span $O(1)$ Myr (Davies & Constable 2014; Driscoll 2016) and can only reach Gyr timescales if very low rotation rates are employed (Wicht & Meduri 2016). In particular, within a single simulation it is impractical to explicitly account for effects arising from slow changes due to growth of the inner core or evolution of buoyancy sources (Anufriev et al. 2005; Davies & Gubbins 2011; Landeau et al. 2017). To apply simulation results over geological time therefore requires a model of long-term core thermal evolution, which is here called a 'thermal history' model.

Another important limitation of the simulations is that they cannot be run with certain parameter values that characterize the properties of Earth's core, in particular the viscous and thermal diffusion coefficients (Jones 2015), though significant recent progress has been made by following a distinguished path in parameter space towards core conditions (Aubert et al. 2017; Aubert 2019). In terms of dimensionless parameters, the Ekman number $E$, the ratio of viscous and Coriolis effects, and the magnetic Prandtl number $Pm$, the ratio of viscous and magnetic diffusivities, are too high while the Rayleigh number $Ra$, measuring the vigour of convection is usually too low. The general approach for using simulation outputs to infer behaviour in Earth's core has been through scaling analysis, where theoretical balances of terms in the governing equations are tested against large suites of simulations (e.g. Christensen & Aubert 2006; Christensen 2010). If a given theoretical scaling collapses the simulation data it gives confidence for using the scaling to extrapolate from conditions in the simulations to those in the core.

A major step forward in using dynamo simulations to model long-term palaeointensity variations was provided by Aubert et al. (2009). They showed that the root-mean-square (RMS) internal field strength in a suite of 43 dynamo simulations was consistent with a theoretical scaling based on the power density $p_A$ provided by buoyancy to drive core convection (Christensen & Aubert 2006) and adopted another empirical scaling to convert this to a dipole field strength at the core surface. They then calculated the true dipole moment (TDM) from two thermal history models, which output $p_A$ over the past 4.5 Gyr given the core-mantle boundary (CMB) heat flow $Q_{cmb}$ and a set of properties that characterize the core material. They found that variations in the predicted and observed field strength were compatible over the whole time period with little long-term change due to the weak dependence of field strength on $p_A$. They also showed that the sharpest change in field strength should occur following inner core nucleation (ICN), but questioned whether this would be observable in the palaeomagnetic data.

In this paper, we revisit the analysis of Aubert et al. (2009), incorporating three important developments from the decade following their study. First, we make use of a much larger suite of simulations that access increasingly realistic physical conditions. Second, we account for the high thermal conductivity $k$ of iron alloys that has recently been obtained by several ab initio studies conducted at core conditions (de Koker et al. 2012; Pozzo et al. 2012, 2013; Gomi et al. 2013; Zhang et al. 2020) and inferred from some (Ohta et al. 2016; Inoue et al. 2020), but not all (Konôpková et al. 2016), experimental works. Thermal history models with high $k$ predict much faster cooling rates and a younger inner core than those with low $k$ (Nimmo 2015; Davies et al. 2015; Labrosse 2015), which influences the predicted field strength as we will show. Third, we use new palaeomagnetic data compilations that now extend back to ~4.2 Ga with improved temporal coverage, particularly during the Archean/Hadean (e.g. Tarduno et al. 2015; Herrero-Bervera et al. 2016; Tarduno et al. 2020), Proterozoic (e.g. Kulakov et al. 2013; Di Chiara et al. 2017; Sprain et al. 2018; Kodama et al. 2019) and Palaeozoic (e.g. Usui & Tian 2017; Hawkins et al. 2019; Veselovskiy et al. 2019).

The objective of this paper is to test whether magnetic field strength predictions from scaling laws can reproduce Earth's modern and palaeofield strength. Our analysis follows the general approach of Aubert et al. (2009), but also differs on three main points. First, we directly compare the dipole CMB field strength and RMS CMB field strength to theoretical predictions as well as the RMS internal field. Second, we consider two plausible theoretical scaling relations for the magnetic field strength based on the theory of Starchenko & Jones (2002)
and Davidson (2013). Both scalings assume a Quasi-Geostrophic (QG) balance of terms in the Navier-Stokes equation at leading order and a second-order balance between Magnetic, Archimedian (buoyancy) and Coriolis (MAC) forces and have hence been named QG-MAC balances (Aubert et al. 2017; Schwaiger et al. 2019); the difference arises in the treatment of the characteristic length scale in the MAC balance. QG-MAC scaling laws are supported by recent high-resolution dynamo simulations (Aubert et al. 2017; Schaeffer et al. 2017; Sheyko et al. 2018; Schwaiger et al. 2019) and match Earth's modern RMS field strength when evaluated at core conditions (Aubert et al. 2017). By comparing predictions from both scalings to geomagnetic and palaeomagnetic data we hope to distinguish the relevant length scale in the QG-MAC balance, which has not yet been fully constrained by simulations (Aubert 2019). We test these scalings against data from 314 simulations and compare the predictions for the internal, CMB and CMB dipole fields against present-day geomagnetic observations before applying them to the palaeofield. Third, we use 275 realizations of core thermal history with high conductivity that span uncertainties in the key parameters (to be defined precisely below).

The paper is organized as follows. In Section 2, we outline two theoretical scaling laws that determine magnetic field strength in terms of the available convective power. Here, we also describe the simulations that are used to test these scaling laws and the thermal history models that are used to apply the scaling results to Earth's palaeofield. In Section 3.1, we compare the scaling law predictions for internal and CMB field strength to the modern geomagnetic field and to empirically derived fits to the simulation data, using various methods to filter the suite of simulations. In Section 3.2, we use both scaling laws to produce synthetic palaeointensity time-series from the 275 core thermal history models. In Section 4, we discuss the implications of our results for the dynamics and evolution of Earth's core.

## 2 METHODS

### 2.1 Theoretical field strength predictions

Much of the theory presented in this section has appeared in various forms in previous work and so only a brief description is given. For more detailed treatment, the reader is referred to King & Buffett (2013), Davidson (2013), Jones (2015) and Aubert et al. (2017). Consider an electrically conducting Boussinesq fluid characterized by its density $ρ$, viscosity $ν$, thermal conductivity $k$, specific heat capacity $C_p$ and magnetic diffusivity $η$. Here and in Section 2.2, these properties will be taken as constants, but in Section 2.3 they will vary with radius $r$. The fluid is confined to a spherical shell of thickness $L = r_o - r_i$ rotating about the vertical $\hat{z}$ direction with frequency $\Omega$. Here, $r_o$ and $r_i$ are the outer and inner boundaries that may be identified with the CMB and inner core boundary (ICB), respectively. For the theoretical considerations conditions on both boundaries are assumed to be spatially uniform.

The goal is to establish the balance of physical effects that determine the characteristic field strength within the dynamo region and on the outer boundary. There are two approaches, based on local and global balances. Since we are interested in both the internal and CMB field it is necessary to use local balances, but useful information can also be gained from the global balance. The Navier-Stokes equation for the local force balance can be written in dimensional form as

$$ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} + 2\Omega \hat{\mathbf{z}} \times \mathbf{u} = -\nabla \bar{P} + \frac{gC' \mathbf{r}}{\rho} + \frac{(\nabla \times \mathbf{B}) \times \mathbf{B}}{\rho \mu_0} + \nu \nabla^2 \mathbf{u}. \quad (1) $$

Here, $\mathbf{u}$ is the fluid velocity, $\mathbf{r}$ is the position vector, $\mathbf{B}$ the magnetic field vector, $C'$ is a density anomaly about a state of rest, $\bar{P}$ the modified pressure (including the centrifugal force), $g$ the acceleration due to gravity at $r_o$ and $\mu_0$ the permeability of free space. The primary balance at leading order is geostrophic in high-resolution simulations (Schaeffer et al. 2017; Aubert 2019; Schwaiger et al. 2019), and possibly in Earth's core (Aurnou & King 2017), and so the vorticity equation, obtained from the curl of eq. (1) is used in the subsequent analysis. Ignoring viscous and inertial effects, which are thought to be very small in the Earth (Davidson 2013; Jones 2015) and have been shown to be small in high-resolution simulations (e.g. Schaeffer et al. 2017; Sheyko et al. 2018; Aubert 2019; Schwaiger et al. 2019) gives a vorticity balance between Magnetic, buoyancy (Archimedian) and ageostrophic Coriolis effects, the MAC balance:

$$ 2\Omega \frac{\partial \mathbf{u}}{\partial z} \sim \frac{g \nabla \times C' \mathbf{r}}{\rho} \sim \frac{\nabla \times [(\nabla \times \mathbf{B}) \times \mathbf{B}]}{\rho \mu_0}. \quad (2) $$

Note that the first term includes only the part of the Coriolis effect that is not balanced by the pressure gradient.

To estimate individual terms we define the characteristic velocity $U$, magnetic field strength $B$ and density anomaly $C$. The theory of Davidson (2013) defines three length scales: $l_u$, the dominant scale of flow structures in the plane perpendicular to the rotation axis; the flow scale parallel to the rotation axis, which is here taken to be $L$; and $l_{Bmin}$, the scale at which magnetic energy is dissipated. With these definitions the terms in eq. (2) can be estimated as

$$ \frac{\Omega U}{L} \sim \frac{gC}{\rho l_u} \sim \frac{B^2}{\rho \mu_0 l_u^2}, \quad (3) $$

where vorticity has been assumed to scale as $U/l_u$.

Eq. (3) is complemented by considering the global kinetic and magnetic energy balance, which can be obtained by taking the scalar product of eq. (1) with $\mathbf{u}$, integrating over the shell volume $V_{oc}$, and using the magnetic energy balance to equate the work done by the Lorentz force to the ohmic dissipation. Averaging over convective timescales (denoted by an overbar) yields an exact balance between buoyant power
$P_A$, ohmic dissipation $D_O$ and viscous dissipation $D_V$: $P_A = D_O + D_V$, or

$$g \int \overline{u_r C'} dV_{oc} = \frac{\eta}{\mu_0} \int (\overline{\nabla \times B})^2 dV_{oc} + \rho \nu \int (\overline{\nabla \times u})^2 dV_{oc}, \quad (4)$$

where $u_r$ is the radial velocity. Assuming Ohmic dissipation dominates, as expected in the core (e.g. Jones 2015; Aubert et al. 2017), the scaling estimate of eq. (4) is

$$g \overline{u_r C'} \sim \frac{\eta \overline{B^2}}{\mu_0 \ell_{Bmin}^2}. \quad (5)$$

To compare to the local balance, multiply eq. (3) by $U$ and assume that $UC = \overline{u_r C'}$, which yields a balance between buoyancy and Lorentz terms given by $gUC/\ell_u \sim B^2U/(\mu_0\ell_u^2)$. This is consistent with eq. (5) provided that

$$\frac{\ell_u}{U} \sim \frac{\ell_{Bmin}^2}{\eta} \Rightarrow \frac{\ell_{Bmin}}{L} \sim Rm^{-1/2} \left(\frac{\ell_u}{L}\right)^{1/2}, \quad (6)$$

where $Rm = UL/\eta$ is the magnetic Reynolds number. This relationship has received support from dynamo simulations (Aubert et al. 2017). Note that it differs from the classical prediction of kinematic dynamo theory where $\ell_{Bmin}/L \sim Rm^{-1/2}$ (Moffatt 1978).

Christensen & Aubert (2006) noted that the large viscosity in current dynamo simulations means that buoyant power is not all dissipated ohmically. In this case eq. (5) can be written (Davidson 2013)

$$f_{\text{ohm}} g \overline{u_r C'} \sim \frac{\eta \overline{B^2}}{\mu_0 \ell_{Bmin}^2}, \quad (7)$$

where $f_{\text{ohm}} = D_O/P_A$. Defining the convective power density $p_A$ as

$$p_A = \frac{g \overline{u_r C'}}{\rho} \approx \frac{gUC}{\rho} \sim \frac{P_A}{V_{oc}} \quad (8)$$

gives a scaling for $B$ as

$$B^2 \sim f_{\text{ohm}} \rho \mu_0 \frac{\ell_u}{U} p_A. \quad (9)$$

Eq. (9) together with the thermal wind balance

$$\frac{U\Omega}{L} \sim \frac{p_A}{U\ell_u} \quad (10)$$

provide two equations to determine the three unknowns $B$, $U$ and $\ell_u$. Starchenko & Jones (2002) assumed that at low $E$ the magnetic field prevents the flow length scale from falling as $E^{1/3}$ and instead sets $\ell_u$ to a fixed fraction of $L$. In this case eq. (10) gives $U^2 \sim p_A/\Omega$ and

$$B^2 \sim f_{\text{ohm}} \rho \mu_0 L \Omega^{1/2} p_A^{1/2}. \quad (11)$$

Alternatively, Davidson (2013) assumed that the field strength is independent of the diffusion coefficients and rotation rate. Dimensional analysis then leads to the result

$$B^2 \sim f_{\text{ohm}} \rho \mu_0 L^{2/3} p_A^{2/3}. \quad (12)$$

Recent high-resolution direct numerical simulations (Aubert 2019) produce behaviour that is more consistent with eq. (12) than eq. (11), however, these simulations still do not entirely adhere to the theory of Davidson (2013). We therefore consider whether the two scalings can be distinguished based on their predictions of modern and palaeomagnetic field behaviour. The scaling laws derived above strictly determine the internal field strength. However, they are in principle valid for describing the field at the CMB if the same balance of terms also holds near the top of the core.

Eqs. (11) and (12) are both QG-MAC balances; the difference arises in the treatment of the convective length scale $\ell_u$. Starchenko & Jones (2002) fix $\ell_u$ to a fixed fraction of $L$ and then use eq. (3) to obtain the unknowns $U$ and $B$ in terms of $p_A$. Davidson (2013) allowed $\ell_u$ to be determined from the vorticity balance, which requires an additional piece of information, in this case that $B$ is independent of the rotation rate and diffusion coefficients. For this reason, we label the scaling (11) as QG-MAC-fixed and the scaling (12) as QG-MAC-free.

## 2.2 Dynamo simulations

We use a total of 314 dynamo simulations, of which 193 employ fixed flux (FF) conditions at the outer boundary as is appropriate for modelling Earth's core. The remaining 121 are driven by a fixed temperature (FT) contrast and are used for comparison purposes since much of the previous work on field strength scaling has employed this setup (Christensen & Aubert 2006). The simulations are from Aubert et al. (2009), Yadav et al. (2016), Christensen et al. (2010), Christensen (2010), Aubert et al. (2017), Schwaiger et al. (2019), Aubert (2019), Davies & Gubbins (2011), Davies & Constable (2014), Sprain et al. (2019) and Meduri et al. (2021). All studies scale length by $L = r_o - r_i$ and define
the Prandtl and magnetic Prandtl numbers as

$$ Pr = \frac{\nu}{\kappa}, \quad Pm = \frac{\nu}{\eta}. \qquad (13) $$

Relations between the different conventions for defining the Ekman number *E*, characteristic velocity *U*, characteristic magnetic field *B* and power density *p* can be established by focusing on the definitions used in Aubert *et al.* (2009), Christensen *et al.* (2010) and Davies & Constable (2014), which are denoted by subscripts A, C and D respectively:

$$ E_A = \frac{\nu}{\Omega L^2}, \quad U_A = L \Omega U_A^*, \quad B_A = \sqrt{(\rho \mu_0) \Omega L B_A^*}, \quad p_A = \rho \Omega^3 L^2 p_A^* $$

$$ E_C = \frac{\nu}{\Omega L^2}, \quad U_C = \frac{\nu}{L} U_C^*, \quad B_C = \sqrt{(\Omega \eta \mu_0 \rho) B_C^*}, \quad p_C = \rho \frac{\nu^3}{L^4} p_C^* $$

$$ E_D = \frac{\nu}{2\Omega L^2}, \quad U_D = \frac{\eta}{L} U_D^*, \quad B_D = \sqrt{(2\Omega\eta\mu_0\rho) B_D^*}, \quad p_D = \rho \frac{\eta^3}{L^4} p_D^* $$

where asterisks denote dimensionless quantities. Here, we use the ‘diffusionless’ units of Aubert *et al.* (2009) and convert all quantities to these units. This choice is suggested by the scaling laws, which do not contain the diffusion coefficients, while Christensen (2010) also found that the choice of units was not critical for the overall results. Converting the various definitions of *p* to diffusionless units requires that

$$ p_A^* = 8 \left( \frac{E_D}{Pm} \right)^3 p_D^* = E_C^3 p_C^*. \qquad (14) $$

The diffusionless measure of field strength is the Lehnert number *Le*,

$$ Le = \frac{B}{\sqrt{(\rho\mu_0)\Omega L}}, \qquad (15) $$

which coincides with the dimensionless $B_A^*$ above. The relevant conversions are:

$$ Le = \sqrt{\frac{4\Lambda_D E_D}{Pm}} = \sqrt{\frac{\Lambda_C E_C}{Pm}}, \qquad (16) $$

where $\Lambda_D = B^2/(2\rho\mu_0\eta\Omega) = \Lambda_C/2$ is the Elsasser number based on the field strength scalings defined above. With these definitions eqs (11) and (12) become

$$ \begin{aligned} Le &\sim f_{ohm}^{1/2} (p_A^*)^{1/4} \quad (\text{QG-MAC-fixed}), \\ Le &\sim f_{ohm}^{1/2} (p_A^*)^{1/3} \quad (\text{QG-MAC-free}). \end{aligned} \qquad (17) $$

Henceforth, we will drop the asterisks on dimensionless quantities.

The simulations are split into groups based on the boundary conditions and heating mode. For simulations that employ homogeneous boundary conditions and standard setups we distinguish between FT, FF and zero flux (0F) conditions on the buoyancy source, which can be thermal, chemical, or a combination of both. Four-letter acronyms such as FTFT denote conditions on the inner and outer boundaries respectively. The final groups are the Coupled Earth (CE) simulations of Aubert *et al.* (2017), Aubert (2019) and Aubert & Gillet (2021) and the ‘mixed’ group of simulations, which both use complex driving modes and boundary conditions. The groups are:

FTFT: Yadav *et al.* (2016) and Schwaiger *et al.* (2019) both consider simulations driven by a fixed temperature contrast, with no-slip and insulating boundary conditions. Yadav *et al.* (2016) report 30 simulations with $Pr = 1, 10^{-6} \le E_C \le 10^{-4}$, $Pm = 1$ at $E_C > 10^{-6}$ and $0.4 \le Pm \le 2$ for $E_C = 10^{-6}$ and $r_i/r_o = 0.35$. Schwaiger *et al.* (2019) report 95 simulations with $Pr = 1, 10^{-6} \le E_C \le 10^{-4}$, $0.07 \le Pm \le 15$ and $r_i/r_o = 0.35$.

FF0F: The Christensen (2010) data set uses no-slip and insulating boundary conditions with a fixed codensity flux at the inner boundary and zero flux at the outer boundary. The simulations span the parameter ranges $Pr = 1, 3 \times 10^{-6} \le E_C \le 10^{-3}$, $0.5 \le Pm \le 40$ and $r_i/r_o = 0.35$.

FTFF: Christensen *et al.* (2010) modelled thermochemical convection and employed fixed temperature on $r_i$ and fixed flux on $r_o$. These simulations span the parameter ranges $Pr = 1 - 3, 3 \times 10^{-6} \le E_C \le 3 \times 10^{-4}$, $0.5 \le Pm \le 33$ and $r_i/r_o = 0.35$.

CE: Aubert *et al.* (2013, 2017) and Aubert (2019) undertook thermochemical simulations with stress-free and electrically conducting upper and lower boundaries. The mass flux is fixed at $r_i$ and there is zero flux at $r_o$, with an internal sink term to conserve mass. In order to match prominent features of the modern geomagnetic field and its secular variation the CE simulations also include: gravitational coupling between the mantle and inner core; magnetic coupling between the liquid and solid cores; and lateral variations in mass anomaly flux at the inner and outer boundaries (Aubert *et al.* 2013). CE simulations follow a path in parameter space that is designed to preserve a constant value of $Rm \sim 1000$ and $\Lambda_C \sim 20$, starting from a simulation that is similar to the original coupled Earth models in Aubert *et al.* (2013). Consequently the simulated field strength follows the prediction of eq. (11).

Mixed: Comprises the simulations from Aubert *et al.* (2009) and a compilation of models which appeared in Davies & Gubbins (2011), Davies & Constable (2014), Sprain *et al.* (2019), Biggin *et al.* (2020) and Meduri *et al.* (2021). Aubert *et al.* (2009) reported 42 simulations of dynamo action driven by thermochemical convection using the codensity formulation. They employed FF conditions on the codensity, no-slip velocity and insulating boundary conditions for the flow and magnetic field respectively, and dimensionless parameters $Pr = 1, 3 \times 10^{-5} \le$
$E_A \le 3 \times 10^{-4}$, $1 \le Pm \le 10$ and $0.1 \le r_i/r_o \le 0.35$. Models from the other studies (Leeds models) all use no-slip boundary conditions and an insulating outer boundary, but use different conditions at the inner boundary (FT or FF, insulating or conducting) and different heating modes (bottom, internal and mixed). Some of these models also include lateral variations in the heat flow at the outer boundary or a stably stratified layer at the top of the fluid domain. The parameter ranges spanned by the Leeds models are $Pr = 1, 1.2 \times 10^{-4} \le E_D \le 10^{-3}$ and $2 \le Pm \le 20$. All except 3 simulations use $r_i/r_o = 0.35$; the others use $r_i/r_o = 0.1, 0.2$.

Overall this large simulation set gives us access to a wide range of physical conditions with which to test the two scaling laws.

## 2.3 Thermal history models

Thermal history models solve equations governing global conservation of energy, entropy and mass, averaged over timescales longer than those relevant to the dynamo process but short relative to the cooling timescale (Nimmo 2015). This averaging is assumed to remove lateral variations in temperature and composition, leaving a state that is adiabatic and chemically well-mixed outside of very thin boundary layers. Convective dynamics enter the model description by preserving the adiabatic state in the bulk of the core and through the CMB heat flow, which is set by mantle convection and will not generally equal the adiabatic heat flow. Detailed descriptions of the modelling process for the convecting core can be found in Gubbins et al. (2003, 2004), Nimmo (2015), Davies (2015) and Labrosse (2015). Here, we use the specific implementation of Greenwood et al. (2021), which models the convecting core in the same way as Davies (2015) and additionally allow regions of stable thermal stratification to develop below the CMB. In these regions, the solution follows a conductive profile, which is matched to the adiabatic and well-mixed bulk at the base of the layer.

Core composition is determined by the core mass and the part of the ICB density jump, $\Delta\rho$, that is not due to the phase change. We use the Fe-Si-O core model of Alfè et al. (2002) and Gubbins et al. (2015) in which Si partitions almost equally between solid and liquid at ICB conditions, while O partitions almost entirely into the liquid. We consider three compositions that are consistent with observational constraints of $\Delta\rho = 0.8 \pm 0.2$ gm cc$^{-1}$ (Masters & Gubbins 2003) defined by mole fractions of 82%Fe–8%O–10%Si, 79%Fe–13%O–8%Si and 81%Fe–17%O–2%Si corresponding to $\Delta\rho = 0.6, 0.8$ and $1.0$ gm cc$^{-1}$, respectively (Davies et al. 2015). The composition determines the melting point depression at the ICB, which anchors the adiabatic temperature. The contributions of all three elements to the gravitational energy and entropy terms, to the entropy of molecular diffusion, and the melting point depression are calculated separately and combined by simple addition as described in Davies (2015).

The global energy balance equates the CMB heat flow $Q_{\text{cmb}}$ to the heat sources within the core. We follow previous work and ignore small effects due to thermal contraction; we also omit radiogenic heating. The energy balance can then be written

$$Q_{\text{cmb}} = -\underbrace{\frac{C_p}{T_{\text{cen}}} \int_{Q_s} \rho T_a dV \frac{dT_{\text{cen}}}{dt}}_{Q_s} - \underbrace{4\pi r_i^2 L_h \rho_i C_r \frac{dT_{\text{cen}}}{dt}}_{Q_L} + \underbrace{\alpha_c \frac{D c_X^l}{Dt} \int_{Q_g} \rho \psi dV_{\text{oc}}}_{Q_g} \quad (18)$$

where $Q_s$ is the secular cooling and $Q_L$ and $Q_g$ are respectively the latent heat and gravitational energy released on freezing. The rate of change of light element X with mass fraction $c_X^l$ in the liquid is

$$\frac{D c_X^l}{Dt} = \frac{4\pi r_i^2 \rho_i}{M_{\text{oc}}} C_r (c_X^l - c_X^s) \frac{dT_{\text{cen}}}{dt} \quad (19)$$

and

$$C_r = \frac{1}{(dT_m/dP)_{r=r_i} - (\partial T_a / \partial P)_{r=r_i}} \frac{1}{\rho_i g_i} \frac{T_i}{T_{\text{cen}}} \quad (20)$$

relates the rate of change of the ICB radius to the cooling rate $dT_{\text{cen}}/dt$ at Earth's centre. Here the density $\rho(r)$, gravity $g(r)$, gravitational potential $\psi(r)$ (referred to zero potential at the CMB), pressure $P(r)$, adiabatic temperature $T_a(r)$, melting temperature $T_m(P)$ and entropy of melting $\Delta s(P)$ are functions of $r$ and are represented by polynomials (Davies 2015). Subscripts i, o and cen refer to quantities that are evaluated at the ICB, CMB and centre of the core respectively, while the subscript oc refers to the outer core. The volume and mass of the whole core are denoted by $V$ and $M$ respectively. In writing eq. (18), the CMB has been assumed to be electrically insulating, consistent with the dynamo simulations, and the specific heat capacity at constant pressure $C_p$ and compositional expansion coefficient $\alpha_c = \rho^{-1}(\partial\rho/\partial c_X)_{P,T}$ are constants. The latent heat coefficient is $L_h = T_a \Delta s$.

The magnetic field appears through the ohmic dissipation $E_J$ in the entropy balance, which reads

$$\frac{1}{\mu_0^2} \underbrace{\int \frac{(\nabla \times \mathbf{B})^2}{T_a \lambda} dV}_{E_J} + \underbrace{\int k \left(\frac{\nabla T_a}{T_a}\right)^2 dV}_{E_k} + \underbrace{\alpha_c^2 \alpha_D \int \frac{g^2}{T_a} dV}_{E_a} = \underbrace{\frac{C_p}{T_{\text{cen}}} \left(M - \frac{1}{T_{\text{cen}}} \int \rho T_a dV \right) \frac{dT_{\text{cen}}}{dt}}_{E_s} - \underbrace{Q_L \left(\frac{T_i - T_{\text{cen}}}{T_i T_{\text{cen}}}\right)}_{E_L} + \underbrace{\frac{Q_g}{T_{\text{cen}}}}_{E_g} \quad (21)$$

Here, $\lambda$ is the electrical conductivity and $\alpha_D$ is defined precisely in Gubbins et al. (2004) and Davies (2015), however it is not important as the entropy $E_a$ produced by barodiffusion is small. $E_k$ is the entropy due to thermal conduction, which depends on the thermal conductivity $k$.
Eqs (18) and (21) can be written in the compact form (Gubbins et al. 2004; Nimmo 2015)

$$ 
\begin{aligned}
Q_{\text{cmb}} &= (\tilde{Q}_s + \tilde{Q}_L + \tilde{Q}_g) \frac{dT_{\text{cen}}}{dt}, \\
E_J + E_k + E_a &= (\tilde{E}_s + \tilde{E}_L + \tilde{E}_g) \frac{dT_{\text{cen}}}{dt},
\end{aligned}
\qquad (22) $$

where the tilde quantities are defined such that $Q_s = \tilde{Q}_s dT_{\text{cen}}/dt$ and similarly for other terms. For given CMB heat flow the energy balance determines the cooling rate $dT_{\text{cen}}/dt$, which is then used in the entropy balance to obtain $E_J$. The ohmic dissipation differs from the ohmic heating $D_O$ by the factor of $1/T_a$ under the integral in eq. (21). We write $D_O \approx E_J T_{\text{mean}}$, where $T_{\text{mean}}$ is the average core temperature (Nimmo 2015). Neglecting viscous heating allows $P_A$ to be obtained from eq. (4):

$$ P_A = D_O + D_V \approx E_J T_{\text{mean}}. \qquad (23) $$

Core properties for the three values of $\Delta\rho$ are listed in table 1 of Davies et al. (2015). The only other model input is the CMB heat flow, which must be specified over the 4.5 Gyr evolution. In principle $Q_{\text{cmb}}$ can be calculated using a parametrized model of mantle convection that is coupled to the core evolution, thus allowing changes in core temperature to alter the heat flow and vice versa (e.g. Nimmo et al. 2004; Driscoll & Bercovici 2014; O'Rourke et al. 2017). However, such a complicated process is not required here, where the goal is to understand long-term variations in magnetic field strength. We therefore use a simple parametrized form

$$ Q_{\text{cmb}} = Q_P \exp^{(4.5-t)/\tau}, \qquad (24) $$

where $Q_P$ is the present-day heat flow at time $t = 4.5$ Gyr and $\tau$ is a timescale. Eq. (24) can approximate a wide range of plausible heat flows including those obtained from coupled core-mantle evolution models (e.g. Driscoll & Bercovici 2014) and 3-D mantle convection simulations (e.g. Nakagawa & Tackley 2014).

Regions of stable thermal stratification can develop if the CMB heat flow becomes subadiabatic (e.g. Lister & Buffett 1998). The thermal conduction equation is solved in the layer, with FF conditions at the CMB and layer base. The layer thickness evolves over time in order to preserve continuity of temperature at the interface. In the models presented here the layers do not grow past 300–400 km and their effect on the bulk evolution is small (Greenwood et al. 2021).

Eqs (22) are time stepped forward from 4.5 Ga to the present with a time step of 1 Myr. At each step, the cooling rate is obtained and used to calculate the temperature and composition at the following step. Initially the core is entirely molten; the inner core begins to grow when $T_a$ drops below $T_m$ at Earth's centre and the ICB radius is tracked over time using the intersection point $T_a = T_m$. The outputs are time-series of $E_J$, temperature at the CMB $T_o$, bulk composition, ICB radius $r_i$ and radius of the stable layer base $r_s$. All reported models are required to satisfy two basic criteria. First, the entropy production $E_J$ must remain positive over the last 3.5 Ga, consistent with palaeomagnetic evidence indicating the persistence of a global field over this period. Second, the model must match the present-day ICB radius to within 10 per cent.

We have conducted 275 thermal history models spanning the parameter space $\Delta\rho = 0.6, 0.8$ and $1.0$ gm cc$^{-1}$, $Q_P = 6-18$ TW (increasing in increments of 1 TW) and $\tau = 2-20$ Gyr (increasing in increments of 1 Gyr). Many of the models fail to produce a dynamo for the whole of Earth's history because $E_J$ falls below zero prior to ICN. This places an upper limit on the allowed value of $\tau$ for fixed $Q_P$. At lower $Q_P$, lower values of $\tau$ are needed to maintain the dynamo, which corresponds to a larger change in CMB heat flow over time.

When determining the TDM time-series for the palaeofield we use the dimensional scaling laws given by eqs (11) and (12) with $\rho = 10^4$ kg m$^{-3}$. Time variations in the shell thickness, $L$, are calculated using the values of $r_i$ and $r_s$ from the thermal history models. A thermal wind flow could arise in the stable layer, in which case it may be more appropriate to calculate $L$ using $r_o$ rather than $r_s$; however, in practice, stable layers rarely emerge in our models and always remain thin, so we do not expect this to significantly affect the results. For $\Omega$, we use the same piecewise linear model as in Aubert et al. (2009) in which the length of day increases from 17 hr at 4.5 Ga to 19 hr at 2.5 Ga to 20.8 hr at 0.64 Ga, and finally to 24 hr today.

# 3 RESULTS

In this section, we first compare the two theoretical scaling laws for Le given by eqs (17) to the results of numerical dynamo simulations. We then present the palaeointensity data set and calculate TDMs for 275 thermal history models that span a wide range of plausible evolutionary scenarios for the core.

## 3.1 Scaling laws for dynamo field strength

We consider the RMS field strength inside the dynamo region, the RMS CMB field strength and the dipole field strength on the CMB, which are defined respectively as

$$ B_{t}^{\text{rms}} = \sqrt{\frac{1}{V_{\text{oc}}} \int \mathbf{B}^{2} dV}, \quad B_{\text{cmb}}^{\text{rms}} = \sqrt{\frac{1}{S} \int \mathbf{B}^{2} dS}, \quad B_{\text{cmb}}^{\text{dip}} = \sqrt{\frac{1}{S} \int \mathbf{B}_{\text{dip}}^{2} dS} \qquad (25) $$
where S is the surface area of the outer boundary and superscript 'dip' refers to the spherical harmonic degree 1 component of the field. All quantities are time-averaged. For each simulation data set, we compute the Lehnert numbers corresponding to these three definitions of the field strength. Yadav et al. (2016) provide the axial CMB dipole field strength, which omits the contributions to the total CMB dipole from spherical harmonic order 1. We do not expect this to influence the results since these terms tend to be much smaller than the axial dipole.

For each individual data set and for the combined data set of 314 simulations we seek the constants c and m that provide the best least-squares fit between the data and an equation of the form

$$Le/f_{\text{ohm}}^{1/2} = c p_A^m. \qquad (26)$$

The theoretically predicted values of m are 1/4 and 1/3 for the QG-MAC-fixed and QG-MAC-free scaling laws respectively (see eqs 17). The prefactors c are not determined by the theory, but should be approximately constant in order for the theory to have captured the dominant parametric dependence of Le. The formal least-squares uncertainty on m is always small and so we also quote the sum of squared residuals (SSR) when comparing results. Following Aubert et al. (2009), we also calculate the vertical standard deviation σ, which is based on the prefactor c using a least-squares fit to the simulation data with the exponent m fixed to the theoretical values determined by the QG-MAC-free and QG-MAC-fixed scaling laws.

It is vital to filter the simulation data set when assessing the fits to theoretical scaling laws. Though eqs (17) do not depend on the topology of the field (Christensen 2010), when applying the results to Earth it is important to focus on dipole-dominated fields. Moreover, the dominant force balance can change significantly as control parameters are varied, with viscous and inertial effects perturbing the expected QG-MAC balance that emerges as more realistic conditions of low E and Pm are approached (Aubert et al. 2017; Schwaiger et al. 2019). In this work, we use two different quantities to filter the simulation data set:

$f_{\text{dip}}$: the time-averaged ratio of the dipole CMB field strength to the RMS strength of all CMB field components up to spherical harmonic degree 12 (Christensen & Aubert 2006). This filter allows to remove simulations that are too dipolar (high $f_{\text{dip}}$) and also multipolar fields (low $f_{\text{dip}}$). Plausible values of $f_{\text{dip}}$ for Earth should exceed 0.4–0.5, which approximately marks the dipole-multipole transition (Christensen & Aubert 2006; Oruba & Dormy 2014). The upper value must include the modern field, for which $f_{\text{dip}} \approx 0.64$ for the CHAOS6 model spanning the last 10 yr (Finlay et al. 2016), and $f_{\text{dip}} \approx 0.70 \pm 0.03$ for the gufm1 model since 1840 (Jackson et al. 2000). Another factor to consider is that weakly driven dynamos, which generally have high $f_{\text{dip}}$, can display significant viscous effects that are not expected to exist in the core. From these considerations Aubert et al. (2009) focused on the range $0.35 \le f_{\text{dip}} \le 0.7$, while Christensen (2010) chose $0.45 \le f_{\text{dip}} \le 0.75$. Here, we report three sets of results: no filter; $f_{\text{dip}} > 0.5$, which conservatively removes multipolar solutions; and the range $0.35 < f_{\text{dip}} < 0.75$.

$E_M/E_K$: the ratio of total magnetic to kinetic energy in the domain. Schwaiger et al. (2019) analysed the force balance in a suite of 95 dynamo simulations and found that the value of $E_M/E_K$ provided a convenient proxy for filtering out dynamos that were not in QG-MAC balance. The critical value of $E_M/E_K$ is around 1 (see Schwaiger et al. 2019, fig. 3) and we test values in the range $E_M/E_K = 0-5$.

Fig. 1 shows fits of m and c to the dynamo simulations for different $f_{\text{dip}}$ and $E_M/E_K$ filters. Quoted c values are calculated by fixing $m = 1/3$, corresponding to the predicted QG-MAC-free scaling. For the RMS internal field, the values of m and c are generally consistent as long as some filtering of the data set has been performed and are tightly clustered for $E_M/E_K \ge 2$. For the CMB dipole field, consistent values of m and c only emerge when $E_M/E_K$ exceeds 2 or 3; indeed, for $E_M/E_K \ge 2$ the variations are at most ~5 per cent for m and ~20 per cent for c. Increasing the critical value of $E_M/E_K$ (below which simulations are filtered out) from 1 to 5 reduces the number of simulations from 225 to 110. In this section, we therefore focus on the case where all simulations with $E_M/E_K < 2$ are filtered out, which produces similar m and c to the more restrictive filters while retaining more data. The resulting data set contains 17 simulations with $r_i/r_o$ that differs from the present-day value; we have verified that retaining these data produces at most a 1 per cent change in the quoted values of m and c.

Fig. 2 shows $Le_{\text{rms}}^m$, $Le_{\text{cmb}}^m$ and $Le_{\text{cmb}}^{\text{dip}}$ computed from eq. (25) as a function of $p_A$ for simulations where $E_M/E_K \ge 2$. For the internal field $Le_{\text{rms}}^m$, the fit to the FTFT data set is close to the QG-MAC-free prediction, which is expected for fixed temperature boundary conditions (Christensen & Aubert 2006). The FF0F simulations fall close to an exponent of $m = 0.25$ as would be expected from a QG-MAC-fixed balance and are not compatible with the QG-MAC-free balance to within the formal uncertainty. The CE simulations also fall close to the $m = 0.25$ scaling as expected because most use a large-scale approximation that fixes the dominant length scale. Notwithstanding these ‘shingling’ effects (Cheng & Aurnou 2016) the best-fitting exponent to the overall data set is $m = 0.32$, in excellent agreement with the QG-MAC-free prediction.

Fits to the RMS CMB field $Le_{\text{cmb}}^{\text{rms}}$ and dipole CMB field $Le_{\text{cmb}}^{\text{dip}}$ (Fig. 2) are similar to the internal field except with more scatter. In both cases the SSR increases by a factor of roughly 2 for all data sets except mixed when compared to the internal field, perhaps in part because of the different spatial averaging. For each simulation grouping the best-fitting exponents are similar between internal and CMB fields, often overlapping within the formal errors. The overall data set displays a clear dependence of $Le_{\text{cmb}}^{\text{dip}}$ on $p_A$, with the vast majority of simulations falling within the $1\sigma$ uncertainty on c (shown by the grey shading in Fig. 2), and SSRs that are comparable to those of the RMS CMB field. The best-fitting exponent to $Le_{\text{cmb}}^{\text{dip}}$ for the overall data set is $m = 0.31$, again in excellent agreement with the QG-MAC-free prediction.

As well as matching simulation data, a viable scaling law should give a reasonable estimate of Earth’s present-day field strength. The ohmic dissipation in the core (which is a proxy for $p_A$) cannot be observed and so we take a wide range of values, $0.1 \le D_O \le 5$ TW, which spans estimates derived from thermal history models (Nimmo 2015; Davies 2015; Labrosse 2015) and scaling analysis (Christensen & Tilgner 2004). For the internal field strength, we use the range 1–10 mT, which spans inferences from satellite field models (Finlay et al. 2016), tidal dissipation (Buffett 2010) and torsional wave periods (Gillet et al. 2010). For the axial dipole field, we take the range 20–40 µT at the surface based on variations observed in the historical (Jackson et al. 2000) and Holocene (Constable et al. 2016) fields.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** SSR versus exponent *m* (top) and prefactor *c* (bottom) for each of the 18 different filters. Squares, circles and triangles show no *f*<sub>dip</sub> filter, *f*<sub>dip</sub> > 0.5 and 0.35 < *f*<sub>dip</sub> < 0.75, respectively, while colours distinguish the filters *E*<sub>M</sub>/*E*<sub>K</sub> = 0–5. Prefactors are calculated by fixing *m* = 1/3, corresponding to the predicted QG-MAC-free scaling. Note that each point is a fit to the (filtered) simulation data set. For the CMB dipole field, the SSR obtained from fitting the unfiltered data set plots above the vertical range shown.</figcaption></figure>

Fig. 3 shows simulation fits and extrapolations for the internal and CMB dipole fields when filtering out all simulations with $E_M/E_K < 2$. For the internal field both QG-MAC-free and QG-MAC-fixed scalings match the modern-day geomagnetic field strength when extrapolated based on the best-fitting *c* value obtained with *m* fixed to the theoretical prediction, though QG-MAC-free provides a better fit to the simulations. For the dipole CMB field the QG-MAC-fixed scaling overpredicts Earth's field strength even given the generous uncertainty bounds, while the QG-MAC-free prediction matches Earth's field strength.

Fig. 3 also shows that simulations with higher *Rm* tend to have lower $Le_{cmb}^{\text{dip}}$ at similar $p_A$, while for $Le_t^{\text{rms}}$ the *Rm* dependence is reduced. To clarify this point Fig. 4 shows $b_{\text{dip}} = Le_t^{\text{rms}}/Le_{cmb}^{\text{dip}}$ as a function of $p_A$ with simulations coloured by *Rm*. There is some dependence of $b_{\text{dip}}$ on the simulation boundary conditions and heating mode as found in Aubert *et al.* (2009), but relatively little dependence on $p_A$. The clear result is that the simulations are systematically biased low, with most $b_{\text{dip}}$ values in the range 4–8 compared to modern Earth values of 10–16. Simulations at higher *Rm* come closer to matching the Earth value of $b_{\text{dip}}$. A potential explanation for this observation is that higher *Rm* reduces the diffusion of field across the outer boundary. The CE simulations come closest to realistic $b_{\text{dip}}$ values because they can reach high *Rm* while remaining at low *E* and *Pm* such that they maintain QG-MAC balance. We will return to this point when comparing synthetic field strength predictions to the palaeofield.

Taken together these results provide support for a relationship between the dipole CMB field and the total power available to drive the dynamo and favour the QG-MAC-free scaling theory of Davidson (2013). In the following sections, we compare both QG-MAC-free and QG-MAC-fixed predictions to the PINT data set to establish whether palaeointensity data can help distinguish between the two predictions.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Field strength as a function of convective power $p_A$ for 225 simulations with $E_M/E_K \ge 2$. The top panel shows the internal field strength $Le_i^{rms}$, middle shows the RMS CMB field strength $Le_{cmb}^{rms}$ and bottom shows the dipole CMB field strength $Le_{cmb}^{dip}$. In each panel, the symbol colour denotes the different simulation types as described in the text. Power-law exponents $m$ for each data set are written in the corresponding colour and the fit for the whole data set is written in black together with the corresponding SSR. The black line is the best fit to the whole data set with $\pm 1\sigma$ uncertainties on the prefactor $c$ shown in grey shading. Symbols are shaded according to the magnetic Reynolds number $Rm$.</figcaption></figure>

We do this by fixing the exponent to the theoretically determined values and using two values of the prefactor as described below. Together with time-series of $p_A$ and $L$ from the thermal history models and the variation of $\Omega$, this completely determines $Le_{cmb}^{dip}$ and hence the TDM from each of the two scaling laws.

## 3.2 Comparing synthetic and observed dipole moment

TDMs obtained from core thermal history models are compared to an expanded version of the PINT data set (Biggin *et al.* 2015), which reports field strength observations at the site-mean (i.e. cooling unit) level. The expanded data set includes new palaeointensity data (Table
<figure><img src="image_5.png" /><figcaption>**Figure 3.** RMS internal field (top) and CMB dipole (bottom) as a function of convective power $p_A$ extrapolated to Earth's core conditions (shaded regions). The data set is filtered by $E_M/E_K \ge 2$. In each panel, the symbol colour denotes the different simulation groupings as in Fig. 2. Symbols are shaded according to the magnetic Reynolds number $Rm$. Power law exponents $m$ and SSRs for each data set are provided with the best fit, $1\sigma$ uncertainty (light dashed black lines) and $2\sigma$ uncertainty (grey shading) for the whole data set. Theoretical predictions based on the $m = 1/3$ and $1/4$ scalings are shown by dashed and dotted green lines with $1\sigma$ uncertainty for the $m = 1/3$ case based on the prefactor $c$ shown by green shading.</figcaption></figure>

1, Supporting Information), the fixes and modifications reported by Kulakov *et al.* (2019), and the removal of select site means which record altered or secondary magnetizations following Smirnov *et al.* (2016) and Bono *et al.* (2019).

We filtered the PINT data set by only including studies that used the following methods to identify laboratory alteration: low-temperature Shaw method ('LTD-DHT-S'; Yamamoto & Tsunakawa 2005), Low-temperature Thellier with partial thermoremanent (pTRM) tail checks ('LTD-T+'; Yamamoto *et al.* 2003), microwave technique with pTRM checks ('M+'; Shaw 1974), Multi-Specimen Parallel Differential Technique ('MSPDp'; Dekkers & Böhnel 2006), Shaw & Thellier ('ST+'), Thellier or variant with pTRM checks ('T+'; Thellier & Thellier 1959), Thellier with pTRM checks and correction ('T+Tv'; Valet *et al.* 1996), Wilson (Wilson 1961) & Thellier with pTRM checks ('WT+'). This yielded a data set containing 2780 field strength observations. We considered further restrictions by requiring $\ge 3$ intensity observations and published $Q_{PI}$ scores $\ge 3$ (Biggin & Paterson 2014), which reduced the data set to 407 observations with most of the exclusions occurring
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Ratio $b_{\text{dip}}$ of the total internal RMS field strength $Le_{\text{rms}}^i$ and the dipole CMB field strength $Le_{\text{cmb}}^{\text{dip}}$ as a function of convective power $p_A$ for all simulations with $E_M/E_K \ge 2$. The magnetic Reynolds number $Rm$ is shown in the colour bar and symbol colours are as in Fig. 2. Values of $b_{\text{dip}}$ for the modern Earth are shown by dashed lines using estimates of the internal field strength from Buffett (2010) and Gillet et al. (2010).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** VADM estimates from PINT observations. Diamonds: all PINT data; blue squares: PINT data meeting additional criteria; black circles: 200 Myr bin median included in our analysis; red circles: 200 Myr bin median not included in our analysis; red crosses: Tarduno et al. (2015) zircon palaeointensity data from single heating step experiments (not included in bin median estimates). Horizontal error bars show minimum and maximum ages for each bin; vertical error bars show inter-quartile range of V(A)DMs. Dotted line shows present day field of ~8 × 10²² Am².</figcaption></figure>

in the last 200 Myr. However, given the overall similarity between the data sets and the large reduction in data (~78 per cent) we chose not to proceed with the more stringent criteria.

Fig. 5 shows the individual data, which are unevenly distributed in time with ~75 per cent of data in the last 200 Ma. We therefore group data into bins that each span 200 Myr, which should sufficiently average secular variation (occurring on timescales of up to 1 Myr), while allowing for the longest term variations (due to secular thermochemical evolution) to be detected. Bins spanning 600–800, 2000–2200, 2800–3000 and 3000–3200 Myr contained no data. Furthermore, bins at 400–600, 800–1000, 1400–1600, 1800–2000, 2200–2400, and 3400–3600 Myr contained only 1, 2, 5, 8, 2 and 7 data points respectively and so these bins (marked by red dots in the figures) were not considered further, leaving a total of $N_b = 8$ bins.

We compare theoretical TDMs, $T_i$, obtained from 275 core thermal history models with the median of the VDM and virtual axial dipole moment (VADM) observations in the ith bin, $V_i$, using the RMS uncertainty:

$$ \text{RMS} = \sqrt{\frac{1}{N_b} \sum_{i=1}^{N_b} (V_i - T_i)^2}. \qquad (27) $$

Using a weighted $\chi^2$ misfit yields similar results to the RMS once the sparsely populated bins (which also have low uncertainties and thus bias the $\chi^2$ estimate) are removed. Misfits for each scaling law are denoted RMS$_j$, where j represents QG-MAC-free or QG-MAC-fixed.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Best-fitting prefactor $c_p$ from PINT data (blue) for the (a) QG-MAC-free and (b) QG-MAC-fixed scalings laws using TDM predictions from 275 thermal history models. The red distribution shows the range of $c_D$ values determined using all simulation data sets with an SSR below 6 (see Fig. 1 for the complete set of prefactors determined for the QG-MAC-free scaling law). Vertical bars show mean (solid), $1\sigma$ (dashed), $2\sigma$ (dotted–dashed) and $3\sigma$ (dotted) bounds based on fitting the dynamo simulation data filtered using $E_M/E_K \ge 2$.</figcaption></figure>

When making direct comparisons, it should be acknowledged that site level palaeomagnetic observations record instantaneous ‘snapshots’ of Earth’s field, which can vary in strength on short timescales (<1 Myr), whereas thermal history TDMs characterize slowly changing core conditions which change on timescales >1 Myr. Synthetic TDMs will therefore provide at best a smoothed representation of the palaeofield behaviour. Both TDM determinations from thermal history models and VDMs grouped in 200 Myr bins should represent a long enough duration that average estimates are robust irrespective of the dynamical state of the core (Driscoll & Wilson 2018).

To specify the scaling prefactor $c$ we compare in Fig. 6 the best-fitting estimates $c_D$ obtained from dynamo simulations to the estimate $c_p$ that minimizes (in a least-squares sense), the RMS error between the binned PINT observations and synthetic dipole moments obtained from the thermal history models. $c_D$ is calculated by fixing the exponent $m$ as determined by the QG-MAC-free or QG-MAC-fixed scaling and fitting to the simulations using all filters shown in Fig. 1 that yield an SSR below 6 (thus removing data sets that are too scattered), while $c_p$ is calculated for each of the 275 thermal histories for both scaling laws. The estimated $c_p$ values fall below $c_D$ for all filters, which is expected because the lower $Rm$ in most simulations compared to Earth’s core leads to higher $Le_{\text{emb}}^{\text{dip}}$ (Fig. 4). For QG-MAC-free, the best-fitting distribution of $c$ values from PINT is between $2\sigma$ and $3\sigma$ below that preferred by the simulations, while for the QG-MAC-fixed scaling the best-fitting PINT distribution sits between the $5\sigma$ and $6\sigma$ bounds. Therefore, for the QG-MAC-free scaling we consider two estimates of the prefactor: $c = 0.23$, a median value among the different filters used in Fig. 1 and corresponding directly to the filter with $E_M/E_K \ge 2$; $c = 0.2$, corresponding to the filter with $E_M/E_K \ge 2$ and $0.35 \le f_{\text{dip}} \le 0.75$ (Fig. 1), which we expect to better fit the PINT data set. For the QG-MAC-fixed scaling, we consider the lowest estimate of $c = 0.0749$ across all filters, which still produces TDMs that far exceed those from PINT as we show below.

Two example thermal history solutions are shown in Fig. 7 together with the predicted TDM. For $\tau < 16$ Gyr, the general behaviour consists of a gradual decline in TDM from 4.5 Ga until ICN, at which time the field strength increases rapidly before peaking and declining towards the present day. The pre-ICN TDM decline arises due to the rapid fall in $Q_{\text{emb}}$ and $D_O$, while the recent decline arises both from the decrease in $D_O$ and the decreasing volume of the liquid core. Changes in $\Omega$ are minor by comparison since it does not vary significantly over time and enters into the scaling laws raised to a low power. For models with $\tau \ge 16$ Gyr, the TDM gradually increases from 4.5 Ga to ICN, at which time it jumps sharply before plateauing. The slow rise in TDM reflects the almost constant $D_O$ before ICN while the recent plateau reflects the balance between increasing $D_O$, which increases TDM, and decreasing core volume and temperature, which decrease TDM. In both cases, the QG-MAC-fixed prediction produces TDMs that are too high to match PINT at all times (Fig. 7). Indeed, Fig. 8 shows that across all 275 models the QG-MAC-free scaling yields the lowest misfit to PINT and so we henceforth focus on this scaling.

Fig. 9 shows RMS misfit for the QG-MAC-free scaling for all $Q_p$ and $\tau$ combinations and the two chosen values of $c$. Here, white regions of the plot denote non-viable models that either failed to generate a dynamo for the last 3.5 Gyr or where the present ICB radius failed to match its seismically determined value. In all cases the models with lowest RMS plot at the interface separating viable and non-viable models. This behaviour arises because the PINT V(A)DM data are relatively flat, which favours high $\tau$, while the predicted present-day TDMs tend to be slightly higher than the PINT average, favouring low $D_O$ and hence low $Q_p$. However, if $\tau$ becomes too large the TDM is too flat and cannot match the general trend of weakening V(A)DM from 3.5 Ga to ~500 Ma observed in palaeomagnetic studies (e.g. Biggin et al. 2015; Bono et al. 2019). As expected, lower $c$ corresponds to lower misfit while also pushing the preferred solution to lower $\tau$ and lower $Q_p$, which corresponds to a lower present-day field strength and a steeper decline in TDM from 4.5 Ga to before ICN.
<figure><img src="image_9.png" /></figure>

**Figure 7.** Two example thermal history calculations together with predicted and observed field strength. The upper panel shows the input CMB heat flow $Q_{\text{cmb}}$ (black) and resulting ohmic heating $D_O$ (red). $Q_{\text{cmb}}$ is defined by $Q_p = 17$ TW and $\tau = 17$ Gyr (dashed lines) and $Q_p = 13$ TW and $\tau = 2$ Gyr (solid lines). The bottom panel shows TDM for QG-MAC-free (blue) and QG-MAC-fixed (orange) scaling laws with $c = 0.20$ and $0.075$, respectively. Diamonds show PINT data, grey shading shows the range of observed field strengths, and the black dotted line denotes the present-day field strength.

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

**Figure 8.** Distributions of log(RMS) obtained from 275 thermal history models for each scaling law, comparing model TDMs with PINT VDMs. Curve shows kernel density estimation. Left-hand (right-hand) panel uses a prefactor of $c = 0.23$ (0.20) for the QG-MAC-free scaling and $c = 0.075$ (0.075) for the QG-MAC-fixed scaling.

In all models ICN occurred between 400 and 1000 Ma (Fig. 10, left), with a median predicted age of 596 Ma. The signature of ICN in the palaeointensity record depends strongly on $\tau$. With $\tau < 16$ Gyr, the minimum predicted TDM always occurs at the time of ICN (Fig. 10, right). With $\tau \ge 16$ Gyr, the minimum TDM occurs at 4.5 Ga. All thermal histories predict a strong increase in TDM directly following ICN.

# 4 DISCUSSION AND CONCLUSIONS

We have considered two power-based scaling laws for determining the strength of the internal and CMB magnetic fields produced by spherical shell convection-driven dynamos. These scaling laws predict exponents $m$ in the relation $Le/f_{\text{ohm}}^{1/2} = cp_A^m$ of $m = 0.25$ (QG-MAC-fixed) and $m = 0.33$ (QG-MAC-free). We have compared these scaling laws to a suite of 314 geodynamo simulations that span over 6 orders of magnitude in the convective power $p_A$ and over 2 orders of magnitude in field strength. We have found that both scaling laws adequately reproduce the
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Contour maps of RMS misfit defined in eq. (27) using the QG-MAC-free scaling laws for all values of *Q*<sub>P</sub> and τ. Magenta lines show thermal history model parameters yielding the lowest misfit; magenta square shows overall best fitting model parameters. Note that our models sample the whole *Q*<sub>P</sub> – τ parameter space; white regions of the plot denote models that either failed to generate a dynamo for the last 3.5 Gyr or where the present ICB radius failed to match its seismically determined value. Top row: prefactor *c* = 0.23 and bottom row: prefactor *c* = 0.20.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Left: histogram of ICN times obtained from thermal history models with kernel density estimate of probability (blue line). Right: time of ICN obtained from the thermal history models plotted against the time of the minimum in TDM using the QG-MAC-free scaling. Colour bar shows variation in RMS for the QG-MAC-free scaling using a prefactor of *c* = 0.23.</figcaption></figure>

amplitude of the present RMS internal magnetic field (Aubert *et al.* 2017); however, only the QG-MAC-free scaling of Davidson (2013) matches the present-day CMB dipole field and provides an adequate fit to the palaeofield over the last 3.5 Gyr.

Fitting individual simulation groups (as determined by differences in boundary conditions and convective driving) reveals variations in empirically derived slopes from *m* = 0.24 to 0.39, with data sets where at least one boundary is held at fixed temperature giving consistently higher exponents than data sets employing fixed-flux conditions. At high *p*<sub>A</sub>, these two groups exhibit similar amplitudes and slopes, but they appear to diverge at low *p*<sub>A</sub>, which may reflect a change in dynamics or the relative sparsity of data at more extreme conditions. The group of simulations using mixed setups is more sensitive to filtering, which perhaps reflects the greater heterogeneity in this data set. At present, the individual groups are too small to separate the role of these different factors and so we have focused on the scaling behaviour of the data set as a whole. However, we do note that predictions from individual simulation groups are broadly consistent with theoretical QG-MAC scaling laws.

To obtain a robust scaling for the CMB dipole field we have found it essential to filter the data set by the magnetic energy to kinetic energy ratio as advocated by Schwaiger *et al.* (2019). Landeau *et al.* (2017) found that changes in the buoyancy distribution can cause the
CMB dipole field behaviour to deviate from the internal field, which follows the QG-MAC-free scaling in their simulations. Our results also suggest a residual dependency of CMB field scaling on the buoyancy source, although the effect is comparable to that seen for the internal field. We also observe similar field amplitudes between data sets with different buoyancy distributions across a wide range of $p_A$. Overall, while the individual simulation groups considered here may show some differences between internal and CMB field scaling behaviour, the combined data set supports the $p_A$-dependence of the QG-MAC-free scaling for both internal and CMB fields.

The majority of our simulations use a modern day aspect ratio of $r_i/r_o = 0.35$. Lhuillier et al. (2019) studied a range of chemically driven dynamos at $E > 10^{-3}$ with a fixed buoyancy distribution and showed that $m$ displays a non-monotonic dependence on $r_i/r_o$ in the range $r_i/r_o = 0.1–0.35$. However, the values of $m$ obtained by Lhuillier et al. (2019) fall below 0.25 for the majority of aspect ratios considered, suggesting that these simulations are not in QG-MAC balance. This raises the possibility that $m$ depends on the choice of control parameters at high $E$, as well as any influence from aspect ratio. In any case, such low values of $m$ will only worsen the fit to the PINT data unless they are associated with much lower values of $c$, which is not suggested by our analysis. Interestingly, for thick shells Lhuillier et al. (2019) obtain $m = 0.33$, which is the QG-MAC-free scaling favoured by our analysis, suggesting that the $m = 1/3$ exponent describes the dependence of dipole moment on convective power over most of Earth's history.

The simulation data sets cannot yet reach the very low $p_A$ values that characterize Earth's core. It is therefore possible that the scaling behaviour changes at more extreme control parameter values (particularly lower $E$ and $Pm$), as arises in non-magnetic rotating convection (Gastine et al. 2016; Long et al. 2020). However, no evidence for a transition from the QG-MAC regime has been found down to extremely low values of $E \sim 3 \times 10^{-10}$ (Aubert & Gillet 2021). The relevant force balance must contain buoyancy (the power source for convection) and the magnetic field (the main product of dynamo action), while rotation breaks reflectional symmetry, which is thought to be crucial for sustaining large-scale magnetic fields (Tobias 2021). At low $E$ and $Pm$ inertia and viscosity become strongly subdominant in the force balance (Aubert et al. 2017; Aubert 2019) and therefore cannot perturb the QG-MAC balance. In principle, the Lorentz force could perturb the large-scale QG balance, though this has not been observed in high-resolution simulations (Schwaiger et al. 2021) and is not expected in Earth's core (Aurnou & King 2017). We therefore believe that the QG-MAC-free and QG-MAC-fixed scaling laws we have considered capture the range of dynamical balances in Earth's core that are plausible given current simulations and theory.

The theoretical scaling laws determine only the exponent of the $Le - p_A$ relation; the prefactor $c$ must be obtained by fitting simulation data. We have assumed a constant prefactor when calculating TDMs, which is clearly an oversimplification because $c$ depends on the time-dependent buoyancy sources and shell thickness. At fixed $p_A$, decreasing the inner core size from its present volume to zero has been found to produce a relative increase in $b_{dip}$ of 30–50 per cent due to the transition from dominantly bottom-driven chemical convection to internally driven thermal convection (Aubert et al. 2009; Landeau et al. 2017). Attributing this change in $b_{dip}$ entirely to the prefactor suggests a 30–50 per cent increase in $c$ from present-day to ICN, which is comparable to our estimated uncertainty on $c$ obtained from fitting all simulation groups together (Fig. 6). Our use of two different constant $c$ values and their associated uncertainties should therefore partly mitigate any effects arising from time variations in the prefactor. We also note that changes in the CMB dipole field due to changes in $p_A$ (with constant $c$) are a factor of two or more (e.g. Fig. 7) and so the main uncertainty in the calculation is the determination of $p_A$ from the thermal history models.

The scaling prefactor obtained from dynamo simulations is generally high compared to an independent constraint obtained by minimising the misfit between TDM predictions from thermal history models and PINT. We do not believe this discrepancy arises from the thermal history models as we have considered a large range of models spanning the plausible range of input parameters. Instead it appears that the available simulations which achieve QG-MAC balance are generally operating at lower $Rm$ than Earth, which promotes diffusion of field out of the core. The path models of Aubert et al. (2017) and Aubert (2019) partially overcome this problem because the effects of inertia and viscosity are sufficiently suppressed to enable high $Rm$ simulations that retain QG-MAC balance and a dipole-dominated field. These models are run along a path where $Rm \sim 1000$; however, $Rm$ in Earth's core could be twice this value if one adopts the higher values of electrical conductivity proposed in some studies (e.g. Pozzo et al. 2013). Future work should investigate whether path-type simulations at higher $Rm$ can improve the fit between simulated and palaeomagnetic field strengths.

The preceding discussion suggests that both the internal and CMB field follow the QG-MAC-free scaling law over the majority of Earth history, with effects due to variations in buoyancy sources, boundary conditions and shell thickness influencing the prefactor $c$. Time variations in CMB dipole field strength are expected to be dominated by changes in convective power rather than the prefactor. Future studies that systematically vary the convective driving modes, boundary conditions and inner core size will provide important tests of these conclusions.

Theoretical predictions of Earth's TDM evolution require coupling dynamo simulations and thermal history models. Our approach utilizes existing simulations and enables a systematic sampling of plausible core evolution scenarios, but assumes a dipole-dominated field. Alternatively, thermal history outputs can be used to set the (interdependent) core geometry and buoyancy sources in a suite of bespoke simulations that represent different stages of core evolution (Driscoll 2016; Landeau et al. 2017). However, while this approach provides the complete field at different epochs, it is restricted to a comparatively small number of simulations and thermal histories and therefore cannot yet definitively constrain long-term TDM evolution and dipole-dominance. Observations suggest that Earth's field has been dominantly dipolar over most of its history (Biggin et al. 2020), but may have undergone periods of 10–100 Myr where the dipole field is weak or absent (Shcherbakova et al. 2017; Hawkins et al. 2019). In principle, it is possible to estimate times of dipole-dominance using theoretical predictions for the dipole-multipole transition; however, the factors that determine the transition in geodynamo simulations are still debated (Christensen
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Distribution of model TDMs compared to binned PINT VDM data (black circles) using a scaling prefactor *c* = 0.23. Black diamonds show the raw PINT data, red circles denote bins that were excluded from the misfit calculation on account of having fewer than 10 data points. The coloured shaded regions show the 1σ uncertainty interval based on the scaling prefactor *c* and the dotted line shows the present day field of 8 × 10²² Am². Top, middle and bottom show Δρ = 0.6, 0.8 and 1.0 gm cc⁻¹ cases, respectively.</figcaption></figure>

& Aubert 2006; Oruba & Dormy 2014; McDermott & Davidson 2019). Further observational constraints and targeted simulation studies extended to broader parameter regimes will shed more light on this important issue.

Figs 11 and 12 compare the binned PINT database shown in Fig. 5 to the best synthetic TDM models (lowest RMS) for each Δρ and the two values of the prefactor *c* = 0.2 and 0.23 obtained from fitting the QG-MAC-free scaling to the simulation data set. Least-squares uncertainties on the TDM, σ, are calculated based on *c* with the scaling exponent fixed to *m* = 1/3. Prior to ICN most solutions show agreement with PINT at just above the 1σ level. In this period the *c* = 0.2 and Δρ = 0.6 gm cc⁻¹ model provides the best fit to the data, even matching to many of the bins that are sparsely sampled by available data (red circles in Fig. 12) and agreeing well with the empirical fit of Bono et al. (2019). Strictly the small differences in misfits between high and low *c* for fixed Δρ mean that is it difficult to differentiate between an overall decline or near-constant field strength on the Gyr timescale preceding ICN. However, given that low *c* solutions are optimal according to our method and that we expect the dynamo simulations to produce anomalously high *c* (see above) we prefer the solutions in Fig. 12 corresponding to a mean decline in field strength before ICN.

All models in Figs 11 and 12 predict field strength for the Brunhes that is compatible with the Holocene field, but is generally at the upper end of the PINT range and cannot reproduce the lowest values in PINT even at the 3σ level. Part of the discrepancy can be explained by the inclusion in PINT of VDMs that may sample a transitional field. For many palaeomagnetic studies on more ancient rocks, it is often unclear whether palaeointensities are sampling a field of stable polarity or in a transitional state. In any case, considering the myriad factors
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Same as Fig. 11, but with *c* = 0.2.</figcaption></figure>

that influence the absolute field strength (discussed above) and the fact that the scaling prefactors are simply fit to simulation data we consider it a success of the overall approach that the theoretical predictions are so close to the observed values for the recent field.

While we do not attempt to fit the VDM low around 0.5 Ga, it is interesting to note that the predicted TDMs around this period vary strongly as a function of Δρ and *c*. For the values of τ favoured by the best-fitting models with the low *c* value (Fig. 12), ICN corresponds to a predicted TDM low around 0.4–1.0 Ga and so the predicted field strength at ~0.5 Ga depends strongly on whether the inner core has nucleated or not. For Δρ = 0.6 gm cc⁻¹, ICN occurs almost contemporaneously with the VDM low in PINT, but models with Δρ = 0.8 and 1.0 gm cc⁻¹ have ICN at earlier times and hence strongly overpredict the field strength at 0.5 Ga. For the high *c* values (Fig. 11), ICN corresponds to a TDM low with high Δρ, while the TDM is basically flat using the lower Δρ values. Following ICN all models predict a steep TDM increase that is not seen in PINT. Indeed the predictions fail to match the PINT bin at ~200 Ma even at the 3σ level.

Figs 11 and 12 clearly mark out a critical period between 400 and 1000 Ma characterized by a relative paucity of palaeointensity data and significant predicted changes in TDM. The large data gap may simply reflect challenges inherent in recovering robust magnetic recorders. With some recent exceptions (e.g. Hawkins *et al.* 2019; Bono *et al.* 2019) the majority of published data in this interval were measured using techniques that cannot detect secondary alteration or the presence of multidomain magnetic carriers, or have been shown to be biased by low unblocking temperatures. Alternatively, intervals of sparse palaeointensity data may reflect the existence of multipolar or dominantly non-dipolar fields (Abrajevitch & Van der Voo 2010; Driscoll 2016; Hawkins *et al.* 2019). In this case, the theoretical TDM would clearly be erroneous since it is derived assuming dipole dominance. Even if the field remained dipole-dominated the simple imposed CMB heat flows used to predict TDM do not capture the rapid dynamical variations seen in global mantle circulation models (e.g. Nakagawa & Tackley 2014) or long-term modulations such as supercontinent cyclicity, which has been suggested to affect the palaeomagnetic record during the Phanerozoic (e.g. Hounslow *et al.* 2018). Landeau *et al.* (2017) suggested an alternative ‘uniformitarian’ scenario in which the dipole field exhibits no significant changes through ICN and declines in strength as the inner core grows. However, this interpretation is not consistent
with the PINT data set, which shows a long-timescale decline in field strength from a high field at the end of the Archean to a dipole field minimum in the Ediacaran (Biggin *et al.* 2015; Bono *et al.* 2019) and, on average, an increase in field strength from post-ICN to present day. The scaling laws predict that the minimum TDM and maximum change in TDM should occur around ICN, which can hopefully be tested with new palaeomagnetic acquisitions. Improved constraints from seismology on the ICB density jump are also crucial for narrowing down the window of inner core formation and hence the low in VDM.

The main conclusions of this study are:

(i) The RMS and dipole CMB field follow scaling behaviour predicted by QG-MAC theory.

(ii) In order to reveal the scaling behaviour of the CMB field it is vital to filter out simulations with a low magnetic to kinetic energy ratio.

(iii) The QG-MAC-free scaling theory of Davidson (2013) yields field strength predictions that are compatible with a suite of 225 geodynamo simulations and both the modern and palaeomagnetic field strength. By contrast the QG-MAC-fixed theory (Starchenko & Jones 2002) overpredicts both the modern and palaeo CMB field. These results further support the application of QG-MAC-free theory to Earth's core dynamics.

(iv) Extrapolating to Earth's core conditions using the QG-MAC-free scaling suggests that the present RMS internal field strength is less than 10 mT (Fig. 3).

(v) For models with a CMB heat flow decay time $τ < 16$ Gyr, ICN corresponds to the lowest TDM value in the last 4.5 Gyr assuming a dipole-dominated field, while for $τ ≥ 16$ Gyr the TDM minimum occurs at 4.5 Gyr.

(vi) TDMs that best fit PINT have $τ ≤ 16$ Gyr and correspond to present-day CMB heat flow of 12–16 TW, increasing to 17–22 TW at 4 Gyr.

(vii) Best-fitting TDMs reproduce binned PINT VDMs before ICN within 1 standard deviation, but PINT does not contain the predicted strong values post ICN.

## ACKNOWLEDGEMENTS

CD acknowledges a Natural Environment Research Council (NERC) personal fellowship, reference NE/L011328/1. CD and AB acknowledge support from NERC grant NE/P00170X/1. RKB, DGM and AB acknowledge support provided by The Leverhulme Trust Research Leadership Award, RL-2016-080. RKB also acknowledges support provided by The Leverhulme Trust Early Career Fellowship, ECF-2020-617. SG acknowledges funding from the NERC SPHERES Doctoral Training Program. We are very grateful to Prof Uli Christensen and Dr Tobias Schwaiger for kindly sharing their dynamo simulation data sets. We also thank Prof Mathieu Dumberry for a constructive and insightful review of our work. CD is grateful to the initial explorations in the undergraduate dissertation projects of Bethany Shenton and Marissa Gonzalez, which provided the impetus for this work.

## DATA AVAILABILITY STATEMENT

Data tables and code are available at https://github.com/scs1cd/Bscaling.

## REFERENCES

Abrajevitch, A. & Van der Voo, R., 2010. Incompatible Ediacaran paleomagnetic directions suggest an equatorial geomagnetic dipole hypothesis, *Earth planet. Sci. Lett.*, **293**(1), 164–170.

Alfè, D., Gillan, M. & Price, G., 2002. Composition and temperature of the Earth's core constrained by combining *ab initio* calculations and seismic data, *Earth planet. Sci. Lett.*, **195**, 91–98.

Anufriev, A., Jones, C. & Soward, A., 2005. The Boussinesq and anelastic liquid approximations for convection in Earth's core, *Phys. Earth planet. Int.*, **152**, 163–190.

Aubert, J., 2019. Approaching Earth's core conditions in high-resolution geodynamo simulations, *Geophys. J. Int.*, **219**, S137–S151.

Aubert, J. & Gillet, N., 2021. The interplay of fast waves and slow convection in geodynamo simulations nearing Earth's core conditions, *Geophys. J. Int.*, **225**(3), 1854–1873.

Aubert, J., Labrosse, S. & Poitou, C., 2009. Modelling the palaeo-evolution of the geodynamo, *Geophys. J. Int.*, **179**, 1414–1428.

Aubert, J., Finlay, C. & Fournier, A., 2013. Bottom-up control of geomagnetic secular variation by the Earth's inner core, *Nature*, **502**, 219–223.

Aubert, J., Gastine, T. & Fournier, A., 2017. Spherical convective dynamos in the rapidly rotating asymptotic regime, *J. Fluid Mech.*, **813**, 558–593.

Aurnou, J. & King, E., 2017. The cross-over to magnetostrophic convection in planetary dynamo systems, *Phil. Trans. R. Soc. Lond. A*, **473**(2199), 20160731.

Biggin, A., Strik, G. & Langereis, C., 2008. Evidence for a very-long-term trend in geomagnetic secular variation, *Nat. Geosci.*, **1**, 395–398.

Biggin, A., Strik, G. & Langereis, C., 2009. The intensity of the geomagnetic field in the late-Archaean: new measurements and an analysis of the updated IAGA palaeointensity database, *Earth Planets Space*, **61**, 9–22.

Biggin, A., Piispa, E., Pesonen, L., Holme, R., Paterson, G., Veikkolainen, T. & Tauxe, L., 2015. Palaeomagnetic field intensity variations suggest Mesoproterozoic inner-core nucleation, *Nature*, **526**(7572), 245–248.

Biggin, A.J. & Paterson, G.A., 2014. A new set of qualitative reliability criteria to aid inferences on palaeomagnetic dipole moment variations through geological time, *Front. Earth Sci.*, **2**, doi:10.3389/feart.2014.00024.

Biggin, A.J., Bono, R.K., Meduri, D.G., Sprain, C.J., Davies, C.J., Holme, R. & Doubrovine, P.V., 2020. Quantitative estimates of average geomagnetic axial dipole dominance in deep geological time, *Nat. Commun.*, **11**, 6100, doi:10.1038/s41467-020-19794-7.

Bono, R.K., Tarduno, J.A., Nimmo, F. & Cottrell, R.D., 2019. Young inner core inferred from Ediacaran ultra-low geomagnetic field intensity, *Nat. Geosci.*, **12**(2), 143–147.

Buffett, B., 2010. Tidal dissipation and the strength of the Earth's internal magnetic field, *Nature*, **468**, 952–954.

Cheng, J.S. & Aurnou, J.M., 2016. Tests of diffusion-free scaling behaviors in numerical dynamo datasets, *Earth planet. Sci. Lett.*, **436**, 121–129.

Christensen, U., 2010. Dynamo scaling laws and applications to planets, *Space Sci. Rev.*, **152**, 565–590.
Christensen, U. & Aubert, J., 2006. Scaling properties of convection-driven dynamos in rotating spherical shells and application to planetary magnetic fields, *Geophys. J. Int.*, **166**, 97–114.

Christensen, U. & Tilgner, A., 2004. Power requirement of the geodynamo from ohmic losses in numerical and laboratory dynamos, *Nature*, **439**, 169–171.

Christensen, U., Aubert, J. & Hulot, G., 2010. Conditions for Earth-like geodynamo models, *Earth planet. Sci. Lett.*, **296**, 487–496.

Constable, C., Korte, M. & Panovska, S., 2016. Persistent high paleosecular variation activity in southern hemisphere for at least 10 000 years, *Earth planet. Sci. Lett.*, **453**, 78–86.

Davidson, P., 2013. Scaling laws for planetary dynamos, *Geophys. J. Int.*, **195**(1), 67–74.

Davies, C., 2015. Cooling history of Earth's core with high thermal conductivity, *Phys. Earth planet. Int.*, **247**, 65–79.

Davies, C. & Constable, C., 2014. Insights from geodynamo simulations into long-term geomagnetic field behaviour, *Earth planet. Sci. Lett.*, **404**, 238–249.

Davies, C. & Gubbins, D., 2011. A buoyancy profile for the Earth's core, *Geophys. J. Int.*, **187**, 549–563.

Davies, C., Pozzo, M., Gubbins, D. & Alfè, D., 2015. Constraints from material properties on the dynamics and evolution of Earth's core, *Nat. Geosci.*, **8**, 678–687.

de Koker, N., Steinle-Neumann, G. & Vojtech, V., 2012. Electrical resistivity and thermal conductivity of liquid Fe alloys at high P and T and heat flux in Earth's core, *Proc. Natl. Acad. Sci.*, **109**, 4070–4073.

Dekkers, M.J. & Böhnel, H.N., 2006. Reliable absolute palaeointensities independent of magnetic domain state, *Earth planet. Sci. Lett.*, **248**(1), 508–517.

Di Chiara, A., Muxworthy, A.R., Trindade, R. I.F. & Bispo-Santos, F., 2017. Paleoproterozoic geomagnetic field strength from the avanavero mafic sills, amazonian craton, Brazil, *Geochem. Geophys. Geosys.*, **18**(11), 3891–3903.

Driscoll, P., 2016. Simulating 2 Ga of geodynamo history, *Geophys. Res. Lett.*, **43**(11), 5680–5687.

Driscoll, P. & Bercovici, D., 2014. On the thermal and magnetic histories of Earth and Venus: influences of melting, radioactivity, and conductivity, *Phys. Earth planet. Int.*, **236**, 36–51.

Driscoll, P. & Wilson, C., 2018. Paleomagnetic biases inferred from numerical dynamos and the search for geodynamo evolution, *Front. Earth Sci.*, **6**, 113.

Finlay, C., Olsen, N., Kotsiaros, S., Gillet, N. & Tøffner-Clausen, L., 2016. Recent geomagnetic secular variation from Swarm and ground observatories as estimated in the CHAOS-6 geomagnetic field model, *Earth Planets Space*, **68**(1), 1–18.

Gastine, T., Wicht, J. & Aubert, J., 2016. Scaling regimes in spherical shell rotating convection, *J. Fluid Mech.*, **808**, 690–732.

Gillet, N., Jault, D., Canet, E. & Fournier, A., 2010. Fast torsional waves and strong magnetic field within the Earth's core, *Nature*, **465**, 74–77.

Gomi, H., Ohta, K., Hirose, K., Labrosse, S., Caracas, R., Verstraete, V. & Hernlund, J., 2013. The high conductivity of iron and thermal evolution of the Earth's core, *Phys. Earth planet. Int.*, **224**, 88–103.

Greenwood, S., Davies, C. & Mound, J., 2021. On the evolution of thermally stratified layers at the top of Earth's core, *Phys. Earth planet. Int.*, **318**, 106763.

Gubbins, D., Alfè, D., Masters, G., Price, G. & Gillan, M., 2003. Can the Earth's dynamo run on heat alone?, *Geophys. J. Int.*, **155**, 609–622.

Gubbins, D., Alfè, D., Masters, G., Price, G. & Gillan, M., 2004. Gross thermodynamics of two-component core convection, *Geophys. J. Int.*, **157**, 1407–1414.

Gubbins, D., Alfè, D., Davies, C. & Pozzo, M., 2015. On core convection and the geodynamo: effects of high electrical and thermal conductivity, *Phys. Earth planet. Int.*, **247**, 56–64.

Hawkins, L., Anwar, T., Shcherbakova, V., Biggin, A., Kravchinsky, V., Shatsillo, A. & Pavlov, V., 2019. An exceptionally weak Devonian geomagnetic field recorded by the Viluy Traps, Siberia, *Earth planet. Sci. Lett.*, **506**, 134–145.

Herrero-Bervera, E., Krasa, D. & Van Kranendonk, M.J., 2016. A whole rock absolute paleointensity determination of dacites from the Duffer Formation (ca. 3.467 Ga) of the Pilbara Craton, Australia: An impossible task?, *Phys. Earth planet. Int.*, **258**, 51–62.

Hounslow, M.W., Domeier, M. & Biggin, A.J., 2018. Subduction flux modulates the geomagnetic polarity reversal rate, *Tectonophysics*, **742-743**, 34–49.

Inoue, H., Suehiro, S., Ohta, K., Hirose, K. & Ohishi, Y., 2020. Resistivity saturation of hcp Fe-Si alloys in an internally heated diamond anvil cell: a key to assessing the Earth's core conductivity, *Earth planet. Sci. Lett.*, **543**, 116357.

Jackson, A., Jonkers, A. & Walker, M., 2000. Four centuries of geomagnetic secular variation from historical records, *Phil. Trans. R. Soc. Lond. A*, **358**, 957–990.

Jones, C., 2015. Thermal and compositional convection in the outer core, in *Treatise on Geophysics Vol. 8, Ch.5*, ed. Schubert, G., pp. 116–159, Elsevier.

King, E. & Buffett, B., 2013. Flow speeds and length scales in geodynamo models: the role of viscosity, *Earth planet. Sci. Lett.*, **371–372**, 156–162.

Kodama, K.P., Carnes, L.K., Tarduno, J.A. & Berti, C., 2019. Palaeointensity of the 1.3 billion-yr-old Gardal basalts, southern Greenland revisited: no evidence for onset of inner core growth, *Geophys. J. Int.*, **217**(3), 1974–1987.

Konôpková, Z., McWilliams, R., Gómez-Pérez, N. & Goncharov, A., 2016. Direct measurement of thermal conductivity in solid iron at planetary core conditions, *Nature*, **534**, 99–101.

Kulakov, E. et al., 2019. Analysis of an updated paleointensity database (QPI-PINT) for 65–200 Ma: implications for the long-term history of dipole moment through the Mesozoic, *J. geophys. Res.*, **124**(10), 9999–10022.

Kulakov, E.V., Smirnov, A.V. & Diehl, J.F., 2013. Absolute geomagnetic paleointensity as recorded by 1.09 Ga Lake Shore Traps (Keweenaw Peninsula, Michigan), *Stud. Geophys. Geod.*, **57**(4), 565–584.

Labrosse, S., 2015. Thermal evolution of the core with a high thermal conductivity, *Phys. Earth planet. Int.*, **247**, 36–55.

Landeau, M., Aubert, J. & Olson, P., 2017. The signature of inner-core nucleation on the geodynamo, *Earth planet. Sci. Lett.*, **465**, 193–204.

Lhuillier, F., Hulot, G., Gallet, Y. & Schwaiger, T., 2019. Impact of inner-core size on the dipole field behaviour of numerical dynamo simulations, *Geophys. J. Int.*, **218**(1), 179–189.

Lister, J. & Buffett, B., 1998. Stratification of the outer core at the core-mantle boundary, *Phys. Earth planet. Int.*, **105**, 5–19.

Long, R., Mound, J., Davies, C. & Tobias, S., 2020. Scaling behaviour in spherical shell rotating convection with fixed-flux thermal boundary conditions, *J. Fluid Mech.*, **889**, A7.

Masters, G. & Gubbins, D., 2003. On the resolution of density within the Earth, *Phys. Earth planet. Int.*, **140**, 159–167.

McDermott, B. & Davidson, P., 2019. A physical conjecture for the dipolar-multipolar dynamo transition, *J. Fluid Mech.*, **874**, 995–1020.

Meduri, D.G., Biggin, A.J., Davies, C.J., Bono, R.K., Sprain, C.J. & Wicht, J., 2021. Numerical dynamo simulations reproduce paleomagnetic field behaviour, *Geophys. Res. Lett.*, **48**(5), e2020GL090544.

Moffatt, H., 1978. Magnetic field generation in electrically conducting fluids, *Cambridge Monographs on Mechanics and Applied Mathematics*, pp. Cambridge University Press.

Mound, J., Davies, C. & Silva, L., 2015. Inner core translation and the hemispheric balance of the geomagnetic field, *Earth planet. Sci. Lett.*, **424**, 148–157.

Nakagawa, T. & Tackley, P., 2014. Influence of combined primordial layering and recycled MORB on the coupled thermal evolution of Earth's mantle and core, *Geochem. Geophys. Geosys.*, **15**, 619–633.

Nimmo, F., 2015. Energetics of the core, in *Treatise on Geophysics*. 2nd Edn, Vol. 8, ed. Schubert, G., pp. 27–55, Elsevier, Amsterdam.

Nimmo, F., Price, G., Brodholt, J. & Gubbins, D., 2004. The influence of potassium on core and geodynamo evolution, *Geophys. J. Int.*, **156**, 363–376.
Ohta, K., Kuwayama, Y., Hirose, K., Shimizu, K. & Ohishi, Y., 2016. Experimental determination of the electrical resistivity of iron at Earth's core conditions, *Nature*, **534**(7605), 95–98.

O'Rourke, J., Korenaga, J. & Stevenson, D., 2017. Thermal evolution of Earth with magnesium precipitation in the core, *Earth planet. Sci. Lett.*, **458**, 263–272.

Oruba, L. & Dormy, E., 2014. Transition between viscous dipolar and inertial multipolar dynamos, *Geophys. Res. Lett.*, **41**(20), 7115–7120.

Pozzo, M., Davies, C., Gubbins, D. & Alfè, D., 2012. Thermal and electrical conductivity of iron at Earth's core conditions, *Nature*, **485**, 355–358.

Pozzo, M., Davies, C., Gubbins, D. & Alfè, D., 2013. Transport properties for liquid silicon-oxygen-iron mixtures at Earth's core conditions, *Phys. Rev. B*, **87**, 014110.

Schaeffer, N., Jault, D., Nataf, H.-C. & Fournier, A., 2017. Turbulent geodynamo simulations: a leap towards Earth's core, *Geophys. J. Int.*, **211**(1), 1–29.

Schwaiger, T., Gastine, T. & Aubert, J., 2019. Force balance in numerical geodynamo simulations: a systematic study, *Geophys. J. Int.*, **219**(Supplement.1), S101–S114.

Schwaiger, T., Gastine, T. & Aubert, J., 2021. Relating force balances and flow length scales in geodynamo simulations, *Geophys. J. Int.*, **224**(3), 1890–1904.

Shaw, J., 1974. A new method of determining the magnitude of the palaeomagnetic field: Application to five historic lavas and five archaeological samples, *Geophys. J. R. astr. Soc.*, **39**(1), 133–141.

Shcherbakova, V., Biggin, A., Veselovskiy, R., Shatsillo, A., Hawkins, L., Shcherbakov, V. & Zhidkov, G., 2017. Was the Devonian geomagnetic field dipolar or multipolar? Palaeointensity studies of Devonian igneous rocks from the Minusa Basin (Siberia) and the Kola Peninsula dykes, Russia, *Geophys. J. Int.*, **209**(2), 1265–1286.

Sheyko, A., Finlay, C., Favre, J. & Jackson, A., 2018. Scale separated low viscosity dynamos and dissipation within the Earth's core, *Sci Rep.*, **8**(1), 1–7.

Smirnov, A., Tarduno, J., Kulakov, E., McEnroe, S. & Bono, R., 2016. Palaeointensity, core thermal conductivity and the unknown age of the inner core, *Geophys. J. Int.*, **205**(2), 1190–1195.

Sprain, C.J., Swanson-Hysell, N.L., Fairchild, L.M. & Gaastra, K., 2018. A field like today's? the strength of the geomagnetic field 1.1 billion years ago, *Geophys. J. Int.*, **213**(3), 1969–1983.

Sprain, C.J., Biggin, A.J., Davies, C.J., Bono, R.K. & Meduri, D.G., 2019. An assessment of long duration geodynamo simulations using new paleomagnetic modeling criteria (QPM), *Earth planet. Sci. Lett.*, **526**, 1157–58.

Starchenko, S. & Jones, C., 2002. Typical velocities and magnetic field strengths in planetary interiors, *Icarus*, **157**, 426–435.

Tang, F. et al., 2019. Secondary magnetite in ancient zircon precludes analysis of a Hadean geodynamo, *Proc. Natl. Acad. Sci.*, **116**(2), 407–412.

Tarduno, J. et al., 2010. Geodynamo, solar wind, and magnetopause 3.4 to 3.45 billion years ago, *Science*, **327**, 1238–1240.

Tarduno, J. et al., 2020. Paleomagnetism indicates that primary magnetite in zircon records a strong Hadean geodynamo, *Proc. Natl. Acad. Sci.*, **117**(5), 2309–2318.

Tarduno, J.A., Cottrell, R.D., Davis, W.J., Nimmo, F. & Bono, R.K., 2015. A Hadean to Paleoarchean geodynamo recorded by single zircon crystals, *Science*, **349**, 521–524.

Tauxe, L. & Yamazaki, T., 2015. 5.13-paleointensities, in *Treatise on Geophysics*, pp. 461–509, Elsevier.

Thellier, E. & Thellier, O., 1959. Sur l'intensité du champ magnétique terrestre dans le passé historique et géologique, *Annal. Geophys.*, **15**, 285–376.

Tobias, S., 2021. The turbulent dynamo, *J. Fluid Mech.*, **912**, P1.

Usui, Y. & Tian, W., 2017. Paleomagnetic directional groups and paleointensity from the flood basalt in the tarim large igneous province: implications for eruption frequency, *Earth Planets Space*, **69**(1), 14–doi:10.1186/s40623-016-0595-x.

Valet, J.-P., Brassart, J., Meur, I.L., Soler, V., Quidelleur, X., Tric, E. & Gillot, P.-Y., 1996. Absolute paleointensity and magnetomineralogical changes, *J. geophys. Res.*, **101**, 25029–25044.

Veselovskiy, R.V. et al., 2019. 1.86 Ga key paleomagnetic pole from the Murmansk craton intrusions Eastern Murman Sill Province, NE Fennoscandia: multidisciplinary approach and paleotectonic applications, *Precamb. Res.*, **324**, 126–145.

Wicht, J. & Meduri, D., 2016. A Gaussian model for simulated geomagnetic field reversals, *Phys. Earth planet. Int.*, **259**, 45–60.

Wilson, R.L., 1961. The thermal demagnetization of natural magnetic moments in rocks, *Geophys. J. R. astr. Soc.*, **5**(1), 45–58.

Yadav, R., Gastine, T., Christensen, U., Wolk, S. & Poppenhaeger, K., 2016. Approaching a realistic force balance in geodynamo simulations, *Proc. Natl. Acad. Sci.*, **113**(43), 12065–12070.

Yamamoto, Y. & Tsunakawa, H., 2005. Geomagnetic field intensity during the last 5 Myr: LTD-DHT Shaw palaeointensities from volcanic rocks of the Society Islands, French Polynesia, *Geophys. J. Int.*, **162**(1), 79–114.

Yamamoto, Y., Tsunakawa, H. & Shibuya, H., 2003. Palaeointensity study of the Hawaiian 1960 lava: implications for possible causes of erroneously high intensities, *Geophys. J. Int.*, **153**(1), 263–276.

Zhang, Y., Hou, M., Liu, G., Zhang, C., Prakapenka, V.B., Greenberg, E., Fei, Y., Cohen, R. & Lin, J.-F., 2020. Reconciliation of experiments and theory on transport properties of iron and the geodynamo, *Phys. Rev. Lett.*, **125**(7), 078501.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

**Table S1.** Additional site level data added to PINT data set. Columns are: description, brief location name; SLAT, site latitude; SLONG, site longitude; AGE, age assigned to magnetization; *N*<sub>INT</sub>, number of samples comprising mean palaeointensity; F, site mean palaeointensity in µT; VDM/VADM, virtual (axial) dipole moment in 1 × 10<sup>22</sup> Am<sup>2</sup>. When VDM/VADM not available, VADM is estimated assuming equatorial site location; reference, published study reference.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.