Constructing stochastic models for dipole fluctuations from paleomagnetic observations

Records of relative paleointensity are subject to several sources of error. Temporal averaging due to grad- ual acquisition of magnetization removes high-frequency ﬂuctuations, whereas random errors introduce ﬂuctuations at high frequency. Both sources of error limit our ability to construct stochastic models from paleomagnetic observations. We partially circumvent these difﬁculties by recognizing that the largest affects occur at high frequency. To illustrate we construct a stochastic model from two recent inversions of paleomagnetic observations for the axial dipole moment. An estimate of the noise term in the stochastic model is recovered from a high-resolution inversion (CALS10k.2), while the drift term is estimated from the low-frequency part of the power spectrum for a long, but lower-resolution inversion (PADM2M). Realizations of the resulting stochastic model yield a composite, broadband power spectrum that agrees well with the spectra from both PADM2M and CALS10k.2. A simple generalization of the stochastic model permits predictions for the mean rate of magnetic reversals. We show that the reversal rate depends on the time-averaged dipole moment, the variance of the dipole moment and a slow timescale that characterizes the adjustment of the dipole toward the time-averaged value. Predictions of the stochastic model give a mean rate of 4.2 Myr (cid:1) 1 , which is in good agreement with observations from marine magnetic anomalies. 2017 Elsevier B.V. All rights reserved.


Introduction
The spectrum of fluctuations in the geomagnetic dipole offers insights into the origin of the magnetic field and the dynamics of Earth's core (Constable and Johnson, 2005). Each distinct timescale bears the fingerprints of the underlying physical processes (e.g. Sakuraba and Hamano, 2007. Paleomagnetic observations are essential for characterizing the long-term behavior, yet no single source of information is sufficient to capture the full range of dynamics. Instead, we require an integrated approach to combine different types of measurements into a composite record that spans a broad range of timescales. One important source of information comes from measurements of relative paleointensity in marine sediments (Valet, 2003). Records are stacked and calibrated using independent estimates of absolute paleointensity to produce models for the virtual axial dipole moment (VADM) over the past two million years (Valet et al., 2005;Ziegler et al., 2011). Sediments acquire a magnetization over several thousand years (Roberts and Winkholfer, 2004), so the true signal is averaged in time. Uncertainties in dating can have a similar affect because paleomagnetic records from different times may be stacked together.
Higher resolution records have been obtained for the past 10 kyr using a combination of archeomagnetic and lake sediment data. These data have improved spatial resolution, so the geomagnetic field can be expanded in low-degree spherical harmonics (e.g. Korte and Constable, 2011). Even higher resolution records are available from historical observations (Jackson et al., 2000). Taken together these records provide a comprehensive description of fluctuations in the dipole field, but the task of combining these results into a single coherent model is a challenge.
Stochastic models are a useful tool because they enable quantitative predictions over a range of timescales. This facility is important for combining different types of data with different levels of temporal resolution. There is also good reason to think that stochastic models can represent the relevant processes in the core. Stochastic models have been constructed from geodynamo simulations with only a few model parameters, yet these models are capable of reproducing most of the variability in these simulations (Kuipers et al., 2009;Buffett et al., 2014;Bouligand et al., 2016 Synthetic studies using geodynamo simulations are an ideal test of the general approach because the simulations have relatively low numerical error and we can control the temporal resolution of the output. None of these advantages apply when we use paleomagnetic observations to construct stochastic models. Significant errors are present in the estimates of the dipole field, which affect the construction of the stochastic model. We also need to deal with temporal averaging because it limits our ability to sample the stochastic process. The goal of this study is to address the practical limitations of dealing with paleomagnetic observations and to devise a strategy for constructing models that can explain both paleomagnetic and historical records. We focus primarily on the power spectrum of dipole fluctuations, although we find that the resulting stochastic models can also account for the observed reversal rate and the duration of polarity transitions.

Stochastic description of dipole fluctuations
Stochastic models were introduced by Langevin (1908) to describe Brownian motion. A small particle in water was assumed to move under the combined influence of viscous resistance and a random force due to collision with (unseen) water molecules. The viscous force was treated as a slowly varying deterministic quantity, whereas the force due to collisions with water molecules was treated as a rapidly fluctuating random process.
Brownian motion serves as a loose analogy for the evolution of the geomagnetic dipole moment. The deterministic part of the dipole moment represents the opposing influences of dipole decay and the time-averaged dipole generation. Rapid fluctuations in dipole generation about the time average can be attributed to (unseen) turbulent flow, which we treat as a random process. We denote the axial dipole moment by xðtÞ and describe its time evolution using a stochastic differential equation (Van Kampen, 1992) dx dt where the drift term, vðxÞ, describes the deterministic part of the evolution and the noise term, DðxÞ, defines the amplitude of the random part. The time dependence of the random process, CðtÞ, is assumed to be Gaussian with a vanishing time average < CðtÞ >¼ 0: We also assume that the correlation time of the noise source is short compared with the sampling of xðtÞ. In this case the autocovariance function of CðtÞ can be approximated by a Dirac delta function, where the factor of two is a common convention (e.g. Risken, 1989).
Estimates for vðxÞ and DðxÞ can be extracted from a realization of the stochastic process (e.g. Friedrich et al., 2011). The drift term is defined by < xðt þ DtÞ À xðtÞ >¼ vðxÞDt þ OðDt 2 Þ ð 4Þ and the noise term can be approximated by where the time averages are taken for a specific value of x ¼ xðtÞ. In practice, the dipole moment is divided into a finite number of bins and a time average is evaluated for each bin. The time increment, Dt, is chosen to be long enough that CðtÞ and Cðt þ DtÞ are uncorrelated, but short enough that higher order terms in Dt are small enough to neglect.
Applying (4) and (5) to the output of a geodynamo model Meduri and Wicht, 2016) shows that the drift term, vðxÞ, is well represented by where < x > denotes the time average and c is a constant that defines the inverse timescale for slow adjustments of the dipole.
A similar representation for vðxÞ has been recovered from VADM estimates (Brendel et al., 2007;Buffett et al., 2013). Very similar values for the constant, c % 0:034 kyr À1 , were reported for the SINT-2000 model of Valet et al. (2005) and the PADM2M model of Ziegler et al. (2011). By comparison, the noise term, DðxÞ, has a weaker dependence on x. It suffices for our purposes to treat D as a constant and denote its value by D eq .
Simple representations for the drift and noise terms permit closed-form solutions for the power spectrum of fluctuations about the time average (e.g. ðtÞ ¼ xðtÞÀ < x >Þ. Defining the Fourier transform of ðtÞ by the power spectrum becomes (Buffett and Matsui, 2015) S ðf Þ ¼ D eq where the power spectrum for a white noise source (with a variance of 2) is The theoretical spectrum in (8) agrees well with a direct calculation of the power spectrum from a geodynamo model (see Fig. 1). Departures at high frequency can be improved by allowing for the influence of correlated noise (Buffett and Matsui, 2015;Bouligand et al., 2016). The spectrum for ðtÞ with correlated noise (denoted by S c ðf Þ) reduces the power at high frequencies, but it does not change the behavior at low frequencies. It is important to note that the drift and noise terms are recovered from the geodynamo model using (4) and (5) with a time difference of Dt ¼ 1 kyr. No longperiod information goes into the estimation of vðxÞ and DðxÞ, yet  , compared to predictions from two stochastic models. One stochastic model assumes a white noise source and the other assumes correlated noise. Both models are capable of predicting the low-frequency fluctuations even though the drift and diffusion terms are constructed from short-period information.
the resulting predictions are in good agreement with the lowfrequency part of the spectrum. This result suggests that a simple stochastic model offers a good description of long-period dipole fluctuations.

Recovering the drift and noise from paleomagnetic models
Several complications arise when the drift and noise terms are computed from paleomagnetic models of the dipole moment. One complication is due to random error and the other is due to temporal averaging of the fluctuations. We explore both of these complications before proposing a possible solution.

Influence of random error
Random error alters the estimates of the dipole moment, so the drift and noise terms are computed from yðtÞ ¼ xðtÞ þ gðtÞ ð 10Þ which includes a time-dependent error gðtÞ. The drift term becomes on substituting for yðtÞ from (10). The presence of random error alters vðyÞ but the time average of the error in (12) is expected to vanish. The same is not true for the noise term. Using yðtÞ to evaluate DðyÞ gives which can be rearranged into the form on introducing Dx ¼ xðt þ DtÞ À xðtÞ and Dg ¼ gðt þ DtÞ À gðtÞ. Even when Dx and Dg are uncorrelated, and gðtÞ represents the effects of white noise, we are left with (Hoze and Holeman, 2015) DðyÞ ¼ DðxÞ þ where r 2 g is the variance of the error. Thus the influence of random error becomes acute when Dt is small. On the other hand, larger Dt causes the higher order terms in (4) and (5) to become more important.
We illustrate the problem using a synthetic example. Consider a stochastic model with a linear drift term (c ¼ 0:034 kyr À1 ) and a constant noise term (D eq ¼ 0:069 Â 10 44 A 2 m 4 kyr À1 ). These numerical values were recovered by Buffett et al. (2013) from model PADM2M of Ziegler et al. (2011). A numerical realization of the stochastic model is run for 2 Myr with values of xðtÞ recorded at 1 kyr intervals. Next we add uncorrelated and normally distributed random error to produce a noisy record, yðtÞ, where the standard deviation of the error is r g ¼ 0:5 Â 10 22 A m 2 . Finally, we recover a constant value for both DðxÞ and DðyÞ from xðtÞ and yðtÞ; we denote these constants by D x eq and D y eq to explicitly identify the input time series. Fig. 2 shows the estimates for D x eq and D y eq as a function of Dt. At the shortest time difference, Dt ¼ 1 kyr, we obtain D x eq ¼ 0:068 AE 0:002 Â 10 44 A 2 m 4 kyr À1 and D y eq ¼ 0:297 AE 0:010Â 10 44 A 2 m 4 kyr À1 , where the uncertainties represent one standard deviation. These results are consistent with expectations from (15). Large deviations from the true value of D eq are found with the noisy record when Dt is small. Smaller deviations occur as Dt increases, although these errors remain relatively large. On the other hand, the value recovered from the error-free record, xðtÞ, is reliable at small Dt but slowly departs from the known value as Dt becomes larger. Consequently, we cannot deal with the influence of random error by arbitrarily increasing Dt.

Influence of temporal averaging
Temporal averaging of the paleomagnetic record can arise in several ways. Errors in dating allows measurements at different times to be stacked. In addition, magnetization is acquired in sediments over several thousand years (Roberts and Winkholfer, 2004). A prolonged acquisition time removes high-frequency variations and affects our ability to sample the stochastic process at short Dt. One way to deal with the problem of averaging is to treat the measured record as a filtered version of the true signal (e.g. Leonard, 1974). We define the measured signal, xðtÞ, as where the filter function, gðtÞ, smooths the true signal, xðtÞ, over some prescribed averaging time (denoted by T). Two popular filter functions are the box-car and gaussian filters (see Fig. 3). The true signal is convolved with a suitable filter function to produce the measured record. The paleomagnetic record, xðtÞ, still obeys a stochastic differential equation, but it is not the same as the differential equation in (1). Applying the filter to (1) gives where we have adopted a constant noise term and a linear drift term. The only difference in (17) is that the random process is driven by CðtÞ rather than CðtÞ. A power spectrum for ¼ xÀ < x > is defined by taking the Fourier transform of (17). Solving for ðf Þ gives Fig. 2. Estimates for the noise term, DeqðDtÞ, computed from exact xðtÞ and noisy yðtÞ time series. Addition of random error to yðtÞ causes D y eq to depart from the known value Deq ¼ 0:069 Â 10 44 A 2 m 4 kyr À1 . Calculations using xðtÞ reproduce the known value to within the uncertainties at Dt ¼ 1 kyr. Discrepancies in D x eq ðDtÞ increase slowly with Dt due to unmodelled contributions from higher-order powers in Dt.
and the power spectrum becomes (Rice, 1954) S ðf Þ ¼ D eq and gðf Þ is the Fourier transform of the filter function. Eq. (20) follows from the convolution theorem (e.g. Bracewell, 1999) because convolution in the time domain corresponds to multiplication in the frequency domain Power spectra for xðtÞ and xðtÞ are the same at low frequencies because gðf Þ ! 1 as f ! 0 (see Fig. 3).
We illustrate the consequences of time averaging using the stochastic model from Section 3.1. A 2-Myr realization is sampled at 1-kyr intervals and a smoothed version is produced using a boxcar filter with an averaging time of T ¼ 3 kyr. Fig. 4 shows the power spectrum of the filtered signal, xðtÞ, compared with the theoretical spectrum from (19). We also show the power spectrum for the original (unfiltered) time series, xðtÞ, versus the theoretical spectrum from (8). Both theoretical spectra are in good agreement with the direct calculations from xðtÞ and xðtÞ. Undulations in the spectrum of xðtÞ is a consequence of the box-car filter, which is oscillatory in the Fourier domain. The main conclusion from this example is that temporal averaging affects only the highfrequency behavior of the record. The filtered dipole moment still obeys a stochastic differential equation and the spectrum is still reliably predicted at low frequencies from the drift and noise terms. Conversely, the low-frequency part of the spectrum constrains the drift and noise terms of the stochastic model. Fig. 5 shows the noise term, D eq , recovered from xðtÞ and xðtÞ as a function of Dt. The most reliable estimate for D eq comes from xðtÞ at the shortest possible Dt (1 kyr in this case). Temporal filtering substantially reduces the estimate of D eq at low Dt, although the recovered value approaches a constant once Dt exceeds the filter width T. A rule of thumb based on the spectrum of the filter (say gðf Þ > 0:9) is that Dt should be roughly twice T. Sampling the process at Dt ¼ 6 kyr gives an estimate for D eq that is nearly independent of Dt. Unfortunately, this estimate is well below the known value (e.g. 0.044 versus 0.069). A similar departure in D eq at Dt ¼ 6 kyr is inferred from xðtÞ (e.g. 0.063 versus 0.069), although the error from the unfiltered time series is much smaller.
The preceding results show that temporal averaging can affect the amplitude of the noise term, particularly when Dt is smaller than the duration of the averaging. Estimates for D eq appear to approach a constant value once Dt > 2T, although this constant can be significantly less than the true value. On the other hand, random noise causes the recovered estimate of D eq to exceed the known value by an amount r 2 g =Dt, where r 2 g is the variance of the error. Both temporal averaging and random noise have the largest affect on the high-frequency part of the spectrum. Averaging removes power at high frequency, whereas random error introduces power across all frequencies, although it is most evident at high frequency. Consequently, the low-frequency part of the spectrum is relatively unaffected by both sources of error. We exploit  this result to construct a broadband paleomagnetic power spectrum.

A composite paleomagnetic power spectrum
We use two sources of information to construct the paleomagnetic spectrum. Model PADM2M of Ziegler et al. (2011) gives the axial dipole moment over the past 2 Myr at intervals of 1 kyr, whereas CALS10k.2 (Constable et al., 2016) gives the axial dipole moment (and other low-degree components of the magnetic field) over the past 10 kyr at 50-year intervals. Fig. 6 shows the power spectrum for each model, calculated using a multi-taper method (function pmtm in Matlab). We also show two theoretical spectra. One spectrum is predicted using the parameters of a simple stochastic model derived from PADM2M (Buffett et al., 2013), which gave c ¼ 0:034 kyr À1 and D eq ¼ 0:069 Â 10 44 A 2 m 4 kyr À1 .
The second spectrum is obtained by applying a gaussian filter to the stochastic model, using an averaging time of T ¼ 2:4 kyr. The sampling used to construct the stochastic model from PADM2M was Dt ¼ 5 kyr, so the filter required to account for the power spectrum of PADM2M is broadly compatible with the proposed rule of thumb Dt % 2T.
CALS10k.2 possesses more power than PADM2M at overlapping frequencies. One interpretation is that temporal averaging has a greater influence on PADM2M, which acts to reduce the power at high frequencies. We might remedy this problem by seeking an independent estimate for D eq using the higher resolution record from CALS10k.2. Fig. 7 shows the resulting estimate for D eq as a function of Dt. The noise term initially increases with Dt, implying temporal averaging and/or correlated noise in the stochastic model. A simple parametric fit of the form D eq ðDtÞ ¼ D eq ð1Þð1 À e ÀDt=T Þ ð 23Þ gives D eq ð1Þ ¼ 0:34 Â 10 44 A 2 m 4 kyr À1 for the asymptotic value of the noise term. We fit (23) through the lower limit of the estimates in Fig. 7 to account for the possible influence of random error (which tends to increase D). A correlation time of T ¼ 120 years suggests that the sampling of the stochastic process should be restricted to Dt > 240 years. (We adopt Dt ¼ 300 years as a lower limit in our subsequent discussion.) Even though the noise term from CALS10k.2 is more than four times larger than that from PADM2M, there are several reasons to think that the estimate from CALS10k.2 is more reliable. First, the random error in CALS10k.2 is lower than that in PADM2M (Ziegler et al., 2011). Second, the short correlation time for CALS10k.2 implies much less temporal averaging. Both of these factors should reduce the errors described in the previous section. An independent assessment of the combined influence of random error and temporal averaging for CALS10k.2 suggests that the net contribution to D eq is less than 25% (see Appendix). On the other hand, CALS10k.2 is too short to reliably recover the drift term, vðxÞ, because we require sampling over a range of x. Instead, we must rely on the longer record from PADM2M to estimate the slope of the drift term.   Recovering a direct estimate of vðxÞ from PADM2M using (4) is not an optimal approach. Increasing the value of D eq without any change in c predicts more power at low frequencies, which is incompatible with the low-frequency part of the PADM2M spectrum. Since we expect random error and temporal averaging to have less affect at low frequencies, we choose to alter c to maintain agreement with PADM2M at low frequencies. In effect, we are using the low-frequency spectrum of PADM2M to estimate c once D eq is inferred from CALS10k.2. The predicted power at low frequency is 2D eq =c 2 , so we take c ¼ 0:075 kyr À1 with D eq ¼ 0:34 Â 10 44 A 2 m 4 kyr À1 to retain consistency with the lowfrequency power in PADM2M. While the slope of the drift term is more than twice the value recovered from PAD2M using (4), it is in rough agreement with the value c ¼ 0:07 kyr À1 estimated for the PISO-1500 model (Channell et al., 2009). (The noise term for PISO-1500 near x ¼< x > is 0:54 Â 10 44 A 2 m 4 kyr À1 , which is somewhat higher than the value recovered from CALS10k.2. Interestingly, the preferred sampling interval for PISO-1500 is 4-5 kyr, which is close to the sampling interval adopted previously for PADM2M and SINT-2000 (Buffett et al., 2013). Thus the time averaging in all three VADM models is roughly the same and may reflect the gradual acquisition of magnetization in sediments). Selective use PADM2M and CALS10k.2 may seem ad hoc, but the underlying motivation is to leverage the relative strengths of these two models. Alternative strategies are possible, such as the maximum likelihood approach of Kleinhans (2012), although we are not aware of applications where these methods have been applied to multiple data sets with different resolutions, durations and accuracies. We now test the stochastic model by assessing the internal consistency and by making predictions for the reversal rate, which can be compared against geological observations. Internal consistency is tested by running a series of 100 realizations of the stochastic model. The model parameters are held fixed but the initial conditions and random component differ between realizations. These realizations are run for 2 Myr and the value of the dipole moment is recorded every 300 years, corresponding to the sampling interval inferred from CALS10k.2. This choice ensures that the noise source, CðtÞ, can be approximated as uncorrelated. The origin of correlation at shorter time steps might reflect the lifetime of convective eddies in the core or possibly the time required to sweep normal and reversed patches of magnetic flux across the surface of the core (Metman et al., 2017). While we could account for correlated noise in the stochastic model (Buffett and Matsui, 2015;Bouligand et al., 2016), we avoid this complication in the realizations by choosing a sufficiently large time step.
A power spectrum is computed for each realization and the results are superimposed on the power spectra computed from PADM2M and CALS10k.2 (see Fig. 8a). We also show a theoretical composite spectrum, based on the parameters of the stochastic model and now including the influences of correlated noise with a correlation time of T ¼ 120 years (e.g. Buffett and Matsui, 2015). A cloud of power spectra for the realizations overlap the low-frequency part of the PADM2M power spectrum and the theoretical spectrum. Comparison of the realizations in Fig. 8a with the power spectrum for CALS10k.2 is not appropriate because the CALS10k.2 spectrum is computed from a 10-kyr time series. A better comparison would rely on 10-kyr realizations (see Fig8b). A series of shorter realizations produces a cloud of power spectra that overlap the computed power spectrum for CALS10k.2, suggesting that the revised stochastic model is broadly consistent with the CALS10k.2 model.
We also test the stochastic model against historical observations (Jackson et al., 2000). A steady decrease in the dipole field has lowered the dipole moment by Dx ¼ 0:68 Â 10 22 A m 2 over a 150-year interval between 1860 and 2010 (Gillet et al., 2013). Such a change is too large to be caused by the drift term, so it must be associated with the noise term. The root-mean-square (rms) variation in the dipole moment due to the noise term is Using the revised estimate of D eq and Dt ¼ 0:15 kyr, we find < Dx 2 > 1=2 ¼ 0:32 Â 10 22 A m 2 . Thus the historical variation is larger than the expected variation, but it is not implausible. A realization of the noise process is described by Risken (1989) Dx where w is a random variable drawn from a standard normal distribution (mean of zero and variance of 1). We require w ¼ 2:13 to account for the recent variation in the dipole field, which would occur about 1.7% of the time. The actual probability could be somewhat lower if the noise source is correlated at Dt ¼ 0:15 kyr (a likely case given our estimate of the correlation time from CALS10k.2). The preceding estimate would then represent a modest overestimate of the probability of occurrence. By comparison, the original value of D eq ¼ 0:069 Â 10 44 A 2 m 4 kyr À1 from PADM2M would require w ¼ 4:73 to account for the historical variation. Such an event would occur less than 0.0001% of the time, so the historical record lends support to the larger value for the noise term. Another useful prediction of the stochastic model is the variance of the dipole moment. We obtain an expression for the variance, r 2 x , by integrating the power spectrum over frequency The revised values for D eq ¼ 0:34 Â 10 44 A 2 m 4 kyr À1 and c ¼ 0:075 kyr À1 give r x ¼ 2:13 Â 10 22 A m 2 . While this value exceeds the estimate r x ¼ 1:48 Â 10 22 A m 2 for PADM2M (Ziegler et al., 2011), it is not too far from the estimate r x ¼ 1:97 Â 10 22 A m 2 for SINT-2000 (Valet et al., 2005) and somewhat smaller than the estimate r x ¼ 2:68 Â 10 22 A m 2 for PISO-1500 (Channell et al., 2009). Thus the predicted variance lies within the range of estimates from recent VADM models.

Geomagnetic polarity reversals
A more general representation of the drift term is needed to describe geomagnetic polarity reversals. The linear approximation in (6) is useful when x varies about < x >, but its utility ceases when x approaches zero during a reversal. The invariance of the magnetic induction equation to a change in the sign of the magnetic field suggests that vðxÞ is an odd function of x. We expect the drift term to adjust x toward the negative value of the time average once x changes sign. One simple extension of the linear approximation is where the expected symmetry is obtained by taking vðÀxÞ ¼ ÀvðxÞ.
The gradient of vðxÞ at x ¼< x > is consistent with the linear approximation in (6), but the value of the drift now vanishes at x ¼ 0. It is convenient to represent the drift as the negative gradient of a potential UðxÞ. Integrating (26) for the UðxÞ gives where the integration constant is chosen to make Uð0Þ ¼ 0. A comparison of UðxÞ with the potential recovered from the PADM2M model of Ziegler et al. (2011) is shown in Fig. 9. The barrier at x ¼ 0 is comparable for both potentials, but the amplitudes of UðxÞ differ at large jxj This is mainly a consequence of increasing c in the revised stochastic model. A larger c produces a narrower potential well and limits the variability of x at a fixed level of noise, consistent with the predicted standard deviation r x ¼ ffiffiffiffiffiffiffiffiffiffiffiffi D eq =c p . We now use the generalization of the drift in (26) to predict the rate of magnetic reversals and the duration of polarity transitions.

Rates of reversals
Random fluctuations in x enable the dipole to jump from one potential well to the other, leading to a magnetic reversal. The average frequency of this transition can be predicted using the stochastic model. Kramers (1940) derived an approximation expression for the reversal rate, r, when the barrier DU ¼ Uð0Þ À Uð< x >Þ is large compared with the noise D. Kramers' formula in our notation gives Substituting for from (27) and using the definition of the variance from (25) gives Remarkably, the rate of reversal depends on the time average, < x >, the variance, r 2 x , and the timescale for slow adjustments of the dipole field, c À1 ; the slow timescale is thought to reflect the decay time of dipole fluctuations (e.g. Gubbins and Roberts, 1987). Geodynamo simulations suggest that the dipole fluctuations can be represented by the first few decay modes . Using < x >¼ 5:3 Â 10 22 A m 2 , r x ¼ 2:13 Â 10 22 A m 2 and c ¼ 0:075 kyr À1 (75 Myr À1 ) gives r ¼ 4:2 reversals per Myr, which is comparable to the observed rate over the past 30 Myr (Lowrie and Kent, 2004). By comparison, a 60-Myr realization of the stochastic process yields 3.9 reversals per Myr when the realization is filtered to a resolution of 30 kyr, comparable to the resolution of marine magnetic anomalies (Gee and Kent, 2015). The need to filter the realization is connected to the complexity of polarity transitions when the noise term is large. We explore this question in the next section.

Duration of polarity transitions
The duration of polarity transitions depends on how the transitions are defined. A definition based on magnetic intensity might depend on the time required for the dipole to recover to the long-term average after a change in sign (i.e. a recovery time). This particular definition is useful for our purposes because it can be computed from the stochastic model. We expect the drift term to be small near x % 0, so the evolution of the dipole during the transition is dominated by the noise term. A useful approximation for the time required for the field to rise above a particular threshold, x t , is (Buffett, 2015) s ¼ where Dð0Þ refers to the value of the noise term at x ¼ 0. The general form of (31) is characteristic of a diffusive process, which includes no contribution from the drift term. A more exact treat- ment of the problem accounts for the drift term as x rises toward the threshold x t . Fig. 10 shows a comparison of the approximation in (31) with the value computed from a numerical solution of the Fokker-Planck equation (e.g. Risken, 1989). Including the drift term shortens the recovery, but the difference is relatively small when we adopt the revised value for D eq . This implies that the recovery of the magnetic field following a reversal is driven mainly by noise (e.g. random turbulent fluctuations in the field generation). We can compute a recovery time from the PADM2M model by interpolating the time when x rises above the time average after a reversal. Each reversal gives a different value for s, but the average and its standard deviation are shown in Fig. 10. The agreement with theory is surprisingly good. We also show the time required for the field to drop from the time-averaged value into a reversal (i.e. a decline time). The mean decline time from PADM2M is 41 kyr, whereas the mean recovery time is 27 kyr. This asymmetry is consistent with previous observations (Valet and Meynadier, 1993). (The decline time was incorrectly reported as the recovery time in Buffett (2015), although the main point in that study was that these short durations require a noise term in excess of D ¼ 0:30 Â 10 44 A 2 m 4 kyr À1 ).
The difference between the recovery and decline times can be attributed to the role of the drift term. The recovery time is shorter than the approximation in (31) because the drift term drives the dipole moment toward the time average, increasing the rate of adjustment after a reversal. Conversely, the dipole must work against the drift term during the decline phase. The approximation in (31) lies roughly midway between the estimates from PADM2M, which suggests that the drift lengthens and shortens the adjustment by comparable amounts, relative to a purely diffusive process with no drift term.
It is reasonable to question whether the PADM2M model can adequately resolve the recovery time when the short-period behavior is not sufficient to compute D eq . A transition that lasts s % 30 kyr would correspond to a frequency of f ¼ 1=2s, assuming the transition represents half a cycle. A nominal frequency of 0.017 cycles kyr À1 lies in the part of the spectrum where PADM2M and the stochastic model are broadly consistent (see Fig. 8). Consequently, there is internal consistency in our argument that the stochastic model is in agreement with both the transition duration and low-frequency power spectrum from PADM2M. It is encouraging that the same stochastic model gives a reasonable estimate for the reversal rate, particularly when no information about the reversal rate is used in the construction of the stochastic model. The dominance of the noise term during a polarity transition has interesting consequences for the complexity of reversals. A process that is driven solely by the noise term is analogous to a random walk. The probability of stepping back and forth across x ¼ 0 increases with the number of steps n during the transition. Dasgupta and Rubin (1998) show that the expected number of zero-crossings is proportional to ffiffiffi n p . As we decrease the step size in a numerical realization, we take a large number of steps through the transition and produce a large number of zero-crossings during a single transition. In practice the time step is limited by the correlation time of the noise source to ensure that realizations of the noise are effectively uncorrelated. Consequently, the number of steps through a transition is not arbitrarily, and a representative number is liable to produce several zero crossings.
Numerical realizations with Dð0Þ ¼ 0:30 Â 10 44 A 2 m 4 kyr À1 (close to the value D ¼ 0:34 Â 10 44 A 2 m 4 kyr À1 proposed here) produced multiple zero-crossings in about 50% of the polarity transitions when Dt ¼ 1 kyr (Buffett, 2015). The average number of zero crossings is 2.8, but this number would go up if Dt ¼ 0:3 kyr. We could expect 3Â more time steps through a transition and roughly ffiffiffi 3 p Â more zero crossings (on average), corresponding to a total of 5 changes in sign during a transition. To make meaningful comparisons with geological observations we would want to remove these short-period polarity changes by filtering the numerical realization to the resolution of the observations. In the previous section we used T ¼ 30 kyr to compare the reversal rate with estimates from marine magnetic anomalies.

Conclusions
Stochastic models have been successfully tested using geodynamo simulations, but their use with paleomagnetic observations requires departures from the standard approach. Two main difficulties are identified. The first is due to random error in the estimates of the dipole moment, which cause the noise term to be over-estimated. The significance of this problem depends on the sampling interval, Dt, and the largest affects occur at short Dt. A second difficulty arises from temporal averaging of dipole fluctuations, either due to errors in dating or gradual acquisition of magnetization in sediment. In either case, temporal averaging reduces the noise term at short Dt, although estimates for D often converge to a constant value as Dt increases. Unfortunately, the noise term does not necessarily converge to the correct value.
An important feature of both random error and temporal averaging is that the largest affects are predicted at high frequency. Because the low-frequency behavior is nearly unaltered, we can use the low-frequency part of the observed power spectrum as a constraint on the stochastic model. We illustrate the approach using the PADM2M model of Ziegler et al. (2011) and the CALS10k.2 model of Constable et al. (2016). An estimate of the noise term is recovered from the high-resolution CALS10k.2 model, while the slope of the drift term, c, is estimated from the lowfrequency part of the spectrum for PADM2M. Realizations of the stochastic model yield a composite power spectrum that agrees reasonably well with both PADM2M and CALS10k.2.
A simple generalization of the stochastic model is needed to allow large deviations from the time-averaged moment. This  (31), where the drift term is assumed to vanish. Discrete estimates from PADM2M are shown for the recovery and decline times. The recovery time agrees well with the theoretical estimate for recovery, whereas the decline time exceeds the recovery time due to contributions from the drift term. modification enables predictions for the mean rate of reversal. A reversal in the stochastic model occurs when a realization jumps between the minima in a double-well potential. Application of Kramers' formula (Kramers, 1940) gives a surprisingly simple expression for the reversal rate. We find that the reversal rate can be defined in terms of the time-averaged dipole moment, the variance of the dipole moment and a slow timescale that characterizes the adjustment of the dipole toward the time-averaged value. Using values from the stochastic model gives a mean rate of 4.2 Myr À1 , which is good agreement with observations (Lowrie and Kent, 2004). Comparable rates are obtained from realizations of the stochastic process, provided we filter the realization to the same resolution as the observations. The need for temporal filtering arises from the importance of noise in driving polarity transitions. Multiple polarity changes can occur within a single transition field, so a quantitative comparison with observations depends on the temporal resolution of those observations.

Acknowledgments
This work is supported by the National Science Foundation -USA (EAR-1644644) and by a Summer Undergraduate Research Fellowship (SURF) from Caltech.

Appendix A. Combined influence of random error and time averaging
Random error and time averaging have opposite effects on estimates for the noise term in a stochastic model. Random error increases the noise term, whereas time averaging tends to reduce D eq when Dt is less than roughly twice the averaging time. The combined influence of both error sources depends on their relative magnitude. As an example, we consider a realization that approximates the CALS10k.2 model. We adopt c ¼ 0:075 kyr À1 and D eq ¼ 0:34 Â 10 44 A 2 m 4 kyr À1 for the parameters of the stochastic model and run a 10-kyr realization with time steps at a 50-year interval. Random error with a standard deviation of r g ¼ 0:4 Â 10 22 A m 2 is added to the realization and a box-car filter is applied with an averaging time of T ¼ 150 years. All of these values are chosen to approximate the statistical properties of CALS10k.2. Fig. A.1 shows the recovered estimates for D eq as a function of the sampling time Dt for both the time-averaged process and the time-averaged process with random error. Estimates for D eq from the time-averaged process increase with Dt and approach a constant value at large Dt. We would obtain D eq ¼ 0:32Â 10 44 A 2 m 4 kyr À1 for the recommended sampling interval of Dt ¼ 2T ¼ 300 years. Addition of random error to the process produces much higher estimates for D eq at small Dt, consistent with expectations for random errors. We obtain D eq ¼ 0:42 Â 10 44 A 2 m 4 kyr À1 at Dt ¼ 300 years and somewhat lower values if we fit a curve through the lower limit of estimates (as done in the text). The discrepancy from the known value is about 24% or less. The time-averaged process with random error looks qualitatively different than the estimates recovered from CALS10k.2 (see Fig. 7). In particular, the estimates from CALS10k.2 are much lower at small Dt. One interpretation is that the actual random errors in CALS10k.2 are lower than r g ¼ 0:4 Â 10 22 A m 2 . Alternatively, we might appeal to a longer effective time-averaging (i.e. larger T), which would lower the amplitude of the time-averaged error. On the other hand, it would be difficult to account for the observed dependence of D eq on Dt in CALS10k.2 when T was much larger than 150 years. It is possible that the reported error for CALS10k.2 includes offsets or biases, which would not affect the recovered estimates for D eq . Assessing the error in paleomagnetic models is a difficult challenge, and it is likely that the error in CALS10k.2 is more complicated than the simple representation considered here. Still, the preceding example suggests reasonable estimates for D eq can be recovered from a shorter, high-resolution model. Recovery of noise term Deq from a 10-kyr realization of the preferred stochastic model. A box-car filter is applied to the realization with and without the addition of random error. Sampling the process at Dt ¼ 2T ¼ 300 years gives a reasonable estimate for Deq, even when a realistic level of random error is added to the process.