Paleomagnetic Biases Inferred from Numerical Dynamos and the Search for Geodynamo Evolution

2 The orientation and intensity of the paleomagnetic ﬁeld is central to our understanding of 3 the history of the Earth. The paleomagnetic signature of the singular most event, inner core 4 nucleation, however, remains elusive. In this study we study numerical dynamo simulations from 5 a paleomagnetic perspective to explore how long observations must be time-averaged to obtain 6 stable virtual geomagnetic pole (VGP) directions and global ﬁeld intensities. We ﬁnd that running 7 averages over 20 − 40 kyr are needed to obtain stable VGP’s with α 95 < 10 ◦ , and over 40 − 120 8 kyr for α 95 < 5 ◦ . We ﬁnd that models with higher heat ﬂux and more frequent polarity reversals 9 require longer time averages, and that obtaining stable intensities requires longer time averaging 10 than obtaining stable directions. Running averages of local ﬁeld intensity and inclination produce 11 underestimates of VDM by factors of 0 . 9 − 0 . 6 and overestimates of VADM by factors of 1 − 1 . 2 as 12 heat ﬂux and reversal frequency increases. We derive a scaling law connecting reversal frequency 13 to paleointensity bias that could be applied to records where reversal frequency is known. Applied 14 to the PINT paleointensity record, these biases produce little change to the overall trend of a 15 relatively ﬂat but scattered intensity over the last 2 Ga. A more careful debiasing applied during 16 periods when the reversal frequency is known could reveal previously obscured features in the 17 paleointensity record. 18


INTRODUCTION
Our knowledge of the history of Earth's magnetic field derives from paleomagnetic signals preserved in

122
The dynamo models are computed using the Rayleigh dynamo code (Featherstone and Hindman, 2016;123 Matsui et al., 2016). All models share the following control parameters: E = 10 −3 , Ra = 10 6 , P r = 1, 124 P m = 10, insulating magnetic boundary conditions, no-slip velocity conditions, inner-outer core radius 125 ratio of 0.35, and fixed temperature gradient at both boundaries (see Table 1). These relatively high Ekman 126 number simulations produce Earth-like large scale magnetic features (see below) and polarity reversals that 127 resemble geomagnetic observations, and are numerically cheap so they can be run extremely long times to 128 produce low frequency statistics. The inner boundary temperature gradient dT /dr i , fixed in time in each 129 model, spans a range of 0.8 − 12 (in non-dimensional units) that produces stable dipolar dynamos on the 130 lower end to unstable, reversing non-dipolar dynamos at the high end (Table 1). The temperature gradient 131 at the outer boundary is set to balance the heat flow at the base so that energy is conserved and there is no 132 internal sink or source. 133 Time is scaled from thermal diffusion times (implemented in the code) to years by multiplying by a factor 134 τ dip /(P m(r o /π) 2 ), where r o = 1.5384 is dimensionless outer core radius and τ dip = 50 kyr is the assumed 135 magnetic dipole decay time of the core. We adopt the same definition of a polarity reversal proposed by 136 Driscoll and Olson (2009): that the dipole co-latitude θ dip spend at least 20 kyr in a stable polarity before 137 and after a reversal.

139
From each model we compute Gauss coefficients over time g l,m (t) and h l,m (t) at Earth's surface from 140 magnetic field spectra at the CMB (e.g. Merrill et al., 1996). Although the dynamo model spectra are 141 resolved out to harmonic degree l max = 64 we only compute Gauss coefficients out to l max = 8 because 142 larger harmonics contribute very little to the surface magnetic field. 143 From the Gauss coefficients we compute the rms non-axial dipole (NAD), (1) This is a provisional file, not the final typeset article and the rms dipole intensity, g D (t) = g 1,0 (t) 2 + g 1,1 (t) 2 + h 1,1 (t) 2 (2) We also compute local magnetic field quantities on Earth's surface that are interpreted as synthetic observations. From the Gauss coefficients we compute the surface vector magnetic field components X, Y , 147 and Z from the magnetic potential Ψ(r, θ, φ) at a point on the surface r = a, and P m l are Schmidt normalized Legendre polynomials (Merrill et al., 1996). From the local field 150 components we compute the local magnetic field intensity F as We will focus in particular on synthetic paleomagnetic observations at an arbitrary point on the equator: 152 θ = π/2 and φ = 0. We will refer to intensity at this location F (θ = π/2, φ = 0) as F eq .

153
Synthetic observations of magnetic direction are also created at the same point on the equator. These inclination I, angular pole distribution α 95 with 95% probability, and Fisher's precision parameter k and the directional cosines are l i = cos D i cos I i , m i = sin D i cos I i , and n i = sin I i (Merrill et al., 1996).

160
The (dimensionless) global dipole moment p DM is where a = 2.8157 is the dimensionless radius of Earth (Merrill et al., 1996). The global axial dipole 162 moment p ADM is the same as in (11) but with only the g 1,0 term in (2). These global dipole moments 163 p DM and p ADM can be considered as the "true" values and will be compared to synthetic or "virtual" 164 observations of these quantities. The "virtual" dipole moment (VDM) is computed from the local magnetic 165 intensity F by (Merrill et al., 1996) 166 and similarly for the "virtual" axial dipole moment (VADM) which reduces to p V ADM = 4πa 3 F at the equator (θ = π/2). 168 We also consider the role of time averaging in producing a paleomagnetic direction or intensity. The  (2014): where x(t) is some output from the dynamo model and t i is time at the i th sampling index within each

RESULTS
We apply the analysis methods described above to the suite of dynamo simulations to investigate how long

185
To demonstrate that these models are in a relevant region of parameter space, we first focus the details 186 of an "Earth-like" model with dT /dr i = 4, which we will refer to as "model 4". These magnetic field statistics give us some confidence that our models produce magnetic fields that are generally "Earth-like" at the largest scales even though they are many orders of magnitude from the Earth   Figure 4 shows the mean and standard deviation of the running mean of four output quantities (i.e. x(τ )) 210 from model 4 for a range of τ . Figure 4a shows that the average of g 1,0 (τ ) (i.e. the average of the running 211 averages from all sub series) converges to zero for all τ , as expected for this reversing model that spends 212 roughly equal time in normal and reversed polarity states ( Figure 2). The relatively large standard deviation 213 of g 1,0 (τ ) reflects the long-term variability that is not averaged out within each sub series. Figure 4b shows 214 that the average rms non-axial dipole field (GNAD) from (1) is also near zero for all τ , implying that all 215 other field harmonics (other than g 1,0 ) individually balance to zero. The standard deviation of GNAD also 216 approaches zero at large τ , implying that non-axial dipole fields more consistently balance out over longer 217 time averaging than g 1,0 . Figure 4c shows that the running mean of rms dipole from (2) and axial dipole 218 coefficients are stationary and non-zero over all τ . Figure 4d shows that the local magnetic amplitude F eq at 219 the equator is similarly stationary and non-zero over all τ and similar in amplitude to the rms axial dipole.

220
A more detailed comparison between observed and true intensities is below. 221 Figure 5 shows the average of running means of four local magnetic pole-related quantities from model 222 4: declination D eq from (6), inclination I eq from (7), α 95 from (8)  The major dynamo transition from dipolar non-reversing models to non-dipolar reversing models occurs 236 around dT /dr i = 3. This transition is apparent in Figure 6a where volume averaged magnetic energy (ME) 237 drops below kinetic energy (KE) due to a weakening of the axial dipole, Figure 6b where rms g 1,0 drops 238 by a factor of ∼ 4, Figure 6c where reversals begin, and Figure 6d where the axial dipolarity (g 1,0 /g rms ) 239 drops below ∼ 0.5.

240
Interestingly, Figure 6a shows that volume averaged ME drops to a minimum at the onset of reversals  non-reversing to reversing models as seen in Figure 6a.

267
Bias in intensity observations can be investigated by comparing the ratio of observed dipole moment 268 p V DM in (12) to true dipole moment p DM in (11), and similarly for the ratio of observed axial dipole moment p V ADM in (13) to true axial dipole moment This is a provisional file, not the final typeset article is systematically low compared to true DM by a factor of ∼ 0.9 for non-reversing dipolar models, and 273 trends down to a factor of ∼ 0.6 for reversing models. The drop in Figure 8c is caused mainly by a drop in 274 F/g D because the running average inclination at τ = 50 kyr is less than 2 • for all models (Figure 5b). The 275 decrease in F/g D is caused by the local intensity F becoming more and more contaminated by non-dipole 276 field components with arbitrary sign as heat flux increases, so that on average F < g D . Figure 8c shows 277 this ratio is the same for τ = 50 kyr, 100 kyr, and 500 kyr, indicating that this bias lingers even with longer 278 time averaging.

279
On the contrary we find an increasing trend in the ratio of VADM/ADM from 1 to ∼ 1.2, which is 280 caused by g 2 1,0 decreasing faster than F with increasing heat flux. This bias can be attributed to the GAD 281 assumption that the field is purely axial when it is not and non-axial dipolar field components contributing 282 to increase F . This trend implies that VADM systematically overestimates the true DM for reversing 283 dynamos. Combining F < g D from Figure 8c and F > g 2 1,0 from Figure 8d gives 284 g 2 1,0 < F < g 2 1,0 + g 2 1,1 + h 2 1,1 implying that local fields of opposite sign combine to decrease F slightly less than g D but to increase it 285 slightly more than g 2 1,0 .

IMPLICATIONS
The intensity biases found above are now applied to the PINT data in Figure 1. In an attempt to remove  is then applied to the unbiased data in Figure 10. This relatively minor biasing effect does not reveal any 297 new long-term trends and no immediate signature of inner core nucleation is apparent.

CONCLUSIONS
We have generated synthetic magnetic observations from numerical dynamos that span a range of dynamical 308 regimes to investigate the time averaged magnetic field orientation and intensity. The range of dynamo 309 regimes found, from stable non-reversing dipolar regimes to reversing non-dipolar regimes, are driven by 310 models that span a factor of 10 increase in bottom boundary heat flux. We find that running averages over 311 20 − 40 kyr are needed to obtain stable VGP's with α 95 < 10 • , and over 40 − 120 kyr for α 95 < 5 • . To 312 obtain stable VGP's we find that models with higher heat flux and more frequent polarity reversals require 313 longer time averages, similar to previous studies (McMillan et al., 2001;Davies and Constable, 2014).

314
However, we also find that obtaining stable intensities requires longer time averaging than directions.
315 Surprisingly we find that running averages of local field intensity and inclination produce underestimates for time averages over 500 kyr. We compute a scaling law connecting reversal frequency to paleointensity 320 bias that could be applied to records where reversal frequency is known.

321
These biases are applied to the PINT paleointensity record, which produce little change to the overall 322 trend of a relatively flat intensity over the last 2 Ga. A more careful debiasing could be applied during 323 periods when the reversal frequency, or some proxy for reversal frequency (such as secular variation or 324 polarity ratio), is known. Correcting for this bias may be important for identifying trends and events (like 325 ICN) in the paleointensity record.

326
Future analysis could be extended to synthetic observations at difference locations on the surface, secular 327 variation, and identifying higher order magnetic components. Models at lower Ekman number and with 328 different boundary conditions should also be analyzed to investigate if these results are sensitive to this 329 region of parameter space.

CONFLICT OF INTEREST STATEMENT
The authors declare that the research was conducted in the absence of any commercial or financial 331 relationships that could be construed as a potential conflict of interest.

AUTHOR CONTRIBUTIONS
The author conceived the project, analyzed the data, and wrote the article.

FUNDING
The author thanks the Carnegie Institution for Science for funding.    Figure 1. Paleointensity (VDM and VADM) from the PINT database (dots) with inverse-distance squared smoothing (solid line) applied to filled circles. Open circles are from paleointensity studies that either contain directional anomalies (Pesonen and Halls, 1983) or low temperature components (Thomas, 1993;Thomas and Piper, 1995 (7). c) Angular spread of magnetic pole location at 95% confidence level, α 95 from (8). d) Fisher's precision parameter k from (9).   Figure 8. Comparison of synthetic observations made at the equator for all models, with τ = 50 kyr. a) True dipole moment p DM from (11) and "observed" p V DM from (12). b) True axial dipole moment p ADM and "observed" p V ADM from (13). c) Ratio of p V DM /p DM from (a). d) Ratio of p V ADM /p ADM from (b).