Global Evolution and Dynamics of the Geomagnetic Field in the 15–70 kyr Period Based on Selected Paleomagnetic Sediment Records

Reconstructions of the geomagnetic field on long time scales are important to understand the geodynamo processes in the Earth's core. The geomagnetic field exhibits a range of variations that vary from normal, dipole‐dominated secular variation to geomagnetic excursions and reversals. These transitional events are associated with significant directional deviations and very low intensities. Here we present a new, global geomagnetic field model spanning the period 70–15 ka (GGFSS70) that includes three excursions: Norwegian‐Greenland Sea, Laschamps, and Mono Lake/Auckland. The model is built from nine globally distributed, high‐resolution, well‐dated, sedimentary paleomagnetic records. The GGFSS70 indicates that the axial‐dipole component changed sign for about 300 years in the middle of the Laschamps excursion (41.25–40.93 ka). The energy comparison at the Earth's surface reveals that the axial‐dipole energy is always higher than the non‐axial‐dipole except over the Laschamps. In the other two excursions, the axial‐dipole is reduced by about one order of magnitude for the Norwegian‐Greenland Sea excursion and less for the Mono Lake/Auckland. At the core‐mantle boundary, the large‐scale non‐axial‐dipole power is comparable to the axial‐dipole power, except over the excursions when the axial‐dipole decreases, though less clearly for the Mono Lake/Auckland excursion. The axial dipole moment over the 15–70 ka varies from 0 to 8 × 1022 Am2, with an average and standard deviation of 5.1 ± 1.5 × 1022Am2. The Laschamps excursion is associated with growth and poleward movement of reversed flux patches and reversed field in the tangent cylinder at the excursion midpoint, which is not the case for the other two excursions.


Data
To derive a model with the highest possible temporal resolution, paleomagnetic sediment records were selected for high paleomagnetic data quality, temporal resolution (high sedimentation rate), and good independent age control, while reasonable global data coverage also was a criterion. An overview of the nine globally distributed records that we considered suitable for the period 15-70 ka is given in Table 1 and plotted in Figure 1. All records have the three components, intensity, and directions, and the total number of data points is 23,197  The average sedimentation rate (SR) is also a mean value over the 15 to 70 ka period from all cores. Details for the cores are given in the text. c Considered only part of the record with ages 45 ka. Details are given in the text.

Table 1
Data Set Used to Construct the GGFSS70 Model Figure 1. (a) Spatial distribution of the nine sediment records used in the study (black triangles). The color map represents how well the core-mantle boundary is sampled by these sites; (b) Temporal distribution of data, by individual components and combined, plotted in 100 years bins. References for the paleomagnetic sediment records are provided in the text and Table 1. 10.1029/2021JB022681 4 of 18 (8,095 declinations, 8,059 inclinations, and 7,043 relative paleointensities [RPI]). Locations of the nine sediment records selected to constrain the model are plotted in Figure 1a. The background color represents the globally summarized data kernel, that is, the averaged sampling of the field at the CMB over the studied period (e.g., Johnson & Constable, 1997) by these data. Although the number of records is limited, no region is completely unrepresented (blue shades), but some regions are poorly sampled. The midlatitude regions in both hemispheres are best-sampled, with maxima in South American and West Pacific regions reflecting the high number of data in the two records from these regions. The whole period of interest is covered by all three components (Figure 1b). Decreased number of data that might influence the reliability of the model for those times are evident in the 62-64 ka interval and toward both ends. The average number of data points per 100 years is 41. Information about the paleomagnetic measurements and chronology are provided in Table 2.
We updated the age models of six records and the references are listed in Table 2. The age models of three records are updated from the GISP2 to GICC05 reference curve using the conversion function in Obrochta et al. (2014). The adjustments to GICC05 introduce offsets that vary in time, prior to ∼38 ka, the changes are smaller than 500 years, ±100 years at the Laschamps midpoint, and increasing after, when the radiocarbon is no longer available, to ∼2,500 years at 50 ka (shift to younger ages) and a maximum shift of ∼4,000 years at 70 ka (shift to older ages). Thompson and Goldstein (2006) calibrated the SPECMAP δ 18 O reference curve (Imbrie et al., 1984;Martinson et al., 1987)   (AICC2012) and further tuning using coincident warming events provide a revised reference for sediment cores from the Southern Ocean (Anderson et al., 2021). The new AICC2012-tuned age model shifts the Last Glacial Maximum (LGM) (18.25 ka) and the MIS-3 to MIS-4 boundary (56.75 ka) (Lisiecki & Raymo, 2005) to ∼330 years later and ∼820 years earlier, respectively. On the other hand, radiocarbon-based age models need to be converted to calendar ages, and this process requires updated values for marine radiocarbon reservoir age (e.g., Reimer et al., 2013). Having in mind these changes in the age reference curves, one realizes the pivotal role of age models in properly recovering the PSV. In the following, we briefly discuss each paleomagnetic record used in the study.
The Black Sea (BS) stack was obtained by averaging 16 individual cores, 10 MSM33 (initially published in Nowaczyk et al., 2018;Liu et al., , 2019, and six M72/5 cores (Nowaczyk et al., 2012(Nowaczyk et al., , 2013. Details about the stack and comprehensive analysis of the geomagnetic excursions as seen in that location can be found in Liu et al. (2020). Age models of M72/5 cores are based on radiocarbon dating, two tephra layers, and correlation of sedimentological parameters to the Greenland ice core (NGRIP) oxygen isotope (δ 18 O) record on the GICC05 age scale (Svensson et al., 2008). Regarding the MSM33 cores, age models were also built by correlating X-ray fluorescence (XRF) logs to the NGRIP record, using the same ratios as for the M72/5 cores. For the time interval between 70 and 15 ka, the sedimentation rates of Black Sea cores range between 10 and 38 cm/ka, with half of the above 20 cm/ka. The average mean sedimentation rate of these cores is 22.7 cm/ka. For the Laschamps interval (39-43 ka), sedimentation rates range between 5 and 31 cm/ka, and 6 out of 12 cores that cover this excursion have sedimentation rates above 15 cm/ka. In general, looking at individual cores, stadial (cold) time intervals have higher sedimentation rates than interstadial (warm) times.
ODP-1233 is a high-resolution record from the southwest margin off Chile in the southeast Pacific Ocean (Lund et al., 2006(Lund et al., , 2007. Magnetic components are obtained from measurements of core halves after a maximum of 25 mT alternating magnetic field (AF) demagnetization. RPI is estimated as the ratio of the natural remanent magnetization (NRM) and magnetic susceptibility. The age model of this record is based on the correlation of magnetic susceptibility and Ca concentration to nearby core GeoB 3313-1 (up to 7 ka, Kaiser et al., 2005), 14 C ages (to 40 ka, Lamy et al., 2004) and correlation of the sea surface temperature to the δ 18 O record of Antarctic Byrd ice core (to 70 ka, Kaiser et al., 2005). We updated the chronology according to Chase et al. (2014) with recalibrated 14 C ages using the Marine13 calibration curve (Reimer et al., 2013). The four individual cores (1233B, C, D, and E) of this record are converted to the new age model. The record exhibits an exceptionally high mean sedimentation rate of 150 cm/ka in the 15-70 ka period.
The JPC-14 record comes from the Blake Outer Ridge, western North Atlantic Ocean (Lund et al., 2005). Samples within the Laschamps excursion were obtained with AF demagnetization up to 100 mT, and AF demagnetization at 20 and 60 mT for the remaining intervals. The directional variations of the discretely sampled paleomagnetic records of JPC-14 and CH89-9P from the Bermuda Rise (1,200 km apart) show excellent agreement . The chronology is based on radiocarbon dates (on GPC cores in the same region) from Keigwin and Jones (1994) that are correlated to JPC-14 through magnetic susceptibility and calcium carbonate variations . Here, we converted this age model, which was based on the GISP2 reference record, to the GICC05 time scale following Obrochta et al. (2014). The average sedimentation rate over the oxygen isotope stage 3 is 35 cm/ka, and ∼28 cm/ka during the period of the Laschamps excursion.
A full-vector paleomagnetic record PLC08-1 has been obtained from Pyramid Lake, USA to study the Laschamps and Mono Lake/Auckland excursions (Lund, Benson, et al., 2017). The GISP2 chronology is based on 20 radiocarbon dates, four tephra layers, and PSV correlation with records from the western North Atlantic, JPC-14, CH89-9P, and CH88-10P . The PSV correlation confirms the radiocarbon and tephrochronological age model, and adding the PSV tie points did not significantly change the latter model. For this study, we converted the GISP2 age model to GICC05 (Obrochta et al., 2014). All samples have been AF demagnetized at max. 60 mT. Additionally, selected samples have been AF demagnetized at max. 100 mT. A few normalizers for obtaining the paleointensity variations had been tested, and the final record was derived as the ratio of NRM to saturation isothermal remanence (SIRM).
The paleomagnetic record MD98-2181 from the Philippines, western Equatorial Pacific Ocean, has an exceptionally high sedimentation rate (on average 53 cm/ka) and represents one of the highest resolution records from the generally low sedimentation region of the Pacific (Lund, Schwartz, & Stott, 2017). The record has been obtained by u-channel sampling with AF demagnetization of max. 20 mT. The chronology reported in Lund, Schwartz, and Stott (2017) is based on radiocarbon and GISP2 oxygen isotope stratigraphy. The number of dates is larger than the originally published RPI record from the same location (Stott et al., 2002). We kept the age model for the period from 28 ka to present, which is based on the revised calibrated radiocarbon dating (Khider et al., 2014), but updated the ages beyond 28 ka to the GICC05 age model following Obrochta et al. (2014). The record on the old and updated age scale is plotted in Figure S4 in Supporting Information S1.
Core MD94-103 from the southern Indian Ocean (Laj et al., 2006;Mazaud et al., 2002), measured using u-channels and covering the period 23-52 ka, provides a pronounced record of the Laschamps excursion in all three components. The preliminary age model has been established by comparing an oxygen isotope record from a neighboring core to the magnetic susceptibility record of MD94-103. The age model was refined with tie points to the NAPIS-75 paleointensity stack (Laj et al., 2000), listed in Sicre et al. (2005), while the younger part of the core, not used here, is based on radiocarbon dating. In this study, we used a revised chronology, namely the AIMtuned age model of Anderson et al. (2021). First, the record is aligned to the Antarctic temperature record on the AICC2012 timescale and further tuned based on the existence of Antarctic Isotope Maximum (AIM) event type warmings (see Figure S5 in Supporting Information S1 for the record plotted on the old and updated age scale).
To complement the Indian Ocean region with a longer record than MD94-103, we selected MD84-528 from the southern Indian Ocean (Tric et al., 1992). The original age model based on oxygen isotope stratigraphy (Imbrie et al., 1984;Martinson et al., 1987) was updated with the help of U/Th calibrated δ 18 O event boundaries from Thompson and Goldstein (2006). Between the MIS stage boundaries, a linear interpolation was applied. This calibration follows Brown et al. (2018) and a shorter version of the record (30-50 ka period) was used in the LSMOD models Korte et al., 2019). Both the original and calibrated MD84-528 records are plotted in Figure S6 in Supporting Information S1. The difference between the records is obvious in all three components. For instance, the difference at the Laschamps excursion is on average about 2.5 ka. The updated model shifts this excursion toward older age, thus making it more consistent with the globally observed age of about 41 ka.
MD07-3076Q from the South Atlantic Ocean is used in this study as originally published by Channell et al. (2017). U-channel samples are AF demagnetized at max. 100 mT, and anhysteretic remanent magnetization (ARM) and isothermal remanent magnetization (IRM) normalizers were tested for obtaining the RPI estimates. The two ratios, NRM/ARM and NRM/IRM, produced similar results for the period of interest here. The age model was built with radiocarbon dates and via correlation of sea-surface temperature proxies to Antarctic ice core (EPICA) records, placed on the AICC2012 chronology (Veres et al., 2013).
The MD04-2822 core from the Rockall Trough (NE Atlantic) covers the past ∼200 ka (Channell et al., 2016). The directional record indicates an excursion at 26 ka and the signals over the other excursions, Laschamps and Mono Lake/Auckland, are not pronounced. When the younger part of this record (<45ka) is considered, the model fails to recover the Laschamps excursion in the region. Additionally, there is twisting of the declination record. Since the North Atlantic region is partly covered by the JPC-14 record, only the older part of the record (45-70 ka) is included in the model. The age model is based on a range of dating techniques: radiocarbon, tephra, correlation to NGRIP δ 18 O record using the updated ages on the GICC05 timescale (Rasmussen et al., 2014) and to the LR04 δ 18 O stack (Lisiecki & Raymo, 2005). For the period we considered (45-70 ka), the age model is constructed by tie points obtained from the correlation to NGRIP δ 18 O, one tephra horizon at 55.38 ka, and one point at 70.15 ka based on correlation to LR04 stack.

Methodology
The new model, named GGFSS70 (Global Geomagnetic Field model from Selected Sediments for the past 70 ka), follows the methodology used to build CALSxk and GGF100k models (e.g., Constable et al., 2016;Korte et al., 2009;Panovska, Constable, & Korte, 2018). The model is built by inversion using spherical harmonic functions in space (up to degree and order six) and cubic B-splines in time (50 years knot spacing). An L2 measure of misfit to the data is employed. For the regularization in time and space, we used a norm based on the second time derivative of the radial magnetic field integrated over the CMB and over the model period, and the Ohmic heating norm, respectively. Temporal and spatial smoothing parameters are selected by trade-off curves, choosing the simplest model that explains the data within a reasonable uncertainty (both trade-off curves are provided in Figure S7 in Supporting Information S1). Although the parameters are chosen subjectively from the knee of the curve, models close to this point have very similar norms and misfits, and exhibit similar structures. It was not necessary to use the smoothing kernels introduced into the forward modeling for the GGF100k model (Panovska, Constable, & Korte, 2018) because of the high-resolution sediment records that constrain the GGFSS70 and their global distribution. Temporal resolution analysis (for more details, see, e.g., Panovska et al., 2012, and one example record in Figure S8 in Supporting Information S1) reveals that 6 of the 9 records have mean smoothing times less than 200 years, with an average of 166 years from all locations ( Figure S9 in Supporting Information S1). The mean and standard deviation of temporal resolution for all records are listed in Table 2.
Relative components, intensity, and declination, are calibrated to absolute values before the modeling. For the calibration of RPI, absolute Virtual Axial Dipole Moment (VADM) data are compiled from the GEOMAGIA50. v3 database (<50 ka; Brown, Donadini, Korte, et al., 2015) and IAGA Absolute Paleointensity (PINT) database (>50 ka; version 2015.05, Biggin et al., 2009Biggin et al., , 2010Perrin & Schnepp, 2004) to cover the period 15-70 ka. The scaling factor is estimated as the ratio of the absolute VADMs average and the average of RPIs from the individual records converted to VADMs. This ratio is further refined by adjusting the scaled VADM distributions of individual records to VADM distributions of absolute paleointensities. Relative declination of the individual records is set to zero mean over the length of the records, excluding transitional directions defined using Fisher's (1953) statistics as directions outside of a 35° circle around the direction expected from a geocentric axial dipole. A detailed description of the calibration process and example plots are available in the Section S1 in Supporting Information S1.
All sediment records are equally weighted with 5 μT uncertainties for the intensities and α 95 of 8.5° for directions. The α 95 values are converted to standard deviations of directional data using the revised equations in Suttie and Nilsson (2019). The inclination uncertainty estimate for all records is 3.5°. Declination uncertainties depend on inclination values (see the example of the Black Sea record for time variations of declination uncertainty in Figure S10 in Supporting Information S1). The large deviations in inclination during excursions produce very large declination uncertainties, for instance, in the Black Sea record, the average declination uncertainty is 6.18° but the uncertainty at few excursional instances gets to about 19°. Therefore, we opted for equally weighting all declination data in one record with the mean value estimated over the studied period, namely 6.18° for BS stack, 5.83° for ODP1233, 5.30° for JPC-14, 6.03° for PLC08-1, 3.54° for MD98-2181, 8.17° for MD94-103, 8.65° for MD84-528, 6.59° for MD07-3076 and 11.18° for MD04-2822. In general, the declination uncertainty increases with increasing latitude ( Figure S11 in Supporting Information S1). The final normalized misfit of the model is 1.03% and 1.3% of the data are rejected based on a 5 standard deviations rejection rule. Model predictions to each record are plotted in Figures S12-S20 in Supporting Information S1.

Time Variations and Transitional Periods
In the following, we analyze several aspects of the model, starting with time variations of the three dipole coefficients, derived dipole moment and dipole axis latitude, as well as the paleosecular variation index (P i , Panovska & Constable, 2017) ( Figure 2). The time evolution of the individual dipole, quadrupole, and octupole terms is plotted in Figure S21 in Supporting Information S1. The GGFSS70 model suggests that the equatorial dipole coefficients are of similar magnitude and vary without large fluctuations from the mean value during the whole 70-15 ka period (mean 1 1 ∶ 0.17 ± 2.12 ; ℎ 1 1 ∶ 0.85 ± 1.56 ) . However, the axial dipole term (Figure 2) varies significantly (mean 0 1 ∶ −19.69 ± 5.76 ) . Moreover, it not only reduces to zero but also reverses for about 300 years at the time of the Laschamps excursion (41.25-40.93 ka). This short reversal is reflected by a double-dip of the dipole moment at the excursion mid-point. Brown and Korte (2016) showed this double-dip situation by exploring simple excursion scenarios. Indeed, seven of the nine records considered in this study show a double-dip in their intensity data over the Laschamps. The three distinct minima of the absolute values of the axial dipole g 0 1 correspond to Norwegian-Greenland Sea (NGS), Laschamps, and Mono Lake/Auckland excursions, denoted in Figure 2. These three excursions are also recognized in the dipole moment, dipole axis latitude, and averaged PSV index over the globe. The Laschamps excursion has the lowest dipole moment of 0.12 × 10 22 Am 2 at 40.97 ka, followed by the Norwegian-Greenland Sea excursion (1.87 × 10 22 Am 2 at 64.87 ka), seen for the first time as a pronounced event globally, and the Mono Lake/Auckland excursion with a relatively higher dipole moment, 3.30 × 10 22 Am 2 at 34.47 ka. Two other periods of comparatively low dipole moment are observed after the NGS excursion at about 59 ka, and after the Mono Lake/Auckland excursion at 29 ka. During the former, the PSV index increases to slightly above the transitional threshold value, whereas only a moderate increase in the PSV index is found for the latter. The dipole moment lows at 60-58 ka result from the very low intensities in the ODP1233 record. Maximum dipole moment values around 8 × 10 22 Am 2 , similar to the present-day, are observed before the NGS and Laschamps excursions, around 67.35 and 50.26 ka, respectively ( Figure 3).
The model dipole moment (DM) agrees well with the high-resolution model of the Laschamps and Mono Lake/ Auckland excursions (50-30 ka), LSMOD.2 , and the general trends in the lower resolution GGF100k model (Panovska, Constable, & Korte, 2018), with some notable differences over the 70-50 ka and 30-15 ka periods (Figure 3). Table S1 in Supporting Information S1 provides minimum, maximum, mean, and standard variations values of the three models. The new model clearly shows the influence of data selection for  reconstructing the geomagnetic field, a strictly selected data set of nine versus a considerable global data compilation of more than 100 records of variable resolution and data quality. Over the period 30-50 ka, the average axial dipole moment of GGFSS70 is in a good agreement with LSMOD.2, but lower than the GGF100k average over the 15-70 ka interval (Table S1 in Supporting Information S1). This difference can be attributed to the different data sets, as 4 out of the 9 records in the GGFSS70 model are new and were not used in GGF100k, updated records and age scales, and different calibration methodologies, with calibration within the inversion and before modeling with the help of archaeomagnetic/volcanic data, for GFF100k and GGFSS70, respectively. In addition, there is on average 7 μT difference between the calibrations of the four intensity records used in both models (Black Sea, JPC-14, MD98-2181, and MD84-528). In particular, GGFSS70 supports the nature of the Mono Lake/Auckland excursion as a double event at 34 and 29 ka (see Figure S22 in Supporting Information S1 for the VDMs of individual records and the GFFSS70 dipole moment over this time interval). Moreover, in contrast to GGF100k, it also has a clear DM low around the time of the postulated GGF-28k event , about half of the present-day value at 27.70 ka, and relatively low DM values (4.7 × 10 22 at 18.70 ka) around the time of the Hilina Pali excursion.
We compared the GGFSS70 dipole moment with independently reconstructed geomagnetic dipole moment variations based on cosmogenic isotopes (Figure 4). The following records are considered: BeDM 20−60 (Simon et al., 2020), constrained by four marine sediment sequences, two from the west-equatorial Pacific, and two from the Portuguese margin, covering the period 20-60 ka; and 36 Cl GRIP and 10 Be stack of NEEM and GRIP ice cores, which span the last glacial period and are "climate-corrected" before converted to dipole moment (Zheng et al., 2021). The conversion is based on theoretical production models, and in case of the BeDM 20−60 stack also statistical calibration using absolute paleointensities measured on lava flows. Both these methods provided similar results (Simon et al., 2020). Besides usual uncertainties in these production rate records, measurement and age errors, climate/environmental influences (transport and deposition processes) on the productions rates are a real challenge for interpreting them as geomagnetic signals (Beer et al., 2012). All broad peaks, lows and highs, observed in GGFSS70 DM coincide well with the peaks in the 10 Be stack. In general, the comparison shows a better agreement during the post-than pre-Laschamps period when looking at the absolute values (Figure 4). To check if this discrepancy comes from individual paleomagnetic records, we estimated the individual virtual dipole moment (VDM) from the nine sediment records. These individual curves show much better consistency among each other in the post-Laschamps period, but a large range in amplitudes over the 55-45 ka period. At this period, the BeDM 20−60 record exhibits relatively high values that are matched only with the Black Sea record, while all others appear to vary about the present-day value (gray dashed line, IGRF-13 model, Alken et al., 2021). The calibration to absolute intensity values should not be a factor since the same paleointensity databases, Geomagia50.v3.3 and PINT2015, are used to scale the normalized 10 Be ratio stack into BeDM 20−60 (see Simon et al., 2020 for details on the selection criteria) and RPIs to absolute intensities in the GGFSS70 model (see Section S1 in Supporting Information S1). The Laschamps dipole low in GGFSS70 is lower than in all cosmogenic isotope reconstructions.  (Simon et al., 2020); 36 Cl GRIP and 10 Be stack of NEEM and GRIP ice cores, which span the last glacial period and are "climatecorrected" before converted to DM (Zheng et al., 2021). Also plotted are the VDMs estimated from the nine individual records that constrained the GGFSS70 model (light gray curves), as well as the present-day DM (gray dashed line) obtained by the IGRF13 model (Alken et al., 2021).
Another striking feature is the maximum at 60 ka observed in the 36 Cl GRIP record (but not in the 10 Be stack) when moderately low DM values are predicted by the GGFSS70 model.
To illustrate the regional differences, we plotted the VDM, virtual geomagnetic pole (VGP), and the PSV index (derived from the previous two quantities) at records locations-as representative of globally distributed points over the Laschamps excursion in Figure 5 (and for the whole period 15-70 ka in Figure S24 in Supporting Information S1). Also marked in the figure is the period when the axial dipole changed polarity, and gained about 5% of the present-day value in the reverse direction. We find an almost synchronous decay of field intensity (expressed as VDM) prior the Laschamps excursion, to the point when axial dipole (AD) reverses, and quite different behavior after when the field recovers in the normal polarity. Note also that the intensity maxima between the double minima do not necessarily occur contemporaneous with the globally reversed dipole direction, due to the regionally different non-dipole influences. The temporally asymmetric behavior is also true for the VGP latitudes, which nearly all synchronously have a strong drop to low values at the start of the AD reversal, while the subsequent changes and the duration of the recovery vary much more. If we consider the PSV index for estimating the start/end of the excursion and its duration (P i > 0.5 for excursional events), then the starting ages for all locations lie within 850 years (the range is 42.13-41.28, mean 41.47 ka). The end ages vary more strongly, translating to regionally different durations from 500 to 2,300 years at these data locations, with an average of 1,400 years. The range becomes even larger when looking at the whole globe, with durations varying between 500 and 3,380 years with a mean of 1,820 years in the time interval 42.3-38.8 ka, in good agreement with values found by Korte et al. (2019). We plotted global maps of the maximum PSV index, starting age, and duration of all excursional events observed in the model in Figure S25 in Supporting Information S1. GGFSS70 always has P i > 0.5 somewhere on Earth in this time interval. Moreover, the transitional PSV index over the Laschamps often has two or more peaks with values above or below the threshold in-between. This suggests that the field varies significantly in intensity and directions during excursional times, which sometimes results in PSV index values characteristic for normal, dipole-dominated field. The picture becomes more complex when looking at the NGS excursion ( Figure S24 in Supporting Information S1). Over a wider interval around this excursion, the field at selected locations is unstable. The PSV index indicates multiple transitional epochs over a 10 ka period, starting at 66 ka. The peaks can be grouped about 65-63 and 60-58 ka. There is no pronounced indication of the Mono Lake/ Auckland excursion in the PSV index, where values shortly exceed the threshold of 0.5 during the 34.37-33.67 ka period at one location, ODP1233. If the threshold is lowered to 0.3, a value never reached by recent and historical models that are representatives of dipole-dominated fields, then the PSV index at seven locations exceeds this threshold at 35.1-34.5 ka or later at 34.0-33.6 ka. For the Mono Lake/Auckland excursion, we analyzed the model predictions to test the presence/absence of this excursion at its eponymous locations ( Figure S23 in Supporting Information S1). The GGFSS70 model produces weakly increased, but not excursional, values at the Mono Lake location. On the other hand, the peak PSV index at the Auckland location results from synchronous lower field intensity and VGP latitude at about 34 ka. Considering the 0.3 threshold, as the 0.5 is not reached, the excursion is confined in the 34.1 ± 0.4 ka period, which is in very good agreement with the age of 34.2 ± 1.2 ka estimated by Laj et al. (2014), where the authors also suggested renaming the Mono Lake excursion to Auckland excursion.

Energy Evolution During Geomagnetic Excursions
We compare the axial-dipole and non-axial-dipole power over 10 kyr periods (5 kyr pre-event and post-event) covering the NGS, Laschamps, and Mono Lake/Auckland at the CMB and Earth's surface in Figures 6a and 6b, respectively. The non-axial-dipole power has to be understood as large-scale non-axial-dipole power (in particular when looking at the CMB) because the effective spatial resolution of the model does not go beyond SH degrees 3-4. The resolution estimates were obtained by comparing the spectra of GGFSS70 with available historical, Holocene, and long-term models ( Figure S26 in Supporting Information S1) and synthetic tests (Section S12 in Supporting Information S1). The comparison to the historical model indicates a spatial resolution of GGFSS70 up to degree 4. Temporally, paleomagnetic models cannot reach the resolution of models obtained from directly observed data. GGFSS70 has lower temporal resolution than LSMOD.2, but higher compared to GGF100k for all degrees. We have also performed synthetic tests to assess the level of resolved features with our limited data set (details are available in the Supporting Information and Figures S1 and S28 in Supporting Information S1). These synthetic tests showed that the effective spatial resolution varies in time and the maximum is degree 4.
The plots in Figure 6 clearly show the difference in the extreme axial-dipole decay over the Laschamps excursion compared to the NGS and Mono Lake/Auckland excursions, which is likely the reason for the global character of the Laschamps compared to the latter two. At the CMB, the large-scale non-axial-dipole energy is comparable to the axial-dipole except over the central excursional intervals, when the dipole decreases, though less clearly for the Mono Lake/Auckland excursion. In all cases, the non-axial-dipole energy fluctuates roughly within the same range over excursional and non-excursional times, in agreement with results found by Brown et al. (2018) and Korte et al. (2019). The Black Sea transformed data in components parallel to the direction expected from GAD and two components perpendicular to it, representing the non-GAD components, showed the same behavior (Liu et al., 2020). Apart from the Laschamps excursion, the axial-dipole energy at Earth's surface is always higher than the non-axial-dipole, even during the dipole minima. In the two less pronounced excursions, the axial dipole reduces about one order of magnitude for the NGS excursion and less for the Mono Lake/Auckland.
In comparison with the LSMOD.2 model ( Figure S1 in Supporting Information S1), over the 5 kyr directly preceding the Laschamps, the GGFSS70 suggests less dipole dominance than LSMOD.2, but some differences in non-dipole power might be due to spatial resolution of the models, as non-dipole power would be higher anyway if we could resolve smaller spatial scales. For the 5 kyrs after the Laschamps, dipole and large-scale non-dipole contributions vary at a similar level, in closer agreement with LSMOD.2. The dipole and non-dipole energies in GGFSS70 and LSMOD.2 agree relatively closely, except for the reduction of dipole power in the excursion midpoint and the non-dipole low 300 years pre-Laschamps in the GGFSS70 model.

Morphology of Geomagnetic Excursions
As the axial dipole moment is diagnosed to be the most variable component during geomagnetic excursions, we look at maps of Zcosθ over the CMB (Figure 7), which represents features that contribute to the dipole moment. Z is the downward vertical component of the magnetic field and θ is co-latitude. These maps are plotted in equal-area projection, centered at the North and South poles, and show the normal (in orange) and reverse (in blue) flux patches that contribute to the axial dipole moment (ADM). Maps of the radial field component at the CMB are shown in Figure S1 Supporting Information S1.
The ADM is proportional to Zcosθ integrated over the CMB (Gubbins, 1987). Its decrease during the Laschamps excursion is a consequence of the growth and poleward movement of two reverse flux patches, one over Siberia and the second over North America in the Northern hemisphere. In the decay period (∼42 − 41 ka) , first the Siberian, and later the North American reverse flux patch cross the tangent cylinder (TC). At the times of the two dipole minima (Figures 7c and 7e), the radial field in half of the TC area in both hemispheres is in the opposite direction, while the whole northern TC has reversed flux at the mid-point when 0 1 has changed sign ( Figure 7d). In the Southern hemisphere, there is reverse flux drifting constantly into and out of the TC, covering only half of the area also at the excursion mid-point (Figure 7d, lower panel). In contrast, for the NGS and Mono Lake/ Auckland excursions, normal flux dominates in both hemispheres also during the lowest ADM times (Figures 7b  and 7f, for NGS and Mono Lake/Auckland excursions, respectively). Even though the NGS excursion has a quite low ADM, ∼22% of the present-day value of 7.6 × 10 22 Am 2 , the normal/reversed flux intensity, and distribution differ from the Laschamps excursion. This is also true for the Mono Lake/Auckland excursion, with its ADM minimum at 42% of the present-day field.
Maps presenting the field intensity at the Earth's surface at intervals of 250 years covering the three excursions, NGS, Laschamps, and Mono Lake/Auckland, show the differences in the global field decrease over the three periods (Figure 8, all maps are plotted with the same intensity scale). During the Laschamps excursion (Figure 8b), the field intensity is globally very low for ∼1,000 years. The intensity is also decreased globally for the NGS excursion (Figure 8a), but the minimum field intensity at Earth's surface is more limited to equatorial and mid-latitudinal regions. The intensity distribution over the Mono Lake/Auckland excursion at 34 ka keeps a dipolar structure at the Earth's surface even more clearly during the global minimum. The pre-, mid-, and post-excursional morphology of the three excursions, however, is quite similar. The strongest decrease in field strength starts over middle America/the equatorial Atlantic and India/the Indian Ocean, and slightly more east in the case of the Mono Lake/Auckland excursion (epochs: 65.75-65.50 ka for NGS; 41.75-41.50 ka for La; and 35.5-35.0 ka for ML). The weakest field during the middle of the excursions is mostly found in midlatitudes (mainly northern) of the Atlantic/Indian Ocean sectors (65.25-64.50 ka for NGS; 34.75-34.50 ka for ML); however, the field is low globally for the Laschamps (41.25-40.50 ka). The recovery to stronger field strength appears to follow a similar geometry for the Laschamps and the Mono Lake/Auckland excursions (40.25-40.00 ka for La; 34.0-33.50 ka for ML). The CMB field morphology associated with the field intensity maps in Figure 8 is plotted in Figure S1 in Supporting Information S1. These figures emphasize the common process that occurs in the three cases related to the field intensity decrease, reverse flux appearing in low-mid latitudes, and moving polewards. In NGS and Laschamps excursions, the patterns when the field intensity decreases and recovers look very similar. Besides the flux positions and motions, the flux intensity is important. What is evident is the intense flux patches that persist throughout the Mono Lake/Auckland excursion, moderate levels for the NGS excursion, but no intense features during the period of the Laschamps' minimum.

Conclusions
We have built a magnetic field model covering the period 70-15 ka based on nine paleomagnetic records, selected to have all three components, high-resolution, good age control (updated when possible), and best possible global distribution. The effective spatial resolution of the model is up to degrees 3 to 4. Relative components have been calibrated before the modeling and all data have been equally weighted. The model includes at least three well-known geomagnetic excursions, the Norwegian-Greenland Sea, Laschamps, and Mono Lake/Auckland. We investigated the differences and similarities of these excursions, globally and regionally, at the Earth's surface and the CMB. The model provides the first global reconstruction of the Norwegian-Greenland Sea excursion (65 ka) with a clear low in the dipole moment and excursional values of the PSV index, which quantifies both the intensity and field directions. The Mono Lake/Auckland (34.5 ka) shows only a slight increase in the PSV index, not reaching the threshold value for a global excursional event, though regionally the threshold is exceeded. Two additional lows in the dipole moment are observed after the NGS and Mono Lake/Auckland excursions, at 59 and 29 ka, respectively. The latter is in good agreement with the second Mono Lake/Auckland feature identified in the LSMOD.2 model . Regarding the GGF-28k event identified in the GGF100k, the better resolved GGFSS70 model suggests no excursional values at 28ka, instead, low dipole moment and dipole energy are observed at 29 ka. Moreover, the limited region of transitional field variations (PSV index 0.5) is found in the south Australian region in contrast to the GGF-28ka observed in the GGF100k model in South America. The GGFSS70 model does not reflect an excursion postulated at 26.5 ka in the MD04-2822 record from the North Atlantic, but this part of the record was not used to constrain the GGFSS70 model. The most recent relatively-low dipole moment at 18.7 ka can be associated with the postulated Hilina Pali excursion, though the PSV index over this interval exhibits no transitional values, above 0.5. GGFSS70 indicates that the axial-dipole component changed sign for about 300 years in the middle of the Laschamps excursion from 41.25 to 40.93 ka, slightly shorter than the 500 years found in the previous version of the model, GGFSS70.1 (Liu et al., 2020). This reversing increases the intensity at the Earth's surface creating a double intensity low over the Laschamps excursion in many locations, as observed in high-resolution records. The Laschamps excursion is characterized by a more substantial decrease of dipole energy compared to NGS and Mono Lake/Auckland excursions. The non-dipole energy varies about the same level in pre-/post-excursional periods and during the excursion for all events. Two minima in dipole energy over the Laschamps coincide with the times when the axial dipole reverses, but the first is less pronounced by 1 order of magnitude due to a stronger contribution from the 1 1 coefficient. The first dipole power decrease is accompanied by a slight non-dipole energy decrease. GGFSS70 suggests that this is not the case in the Norwegian-Greenland Sea or Mono Lake/ Auckland excursions. It remains unclear if this non-dipole decrease accompanying the dipole decay might be a  characteristic feature for major global geomagnetic excursions, during which the axial dipole changes sign. This will have to be confirmed or rebutted with analyses of more excursions of this type. Clearly, for robustly resolving the characteristics of the geomagnetic excursions, high-quality, high-resolution paleomagnetic records with good, independent age control are needed. The new model once more confirms earlier findings that excursions may appear quite different in data records from different regions, which has to be kept in mind when considering geomagnetic excursions as stratigraphic tie points or inferring global field properties from individual records.

Data Availability Statement
Geomagnetic field maps are plotted with the programs Magmap and Color by Robert L. Parker (https://igppweb. ucsd.edu/∼parker/Software/). The manuscript includes Supporting Information file that contains all the supplementary sections and figures. An animation of the GGFSS70 model is available at the website https://earthref. org/ERDA/2471/. The model coefficients, FORTRAN codes for producing field predictions, and Gauss coefficients for a particular epoch from the time-dependent coefficients, and a file with dipole moment, dipole-axis coordinates, and PSV index can be found at https://earthref.org/ERDA/2472/.