Paleosecular Variation and the Time‐Averaged Geomagnetic Field Since 10 Ma

Investigations into long‐term geomagnetic variations provide useful information regarding paleomagnetic field behavior. In this study, we assess the latitudinal structure of paleosecular variation (PSV) and the time‐averaged field (TAF) for the Brunhes normal and Matuyama reverse chrons, and for the 0–10 Ma period, from an updated and reviewed paleodirectional database spanning the past 10 Myr. The new database comprises 2,543 paleomagnetic sites from igneous rocks, providing improvements in the geographic and temporal distributions of high‐quality data relative to previous compilations. In addition, the new data collection differs considerably in application of strict selection criteria. Statistical analysis of the virtual geomagnetic pole (VGP) dispersion curve of Model G reveals a low latitudinal dependence of PSV for the last 10 Myr. For this period, we present a zonal TAF model based on the latitudinal distribution of inclination anomaly data. The best estimates found for axial quadrupole and octupole components were about 3% and 1% relative to axial dipole component, respectively. The new statistical models for the Brunhes and Matuyama chrons have different patterns in both PSV and TAF, in compliance with earlier studies. Our quantitative assessments indicate an apparent hemispheric PSV asymmetry, particularly in the Brunhes chron, with a stronger latitudinal signature in the southern hemisphere compared to the north. These findings suggest that equatorial PSV asymmetry, that has previously been found in modern, historical and millennial scale geomagnetic models, has persisted over the past 0.78 Ma.

2 of 23 understand paleofield behavior. Over the past two decades, significant progress has been achieved in obtaining high-quality data for the last 10 Myr (Cromwell et al., 2018;Johnson et al., 2008;Opdyke et al., 2015). These studies imposed rigorous selection criteria for obtaining acceptable data, requiring improved statistical parameters and determined from modern laboratory methods. In essence, these high-quality data compilations provide support for statistical investigations of non-GAD field contributions, detection of distinctions between normal and reverse polarity fields, and evaluation of possible hemispheric asymmetries in PSV and TAF estimates.
PSV is usually assessed by the angular dispersion of virtual geomagnetic poles (VGPs) assuming a GAD field (Merrill & McFadden, 2003). The latitudinal variation in VGP dispersion has been described by two main types of statistical PSV models: Model G (McFadden et al., 1988) and Giant Gaussian Process (GGP) models (Constable & Parker, 1988). These models were designed based on observable modern geomagnetic field properties, but use different statistical approaches. Both models consider the surface magnetic field in spherical harmonic form, where Gauss coefficients m l E g and m l E h (l is the degree and m is the order) define the field geometry. In the phenomenological Model G of McFadden et al. (1988), the overall VGP angular dispersion (S) is composed of two independent dynamo families, known as primary (or antisymmetric) and secondary (or symmetric) families, expressed by: The contribution relative to the primary family (S p ) comprises spherical harmonic terms with odd-numbered l − m (i.e., equatorially antisymmetric terms), which was assumed to vary linearly with latitude (λ). The secondary family (S s ) that contributes to the VGP dispersion is given by the equatorially symmetric terms for which l − m is even and is assumed to be approximately constant for all latitudes. In this respect, Model G can be described by a latitudinal variation curve fitted to the dispersion data, considering two shape parameters a (= S s ) and b (= S p /λ), respectively, which are related to the secondary and primary families (McElhinny & McFadden, 2000).
GGP-type models offer a different statistical perspective capable of predicting the distribution of geomagnetic field vectors in any geographical location, by describing the variances of Gauss coefficients. Initially defined by Constable and Parker (1988) who considered each Gauss coefficient as a normally distributed random variable, and the prescribed variances that reproduce the spatial power spectrum (Lowes, 1974) of the modern field. All non-dipole spherical harmonic terms have zero mean except for the axial quadrupole component. Prior to 2020, four GGP models: CP88 (Constable & Parker, 1988), QC96 (Quidelleur & Courtillot, 1996), CJ98 (Constable & Johnson, 1999), and TK03 (Tauxe & Kent, 2004) were presented, and are differentiated by the analytical methods and paleomagnetic 0-5 Ma database used (see Johnson & McFadden, 2015). Recently, the GGP models of BCE19 (Brandt et al., 2020), BB18 and BB18.Z3  have been used to evaluate the latitudinal behavior of PSV spanning the last 10 Myr using the 0-10 Ma database (hereafter PSV10) of Cromwell et al. (2018). The BCE19 model was constructed from analysis of paleodirectional data, taking into account the shape and dispersion of directional distributions at any latitude. It assumes that only the axial dipole term 0 3 of 23 contributions were reported for reverse polarity periods compared to normal polarity periods (Johnson et al., 2008;McElhinny & McFadden, 1997).
A further debate concerns potential hemispheric paleomagnetic field asymmetries Engbers et al., 2020). In particular, the latitudinal PSV structure for the past few million years suggests higher southern hemisphere VGP dispersion than northern hemisphere (Cromwell et al., , 2018Lawrence et al., 2009). Predictions from numerical geodynamo simulations indicate that lateral thermal variations at the core-mantle boundary (CMB) could influence geomagnetic field morphology and promote asymmetries (Aubert et al., 2013;Davies et al., 2008;Terra-Nova et al., 2019). However, our knowledge of long-lived PSV and TAF asymmetries is limited by the uneven temporal and geographic sampling of paleomagnetic records. In this context, improvements and expansion to the current 0-10 Ma database are important to better assess the TAF structure, and to constrain Earth-like geodynamo models.
Here we review and update the paleodirectional database derived from igneous rocks for the past 10 Myr using stringent selection criteria compared to the previous compilations. We then evaluate the latitudinal PSV and TAF structure for three age intervals corresponding to the Brunhes and Matuyama chrons, and the 0-10 Ma period. VGP dispersion estimates from Model G are compared with recent PSV studies and with the BCE19 and BB18 model family. We also quantitatively assess VGP dispersion patterns for the northern and southern hemispheres and extend investigations for the historical period from the COV-OBS geomagnetic field model (Gillet et al., 2015). We present new zonal TAF models for each age interval based on robust inclination anomaly estimates, and discuss the presence of non-dipole components relative to the GAD term.

Paleomagnetic Database
In order to investigate PSV behavior and TAF structure in the 0-10 Ma interval, we compiled an updated database of high-quality paleodirections from published studies. They include: the Magnetics Information Consortium (MagIC) database (https://www2.earthref.org/MagIC; Tauxe et al., 2016), academic search engines (e.g., Scopus at https://www.scopus.com/home.uri), and the PSV10 compilation (Cromwell et al., 2018). Only paleomagnetic data derived from volcanic rocks and thin dykes were accepted because these are considered most suitable for PSV and TAF analysis, and they provide instantaneous paleomagnetic field records in contrast to smoothed sedimentary records (Johnson & McFadden, 2015). We further employ the following selection criteria.
1. All paleomagnetic data must be within the age interval 0-10 Ma. 2. The characteristic remanent magnetization (ChRM) directions must have been acquired using modern demagnetization procedures and processing techniques. Additionally, site-level data must have an associated demagnetization code (the number of demagnetization site used for ChRM determination) equal to four or five (DC ≥ 4, McElhinny & McFadden, 2000). 3. Studies that did not provide site-mean directions and site coordinates were rejected. 4. Studies where the sampling region has been subjected to significant post-emplacement tectonism (e.g., tilting) were not included, based on information provided in the original publications. 5. All studies must comprise paleodirectional data from at least 10 sites (i.e., N ≥ 10) related to independent geomagnetic field records, for a given study location. 6. We require a minimum of five samples per site (n ≥ 5), with an estimated Fisher precision parameter (Fisher, 1953) of k ≥ 50 for each site-mean direction.
Detailed discussion of these six selection criteria can be found in Text S1 in Supporting Information S1.
The paleomagnetic database contains 2,543 directional sites from 80 paleomagnetic studies published between 1989 and 2020 that meet our selection criteria. Selected data sets associated with their respective geographic locations, average ages, DC values, and references are listed in Table S1. Moreover, our compilation supplies additional information for all site-level data (summarized in Table S2), including site coordinates, paleolocations (computed using the NNR-MORVEL 56 model (Argus et al., 2011) for plate motion corrections), Fisher site-mean directions, paleomagnetic statistical parameters, ages, and VGP coordinates. We 4 of 23 consider here site locations corrected for plate tectonic motions in the PSV and TAF statistical analyses (as adopted by Cromwell et al. [2018]), as discussed below.
Regarding the number of paleomagnetic sites and quality criteria, the new database supersedes the PSV10 data set (as discussed in Section 5.1), and offers a latitudinal coverage from 78°S to 78°N (Figure 1). There is an uneven spatial distribution of paleodirectional data between the northern and southern hemispheres; the latter consists of 23 data sets which corresponds to 29% of the collection.
Most of paleodirectional data are associated with high-level demagnetization procedures (DC = 5; 72% of the paleomagnetic sites; Table S2, Figure S1). The number of samples per site are mainly within the 5-10 range (∼89% of sites), whereas the Fisher precision parameter, k is concentrated between 50 and 200. Only 844 paleodirectional sites (∼33%) are associated with radiometric dating based on information from the original references. For other sites, average ages were inferred from stratigraphic and historic records discussed in selected studies. Paleomagnetic sites for last 10 Myr are mainly from the Brunhes (39%) and Matuyama (29%) chrons.

Methods
We examined the statistical behavior of the paleofield during three time periods similar to Johnson et al. (2008, hereafter J08): (a) the entire 0-10 Ma data sets, (b) Matuyama (0.78-2.58 Ma) reverse polarity data, and (c) Brunhes (0-0.78 Ma) normal polarity data. We aim to ascertain whether there are differences between normal and reverse polarity intervals. Data sets assigned to the Gauss and Gilbert chrons were not assessed due to insufficient data to perform a separate analysis for these time intervals.

Estimate of Paleosecular Variation
To assess PSV behavior, the between-site dispersion (S B ) for each data set was calculated according to (McElhinny & McFadden, 1997):  Table S1. 9 3 2 cos 2 cos 3 10 cos 6 cos 2 2 tan . 15 3 sin 2 3 cos sin 3 cos sin sin 2 2 We determined the zonal quadrupole ( G g g 2 2 0 1 0  / ) and octupole terms ( G g g 3 3 0 1 0  / ) by weighted least squares fittings between the observed ∆I anomalies and the predicted inclination anomaly (from Equations 6-8), where the weighting factors correspond to 95% uncertainties in the observed ∆I data. The best-fit TAF model was obtained based on the Levenberg-Marquardt method (see Section 3.1).

Transitional VGPs
Anomalous VGP data (characterized by large deviations from the mean pole) are generally excluded in PSV and TAF studies to ensure that data reflect stable field regimes. These anomalous VGPs (also referred to as outliers; Biggin et al., 2008) can be associated with geomagnetic excursions, polarity transitions, or spurious data caused by experimental measurement errors or secondary magnetizations (Johnson & Mc-Fadden, 2015). To evaluate the effects of transitional data on PSV and TAF estimates, we calculate the dispersion S B and ∆I after applying three approaches for each time period. The first applied no VGP cutoff and considered all directional data. The second used a fixed 45° cutoff angle, removing paleodirectional sites with VGPs that deviate >45° from the mean paleomagnetic pole. The third used the criterion of Vandamme (1994), where the cutoff angle (A) relative to mean pole is calculated as a function of dispersion S, by: An iterative method is employed according to the A value for a given data set; if VGP data are larger than the A estimate, they are excluded, S is recalculated, and the procedure is repeated until there are no larger deviations than the value determined from Equation 9.
Approximately 3% of paleodirectional sites from the entire 0-10 Ma data set are regarded as anomalous or transitional after applying a fixed 45° cutoff; this increased to 5% for the Vandamme (1994) criterion. The statistical distribution of paleomagnetic sites retained after applying a VGP cutoff for the three time periods is shown in Figure S2.

Analysis of Latitudinal Variation of VGP Dispersion
S B estimates for three temporal subsets: 0-10 Ma, Matuyama, and Brunhes age intervals are presented in Tables S3-S5, respectively. Tables S3-S5 also include, for each age group, S B results with lower and upper 95% confidence limits determined by the bootstrap method (Efron & Tibshirani, 1993), alongside values of average site latitude and longitude, mean direction, corresponding paleomagnetic pole, and paleomagnetic statistical parameters. Some data sets have substantially higher or lower S B estimates, so we apply an additional criterion proposed by Deenen et al. (2011) to assess whether the data adequately sample PSV, which is defined by an envelope A 95 limited by a range of values between A 95min and A 95max (see Text S2 in Supporting Information S1). This criterion was used to identify data with extremely different behaviors from that generally observed in paleomagnetic studies; these data were not considered in our PSV and TAF analyses. Furthermore, this criterion was used here because it was designed for datasets and models from targets with ages of a few million years and with thermoremanence acquisition similar to those found here, that is, data derived from rapidly cooling igneous rocks. In the 0-10 Ma period, we identify a maximum of 11 data sets that do not meet the criterion of Deenen et al. (2011; see Table S3). These possible PSV estimate biases may be related to transitional data that are over-represented in some data sets, especially when no VGP cutoff or a variable cutoff is employed (data sets # 5, 12, 22, 31, 76, and 80). In addition, data that are serially correlated or that are temporally underrepresented may contribute to anomalously low S B values in three paleomagnetic studies (data sets # 32, 39, and 48). Another relevant factor might be undetected regional tectonic effects that can enhance S B estimates.  Table 1. Furthermore, the S B data distribution is compared with three recent GGP models: BCE19, BB18, and BB18.Z3. The method employed for predicted dispersions is described in Text S3 in Supporting Information S1.

of 23
S B estimates for the 0-10 Ma interval (80 data sets, Table S3), defined using three VGP cutoff approaches, are displayed in an interhemispheric representation in Figure 2a. Overall, an increasing dispersion S B trend can be observed as a function of latitude for both hemispheres. Differences in PSV estimates are identified when the fixed (45°) cutoff and Vandamme (1994) criterion exclude transitional data. Model G fits most S B data reasonably well, regardless of the VGP cutoff employed (Figures 2b-2d). S B data with no VGP cutoff (69 data sets, Figure 2b) produces a best-fit Model G curve with parameters  Vandamme (1994) cutoff. In (d), the yellow line denotes the Model G curve fitted to the PSV10 data compilation (Cromwell et al., 2018). Also shown are VGP dispersions predicted for three GGP models: BCE19.GAD model of Brandt et al. (2020); light green dash-dot line; BB18 and BB18.Z3 models  purple and pink lines, respectively.
OLIVEIRA ET AL.

10.1029/2021GC010063
8 of 23 For the Matuyama reversed polarity chron (19 data sets, Table S4), a pattern of increased S B values is observed as a function of latitude for both hemispheres (Figure 3a). Nevertheless, there is a data gap at low (0-25°N) and high (>65°N) northern latitudes, and data coverage is even more restricted in the southern hemisphere. When a VGP cutoff (variable or fixed) is employed, S B estimates are especially different for data sets from Antarctica at 78°S (Lawrence et al., 2009) and North America at latitudes 60°N (Coe et al., 2000), and 46°N (Lhuillier et al., 2017). Model G for data sets with no cutoff (13 data sets, Figure 3b) yields , but these estimates are not statistically different considering the 95% confidence intervals. When the Vandamme (1994) criterion (16 data sets, Figure 3d) Figure 3d, predicted VGP dispersion values for BB18-family models are lower at low and high latitudes compared to the Model G fit, whereas the BCE19 model has lower dispersion along the entire latitudinal range relative to the Model G. BB18 and BB18.Z3 models fit the S B data better as a function of latitude compared to the BCE19 model.
The interhemispheric distribution of Brunhes normal polarity data (46 data sets, Table S5) is illustrated in Figure 4a. The northern hemisphere has good data coverage to 70°N, while the southern hemisphere presents a sparse data coverage, especially in mid-to high-latitudes (40°-60°S). Most studies have small differences in VGP dispersion when fixed or variable cutoffs are used. The best-fit curve for S B estimates when no , respectively. All Model G parameters are similar regardless of cutoff employed. Regarding the GGP models, predicted dispersions for BB18 and BB18.Z3 models are compatible with the observed S B data in a broad range of latitudes, while the BCE19 model is compatible with some data at low to mid-latitudes. When GGP models are compared to the Model G curve, the expected dispersions for BB18-family models are slightly higher for latitudes from 20° to 70°, whereas the BCE19 Model provides lower dispersion estimates for all latitudes.

Analysis of Latitudinal Structure of Inclination Anomaly
In general, ∆I data for the 0-10 Ma interval (Table S3 and Figures 5a-5c) are not statistically distinguishable from the GAD field model (dashed line) at mid-to high-latitudes for both hemispheres, which contrasts with the large negative inclination anomalies at low latitudes (0-30°N and S). The best-fit field model for 0-10 Ma data without VGP cutoff (Figure 5a) Figure 5c). Some data differ significantly from the predicted zonal field model within the 95%  (Table 1), for instance, the high negative or positive ∆I values in Norway (78°N;  and Saint Helena Island at 18°S, respectively (Engbers et al., 2020).
ΔI estimates for the Matuyama chron (Table S4 and Figures 5d-5f) are more limited, with negative and positive values in mid-to high-northern latitudes, while the southern hemisphere has positive (<10°) anomalies; exceptions are low negative ΔI values at low latitudes from Ecuador at 0.5°S (Opdyke et al., 2006) and French Polynesia at 18°S (Yamamoto et al., 2002). The distribution of ∆I estimates with no VGP cutoff ( Figure 5d) Tables S3-S5). u Upper 95% confidence limit. l Lower 95% confidence limit.

of 23
Israel (Behar et al., 2019), and Saint Helena Island (Engbers et al., 2020). In terms of temporal distribution, there was a 5% increase in data coverage for the age interval older than 5 Ma.
Moreover, our stricter selection criteria incorporate a minimum number of paleodirectional sites per study (N ≥ 10) and number of samples per site n ≥ 5, which differs from PSV10 database (with n ≥ 4). These criteria exclude 15 of the original 81 data sets from the PSV10 compilation. Through inspection of the S B distribution for individual PSV10 data sets by applying the Vandamme (1994) criterion (hereafter PSV10 vcut ; Table S6), higher S B estimates are obtained for some data than are expected from Model G (suggested by Doubrovine et al. [2019]) at mid northern latitudes ( Figure S3). Most of these overestimates may be associated with data undersampling, with N < 10, in accordance with the analyses of J08. Large N and n values are recommended to reduce bias in PSV and TAF estimates (Biggin et al., 2008;Johnson & McFadden, 2015). Cromwell et al. (2018) demonstrated that intrasite directional variance is reduced by ∼7% for n = 5 compared to n = 4. These assessments indicate the high-quality of data sets in the present study that can be used to examine paleofield latitudinal structure.
By comparing the Model G curve fitted to the PSV10 vcut subset (yellow line in Figure 2d) and the best-fit curve for the new 0-10 Ma data sets using a variable cutoff, a lower latitudinal dependence of the VGP dispersion curve is observed here (Figure 2d) Table S7). The small differences in the VGP dispersion curves are probably associated with the method used to estimate PSV. We determined the VGP dispersion for each data set individually, in contrast to the mean S B calculated for 10° latitude bands by Cromwell et al. (2018). We choose a statistical approach that allows us to investigate PSV and TAF at the study level, as has been widely employed for studies over the Phanerozoic and Archean time scales (e.g., de Oliveira et al., 2018;Doubrovine et al., 2019;Franco et al., 2019;Veikkolainen & Pesonen, 2014). In addition, it is possible to analyze differences between S B estimates that may be caused by the paleomagnetic study location where binning is not viable.

Evidence for Hemispheric Asymmetry of Paleosecular Variation
In order to identify possible VGP dispersions pattern differences between the northern and southern hemispheres, we tested PSV equatorial symmetry by evaluating the latitudinal behavior of the Model G curves fitted to reliable S B data for the 0-10 Ma interval ( Figure 6a) and Brunhes chron (Figure 6b) for the northern and southern hemispheres (closed and open symbols, respectively). The Matuyama subset was not considered separately due to poor S B coverage data over a broad latitudinal range. For the past 10 Myr, the overall patterns of the VGP dispersion curves for the northern and southern hemispheres are similar, regardless of the VGP cutoff filter used (see also Figures S4 and S5 in Supporting Information S1). The best-fit curve for S B data using the Vandamme (1994) criterion (Figure 6a) Figure S4.
For the Brunhes subsets, the best-fit Model G curves suggest northern-southern hemisphere asymmetry. When the variable cutoff is employed (Figure 6b) ). However, the differences between the Model G parameters are both within the 95% confidence interval (see also Figure S5 for Brunhes subsets without a VGP cutoff and using a constant 45° cutoff).
Based on comparative assessments among VGP dispersion curves for the northern and southern hemispheres for 0-10 Ma and Brunhes data, the latter provides evidence of a hemispheric PSV asymmetry, especially when the Vandamme (1994) criterion is applied. These findings support the hypothesis of higher (lower) VGP dispersion in the southern (northern) hemispheres for at least the past 0.78 Ma. In addition, our quantitative assessment differs from previous observations (Cromwell et al., , 2018. Recent geomagnetic field models that evaluate PSV indices (P i , a measure of regional field variability) for present, historical (Panovska & Constable, 2017), multimillennial (0-10 ka; Constable et al. [2016], and 0-100 ka 13 of 23 [Panovska, Constable, & Brown, 2018;Panovska, Constable, & Korte, 2018]) timescales yield higher P i values in the southern hemisphere than the northern hemisphere. These differences are associated with low field intensities, which suggest that the equatorial PSV asymmetry is a long-period feature of the geomagnetic field.
Hemispheric field structure differences have been attributed to lateral CMB heat flux heterogeneity, based on numerical dynamo simulations (e.g., Aubert et al., 2013;Christensen & Olson, 2003;Terra-Nova et al., 2019). Davies et al. (2008) reported a higher latitudinal variation of synthetic VGP dispersion in the southern hemisphere than the north, from a dynamo model run under thermal heterogeneous CMB conditions. However, VGP dispersion estimates observed from the geodynamo models were significantly lower when compared to PSV models (e.g., Johnson et al., 2008;McElhinny & McFadden, 1997).  Vandamme (1994) cutoff. Light (dark) red lines represent the Model G curve fitted to northern (southern) hemisphere data, with its 95% confidence intervals (dashed lines).

Comparison With Historical Equivalent VGP Dispersion
It is well known that snapshots of geomagnetic field models are not equivalent to TAF models (Hulot & Gallet, 1996;Kono & Tanaka, 1995;Merrill et al., 1998). Even so, McFadden et al. (1988) showed that VGP dispersion curve for the present geomagnetic field (computed from the IGRF65 model) is similar to the 0-5 Ma paleomagnetic field. Nevertheless, their results only considered the average VGP dispersion over both southern and northern hemispheres. Furthermore, Hulot and Gallet (1996) showed that historical VGP dispersion is highly time-dependent; for earlier historical periods, the similarity observed by McFadden et al. (1988) no longer holds. They also showed that the Gauss coefficients of degree and order (1,1), (2,1), (3,1), and (4,1) are the main terms that influence the shape of the VGP dispersion curve. Here, we further investigate the time dependence of the equatorial asymmetry in VGP scatter curves.
We evaluate the equivalent VGP dispersion for the historical period using the time-dependent geomagnetic field model COV-OBS based on stochastics methods (Gillet et al., 2015). Distributions of dispersion data were computed for each 5° latitude band, with 5° longitude spacing. Equivalent VGP dispersion curves from 1840 to 2015 in five-year intervals are shown in Figure 7a, considering VGP dispersion estimates between northern and southern hemispheres. As expected from previous observations (e.g., McFadden et al., 1988;Hulot & Gallet, 1996), the VGP dispersion tends to increase from the equator to the pole. However, for earlier epochs (before 1880) a decreasing dispersion from 60° latitude is observed. By comparing the Model G best-fit curves from this study with those for previous PSV studies for the 0-5 Ma (Opdyke et al., 2015) and 0-10 Ma (Cromwell et al., 2018) intervals, we detect an incompatibility between paleomagnetic VGP dispersions and equivalent VGP dispersions from a snapshot of historical geomagnetic fields, in contrast to the observations of McFadden et al. (1988). Nevertheless, latitudinal geomagnetic field variation for older periods (yellow curves in Figure 7a) and the 0-10 Ma paleomagnetic field (from this study using the Vandamme [1994] criterion) capture a similar pattern within the 0°-40° latitudinal range.
When historical equivalent VGP dispersions are examined separately for each hemisphere, there is a pronounced equatorial PSV asymmetry especially for recent intervals of geomagnetic field behavior (Figure 7b). Neither the northern nor the southern hemisphere equivalent dispersions have any striking similarity compared to VGP dispersion curves from PSV data (Figure 7a). Southern hemisphere scatter is higher with a bump around 20°, which probably relates to the historical latitudinal position of the South Atlantic Anomaly (e.g., Amit et al., 2021;Hartmann & Pacca, 2009;Terra-Nova et al., 2017). In the northern hemisphere, there is a VGP dispersion decrease with latitude, with a bump at around 50-60°N. A similar feature was found by McFadden et al. (1988) for the IGRF65 model. These results highlight that the historical equivalent VGP dispersions produce higher equatorial asymmetry than that discussed in PSV studies. This could relate to the short time span of the geomagnetic field model, but should not be over-interpreted. Curiously, equivalent southern hemisphere dispersion better describes the paleomagnetic data than the equivalent northern hemisphere dispersion.
We also investigate the spherical harmonics responsible for the observed time evolution of equivalent VGP dispersion curves. First, we calculate dispersion distributions in terms of specific Gauss coefficients from the COV-OBS model for 1965 (Figure 7a). Following Hulot and Gallet (1996), we expand the investigation of how non-zonal harmonic terms of low degree (l = 1, 2, 3, 4) and order m = 1 alter the shape of the VGP dispersion curve. The nullifying of certain non-zonal harmonics changes significantly the latitudinal behavior of VGP dispersions. For instance, when the four non-zonal terms are reduced to zero, the dispersion distributions tend to decrease at high latitudes (>60°; brown curve in Figure 7c).
To quantify hemispherical differences in the time evolution of historical equivalent VGP dispersions, we define the interhemispheric variance (∆S) expressed by: where S NH and S SH are dispersion estimates for the northern and southern hemispheres, respectively. The higher ∆S, the more asymmetric are the northern and southern equivalent VGP dispersion curves. By evaluating the average energy of non-zonal harmonics as a function of ∆S estimates (Figure 7d), ∆S is observed to increase progressively with dipole field decay (including the axial 0 1 E g term) and increased non-dipole field contributions over historical time. When establishing a linear regression to the data, a best-fit was found for OLIVEIRA ET AL.

10.1029/2021GC010063
15 of 23 the (2,1) term. Further investigation of ratios of non-dipolar fields at the CMB might clarify the source of this asymmetry but is beyond the scope of this study.

Evaluations of Non-Dipole Zonal Terms Over the Last 10 Myr
The new zonal TAF models fitted to ∆I data indicate a small non-dipole field contribution. Specifically, for the 0-10 Ma period, the best estimate for the axial quadrupole 0 2 E g term is 3.2%-3.7% of 0 1 E g , with a smaller axial octupole contribution g 3 0 ranging from 1.2% to 1.8% of 0 1 E g (Table 1). Our results suggest that small departures from the GAD model have persisted over long timescales for the 0-100 ka and 0-5 Ma periods. From their GGF100k model (Global Geomagnetic Field for the 0-100 ka interval), Panovska, Constable, and Korte (2018) revealed the presence of the zonal terms G2 = 4.2% and G3 = 2.5%. For a statistical TAF model for the last 5 Myr, J08 suggests estimates of G2 = 3% and G3 = 3% from high-quality paleomagnetic data.
The BB18.Z3 model  suggests values of G2 = 2.9% and G3 = −1.3% for the past 10 Myr. Thus, small non-GAD contributions seem to be a common feature of paleomagnetic field models.

Examining Paleomagnetic Inclinations Relative to the GAD Approximation
The question of the paleomagnetic field geometry found to be predominantly dipolar over the past million years has been debated (Johnson & McFadden, 2015;McElhinny, 2004). For the three time periods investigated in this study, our TAF models reveal low values of the G2 and G3 terms (Table 1). Following the previous TAF studies (e.g., Opdyke & Henry, 1969;Opdyke et al., 2015;Schneider & Kent, 1990), we also examine the latitudinal distribution of mean inclination data (applying the Vandamme [1994] cutoff) to test compliance with the expected GAD inclination (determined by Equation 7), as shown in Figure 8.
Most of mean inclination data for the 0-10 Ma interval agree with the latitudinal variation curve from a GAD model (black dashed line in Figure 8a), and are not statistically different at the 95% confidence level relative to predicted GAD inclination. Similarly, Matuyama reverse and Brunhes normal polarity data (Figures 8b and 8c,respectively) show a close correspondence with the GAD inclination curve, considering the 95% uncertainties of the mean inclinations. Nevertheless, the Matuyama estimates at northern latitudes produce large deviations from the GAD prediction compared to the Brunhes chron.
For each age group, we calculate the weighted root mean square (RMS) deviation of the difference between the mean paleomagnetic inclination and the predicted GAD inclination. The weights correspond to the 95% uncertainties of the mean paleomagnetic direction. Our analyses show that the three temporal data sets are best-fit for the corresponding TAF models presented in this study (summarized in Table 1), with slightly lower RMS values than for GAD model. The inclination curves for TAF models (red line) are shown in Figure 8. These results suggest that the inclination estimates are compatible with the GAD field, however, the paleofield models with small G2 and G3 contributions provide a statistically acceptable fit to the observed inclination data. Furthermore, the non-dipole field components do not yield significant changes for the last 10 Myr.

Distinct PSV and TAF Patterns: Brunhes Normal and Matuyama Reverse Polarity Data Sets
Statistical analysis of PSV and TAF latitudinal behavior reveals differences between the Brunhes normal and Matuyama reverse polarity data, which partially supports the observations of J08. For PSV analysis, we examine Model G parameter estimates fitted to S B data for both age groups, which differs from that presented by J08. Best-fit curves for Matuyama subsets yield higher S B estimates at low latitudes compared to the Brunhes subsets, especially when a 45° cutoff angle is used to exclude outlier data (Figures 3 and 4). When compared to J08, our Matuyama estimates (using a fixed 45° cutoff) produce less latitudinal variation in the VGP dispersion curve (Figure 3c). However, the difference is statistically significant with respect to Model G parameter b (Table 1 and Table S7). Differences between the two VGP dispersion curves are probably  Vandamme (1994) criterion. The uncertainties of the inclination data correspond to the 95% confidence cone (α 95 ) of the mean direction. Black curves represent the expected GAD inclination. Red curves are the inclination predicted for TAF models proposed in this study (see Table 1). Vand. = Vandamme (1994); RMS = root mean square.
OLIVEIRA ET AL.

10.1029/2021GC010063
17 of 23 related to inclusion of new mid-to high-latitude (±50-65°) records that were not included in J08 and that may have influenced the Model G fit. For the Brunhes chron (Figure 4c), our predicted S B values for Model G are lower for J08, but Model G parameters in both studies are not statistically distinguishable at the 95% confidence level. The slight observed differences may be associated with the data sets distribution and size in both studies (Table 1 and Table S7).
According to some studies (Cromwell et al., 2018;Valet & Herrero-Bervera, 2011), high (low) S B estimates for reverse (normal) polarity fields might be attributed to lower (higher) average paleointensity during the Matuyama (Brunhes) chron. Ahn and Yamamoto (2019) analyzed 0-5 Ma absolute paleointensity data from the PINT database (Biggin et al., 2009) and concluded that the Matuyama had a lower dipole field strength compared to the Brunhes period. Recent studies Meduri et al., 2021) have also pointed out an inverse relationship between the a parameter from Model G with the degree of dipole dominance, inferred from numerical geodynamo models. Here, the shape parameter a is higher for Matuyama data sets than for the Brunhes, whichever VGP cutoff is applied (Table 1). These findings suggest that the TAF was less dipolar during the Matuyama than the Brunhes chron. However, observed differences between Model G parameters for these two age groups are not statistically significant at the 95% confidence level.
Zonal TAF model fits for the Matuyama ∆I data result in higher quadrupole (G2 of 3.1%-3.6%) and octupole (G3 of 2.0%-2.9%) contributions compared to Brunhes estimates (G2 = 2.0-3.1% and G3 = 1.3-1.8%). The exception is the low value obtained for the G2 term by applying a fixed 45° cutoff (see Table 1). J08 reported estimates of G2 = 4% and G3 = 5% and G2 = 2% and G3 = 1% for the Matuyama and Brunhes chrons, respectively. The main differences in TAF estimates between the studies can be assigned to two factors: (a) incorporation of new data produced after J08 and (b) the formalism used for ∆I calculations (here ∆I estimates were calculated for each paleomagnetic study instead of latitude band averages, as adopted by J08). Furthermore, high estimates of non-dipole components detected for the Matuyama epoch could be linked to the lower dipolar paleofield dominance in this period than during the Brunhes chron.
Thermal and compositional heterogeneities at the CMB have been suggested as a possible explanation for differences between normal and reverse polarity fields, as debated in the literature (Johnson & McFadden, 2015). McElhinny and McFadden (1997) suggested that distinctions between these two polarity states, ascertained from the latitudinal behavior of 0-5 Ma VGP dispersions, could rather be caused by incompletely removed viscous overprints for reverse polarity data (with higher dispersion compared to normal polarity data). However, considering the high-quality data assessed in this study, we assume that this effect is minimized. Thus, polarity asymmetries in the PSV and TAF models for Matuyama and Brunhes subsets appear to be an empirical feature of the paleomagnetic field, which corroborate earlier observations (e.g., Cromwell et al., 2018;Johnson et al., 2008).

Conclusions
1. An upgraded database of high-quality igneous paleodirections produced using rigorous selection criteria, provides robust PSV and TAF assessments for the Brunhes and Matuyama chrons over the 0-10 Ma period. The new data compilation contains 2,543 paleomagnetic sites covering the past 10 Myr, which improves the temporal and spatial distributions of paleodirectional data relative to the previous PSV10 database (Cromwell et al., 2018). 2. VGP dispersion patterns from the 0-10 Ma data set, and from Matuyama and Brunhes subsets, are well-fitted by an adapted Model G (McFadden et al., 1988) that differs slightly from the dispersions predicted by the GGP models BCE19 (Brandt et al., 2020) and BB18 and BB18.Z3 . It is noteworthy that the BCE19 model and BB18-family models were designed to fit, respectively, paleodirectional and VGP dispersion estimates in the PSV10 compilation. 3. The VGP dispersion curve fitted for new 0-10 Ma data sets using Model G produces lower latitudinal dependence of S B than the PSV10 data set. Analyses of the S B data as a function of latitude suggest that differences in the shape of Model G curves are associated with data selection and calculation methods for PSV estimates. Our results differ by including more stringent quality criteria (with N ≥ 10, n ≥ 5), and assess latitudinal PSV variation by locality using the Deenen et al. (2011) criterion rather than mean S B values by latitude bands.
4. New zonal TAF model for the 0-10 Ma interval indicate the presence of small non-dipole field contributions, with an axial quadrupole component G2 = 3.2-3.7% and a smaller axial octupole component G3 = 1.2-1.8%. 5. Assessments of the latitudinal PSV behavior do not provide clear evidence of a north-south hemispheric asymmetry in VGP dispersion for 0-10 Ma data sets. Brunhes data have differences between hemispheres that are characterized by a stronger latitudinal S B dependence in the southern relative to the northern hemisphere, especially when the Vandamme (1994) cutoff is used. This finding suggest that equatorial asymmetry of geomagnetic secular variation has persisted over the last 780 ka. 6. Statistical simulations from the COV-OBS model for the 1965 epoch indicate that the reduction or combined effect of non-zonal harmonic terms (1,1), (2,1), (3,1), and (4,1) to zero implies significant modifications in the shape of the equivalent VGP dispersion curve. Furthermore, investigations of the historical evolution of dispersion estimates for the northern and southern hemispheres indicate an equatorial asymmetry of VGP scatter that gradually increases with time, which can be associated with dipole field decay and increased non-dipole field contributions. 7. We report differences between Brunhes normal and Matuyama reverse polarity data in both PSV and the TAF analysis. The Matuyama subset produces higher S B patterns at low latitudes than the Brunhes chron. In terms of non-GAD TAF structure, zonal quadrupole and octupole contributions are larger for the Matuyama chron (with G2 = 3.1-3.6% and G3 = 2.0-2.9%) compared to the Brunhes (G2 = 2.0-3.1% and G3 = 1.3-1.8%), which support observations of lower and higher average paleointensities for the respective periods. 8. Finally, additional high-quality data are essential to enhance temporal and geographic sampling of the 0-10 Ma database (especially data coverage older than 1 Ma and new southern hemisphere records), and to better understand long-period geomagnetic field asymmetries.