Rapid Changes in Strength and Direction of Earth's Magnetic Field Over the Past 100,000 Years

Previous studies of rapid geomagnetic changes have highlighted the most extreme changes in direction and field strength found in paleomagnetic field models over the past 100 ky. Here we study distributions of rates of change in both time and space. Field models based on direct observations provide the most accurate values for rates of change, but their short duration precludes a complete description of field behavior. Broader representation is provided by time‐varying paleofield models, here including GGF100k, GGFSS70, LSMOD.2, CALS10k.2, HFM.OL1.A1, pfm9k.2, and SHAWQ‐iron age although variability across models and lack of temporal and spatial resolution of fine scale variations make direct comparisons difficult. For the paleofield we define rapid changes as exceeding the peak overall value of 0.4° yr−1 for directional changes and 150 nT yr−1 for intensities as established by the gufm1 model spanning 1590–1990 CE. We find that rapid directional changes are associated with low field strength and can spread across all latitudes during such episodes. Distributions of directional rates of change exhibit high skewness for models that include excursions. Rates of change in field intensity exceeding 150 nT yr−1 arise in brief intervals during the Holocene particularly associated with the strong field Levantine Iron Age Anomaly. Around the Laschamp excursion there are also rare localized occurrences of rapid intensity change. Limitations in current models make it difficult to define absolute rates for past changes, but we see that rapid changes are essential field characteristics not observed in the modern field that should nevertheless be regarded as an essential for Earth‐like dynamo simulations.


Introduction
Extreme changes in Earth's magnetic field are found in the paleomagnetic record and the largest of these are associated with geomagnetic excursions and reversals.However, we still do not know how rapidly these can occur or the details of the core dynamics that give rise to them.Efforts to improve understanding have followed two parallel tracks: global field models and dynamo simulations.Recently developed time-varying spherical harmonic field models based on paleomagnetic data now span several excursions over the last 100 ky (Brown et al., 2018;Korte et al., 2019;Panovska et al., 2021), with the prospect for improving these with better data and extending further back in time as recently done for the Matuyama-Brunhes reversal (Mahgoub et al., 2023).In parallel, numerical dynamo simulations have pushed to increasingly extreme physical conditions (Aubert et al., 2017;Schaeffer et al., 2017) and have been compared in detail to magnetic field observations spanning a wide range of timescales (Christensen et al., 2010;Davies & Constable, 2014;Sprain et al., 2019, for example), including exploring the fundamental core dynamics associated with rapid geomagnetic events (Davies & Constable, 2018, 2020).Each of these approaches has important limitations that need to be addressed to facilitate understanding of the full spectrum of geomagnetic field changes.The work presented here focuses on garnering information from paleofield modeling and determining how it might be related to actual rates of geomagnetic change.
The power spectrum of variations in the axial dipole moment as a function of frequency based on a selection of paleofield and observational models (Constable & Constable, 2023) is shown in Figure 1.The increasing general falloff rate with frequency clearly demonstrates that the largest changes in the axial dipole occur at the longest periods.Black arrows on this figure correspond to indicated transitions in the approximate exponent of the power law in frequency.In contrast, the spectrum of the time derivative in Figure 2 reveals that the largest rates of change in dipole moment occur at periods between 100 years and 50 ky, and that there may be more structure than the simple f 2 falloff across this range that is depicted in Figure 1.At much higher frequencies the rates of change in axial dipole moment contribute only a fraction to the overall spatial power spectrum of secular variations for the global field model gufm1 covering the past 400 years as discussed further in Section 2.3.
For the modern magnetic field, direct observations including those at geomagnetic observatories and low earth orbiting satellites provide excellent data used to construct global time varying field models of low-degree poloidal field B, from which it is possible to calculate rates of change and even accelerations, and to study wave motions and/or geomagnetic jerks.Over the past 400 years the rates of intensity change, ∂B⁄∂t, display modest peak values of ∼150 nT yr 1 , while directional changes ∂ B/∂t yield maximum values less than 0.4°yr 1 (Finlay et al., 2020;Jackson et al., 2000).By contrast, intensity changes exceeding 750 nT yr 1 have been reported in the Levantine region around 1000 years BCE (Ben-Yosef et al., 2017;Korte & Constable, 2018;Shaar et al., 2016), while localized directional changes of 1°yr 1 or more have been recovered in a number of mid-latitude regions, often (but not always) coeval with the most recent excursions and reversals as summarized in Figure 1 of Maffei et al. (2021).The spatial scale over which such rapid changes can be supported is unclear, but it must be fairly broad in order for the signal to originate from the core.For example, Davies and Constable (2017) argued that the Levantine Iron Age Anomaly (LIAA) must have spanned at least 60°at Earth's surface.Davies andConstable (2018, 2020) identified maximum rates of change in geomagnetic field strength and intensity and highlighted the fact that there are broad similarities between results from both paleofield models and numerical simulations.They found that the largest directional changes typically occur at low to moderate latitudes in both.In numerical simulations extremal changes in intensity are more often found at high latitude, while the information about such changes in paleofield models remains somewhat uncertain.From the numerical simulations, Davies and Constable (2018) argued that peak intensity changes occur due to rapid intensification of normal polarity patches migrating across the core surface, while Davies and Constable (2020) found that peak directional changes were caused by movement of regionally reversed flux across the core surface.
Here we go beyond the initial analyses of extremal events (which uncovered only peak rates of change over time and space within each studied model) and provide a more comprehensive analysis of distributions of rates of change in both time and space, focusing on 8 global time-dependent field models that span different parts of the past 100 ky.In particular, we focus on "rapid changes" which we define as faster than what is found in the historical geomagnetic field over the last 400 years.This provides a perspective on the broad issue of the degree to which temporal field variations over the historical era reflect variations over longer timescales.We address the following questions: (a) How do rapid changes in direction and intensity vary with time and latitude?(b) Are there spatio-temporal relations between rapid changes in direction and intensity?(c) Do rapid changes reflect specific features of the core field, for example, regions of anomalously strong/weak field?

Time-Varying Paleofield Models and Rates of Change
Our analyses are based on a collection of time-varying spherical harmonic geomagnetic field models which we loosely divide according to three time intervals: (a) the past 400 years where direct observations are available; (b) the Holocene/archeomagnetic time interval, and (c) spanning all or some substantial fraction of the past 100 ky.In general the spatio-temporal resolution of the underlying data decreases as longer timescales are considered; moreover, within each time interval there are distinct limitations on our ability to detect rapid field changes leading to some difficulties in making an explicit numerical definition of the rate of change that defines an unusually rapid geomagnetic event.
Details about the 8 global time-varying field models have been reported elsewhere so we only provide a brief description here.All except 2 are parameterized spatially in global spherical harmonics and used cubic splines to represent the time variations.The regularized inversion technique used in construction is intended to minimize unnecessary structure in the resulting model.

Direct Observations
As a representation of the modern field we use the canonical gufm1 historical field model (Jackson et al., 2000), which covers the period 1590-1990 CE, and is derived from directional measurements recorded in ships logs, land surveys, and geomagnetic observatories.Full vector data are only available post 1840 CE and the g 0 1 axial dipole coefficient prior to this time is based on a linear extrapolation back in time.In Figures 1 and 2 only the post-1840 part of the record was used in the spectrum calculations of Constable and Constable (2023).Our subsequent calculations use the whole time range of the model from 1590 to 1990 CE, and we note here that the rate of change for g 0 1 prior to 1840 CE is not well constrained by direct historical observations (Finlay, 2008).Average values are around 2.3 ± 2.7 nT yr 1 in reasonable agreement with values for Holocene field models.gufm1 extends to spherical harmonic degree and order 14, and the cubic splines have temporal knots every 2.5 years.

The Past 10 ky
Two of our models span the past 10 ky, CALS10k.2 and HFM.OL1.A1 (Constable et al., 2016), pfm9k.2(Nilsson et al., 2022) extends to 9 ka, while SHAWQ-iron age (Osete et al., 2020) which we denote SHAWQ-IA, covers a much shorter time frame from 1300 BCE to 0 CE.Because they include a number of more recent data SHAWQ-IA and pfm9k.2 are each expected to have higher resolution results than CALS10k.2 and HFM.OL1.A1 during the Levantine Iron Age Anomaly, a time of particularly rapid growth and decay in regional field strength.The CALS10k.2 and HFM.OL1.A1 and SHAWQ-IA models extend to spherical harmonic degree and order 10, while pfm9k.2 is limited to degree and order 5. CALS10k.2 and HFM.OL1.A1 are both built from globally distributed archeomagnetic and sediment records and use cubic splines with knots every 40 years.SHAWQ-IA uses only archeomagnetic data and is thus strongly influenced by mid-latitude northern hemisphere records.And, although it was built using a cubic spline parameterization in time (Campuzano et al., 2019), the available public record (https://earthref.org/ERDA/2453/) used here is provided in the form of snapshot models every 25 years.The strategy for constructing pfm9k.2 is quite different as it adopts a probabilistic Bayesian approach to jointly model ages and field (Nilsson & Suttie, 2021) and uses archeomagnetic data together with a limited selection of sediment records.We have used the posterior mean model time series for the Gauss coefficients every 50 years available at https://earthref.org/ERDA/2522.

The Past 100 ky
On longer time scales GGF100k (Panovska et al., 2018), GGFSS70 (Panovska et al., 2021), and LSMOD.2 (Korte et al., 2019) provide time-varying models spanning 0-100 ka, 15-70 ka, and 30-50 ka, with knot spacings of 200, 50, and 50 years, respectively.Both GGF100k and LSMOD.2 extend to spherical harmonic degree and order 10, while GGFSS70 stops at 6.The models have a range of limitations and inconsistencies related to data selection, quality, and coverage in space and time.We take gufm1 as our standard for comparison and representative of a snapshot of 400 years of variation, and note that in the paleofield models lack of data coverage leads to non-uniformity in both temporal and spatial resolution.In general, the best coverage is at mid-latitudes in the northern hemisphere.There is a particular paucity of southern hemisphere records for SHAWQ-IA because of the restriction to archeomagnetic contributions, and much of the southern hemisphere structure may reflect more about the regularization imposed in modeling than actual field structure.The same is true for regions with minimal data in other paleofield models.
Figure 3 provides a partial guide to the comparative temporal and spatial resolution across the 8 models.Solid lines represent the average spatial power in the field, R l , defined (Lowes, 1974) in terms of the Schmidt normalized Gauss coefficients g m l , h m l for spherical harmonic degree l and order m The dashed lines are average power in the "instantaneous secular variation," S l , 2 with ġm l , ḣm l representing the time derivative of the corresponding Gauss coefficients.We must recognize that the term "instantaneous" is used loosely here and depends on the temporal resolution inherent to each model, which ranges from about 2.5 years for gufm1 to a few hundred years for GGF100k and imparts an implicit temporal smoothing on our results.While gufm1 is often taken as the standard to represent the geomagnetic field, and certainly has higher spatial resolution, it is important to note that the instantaneous secular variation does not provide an adequate sample of field variability.This is especially clear in (but not limited to) the l = 1, dipole term.In Figure 3, we also see that despite the dipole dominance in the spatial power spectra at Earth's surface (solid lines), the black dashed line shows that the current contribution to spatial power in the secular variation is Figure 3. Average spatial power spectra, R l , of magnetic field (solid lines) and average spatial power in "instantaneous secular variation," S l (dashed lines) at Earth's surface for gufm1 and seven paleofield models considered in this paper.largest at spherical harmonic degree 2, and still significant for the higher degree terms which can be difficult to resolve in much longer term field models.
For each field model we calculate the local magnetic field vector B and its rates of change on a 5°× 5°spatial grid at 10 years intervals.The cubic spline parameterization allows for an analytical calculation of temporal rates of change in the Gauss coefficients in all cases except for SHAWQ-IA and pfm9k.2where we assumed a constant rate of change based on linear interpolation between Gauss coefficients provided for each snapshot.All the rates of change were propagated into rates of change in intensity, local field directions, virtual dipole moments (VDM), and angular virtual geomagnetic pole (VGP) positions.This analytical calculation differs slightly from the strategy used in our previous papers (Davies & Constable, 2020) where we used first differences in time for each element.We find that the maximum difference between the two methods is a few percent and is generally less, producing no impact on the overall interpretations.Changes in field intensity and direction are denoted ∂B⁄∂t and ∂ B/∂t respectively.Changes in the equivalent Virtual Dipole vector P V are denoted ∂P V ⁄∂t and ∂ PV /∂t respectively.

Global Rates of Change in the Geomagnetic Field
Table 1 presents summary statistics of the globally sampled rates of change for field strength for the 8 global timedependent field models considered.Both local field intensity (nT yr 1 ) and VDM (ZAm 2 yr 1 ) are provided and we distinguish mean absolute rates of change from signed values.Directional rates of change are given in °yr 1 for both local field vectors and VGP positions and are intrinsically positive.
For both ∂B⁄∂t and ∂P V ⁄∂t, the overall mean rates of change (RoC) in all models are generally close to zero, reflecting an overall balance between growth and decay of field strength.The exception is gufm1 which despite its better temporal resolution lacks a long enough time span to achieve a robust average.The average RoC of 22 nT yr 1 reflects the short term more or less continuous decrease in field strength since 1830 CE with no corresponding recovery.Looking at the most extreme values we find that the signed maximum and minimum rates of intensity and VDM changes (Table 1) show values fairly comparable to the modern field in all models except for GGF100k (lower by a factor of 2) and SHAWQ-IA (larger especially on the growth side).Global Holocene field models do not yet have sufficient resolution to fully capture rapid localized intensity variations like the LIAA (Davies & Constable, 2017) and so we suspect that the actual maximum rates of change were higher than those seen in current paleofield models.Interestingly the peak values for pfm9k.2 are quite similar to HFM.OL1.A1 despite the updated data sets in the former.Differences in RoC between HFM.OL1.A1 and CALS10k.2 reflect tradeoffs between spatial and temporal regularization in model construction and similar effects might also explain the larger extreme RoC values in SHAWQ-IA than in pfm9k.2.
For both local field and VGP positions global mean directional RoC range from 0.01 to 0.07°yr 1 across all models, and even the peak values remain low in time intervals that do not contain excursions.However, the maximum directional changes prior to the Holocene far exceed those found in the modern field, with peak values of ∂ B/∂t in GGFSS70 and LSMOD.2 exceeding 60°yr 1 and 20°yr 1 respectively.GGF100k is smoother in time than these models, but still produces angular rates that exceed 1°yr 1 , which is over twice as large as the fastest variation in the modern field.
The spatial distributions of these changes are summarized in Figure 4 which shows mean rates of change as a function of latitude for the 8 models.Rates of change of P V generally display weaker latitudinal variations than   timescales.The right panels (b, d, f, h) plot the actual values of the anomalous rates as a function of time.For directional changes (a, b, c, d) we plot only times where rates exceed the fastest variations in gufm1 0.4°yr −1 for ∂ B/dt, and 0.25°yr 1 for ∂ PV /dt.Rapid directional changes are generally coeval with the three excursions bands in gray showing Norwegian/Greenland Sea, Laschamp, and Mono Lake/Auckland in successively decreasing age and most prominent in regions where the field is anomalously weak (Davies & Constable, 2020).Interestingly, rapid low-latitude directional changes occur toward the end of each of the excursions, reminiscent of a postcursory phase of activity arising as the dipole field gradually recovers toward its pre-excursion strength.A comparison across panels (a-d) in Figure 5 reveals broad similarities in the temporal distribution of anomalous values of ∂ B/dt and ∂ PV /dt together with some variations in the latitudinal distributions, most notably in the higher latitude variations in ∂ PV /dt, as could be expected from the dipole transformation.
For the field strength (Figures 5e-5h) we retain rates with ∂B⁄∂t > 150 nT yr 1 and ∂P V ⁄∂t > 0.2 ZAm 2 yr 1 .These lower cutoffs were chosen for the data in order to see representative large values which generally lie below the peak values recorded for gufm1 in Table 1 given by the dashed purple lines.Intensity changes are not particularly fast when compared to the modern field (purple dashed line, in (f, h)), except in SHAWQ-IA.HFM.OL1.A1 has some moderately rapid changes at high latitudes near the beginning of the model and pfm9k.2 has a peak rate of increase of 170 nT yr 1 , larger than the peak positive value found in gufm1.Both pfm9k.2 and SHAWQ-IA exhibit the peak RoCs in times around the LIAA.All the Holocene models reflect changes above their average rates during recent times, but these remain lower than in gufm1.The largest pre-Holocene ∂B⁄∂t values tend to occur near the start and end of the excursions when the field is weakened but still substantially stronger than its local minimum.In SHAWQ-IA the largest ∂B⁄∂t values occur between 1300 and 900 BCE around the time of the LIAA; surprisingly these rapid changes are not confined to the northern hemisphere where the data record originates.There are also some intensity changes faster than gufm1 around 100-0 BCE in the SHAWQ-IA model, which occur at low-to-mid latitudes.The story for ∂P V ⁄∂t given in (g, h) is broadly similar, but with anomalous values in GGFSS70 during the Laschamp excursion that do not appear in the local field intensity RoC.
Our reported rates of change will obviously depend on the spatial and temporal sampling of the field models.
The 10 years sampling used is denser than the temporal resolution of the paleofield models.To investigate how spatial sampling affects the reported peak rates we resampled LSMOD.2 on spatial grids of 10°× 10°and 20°× 20°.Compared to the peak value of ∂ B/∂t = 27°yr 1 for the 5°× 5°grid reported in Table 1, the coarser grids yield peak rates of ∂ B/∂t = 20°yr 1 and ∂ B/∂t = 10°yr 1 respectively.A similar reduction is observed for ∂ PV /∂t.We also note that the range of values seen in Table 1 for rates of change among the various models considered here clearly highlight their differences in apparent temporal and spatial resolution and demonstrate a need for better characterization of geomagnetic RoC over the various time scales sampled.Current attempts to characterize the uncertainty in field prediction across the various models remain difficult to evaluate, and at present it is our view that the variability across models may well provide the best estimates for such uncertainties.We therefore do not ascribe particular significance to the specific values of peak rates of change.
Nevertheless, the peak rates remain significantly faster than anything seen in the modern geomagnetic field for all model resolutions considered and so we view this result as a robust feature of the multi-millennial field models.
In Figures 6a and 6b we again select rates of change in GGFSS70, LSMOD.2,GGF100k, and SHAWQ-IA that exceed 0.2 ZAm 2 yr 1 for |∂P V ⁄∂t| and 150 nT yr 1 for |∂B⁄∂t| to look for patterns as a function of B and of P V normalized by P 0 V = 70 ZAm 2 .The fastest intensity changes in (b) are associated with strong local field strength.Somewhat rapid >0.2 ZAm 2 yr 1 changes in P V are seen in (a) that occur at weak values of the local VDM, but all are <0.25 ZAm 2 yr 1 so these changes remain below the definition of anomalous based on gufm1.In contrast to RoC in field strength, the vast majority of directional changes in (c) that exceeded the fastest variations for the modern field occurred when the local VDM was less than half the 100 ky mean.Larger directional changes in (d) generally occur at lower field strengths and the very fastest changes occurred at the weakest intensity values found in the model.Taken together, these results suggest that rapid changes in direction and intensity reflect differing evolution of features of the underlying core field.
To further illustrate this point, Figure 7 shows contrasting maps of high directional RoC, ∂ B/∂t, compared with coeval intensity in LSMOD.2 (a, c) and similar maps of high intensity RoC, ∂B⁄∂t, compared to intensity in SHAWQ-IA (b, d).The time point for LSMOD.2 is chosen such that the rapid changes are located in regions that are closely associated with the underlying data locations, while the time point for SHAWQ-IA is chosen near the point of maximum ∂B⁄∂t, which is well represented by data from the LIAA.These maps provide sample snapshots in time, but similar results are obtained at other times where there are rapid changes.The fastest directional changes in (a) are indeed well correlated spatially with the minimum field strength in (c).However, the fastest intensity changes in (b), are not correlated with the strongest intensities in (d).This reflects the fact that we expect rapid changes in both intensity and direction to reflect the emergence and/or movement of magnetic features at the core surface.Detailed examination of a series of snapshots of B from LSMOD.2 indeed reveals that the maximum ∂ B/∂t is associated with the rapid growth of a region of very weak field (<50 nT) around 80°S and 30°W.The maximum changes in ∂B⁄∂t in SHAWQ-IA are associated with the westward migration of an intense normal polarity flux patch over the LIAA in the northern hemisphere.

Conclusions and Discussion
In this paper we have identified spatial and temporal geomagnetic field variations over the past 100 ky that are not present in the modern field spanning the most recent few hundred years.Previous studies have treated changes in direction and intensity separately (Davies & Constable, 2017, 2020;Korte & Constable, 2018;Maffei et al., 2021) and have generally focused on specific events such as the Matuyama-Brunhes reversal (Sagnotti et al., 2014) or the LIAA (Ben-Yosef et al., 2017;Osete et al., 2020;Shaar et al., 2016).Here we have synthesized the temporal and latitudinal variations in rapid magnetic field changes, exploiting the recent advances provided by global timedependent representations of the field over various portions of the last 100 ky.Over this time period we find that the fastest changes for direction and intensity occur under different circumstances.The fastest directional changes over the past 100 ky are associated with the recovery phase of excursions in which the dipole field slowly increases from its local minimum.These directional changes occur at all latitudes and can be up to 100 times faster than the fastest changes seen in the modern field.By contrast, the fastest intensity changes over the last 100 ky as represented in global field models are generally slower than for the modern field except in the SHAWQ-IA model that includes data from the LIAA.The fastest directional changes occur when the field strength is weak, while the fastest intensity changes occur when the field strength is comparable to or stronger than its average value over the last 100 ky.Some care is needed when interpreting the results from the global field models because of the uneven spatial and temporal sampling of the underlying data distribution.The fastest directional changes are found in GGF100k, LSMOD.2, and GGFSS70.Sensitivity kernels show that all of these models sample all parts of the CMB over their total time-spans (Panovska et al., 2019(Panovska et al., , 2021)), though some regions are poorly sampled, with the best spatial sampling generally in the northern hemisphere.Temporal sampling improves around excursions, particularly the Laschamp (Panovska et al., 2019).We therefore expect that modeled directional changes in the northern hemisphere around the Laschamp excursion are a consistent representation of the underlying data and indeed we see very rapid changes at this time in the northern hemisphere in all 3 models.Reducing the spatial sampling of the models does not change this result.The LIAA as represented in SHAWQ-IA and pfm9k.2displays a smoother temporal variation than what is seen in detailed records from the region (e.g., Ben-Yosef et al., 2017;Davies and Constable, 2017) and so we expect the intensity rates of change around this time to be under-estimates.Therefore, within the context of the available data and modeling strategies, we conclude that the rapid field changes we have found are robust.
Numerical geodynamo simulations provide independent and complementary information on rapid magnetic field changes.Davies and Constable (2020) argued that rapid directions change of >1°yr 1 were associated with migration of reversed flux patches across the CMB in their simulations.The link to reversed flux is not so clear in the field models, partly because of their reduced spatio-temporal resolution and partly because the fastest changes generally occur in and around excursions where the distinction between normal and reversed polarity field is less obvious.However, the field models do show that the fastest directional changes occur when the field is anomalously weak, as would be the case within a reversed flux patch, and when the field morphology is more "patchy."The link between excursions, field "patchiness," and rapid directional changes can be tested by future numerical simulations.Davies and Constable (2018) found rapid intensity changes of >0.5 μT yr 1 in numerical dynamo simulations that were caused by migration of normal polarity flux at the CMB.This picture is consistent with reconstructions of the LIAA (Davies & Constable, 2017;Osete et al., 2020), but on longer timescales the situation is less clear because the available field models do not produce rapid intensity changes.We suspect that intensity changes over the last 100 ky are under-estimated in the models because they are associated with relatively small scale features on the CMB that are hard to resolve with the available data.The development of more detailed models of regional field changes may help to address this issue.
Despite the aforementioned limitations, global field models spanning the past 100 ky show that rapid changes in intensity and direction are associated with strong and weak field respectively, strongly suggesting that they reflect different features on the CMB.The expectation from the field models that rapid intensity changes occur predominantly at high latitudes and rapid directional changes occur mainly at low latitudes can hopefully spur future paleomagnetic data acquisitions to test these relations, while links to core dynamics can be probed using numerical simulations.Indeed, we suggest that these rapid changes, which are not seen in the modern magnetic field, should be regarded as an essential ingredient for geodynamo simulations to reproduce in order to be considered "Earth-like."Existing criteria that assess temporal variability in simulations either rely on high-resolution observations of the secular variation (Mound et al., 2015), frequency-dependence of the dipole moment (Davies & Constable, 2014), or the My timescale variations recovered by syntheses of paleomagnetic data that provide relatively sparse spatio-temporal sampling (Sprain et al., 2019).However, the global field models considered in this paper that access intermediate timescales reveal temporal fluctuations, including local rapid intensity and directional changes, that are distinct from those seen on centennial or Myr timescales and cannot be assessed using the aforementioned criteria.We speculate that the PSV activity index (Panovska et al., 2018), which measures both intensity and directional variations in the field and can be calculated locally and globally, could be used to consistently define temporal geomagnetic field variability across these broad ranges of timescales, enabling systematic comparisons with numerical geodynamo simulations.

Figure 1 .
Figure 1.Power spectrum showing variations in the paleomagnetic axial dipole component as a function of frequencymodified from Figure 2 of Constable and Constable (2023) who provide details about the individual contributing models.The black arrows on this figure correspond to transitions in the approximate exponent of the power law in frequency given by the dashed lines.

Figure 2 .
Figure 2. Power spectrum of time rate of change in the axial dipole component of the geomagnetic field computed from the grand spectrum of Constable and Constable (2023) available at https://earthref.org/ERDA/2719/.Gray bars give one sigma uncertainty estimates.

Figure 5
Figure5exhibits the temporal distributions of anomalous rates of change.The left panels (a, c, e, g) give the latitudes at which anomalous values occur as a function of time for the 7 models that span Holocene and longer

Figure 5 .
Figure 5. Latitudinal distributions of anomalous rates of change as function of time (left) and anomalous values as function of time (right).For directions (a, b) shows B;(c, d), PV .For directional rates only locations and times exceeding the corresponding maximum value for the modern field (dashed purple lines on right-hand plots, at 0.4°yr 1 for d B/dt, and 0.25°yr 1 for d PV /dt) are plotted.For field intensity (e, f) shows B and (g, h) shows virtual dipole moment, P V .We used different intensity criteria, retaining rates with ∂B⁄∂t > 150 nT yr 1 and ∂P V ⁄∂t > 0.2 ZAm 2 yr 1 (lower than the peak values in gufm1-indicated by dashed purple line in f, g).The Norwegian/Greenland Sea, Laschamp and Mono Lake/Auckland excursions are shown by gray shading with start and end times as determined byPanovska et al. (2019).Negative years for time reflect BCE and positive are CE.

Figure 7 .
Figure 7. Mollweide projections of ∂ B/∂t in (a) compared with intensity, B in (c), at snapshot of high ∂ B/∂t in LSMOD.2;For SHAWQ-IA intensity RoC (∂B⁄∂t) in (b) compared to intensity, B, in (d) at the time of maximum ∂B⁄∂t.

Figure 6 .
Figure 6.Large values of RoC for (a) VDM and (c) VGP direction versus normalized VDM, P 0 V .Black vertical lines denote half the mean VDM of 70 ZAm 2 used as P 0 V .Over the past 100 ky lower values have only occurred during excursions.(b) |∂B/ ∂t| and (d) ∂ B/∂t versus B at the same locations.For directions, only rates of change larger than those in gufm1 are shown, while for the intensities the cutoffs are ∂B/∂t > 150 nTyr 1 and ∂P V /∂t > 0.2 ZAm 2 yr 1 .

Table 1
Rates of Change Statistics in Geomagnetic Field Models: Mean Is the Time-Averaged Rate of Change (RoC), av abs Value Is the Time-Averaged Absolute Value of RoC, Mad Is Value of the Mean Absolute Deviation, Max (Min) Is the Maximum (Minimum) Value of RoC Over the Model Time Interval, Skew Gives the Skewness Value for the RoC Distribution Bold face text indicates values larger than the thresholds defined based on gufm1.
If we consider the unsigned RoC (listed as av abs value in the table) the smallest mean values in |∂B⁄∂t| and |∂P V ⁄∂t| occur for GGF100k which has the smoothest representation of field behavior considered here.Rapid intensity changes in that model are erased by smoothing arising from limitations in the available chronological constraints.GGFSS70, LSMOD.2,CALS10k.2,HFM.OL1.A1, and pfm9k.2exhibit intermediate mean RoCs, while SHAWQ-IA and gufm1 are both considerably larger.It seems probable that the values seen for SHAWQ-IA remain underestimates because of the lack of spatial and temporal data coverage during the LIAA.Both CALS10k.2 and HFM.OL1.A1 data compilations predate the LIAA data acquisition -and have much lower mean RoC for both |∂B⁄∂t| and |∂P V ⁄∂t| as does pfm9k.2which does include these data, partly a reflection of the short time interval with high RoC spanned by SHAWQ-IA.