One Hundred Thousand Years of Geomagnetic Field Evolution

Paleomagnetic records from sediments, archeological artifacts, and lava flows provide the foundation for studying geomagnetic field changes over 0–100 ka. Late Quaternary time‐varying spherical harmonic models for 0–100 ka produce a global view used to evaluate new data records, study the paleomagnetic secular variation on centennial to multimillennial timescales, and investigate extreme regional or global events such as the Laschamp geomagnetic excursion. Recent modeling results (GGF100k and LSMOD.2) are compared to previous studies based on regional or global stacks and averages of relative geomagnetic paleointensity variations. Time‐averaged field structure is similar on Holocene, 100 ky, and million‐year timescales. Paleosecular variation activity varies greatly over 0–100 ka, with large changes in field strength and significant morphological changes that are especially evident when field strength is low. GGF100k exhibits a factor of 4 variation in geomagnetic axial dipole moment, and higher‐resolution models suggest that much larger changes are likely during global excursions. There is some suggestion of recurrent field states resembling the present‐day South Atlantic Anomaly, but these are not linked to initiation or evolution of excursions. Several properties used to characterize numerical dynamo simulations as “Earth‐like” are evaluated and, in future, improved models may yet reveal systematic changes linked to the onset of geomagnetic excursions. Modeling results are useful in applications ranging from ground truth and data assimilation in geodynamo simulations to providing geochronological constraints and modeling the influence of geomagnetic variations on cosmogenic isotope production rates.


Introduction
The internal magnetic field generated by dynamo processes in the fluid outer core is one of Earth's fundamental properties. It sets the stage for interactions with the solar wind plasma and likely protected the Earth from early volatile losses suffered by other planets in the solar system, thereby preserving the atmosphere necessary for life on Earth. It shields the planet from charged particles originating mostly from the Sun and moderates the flux of cosmic rays and ultraviolet radiation, which control cosmogenic isotope production in the upper atmosphere. Relatively high modern geomagnetic field strength supports these benefits, but very large changes in the form of past geomagnetic reversals and excursions accompanied by drastic reduction in field strength may have produced large changes in the structure of Earth's magnetic environment.

10.1029/2019RG000656
Aside from intrinsic scientific interest in the origin and physics of the geodynamo it is thus important to understand in as much detail as we can the nature and origin of geomagnetic variations.
This review, covering current knowledge of geomagnetic field evolution over 0-100 ka, is intended to target a wide range of audiences. For paleomagnetists and geomagnetists, a global synthesis of knowledge of 0-100 ka geomagnetic variations provides a more complete view of spatial and temporal magnetic variations at Earth's surface, and if downward continued to the core-mantle boundary (CMB), it provides a closer view of the core processes responsible for secular variation. Statistical representations of paleosecular variation (PSV), including average field, variance, power spectra, and symmetry properties, are needed by the geodynamo community for validating numerical simulations. In geomagnetic data assimilation studies (e.g., Fournier et al., 2010;Gillet, 2019;Sanchez et al., 2019), output from time-varying field models based on paleomagnetic data can be used to provide additional constraints for numerical geodynamo models. Local predictions of time-varying geomagnetic elements can be used for geochronological purposes such as paleomagnetic dating of sediments, lava flows, archeomagnetic artifacts, and stalagmites by matching the paleomagnetic record with model predictions. Potential applications also extend to cosmogenic isotope studies. 14 C, 10 Be, and 36 Cl, are produced in the Earth's upper atmosphere by nuclear interactions between energetic cosmic ray particles with target elements. This cosmogenic isotope production is directly modulated by the geomagnetic field as was initially inferred from archeomagnetic absolute paleointensities (Elsasser et al., 1956) and by the solar activity (Beer et al., 1990). Separating these influences remains an active area of study.
Paleomagnetic records from several sources (volcanics, archeological artifacts, stalagmites, and sedimentary materials) that serve as proxy magnetometers provide access to geomagnetic field evolution before the age of systematic ground and satellite measurements or historical observations of Earth's magnetic field. The spatial and temporal coverage, accessible detail, and temporal resolution attainable are diminished as one goes further back in time. Although several excellent review articles describe the state of geomagnetic field evolution during the Holocene (see, e.g., , and others have focused explicitly on geomagnetic excursions (Laj & Channell, 2015), or intervals with polarity transitions (Valet & Fournier, 2016), the longer 0-100 ka time interval has not received comprehensive attention. In part this is because the paleomagnetic data basis for 0-100 ka with associated chronological constraints has only recently become sufficient for global field modeling. However, the time interval is an important one, as it is long enough to span an interval containing several geomagnetic excursions whose global extents are not yet fully understood, but not so long that poor data coverage and temporal resolution greatly reduce the capacity to make inferences about global properties of the magnetic field. Previous discussions of global/regional paleointensity stacks (e.g., Laj et al., 2000;Laj et al., 2004;Stoner et al., 2002) have focused more on the data collected than the current review, which will deal with presently available compilations of high quality paleomagnetic records (more than 150 sediment records) and a new global, time-dependent geomagnetic field model dominated by sediment records and extending over the past 100 ka.
The layout of this review is as follows. Bearing in mind the breadth of our anticipated audience, we introduce a number of paleomagnetic terms, commonly used abbreviations, and concepts in section 2. Then we describe the spatial and temporal coverage provided by the various available data and compilations that provide information about 0-100 ka geomagnetic field evolution (section 3). Data syntheses are an important part of geomagnetic field research providing insight into global field characteristics and dynamo processes in the outer core. We give an overview of synthesis methods in section 4 and describe results from published regional and global syntheses of the observations in section 5. These have evolved over time from stacking and averaging regional and global records into spherical harmonic (SH) reconstructions of the global geomagnetic field evolution over the past 100 ka, using variations on modeling strategies developed for modern field and the Holocene interval (roughly 0-10 ka; section 5.2). Field properties are discussed in the global context in section 6 in terms of various field characteristics at Earth's surface and the CMB. In section 7, we focus on excursions, with special attention to the Laschamp excursion, and to the signatures of other suggested excursional events in the global model over the past 100 kyr. We discuss characteristic geomagnetic field properties, which can be retrieved from global empirical models and may be used in the validation of geodynamo simulations or in data assimilation (section 8). Finally, in section 9 we summarize recent progress in understanding geomagnetic field evolution over the past 100 ka and suggest perspectives and implications for other general applications.
If the geomagnetic field were a pure dipole, individual local field observations of direction and intensity would provide information about global field strength and dipole tilt. Thus, transforming local observations from a site with geographic latitude , longitude to equivalent dipole properties is commonly used to remove the largest geographic variations in these properties. Geomagnetic colatitude, g , and latitude g (with g = 90 • − g ) are referred to an equivalent dipole axis instead of Earth's rotation axis, and can be obtained from an inclination measurement, I, via tanI = 2tan g . (1) A virtual dipole moment (often abbreviated as VDM) is obtained from an intensity (F) measurement at geomagnetic colatitude g by with magnetic vacuum permeability 0 and Earth's radius a. If inclination is not known, a virtual axial dipole moment (VADM) can be obtained using geographic latitude, , instead of g in equations (1) and (2). Here an even stronger assumption (and correspondingly less accurate approximation at any given instant) would be that the field is a pure dipole oriented along the rotation axis. Geographic latitude, V , and longitude, V , of the virtual geomagnetic pole (VGP) can be obtained through geometrical considerations from an observation of declination, D, and inclination, I, at a known site location ( , ) (see, e.g., Merrill et al., 1996, for more details).
Individual VADMs, VDMs, and VGPs can deviate substantially from the DM and geomagnetic pole, depending on the influence of non-axial-dipole or nondipole field contributions. A common paleomagnetic assumption is the geocentric axial dipole (GAD) hypothesis, which says that over long times the geomagnetic field variations average to a centered dipole aligned with Earth's rotation axis. It is thus often assumed that nondipole field contributions average out if V(A)DMs or VGPs can be averaged over a range of locations and long times. However, the validity of this assumption is unclear. Merrill et al. (1996) report that time averages of 10 4 to 10 5 years have been considered sufficient, while studies of the time-averaged field (TAF) frequently find deviations from GAD on timescales ranging from millennia up to 5 Myr (see, e.g., Constable et al., 2016;Cromwell et al., 2018, and Sections 5.3, 6.4).

PSV and the PSV Index P i
The continuous variation of the paleomagnetic field is called PSV. Panovska and Constable (2017) introduced the PSV index, P i , that can be used to measure PSV strength and variability in either a geographical or temporal context, or both. At a given location ( , ) and time t the nondimensional index depends on VGP latitude, V , and the VDM, and is scaled by the present-day DM of DM 0 = 80 ZAm 2 : As defined in equation (3), P i is a combined measure of local deviation from an axial dipole field direction and from modern dipole strength. The standard deviation of the index over time is a measure for field variability in the time interval and is called PSV activity. The PSV index can be determined for data records, empirical geomagnetic field models or numerical dynamo simulations, and can be used to compare field variability regionally, globally, and in time .
Acquiring better understanding of excursions is complicated by the fact that a directional signal of an excursion might not be detected even in a sedimentary record that is supposed to be continuous, due to inherent smoothing of the signal from low sedimentation rates or from postdepositional remanent magnetization acquisition (Roberts & Winklhofer, 2004). Likewise the decrease in paleointensity may be muted in sedimentary environments. Over the past few decades a lot of individual excursional observations have been reported, both from sediments and volcanic sequences (see Laj & Channell, 2015). In order to be commonly accepted, excursions need to be documented from different regions and with good age control. In fact, uncertainties in ages can be the major obstacle to determine if two independently observed excursional signatures belong to the same or two different events. Even the ages of well-documented excursions such as the Laschamp (now ∼41 ka) have been notably revised over time. In that case, age estimates have ranged from 50 ka to between 8 and 20 ka (Laj & Channell, 2015;Singer, 2014), and it was disputed whether the Laschamp and Mono Lake (now ∼33 ka) excursions were indeed two different events until both were found in the same records. A number of excursions are generally accepted, and several additional ones have been proposed but await confirmation by additional data or rejection through age revisions or for reasons of recording mechanism. The review by Laj and Channell (2015) lists seven well-constrained global excursions plus several less certain ones for the Brunhes Chron, that is, the normal polarity epoch lasting since the most recent reversal that occurred ∼780 kyr ago. We discuss the excursions falling in the recent 100 ka interval in section 7 (see Figure 16 for an overview of their names and ages).
Little is known about the driving mechanism and the physics of excursions and their relation to field reversals. Excursions might be due to particularly strong (dipole) secular variation (e.g., Lund, 2018;Lund et al., 2005), or they, perhaps together with reversals, might individually have distinct driving mechanisms . It has often been argued that excursions are aborted reversals (e.g., Cox et al., 1975;Hoffman, 1981;. Gubbins (1999) suggested that this happens when the field reverses polarity only in the outer core, without reversing in the solid inner core. Three main descriptions of the excursion and reversal process have been proposed (Merrill & McFadden, 1994): (1) a significant rotation of the dipole axis away from Earth's rotation axis, accompanied with a decrease in nondipole power to explain the observed intensity low (e.g., Lund et al., 2006;Tric et al., 1991); (2) a decrease of the (axial) dipole field so that nondipole field becomes dominant (e.g., Constable, 1990;Courtillot et al., 1992;; and (3) localized perturbations in Earth's outer core that cause regional excursions with strong nondipole field, but no dipole reversal. It is widely discussed whether (all) excursions show preferred transitional VGP paths or not (see Laj & Channell, 2015;Roberts, 2008), which could indicate the respective relevance of dipole and nondipole field contributions to the excursion and/or an influence of lower mantle structure on the geodynamo process (Gubbins, 1994;Gubbins & Coe, 1993;Laj, 1991). Comparisons with numerical dynamo simulations showing reversals and excursions  can be instructive but require a robust empirical understanding of the global excursion dynamics, a topic discussed in section 7.

Data on the Past Geomagnetic Field
Spatial and temporal variations of the geomagnetic field have been observed directly over the past 500 years via satellites, geomagnetic observatories, and historical measurements (see, e.g., Hulot et al., 2015;Jackson & Finlay, 2015). Over geological timescales, knowledge about past field evolution relies on indirect measurements based on remanent magnetization in lake and marine sediments, volcanic rocks and archeological artifacts (see, e.g., Laj & Channell, 2015). Volcanic and archeomagnetic materials acquired their magnetization via a process called thermoremanent magnetization, TRM (see, e.g., Tauxe, 2002). Basically, individual magnetic moments that freely fluctuate above the Curie temperature become statistically aligned with the geomagnetic field below the Curie temperature, when the material cools down. If the sample orientation is known, archeomagnetic and volcanic material can provide information about three absolute field components: declination, inclination, and paleointensity. The directional angles are commonly determined by principle component analysis after thermal or alternating field demagnetization. Paleointensities have to be determined by experiments that assume that a TRM acquired in the laboratory is proportional to the remanent magnetization from the past magnetic field. Several methods to determine paleointensities exist, and the failure rates of all the experiments are high. Moreover, in order to verify that paleointensity measurements are reliable, a range of tests need to be performed to check the thermal stability and detect any possible chemical alterations of magnetic remanence carrying minerals (see, e.g., Tauxe & Yamazaki, 2015;Dunlop & Özdemir, 2015). Records from archeological artifacts and volcanic rocks consist of individual, often sparse snapshots for particular locations and specific times, where the latter must be independently identified.
Sediments acquire their magnetization during deposition when the magnetic moments of grains are statistically aligned with the geomagnetic field, the process called (post) depositional or detrital remanent magnetization, DRM (Johnson et al., 1948;King, 1955;Tauxe, 2002). Sediments offer the possibility for long and continuous paleomagnetic records and sufficient global distribution and are therefore highly valued in reconstructing long-term geomagnetic field variations. Typically, they provide two components in relative form, declination and relative paleofield intensity (RPI), while inclination is commonly supposed to be absolute, unless there is evidence for nonvertical entry during coring. The intensity of the natural remanent magnetization must be normalized using some parameter based on magnetic mineralogy, to remove variations unrelated to changes of Earth's magnetic field. Tauxe (1993) suggested criteria that paleointensity measurements need to fulfill to be considered reliable. Because most cores are not oriented during drilling, declinations are generally relative and usually presented as zero-mean values. The two relative components can be calibrated with the help of absolute data derived from volcanic rocks and archeological artifacts, or using predictions from existing global or regional models based on absolute measurements (Korte & Constable, 2006). Sedimentation rate plays an important role in the geomagnetic signal recorded in sediments. Low sedimentation rates can significantly smooth the signal, and some features may be completely absent from the record (Roberts & Winklhofer, 2004;. The resolution also depends on whether the sampling and measurement strategy uses discrete or u-channel samples. Sampling using u-channels can provide fast and dense measurements, but the signal is attenuated by the magnetometer response function that suppresses recording of large paleomagnetic changes (e.g., Philippe et al., 2018). Additionally, it can be difficult to deal with the end effects from sectioning cores. Deconvolution may be employed on u-channel data to produce higher-resolution records (e.g., Constable & Parker, 1991;Guyodo et al., 2002;Jackson et al., 2010;Oda et al., 2016). Roberts (2006) discussed the strengths and limitations of u-channels and suggested strategies for mitigating the limitations of these measurements. In field modeling, however, it is preferable to take account of the instrument convolution in the forward modeling thereby directly predicting the smoothed observations .

Data Uncertainties
Many factors, ranging across recording processes in specific media, sampling, and lab protocols, influence the accuracy of paleomagnetic and volcanic/archeomagnetic data, and these data have much higher uncertainties than direct field measurements. It is crucial for many applications (geochronology, field modeling, etc.) to acquire a reliable estimate of the uncertainty in the magnetic field record and a variety of approaches has been used. Often, minimum uncertainties are assigned for different classes of data, based on broad inferences. For the Holocene period, sediment data have been compared with historical models (Jackson et al., 2000) over the past 400 years to provide an assessment of uncertainties (Constable et al., 2000). The sediment records used in the CALSx Holocene models propagated uncertainties from the 95 confidence circle of the direction (Fisher, 1953) of 6 • for directions into declination and inclination (see, e.g., Donadini et al., 2009), and estimated 5 μT for intensities (e.g., . Suttie and Nilsson (2019) recently found an error in the relationship between 95 and directional uncertainties. Apparently, the commonly used conversion produces a higher uncertainty by a factor of √ 2 (for example, 95 = 6 • should convert to an uncertainty in inclination of 2.5 • instead of 3.5 • ). However, as paleomagnetic uncertainties mostly are not known in a strict statistical sense, this only translates to a slight shift in relative weighting of directional versus intensity data in most of the currently applied global modeling methods. Individual uncertainty estimates have been derived for sediment records covering the past 10 and 100 ka by analyzing smoothing spline fits to the records (Panovska et al., 2012;Panovska, Constable, & Brown, 2018). These results show a wide range of uncertainties for the three components, with inclination being the most reliable. The analyses also underline the range of quality across the records and the importance of appropriate treatment when used to reconstruct the geomagnetic field.
Volcanic and archeomagnetic data are often obtained by averaging measurements from several samples or specimens and the uncertainties are reported in 95 values for directions and as standard deviations or standard error for intensities. Archeomagnetic/volcanic data usually have a priori uncertainties assigned and published in the original studies, in contrast to the majority of sediment studies. In all cases when uncertainties are not available, or are smaller than a certain plausible threshold, it has become customary to allocate minimum uncertainties (e.g., Panovska, Constable, & Korte, 2018). Other methods include introducing a modeling error and adding it in quadrature to the original data uncertainty (Licht et al., 2013), assigning minimum errors based on the number of samples/specimens used to calculate the mean (Nilsson et al., 2014), grouping the data in different categories based on paleomagnetic quality criteria (Korte et al., 2005;, and via comparison with historical data, considering error estimates for the difference (Arneitz et al., 2017).
Dating uncertainties are another important source of uncertainty in sediment, volcanic, and archeomagnetic records. Depending on the dating method, age uncertainties span a wide range of values, from a few decades to a few centuries, or even millennia. (Panovska, Constable, & Brown, 2018) summarized all dating methods used for sediment, volcanic, and archeomagnetic data spanning the past 100 ka. The following methods are most commonly used for sediments: radiocarbon ( 14 C), oxygen isotopes ( 18 O), and correlation based on magnetic properties, for example, magnetic susceptibility, for transferring the age scale of a dated record to parallel cores. When only a few tie points are used to assign ages for a whole core, then constant sedimentation rates are assumed. This may not always be correct, and data lying between tie points can have significantly larger age uncertainties. Geochronological methods for archeomagnetic data are discussed in (Brown, Donadini, Korte, et al., 2015). Limitations of radiocarbon and Ar/Ar dating can diminish accuracy of dating for volcanic records older or younger than 50 kyr, respectively, making it difficult to accurately tie them into the global record. (Panovska, Constable, & Brown, 2018) recently analyzed a global set of all available published paleomagnetic data for the past 100 kyr, considering two types of records: continuous time series from lake and marine sediments, and discrete points in time from lava flows and archeological artifacts. Sediment data were compiled from the following sources: the SEDPI06 collection (Tauxe & Yamazaki, 2015), GEOMAGIA50.v3 (sediment part) , Pangaea (Diepenbroek et al., 2002), the MagIC database https://www.earthref. org/MagIC, records longer than 10 ka from the Holocene compilation of , and new sediment records provided directly by several authors. Volcanic and archeomagnetic data were extracted from the GEOMAGIA50.v3 database Brown, Donadini, Korte, et al., 2015) for the period 0-50 ka and the global 0-10 Ma PSV10 data set (Cromwell et al., 2018) for the period 50-100 ka.  Panovska, Constable, and Brown (2018), and archeomagnetic and volcanic data from the GEOMAGIA50.v3 (Brown, Donadini, Korte, et al., 2015) and PSV10 (Cromwell et al., 2018) databases. Distributions of studied regions in publications from the GEOMAGIA50.v3 and PSV10 databases (b) and the sediment data compilation (c).

Data Compilations
In Figure 1a we present a summary of numbers of publications and database sources for sediment, volcanic, and archeomagnetic data covering 0-100 ka. The data collection includes 121 studies from the sediment data compilation (listed in the supporting information of Panovska, Constable, and Brown (2018)), 728 from the GEOMAGIA50.v3 and 81 from PSV10. In general, the number of studies increases toward the present, with an apparent decrease in recent years because many newer studies were not available in the data compilation. If these were all included, the general trend would reflect increased attention to the production of paleomagnetic records, although PSV is not always the primary goal of every study. The earliest study of archeomagnetic data used is the pioneering work of Thellier and Thellier (1959), while in the sediment compilation Kent and Opdyke (1977) presented intensity variations in the core RC10-167 from the North Pacific over the Brunhes epoch.

Data Coverage
Our review is focused on the global view on geomagnetic field evolution, and all findings are based on global or regional syntheses of data, and not individual records. Therefore, only records that have previously been used in such syntheses are considered here, and some more recently published records are not included. The geographic distribution of the GEOMAGIA50.v3 (Brown, Donadini, Korte, et al., 2015) and PSV10 (Cromwell et al., 2018) databases and of the sediment data compilation of (Panovska, Constable, & Brown, 2018) are given in Figures 1b and 1c. The European region is covered best by volcanic and archeomagnetic data (60%), followed by Asia and North America each with 13%. Other regions are poorly represented, 5% of the data come from Central America, and 3% from each of Africa and South America. Deep-sea sediment records from the North Atlantic and North Pacific (58% combined) prevail over the South Atlantic, South Pacific, Southern, and Indian Ocean (13%). Twenty-four percent of the data are lake sediment records that cover the Holocene in addition to part or all of the 100 ka interval, for example, Lac du Bouchet, France (Thouveny et al., 1990). The remaining marine sediments segment ( Figure 1c) contains records from the Black Sea (Nowaczyk et al., 2012;Nowaczyk et al., 2013), the Mediterranean (Tric et al., 1992), and Adriatic and Ionian Seas (Vigliotti, 2006). A single record (Pan et al., 2001) in the present compilation comes from  (Panovska, Constable, & Korte, 2018). Record names are omitted. Sediment cores are plotted in an alphabetic order, as listed in supporting information Table S1 in Panovska, Constable, and Brown (2018). Lowermost, shorter records of each component are drawn from the Holocene sediment compilation. Top plot shows the total number of data in 1,000-year bins. a different kind of sediment archive for paleomagnetic secular variation, namely, the well-known Chinese loess deposits. Chinese loess-paleosol sequences are mainly used in paleoclimate studies, although a few long paleomagnetic records exist, which cover the Matuyama-Brunhes transition (e.g., Yang et al., 2010).
More details on the temporal sampling provided by the sediment collection are provided in Figure 2. Around ∼60% of the records span 50 ka or longer and 45% and 40% of records cover ≥80 ka and ≥90 ka, respectively. The number of data points in the individual records varies, influencing the temporal resolution or smoothing time of the time series. (Panovska, Constable, & Brown, 2018) show that more than 70% of the records have a temporal resolution better than 1 kyr, but a few low-resolution records have smoothing times of up to 6 kyr. Over the whole time interval, the highest numbers of directional data are found in the Holocene period and during the Laschamp excursion (∼41 ka BP). The number of intensity data increases steadily toward the present, with a small peak over the Laschamp and a decrease over the past 10 ka. Many of the intensity records come from ocean environments rather than lakes, which dominate the Holocene, and are much longer (sometimes million of years), with lower sedimentation rates compared to Holocene records. The absence of the most recent part of the record in many of these cores also makes their contribution to the past 10 ka smaller.  Panovska, Constable, and Brown (2018). Maps represent probability distribution functions on a log scale. The more data, there are the higher the concentration is (red colored locations). (b) Data kernels representing averaged sampling of the field at the core-mantle boundary over the past 100 ka. These will vary for different time intervals across the data distribution. Note different color scales for the plots in (b). Combined kernels depend on the number of data since they sum the magnitude of all contributions from the available sampling sites. Differences in the color scale between kernels of sediments and volcanic/archeomagnetic data reflect the small number of data in the second data set. See Figure 4 for an example of data kernels for a single location. Spatial data distributions of the sediment, volcanic, and archeomagnetic data covering the past 100 ka are presented in Figure 3a. The Northern Hemisphere is better covered than the Southern Hemisphere for both types of data. Africa, Southern Atlantic and Pacific latitudes, Indian and Southern Ocean are poorly sampled. However, new sediment records have been published since this data set was compiled, for example, from the South Atlantic mid-ocean ridge (Channell et al., 2017), or the Congo deep-sea fan , and the geographical distribution is continually improving.
The consequences of the spatial distribution of observations can also be understood in terms of data kernels which quantify, via the relevant Green's function, the extent to which a measurement at the Earth's surface is sensitive to changes in the radial magnetic field at the CMB. The sampled region on the CMB is much broader, because of the effects of upward continuation of the field through Earth's mantle. Since declination, inclination, and intensity components are nonlinearly related to the radial magnetic field at the CMB, the kernels are linearized about the axial dipole field (Johnson & Constable, 1997). Details about the Green's function for surface magnetic field from the radial CMB field can be found in Constable et al. (1993). An example of data kernels (which vary with location) is provided for a vector field observation in the Black Sea in Figure 4. Basically, a single declination observation in the Black Sea samples a broad longitudinal area of the CMB (peaking at a distance of 22 • ), and inclination and intensity are most sensitive to areas beneath the measurement point, but at slightly lower and higher latitudes, respectively. When all three components of the field are available the summed data kernel is as given in the rightmost panel of Figure 4.
When summed kernels for the entire data distribution are given by type in Figure 3b, we can see the overall impact of the geographic sampling in part (a) of that figure. The temporally averaged sensitivity kernels are estimated for the different data types and for the overall data set. No region is completely unrepresented, but if a model covering the past 100 ka is based only on volcanic and archeomagnetic data set, the entire Southern Hemisphere is insufficiently sampled. The maximum value of the combined summed kernels and thus the best-sampled region on the CMB occurs at midlatitudes in the Atlantic. Sediment records in the Pacific and Southern Hemisphere contribute significantly to better constrain a model in these areas. Note that one should keep in mind that these data kernels correspond to the whole 100 ka period, and the sampling varies from epoch to epoch.

Data Synthesis Methods
Geomagnetic field evolution is governed by global processes in the fluid outer core, yet the field variations have regionally distinct signatures at the Earth's surface. Long-term (multicentennial to multimillennial) trends in paleointensity are generally supposed to follow global DM variations, but shorter-term modulations and directional changes are largely dominated by the nondipole field and thus may not be clearly correlated across different regions. Syntheses of data covering larger regions or the whole globe are used to gain insights into the underlying global processes from localized data records. In this section we describe three commonly used methods to derive such syntheses. These are global or regional stacks of records from distinct locations, global geomagnetic field models, and time-averaged global field models. With all methods one has to keep in mind the underlying data distribution. The global dipole field variation might be well recovered in global averages, if regional nondipole field variations cancel out with a good data coverage. Geographical biases in data distribution, however, might influence both stacks and global models and lead to leakage of nondipole field into the dipole reconstruction (e.g. Genevey et al., 2008). We provide an overview of the results using these methods for time intervals ranging from 10 kyr to 5 Myr in section 5.

Stacks
As outlined in section 2 local field data are frequently translated into equivalent virtual (axial) dipole moments (VADM or VDM) or VGP locations ( V , V ) to remove gross geographic variations represented by a dipolar field. Because directional and intensity data are not necessarily linked in the collection process, data stacks for the DMs and VGP locations are usually independent entities, which do not take account of any covariance between the two properties of the field. 4.1.1. VADM Reconstructions Using Absolute Paleointensities Global (axial) dipole moment reconstructions have been obtained by scalar averaging of V(A)DMs over the globe and over certain time intervals to ensure reasonable spatial data coverage (e.g., Korhonen et al., 2008;McElhinny & Senanayake, 1982). It is generally assumed that nondipole field contributions from the different locations are thereby averaged out and the global reconstruction of VADM(t) or VDM(t) evaluated at distinct time intervals t i gives a reasonable proxy for the actual geomagnetic DM variations.

RPI Stacks
A similar approach has been adopted for RPI records where is assumed that multimillennial RPI records mainly reflect axial dipole strength variations, so that global stacks represent them well. Regional stacks have been constructed for areas with particularly good data coverage or better temporal resolution than is globally obtainable. Since regional stacks contain nondipole contributions, a comparison between global and regional reconstructions must be considered with care (Genevey et al., 2008). Individual RPI records have specific temporal sampling intervals, depending on sedimentation rate and experimental technique (discrete or u-channel sampling), and scaling. In particular, if individual records have no age model assigned yet, or there are large uncertainties in existing age models, records to be stacked are often temporally correlated first, for example, by using the measured susceptibility or 18 O (e.g., Guyodo & Valet, 1999) but sometimes also directly using the RPI variations (e.g., Laj et al., 2004) or specific events such as excursions and reversals. The records then are interpolated and resampled in common intervals, and they have to be normalized to unify them in terms of mean and amplitude of variations (e.g., Guyodo & Valet, 1999;Laj et al., 2004). A stack is commonly obtained as arithmetic mean of the resampled values. Uncertainty estimates may be given in the form of standard deviation or standard error from the averaging (e.g., Yamamoto et al., 2007;Yamazaki & Oda, 2005), or by bootstrap resampling of variations in the data set (e.g., Laj et al., 2000;Xuan et al., 2016). Occasionally, methods like iterative outlier rejection are used to reduce potentially unreliable variability in the stack (e.g., Laj et al., 2004). Several RPI stacks have been calibrated to absolute values by comparison to globally averaged volcanic and archeomagnetic VADMs for some time interval (often not spanning the full length of the stack) (e.g., Guyodo & Valet, 1999;Valet et al., 2005) and sometimes including additional scaling by minimum RPI values, for example, during the Laschamp excursion (e.g., Channell et al., 2009;Laj et al., 2004). Ziegler et al. (2011) developed an inverse modeling approach, intended to minimize the bias due to non-axial-dipole contributions, first for global and later (Ziegler & Constable, 2015) for regional axial dipole moment reconstructions that are continuous in time. The method combines information from sediment RPI and igneous rock absolute paleointensities and uses a penalized maximum likelihood approach to trade-off data misfit against complexity in an axial dipole field reconstruction. The temporal parameterization used for the resulting model is in the form of cubic B-splines. The absolute paleointensities are first fit by a penalized spline of very low temporal resolution, and this model is used to calibrate all RPI records. All data are then refit iteratively until convergence is reached, using a penalized maximum likelihood method based on successively refined RPI calibrations and empirical estimates of noise distributions for the observations to find the appropriate cubic B-spline representation. As in the SH inverse methods discussed in section 4.2 the resulting model is regularized by trading off an appropriate target misfit based on noise in the observations against minimizing the second time derivative of VADM(t) to avoid spurious oscillations.

VGP Stacks
For directional information averaging or stacking VGPs globally to estimate past motions of the dipole axis has been done for (sub-)Holocene times Nilsson et al., 2011;Merrill et al., 1996), but to our knowledge has not been widely applied on longer timescales. The reasons probably are that there are too few absolutely oriented igneous rock results to determine robust global averages over time, and the lack of absolute declination values prevents the determination of robust VGP longitudes from sediments. We present a global stack of VGPs from the data compilation by (Panovska, Constable, & Brown, 2018) in section 6.

Global SH Field Reconstructions
The commonly used method to obtain the full planetary picture of geomagnetic field distribution and evolution is inverting globally distributed data for global models based on SH basis functions. The method was devised by Carl Friedrich Gauss (Gauss, 1839) in the early 19 th century and is still widely used to map and investigate the present-day field (see, e.g., Lesur et al., 2011;Gillet et al., 2010;Hulot et al., 2015) from ground and satellite observations, including the International Geomagnetic Reference Field (IGRF) (Thébault et al., 2015), which goes back to 1900 in the form of discrete snapshots in time published every 5 years. The past two decades have seen the development of a growing number of continuous longer-term global field reconstructions, where the temporal evolution is described based on continuous functions. Model gufm1 (Jackson et al., 2000) covers 400 years (1590 to 1990), relying on a large number of historical observations of declination and a few inclination measurements mainly carried out on ships for navigational purposes before systematic geomagnetic field observations became routine starting in the early nineteenth century. A growing number of SH models now spans the past 2 to 14 millennia (see, e.g., reviews by Korte et al., 2018;, using global compilations of archeomagnetic, volcanic and mainly lacustrine paleomagnetic sediment data (see, e.g., Brown, Donadini, Korte, et al., 2015;. Longer-term global field reconstructions have long remained a challenge for reasons given in section 4.2.2.

Mathematical Formulation
The fitting of global geomagnetic or paleomagnetic data by SH functions takes advantage of the fact that in a source-free region a magnetic field B can be described as the negative gradient of a scalar potential V, where t is time and r is a position vector (r, , ): that satisfies Laplace's equation This equation holds everywhere outside Earth's core when only the main field produced in the liquid outer core is considered. Neglecting the lithospheric, ionospheric, and magnetospheric field contributions is justified in paleomagnetic applications as these are too small to be resolved from the paleomagnetic signal recoverable from sediments or volcanic rocks. In spherical polar coordinates (r, , ), with r the radius from Earth's center, the colatitude and the longitude, this potential can be described by a series of associated Legendre functions P m l (cos ) of degree l and order m: where the reference radius, a, is usually taken as Earth's mean radius with a value of a = 6, 371.2 km in geomagnetism. Equation (6) and the same for h m l (t). The functions P m l (cos ) cos(m ) and P m l (cos ) sin(m ) are the SH basis functions, and a model is given by the set of {g m,k l , h m,k l } that determine the contributions of the various basis functions to describing the field, and which are determined by fitting the data using inverse methods. Predictions for geocentric field components at any location (r, , ) can be obtained from model coefficients as the component derivatives from equations (6) and (4), and converted to the commonly used components declination, inclination, and field intensity by geometrical considerations. For more details see, for example, Merrill et al. (1996), Backus et al. (1996). Under the assumption of an electrically insulating mantle the region outside Earth's core is free of magnetic field sources. The SH formulation thus offers straightforward upward and downward continuation from the Earth's surface to higher altitudes or down to the CMB by simply changing the reference radius r from a to the desired value, c, or about 3,485 km for the CMB. Continuous global SH models therefore can provide time series of field evolution anywhere on Earth, and global maps of any field component and its evolution at any altitude, in particular for our purposes at Earth's surface and the CMB.
The geometry of the spatial basis functions allows for convenient study of some field characteristics. The basis function for degree l = 1 is a pure dipole, with g 0 1 describing the contribution from a GAD, and g 1 1 and h 1 1 from two orthogonal equatorial dipoles. The dipole moment (DM), that is the strength of the best fitting, geocentric, tilted dipole is given by with 0 = 4 · 10 −7 Vs∕(Am) the permeability of free space. The dipole tilt is given by the colatitude of the dipole axis, D , via and the latitude of the dipole axis is D = 90 • − D . The full axis location including longitude, D is obtained by further geometrical considerations, tan D = h 1 1 ∕g 1 1 , see, for example, Merrill et al. (1996), Lanza and Meloni (2006).
The coefficients of SH degree l = 2 describe the contributions from a quadrupole, l = 3 from an octupole and so on. Higher degree basis functions thus describe ever smaller-scale field contributions in the spatial domain. The spatial geomagnetic power spectrum (Lowes, 1966;Mauersberger, 1956), R l , describes the average squared magnetic field strength contribution at each SH degrees or spatial wavelengths, l and is obtained by It is common to study the spectrum at the CMB (Lowes, 1974), where it is approximately white at present, except for the higher dipole contribution. The relation of energy from the (axial) dipole to that from the non-axial-dipole or nondipole contributions can be studied by summing only over the relevant degrees and orders in equation (10). The temporal basis functions M k (t) are differentiable, and it is therefore possible to determine secular variation coefficients as first derivatives of the main field coefficients. A secular variation power spectrum is defined by these along the same lines as equation (10).
The SH representation also allows the study of symmetries in the field. Coefficients of order m = 0 describe zonal, that is, rotationally axis symmetric, field contributions. The functions with even numbers for l−m are hemispherically symmetric about the geographic equator, while the ones with odd l−m are hemispherically antisymmetric. When l = m the basis function is sectorial (see, e.g., Merrill et al., 1996, for an illustration).
The spatial and temporal complexity that a model can describe is in general limited by its maximum SH degree, L, and the number of splines K. Knowing that the measurements are uncertain, the aim is not to fit the data exactly, as this will lead to a model containing unrealistically small-scale structure in time and space. However, unless the data and age uncertainties are very well known (which they in general are not), it is impossible to know how closely to fit the data. Most paleomagnetic global SH models are parameterized so that they can describe smaller scale and faster variations than the data can be expected to resolve and include regularizations in space and time, which are additional conditions that allow a trade-off between fit to the data and recovering a smooth, simple model (see, e.g., Jackson et al., 2000;Korte et al., 2009). This trade-off is accomplished by using a regularization or damping parameter on some measure of model complexity. In its usual implementation, if the spatial regularization is made extremely strong, the resulting model will tend to be a pure dipole, and if the temporal damping is very strong, this dipole will only change linearly in time. The spatial regularization is generally designed so that smaller scales are damped more than larger scales. A range of possible solutions can be explored with different damping parameters, and the choice made of a preferred model, that is considered to give the most reliable representation of the global paleomagnetic field given the information contained in the available data. This choice is based either on goodness of fit to the data or on physical considerations, whether the spectral energy distribution appears realistic in space (i.e., close to white for the large scales at the CMB, with a dropoff for the unresolved smaller scales) and not higher than seen in the present-day field in time. A simple truncation at low SH degree and small number of splines is a particularly severe and inflexible form of regularization and carries the danger of aliasing, that is, smaller-scale and faster variations that cannot be accommodated by the model are mapped into the larger scales. It has to be noted, however, that any regularization has a strong influence on the resulting paleomagnetic model and the specific method and strength of imposed regularization remain subjective choices. Other methodologies include alternative spatial constraints, for example, a dynamo norm based on the statistics of a numerical dynamo simulations with Earth-like features (Sanchez et al., 2016) or constructing the prior information from spatial and temporal statistics of the geomagnetic field, based on satellites, ground observatories, and paleomagnetic measurements, and creating temporal cross-covariance functions in which Gauss coefficients are projected (Hellio & Gillet, 2018).

Challenges in Modeling Paleomagnetic Data
The use of paleomagnetic data poses additional challenges compared to deriving SH models of the recent field, which limit the spatial and temporal resolution and also the reliability of paleomagnetic SH reconstructions compared to models from direct field observations. These are 1. Limited global data coverage 2. Lack of full vector information 3. Paleomagnetic data uncertainties 4. Age uncertainties 5. Variable resolution of time series.

Global Data Coverage
Data coverage has already been discussed in section 3.3 so we will not belabor it here but simply reiterate that despite the significant increase over time in the number of studies ( Figure 1a) several areas of the world remain poorly sampled ( Figure 3). For the Holocene, these are mainly the oceans in general, where sedimentation rates are too low to provide the desired centennial or better resolution, and the Southern Hemisphere continents. On longer timescales data mainly come from marine sediments and volcanic regions. Large parts of the Southern Hemisphere, the Asian continent and the high latitudes are scarcely covered. Predictions from SH models are clearly less reliable in areas were data coverage is sparse (see, e.g., Sanchez et al., 2016;Korte et al., 2018;Hellio & Gillet, 2018).

Lack of Full Vector Information
When the full field vector information is known the measured components can be converted to the three geocentric components B r (radial field), B (equivalent to geodetic East component), and B , which are linearly related to the potential V. The field components declination, inclination, and intensity are nonlinearly related so that the inverse problem has to be linearized and solved iteratively from a starting model. By itself this is generally no problem if the data coverage is good enough. An additional complication with paleomagnetic data is that sediment declinations and intensities are only known as relative values, and the absolute calibration has to come from comparatively few available volcanic data, especially beyond the past three millennia where a relatively large number of archeomagnetic data exist at least for parts of the world. Declination and RPI records can either be calibrated prior to modeling or during the inversion process . For prior calibration, ideally results from nearby volcanic data are used. If such are not available, previous global models can be used, keeping in mind that these too might not be well constrained for the particular region. In that case it might be as good a compromise to use a global VADM reconstruction and assume zero-mean declination over each record, that is, scaling and orienting the records under the GAD assumption. While in general an iterative solution should be independent of the starting model,  showed that this is not always the case if declination and RPI calibrations are needed as part of the inversion and absolute declination and intensity data are comparatively sparse. One particular consequence is that the dipole strength is not always as reliably determined by global SH models as one might hope.

Data Uncertainties
In contrast to studies of the present-day magnetic field, contributions from field sources other than the core are negligibly small compared to the experimental paleomagnetic data uncertainties. The latter, however, must be evaluated using one or more of the methods described in section 3.1. In many cases it may be desirable to specify a minimum expected misfit value if more specific information is unavailable. The data

10.1029/2019RG000656
are commonly weighted by their uncertainty estimates, thereby encouraging a model to fit reliable data more closely than those of lesser reliability.

Age Uncertainties
Dating uncertainties as discussed in section 3.1 are one of the reasons that temporal variability of paleomagnetic SH models (or stacked curves) is limited. If the same magnetic field variation feature is recorded by two adjacent sediment cores, but assigned slightly different ages due to uncertainties in the applied dating methods, any interpolation or modeling method will tend to smear out the feature. Conversely, if a concurrent variation in different locations is assigned different ages the difference might be described incorrectly by the model as rapidly changing regional secular variation. Prior age alignments of records as frequently done for stacks (e.g., Channell et al., 2009;Laj et al., 2004;Stoner et al., 2002) or applied regionally by Brown et al. (2018), or alignment by iterative modeling (Leonhardt et al., 2009) can mitigate the effect and can also be prone to reflect potentially wrong assumptions about contemporaneity of magnetic field variations at different locations. Several attempts at tackling age uncertainties have been made in paleomagnetic models, for example, by adjusting the sediment record chronologies by randomly stretching and compressing the individual timescales while preserving the stratigraphic relation (Nilsson et al., 2014), using a probabilistic approach (Hellio et al., 2014;Lanos, 2004) or bootstrap resampling to create an ensemble of models . All these methods can potentially be applied for global studies and beyond archeomagnetic data and Holocene timescales.

Time Series Resolution
The second reason for limited time resolution in models lies in the resolution of sediment time series (see section 3). However, advances to take into account the sediment smoothing effects are in progress. Nilsson et al. (2018) presented a method to consider postdepositional remanent magnetization acquisition effects by modeling the lock-in depth (following Roberts & Winklhofer, 2004) in the construction of age-depth models using a Bayesian method and information from high-resolution archeomagnetic field models. In order to increase the resolution of global field models in areas where high-and low-resolution sediment records exist in proximity, (Panovska, Constable, & Korte, 2018) implemented smoothing kernels adapted to each sediment record's characteristic smoothing time in the forward modeling when determining the fit to the data in the inversion. The smoothing time is determined by fitting a spline function to the record (Panovska et al., 2012) and the method was applied in constructing the GGF100k model (see section 5.2).

Time-Averaged Global Field Models
On even longer timescales, particularly the past 2.5 and 5 Myr, time-averaged global field models (TAF) also based on SHs are used to test the GAD hypothesis, seek out any persistent departures from GAD, and determine statistical properties of the geomagnetic field. The data for these models are usually global compilations of directional lava flow data (e.g., Cromwell et al., 2018;Johnson & Constable, 1996;Merrill & McElhinny, 1977;Quidelleur et al., 1994), of sediment inclinations (e.g., Schneider & Kent, 1990) or a combination thereof. Sediment-only TAF models based on inclination information alone have not proved adequate to investigate longitudinal structure (Johnson & Constable, 1995). Corrections for plate tectonic movements have sometimes been carried out (e.g., Schneider & Kent, 1990) but have mostly been considered negligible (e.g., Gubbins & Kelly, 1993;Johnson & Constable, 1995). Cromwell et al. (2018) implemented plate movement corrections on their recent PSV10 0-10 Myr lava flow data compilation and found only a small effect. The available data are usually split into normal and reverse intervals, for which separate models are derived. The available normal polarity data are mostly more numerous than the reverse polarity ones. Transitional data are generally excluded.
Two different approaches have mainly been used, either fitting a few large-scale zonal coefficients (typically dipole to octupole, that is, g 0 1 , g 2 2 and g 0 3 ) to the data (e.g., McElhinny et al., 1996;Merrill & McElhinny, 1977;Schneider & Kent, 1990), or applying a regularized inversion, typically up to SH degree and order 10, as described above (equation (6) for static coefficients without time dependence). Apart from this systematic difference the models mostly vary in terms of underlying data, where different quality selection criteria are applied and newly published data have been added over time. Differences in data weighting result from whether data are averaged temporally or binned spatially prior to modeling or not. . Data distribution of sediment records in global and regional RPI stacks, and global, spherical harmonics field models. All these reconstructions cover the whole or part of the 100 ka period. Studies are ordered by the year of publication. One of the records in the IMOLE model is from volcanic origin: Skalamaelifell record from Iceland. The Holocene period is represented by the most recent models, CALS10k.2 and HFM.OL1.A1. Several of the studies incorporate archeomagnetic and volcanic data as well: PADM2M, RADM, most Holocene models, LSMOD.1 and 2, and GGF100k. Archeomagnetic and volcanic intensity data were also used to convert some RPI stacks to VADM (or VDM). References for the stacks and models are provided in Figure 6.

Data Synthesis Results: Stacks and SH Models From 10 kyr to 5 Myr
In this section we provide an overview of results available from the various stacking and modeling approaches described in the previous section, before moving on to discuss characteristics of the paleomagnetic field over the past 100 kyr and comparisons with those from both shorter and longer intervals in section 6. Figures 5 and 6 provide the basis for our discussion, indicating the spatial coverage for sediments and temporal coverage for all types, respectively, for the data used in a wide range of different studies. Figure 6 includes models, absolute and RPI stacks longer and older than 10 ka.
The first global RPI stack (Sint-200) was produced by Guyodo and Valet (1996) for the past 200 kyr, based on 17 records from the North Atlantic, the Mediterranean, the Indian Ocean, and the western equatorial Pacific (Figure 5a). It was later extended to 800 kyr including 33 records (Sint-800, Figure 5b) (Guyodo & Valet, 1999) and 2 Myr with 15 records (Sint-2000, Figure 5d) (Valet et al., 2005). Both stacks are calibrated, with Sint-800 using volcanic absolute intensities over 0-40 kyr, and Sint-2000 the 0.8 Myr time-averaged value of 74.6 ZAm 2 . PISO-1500 (Channell et al., 2009) is a higher-resolution 1.5 Myr global stack based on 13 records (dominated by north Atlantic ones, Figure 5e) and calibrated to absolute values similarly to Sint-2000 (but note the age correlations across cores are partially based on RPI records). The shorter (75-10 ka BP) GLOPIS-75 stack (Laj et al., 2004) is based on 24 high-resolution records (>7 cm/kyr sedimentation rate with time resolution between 110 and 500 years) from North and South Atlantic, Mediterranean and Indian Ocean (Figure 5c). It is calibrated by archeomagnetic and volcanic data for 20-10 ka and Laschamp minimum intensities from French and Icelandic lavas. The 2 Myr model PADM2M (Ziegler et al., 2011), derived by a penalized maximum likelihood spline fit, is the only global axial dipole reconstruction directly incorporating absolute and RPI data. Seventy-six RPI records, 1,800 paleointensities from igneous rocks, and 3,300 archeointensities contribute to this model (Figure 5f). A recent stack (C2018-Overall) for 2-45 ka (Channell et al., 2018) combines 12 records from the Iberian margin with 12 records from other parts of the world ( Figure 5g) and was scaled to absolute VADM by Holocene global models and a value of 30 ZAm 2 for the Laschamp minimum.
Recognizing the uneven global data coverage some long-term stacks are regionally confined. NAPIS-75  is a North Atlantic stack spanning the interval 10-75 ka based on six high-resolution North Atlantic records with sedimentation rates >10cm/kyr (Figure 5h). Five sub-Antarctic South Atlantic records with 15 to 25 cm/kyr sedimentation rates are stacked for the 0-80 ka South Atlantic stack SAPIS (Stoner et al., 2002) (Figure 5i). Updated and extended stacks for the central South Atlantic (SAS-300, SACS-300, Figures 5k and 5l) and the past 300 kyr were given by Hofmann and Fabian (2007), Hofmann and Fabian (2009), who used 8 RPI records and paid special attention to improving RPI normalization for grain size and reductive diagenesis to remove environmental influences for SACS-300. A new high-resolution North Atlantic stack spanning 1.5 Ma (HINAPIS-1500, Figure 5n) based on four RPI records of >10 cm/kyr sedimentation rate was devised by Xuan et al. (2016), and a stack for 1.2-2.2 Ma (NARPI-2200) also exists for that region (Channell, Hodell, & Curtis, 2016). Yamamoto et al. (2007) presented a 250 kyr RPI stack for the northwest Pacific using 10 sediment records (NOPAPIS-250 Figure 5j), and six equatorial Pacific cores were stacked to form EPAPIS-3Ma, spanning 0.8-3.0 Ma (Yamazaki & Oda, 2005).
Regional axial dipole moment (RADM) penalized maximum likelihood spline models by Ziegler and Constable (2015) for the past 300 kyr combine absolute and RPI data and have been derived for several regions where all times are constrained by at least five RPI records (Figure 5m). In order to study various potential differences, RADMs were obtained for three longitudinal bands (Indian, Pacific, and Atlantic Oceans, respectively), two latitudinal bands (middle-high northern and equatorial latitudes, respectively, with too few data for a third southern latitude band) and three regions, two of them approximately covering the African and Pacific seismologically identified large low shear wave velocity provinces (LLSVP) (e.g., Lekic et al., 2012), respectively, and the third one assumed not to be influenced by an LLSVP. The calibration by absolute data had to be done globally due to a lack of sufficiently distributed regional absolute intensity data, and the individual RADMs have different resolution depending on the input data.

SHA Models
The most complete global picture is gained from continuous global geomagnetic field models (section 4.2). However, a large number of data, well-distributed over the globe and with sufficient age control is required to derive continuous global models. Up to now, such models have mainly been obtained for millennial-scales up to 14 kyr and a few individual excursions and reversals ( Figure 6) in the paleomagnetic context. Holocene models, their applications, limitations, and inferences are discussed in the reviews by Constable and Korte (2015), Korte et al. (2018). The latest models spanning (most of) the Holocene are named CALS10k.2 and HFM.OL1.A1 (Constable et al., 2016), based on both archeomagnetic /volcanic and sediment data ( Figure  5o), and SHA.DIF.14k , which is constrained by a geographically more limited archeomagnetic /volcanic data set only.
Due to the challenges listed in section 4.2, which increase with paleomagnetic age, only a few studies have so far attempted global SH models beyond Holocene times. The first continuous global model spanning 100 kyr, GGF100k, was recently presented by Panovska, Constable, & Korte, (2018). The model is constrained by more than 100 sediment records ( Figure 5r) and all available volcanic and archeomagnetic data (Panovska, Constable, & Brown, 2018). Uncertainties for all sediment records were estimated from the random variability in each record. Large differences in temporal resolution of the sediment records are considered through the implementation of smoothing kernels in the forward model when assessing the misfit, which increases the model resolution. Nevertheless, the temporal resolution of GGF100k is limited, not least due to potential chronological misalignment of the records. Centennial variations and in particular field excursions are not fully resolved.
Understanding the global processes driving the extreme geomagnetic field variations during reversals and excursions has motivated the derivation of SH models spanning relatively short intervals (around 10 to 20 kyr) around such events ( Figure 6). Mazaud (1995) reconstructed the upper Olduvai reversal (∼1.66 Ma) using five sedimentary records. The data are fairly well distributed over the globe, but only directional data were used and strong assumptions were made about field energy and timing of the directional changes. Shao et al. (1999) and Ingham and Turner (2008) reconstructed the Matuyama-Brunhes reversal for selected Reviews of Geophysics 10.1029/2019RG000656 snapshots in time using SH basis functions up to degree and order 3, and temporal alignment of the reversal across the data records. Leonhardt and Fabian (2007) devised a Bayesian inversion method, expanding the SH basis to degree and order 5 and implementing splines for temporal continuity. The method was applied to the Matuyama-Brunhes reversal (∼780 ka, model IMMAB4), the Iceland Basin excursion (∼188 ka, model IMIBE) (Lanci et al., 2008) and the Laschamp excursion (∼41 ka, model IMOLE) (Leonhardt et al., 2009). These three models are each based on four to five full vector records (Figure 5p), which come from different, rather uneven longitudinal distributions, and with no more than one from the Southern Hemisphere in each case. The models were obtained by an iterative procedure where a first model was devised from a single record and further records were subsequently added. The ages within each record were adjusted to the model from the previous iteration step.
A model (LSMOD.1) for the time interval 50-30 ka, including both the Laschamp and Mono Lake (∼34 ka) excursions, has been presented by Brown et al. (2018) and recently slightly updated by correction to one of the underlying data sets to version LSMOD.2 . All available sediment and volcanic data were checked carefully. Several sediment age-depth models and volcanic data ages were updated (e.g., when the original chronology was based on now outdated radiocarbon calibration or oxygen isotope reference curves). Paleomagnetically dated and regionally inconsistent sediment records were rejected. Closely adjacent records were then stacked, in a few cases including an alignment of the excursion signal to a master record for the region. The model thus incorporates directional and RPI records from 10 and 12 locations, respectively, comprising 18 declination, 20 inclination, and 35 RPI individual records (Figure 5q).  investigated the robustness of model features by comparing this preferred solution with a large number of similar models based on variations of the LSMOD.1 data set, in particular including end members with no and maximum alignment of the excursion signal. It was found that several general excursion characteristics are robust in all these models (these will be discussed in section 7), while details, in particular the power in individual nondipole SH degrees, should be interpreted with caution.

TAF Models
For even longer times, continuous SH models have not yet been attempted, but a number of SH models of the TAF do exist (Figure 6). Schneider and Kent (1990) and Gubbins and Kelly (1993) each give TAF models for the Brunhes and the Matuyama chrons (past 2.5 Myr), the first using only sediment inclinations for purely zonal models, and the latter using lava flow and sediment data in regularized inversions. Zonal models for the normal and reverse TAF of the past 5 Myr are given by Merrill and McElhinny (1977), McElhinny et al. (1996). Kelly and Gubbins (1997) were the first to include lava flow intensities in addition to directions and to invert the individual measurements rather than temporal (and sometime spatial) averages. Their regularized full SH models for the normal and reverse 5 Myr TAF also includes sediment inclinations. Johnson and Constable (1995) tested whether zonal models could explain their set of 2187 lava flow records from 104 locations (1528 normal and 659 reversed) and concluded that nonzonal structures as found in their regularized full SH models LN1 and LR1 are required. However, they also found from jackknife estimates for different subsets of data that a few, probably unreliable, site data contribute significant structure to the models. Carlut and Courtillot (1998) used the same data and modeling method as Johnson and Constable (1995) to test the robustness of TAF models and investigate effects of data errors and site distribution. Due to limitations from these two factors, the only robust terms identified in TAF models are the axial dipole and quadrupole. The 5 Myr TAF models MF-1 (reversed polarity), MF-2 (normal polarity) and MF-3 (Brunhes chron only) (Shao et al., 1999) make use of the Johnson and Constable (1995) lava flow data but are methodologically different as VGP and equatorial virtual pole (offset 90 • from VGP on the site meridian) distributions are fitted instead of mean vectors. The method has not (yet) been applied to an updated data set.
Models LSN1 and LSR1 (Johnson & Constable, 1997) might be considered as updates of LN1 and LR1, including the sediment inclination compilation by Schneider and Kent (1990), although the overall quality of the sediment data set is poor. The findings by Johnson and Constable (1997), that the lava data require zonal structure while sediment inclinations do not and that models LSN1 and LSR1 indicate insignificant differences between normal and reversed polarity in contrast to their earlier LN1 and LR1, illustrate the strong influence of the data quality and distribution on model results. The much more recent models LN3 and LN3-SC by Cromwell et al. (2018) are direct updates with greatly improved data over the normal polarity LN1 model. They are based on the 5 Myr subset from the new PSV10 0-10 Ma lava flow compilation (Cromwell et al., 2018), with the data for LN3-SC corrected for serial correlation by averaging the directions from all sites from the same stack of lava flows to a single site mean. LN3 and LN3-SC confirm the earlier result that nonzonal structure is required to adequately fit the data.

Global Magnetic Field Characteristics for the Past 100 kyr
The global GGF100k field reconstruction offers the possibility to study a large range of geomagnetic field characteristics over the past 100 kyr. These are discussed in this section and put into perspective with findings from other data syntheses. Many of the geomagnetic field characteristics exhibited by GGF100k and the shorter global reconstructions LSMOD.2 (50-30 ka) and IMOLEe (44-36 ka) that are discussed here are illustrated in comprehensive animations available from the EarthRef digital archive (http://earthref.org/ ERDA/2384/).

DM
Several estimates of DM evolution over the past 100 kyr are shown in Figure 7. The GGF100k DM closely tracks the axial dipole moment (ADM) for the model (Figure 7a), indicating that in general the equatorial dipole contributions are small at this temporal resolution. Also shown is a VADM stack built from all the sediment records that were used in the GGF100k model. The large standard deviation range obtained from bootstrapping is due to both the considerable amount of data that exhibit some variations that are not aligned in age, and the RPI calibration. Many of the previous VADM and scaled RPI stacks (Figure 7b) fall within one standard deviation most of the time, and there is a good general agreement of trends with the GGF100k DM.
Mean, minimum, and maximum values of VADM (or ADM) and their corresponding ages are listed in Table 1 for all the available stacks and models. The values are only given for stacks which provided absolute values, rather than just RPIs, and if records are shorter than 100 ka then the available period is used. Consequently, direct comparisons are not possible among all entries. Furthermore, the absolute values are strongly influenced by the calibration method used.
ADM and consequently also DM were highly variable over the last 100 ka (Figure 7). The earliest intensity low (∼100-90 ka) corresponds to the post-Blake excursion. After some oscillations about higher average DM values, another pronounced minimum (actually appearing as two minima in GGF100k) occurs between ∼66 and 60 ka. This is followed by an increase to a broad maximum that lasts more than 10 ka, with DM values similar to those for historical times and the twentieth century. The highest values appear at 47 ka in GLOPIS-75, and at 49.84 ka in LSMOD.2. GLOPIS-75 actually has two distinct maxima at 53 and 47 ka that  are not found in the other reconstructions. Although the GGF100k model shows more than one peak over the same period, these are not as pronounced as in GLOPIS-75 nor are the lows in between very significant.
The dramatic drop to a very low value at 41 ka, the lowest value seen during the period studied, is the Laschamp excursion, which occurs within the age range 37-42.45 ka in the models and stacks. However, the latest occurrence of the low at 37 ka (29.0 ± 3.2 ZAm 2 ) in Sint-800 (Table 1) is one of the two local minima occurring during the Laschamp excursion, the first one being at a similar time with other models/stacks, at 41 ka with VADM of 31.3 ± 3.0 ZAm 2 . The minimum values for ADM and VADM range from 2.9 ZAm 2 in LSMOD.2 (a dedicated Laschamp excursion model, see also section 7) to 45.1 ZAm 2 in sint-2000, the long-term paleointensity stack spanning the past 2 Ma, which, in general, has higher VADM values and is very smooth over the whole 100 ka interval (Figure 7b).
Most reconstructions suggest that the overall trend since the Laschamp event has been a gradual increase toward another maximum at around 2 ka. About 25 kyr are needed in order for the average post-Laschamp DM to achieve the average pre-Laschamp value (over the period 100-41.6 ka). For 20 kyr before and after the excursion, excluding 0.4 kyr around it (when the modelP i exceeds 0.5), the average pre-and post-Laschamp DM values from GGF100k are 71 ZAm 2 and 64 ZAm 2 , respectively. However, if the whole 100 ka period is considered, pre-Laschamp (100-41.1 ka) and post-Laschamp (40.7-0 ka) values are 66 ZAm 2 and 71 ZAm 2 , respectively. Considering the standard error in the mean of the bootstraps on data VADMs of 1.7 ZAm 2 as a rough uncertainty estimate (though likely too small due to incomplete global coverage), the latter result suggests that the mean DM indeed was slightly lower on average before than after the Laschamp over the time interval considered.
During the period after the excursion, when the DM has not yet recovered to values as high as pre-Laschamp, a few local intensity lows appear (Figure 7), some of which correspond to reported ages of less well constrained excursional events, for example, Mono Lake/Auckland and Hilina Pali (see section 7). Sediment records that constrain the GGF100k in this period show a wide range of ages for these events (see Figure 16), which is reflected in the predicted DM curve. Similar to the 55-45 ka interval, GLOPIS has more pronounced variations than most other models (except for the C2018-Overall stack) between 35 and 30 ka. This might be due to a predominance of Atlantic records in these two stacks (Figures 5c and 5g) and reflect higher than global average PSV in this region (see section 6.3).
A short maximum at about 15 ka is followed by a rapid drop. This feature is much more pronounced in the C2018-Overall stack, with a higher and earlier maximum than seen in any of the other reconstructions. This is one of the few features falling outside the one sigma range of the GGF100k data VADM stack. The average value for the past 100 ka from the five full length models/stacks in Figure 7b is 66.7 ZAm 2 , which is in excellent agreement to the GGF100k ADM of 67.8 ZAm 2 considering the VADM 1.7 ZAm 2 standard error. Interestingly, the two sint-800 and sint-2000 series of stacks predict the lowest and the highest averages over the past 100 ka: 58.6 and 78.6 ZAm 2 , respectively, perhaps reflecting differences in calibration method or data set.
The VADM stack as plotted in Figure 7a provides global averages over 400 year bins, intended to average out the contributions of the non-axial-dipole field in both space and time and thus does not have the same temporal resolution as the DM for GGF100k. The spectrum of GGF100k variations (see supporting information Figure S1) shows several estimates of power spectral density with average frequency resolution ranging from 0.16 to 10 kyr (see also Figure 8 of Panovska, Constable, & Korte, (2018) for the context of the complete spectrum from 0.05/Myr to 10 5 /Myr). Above 2 kyr the spectrum reveals the very strong impact of the spline knot spacing which is 200 years and the temporal damping greatly reduces the power in axial dipole variations. Above 1 kyr we already see the characteristically increasing falloff rate associated with temporal regularization. Between 0.1 and 1.0 kyr (10-1 kyr periods) it is likely that the spectral power is also underestimated because of limited resolution and noise in the data. What is clear is that in addition to longer period trends there is increasing power in the axial dipole variation over this frequency interval and some suggestion of aperiodic power concentrations. Improvements in age control across the various records will be needed to determine whether any characteristic timescale for axial-dipole oscillations can be robustly identified.

Pole Positions and VGPs
Using the 110 inclination sediment records that constrained the GGF100k model, we produced a scalar global VGP latitude stack by averaging individual VGP latitude series (Figure 8). Declination is assumed 0 in the calculation of the VGP. This stack is compared with the dipole latitudes predicted from the GFF100k model (using the dipole coefficients), as well as the global averages of VGPs estimated on an equal area grid (with 1 • spacing at the equator) at 200 year intervals from the model. The latter also contains nondipolar field contributions. The VGP latitude stack from the data is in good agreement with the average VGPs from the model, but the geomagnetic pole latitude from GGF100k is somewhat different. The difference is especially noticeable during transitional periods, for example, over the Laschamp excursion when the field is not dipole dominated. However, from this analysis of VGP latitude over the past 100 ka, it can be concluded that the field was predominantly dipolar. Occasional outliers of 30 • or 40 • reflect times when excursions have been reported, although the excursional deviations are not as large as those predicted from models with higher resolution that are discussed further in section 7. In general, the two curves obtained from the GGF100k model show that the contribution of the nondipole field that is mapped into VGPs is not negligible. As can be seen in the animations of the GGF100k model available at the website (http://earthref.org/ERDA/ 2384/), VGPs are strongly clustered about the geomagnetic dipole axis at 1.5 ka BP when the highest DM is observed over the past 100 ka. VGPs deviate to middle and equatorial latitudes at 41 ka BP and 28.5 ka BP, both times when more complex structure dominates the field. Model results confirm earlier statements, for example, Lund (2018), that high VGP dispersion is associated with low field strength and vice versa.

PSV Index Over the Past 100 ka
The PSV index P i suggested by Panovska and Constable (2017) is used to quantify the geomagnetic field activity over the past 100 ka using both the GGF100k model and by directly stacking results from the individual sediment data (Figure 9). The globally averaged PSV index is generally low over the past 100 ka following the dipole dominated structure seen in DM and VGP, the two components that are used to build the index (previous two subsections). The pronounced peak at 40.6 ka and the smaller side peak at 37.0 ka are associated with the Laschamp excursion, with the latter peak visible only in sediment records from the Southern Hemisphere (Panovska, Constable, & Brown, 2018). Increased activities seen only locally or regionally and/or at different ages are averaged out in the global picture. Apart from the Laschamp excursion, the second largest peaks can be associated with the post-Blake (∼100 ka) and Norwegian-Greenland Sea excursions (∼61 ka). These are the two occasions during the past 100 kyr when the ADM dips below half of the ∼80 ZAm 2 present-day value. A peak around 29 ka appears slightly younger than the Mono Lake excursion and is discussed further in section 7.
The distribution of the P i can also be evaluated separately by hemispheres over the past 100 ka. Figure 10 indicates variations in P i indices from GGF100k averaged every 2 • by latitude in (a) and longitude in (b)  from estimates at 500 year intervals. In general, the indices are higher in the Southern Hemisphere than the Northern Hemisphere and in the Atlantic compared to Pacific Hemisphere, a characteristic already observed for present, historical (Panovska & Constable, 2017) and Holocene geomagnetic field models (Constable et al., 2016). In the regional VADM models produced by Ziegler and Constable (2015), they also found higher average ADM in the Pacific than in the Atlantic, associated with lower field variability over the past 300 ka. These results support the idea that the hemispherical asymmetry is a long-term pattern of the field behavior. Large-scale heterogeneities in conditions at the top or the bottom of the outer core that are considered responsible for the asymmetry observed in direct and historical observations (e.g., Aubert et al., 2013) and also found in numerical dynamo simulations with heterogeneous CMB heat flow (e.g., Terra-Nova et al., 2019), may also play an important role at multimillennial timescales.

From the Holocene to 5 Ma TAF
The first-order approximation of the geomagnetic field as a GAD is widely used, but it is important to characterize second-order contributions and their persistence at various timescales. Thus, we have compared in Figure 11 the TAF determined from the mean Gauss coefficients over the past 10 ka (from the CALS10k.2 model), 100 ka (GGF100k), and 5 Ma normal polarity field (LN3 model) at both the CMB and at Earth's surface. We look at both the radial field component B r and nonaxial dipole (NAD) part of B r at the CMB, and at inclination anomaly (differences between inclination of the total field and the GAD prediction) and declination at the Earth's surface. It is clear from the time-averaged B r at the CMB that the axial dipole component is dominant over all three intervals. Intense normal flux lobes appear in the Northern Hemisphere, with two clear regions in the Holocene model, over North America and Asia (for more details of the TAF in Holocene models see Panovska et al., 2015;Constable et al., 2016). The averaged B r NAD field shows similar structure in all three models: positive anomalies over the west Atlantic Ocean/Gulf of Mexico, north-central Africa, and western equatorial Pacific. Negative anomalies appear in the North and South Atlantic Ocean,

Reviews of Geophysics
10.1029/2019RG000656 Figure 12. Ratios of the time-averaged NAD terms to the axial dipole coefficient (g 0 1 ), up to degree 4. Time-averaged fields (TAF) are estimated from the following models: 10 ka from the CALS10k.2, 100 ka from the GGF100k, and 5 Ma normal polarity field from the LN3 model. while the mildly positive B r NAD anomaly at high southern latitudes, over Antarctica, is present only in the 100 ka and 5 Ma averages. In general, the overall structure of latitudinal/longitudinal variations of the average NAD radial field observed over the past 100 ka is also present in the several million year averages (e.g., Cromwell et al., 2018;Johnson & Constable, 1997).
The TAF inclination anomalies in all three models have larger negative than positive values, the largest being observed in the GGF100k model. The structure looks similar, negative inclination anomaly in the equatorial region, with peaks over Africa and/or Indonesia region (western Pacific), manifested in the longitudinal variations of the magnetic equator. Positive inclination anomalies span all latitudes in the eastern Pacific in GGF100k, but only northern parts for CALS10k.2 and mostly southern latitudes in LN3. Time-averaged declination shows large regions of positive and negative anomalies and similar patterns in the 10 and 100 ka averages in contrast to the 5 Ma. This probably results from the significantly larger data sets that constrain the CALS10k.2 and GGF100k models.
The SHA decompositions of the field for the three models allow us to analyze the terms responsible for the TAF structure explained hitherto. NAD coefficients up to degree 4 as a percentage of g 0 1 are presented in Figure 12. g 1 1 and h 1 1 have notably higher percentages in the 100 ka TAF compared to the 10 ka and 5 Ma averages. The axial quadrupole term (g 0 2 ) contributes significantly on all timescales: 4.8%, 4.2%, and 3.0% of the axial dipole in 10 ka, 100 ka, and 5 Ma TAF, respectively. These values lie in the range previously found by other TAF studies (2-5% of g 0 1 ) that recognize the axial quadrupole term as one of the largest contributions to persistent deviation from the GAD (Johnson & McFadden, 2015). The axial octupole term (g 0 3 ) of 2.5% and 1.3% of g 0 1 in 100 ka and 5 Ma also make a substantial contribution, although g 0 3 is negligible in the 10 ka TAF. Other nonzonal dipole, quadrupole and octupole terms (up to degree and order 3) that have nonzero values are also responsible for persistent deviations of the geomagnetic field from the GAD over the periods studied, but specific persistent nonzonal patterns are difficult to identify. Some of these terms, for example g 2 2 , g 2 3 , h 3 3 , show clear differences in sign and magnitude for the different timescales. Nevertheless, all three models suggest nonzonal (longitudinally varying) structures explain the variations in the TAF (as seen in Gubbins & Kelly, 1993;Johnson & Constable, 1995;Johnson & Constable, 1997;Cromwell et al., 2018).

Field Morphology
An interesting observation from the GGF100k model is the similarity of the present-day Earth's surface and CMB field to the field observed at several points over the past 100 ka (Figure 13). The field morphologies seen at 85.5, 69.2, 49.2, 45.9 (before the Laschamp excursion), and 33.2 ka at Earth's surface show low field intensity in areas over the South Atlantic Ocean or South America.
During these periods, the location of the reverse flux patch on the CMB is similar to the current position of the reverse patch which causes the South Atlantic Anomaly (SAA) (Figure 13b). Brown et al. (2018) noted this similarity in a SH model spanning 50 to 30 ka BP and concluded that this seems to be a recurring feature, but that it does not indicate the beginning of an excursion or field reversal. The periods identified in the GGF100k when the field resembles the SAA (Figure 13) in most cases have ADM values in the range of 63 ZAm 2 (with the exception of 80 ZAm 2 at 49.2 ka) and are related to weakly pronounced ADM (or g 0 1 ) minima. The field starts to increase after or shortly after these times (150 years for the 69.2 ka feature). Terra-Nova et al. (2017) showed that the SAA clearly results from a combination of reversed and normal

10.1029/2019RG000656
flux regions below the South Atlantic, and the axial dipolarity of the field influences the position of the minimum field intensity at the Earth's surface. Unusual CMB composition defined by an LLSVP (large low shear wave velocity province) beneath South Africa has been suggested as the reason for the SAA and the LLSVP's longevity could play a role in recurrent weak field or persistent triggering of transitional events in this region (Bloxham & Gubbins, 1987;Gubbins, 1987;Tarduno et al., 2015;Terra-Nova et al., 2019).
GGF100k supports an assessment that SAA-like features do not necessarily directly or always lead to transitional events. Times when transitional field (P i >0.5) patches are found at the Earth's surface are usually characterized by more than one intensity minimum at the Earth's surface and more than one reversed flux patch at the CMB. Figure 14 presents the field morphology at Earth's surface and the CMB for three periods in the early phases of proposed excursions.
In all cases, the field presents wider areas of low field intensities than at present day and reverse flux entering both Northern and Southern Hemispheres (see also the animations available from the website (http:// earthref.org/ERDA/2384/)).
In Figure 15 we show the variation in power in the dipole compared to the sum of the power in all ND degrees up to SH degree and order 5, approximately the effective spatial resolution of these models, from GGF100k together with that for the 50-30 ka model LSMOD.2 and the Holocene CALS10k.2. Some clear disagreements in the comparison of GGF100k with these two higher-resolution models shows that the ND contribution, and in particular power in individual SH degrees might not always be robustly resolved. Differing data and age constraints play a role here. Note, however, that with incomplete global data coverage the available data at the Earth's surface might in general be fit similarly well by different distributions of power among SH coefficients (see, e.g., . Surface field characteristics therefore will be more robustly resolved than CMB structure and individual SH degree power, and these should be interpreted with care.

Excursions
Three excursions are listed as confirmed for the past 100 ka in a recent comprehensive review by Laj and Channell (2015). The best documented is clearly the Laschamp (∼41 ka), which has been found in records distributed around the world. Note that it was recently suggested to revise the spelling to Laschamps, because apparently the original "s" at the end of the name has been forgotten by cartographers since the eighteenth century (Kornprobst & Lénat, 2019). However, although the adajenct village indeed is spelled "Laschamps" on many modern maps, the volcano after which the excursion is named is still mostly labeled "Puy de Laschamp," and we here maintain the traditional spelling "Laschamp" for the time being. The Mono Lake (∼33 ka) and the Norwegian-Greenland Sea (∼60 ka) excursions are less well documented globally. Laj et al. (2014) and Singer et al. (2014) argue for renaming the Mono Lake to the Auckland excursion because of recent dating results indicating that the excursion originally recorded at the Mono Lake type locality might be the Laschamp (see also Marcaida et al., 2019). We will use the term Mono Lake/Auckland from here on. Additional, less established events have been reported. Nowaczyk et al. (1994), Singer et al. (2014) and Ahn et al. (2018) list several postulated events younger than 20 ka, which, given the age uncertainties, might be manifestations of the Hilina Pali excursion now mostly dated at 17 ka (Singer, 2014). Individually reported excursion signatures between ∼106 and 91 ka and probably belonging to the same event have been termed Fram Strait (Nowaczyk et al., 1994), Post-Blake  or Skálamaelifell excursion (see Jicha et al., 2011). The term most widely used in present literature seems to be Post-Blake, which we adopt in the following. Figure 16 gives an overview of the locations where excursional data have been found in the respective time intervals. It is difficult to establish the global magnetic field behavior during excursions from individual records, so in the following we take a global view based on data syntheses. The syntheses give a regional or global interpretation of the signals found in individual records. To some degree they can be able to take physical consistency of information from different regions into account and disregard inconsistent signals. However, they also might be affected by erroneous or missing information, such as incorrect chronologies or unresolved fast variations. In the following we describe conceivable magnetic field behavior during excursions. We acknowledge that other interpretations, in particular if assuming that many excursions are often not recorded due to brevity, are possible.
With the possible exception of Hilina Pali all these excursions are indicated, more or less clearly, by minima in global intensity stacks (see Figure 7), the ADM of GGF100k (Post-Blake, Norwegian-Greenland Sea, and Laschamp, see Figure 7) and/or the VGP latitude stack in Figure 8 (for Post-Blake, Norwegian-Greenland Sea, Laschamp, and Mono Lake/Auckland). For the rest of this section we use the P i threshold of 0.5 whenever estimating beginning end or duration of an excursion. Although the globalP i from GGF100k only exceeds the 0.5 threshold for excursions and reversals suggested by Panovska and Constable (2017) during the Laschamp excursion, Figure 9 shows mildly elevated values at ∼28 ka (close to Mono Lake/Auckland), compared them, analyzed the robustness of their excursion features, and updated the 50-30 ka model to LSMOD.2. A general conclusion with relevance for the interpretation of the other excursions in GGF100k is that, as for the Holocene, the long-term model provides a strongly smoothed picture of the events and several characteristics appear subdued. The IMOLEe model appears strongly affected by the sparse data distribution ( Figure 5p) so that findings regarding dipolar and nondipolar variations and field morphology during the Laschamp excursion from that model seem outdated. Animations of all three models are provided at the website (http://earthref.org/ERDA/2384/), where there are many more details than are discussed below.
LSMOD.2 and GGF100k disagree about the Mono Lake/Auckland excursion age: LSMOD.2 recovers it as a series of two regionally manifest events at around 34 and 31 ka. GGF100k has no indications for an excursion at these ages, but does show regional midlatitude VGPs and mildly elevatedP i at about 28 ka (an age not covered by LSMOD.2 or any other global model at present). A pure age discrepancy might be due to age updates and regional age alignments in LSMOD.2, but general characteristics of these events also differ, as will be seen in the following. From the comparison of the Laschamp excursion manifest in both models it seems conceivable that the two events seen in LSMOD.2 are suppressed by smoothing in GGF100k and a third (possibly stronger) event occurred around 28 ka. The existence of an excursion around that time is supported by a recent record from the NE Atlantic but not yet included in the compilation. Reversed VGP Figure 15. Dipole (D), nondipole (ND), power at the CMB from models GGF100k (100-0 ka; black, blue), LSMOD.2 and CALS10k.2 (50-30 ka and 10-0 ka, respectively; gray, cyan). The ND power is calculated from SH degrees 2 to 5, higher degrees contribute little to these models' ND power. . Sediment records are in blue/light blue and volcanic data are in red/light red colors. Note that on the histogram there are a few events that are not in the time ranges considered for the maps. If a time interval is reported for an event instead of a fixed age then the average of that interval is taken for the histogram plot. Excursions on the maps are grouped according to their age. Sometimes, especially in older publications, an excursion can be named differently from the names given on these maps. If excursional features and excursions are reported at the same location, the excursion symbol is plotted on top and the light color symbols are not visible. Supporting information Table S1 lists all locations, references, ages, and names of excursions used to produce these plots.
latitudes were found in there around 26.5 ka . We keep the name Mono Lake/Auckland for the two events seen in LSMOD.2 and label the GGF100k event "GGF-28k" in the following discussion of excursion characteristics.
Other data can also be brought to bear on the question of how to reconstruct the past geomagnetic field. Specifically, cosmogenic nuclides, such as 14 C, 10 Be, and 36 Cl, are produced in the Earth's upper atmosphere by nuclear interactions between energetic cosmic ray particles with target elements. A nonlinear inverse relationship has been found between paleomagnetic dipole strength and production of cosmogenic isotopes (e.g., Masarik & Beer, 1999). Isotopes from any geographic location are in general considered to give an indication of global magnetic dipole strength under the assumptions that the cutoff rigidity (defined as the momentum per unit charge of a particle that shows its ability to penetrate the geomagnetic field) is not strongly influenced by higher-degree geomagnetic terms, that there is rapid global mixing in the atmosphere, and adequate corrections can be made for any environmental disturbances such as weathering effects on 10 Be estimated production rates (e.g., Carcaillet et al., 2004;Frank et al., 1997;Robinson et al., 1995). In  Figure 17 the (axial) dipole moment variations from the GGF100k model are compared to corresponding variations derived from two sets of 10 Be records. The 50-20 ka 10 Be-derived VDM by Ménabréaz et al. (2012) is based on an authigenic 10 Be/ 9 Be marine stacked records from equatorial Pacific and northeast Atlantic, and VADM from Simon et al. (2016) uses two cores from the equatorial Pacific covering the past 100 ka. The calibration to VDM values had been carried out with the help of absolute paleointensities from the GEO-MAGIA50.v3 database in both cases. Ménabréaz et al. (2011) linked the minimum in VDM from the highest production rate of cosmogenic nuclides measured in the Portuguese Margin sediment core to the Laschamp excursion. Indeed, there is a clear agreement to the Laschamp DM minimum in GGF100k in both reconstructions from the isotope records, and minima related to other excursion are also seen. The longer-term record in Figure 17b shows very similar minima related to the Post-Blake and Norwegian-Greenland Sea excursions as GGF100k, supporting the inherent assumption that the isotope record is primarily responsive to global dipole changes. Additional minima not seen in GGF100k Figure  17b). It seems quite conceivable that the 34 ka and a weaker 31 ka minimum in the Ménabréaz et al. (2011) record are expressions of the contemporaneous events found in LSMOD.2, but not resolved in GGF100k and the Simon et al. (2016) record. The stronger ∼26a and ∼28 ka minima in the two isotope-based records, respectively, might represent the GGF-28k event within the different age uncertainties of the records, and the minimum at ∼19 ka might be linked to the Hilina Pali excursion.

Excursion Duration and Ages
In Figure 18, the maximum P i values found during the suggested excursion intervals (see left panels of the figure) are mapped together with duration (middle) and earliest occurrence (right panels) of P i ≥0.5. The Laschamp is clearly the strongest excursion, and is seen globally in model LSMOD.2, with lower maximum P i over the Pacific and Indian Ocean and highest values over the Atlantic and in parts of the Southern Hemisphere. The comparison between Figures 18a and 18b illustrates the influence of stronger regularization in GGF100k, where P i exceeds 0.5 only in the regions with the strongest PSV, and not globally. The regional duration reaches up to 3.5 kyr in LSMOD.2, but only 2.5 kyr in GGF100k. It is likely that GGF100k also underestimates the regional duration and spatial extent of the other excursions. In particular, P i never exceeds 0.5 anywhere on Earth during the Norwegian-Greenland Sea and Hilina Pali excursions in that model. Many data and model predictions indicate that fast and strong directional changes occur during the intensity minimum, and excursion durations appear shorter based purely on directional changes than on the intensity low and the strong directional changes. For instance, Laj et al. (2014) report durations of 640 years from directional and 1,500 years from intensity results of the Laschamp from the French Chaîne de Puys, and Nowaczyk et al. (2012) found 440 years for the excursion directional swing in Black Sea sediments, while the intensity low at the same location (Nowaczyk et al., 2013) has a width of about 2 kyr. Rapid directional field changes will occur naturally if the dipole field is weak but the ND secular variation continues unchanged (e.g., Brown & Korte, 2016;Constable, 1990).
For the Mono Lake/Auckland excursion, LSMOD.2 predicts strongest P i over the southwestern Pacific and Indian Ocean, and the northeast Atlantic. The plots of start age and duration indicate that the excursional signals from these two regions are separated in time by ∼2 kyr, so that they might be seen as two different, regional events of short, a few centuries to ∼1 kyr, duration. Some high-resolution cosmogenic isotope records also show more than one peak that might be related to the Mono Lake/Auckland excursion, for example, the 36 Cl record presented by Wagner et al. (2000) has maxima at 31.5 and 34 ka. Interestingly, none of the models predicts P i ≥0.5 around the type locality in North America. This might further support the hypothesis that the Mono Lake type locality results belong to the Laschamp excursion, but it also has to be kept in mind that the smoothed models might fail to fit the full amplitudes of regional field variations. However, since P i applied locally to sediment records show considerably lower values and at fewer locations for the Mono Lake compared to the Laschamp excursion (Panovska, Constable, & Brown, 2018), weaker P i signals and regional variability are expected in the models. GGF100k predicts an excursional field over south America corresponding to GGF-28k, also evident as a dipole low in the 10 Be/ 9 Be records in Figure 17. The differences noted above between the two 10 Be/ 9 Be records and also from GGF100k might be due to either differences in resolution or age discrepancies.
Several patches of regional transitional field behavior according to P i ≥0.5 go along with the smaller maxima in the global PSV index from Figure 9 related to the Norwegian/Greenland Sea and Post-Blake excursions. The peak activity appears in different regions at different times, that is, over the West Indian Figure 18. Maximum P i , duration, and starting age of excursion with P i ≥ 0.5 in global models within the past 100 ka. Results from LSMOD.2 (a) and GGF100k (b) for the Laschamp excursion, from LSMOD.2 for the Mono Lake/Auckland (c) and GGF100k for the GGF-28k excursion (d), and from GGF100k for the Post-Blake excursion. White contour line in maximum P i plots is 0.5. That value is not exceeded in GGF100k for the Norwegian-Greenland Sea and Hilina Pali excursions.

Reviews of Geophysics
10.1029/2019RG000656 Figure 19. Minimum VADM and its occurrence age and minimum VGP latitude and its occurrence age from left to right for the Laschamp from models LSMOD.2 (a) and GGF100k (b), the Mono Lake/Auckland from LSMOD.2 (c), and the GGF-28k (d), the Norwegian-Greenland Sea (e), and Post-Blake (f) excursions from GGF100k. Black and white patches indicate older (black) or younger (white) ages outside the shown time interval. White contour lines in minimum VGP plots are 45 • .  Figure 19) but P i <0.5). Areas of stronger PSV broadly agree with the locations, for which excursions or excursional field features have been reported (compare Figure 16), while P i and local field variation amplitudes and regional excursion durations are likely to be underestimated by the smoothed model.
Although P i exceeds 0.5 (in LSMOD.2) all over the globe in the course of the Laschamp, we find no epoch where this occurs simultaneously everywhere, and similarly there is no epoch where transitional or reversed directions are seen at the same time all over the globe. As seen in Figure 18 and from Table 2 the regional duration of the excursions varies, from a few centuries to ∼3.5 kyr. This is in good agreement with various excursion durations reported from individual records (see Laj & Channell, 2015;Roberts, 2008). Global duration according toP i ≥0.5 can only be determined for the Laschamp, and is 1.8 kyr according to LSMOD.2, somewhat shorter than an independent estimate of 2.5 kyr from full width at half maximum of the peak in low-pass filtered 36 Cl production from the Greenland Ice Core Project ice core (Wagner et al., 2000).
The ages of the excursion midpoints according to maxima in theP i timeseries of GGF100k are 98.8 ka, 61.6 ka, and 17.3 or 14.5 for the Post-Blake, Norwegian-Greenland Sea and Hilina Pali excursions. For the Laschamp, both GGF100k and LSMOD.2 give 41.1 ka, and for the Mono Lake/AucklandP i occur at 34.7 and 31.6 ka (LSMOD.2) and for GGF-28k at 28.5 ka (GGF100k). These ages agree well with the ages obtained from compilations of radioisotope dating of volcanics from different locations for each excursion by Singer et al. (2014). He obtained ages of 100 ka for the Post-Blake and 32 ka and 17 ka for Mono Lake/Auckland and Hilina Pali, considering each of these as only one event. His age determination of 40.7 ± 0.9 ka for the Laschamp, the Laschamp excursion midpoint at 41 ka in GLOPIS-75, and an independent result from the 10 Be flux peak in the NorthGRIP ice core (Singer et al., 2009) are all in excellent agreement with the results from the global field models.

VGP Paths and Minimum Field Intensity
Geomagnetic excursions have mostly been defined by VGP deviation from the GAD, for example an arbitrary 45 • threshold (e.g., McElhinny & Merrill, 1975;Wilson et al., 1972). It was often noted that some records exhibit reversed VGPs in high southern latitudes (e.g., Channell, 1999;Laj et al., 2006;Nowaczyk et al., 2012) and others reach only low to middle latitudes (e.g., Horng et al., 2003;Meynadier et al., 1992;Tric et al., 1992), which may be attributed to sediment recording fidelity. Low sediment accumulation rate and inefficient recording have been suggested as a reason for the latter (Channell & Guyodo, 2004;Roberts & Winklhofer, 2004). However, equally strong deviations from axial dipole direction should not be expected in all regions if the field is not dipole dominated during excursions Brown & Korte, 2016).
Very few predicted VGPs from GGF100k reach high southern latitudes since the model has smoothed the fully reversed directions present in some individual sediment records and not fully captured their complexity;compare Figures 19a and 19b. This figure includes minimum VADM and minimum VGP latitude with their ages of occurrence in the respective time intervals. However, the VGP predictions of LSMOD.2 ( Figure 19a) do not reach southern latitudes all over the globe either, suggesting that strong regional variations in minimum VGP latitudes can be real.
The regional distribution of minimum VADM and VGP latitudes individually broadly agrees with that of maximum P i as expected by the P i definition (cf. Figure 18). The models accommodate regional differences in minimum VADM and VGP latitude of up to a few kiloyears. Although we cannot rule out that insufficient age control might distort excursion signals in a model, this result indicates that regional age differences in

10.1029/2019RG000656
paleomagnetic excursion signatures can be real and care should be taken when using them as stratigraphic markers. The sparsity of black and white patches in Figure 19, which indicate that smaller VADM or VGP latitude values are found shortly before or after the displayed time interval in Figure 18, indicates that minima are seen in time series of these quantities even if they do not qualify as transitional by standard definitions.
Field intensity might have produced VADM values below 30 ZAm 2 everywhere during the Laschamp excursion, with minimum values < 1 ZAm 2 (Figure 19a). This is slightly lower than very low intensities found in lavas from France (Laj et al., 2014) and New Zealand (Mochizuki et al., 2006). Although we have to consider the more limited resolution of GGF100k, it seems clear that none of the other excursions within the past 100 ka reached similarly low values over a comparably large portion of the Earth.
A clockwise loop of VGPs moving south over Asia and the west Pacific and then recovering northward through Africa and west Europe has been found in a compilation of high-resolution sediment paleomagnetic records of the Laschamp excursion by Laj et al. (2006). However, Laj and Channell (2015) note that VGPs of other records of the same excursion make the transition at different longitudes. A general caveat when interpreting VGPs is the need for orientation of relative declinations, that has an influence in particular on VGP longitude.
The investigation by  on Laschamp VGP paths from a range of global models supports the results by Laj et al. (2006) that a higher than average number of paths loop down over Pacific and up over African/European longitudes. This is shown in Figure 20a, together with temporal distribution of VGPs below 45 • , where the maximum is found close to the excursion midpoint. The smoother GGF100k does not fully reproduce the longitudinal preference ( Figure 20b) but does give preferred longitudes roughly around 90 • W to 120 • W for all excursions (Figures 20d-20f), although they are observed from very different locations (Figures 21d-21f). LSMOD.2 has a similar longitudinal VGP path distribution for the stronger 31.5 ka part of its bimodal Mono Lake/Auckland distribution (Figure 20c) as for the Laschamp. Despite a clear preference for certain VGP longitudes not all VGP paths follow the same pattern, and a simple excursion geometry of a rotating dipole with reduced nondipole contributions (Laj et al., 2006) is not supported by the global field models (see also section 7.3).

Field Morphology During Excursions
Maps of surface field intensity, VGP latitude distribution, and B r at the CMB from GGF100k and LSMOD.2 are shown in Figure 21 for the excursion midpoints according to maximum P i .  concluded from their robustness analysis that the complex patterns of nondipole dominated CMB flux evolution predicted by SH models during excursions depend notably on data selection, treatment for regional consistency and assumptions about temporal alignment of the observed intensity minima. Therefore, it is not surprising that GGF100k and LSMOD.2 show quite different B r patterns, at concurrent epochs during the Laschamp as well as at the respective excursion midpoints as shown in Figures 21a and 21b. In particular, LSMOD.2 has significant amounts of reverse flux entering into the tangent cylinder in both hemispheres during the Laschamp. At first sight, the field morphology looks significantly different in the middle of the different excursions, but reverse flux seems to appear preferentially around south Africa, western South America, and eastern North America. This is probably linked to the preferred VGP paths, a feature that might not be expected from the various distributions of VGP latitudes observed during the excursion midpoints. Surface intensities also have different patterns for the different excursions. LSMOD.2 suggests that the field was weakest in the Atlantic, south America, and Indian Ocean in the middle of the Laschamp, but the lower-resolution GGF100k gives a different picture. The 28 ka event ( Figure 21) appears extremely localized around central and south America with unusually strong dipole field structure and high northern latitude intensities compared to the other events.

Excursion Mechanism
The model observations resemble the excursion scenario of strong axial dipole variation with persistent normal secular variation in nondipole contributions Brown & Korte, 2016). In particular, transitional and reversed directions do not occur all over the globe at the same time. This would require the axial dipole contribution to at least reverse sign, together with a (slight) increase of DM during the central phase of the excursion, which is not predicted by any of the models. Both the GGF100k and LSMOD.2 models suggest that excursions might be produced by stronger than usual axial dipole strength variations, perhaps including some transfer of energy from axial dipole to nondipole, but

Reviews of Geophysics
10.1029/2019RG000656 without significant changes in NAD secular variation or contributions from equatorial dipole contributions . Figure 15 confirms the hypothesis by Brown et al. (2018),  that excursions occur when dipole power drops to nondipole levels at the CMB in field models with sufficient resolution to around SH degree 5: D and ND have similar values around 100-93 and 60 ka, ND clearly dominates during the Laschamp, the strongest of the excursions, the levels remain similar after the Laschamp from 36 to about 28 ka, and again about 21-18 ka. There is a clear difference between the Laschamp, where dipole power drops far below ND power for ∼5 kyr, with the axial contribution getting close to 0 and the other excursions. For these, dipole and nondipole fields remain at similar levels for several millennia, likely causing series of two or more regionally confined excursions at Earth's surface during all these times as found in LSMOD.2 for the Mono Lake/Auckland at 34 and 31 ka. Care has to be taken when using such excursions as stratigraphic markers, as similar geomagnetic excursion signatures from different regions might be offset by a kiloyear or more.  suggested three phases of geodynamo behavior: (1) a stable, strongly dipole dominated state, (2) strong axial dipole variations causing globally observed excursions such as the Laschamp, and (3) a state where D and ND are of similar strength and regional excursions occur at Earth's surface. It is unclear if the second state is the same as occurs during reversals. All three states might be part of a continuum of secular variation , with the apparent differences caused only by the amplitude and possibly also speed of axial dipole variation. In that case, there might be no substantial difference between excursions and reversals, but it might be just a matter of chance whether the dipole recovers to normal or reverse polarity once it has dropped to zero or reversed sign. Similar field behavior with stable phases of high DM and phases of low DM with excursions and reversals has been found in a long numerical dipole simulation by Wicht and Meduri (2016). There is one exception in GGF100k from the scenario of excursions being generated by strong axial dipole variations in GGF100k, and that is the GGF-28k event. Here, ND power rises strongly, clearly surpassing the moderately varying dipole power for a short interval ( Figure 15). It remains to be seen whether similar field behavior will be found for other epochs, or whether some problematic or inconsistent data within this time interval have caused an unrealistic description of spatial spectral power distribution for this feature in GGF100k. The finding of clear V(A)DM minima obtained from 10 Be production rates around this time might indicate that more dipole variability than predicted by GGF100k might be involved in the generation of this event.
In general, the total nondipole power did not drop during any of the excursions. On the contrary, it increased during the Laschamp, a result that is in good agreement across both models. Transfer of energy from dipole to higher degree structure and back (Amit et al., 2018;Amit & Olson, 2010;Huguet & Amit, 2012;Williams et al., 1988;Williams & Fuller, 1981) might play a role. The increase in ND power across the Laschamp in GGF100k and LSMOD.2 contrasts with the earlier conclusion by Leonhardt et al. (2009) from the IMOLEe model, which predicted a ND decrease along with the dipole drop.

Linking Earth Properties to Dynamo Studies
The growing number of numerical geodynamo simulations now available has led to efforts to provide statistical criteria that can be considered Earth-like (Christensen et al., 2010;Davies & Constable, 2014), and can be readily computed for any numerical simulation. To date these Earth-like criteria have been reliably estimated from modern field models, which lack a sufficient time sample but have high spatial resolution, and from early Holocene models with longer time spans, but poor temporal and spatial resolution. Statistical models describing PSV over the past few million years have received limited attention. The global GGF100k model allows a comparison of these time-varying geomagnetic field properties over an extended timescale of 100 kyr.
We use four quantities proposed by Christensen et al. (2010) and determine them from both GGF100k and two state of the art Holocene models to provide more reliable long-term constraints for testing geodynamo simulations. Each one is evaluated at the CMB, a straightforward calculation for a numerical simulation or a SH model, but requiring some care when comparisons are made between highly detailed and lower-resolution paleofield models. The properties used are (1) AD/NAD, the ratio of power in the axial dipole to that in the rest of the field; (2) O/E, the ratio of the power in equatorially antisymmetric (l − m odd) to equatorially symmetric (l-m even) nondipole components (SHs of degree 2 and higher); (3) Z/NZ, the ratio  Table 3. All models are truncated to degree 5. Mean and standard deviations values are estimated over the model validation periods. Gray dashed lines are 1.0 for the AD/NAD, 0.78 and 0.143 obtained for the purely random equipartioned nondipole field truncated at degree 5 for the O/E and Z/NZ, respectively, and 0.8 for the FCF of a pure dipole field.
of power in zonal components to nonzonal components of the nondipole field; and (4) the degree of spatial concentration of magnetic flux (flux concentration factor, FCF). The first three quantities are derived from the spatial power spectrum at the CMB, see for example, Equations 1 and 2 in Christensen et al. (2010) and Equation 6 in Davies and Constable (2014). FCF = [⟨B 4 r ⟩ − ⟨B 2 r ⟩ 2 ]∕⟨B 2 r ⟩ 2 , where B r is the radial component of the magnetic field and angle brackets stands for the mean value over the spherical surface.
We have estimated these quantities at full model resolution (l = 10) and also after truncation at l = 5. Figures 23 and 24 demonstrate that this makes hardly any difference in the 100 ka model and only a small difference in Holocene models, which confirms the models have little power in terms with l > 5. Christensen et al. (2010) based their analysis on models truncated at degree 8. The influence of the truncation on the individual quantities is summarized in table 3, where we compare them with the IGRF12 up to degree 13 and truncated to degrees 8 and 5, together with the values averaged over their respective time intervals for GGF100k and Holocene models CALS10k.2 and HFM.OL1.A1 (all for maximum l = 5). The model averages together with standard deviations from the variability over their time spans are shown in Figure 22. Supporting information Figures S2 and S3 illustrate how the quantities vary in IGRF12 with time and truncation level. Note that direct comparison to the CAH10 nominal values (Christensen et al., 2010) is not possible due to the different truncation degree.

Dipolarity
The AD/NAD ratio obviously depends strongly on the truncation degree of the model, and whether excursions are present (Table 3). D/ND, the ratio of power in the dipole to the nondipole part of the field has similar variations to AD/NAD, and therefore, it is not considered separately here.
Most of the time AD/NAD for GGF100k clearly exceeds the present-day (IGRF12) ratio truncated at the same degree. The main exceptions arise during the previously discussed excursions when NAD power occasionally exceeds axial dipole power even within these low degrees. This does not happen during the Holocene according to either CALS10k.2 or HFM.OL1.A1: Note that the two models show some overall differences related to variations in their spatial and temporal resolutions (Figures 24a and 24b). Particularly, large values of AD/NAD are seen in GGF100k at around 2 and 10 kyr ago, at about 52-54 ka and just before the Laschamp excursion, and highly variable values occur between 70 and 92 ka, that is, in the time interval between the Norwegian-Greenland Sea and Post-Blake excursions ( Figure 23). The Holocene models additionally suggest unusually high dipolarity around 5 ka. Interestingly, as noted earlier the low AD/NAD ratio associated with the GGF-28k event seems to reflect an increase in NAD rather than a decrease in AD.

Equatorial Symmetry
O/E depends much less on truncation degree (Table 3), and in fact, the 100 ka mean is close to the average IGRF12 value over the time interval 1900-2015. The Holocene values are slightly higher on average, but with a large overlap in standard deviations (Figure 22) This ratio has been used to assess the stability of simulated geodynamo fields (Coe & Glatzmaier, 2006) as well as in a parameter that represents the variation of VGP dispersion with latitude (Tauxe & Kent, 2004). Christensen et al. (2010) estimated a value of 0.833 for a purely random equipartioned nondipole field extending to l = 8 (i.e., a white spectrum at the CMB). The equivalent value for truncation at l = 5 is 0.78, while l = 13 yields 0.882. The average IGRF values for truncation to l = 13 and l = 5 are the same with decreasing variability (standard deviation). Significant variability can be present in the various contributions at any instant in time (see supporting information Variations of field dipolarity given by axial dipole to non-axial-dipole ratio (AD/NAD), equatorial symmetry given by the ratio of odd to even (O/E) spherical harmonics for the nondipole field, zonality by zonal to nonzonal power ratio (Z/NZ) of the nondipole field, and the flux concentration factor (FCF) at the CMB, all from the GGF100k model. Red dashed lines represent the model with all harmonics up to degree and order 10, and blue solid lines are for the model truncated at degree 5. Mean values are plotted with dotted lines in the same colors, respectively, for the total and truncated model. Yellow and gray areas indicate ratios above and below 1 for the AD/NAD, 0.78 for O/E and 0.143 for Z/NZ (values for purely random equipartioned nondipole field truncated at degree 5) and 0.8 for FCF (a pure dipole field). Figure S2). The mean O/E ratio for GGF100k is 0.85, and there are only a few intervals where it shows strongly dominant equatorially antisymmetric contributions in the ND field ( Figure 23).
The mean O/E ratios of CALS10k.2 and HFM.OL1.A1 for the past 10 ka are 1.15 and 1.21, respectively, which suggests a stronger antisymmetry and a more stable paleomagnetic field. Note that the Holocene mean calculated for GGF100k (0.60) is significantly lower than for the Holocene models because the temporal variability is more strongly damped in GGF100k. The Laschamp excursion does not have a characteristic signature in the O/E ratio. With a value of 0.80 the excursion time interval is comparable to the 100 ka average. The unusual maximum at 18.5 ka BP arises from low power in the symmetric components at a time when the antisymmmetric components are high. Christensen et al. (2010) note that the ratio Z/NZ would be 0.10 for a purely random equipartioned ND field extending to degree and order 8. This must be modified to Z/NZ=0.143 for truncation at l = 5, but all the GGF100k and Holocene average values are still quite a bit larger than that, indicating a greater dominance of zonal components in the nondipole field than expected for a white spectrum at the CMB. Of course the sum of the nondipole power in the nonzonal components at any given time is almost always larger than the power of the zonal components over the past 100 ka (Figures 23 and 24) as there are more nonzonal harmonics. The mean Z/NZ ratio of GGF100k is 0.28 and agrees well with the CALS10k.2 and HFM.OL1.A1 means (0.23 and 0.33, respectively) and the present-day IGRF12 value (Table 3). This quantity also does not depend much on truncation degree in the various paleomagnetic models. However, GGF100k does show a difference between 0 and 20 ka and longer-term behavior in Z/NZ ratio, with fewer fluctuations and lower average in the more recent time interval. This is in general confirmed by the Holocene models ( Figure 24), with the exception of HFM.OL1.A1 showing variations similar to earlier times in GGF100k between 7 and 9 ka BP.

Flux Concentration
FCF describes radial field concentration, and spatial model resolution clearly has to be considered in its interpretation (Table 3). A pure dipole field has FCF 0.8 (Christensen et al., 2010), while the average IGRF from 1900 to 2015 truncated to degree 5 gives 1.12 ± 0.05. The mean FCF is 1.70 (1.64) over the past 100 ka (10 ka) from GGF100k, 1.51 from CALS10k.2 and 1.44 from HFM.OL1.A1. More complex morphologies at the CMB and associated increased FCF occur over the past 100 ka during times of instabilities (geomagnetic excursions): 29 ka BP, ∼42 ka BP (before the Laschamp excursion), 64 ka BP and toward the oldest part of the model (92-100 ka BP). FCF has very similar trends to the globally averaged P i variations (cf. to Figure 9), because both can be seen as a measure for field complexity.

Summary and Perspective
In this final section we attempt a brief summary highlighting recent progress in our knowledge of geomagnetic field evolution over the 0-100 ka interval, and suggest some likely future developments.
There is no question that a steadily improving number of available paleomagnetic secular variation data is able to provide a greatly enhanced global view of geomagnetic field variations and useful comparisons with both longer and shorter-term field descriptions. Several different representations of the field now exist in the form of global data stacks and models, but GGF100k is the first global time-varying SH field model spanning 0-100 ka. The ADM exhibits a broad range of values from a low of 2.9 ZAm 2 during the Laschamp excursion for the high-resolution LSMOD.2 model to around 98.9 ZAm 2 for GGF100k about 1,500 years ago.  Christensen et al. (2010) as "Earth like" within a factor of * , numbers noted with a star. These factors are not standard deviations as for the models, but they denote a tolerance range (from E / * to E × * ). Model periods: 1900-2015 for IGRF12, past 10 ka for CALS10k.2 and HFM.OL1.A1 and past 100 ka for GGF100k.
When averaged over the 0-100 ka time interval (coefficients available at https://earthref.org/ERDA/2382) the GGF100k has zonal quadrupole, g 0 2 , and octupole, g 0 3 , contributions that are not dissimilar to previous Holocene and 0-5 Ma averages. At Earth's surface, zonal contributions at each SH degree, l, dominate higher order terms m > 0, but in agreement with earlier studies the GGF100k average clearly exhibits significant (if individually poorly determined) nonzonal contributions to the field. Large equatorial dipole contributions are a somewhat surprising feature of the 100 ky average. PSV activity is generally higher in the Southern Hemisphere relative to the Northern Hemisphere, and is greater on average in the Atlantic Hemisphere than the Pacific Hemisphere. At several times during the past 100 ky the field morphology is rather similar to that at present, with low surface field strength over the South Atlantic or South America and similarly located reverse flux patches at the CMB.
Comparisons with higher-resolution results are available. While the recent end of the GGF100k model is generally consistent with available Holocene models, some notable differences exist. Since the Holocene models have higher resolution and are constrained by more paleomagnetic records, these models remain preferred for those timescales. Likewise LSMOD.2, the most recent effort to reconstruct the field during the 50-30 ka interval, indicates the potential benefits of better age control and high-quality data. However, there remains considerable scope for improvement in both spatial and temporal coverage by the data. It is to be expected that continued effort in gathering high quality paleomagnetic data and the development of improved, up to date, and accessible age databases will provide continued enhancements to our views of the 0-100 ka field. The fact that global mapping is now possible for the evolution of all magnetic field elements at both Earth's surface and the CMB allows visualization of continuous changes in the field (download animations from http://earthref.org/ERDA/2384/). New paleomagnetic models spanning the past 100 ka like GGF100k can place constraints on the processes that generate the geomagnetic field and give insight into geomagnetic field behavior during transitional times. The synthesis presented here suggests that a single global excursion (the Laschamp) occurred in the interval studied. PSV activity is heightened at various locations during several other time intervals but only exceeds the critical value of 0.5 globally during a limited time around 41 ka. Other recorded events appear in limited locations and are likely regional and/or transgressive in nature. Improved chronological information will be a key factor in acquiring a better determination of the progression of regional geomagnetic excursional activity.
Cosmogenic nuclide records also provide a proxy measurement for past geomagnetic DM after calibration by absolute paleointensity data. The agreement between ADM variations predicted from GGF100k model and the cosmogenic nuclides production record is generally good. Some mismatches in timing and amplitude of the variations further highlight intervals where new and improved data could help resolve inconsistencies.
We have discussed several quantities derived from GGF100k that can describe general characteristics of the geomagnetic field over multimillennial timescales and might serve to assess whether numerical geodynamo simulations appear Earth like. Our results are compatible with the broad intervals suggested by Christensen et al. (2010) for Earth-like properties, indicating strong average dipolarity with low values of AD/NAD during excursions; a tendency toward equatorial asymmetry in power in the nondipole field with occasional spikes in O/E values; larger values of the ratio of zonal to nonzonal power than expected from a white spectrum at the CMB. The FCF is probably strongly affected by limited spatial resolution, but varies over time and is positively correlated with the PSV index. The new models present the possibility of reassessing the criteria for numerical dynamos to be Earth like in the context of a much longer time interval than has been previously available.
The emerging global view on geomagnetic field variations over the past 100 kyr provided by models like GGF100k can also be useful beyond geomagnetism and paleomagnetism. Magnetic field variations on all timescales can aid in assigning chronologies to sediment cores, for example, for climate or cosmogenic isotope production studies, or samples from lava flows to study the volcanic eruption history. Global models and stacks can be useful tools in this regard, although some caveats concerning spatial correlation length of geomagnetic variations and temporal resolution exist (see recent reviews by Roberts et al., 2013;Korte, Brown, Gunnarson, et al., 2019). SH models allow estimation of geomagnetic shielding for cosmogenic and in situ isotope production, or determination of the paleomagnetosphere for space climate studies. Accuracy and reliability of empirical models depend on density of data coverage, reliability of recovered paleomagnetic signal from sediment and volcanic data, and temporal resolution and quality of independent age control in the data. Future improvements to models like GGF100k can be expected, contributing to further progress in all these fields.

Glossary
ADMAxial dipole moment, the magnetic moment of the contribution from a pure dipole field aligned with Earth's rotation axis to the geomagnetic field. This cannot be determined from individual observations, but is easily determined from the first coefficient in SH field representations. Global stacks of VADM records are considered to give a reasonable approximation of the ADM; see sections 5 and 6. CMB Core-mantle boundary; relevant as the upper boundary of the geodynamo. SH global field models can conveniently be downward continued under the assumption of an insulating mantle, to provide an indication of geomagnetic field morphology at the top of its region of origin. DM Dipole moment of the geomagnetic field, in general determined from the first three SH coefficients of empirical field models. The dipole moment cannot be directly determined from individual observations. Global stacks of VDM records are considered to give a reasonable approximation of the true DM, see sections 5, 6. DP Dipole power; energy in the (tilted) dipole field contribution. DRM Depositional or detrital remanent magnetization; sediments acquire their magnetization via this process, when magnetic moments of grains are statistically aligned with the geomagnetic field during deposition. F Abbreviation used for geomagnetic field intensity (Force). GAD Geocentric axial dipole; mostly referring to the hypothesis that over long time intervals the geomagnetic field averaged to this simple field geometry. The validity of this hypothesis is unclear, see section 2 GEOMAGIA50 Database for archeomagnetic and paleomagnetic data and metadata for the past 50 kyr (http://geomagia.gfz-potsdam.de/index.php) Geomagnetic coordinates Geomagnetic latitude ( g ) and longitude ( g ) in a coordinate system aligned with the geomagnetic dipole axis instead of Earth's rotation axis. The geomagnetic equator thus is (presently) tilted with respect to the geographic equator. Geomagnetic dipole field Over much of its geological history the geomagnetic field is well approximated by a magnetic dipole field (at present accounting for >93% of the field observed at Earth's surface), that is tilted slightly with respect to Earth's rotation axis (at present ∼10 • ). The geomagnetic dipole field is described by DM (see above) and tilt of the dipole axis, given by the geographical position of its geomagnetic pole ( D , D ). Geomagnetic field model This term generally refers to an SH model, either for one epoch or nowadays often spanning a certain time interval using temporally continuous basis functions. An overview of existing geomagnetic field models from 10 ka to 5 Ma is given in section 5.

10.1029/2019RG000656
Geomagnetic pole Pole of the approximated geomagnetic dipole field (see above). Geomagnetic North and South Poles lie diametrically opposite on Earth. Their locations ( D and D ) cannot be determined from individual observations, but are readily obtained from the first three coefficients of SH field models (see "geomagnetic field model" and "SH"). Globally averaged VGPs (see below) from individual observations are generally considered to approximate the geomagnetic pole locations. See also section 2. GGF100k The first and at the time of writing only global geomagnetic field model spanning the past 100 ky. See section 5.2. IAGAThe International Association of Geomagnetism and Aeronomy. http://www.iaga-aiga.org. IGRF The International Geomagnetic Reference Field updated every 5 years by a Working Group of IAGA (https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html) MagICComprehensive database and archive by the Magnetics Information Consortium for all kinds of archeomagnetic and paleomagnetic data and metadata. https://www.earthref.org/MagIC Magnetic Pole Magnetic pole or dip pole is the term used for the location where the magnetic field lines penetrate the Earth vertically. In contrast to the geomagnetic poles (see above), the locations of the magnetic poles are influenced by nondipole field contributions, and northern and southern magnetic pole do not lie diametrically opposite on Earth. NAD Non-axial-dipole; referring to all field contributions apart from the axial dipole. ND Nondipole; referring to all field contributions apart from the tilted dipole. NDP Nondipole power; energy in the nondipole field contributions. P i Index to characterize the strength of PSV, considering dipole field strength and VGP latitude. See section 2. PINT Database for paleointensities older than 50 ka. http://earth.liv.ac.uk/pint/ PSV Paleosecular variation; change of the geomagnetic core field over geological times. RPI Relative paleointensity; the remanent magnetization that can be measured in sediments is influenced not only by geomagnetic field strength, but also by several environmental factors. These might be normalized for, but the result only produces relative paleomagnetic field intensity variations, the RPI, not absolute field values. SH Spherical harmonic; referring to the basis functions of the most widely used empirical global field modeling method. See section 4.2. SH coefficients The coefficients determined by fitting data by SH functions; they define a global SH field model, from which predictions for all field components at any location on Earth can be obtained. They also allow investigation of global geomagnetic field characteristics like DM or symmetries. See section 4.2. TAF Time-averaged field; mostly referring to empirical global models of the geomagnetic field averaged over a certain, usually geological, amount of time. TRM Thermoremanent magnetization; volcanic and archeomagnetic materials acquired their magnetization via this process, when individual magnetic moments that freely fluctuate above the Curie temperature, become statistically aligned with the geomagnetic field below the Curie temperature, when the material cools down. VADM Virtual axial dipole moment; this can be determined from individual intensity field measurements under the assumption that the geomagnetic field was a pure dipole aligned with Earth's rotation axis. When comparing field intensities from different locations, VADMs the systematic intensity latitude dependence from the dipole field contribution. See section 2. Global and regional stacks of VADM records are described in section 5. VDM Virtual dipole moment; this can be determined from individual intensity and inclination measurements under the assumption that the geomagnetic field was a pure dipole, possibly tilted against Earth's rotation axis. When comparing field intensities from different locations, VDMs eliminate the systematic intensity latitude dependence from the tilted dipole field contribution. See section 2. VGP Virtual geomagnetic pole; this can be determined from individual declination and inclination measurements as the pole of a pure, tilted dipole field. Geographic latitude and longitude of a VGP here are labeled V and V , respectively. See section 2.