Numerical Dynamo Simulations Reproduce Paleomagnetic Field Behavior

Numerical geodynamo simulations capture several features of the spatial and temporal geomagnetic field variability on historical and Holocene timescales. However, a recent analysis questioned the ability of these numerical models to comply with long‐term paleomagnetic field behavior. Analyzing a suite of 50 geodynamo models, we present here the first numerical simulations known to reproduce the salient aspects of the paleosecular variation and time‐averaged field behavior since 10 Ma. We find that the simulated field characteristics covary with the relative dipole field strength at the core‐mantle boundary (dipolarity). Only models dominantly driven by compositional convection, with an Ekman number (ratio of viscous to Coriolis forces) lower than 10−3 and a dipolarity in the range 0.34–0.56 can capture the observed paleomagnetic field behavior. This dipolarity range agrees well with state‐of‐the‐art statistical field models and represents a testable prediction for next generation global paleomagnetic field model reconstructions.

• We present the first numerical geodynamo simulations known to reproduce the main features of paleomagnetic field variability since 10 Ma • All simulated characteristics of paleomagnetic behavior covary with the degree of dipole dominance (dipolarity) • Only chemically driven dynamos at sufficiently low Ekman numbers in a specific dipolarity range capture paleomagnetic field behavior Supporting Information: • Supporting Information S1 • Table S1 • Table S5 10 Myr, the authors found that none of the 46 simulations explored could capture the main aspects of the observed variability. In fact, the large majority of the simulations performed poorly; only a few passed three out of the five Q PM criteria with large total misfits.
Here we present a new set of simulations reproducing paleomagnetic field behavior of the last 10 Myr according to the Q PM criteria. First, we show that the relative strength of the dipole field to the total field up to spherical harmonic (SH) degree and order 12 at the core-mantle boundary (CMB) can be used as a proxy for all paleomagnetic observables considered. We then examine the conditions for obtaining paleomagnetic-like simulations and discuss implications for the Earth's core.

Model Formulation
The setup and solution method for the geodynamo models are standard and extensively documented elsewhere (Davies & Constable, 2014;Wicht, 2002;Wicht & Meduri, 2016, WM16 hereafter;Willis et al., 2007). We therefore provide only a brief description here (see also Section S1). We consider a convection-driven magnetohydrodynamic flow under the Boussinesq approximation with the fluid confined to a spherical shell of thickness d = r o − r i rotating at a constant angular velocity Ω. Here r i and r o are the inner and outer boundary radii, which are identified with the inner core radius and the CMB radius, respectively.
All models assume no-slip mechanical boundary conditions and an electrically insulating mantle. We employ the codensity approach where density perturbations due to compositional and temperature differences are described by only one variable. Different convective driving scenarios are modeled via the boundary conditions and homogeneous volumetric codensity sinks. Thermal dynamos are bottom heated with either fixed flux or fixed temperature at r i . All the heat entering at r i leaves the system through r o where a fixed flux condition is imposed. Some models employ lateral variations in the outer boundary heat flux in the form of a recumbent SH of degree ℓ = 2 and order m = 0 as an approximation of the observed lower mantle seismic shear-wave structures (Dziewonski et al., 2010).
Chemical dynamos are driven by either a fixed light element concentration or concentration gradient at r i , which is balanced by a volumetric sink. The flux through r o is set to zero. While the chemical dynamos assume an electrically conducting inner core, the thermal dynamos use an insulating inner core for simplicity. Wicht (2002) showed that the impact of inner core conductivity on the magnetic field and its variability is minor in thermal dynamos, although this may depend on the details of the convective driving and mechanical boundary conditions employed (Dharmaraj & Stanley, 2012;Lhuillier et al., 2013).
The dimensionless parameters controlling the system are the Ekman number Ek, the Prandtl number Pr, the magnetic Prandtl number Pm, and the shell aspect ratio χ: Here ν, η, and κ are the kinematic viscosity, magnetic diffusivity, and thermal (or compositional) diffusivity of the fluid, respectively. The Rayleigh number controls the vigor of convection and is defined in Section S1. Ek varies between 3 × 10 −4 and 2 × 10 −3 , and Pm spans the range 3-10. These ranges are constrained by the need to perform long temporal integrations.
Our data set is summarized in Table S1 and consists of 50 simulations: 21 from S19, 7 from WM16, and 22 new runs. From S19, we excluded thermal dynamos which use specific buoyancy profiles , large amplitudes of the CMB heat flux anomalies (ε = 1.5; see Table S1 for the definition of ε), and low Rayleigh numbers (regime of locked dynamo action). We also excluded the cases at Pm = 20. All these simulations poorly comply with Earth, having total Q PM misfits larger than 5 and total Q PM scores of 3 at most (see Section 2.2). The new thermal runs complement the Rayleigh number range explored by S19 and additionally include cases at Ek = 3 × 10 −4 . The new chemical runs are similar to those of WM16 but focus on reversing dipolar solutions.

Paleomagnetic Criteria for Geodynamo Simulations
The Q PM framework is described in detail in S19 and we recall only the essentials here. S19 identified five paleomagnetic observables that describe the long-term paleosecular variation (PSV) and time-averaged field (TAF) behavior. The first two observables characterize the virtual geomagnetic pole (VGP) angular dispersion S by estimating its equatorial value and latitudinal dependence. They are the parameters a and b of the empirical quadratic fit with paleolatitude λ introduced by McFadden et al. (1988), The third Q PM observable is the maximum of the absolute inclination anomaly Here I is the Fisher mean inclination (Fisher, 1953) and I GAD is the inclination expected under a geocentric axial dipole field. The fourth observable, V%, is the ratio of the interquartile range to the median of the virtual dipole moment (VDM) distribution. The last observable is the relative transitional time τ T , defined as the fraction of time spent with an absolute true dipole latitude lower than 45°, which is associated with a criterion on the presence of reversals.
Using the most recent compilation of paleomagnetic directional data PSV10 (Cromwell et al., 2018) and the paleointensity database PINT (Biggin et al., 2009(Biggin et al., , 2015, S19 estimated values and uncertainties of these five observables for the last 10 Myr (see Table S2). The sum of normalized linear misfits between simulated and observed values for each Q PM observable is ΔQ PM . If the normalized misfit of a given observable is ≤1, the observed and simulated paleomagnetic characteristics overlap within the respective estimated uncertainties.
Together with the misfits, S19 defined binary scores. The score of a given Q PM observable is 1 if the normalized misfit in that observable is ≤1 and is zero otherwise. The total score Q PM , obtained by summing the single scores, thus ranges from 0 to 5. By definition, a paleomagnetic-like simulation with the maximum score Q PM = 5 has a total misfit ΔQ PM ≤ 5. Even when Q PM < 5, however, the total misfit can be smaller than 5. While a large Q PM signifies a good compliance with Earth, a large ΔQ PM means the opposite.

Evidence for Paleomagnetic-like Geodynamo Simulations
The magnetic fields obtained in geodynamo simulations are generally characterized by their degree of dipole dominance, which is often measured by the dipolarity D 12 , defined as the time-averaged ratio of the root mean square (RMS) dipole field strength to the total RMS field strength up to SH degree and order ℓ = m = 12 at the CMB (Christensen & Aubert, 2006). Multipolar solutions generally have D 12 ≲ 0.35, while dominantly dipolar ones like the present geomagnetic field have D 12 ≳ 0.7 (Christensen, 2010;Christensen & Aubert, 2006). Dipolar reversing solutions, that is, dynamos which are dipole dominated most of the time but occasionally undergo polarity reversals, occur in a narrow dipolarity range sandwiched between the stable dipolar and the multipolar regimes (Driscoll & Olson, 2009;Wicht & Tilgner, 2010;Wicht et al., 2015).
Figures 1a and 1b demonstrate that D 12 is a good proxy for the total misfit ΔQ PM . When the Rayleigh number increases (in the direction indicated by the arrows in the connected lines in Figure 1) the dipolarity systematically decreases together with ΔQ PM until D 12 ≈ 0.5. For smaller values of D 12 , ΔQ PM increases again and the simulations roughly describe parabolic paths in the D 12 -ΔQ PM plane. These paths show no apparent dependence on the codensity boundary conditions and on Pm (see also Figure S1), but depend strongly on the convective driving mode and on the Ekman number. While the thermal dynamos barely reach misfits of ΔQ PM ≈ 4 with a score Q PM = 2 (Figure 1b), the chemical runs show ΔQ PM as low as 2.7 with Q PM = 4 (Figure 1a). In fact, these latter runs come close to a score of Q PM = 5, having either a moderately low relative transitional time or an equatorial dispersion only a few degrees higher than Earth (Section 3.2.1). These paleomagnetic-like dynamos combine chemical driving with the lower Ek = 3 × 10 −4 . The chemical runs at Ek = 10 −3 barely reach ΔQ PM ≈ 5 with Q PM = 3.
Remarkably, the optimal field solutions, that is, those which yield the lowest ΔQ PM and the highest Q PM in each Rayleigh number track, lie in a well-defined D 12 range. The quadratic function well describes the behavior of the simulations in both types of convective forcing with a high coefficient of determination R 2 (gray curves in Figures 1a and 1b; see Table S3 for values of the regression coefficients and R 2 ). The minima of the quadratic fits occur at D 12 = 0.45 and D 12 = 0.48 for the chemical and thermal runs respectively. Small departures from these values cause a large increase of ΔQ PM and a decrease in Q PM . The values of D 12 where the quadratic fit of the chemical runs at Ek = 3 × 10 −4 intersects the threshold ΔQ PM = 5 below which Earth-like models are expected define the dipolarity interval δD 12 = [0.34, 0.56].
The dipolarity values of our paleomagnetic-like dynamos are compatible with estimates obtained for Earth from global paleomagnetic field model reconstructions. Since these models have spatial power spectra that are generally considered to be well resolved only for SH degrees ℓ ≤ 4 (Korte & Constable, 2008;Nilsson et al., 2014;Wardinski & Korte, 2008, see also Figure S2), we analyze the modified dipolarity D 4 which includes SH contributions up to ℓ = m = 4. Note, however, that even degrees ℓ ≤ 4 may suffer from spatial and MEDURI ET AL.
10.1029/2020GL090544 4 of 10  (a) and (c) present the GGP models TK03 (Tauxe & Kent, 2004), BCE19 (Brandt et al., 2020), and BB18  (error bars denote one standard deviation above and below the dipolarity values). Capital letters A-F mark the six simulation runs discussed in the main text (see Table S1 for additional information). Ek, Ekman number; Pm, magnetic Prandtl number; GGP, giant Gaussian process. temporal regularization (Hellio & Gillet, 2018;Sanchez et al., 2016) and the true D 4 values could be somewhat smaller. Figures 1c and 1d show that the behavior of D 4 can also be described by a simple quadratic dependence. The paleomagnetic-like dipolarity interval predicted by our chemical dynamos at Ek = 3 × 10 −4 is δD 4 = [0.57, 0.80].
GGF100k (Panovska et al., 2018), the longest global field reconstruction to date, spans the last 100 kyr and provides D 4 = 0.84 ± 0.08, in relatively good agreement with our numerical prediction (Figure 1c; see also Table S4). LSMOD.2 (Brown et al., 2018) has a lower value of D 4 = 0.74 since it deliberately models the field during the Laschamp and Mono Lake excursions. This lower estimate falls within δD 4 and is in excellent agreement with the value obtained for run B, our most paleomagnetic-like simulation (Figure 1c and Table S4).
For field reconstructions covering shorter time intervals, the Holocene CALS10k.1b (Korte et al., 2011) model provides D 4 = 0.92, and the historical gufm1 (Jackson et al., 2000) and IGRF-13 (Thébault et al., 2015) models give D 4 = 0.88 and 0.82, respectively. Such high values of D 4 are likely due to the short timespans sampled by these models. It is encouraging that the differences with our numerical predictions of D 4 reduce for the longer time averages obtained from GGF100k and LSMOD.2.
Estimates of D 12 for Earth on timescales of million years can be obtained from statistical field models based on a giant Gaussian process (GGP). Here we consider the GGP models TK03 (Tauxe & Kent, 2004), BCE19 (Brandt et al., 2020) and BB18 , which are explicitly constructed to reproduce the PSV over the last 5-10 Myr, with BB18 also capturing the observed VDM distribution. TK03 and BCE19 differ only in the assumed variances of their independent and normally distributed Gauss coefficients, while BB18 additionally employs a covariance pattern for degrees ℓ ≤ 4 inferred from dynamo simulations. These GGP models have 0.54 ≤ D 12 ≤ 0.56 and thus fall within the paleomagnetic-like interval δD 12 predicted here (see Figure 1a and also Table S4). Figure 2 presents the five Q PM observables as a function of D 12 for all chemical runs. All observables increase as D 12 decreases, a trend which is also observed for the thermal dynamos ( Figure S3). Though our three most paleomagnetic-like simulations (runs B-D) only reach a total score of 4 (Figure 1a), they still reproduce the single missed Q PM observable to a reasonable level.

Variation of the Paleomagnetic Observables With D 12
Run B closely captures all observables except the relative transitional time τ T (Figure 2; Table S4), which is too small at 0.018, about half the Earth's lower bound value (Table S2). While this run showed one reversal and three excursions in 35 magnetic diffusion times, we cannot exclude it may yield Earth-like τ T when a robust statistic is obtained for longer integrations. Runs C and D have Earth-like transitional times (τ T = 0.046 and 0.065 respectively; Figure 2e) but a median equatorial dispersion a about 4° and 6° too high with misfits <1.6 (Figure 2a; Tables S2 and S4). Relaxing the uncertainties on a by <2° would yield a total score of 5 for these two runs. On this basis, we consider the simulated paleomagnetic behavior of runs B-D to be an excellent approximation to that of Earth in the last 10 Myr, while acknowledging it does not quite meet the full requirements for being classified as "Earth-like" by the current Q PM criteria.

Influence of the Ekman Number in Chemical Dynamos
As well as intermediate values of D 12 , chemical dynamos can reach low misfits and high scores only if the Ekman number is low enough (Section 3.1). This dependency on Ek results from the behavior of a and τ T . In the paleomagnetic-like dipolarity interval δD 12 , chemical dynamos at Ek = 10 −3 have comparable or higher a and lower τ T than the cases at Ek = 3 × 10 −4 (Figures 2a and 2e). Misfits in a and τ T are up to three times smaller in these low-Ek runs compared to the high-Ek cases, while misfits in the other Q PM observables are similar (for example, in Table S4 compare run D with run E, the simulation with the lowest ΔQ PM at Ek = 10 −3 ).
In the simulations at Ek = 10 −3 , more frequent polarity transitions leading to Earth-like values of τ T start at D 12 ≈ 0.3 where a is already far too high (Figure 3a; white symbols show Earth-like τ T ). At the lower Ek of 3 × 10 −4 , on the other hand, reversals start to appear at D 12 ≈ 0.45 where a is still relatively Earth-like. Such a dependency on Ek for the onset of reversals arises because the dipole field variability, measured by the relative standard deviation  d d / B (the ratio of the standard deviation of the dipole field strength at the CMB to its time-averaged value), increases with decreasing Ek in our simulations (Figure 3b). These larger dipole fluctuations naturally lead to an increased likelihood of both transitional periods and polarity reversals (Driscoll & Olson, 2009;Meduri & Wicht, 2016, WM16). We note that a remains Earth-like in the low-Ek runs at D 12 ≳ 0.45 since the equatorially symmetric CMB field, which determines a (Coe & MEDURI ET AL.   Table S2 for further details). Dark gray shaded regions highlight the predicted paleomagnetic-like dipolarity interval δD 12 . Glatzmaier, 2006;McFadden et al., 1988), is weaker than in the high-Ek cases in the same dipolarity range ( Figure S4).

Influence of the Convective Driving Mode
Thermal dynamos at Ek = 3 × 10 −4 do not reach ΔQ PM as low and Q PM as high as the chemical dynamos at the same Ekman number (Section 3.1). In these runs, the Q PM observables vary similarly with D 12 , with the exception of the latitudinal VGP dispersion b and τ T ( Figure S3). The thermal runs, unlike the chemical ones, present values of b which remain low and non-Earth-like for all D 12 explored (Figure 3c). The weak variation of b with D 12 in the thermal runs can be attributed to the almost unchanged odd and even CMB field contributions ( Figure S4). The chemical runs, unlike the thermal cases, already show reversals at an intermediate D 12 . At D 12 = 0.53, for example, run F (thermal) shows a stable dipolar solution with τ T = 0, whereas run A (chemical) undergoes few excursions so that τ T = 0.01 (Figure 3c; see also Figure S3e). Chemical dynamos then reach Earth-like τ T at slightly lower values of D 12 where also the other Q PM observables are captured.

Discussion and Conclusions
We tested whether a new set of numerical dynamo simulations reproduces the paleomagnetic field behavior of the last 10 Myr by applying the Q PM criteria of Sprain et al. (2019). These criteria examine the equatorial and latitudinal VGP dispersion, the inclination anomaly, the VDM distribution, and the relative time spent by the true dipole pole in transitional latitudes, along with the presence of reversals.
We reported the first numerical simulations known to reproduce these fundamental characteristics of the paleomagnetic field since 10 Ma. The dipolarity D 12 , which measures the degree of dipole dominance at the CMB, appears to be a good proxy for all five Q PM observables across a variety of simulations differing in control parameters, boundary conditions and convective driving mode, and it allows predictions of the total Q PM misfit and score. Simulations capturing the observed field behavior are characterized by (1) a compositional driving, (2) an Ekman number Ek below 10 −3 , and (3) a dipolarity D 12 in the interval δD 12 = [0.34, 0.56]. Previous numerical studies exploring long-term geomagnetic field behavior do not generally employ simulations that fulfill all such conditions; for example, they often use Ek ≳ 10 −3 due to computational reasons.
Our best performing simulations employ a setup where buoyancy is released at the inner core boundary and absorbed by the outer core. This seems an appropriate scenario for the geodynamo after the inner core started crystallizing and the light elements, emanated from the growing inner core front, may have dominated convective driving (Labrosse, 2015;Nimmo, 2015). Taken at face value, our analysis appears to favor compositional over thermal convection as the dominant driving mode of the geodynamo over the  last 10 Myr, in agreement with both thermal history calculations (see, e.g., Nimmo, 2015) and numerical dynamo studies exploring the influence of the two convective drivings on the magnetic field morphology (Kutzner & Christensen, 2000).
The dipolarity interval δD 12 where paleomagnetic-like simulations are found is bounded above by the modern field. Current GGP models provide an estimate of D 12 for Earth of about 0.55, which falls within δD 12 and confirms the robustness of our numerical results. The suggested dipolarity interval represents a specific, testable prediction for next generation paleomagnetic field models, once they reach higher spatial resolutions and cover longer time intervals than the models currently available. Earth's core may lie at the transition between the dipolar and the multipolar dynamo regimes (Christensen & Aubert, 2006;Oruba & Dormy, 2014) and our results are compatible with this finding. This transition indeed occurs at D 12 in the range 0.35-0.50, which is included in the paleomagnetic-like interval δD 12 predicted here.
Our results suggest the possibility of constructing a path toward Earth's core conditions which preserves paleomagnetic-like dynamo characteristics. According to Oruba and Dormy (2014), the parameter combination Re Ek 2/3 (Re = Ud/ν is the Reynolds number, where U is the time-averaged RMS core flow velocity) defines the dipolar-multipolar transition at D 12 ≈ 0.5, which is close to the optimal value D 12 = 0.45 inferred from our analysis. To maintain D 12 constant while reducing Ek, Re needs to increase and this can be achieved by increasing the Rayleigh number and by decreasing Pm. Following this path to lower Ek and Pm and higher Rayleigh numbers, the relevant balance for the Earth's core between magnetic, Coriolis and buoyancy forces is expected to emerge naturally, as also suggested by recent high-resolution numerical simulations (Aubert, 2019;Aubert et al., 2017;Schaeffer et al., 2017;Yadav et al., 2016).
It has recently been argued that the geomagnetic field displayed similar average degrees of surface axial dipole dominance over large swathes of geological time . Insofar as dipolarity and surface axial dipole dominance may be assumed to be related, we expect that many of the conclusions reached here are also valid at certain earlier times in Earth's history.

Data Availability Statement
The simulation output data used to produce Figures 1-3 can be found in the supplementary information and can be accessed via DOI 10.17605/OSF.IO/PY85T.