Characteristics of the Matuyama‐Brunhes Magnetic Field Reversal Based on a Global Data Compilation

Magnetic field reversals are irregular events in Earth's history when the geomagnetic field changes its polarity. Reversals are recorded by spot and continuous remanent magnetization data collected from lava flows and marine sediments, respectively. The latest field reversal, the Matuyama‐Brunhes reversal (MBR), is better covered by paleomagnetic data than prior field reversals, hence providing an opportunity to understand the physical mechanisms. Despite the quantity of data, a full understanding of the MBR is still lacking. The evolution of the MBR in time and space is explored in this work by compiling a global set of paleomagnetic data, both from sediments and volcanic rocks, which encompass the period 900–700 ka. After careful evaluation of data and dating quality, regional and global stacks of virtual axial dipole moment (VADM), virtual geomagnetic pole (VGP), and paleosecular variation index (Pi) are constructed from the sediment records using bootstrap resampling. Individual VADMs and VGPs calculated from lavas are compared to these stacks. Four phases of full‐vector field instability are observed in these stacks over the period 800–770 ka. The first three phases, observed at 800–785 ka, reflect a rapid weakening of the field coupled with low VGP latitude, after which the field returned to the reverse polarity of the Matuyama chron. The fourth phase, lasting from 780 to 770 ka, is when the field reversal process completed, such that the field entered the Brunhes normal polarity state. These findings point to a complex reversal process lasting ∼30 Kyr, with the reversal ending at ∼770 ka.

The geometry of the MBR was investigated by tracking VGPs, known as VGP paths, calculated from lavas and sediment data of transitional polarity. Different assumptions were proposed from these paths, which include an axisymmetric transitional field (Clement & Kent, 1984;Hoffman & Fuller, 1978;Valet et al., 1988), longitudinally preferred VGP paths and VGP clusters (Clement, 1991;Laj, 1991;Valet et al., 1989), and a long-lived transitional field state (Hoffman, 2000). The validity of these assumptions were, however, questioned (Langereis et al., 1992;Prévot & Camps, 1993;Valet & Fournier, 2016). Also, Hartl and Tauxe (1996) have found a period of decrease in paleointensity ∼15 Kyr (at ∼795 ka) before the main MBR, known as precursor, from several sediment cores. This precursory event was found to be accompanied by significant directional change, from a series of volcanic records span from 180 to 0.78 Myr ago (Valet et al., 2012). The field precursor was also seen in records from the North Atlantic (Channell et al., 2004), North Pacific (Korff et al., 2016), Antarctica (Macrì et al., 2010), Mediterranean (Just et al., 2019;Sagnotti et al., 2014), and Indian Ocean (Valet et al., 2014). Following the polarity transition, Valet et al. (2012) observed another unstable geomagnetic phase with low-latitude VGPs, which they called a rebound. Unlike the precursor, the post-reversal rebound state was only documented in a few records (Singer et al., 2019). According to Valet et al. (2012), the durations of the precursory and rebound phases were estimated at 2.5 Kyr. The global magnetic field during MBR and its evolution in space and time were reconstructed using spherical harmonic (SH) models (Ingham & Turner, 2008;Leonhardt & Fabian, 2007;Shao et al., 1999). The iterative Bayesian inversion model of Leonhardt and Fabian (2007), called IMMAB4, indicates the existence of equatorial radial magnetic flux patches that move poleward during the MBR at the core-mantle boundary. Similarly, the Ingham and Turner (2008) SH model predicts low-latitude flux patches at ∼795 ka, which can be linked to the MBR precursor (Hartl & Tauxe, 1996). However, because the input data in these models are limited and may have certain issues with resolution and timescale reliability (see, e.g., Valet & Fournier, 2016), further global analyses on high quality data are required.
In order to precisely capture the field reversal behavior, lava samples should have precise age data (e.g., K-Ar and 40 Ar/ 39 Ar) with uncertainty values lower than duration of the field reversal itself. Also, well-dated sediment cores with high sedimentation rate (SR) and temporal resolution better than the reversal duration are needed. Ideally, the age uncertainty of lavas should not exceed few hundreds of years, and the temporal resolution of sediments should be lower than this limit to properly define the reversal. For sediments, in order to have a record with 200 year resolution, the SR should be greater than 10 cm/Kyr (Roberts & Winklhofer, 2004), but some additional factors, such as presumption of constant and continuous SR and type of sample discrete or U-channels, (e.g., Roberts, 2006;Sagnotti et al., 2016) are fundamental as well. A second issue is the ability of a paleomagnetic sample to reliably retain paleomagnetic field components (declination, inclination, and intensity) during the reversal when field intensity reaches minimal values (Singer et al., 2019). In this context, lava samples are preferable because they acquire the paleomagnetic field information by thermal remanent magnetization (TRM), for which the theoretical background is well established (see Néel, 1955). Sediment-derived detrital, depositional, and post-depositional remanent magnetization (PDRM), on the other hand, has some limitations in terms of, among others, signal smoothing (Lund & Keigwin, 1994), potential inclination shallowing (Arason & Levi, 1990;Deamer & Kodama, 1990), and lock-in depth (Channell et al., 2004;Lund & Keigwin, 1994;Sagnotti et al., 2005), which can significantly impair their dependability. Moreover, the magnetization acquisition processes in sediments decrease their efficiency in recording the field during transitions (Valet & Fournier, 2016).
However, (P)DRM data are continuous while lavas only give scattered TRM data, therefore sediment cores are more efficient for investigating the temporal evolution of the field reversal. The best solution is an integration of different inputs, as done by Singer et al. (2019).
In this study, we evaluate the MBR by combining the time series of 38 sediment cores with independent age models, selected from a total of 68 records distributed over 16 regions of the globe. Half of these cores have mid-to high SR in the range of 5-16 cm/Kyr, so we expect to have a reasonable resolution. The data, their selection and treatment are described in Section 2, including regional consistency checks. Some global characteristics of the magnetic field during the MBR are inferred from regional and global stacks of VGP and virtual axial dipole moments (VADM) of these data in Section 3. Volcanic data are included in this analysis. We estimate the time and duration of the last field reversal, using a global stack of the paleosecular variation index (Panovska & Constable, 2017) in Section 4. The relation between duration and local site latitude is also discussed there.

Data Compilation and Treatment
Paleomagnetic data from lavas and marine cores have been compiled and analyzed in this study. We included all data, accessible to us, of the age period 900-700 ka. The compilation process was completed by carefully assessing the original articles, and by using open-access global databases, such as the Absolute Paleointensity database (PINT; Bono et al., 2022); PANGAEA (https://www.pangaea.de/) and Magnetics Information Consortium (MagIC; https://www2.earthref.org/MagIC). Also, we consider the compilations SEDPI06 (Tauxe & Yamazaki, 2007 and PADM2M (Ziegler et al., 2011), where the relative paleointensity (RPI) data of sediment records were provided.
The next sections summarize how we compiled and examined age and paleomagnetic data of lavas and sediments records. More details about the compilation and evaluations processes are available in the Texts S1 and S2 in Supporting Information S1, respectively.

Lava Data
We compiled paleomagnetic data from lavas with only absolute radioisotope, 40 Ar/ 39 Ar and K-Ar, age data. Lavas that have been dated relatively using magnetic polarity and stratigraphy are not included in this study. Considering this dating constraint, a number of 108 lava sites have been compiled, and 107 site-mean directions and 42 site-mean absolute paleointensities (PI) were calculated from them. Table 1 contains a list of these sites in terms of geographic location, ages, paleomagnetic directions, and PIs. The distribution of compiled lavas is illustrated in Figure 1a, which reveals that lavas were reported from 11 regions (Tahiti, Hawaii, Guadalupe, Chile, Canary Islands, Azores Islands, Iceland, NW-USA, Mexico, Germany, and Antarctica). Six of the 11 lava regions are located in the western hemisphere, three are located along the Atlantic (longitude ∼335-343°E), and only two are located in the eastern hemisphere. Of the 108 sites, 29 have K-Ar ages while the rest were 40 Ar/ 39 Ar dated.
Our compilation is based on the latest version of the PINT database (http://www.pintdb.org/), which contains 200 PI data points for the period 700-900 ka, 198 of them from lava flows. In terms of age, these 198 data points were dated as follows: 91 have 40 Ar/ 39 Ar ages, 26 with K-Ar, and the remainder have been dated by relative means (magnetic polarity and stratigraphy). By considering only absolute ages, only 117 lava sites remain. We found that certain sites in PINT include absolute 40 Ar/ 39 Ar age data, but the original papers where these data were cited had no such absolute ages. For instance, we compiled 10 40 Ar/ 39 Ar dated lava sites from Chile, whereas PINT reported 15 sites. After examining them, we noticed that the five extra sites in PINT had not been assigned 40 Ar/ 39 Ar age data. Moreover, we consider data updates; when a lava flow is re-sampled for new age estimates or PI measurements, we include only the most recent updates if their quality is superior. An example for that is PI data reported from La Guadeloupe Island. Carlut et al. (2000) obtained PI data from three La Guadeloupe Island lava sites, which were then resampled for new 40 Ar/ 39 Ar age estimations (M. C. Brown et al., 2013) and PI measurements (M. C. Brown et al., 2009). We take these updates into account in terms of ages (M. C. Brown et al., 2013) and average the PI data of the two studies (M. C. Brown et al., 2009;Carlut et al., 2000). Also, we note that some of the PI data from previous works were rejected by us since further analysis revealed that their ages were not in logical sequence with their stratigraphic context; see Text S1 and Figure S1 in Supporting Information S1 for more information. This explains why the PINT data set is larger than our final data set.
It is known that conventional K-Ar ages are not as accurate as 40   Quidelleur and Valet (1996) Quidelleur and Valet (1996) Quidelleur and Valet (1996) Quidelleur and Valet (    Ar/ 39 Ar, Argon/Argon; K-Ar, potassium/Argon; VGP lat, virtual geomagnetic pole latitude; A95, radius of the 95% confidence circle about the calculated mean pole. We estimated it as the square root of dp × dm (see Khramov (1987)) which are the semi axes of the ellipse at 95% confidence level (calculate from equations 7.8 and 7.9 in Butler (1992)), VADM, virtual axial dipole moment; SD, standard deviation; Am 2 , Ampere/meter 2 ; PI, paleointensity; na, not applied; Q PI , quality of paleointensity data.
All 40 Ar/ 39 Ar age errors are provided with 2σ, but K-Ar age errors are provided with 1σ. conv and CG refer to conventional and Cassignol-Gillot (Gillot & Cornette, 1986;Quidelleur et al., 2001) techniques that were use during K-Ar dating, respectively. a Refer to VADM calculated from RPI data (Azores) after scaling them to absolute values (See Section 2.2.5. We use the Q PI values as listed in the PINT database (Bono et al., 2022). The Q PI values of the four Guadeloupean sites (G03, G02, G01, and GD01) and the Chilean site (QTW11-20), which are noted with an asterisk *, were newly assigned for this study. Also, we assigned Q PI values for the Azores sites (designated with a **)). McDougall, 1967). In addition to the traditional K-Ar approach, the K-Ar Cassignol-Gillot technique (Gillot & Cornette, 1986;Quidelleur et al., 2001), which applied to separated groundmass, is mentioned here. This technique is suitable for dating Quaternary volcanic rocks (Guillou et al., 1996;Quidelleur et al., 2001Quidelleur et al., , 2003Samper et al., 2007), with a detection limit being of 10% of the radiogenic Ar content (Quidelleur et al., 2001). For K-Ar age data assembled here (Table 1), four studies (Abdel-Monem et al., 1972;Armstrong, 1978;Nelson & González-Caver, 1992;Wilch, 1991) provided K-Ar ages using the conventional approach, whereas the remaining K-Ar ages were determined using the refined Cassignol-Gillot technique. However, we accept all K-Ar age data in this study because the available database is already sparse. Reported age uncertainties range from 3 to 660 ka, with 22 sites having age errors ≤5 ka, 30 in the range 5-10 ka, and 56 sites have errors greater than 10 ka. Age errors are larger in sites with K-Ar ages than in 40 Ar/ 39 Ar dated sites. This situation prevents from obtaining crucial information regarding the fast changes in the magnetic field during the reversal. Nevertheless, the age data are carefully checked, and any updates are taken into account. The 40 Ar/ 39 Ar method requires standards that are used to calibrate the quantity of 40 K that can be attributed to a specific amount of 40 Ar in an analysis. Five standards are typically used in the 40 Ar/ 39 Ar geochronology: FCs, ACs, Mklhb-I, GA-1550, and TCs (for details of these standards see, e.g., Renne et al. (1998)). In this study, all 40 Ar/ 39 Ar ages have been recalibrated using the equivalent ACs (1.1851 ± 0.0004; Schaen et al. (2020)) or FCs (28.201 ± 0.046; Kuiper et al. (2008)) standards using the Min et al. (2000) decay constant (λ = 5.463 ± 0.054 × 10 −10 /year). This is another source of difference between our compilation and PINT, as the latter did not perform any recalibration on the 40 Ar/ 39 Ar ages. The original ages are listed alongside the calibrated ages in Table 1. The calibration is done using the following equation given by Renne et al. (1998): Paleomagnetic directions amount to nearly 2.5 times the PI values. The scarcity of PI data is likely due to the challenges of obtaining reliable PI from lavas (see, e.g., Tauxe, 2010, for details). We do not set threshold values for statistical parameters of the gathered site-mean directions, which are the number of specimens (N) utilized to compute the mean direction, the 95% angle of confidence (α 95 ), and the precision parameter (k). A significant portion of volcanic data does not fulfill high statistical standards (see Table 1), which is to be expected when dealing with paleomagnetic data in a transitional field state. The 42 site-mean PI values were obtained by different methods, including: the double heating Thellier technique (Thellier & Thellier, 1959) with various modifications (Aitken et al., 1988;Coe, 1967;Tauxe & Staudigel, 2004); the Microwave approach (Hill & Shaw, 1999); Shaw method (Shaw, 1974) and its variants (e.g., Tsunakawa & Shaw, 1994); and the multispecimen technique (Dekkers & Böhnel, 2006). Beside absolute PI data, five relative PI data are available from the Azores (Ricci et al., 2020). These data are scaled to absolute values in this study, following the approach described in Section 2.2.5. The quality of the current PI data was evaluated using the quality of PI (Q PI ) metrics, which were first introduced by Biggin and Paterson (2014). The updated PINT database (Bono et al., 2022) increased these metrics from 8 to 10 criteria. In general, the criteria are based on age information, type and stability of remanent magnetization, rock magnetic properties, PI data (number of specimens, uncertainty, and reliability parameters), as well as its availability. The PINT database was used to score the sites of this study by Q PI , with a notice that the Q PI values for four sites from Guadeloupe and one site from Chile (see Table 1) were adjusted in this study. Additionally, Q PI of the scaled intensity data from the five Azores sites was assigned in this study. These sites get a Q PI score of 1/10 because, aside from the age-related information, they lack information on the other nine requirements. The Q PI in the 42 sites have a range of 1/10-6/10, with a median value of 3/10. This suggests that the majority of PI data recording MBR have low Q PI values. The highest Q PI values for PI data come from Antarctica (Lawrence et al., 2009) and NW-USA (Lhuillier et al., 2017) sites, which have 6/10 and 5/10, respectively.
We use VADM and VGP parameters to represent intensity and directional changes, respectively, from our compilation of lava data over the past 900-700 ka. Here we consider the transitional polarity state by the commonly used range of VGP latitude values between +45° and −45°. Most of the sites between 900 and 800 ka have reverse  ( (Alken et al., 2021). Only nine VADM data exist between 762 and 700 ka, six of which have a mean of 6.0 × 10 22 Am 2 , showing that the magnetic field has most likely stabilized. This can also be seen in the VGP data, where the majority of the data during this time period is of typical normal polarity. One site from the West Eiffel, Germany (Singer et al., 2008), with a calibrated 40 Ar/ 39 Ar age of 726 ka, exhibits unusual behavior, with transitional polarity state (VGP = 35°N) and low VADM (2.5 × 10 22 Am 2 ). This anomalous behavior could, however, be connected to the site's significant 2σ age uncertainty of 38 ka. Due to large age errors and still insufficient data set, it is impossible to pin down the beginning and end of the MBR using solely lava locations. This necessitates the use of continuous sediment records.

Compilation and Spatio-Temporal Distribution
During the data compilation procedure for sediments, we focus on records with magnetic field components on a time scale. Records with paleomagnetic data on a depth scale that clearly specify the SR and the depth at which the MBR occurred have also been compiled. In this case, linear interpolation was used to convert to a time scale considering constant depositional rate. Other cases in which data is only available on a depth scale are not taken into account. Therefore no paleomagnetic data from the MBD97 (Love & Mazaud, 1997) database has been considered. Paleomagnetic data that have been reported from deep-sea sediments, shallow-sea sediments, and lake sediments (including lacustrine deposits) are only considered during data compilation. We do not use paleomagnetic data reported from Loess deposits, because doubts have been raised about the ability of these sediments to record the magnetic field signal during reversals (Wang et al., 2014;Wu et al., 2016;Zhou et al., 2014).
Under these selection criteria, we compiled 68 sediment records, of which 56 are deep-sea sediments and the remaining 12 are shallow-sea and lake sediments. Table 2 contains a list of these sediment data and the geographic distribution is shown in Figure 2a. The majority are located in three regions from the north-eastern hemisphere (North Atlantic, Western Equatorial Pacific, and North Pacific), which account for over 92% of the total. Note that this is in contrast to the lava sites distribution (see Figure 1a), where most of data are distributed in the western hemisphere. As illustrated in Figure 2a, the southern hemisphere contributes just a modest amount of data; data are not available in the southern half of South America, Africa, and Australia. A number of 28 sediment cores have full vector paleomagnetic data, 26 have only RPI, and the other cores have one or two magnetic field components.
Temporal distributions of sediment data are plotted in Figure 2b where the number of data are shown in 2 ka years bins. There are 13,068 declination, 13,454 inclination, and 16,897 intensity data points in total. Note that we only display data that have been accepted in this study. Between 760 and 800 ka, the period when the MBR was supposed to occur, almost 30% of the full-vector (declination, inclination, and RPI) field data (5,177) are concentrated. The North Atlantic region holds the majority of these data, with 10,067, 10,067, and 12,262 data for declination, inclination, and intensity, respectively, equally distributed over the studied period. These data account for ∼74% of the whole data set. Most of the N. Atlantic data are distributed between 800 and 780 ka.  (1977)  The Western Equatorial Pacific is the second largest data-rich region, with data highly concentrated between 790 and 760 ka. It is worth noting, however, that it only amounts to 13% of the North Atlantic data and 9% of the whole data set. In the North Pacific, like the Western Equatorial Pacific, the majority of the data are centered between 780 and 760 ka. The rest of the regions (Antarctica, Eastern Northern Pacific, Mediterranean, Equatorial Indian Ocean, Eastern Equatorial Pacific, Caribbean Sea, and Russian Arctic) together only have 947 declination, 1,251 inclination, and 2,137 intensity data points. The larger part of these data fall into the Brunhes chron, more precisely between 800 and 700 ka.

Timescale Reliability and Updates
The gathered sediment records were dated using a variety of methods. We did not consider records that were dated purely by comparing the magnetic field direction to a reference record of known age derived from global or local data sets, because this method is not independent of the magnetic field variations that we want to investigate here. Polarity reversals (or magnetostratigraphy) and RPI correlation, are two other dating methods that depend on the magnetic field. In the polarity reversal method (POL), conversion from the depth-scale to the age-scale is generally performed by comparing the directional data of a specific record to one of the geomagnetic polarity timescale references (Cande & Kent, 1995;Gradstein et al., 2012;Hilgen, 1991;Shackleton et al., 1990). The conversion is carried out by assuming a constant SR and employing one or a few tie points. The reference in the RPI dating method is a regional or global RPI stack or model, such as SINT-800 (Guyodo & Valet, 1999), PISO-1500 (Channell et al., 2009), or PADM2M (Ziegler et al., 2011). In RPI dating, the timescale is derived from a large number of tie points, which is typically greater than in the POL approach, and it does not assume constant SR along the record. A total of 46 records were dependently dated, with 37 using POL and 9 using RPI dating ( Table 2). Seventeen of the 37 POL age cores have only employed the POL method. Knowing the age of the MBR is required for POL dating, as this age is used as a primary (and in some cases unique) tie point. Therefore we do not include these 17 records in our following analyses, where we aim for independent estimates of MBR timing and duration, but keep them in our compilation for potential future use. The remaining 20 records combine POL with other dating methods and we include them. Moreover, we retain the 9 RPI-dated records since they do not require the MBR age and use more tie points. Although RPI does not vary synchronously over the globe, it is clearly influenced by dipole changes (see e.g., M. C. Brown et al., 2018;Korte et al., 2019) which might justify the age determination by this method dependent on magnetic field changes.
The preferable dating methods use features that are not tied to the Earth's magnetic field, which we call independent properties. The oxygen isotope (δ 18 O) method, which was used on 29 of the 68 sediment cores, is the most common independent method. Global δ 18 O stacks can be used to establish the paleoclimate by combining the mean time series of δ 18 O data from numerous places throughout the world. SPECMAP (Imbrie et al., 1984) that covers the last 750 Kyr, S95  that covers the last 6 Myr, and LR04 (Lisiecki & Raymo, 2005) that covers the last 5.3 Myr, are examples of such stacks. The age model of these stacks was built through orbital tuning; consequently, one may date a record by comparing δ 18 O data to one of these stacks. The highest resolution is given by LR04, because it was created using 57 benthic δ 18 O records published from globally well-distributed (in latitude, longitude, and depth) sites from the Atlantic, Pacific, and Indian Ocean. Its age model was determined by aligning the δ 18 O stack with an ice model based on summer insolation at 65°N on the 21st of June (Lisiecki & Raymo, 2005). No modifications were done to sediment cores having LR04-assigned age data. For records that used SPECMAP and S95, conversion of their age to the LR04 age model can be done using the conversion relation available online (https://lorraine-lisiecki.com/LR04_age_conversion.txt). It is worth noting that a 5-Myr probabilistic stack was created from 180 δ 18 O records, and can be used for dating (Ahn et al., 2017). However, because it was not used in any of the assembled records, we used LR04 and did not convert all data to this probabilistic stack. In addition to δ 18 O, a total of 13 cores have been dated by comparing their rock magnetic properties, such as MS, S ratio , ARM, and ARM/MS, to one of the δ 18 O reference records or stacks. Other independent dating techniques, such as palaeontology (2 records), electron spin resonance (2), stratigraphy (1), tephrostratigraphy (4), diatom studies (1), gamma-ray attenuation porosity evaluator (1), and barium/titanium (Ba/Ti) ratio (1), have been utilized. We cannot go into details here, but we accepted all of these records.
During data compilation, published updates on the original age scale have been applied for four records (MD97-2143, OB, RC10-167, and HS), see supplement for details.

Temporal Resolution
In different aquatic environments, sediments are deposited at various SR. Deep-sea sediments have generally the poorest resolution, with typical SR on the scale of 1 cm/Kyr. Coastal ocean and shallow sea sediments have resolutions of a century to a millennium (10 cm/Kyr), whereas lake cores have resolutions of a decade to a century (100 cm/Kyr). The SR in the compiled records (Table 2) ranges from 0.5 to 89 cm/Kyr, with 40 records having SR < 5 cm/Kyr, 4 having SR range 5-10 cm/Kyr, and 22 records having SR ≥ 10 cm/Kyr, that is, more than half of all the sediments data are of low resolution. There are two records (807A and KS752) without SR information provided. Regarding only the accepted records (which are 38): 22 records have SR < 5 cm/Kyr; 3 in the range 5-10 cm/Kyr; and 13 records have SR ≥ 10 cm/Kyr. This indicates that there are no systematic differences in SR between accepted and non-accepted records. The SR distribution is shown in Figure 2a, where we can see that North Pacific and Mediterranean records have the highest SR, and the maximum SR value of 89 cm/Kyr is found in three shallow marine-records (YGC, YT, and CBCS) from the North Pacific. The second largest SR, in the range of 25-60 cm/Kyr, are found in three Mediterranean records: the lacustrine section HS; shallow-marine sediment record CB, and lake record International Continental Scientific Drilling Program (ICDP) 5045-1. North Atlantic records have semi-uniformly high SRs, mostly around ≥10 cm/Kyr. The data from the Western Equatorial Pacific, on the other hand, show the lowest SR, which does not surpass 3 cm/Kyr, with the exception of three cores (ODP 769B, ODP 769A, and ODP 769B, SR = 10 cm/Kyr).
The temporal resolution of a sediment record is also influenced by the sampling method. U-channel and pass-through measurements lead to stronger smoothing of the signal than measurements on discrete samples. We assessed the temporal resolution of the compiled records by performing smoothing spline analysis (Panovska et al., 2012). This analysis provides estimates of the smoothing time (Ts), which is the length of time that can be resolved in a sediment record. It is performed on each magnetic field component separately and can also help to assess the random uncertainty in the records. To perform the analysis, a prior Ts value, which represents a threshold value for the time resolution during the spline analysis is required. We determine it from the SR and lock-in depth; SR values are used as listed in Table 2 while the lock-in depth is set to 10 cm for all sediment records.
After obtaining a spline model of a given component, the Ts is determined by estimating the full width at half maximum height of a resolving kernel that diagnoses the degree of smoothing inherent in the spline model (for more information see Panovska et al., 2012). We calculated an average Ts for each record using smoothing spline analysis on 79 individual components established from the 38 acceptable sediment records. An example of the smoothing spline analysis is given in Supporting Information S1 ( Figure S2). The distribution of the Ts values is presented in Figure 2c. The resultant Ts values range from 0.2 to 18.8 ka, with the median being ∼1.0 ka. The longest smoothing time (18.8 ka) was calculated from the Caribbean Sea record (CAS16-24PC). We note that, for the time range 900-700 ka, there are only 18 data points in this record. We found no significant changes in the resultant Ts values when we decreased the lock-in depth parameter for this record to 0 or 5 cm. Except for this record, all analyzed records have Ts < 10 ka, 25 have Ts < 2 ka, and 18 records have Ts ≦ 1 ka. We may deduce from the temporal resolution analyses that ∼50% of the sediment records selected in this study can capture magnetic field behavior during the MBR with a resolution better than 1 ka. It should be noted that the majority of the high resolution records are from the North Atlantic region.

Declination Orientation
We ascertain if the declination component in the records is azimuthally oriented. Since we are dealing with two antipodes polarities, median declination has to be south and north during the Matuyama reverse and Brunhes normal polarity chrons, respectively. Records that do not match this requirement are rotated in this study. In order to achieve this, mean average of the declination data over the two chrons were calculated, eliminating data of transitional polarity (VGP latitude between +45° and −45°). After that, the declination values are rotated, so as to have ∼180° and ∼360° median declination over Matuyama and Brunhes chrons, respectively. This approach is quite comparable to earlier ones (M. C. Brown et al., 2018;Panovska et al., 2021). Here, four records (ODP 980, ODP 984, MD97-2143, and CADO) were rotated. Figure S3 in Supporting Information S1 provides an example of the rotation process for record ODP 984.

Scaling of RPI Data
RPI values were estimated using different normalization factors: ARM, MS, and IRM. To scale RPI data, we must compare them, ideally to PI values from nearby volcanic data of known age, or to a calibrated stack or model constructed from RPI records (e.g., SINT-800, PISO-1500, PADM2M). The volcanic data are too sparse in the past 900-700 ka, as seen in Figure 1. We chose to use the 2,000 ka paleomagnetic axial dipole moment model PADM2M (Ziegler et al., 2011) that is built from 96,032 absolute and relative PI data obtained from archeomagnetic, igneous, and sedimentary archives.
Scaling was done by the approach described by M. C. Brown et al. (2018), using: where F rescale is the rescaled RPI data at each time step of a record, F RPI is the RPI data at each time step, M RPI is the median RPI over the full time range of the input sediment record, and M p is the median absolute PI over the time range of the sediment record, obtained from PADM2M. In order to convert VADM to intensity at each sediment core geographic location, we use the equation (Tauxe, 2010) = VADM 0∕ ( where a is the average radius of the Earth, μ 0 is the permeability of free space, F is absolute PI, and θ is the geographic co-latitude. We investigated scaling the sediment records across a 900-700 ka segment in addition to the full length. With the exception of one Balkan record (ICDP5045-1), the scaling results were nearly identical in both timeframes ( Figure S4 in Supporting Information S1). When employing the full-length timescale, the M P /M RPI ratio for this record is 92.1, whereas using the 900-700 ka segment yields a ratio of 19.5. Except for ICDP5045-1, RPI data of all records are scaled over the full period covered by the records.

Regional Consistency of Sediments
After the dating check and RPI scaling, we have further selected the sediment records to ensure a high-quality, consistent data set for our investigation of the MBR. Three aspects are considered in this evaluation process: (a) regional consistency, (b) timing of the polarity transition, and (c) range of the scaled intensity values. Several regions contain more than one record, located within relatively close distances, so that similar geomagnetic variations are expected (e.g., Korte et al., 2019). We rejected records that display clearly different behavior from the others in regional proximity, as the credibility of the recovered paleomagnetic signal and/or the ages have to be questioned. Rejection was made also based on the polarity transition's timing: records that show a polarity reversal in VGP latitude at periods that are strongly different from all others were not considered in the following analysis. Records with very unusual intensity values are also discarded since they may indicate that the normalization might be insufficient and, for example, some climatic influences were not removed. It is not possible to define strict criteria because we do not know how irregular the magnetic field was around the MBR. Our decisions therefore are somewhat subjective based on current understanding of the reversal. Our choices are described in detail below, illustrated by plots of the time series of intensity ( Figure 3) and direction (Figure 4), expressed by VADM and VGP latitude, of the records combined according to their regional distribution (Table 2).
In comparison to the rest of the North Atlantic records (Figures 3a and 4a), we observe that sites U1305 (Mazaud et al., 2012) and U1302-03  show different VADM and VGP behavior. When three cores of U1305 were examined in detail (see Channell, 2017), hole-to-hole variations were noted. Two other records, U1304 and U1306, have similar hole-to-hole differences (Channell, 2017), but their VGP curves are not substantially different from the rest of the region's records. Note that Site U1305 was compared with Sites U1306 and U1307 (Text S2 in Supporting Information S1), as all three are from Eirik Drift and are ∼70 km apart. The only instance for the Labrador Sea is Site U1302-03, which spans 764 ka . Its paleomagnetic results are inconsistent with the North Atlantic record, but agree with Site U1305. Nevertheless, we reject these two sites from the North Atlantic data collection. Record ODP1063 from the Mid-Latitude North Atlantic area is also not accepted, because the primary polarity transition is seen at ∼762 ka (Figure 4j). From the Southern Atlantic, RPI data of core KS752  had a significant problem that the natural remanent magnetization (NRM)/(ARM) and NRM/MS ratios are not the same, implying that it is characterized by non homogeneous magnetic mineralogy content. Therefore, we do not consider the VADM data ( Figure 3k) derived from this core.
For the Pacific, according to the regional consistency and intensity range criteria, records ODP 767B (Oda et al., 2000) and MD012380 (Huang et al., 2009) from the Western Equatorial Pacific have distinct behavior and/or abnormal intensity values at 770-760 ka and 820-810 ka, and have been rejected (Figure 3b). We chose 769A from the two Sulu Sea holes (ODP 769A and 769B; Oda et al., 2000) because it is more consistent with other records of the region. The VADM curve of records NGC69 (Yamazaki, 1999) and RC10-167 (Kent & Opdyke, 1977;Meynadier et al., 1995) are different from rest of North Pacific records ( Figure 3c) and therefore not considered. Record ODP1021.B  from the Eastern North Pacific has two main intensity lows, at ∼825 ka and ∼730 ka (Figure 3e). From the VGP latitude (Figure 4e), the main polarity transition occurred at ∼730 ka. We reject this record since the time is very different from the widely accepted ages of the MBR (Channell et al., 2010;Shackleton et al., 1990;Singer et al., 2019;Valet et al., 2014). Similarly, record ODP851 from the Eastern Equatorial Pacific is discarded because its VADM and VGP data (Figures 3h and 4h) indicate the MBR at ∼730 ka .
In the Mediterranean, we notice that record ICDP5045-1 from the Balkan (Just et al., 2019) (Figure 3f) has distinct behavior in comparison to records HS and LC07 (from Italy), which are located ∼700 km west from it. VGP latitude of the ICDP5045-1 is also different from HS and LC07 (Figure 4f), hence we do not accept it. The long Integrated Ocean Drilling Program U1396 record from the Caribbean Sea (Hatfield et al., 2021) is not accepted due to the poor paleomagnetic data quality of unit 2a, which includes the time of interest, as described in the original paper. This might be evident in the VGP curve derived from this record (Figure 4i). The main polarity transition for record ICDP-5011 from Lake El'gygytgyn, East Russian Arctic, is at ∼780 ka (Nowaczyk , which is consistent with other investigated regions. Although Nowaczyk et al. (2013) pointed out that there is no evidence for geomagnetic excursions during the Brunhes chron, we accept this record because the MBR appears to have been faithfully recorded.

VADM Stack
To explore the stability of the magnetic field strength during the MBR, we create a set of regional VADM stacks and one global stack. In the regional stacks, our input curves are accepted records (which are 35) while the global stack was then built from these regional stacks. We should point out that individual records were not temporally aligned prior to establishing regional stacks; instead, we used their individual ages (see Section 2.2.2). This strategy differs from the commonly used approach (Channell et al., 2009;Guyodo & Valet, 1999;Laj et al., 2004;Valet et al., 2005) to align all records with a master record, which results in an average curve with lower uncertainty. However, some true variability of the non-dipole field may be lost by this process, and we do not employ this technique here.
VADM regional and global stacks were calculated by a bootstrap resampling technique where the data are resampled 5,000 times. VADM values of individual records from same region are combined into one data set. From this data set, 5,000 realizations are obtained by bootstrapping. For each of these realizations, a smoothing spline fit is obtained. The mean of these 5,000 fits is then calculated and produced at 500-year time steps, and this mean represents regional VADM stack. Confidence interval of each regional VADM stacks is provided at 2σ level. Seven regional VADM stacks were constructed using this approach (Figure 5a). Uncertainties for the regional stacks are not shown to avoid cluttering the figure. Indications for the uncertainties come from the comparison of the individual curves in Figure 3. Note that the bootstrap uncertainties may depend on the number of records and might not be fully reliable. In the western hemisphere, two regions (the Eastern Equatorial Pacific and the Caribbean Sea) each have only one acceptable record (Figures 3h and 3i). These two regions were not subjected to bootstrap analysis; instead, we kept their original VADM data and treated them as representative of the two regions. We maintain them because sediment records from the western hemisphere are few (Figure 2a). A global VADM stack with its 95% confidence interval was then generated by bootstrap resampling of the regional VADM stacks and the two representative regions, as shown in Figures 5a-5c. Note that the resampling does not take the age uncertainties of individual records into account, so that the confidence interval likely underestimates the true uncertainty. A global stack can also be build from all individual VADM data points of the 35 accepted records. However, the resulting curve appears to be biased toward the regional North Atlantic and Western Equatorial Pacific stacks, which are the most intensively covered regions (see Figure S5 in Supporting Information S1). We consider the stack based on regions as more representative for the global VADM variations.
The comparison of the regional stacks ( Figure 5a) shows that the North Atlantic, North Pacific, and Antarctic stacks have the highest resolution, while the Western Equatorial Pacific, Eastern Equatorial Pacific, and Caribbean sea have the lowest resolution. The resolution is determined by the SR and the quantity of input records, among other factors. At ∼820 ka, most of the regional stacks reveal a gradual decline in VADM. In addition, all regions experienced two major VADM drops at ∼800-790 ka and 780-775 ka. Higher field intensity of different magnitude is recorded between these two minima, with the highest intensity being observed in the East Indian Ocean stack. The rates of magnetic field decay and growth vary between the curves during the entire 900-700 ka interval. For example, regions of the North Pacific, Mediterranean, and Eastern Equatorial Pacific display rapid intensity declines between 810 and 800 ka, whilst the North Atlantic stack shows more gradual decay. The magnetic field intensity increased rapidly after 775 ka, reaching generally higher levels than prior to 775 ka. This quick growth may be seen in the Eastern Equatorial Pacific curve, where the magnetic field intensity doubled in less than 1 ka. The rest of the regions display rapid growth as well.
The global VADM stack (Figures 5a-5c) shows a small VADM drop at ∼820 ka (to 4.2 × 10 22 Am 2 ), but the field increased back to its pre-820 ka VADM after 10 Kyr. Then, between 810 and 795 ka, the magnetic field intensity gradually decreased until it reached 3.7 × 10 22 Am 2 . After a short increase, the decrease to the lowest values observed in our stack is seen at 780 ka, with the VADM dropping to ∼2 × 10 22 Am 2 . After 780 ka, the axial dipole moment rises steeply and only ∼23 Kyr later (at 757-750 ka) the field reached ∼10 × 10 22 Am 2 . This broad peak is about double as high as the values observed in the period 900-810 ka. A comparison of our new stack to some Figure 5. (a) Seven regional virtual axial dipole moment (VADM) stacks, constructed from 35 accepted records by bootstrap resampling, and two individual VADM curves. The global VADM stack (thick black line) with its 95% confidence interval (thin black lines) constructed from these regional stacks is shown. (b) The global VADM stack is compared to previous stacks: SINT-2000(Valet et al., 2005; PISO-1500 (Channell et al., 2009); and the PADM2M model (Ziegler et al., 2011). (c) Comparison with the volcanic VADM data compiled and evaluated in this study, indicated by the same symbols as in Figure 1.
previous stacks that include the 900-700 ka period, SINT-2000, PSIO-1500, and PADM2M, is given in Figure 5b. The SINT-2000 stack (Valet et al., 2005) was built from 43 records, 33 of which are younger than 800 ka while the remaining 10 are between 2,000 and 800 ka. The PISO-1500 RPI stack (Channell et al., 2009) that covers the past 1,500 ka was created by tandem matching and stacking of 13 paired RPI data and δ 18 O records, with a note that most of these records are located in the North Atlantic. PADM2M (Ziegler et al., 2011) is a regularized cubic B-spline VADM model, constructed from global absolute and scaled intensity data. The VADM values in our stack range from ∼2 to ∼10 × 10 22 Am 2 , which is comparable to SINT-2000. PADM2M, on the other hand, has relatively lower axial dipole moment strength, which does not exceed ∼8 × 10 22 Am 2 . Between 900 and 780 ka, PISO-1500 has greater amplitudes than our stack, but after 780 ka, they are comparable. We note that our global stack is substantially more consistent to SINT-2000 and PADM2M, while the peaks between 860 and 800 ka in PISO-1500 are similar to intensity maxima found when stacking the individual records dominated by Atlantic data ( Figure S6 in Supporting Information S1). The main VADM decrease that occurred at 780 ka in our stack is also seen in PADM2M at the same time, while SINT-2000 and PISO-1500 show the main drop about 5 Kyr later (at 775 ka). The previous stacks drop to a minimum value of ∼1 × 10 22 Am 2 , which is half of the value observed in our stack. Hence, either the scaling of our RPI is systematically too high, or, more likely, our global VADM stack contains non-dipole field contributions, which cannot cancel during the transition as intensity or VADM are always positive values regardless of field polarity. Note that the global, pure axial dipole moment has to be zero at some point in time during the MBR.
When the global VADM stack is compared to the absolute PI-derived VADM data of the 42 lava sites (Figure 5c), reasonable consistency is observed for most regions, considering the uncertainties in age data, PI results, and in the stack itself. Several data from Tahiti, Guadeloupe, Chile, the Canary Islands, Azores Islands, and NW-USA match the stack at 875 ka, ∼790-770 ka, and 750-700 ka. The low VADM values obtained at two sites from Canary Islands and Chile in ∼820 ka correspond to the global VADM drop indicated in the global VADM stack. Because absolute PI data are more reliable than relative intensities, this consistency demonstrates that the conversion of RPI to absolute intensity was done with a suitable scaling factor.
VADM values for some data are higher than the stack values. Three Mexican VADM data, for instance, are much higher than stack values at 800 ka and 740 ka, and the intensity at one site in Guadeloupe and one site on the Canary Islands are higher than the stack value at 785 ka. However, when the age inaccuracy is considered, the aforementioned disparities are mostly insignificant. Five sites from Hawaii and the Azores exhibit very low VADM values (0.4-1.4 × 10 22 Am 2 ) between 770 and 760 ka, which are not reached by our stacked curve nor in previous stacks. In addition, VADM values in two Hawaiian sites with ages 725 and 748 ka are statistically lower than the stack's values. In this comparison we have to keep in mind that each individual VADM value translates non-dipole field contributions to axial dipole moment, which might account for some of the observed discrepancies.

VGP Latitude Stack
Although directional field variations are expected to differ even more strongly regionally than VADMs, and this is particularly true when the dipole dominance is weak (M. C. Brown et al., 2018;Korte et al., 2019), we construct regional and global VGP stacks to characterize the global field direction change during the MBR. As in VADM, the VGP data of individual records were not temporally aligned and six regional stacks were created by bootstrap resampling. Two individual records (Russian Arctic and Caribbean Sea) were maintained as representative records for the region, as shown in Figure 6a. A global VGP stack was created from these stacks and example curves ( Figure 6). As in VADM, a global VGP stack built from all individual records rather than the regional stacks is similar to the North Atlantic regional stack (see Figure S6 in Supporting Information S1) and is not considered here.
All curves indicate reverse polarity directions from 900 to 835 ka (Figure 6a). The North Atlantic curve does not reach as high southern VGP latitudes as the other ones. This observation can be linked to the SR and Ts of regional data. Records with a high SR and a low Ts (as in the case of the North Atlantic) should have large amplitudes of directional changes, whereas records with a low SR and a high Ts should have smooth directional variations. Moreover, as seen in Figure 2b, the North Atlantic region has the most input data. Looking at the VGP distributions of the various regions ( Figure S7 in Supporting Information S1) we see that records with the most data points and a high SR have the widest range of VGP latitudes, including rather low values. Sediment cores with high sedimentation rates more likely record short-lived directional anomalies, for example, excursions, so their VGPs are more dispersed. Therefore, this regional difference in the VGP latitudes likely does not have magnetic origin. All regions indicate that a sequence of transitional field states developed after 835 ka and lasted until ∼770 ka, when the polarity state turned to normal. Unstable field behavior is first observed in the Russian Arctic curve between 831 and 804 ka. However, more records from this region are needed to corroborate this result, which is based on only one record.  (Figures 6a and 6b) exhibits the majority of the previously mentioned characteristics. The field polarity was reverse between 900 and ∼800 ka, followed by a transitional state between 800 and 773 ka, when the field became normal polarity. The VGP latitudes obtained from 107 lava sites are also shown (Figure 6b). In general, there is good agreement, particularly between the transitional direction seen in the stack at ∼784-778 ka and the low VGP latitude values seen in Guadeloupe, Chile, and NW-USA. Low VGP latitude values from Tahiti and Hawaii accompanied the transitional stage noted at 795 ka and 785 ka in the global stack (Figure 6b). The two stable polarity epochs in the stack in general agree well with the lava VGPs. However, some differences between the stack and the lava data exist. For example, between 870 and 770 ka, the normal polarity directions seen at Tahiti, Guadeloupe, the Canary Islands, Iceland, Mexico, and Antarctica differ from the reverse polarity direction obtained from the stack. Between 795 and 780 ka, four regions of lavas (Tahiti, Guadeloupe, Canary, and Azores Islands) have reversed directions, whereas the stack shows transitional polarities. The lava VGP data also reveal a transitional polarity state around ∼775-770 ka, where the stack has reached stable polarity. Most of the disparities are small when the age uncertainties are taken into account. Moreover, like VADMs, individual VGPs contain non-dipole field contributions and represent the local directions rather than a global field property. The field in the diverging regions could deviate particularly strongly from the global average, which should be given by the stack, but also might be somewhat biased.

Paleosecular Variation Index
We use the paleosecular variation index (P i ) suggested by Panovska and Constable (2017) to investigate periods of unstable field during the past 900-700 ka quantitatively. This index measures the departure of the VGP latitude and virtual dipole moment (VDM) from the geographic pole and the present day dipole moment strength, respectively, that is, takes both field direction and intensity information into account to characterize field complexity and variability. Low index values represent a stable, strong, dipole-dominated magnetic field, while high index values characterize unstable and weak magnetic field. To appropriately quantify a threshold value dividing stable from unstable field states, Panovska and Constable (2017) tested this index on recent field models and paleomagnetic data and suggested a threshold of 0.5 for transitional events.
We calculated P i on local, regional and global scales. The time series of three P i curves calculated from three sediment records are illustrated in Figure 7a, as examples from which the temporal evolution of the local magnetic field state can be investigated. P i curves of all records with accepted RPI and directional results are shown in Supporting Information S1 ( Figure S8). Most of the analyzed records show transitional field state (P i ≥ 0.5) at more than one period. For instance, the P i curve of site ODP 983 from the North Atlantic shows unstable magnetic field behavior at 859-855 ka, 853-850 ka, and 776-766 ka (Figure 7a). Similarly, record MD90-961 from the Equatorial Indian Ocean shows transitional states at 800-793 ka and 779-771 ka. Moreover, the North Pacific record SO202-1 has transitional states at 900-893 ka and 800-796 ka. From 800 to 796 ka, P i ≥ 0.5 was noted in only three data points, with a high value at 796 ka. The large smoothing time of ∼8 ka of this record (Section 2.2.3) should be noted. Bootstrap resampling is performed on P i values of individual records to construct regional P i stacks using the same methodology as for the VADM and VGP stacks. However, doing so results in unrealistic fits for the North Pacific region, which is caused by one very high P i value observed in record SO202-1 (Figure 7a). In order to obtain a more realistic regional P i stack for the North Pacific region, we consider only smoothing spline fits with positive P i values, and apply this procedure to all regions. We note that it is also possible to use VGP ( Figure 6b) and VDM stacks to calculate regional P i stacks, which gives comparable results (see Figure S9 in Supporting Information S1). Except for the Caribbean Sea, all of the analyzed regions show high P i values between 800 and 775 ka that reflects the global transitional field state. The lack of transitional field in record CAS16-24PC (see Figures S8 and S9 in Supporting Information S1) could be due to its low SR (1.7 cm/Kyr), which makes it difficult to record quick field changes. The transitional state is first seen in the Equatorial Indian Ocean, Western Equatorial Pacific, North Pacific, and Antarctic regions at ∼803-795 ka. Five regions (North Atlantic, Western Equatorial Pacific, Antarctica, Mediterranean, Equatorial Indian Ocean) then have P i ≥ 0.5 at 795-790 ka. Between 785 and 770 ka, transitional field states are found almost simultaneously in six regions (except for the Caribbean Sea) which, as will be discussed below, reflects the main phase of the MBR. Our results suggest that the field in the Mediterranean region has transitional behavior between 792 and 781 ka, which is longer than the duration documented in different individual places in this region (Just et al., 2019;Macrì et al., 2018;Sagnotti et al., 2014). We note here the Maffei et al. (2021) deduction that the extreme variations recorded in the Sulmona site (Sagnotti et al., 2014) are spatially localized and only represent a transient feature during a longer transition.

Timing and Duration
A global P i stack was established (Figure 7c) from the five regional P i stacks and two P i curves ( Figure 7b). We propose that the time and duration of the MBR can best be estimated using this stack because it considers the full-vector magnetic field information from 18 globally distributed, high-quality sediment records. The field polarity was stable reverse between 900 and 801 ka, followed by a transitional state between 800 and 770 ka, when four intervals of strong field instability occurred, indicated by four P i peaks. The first three peaks are observed in the interval 800-785 ka, and they are characterized by small drops in the axial dipole moment ( Figure 5c) and transitional VGP latitudes (Figure 6b) in some regions, but the field resumed its previous orientation. These short intervals of magnetic field instability can be considered as precursors (Hartl & Tauxe, 1996;Kent & Schneider, 1995;Valet et al., 2012), which are seen regionally, similar to some excursions (see, e.g., Panovska et al., 2019). Although we cannot fully rule out some distortion from age uncertainties in individual records, that facts that the peaks are found in different regional stacks coming from independently dated records and occur at times of moderately low VADM ( Figure 5) supports our interpretation. The fourth major phase, which lasts between 780 and 770 ka, is the main polarity transition period, after which the polarity quickly becomes stable in the present day direction. Considering all four peaks, we deduce that the MBR evolved over ∼30 Kyr (800-770 ka). We do not find indications for a rebound state (Valet et al., 2012) with our data selection, which is in agreement with the results by Singer et al. (2019).
From the above, the end stage of the Matuyama-Brunhes polarity reversal is around 770 ka. This estimate is consistent with previous findings (Channell et al., 2020;Jouzel et al., 2007;Raisbeck et al., 2006;Singer et al., 2019;Valet et al., 2014) where the mid-point age of the MBR was determined at 773 ka. We find an even longer duration than Singer et al. (2019), who suggested that the MBR evolved over ∼22 ka, from 795 to 773 ka, with main instabilities in three successive periods around ∼795, ∼784, and 773 ka. This is probably due to the substantially lower number of data used by Singer et al. (2019). The inferred duration is much longer than the range of 1-10 ka observed regionally in the IMMAB4 geomagnetic field model (Leonhardt & Fabian, 2007), which spans the interval 795-764 ka. The discrepancy can probably be attributed to IMMAB4's limited data set, which included only three sediment records (609B, 664d, and V16-58) and one volcanic region (La Palma). Actually, these three sediment records are purely POL dated, which aligns the MBR signal, and they were not accepted in our study.

Relation Between Local Duration and Latitude
The local MBR duration of 18 individual records, defined from P i curves ( Figure S8 in Supporting Information S1), is displayed versus geographic latitude to see if the duration computed in this study depends on latitude. A latitudinal dependence of duration had been suggested by Clement (2004) on four field reversals, including the MBR which was studied in 20 records. Clement (2004) found that the duration increased with latitude, ranging from 2 Kyr at low latitude sites to 12 Kyr at latitudes 50-60° north and south. The duration in Clement (2004) was calculated using only directional data, by estimating the thickness of the stratigraphic section where the polarity switch was observed.
Using the same methodology, Leonhardt and Fabian (2007) evaluated the duration for identical locations (as in Clement (2004)) as predicted from their global IMMAB4 model. The duration estimated from the model is broadly similar with the duration derived from paleomagnetic data, confirming the suggestion made by Clement (2004). Latitudinal dependence of duration has also been found in the numerical dynamo model T4 of Wicht et al. (2009). In order to compare the durations from our full-vector field information with previous values, all these different estimates are plotted in Figure 8. To distinguish between excursions and reversals, Wicht et al. (2009) stipulated that field reversals be separated by a stable period of duration T s that does not include any excursion. They used three values of T s to investigate the latitudinal dependence of duration; we include the duration obtained from T s = 10 Kyr in Figure 8 because it is closest to the durations determined from paleomagnetic data. As seen there, the durations generated by the numerical model are longer at higher latitudes and at the equator, and shorter at mid-latitude regions.
We first remark that the duration range that we find is different, with durations ranging from 1 to 35 Kyr in our study, which is roughly three times longer than the 2-10 Kyr durations of Clement (2004) and Leonhardt and Fabian (2007). The dynamo model-derived durations (Wicht et al., 2009), on the other hand, are around two times longer than ours and seven times longer than the Clement (2004) and Leonhardt and Fabian (2007) estimates. Note that in Figure 8, we have plotted the paleomagnetic sites used in this study and Clement (2004) and the duration derived from IMMAB4 model (Leonhardt & Fabian, 2007) scaled by SR to see if it has an effect on the calculated durations. We observe short and long durations with distinct SR in our data set, indicating no SR affect. The MBR durations of our high-latitude records (North Atlantic and Antarctica) are generally longer than those of low-to middle-latitude records (Western Equatorial Pacific, Equatorial Indian Ocean, North Pacific, and Mediterranean regions). The same applies to durations of Clement (2004) and Leonhardt and Fabian (2007). This indicates that duration increases with latitude, regardless of how it was calculated. Also, we note that one record from the Equatorial Indian Ocean (lat ∼ 5°N) has a duration of 35 Kyr, while another from the Western Equatorial Pacific (lat ∼ 9°N) has a local duration of 15 Kyr. These durations are longer than those reported in six mid-latitude records (lat ∼35-39°N), suggesting that low-and mid-latitude regions might have different durations. Indeed, this result is similar to the numerical model-derived curve ( Figure 8). This observation was believed to be an evidence of the presence of low-latitude reversed flux (Amit et al., 2010;Wicht, 2005;Wicht et al., 2009). However, as we also have two records of latitudes ∼16°N and ∼6°S with durations around one Kyr each, this hypothesis requires further paleomagnetic data to be confirmed. We also notice that two North Atlantic records give durations of ∼5 Kyr and ∼11 Kyr, which are comparable to low and mid-latitude results. More data, in particular from the southern hemisphere, are still needed to fully investigate any latitudinal dependence of MBR duration.

Conclusions
High quality paleomagnetic sediment records have been used to produce regional and global stacks of VADM, VGP latitude, and of the paleosecular variation index (P i ) for the past 900-700 ka, that is, the time interval around the Matuyama-Brunhes geomagnetic field reversal. Only well-dated records that have reasonable consistency with surrounding records (if any exist) and that also have no aberrant relative intensity values were accepted. As a result, only 38 sediment records were chosen from a total of 68. The accepted records have a reasonable spatial distribution, though the southern hemisphere is underrepresented, and there are only few records from mid-latitude regions. We calculated the MBR's timing and duration, and investigated the relationship between the local reversal duration and geographic latitude. Unlike previous approaches, we analyze these features by evaluating the full vector field data of the accepted records by use of the P i , which is computed from both field intensity and direction. At 800 ka, the magnetic field became unstable, indicating the start of the MBR, which ended at 770 ka. It has been shown, from regional and global VADM stacks, that the global axial dipole moment began to steadily drop ∼10 Kyr before the onset of the reversal. The 770 ka age for the end of the polarity reversal at a global scale determined in this work corresponds well with the mid-point age (773 ka) reported previously from high-quality and high SR marine cores. The ∼30-kyr MBR duration is even longer than the 22-kyr evolution time proposed recently based on a compilation of well-dated lava flows, sediments, and ice cores, according to which the MBR evolved between 795 and 773 ka. Our estimate's longer duration may be explained by the fact that it is based on a larger data set. The prior idea that local reversal duration is latitudinal dependent is generally confirmed in this study, as high-latitude records from the north and south hemisphere have longer durations than those found at low-to-medium latitudes. We also found indications for longer duration at very low latitudes similar to earlier results from a numerical dynamo simulation, but these results might not be robust as nearby records suggest short duration. More data from a wider range of locations is still required to establish the relation between MBR duration and latitude. The new data compilation presented here will be used to build a new global SH model of the Matuyama-Brunhes reversal.

Data Availability Statement
Data of regional and global stacks of the virtual dipole moment, virtual geomagnetic pole latitude, and paleosecular variation index, which were calculated in this study, as well as the paleomagnetic data of accepted sediment records compiled from previous studies can be found at https://earthref.org/ERDA/2545/.

Acknowledgments
The Alexander von Humboldt Foundation awarded A.N. Mahgoub a research fellowship, which he expressed thanks for. S. Panovska gratefully acknowledges the Discovery Fellowship at the GFZ Potsdam, Germany. We would like to express a deep appreciation to all authors who shared their data with us personally, or made the data available through supporting information and public databases. We also thank Leonardo Sagnotti and an anonymous reviewer for their constructive comments and suggestions. Open Access funding enabled and organized by Projekt DEAL.