A Late Cretaceous‐Eocene Geomagnetic Polarity Timescale (MQSD20) That Steadies Spreading Rates on Multiple Mid‐Ocean Ridge Flanks

Magnetic anomalies over mid‐ocean ridge flanks record the history of geomagnetic field reversals, and the width of magnetized crustal blocks can be combined with absolute dates to generate a Geomagnetic Polarity Timescale (GPTS). We update here the current GPTS for the Late Cretaceous‐Eocene (chrons C33–C13, ~84–33 Ma) by extending to several spreading centers the analysis that originally assumed smoothly varying spreading rates in the South Atlantic. We assembled magnetic anomaly tracks from the southern Pacific (23 ship tracks), the northern Pacific (35), the southern Atlantic (45), and the Indian Ocean (51). Tracks were projected onto plate tectonic flow lines, and distances to magnetic polarity block boundaries were estimated by fitting measured magnetic anomalies with a Monte Carlo algorithm that iteratively changed block model distances and anomaly skewness angles. Distance data from each track were then assembled in summary sets of block model distances over 13 ridge flank regions. We obtained a final MQSD20 GPTS with another Monte Carlo algorithm that iteratively perturbed ages of polarity chron boundaries to minimize the variability of spreading rates over all ridge flanks and fit an up‐to‐date set of radioisotopic dates. The MQSD20 GPTS highlights a major plate motion change at 50–45 Ma, when spreading rates decreased in the Indian Ocean as India collided with Eurasia while spreading rates increased in the South Atlantic and Northern Pacific and the Hawaii‐Emperor seamount chain changed its orientation.


Background
Accurate timescales are crucial to establish the history of tectonic plate motion and to determine past rates of change documented in the rock record. Magnetic measurements and radioisotopic dating of volcanic rocks led to the discovery of globally synchronous reversals of the Earth's magnetic field (see the historical overview of Glen, 1982). The history of these reversals was recorded by magnetic anomalies created by seafloor spreading on mid-ocean ridge flanks, allowing for the development of a geomagnetic polarity timescale (GPTS; acronyms used in this paper are listed in Table 1).
The first GPTS was based on the C-sequence (Late Cretaceous-present) magnetic anomalies in the South Atlantic with the assumption of a constant spreading rate determined by a single age tie at 3.35 Ma (Heirtzler et al., 1968). Kent (1992, 1995) derived an improved CK95 GPTS for the last 83 Ma that accounted for long-term variations in spreading rate. CK95 was based on a reference South Atlantic plate tectonic flow line defined by a detailed set of finite plate rotations (Cande et al., 1988). The finite rotations provided nine key magnetic polarity block model distances (BMDs); more detail was inserted from the relative widths of polarity blocks representative of uniform seafloor spreading in the South Atlantic and in faster spreading ridges of the North Pacific and Central Indian Oceans. The CK95 GPTS was obtained by interpolating BMDs between nine radioisotopic dates with a cubic spline, thus ensuring a smooth temporal variation in South Atlantic spreading rates. The South Atlantic BMDs of CK95 have been the key source of information from marine magnetic anomalies for timescale development over more than two decades and have been repeatedly used to construct GPTSs with different age constraints (e.g., Ogg & Smith, 2004;Vandenberghe et al., 2012;Wei, 1995).
The last few decades also saw the development of astrochronology, based on sediment cycles that record Milankovitch periodicities in the Earth's orbit and the orientation of its spin axis (Hinnov & Hilgen, 2012). In the Neogene portion of the most recent GTS12 GPTS , information from marine magnetic anomalies has been mostly replaced by astronomical dating of sedimentary sections with a reliable magnetostratigraphy.
Astronomical dating of older intervals, however, is more challenging, and the marine magnetic anomaly record remains an independent source of information for timescale development. For example, in the Paleogene, the GTS12 GPTS uses a combination of astrochronology (66-47 Ma) and a GPTS based on the CK95 BMDs , and there are significant discrepancies in the Eocene and Paleocene (Vandenberghe et al., 2012). Many astrochronology studies used the CK95 GPTS to initially match sedimentary cycles with astronomical periods (Herbert et al., 1995;Röhl et al., 2003), to decide between alternative tuning options (Röhl et al., 2003), to estimate the duration of hiatuses (Pälike et al., 2001), to provide age constraints to floating timescales (Jovane et al., 2010), and to compare ages and durations of magnetic polarity chrons (Billups et al., 2004;Francescone et al., 2019;Herbert et al., 1995;Husson et al., 2011;Westerhold et al., 2008, Westerhold et al., 2017.

Updating Late Cretaceous-Eocene BMDs and GPTS
The two objectives of our work are (1) to generate a global set of BMDs obtained independently over a number of spreading centers and (2) to construct an updated GPTS that minimizes global spreading rate variations. The global record of magnetic anomalies has grown since it was analyzed by Cande and Kent (1992), and detailed finite plate rotations have been determined over a number of spreading centers besides the South Atlantic. We examine here a large data set of 154 marine magnetic anomaly ship tracks projected onto plate tectonic flow lines computed from finite rotations. The resulting up-to-date global BMDs will inform astrochronology interpretations and constrain GPTS construction.
A drawback of the CK95 GPTS is that the calculated spreading rates in other ocean basins vary more erratically than in the South Atlantic (e.g., Figure 42 of Cande & Kent, 1992). Huestis and Acton (1997) pointed out that there is no reason to expect that spreading rates change smoothly at one particular location while they are more variable at all other ridges. They argued for a "least favoritism" approach where a GPTS is constructed by minimizing the variation in spreading rates in all the magnetic anomaly profiles examined. Independent astronomical dating supports this suggestion: Wilson (1993) and Krijgsman et al. (1999) found that the spreading rates implied by astrochronology in the last~10 Ma were less variable than rates computed from the CK95 GPTS. Astronomical age control steadies the spreading rates (see also Baksi, 1994;Gordon, 1993;Langereis et al., 1994), indicating that minimizing spreading rate fluctuations over several spreading centers will provide a more reliable GPTS.
We concentrate here on the Late Cretaceous-Eocene (chrons C33-C13,~84-33 Ma), which is the interval where the CK95 GPTS is most uncertain. Cande and Kent (1992, p. 13947-13948) recognized that the "largest uncertainties in our new time scale are probably for the anomaly spacings in the late Cretaceous and early Cenozoic. This corresponds to the time interval of the most rapid change in spreading rate […] in the South Atlantic […] when there is the most potential for error." This time interval is also of great interest in paleoenvironmental studies, as it records a long-term warming trend in the Paleocene, a sudden warming event at the Paleocene-Eocene thermal maximum followed by a series of hyperthermals, the greatest Cenozoic warmth in the Early Eocene, and the onset of the cooling trend that resulted in the formation of the Antarctic ice sheets at the Eocene/Oligocene boundary (Zachos et al., 2001(Zachos et al., , 2008. An accurate timescale is critical to advance our understanding of these climatic changes. As noted earlier, astronomical calibration of the Late Cretaceous-Eocene is still in flux, and an updated GPTS would be most useful to constrain timescale development in this critical interval. The GPTS we generate uses only an up-to-date set of radioisotopic dates and includes no astrochronological constraints. Our goal is to provide independent information to astrochronology that can help resolve timescale inconsistencies, for example, around chrons C23 (Vandenberghe et al., 2012;Westerhold et al., 2017) and C28 (Hilgen et al., 2010;Westerhold et al., 2008). A GPTS that fully integrates information from the global marine magnetic anomaly record, radioisotopic dates, and astrochronology will be a future development.
In this paper, we apply a Bayesian formulation to integrate information from sets of diverse, uncertain data (e.g., Agrinier et al., 1999). We first introduce the Monte Carlo procedure used to obtain BMDs that match measured magnetic anomalies and a GPTS that fits radioisotopic dates and minimizes global spreading rate variations. We then describe the fundamental data used to build a GPTS: BMDs obtained over 13 mid-ocean ridge flank regions and radioisotopic dates tied to magnetostratigraphy. A new MQSD20 GPTS is constructed following the approach previously applied by Malinverno et al. (2012) to the M-sequence magnetic anomalies (Late Jurassic-Early Cretaceous). We conclude by comparing the ages and chron durations in MQSD20 with those in existing GPTSs, exploring the impact of the newly obtained BMDs on testing astrochronology interpretations, and describing a global change in spreading rates at 50-45 Ma, when India collided with Eurasia and the Hawaii-Emperor seamount chain changed its orientation.

Nomenclature
We use here the CK95 sequence of chrons to define a series of magnetic polarity crustal blocks that recorded field reversals and do not consider magnetic field excursions such as tiny wiggles (Cande & Kent, 1992). Following general nomenclature (Gee & Kent, 2007;Opdyke & Channell, 1996), C-sequence magnetic anomalies are named 13n, 13r, and so forth, and the corresponding polarity chrons are C13n, C13r, and so forth. Magnetic polarity is denoted by "n" for normal and "r" for reversed. Boundaries of magnetic polarity blocks and chrons are denoted by appending "y" for the young boundary and "o" for the old (e.g., C13ny is the young end of chron C13n). We denote years of age as "a" and years of duration as "yr," with the usual prefixes (e.g., 1 Ma = 1 million years ago). We report uncertainties as one or two standard deviations (1s or 2s, respectively).

MCMC Sampling
To determine the best values and uncertainties of BMDs estimated from ship track magnetic anomalies and of a GPTS that minimizes the global variation of spreading rates, we use here a Markov chain Monte Carlo (MCMC) algorithm (Brooks et al., 2011;Gilks et al., 1996). The algorithm asymptotically generates a sample of model parameter vectors m (e.g., distances to polarity block boundaries) that are distributed as in a target probability density function (PDF). We follow a Bayesian formulation, where the target distribution is a posterior PDF proportional to the product of a prior PDF p(m) that quantifies prior information on the parameters and a likelihood function p(d|m) that quantifies how closely the measured data in a vector d (e.g., magnetic anomalies) are fitted by data predicted by the parameters in m. The posterior PDF is given by Bayes' rule as where k is a normalizing constant. In our application, the parameter vector m also contains "hyperparameters," which are additional variables that affect the solution and are not closely constrained a priori (Gelman et al., 2004;Malinverno & Briggs, 2004). An example is the variance of the misfit between predicted and measured data, which is needed to define the likelihood function but is not known beforehand.
We apply here the Metropolis-Hastings MCMC algorithm (Chib & Greenberg, 1995;Metropolis et al., 1953), which performs a random walk in the space of the parameter vector m. The algorithm first proposes a candidate parameter vector m* that is a small perturbation on the current value of the parameters. The construction of candidate parameter vectors needs to perform a random walk that can go from any point to any other point in the space of the parameters. At each step of the random walk, the proposed candidate is accepted or rejected with a probability that depends on the ratio of the posterior PDF for the current and the candidate parameter vector (so that the normalizing constant k in Equation 1 is irrelevant). By applying this simple accept/reject rule, it can be shown that the resulting Markov chain will asymptotically sample the posterior PDF. This MCMC procedure has been widely used in geophysical inverse problems (e.g., Malinverno, 2002;Malinverno & Leaney, 2005;Piana Agostinetti & Malinverno, 2010;Sambridge & Mosegaard, 2002;Sen & Stoffa, 2013).
We implement here the simple MCMC scheme proposed by Mosegaard and Tarantola (1995). If the random walk that samples candidate model vectors is designed to sample the prior PDF of the parameters, then the probability of accepting the candidate simply depends on the ratio of the likelihoods. The sampling algorithm starts by defining an initial model parameter vector m and calculating its likelihood p(d|m). Each sampling iteration is as follows: 1. Sample a candidate m* following a random walk that samples the prior PDF of the model parameter vector; 2. Compute the candidate likelihood p(d|m*); 3. Accept the candidate with probability (2) 4. If the candidate is accepted, set m = m* and p(d|m) = p(d|m*); if not, m and the likelihood stay the same.
Finally, to ensure that the parameter values output by the MCMC algorithm converge to sampling the posterior PDF, we compare the results of multiple independent sampling chains as suggested by Gelman et al. (2004). The detailed implementation of the MCMC sampling algorithm is described in the supporting information.

Marine Magnetic Anomaly Data
This section describes how magnetic anomaly ship tracks were processed to obtain distances to the boundaries of modeled magnetic polarity blocks. These BMDs are the input to construct a GPTS that minimizes the global variability of spreading rates.
Ship tracks were selected over 13 mid-ocean ridge flank regions that recorded anomalies between 13ny and 34ny at intermediate and fast spreading rates and that had previously determined sequences of finite rotation poles ( Figure 1 and Table 2). The main source of magnetic anomaly data was the NOAA-NCEI archive of trackline geophysical measurements (NOAA National Geophysical Data Center, 1977), supplemented by additional surveys in the Indian Ocean (Yatheesh et al., 2019). We chose tracks that approximately followed the direction of plate motion, did not go over seamounts, and did not cross the fracture zones and ridge axis discontinuities listed in the Global Seafloor Fabric and Magnetic Lineation Data Base (Matthews et al., 2011). Each track was projected on a nearby plate tectonic flow line defined by the finite rotation poles. This projection ensured that the BMDs determined on each projected track were measured along the direction of plate motion.
Initial BMDs in each projected track were determined from the position of key anomaly picks taken from the Global Seafloor Fabric and Magnetic Lineation Data Base (Seton et al., 2014). These key anomaly picks were then interpolated assuming piecewise constant spreading rates between the polarity blocks defined by the CK95 GPTS to obtain an initial set of BMDs. Anomalies predicted in each track by the initial BMDs and an initial anomaly skewness angle (Schouten & McCamy, 1972) were compared to the measured anomalies. The initial BMDs and skewness were then manually adjusted to improve the overall fit to the measured data. While this manually adjusted fit did not have to be exact, the adjusted BMDs and skewness had to result in an approximate match between modeled and measured anomalies for the Monte Carlo algorithm in the next step to be successful (e.g., Figure 2a). During this process, tracks or track segments that did not display an

10.1029/2020JB020034
Journal of Geophysical Research: Solid Earth unambiguous correlation to the overall predicted anomaly pattern were discarded, resulting in a data set of 154 original and projected tracks (location maps are in the supporting information).
The manually adjusted skewness angle and BMDs for each track were then iteratively modified by a MCMC algorithm to obtain a sample of BMDs that fit the magnetic anomaly data. The MCMC algorithm also adjusted a set of nodes that define an interpolated multiplier for the magnetic anomaly amplitudes (Figure 2). Accounting for the variation of anomaly amplitudes along tracks proved necessary to prevent the Monte Carlo algorithm from occasionally sampling unrealistically narrow polarity blocks. If the amplitudes of the predicted anomalies were not adjusted, a narrow, low-amplitude anomaly could only be fitted by making the corresponding polarity block narrower than it should be. The MCMC algorithm was generally successful in sampling values of skewness angle and BMDs that gave a close match between predicted Figure 2. Comparison of measured-predicted magnetic anomalies for starting BMDs (a, b) and best fit BMDs sampled by a MCMC algorithm (c, d). The projected ship track is nbp9604_PAC_a in the PACBAN-PAC ridge flank region. The vertical dashed lines in (a) show the starting BMDs manually adjusted to approximately fit measured anomalies. The MCMC algorithm modifies iteratively the starting BMDs, the anomaly multiplier nodes, and the skewness angle to maximize the fit between measured and predicted anomalies. and measured anomalies (e.g., Figure 2c). In a few instances, it was necessary to manually adjust further the starting skewness and BMDs and repeat the MCMC sampling to obtain a good fit.
The final product of Monte Carlo sampling is the best value (the sample average) and uncertainty (the sample standard deviation) of the BMDs in each projected track. Details on the MCMC procedure and figures with the geographic locations of polarity block boundaries are in the supporting information.
We then obtained summary BMDs over each of the 13 ridge flank regions listed in Table 2 by rescaling and averaging the BMDs estimated in each projected track. Rescaling the BMDs accounts for the systematic variation of local spreading rates as a function of the distance between the track and the plate rotation poles (Figure 3). The rescaling was based on a reference projected track in each ridge flank region. This reference track was either a single long track or multiple tracks that spanned the full range of anomalies recorded over the ridge flank region. The BMDs estimated in each of the other projected tracks were rescaled with a least squares fit so that the average spreading rate in the interval covered by the track was the same as in the matching interval in the reference track. This simple rescaling resulted in a consistent set of BMDs in each ridge flank region (see Figure 3 for an example). The best value for the ridge flank region BMDs was set to the median of the rescaled BMDs on each track. We used the median because it is a central value estimator that is not affected by occasional outliers. The uncertainty of the ridge flank region BMDs was quantified from the standard deviation of the rescaled BMDs for each track. The standard deviation will be increased by outliers, but this is a useful characteristic as it will decrease the influence on the final GPTS of BMDs that contain outliers and are poorly determined.
Finally, we deleted from the summary BMDs polarity blocks narrower than 2 km. Narrow polarity blocks are poorly recorded by magnetic anomalies generated by a source layer whose top is at~4.5 km water depth in 33 Ma crust (Parsons & Sclater, 1977). The width of these narrow blocks is mostly constrained by the distances to adjacent polarity block boundaries and does not provide independent information. The summary BMDs for 13 ridge flank regions are illustrated in Figure 4. Files listing the 154 tracks used here, the skewness angle in each track, the geographic positions of polarity block boundaries and BMDs in the original and projected tracks, and the summary BMDs for each of the 13 ridge flank regions are available in two open access data publications (Malinverno et al., 2019a(Malinverno et al., , 2019b.

Radioisotopic Dates
The radioisotopic dates used here (Table 3) are from Table 28.1 of GTS12, with a few modifications. All radioisotopic dates in Table 3 are 40 Ar/ 39 Ar for a Fish Canyon sanidine age of 28.201 Ma (Kuiper et al., 2008), except for the 206 Pb/ 238 U 55.785 Ma date of C24n.3r. Further details and source references are in Vandenberghe et al. (2012).
GTS12 reports both a radioisotopic dating uncertainty and a stratigraphic position uncertainty expressed as a fraction of the stratigraphic thickness of the polarity chron (assumed here to represent two standard deviations). For GPTS construction, we then assigned to each radioisotopic date a total uncertainty from the sum of the variances due to radioisotopic dating uncertainty and stratigraphic uncertainty. The temporal uncertainty due to stratigraphic uncertainty was calculated as the product of the dimensionless stratigraphic uncertainty times the duration of the respective chron in CK95. The two radioisotopic dates in chron C13r listed in Table 28.1 of GTS12 have been averaged to a single date (first row of Table 3). These dates were obtained in the Massignano quarry stratotype section from biotite-rich clayey layers of possible volcanic origin (Odin et al., 1991). GTS12 used a Fish Canyon sanidine age of 28.201 Ma to recalculate the original Ar/Ar radioisotopic dates as 34.4 ± 0.2 Ma (14.7 m stratigraphic height in Massignano quarry section, 0.4 up from base of C13r, 2s uncertainty) and 35.2 ± 0.2 Ma (12.7 m stratigraphic height, 0.14 up from base). We referred these dates to the stratigraphic framework of a drill core located about 110 m south of the Massignano stratotype section (Lanci et al., 1996). The drill core samples were taken every 12-15 cm, a much more detailed sampling interval than that possible in the Massignano outcrop, where strong weathering makes it difficult to obtain closely spaced, pristine samples. The C13r interval in the drill core was clearly established between stratigraphic depths of 14.2 and 24.8 m ( Figure 6 of Lanci et al., 1996). Using the conversion in Figure 3 of Lanci et al. (1996), the predicted stratigraphic height of C13r in the quarry section is 11.7 to 22.3 m. This estimate is in close agreement with the location of reversely magnetized samples in the outcrop (Figure 9 of Lowrie & Lanci, 1994). Based on the high-resolution magnetostratigraphy in the drill core, the stratigraphic positions of the two dates in C13r become 0.28 and 0.09 up from base.
Using directly the two C13r dates listed in GTS12 with the stated small uncertainties of 0.2 Ma will unduly bias the duration of C13r in the GPTS. To avoid this, we simply averaged the two dates (34.8 Ma) and located the average age midway between the two stratigraphic positions (0.19 up from the base of C13r). We did not change the radioisotopic date uncertainty but increased the 2s stratigraphic position uncertainty (0.05 in GTS12) to half the range of the two original dates, which is (0.28-0.09) / 2 ≈ 0.1.

The MQSD20 GPTS
We obtain a new MQSD20 GPTS with a MCMC algorithm that generates a large ensemble of GPTSs. This ensemble will be asymptotically distributed as in the posterior PDF of Equation 1 and is used to calculate a reference GPTS and quantify its posterior uncertainty. Full details of the MCMC procedure are in the supporting information, and here we illustrate how the sampling proceeds when the GPTS is constrained by different types of information ( Figure 5).  (Table 3). The sampling starts from a perturbed version of the CK95 GPTS and proceeds by iteratively changing chron boundary ages, accepting or rejecting such changes as in the Metropolis acceptance probability of Equation 2. The sampled GPTSs fit the radioisotopic dates within their uncertainty, but the chron boundary ages and durations in intervals between radioisotopic age ties are unconstrained and extremely variable (e.g., between about 66 and 80 Ma).  (Figure 4). In this case, the likelihood function is greater if the GPTS implies a smaller variation of spreading rates (which depend on the BMDs). The variability of the sampled chron ages between radioisotopic age ties is much less than in Figure 5a, as it is constrained by the need to minimize spreading rate fluctuations.  (Swisher et al., 1993) and was omitted. d Stratigraphic uncertainty was not reported in Table 28.1 of Vandenberghe et al. (2012) and was set conservatively as the maximum given elsewhere (0.1).

Journal of Geophysical Research: Solid Earth
Figure 5b also illustrates that the MCMC algorithm does not converge to an optimum GPTS but rather performs a random walk in the space of possible GPTSs. After an initial phase where the random walk moves to the region of high posterior probability (the first~1,500 iterations in Figure 5b), the algorithm samples the posterior PDF, exploring the space of possible GPTS ages that fit the radioisotopic dates and minimize spreading rate variability within the data uncertainties.
The final sample of GPTSs used to derive MQSD20 is constructed by combining the results of 10 independent sampling chains such as that in Figure 5b (see supporting information). The MQSD20 chron boundary ages, chron durations, and respective uncertainties are in Table 4 and Figure 6. The reference GPTS chron boundary ages and chron durations are the average values of the ensemble obtained by MCMC sampling. GPTS uncertainties are quantified from the standard deviations of chron ages and durations in the ensemble. Near age ties, the uncertainties in the GPTS chron boundary ages are smaller than the uncertainties of the  (Table 3). Sampling is constrained to fit only the radioisotopic dates (a) or to fit radioisotopic dates while minimizing spreading rate variability over 13 ridge flank regions (b). Journal of Geophysical Research: Solid Earth radioisotopic dates themselves (Figure 6), reflecting the additional constraints given by spreading rate information.

Journal of Geophysical Research: Solid Earth
The chron boundary ages and chron durations of MQSD20 are compared to those of CK95 (Cande & Kent, 1992 and of GTS12 (Ogg, 2012) in Figure 7. The differences in chron boundary ages with CK95 are due to the use of a different set of BMDs and radioisotopic dates in MQSD20. These differences reach~1 Myr for chrons C24 and earlier (ages ≥ 53 Ma), where they are mostly due to radioisotopic age recalibration.
As MQSD20 used the same radioisotopic age ties as GTS12, the age differences are smaller (~0.5 Myr or less). The Paleogene GTS12 GPTS was based on a spline fit of radioisotopic dates and the CK95 BMDs for chrons C16r-C23r and an astrochronology age model for chrons C24n-C29r (Table 28.3 of Vandenberghe et al., 2012). In contrast, MQSD20 used a different set of BMDs from CK95 and did not include astrochronology information. The largest age difference is for C27n, which is midway within an~5.3 Myr interval Figure 6. MQSD20 GPTS and its 1s uncertainties in chron boundary ages and chron durations. Red error bars show the 1s uncertainties of radioisotopic dates (Table 3).

10.1029/2020JB020034
Journal of Geophysical Research: Solid Earth spanned by two radioisotopic dates (in C26r and C28r). The~0.5 Myr age difference is explained by the use of different data in the two GPTSs: magnetic anomaly information from a global set of BMDs in MQSD20 and astrochronology in GTS12. Finally, chron duration differences are at most 0.3 Myr with respect to both CK95 and GTS12 and are generally within the 2s uncertainty of the chron durations in MQSD20.
The half-spreading rates implied by MQSD20, CK95, and GTS12 for the summary BMDs in the 13 ridge flank regions considered here are compared in Figures S30-S42. The overall variability of spreading rate in each ridge flank region can be quantified by a coefficient of variation (CV) that equals the standard deviation of spreading rate divided by its mean value. Table 5 lists the CVs of spreading rate computed using different GPTSs and the summary BMDs in the 13 ridge flank regions. As MQSD20 minimizes spreading rate variability, the CVs of spreading rate are less than those implied by CK95 and GTS12. The exception is PACVAN-PAC in the North Pacific (Figure 1), where MQSD20 has more variable spreading rates around chrons C17-C18 than CK95 and GTS12 (see Figure S29). The likely reason is that the duration of chrons C17-C18 in CK95 and GTS12 is mostly controlled by magnetic anomaly records in this area of the North Pacific, whereas MQSD20 includes information from other mid-ocean ridge flanks.

Addressing the 50 Ma Discrepancy
As it has been obtained independently from astronomical dating, the MQSD20 GPTS can be used to address conflicting results in astrochronology interpretations. We discuss here as an example the "50 Ma discrepancy" noted by Vandenberghe et al. (2012), which centers on the duration of chron C23n.2n. Whereas this chron lasts 696 Kyr in the CK95 GPTS, it has been estimated to be less than 400 Kyr in several astrochronology studies. From an analysis of the sedimentary record at ODP Site 1258, Westerhold and Röhl (2009) obtained a C23n.2n duration of 379 to 399 Kyr, depending on the astronomical cycle chosen (long eccentricity or precession, respectively). Westerhold et al. (2015) confirmed this interpretation in a study that included other drill sites and listed a duration of 377 Kyr for C23n.2n. From the ODP Site 1263 record, Lauretano et al. (2016) reported an even shorter C23n.2n duration of 295 Kyr in their preferred age model (with an alternative estimate of 395 Kyr). When spreading rates are calculated from the South Atlantic BMDs of CK95, these astronomically determined durations result in spreading rates that are more than twice as fast during chron C23n.2n than in adjacent chrons (Figure 6a of Westerhold et al., 2017). A possible explanation offered for this discrepancy is that the CK95 BMDs may be poorly determined around chron C23 (Westerhold et al., 2015;Westerhold & Röhl, 2009), as the width of the chron C23 block has the largest uncertainties reported in CK95 (17.3% of its width; see Table 4 of Cande & Kent, 1992).
Distances to both the young and old end of C23n.2n were estimated in 85 of the 154 ship tracks examined here. All the ridge flank region BMDs record the young and old end of C23n.2n except for PACFAV-PAC, which only spans chrons C24n and older. When spreading rate variations over all ridge flank regions are considered, the duration of C23n.2n in MQSD20 is 699 ± 180 ka (2s; Table 4). Although we examined a larger magnetic anomaly  Journal of Geophysical Research: Solid Earth data set and used a different set of radioisotopic dates, the duration of C23n.2n we obtain is effectively the same as that in CK95 (696 Kyr). This supports the reliability of the CK95 BMDs and implies that a 300-400 Kyr duration requires a doubling of spreading rates during C23n.2n in multiple mid-ocean ridges, which is implausible. The duration of C23n.2n in MQSD20 agrees with the solution of the 50 Ma discrepancy put forward by Westerhold et al. (2017), who concluded that chron C23n is too short in the magnetostratigraphic interpretation of Site 1258. Their revised astronomical timescale gives a C23n.2n duration of 712 ± 123 Kyr, which is consistent with the MQSD20 results. Figure 8 plots the summary BMDs in each of the 13 ridge flank regions as a function of MQSD20 age. The minimization of global spreading rate fluctuations highlights a major spreading rate change in the interval 50-45 Ma (chrons C22r-C20r). During this time, spreading rates decreased by a factor of 2-3 in the Indian Ocean while they approximately doubled in the North Pacific and in the South Atlantic (Figure 9). These spreading rate changes coincide with a previously noted set of global tectonic events (Francescone et al., 2019;Sutherland et al., 2020;Wessel et al., 2006;Whittaker et al., 2007) that we summarize below.

Spreading Rate Changes and Global Tectonic Events at 50-45 Ma
The prominent decrease in Indian Ocean spreading rates at C22-C20 time is related to the onset of the collision between the India subcontinent and Eurasia (Copley et al., 2010;Molnar & Stock, 2009;Patriat & Achache, 1984). To the west, the contemporaneous spreading rate increase in the South Atlantic confirms a general pattern of spreading rate anticorrelation between the Indian and South Atlantic Ocean observed throughout the period 80-30 Ma (Figure 14 of Cande & Patriat, 2015). Moving eastward in the Indian Ocean, in the interval C22-C20 spreading rates decreased substantially in the Wharton Ridge separating the Indian and Australian plates (Jacob et al., 2014). At the same time, a major Australia-Antarctic plate reorganization took place (Whittaker et al., 2007), and spreading rates in the Southeast Indian Ridge markedly increased immediately after C21 (Figure 2h of Cogné & Humler, 2006).  Figures S17-S29). The slope of these distance-age plots is the half-spreading rate. Red symbols are radioisotopic ages (Table 3)

10.1029/2020JB020034
In the southeast Pacific Ocean, another plate boundary reorganization took place around C22-C21, when the Pacific-Antarctic ridge propagated northward breaking off a large fragment of the Pacific plate that became attached to the Antarctic plate (Cande et al., 1982). In the western Pacific, the onset of subduction-related volcanism in the Izu-Bonin-Mariana and Tonga-Kermadec arcs has been dated to 52-45 Ma (Arculus et al., 2015;Bloomer et al., 1995;Ishizuka et al., 2011).
In the northern Pacific, an about twofold increase in spreading rate at C22-C20 time has been noted by others (Barckhausen et al., 2013;Wright et al., 2015). This spreading rate increase coincides with the prominent change in orientation in the Hawaiian-Emperor seamount chain (Hawaiian-Emperor bend [HEB]), whose date has been recently updated to~47 Ma (O'Connor et al., 2013;Torsvik et al., 2017;Wessel et al., 2006). The HEB was originally explained by a change in absolute motion of the Pacific plate over a fixed hot spot, but later on, several authors argued that it resulted from a slowing southward motion of the Hawaiian hot spot (Norton, 1995;Tarduno et al., 2003;Wright et al., 2015). Torsvik et al. (2017), however, recently concluded that both a southward shift of the Hawaiian hot spot and a change in Pacific plate motion direction are necessary to explain the HEB. The substantial increase in northern Pacific spreading rates at the same time of the HEB strongly suggests a connection. We conjecture that even if the direction of absolute Pacific plate motion did not change, a substantial acceleration in Pacific spreading rate at the time of the HEB over a southward drifting Hawaiian hot spot may have turned the orientation of the Hawaiian-Emperor seamount chain closer to an E-W direction.
Contemporaneous worldwide spreading rate changes and plate boundary reorganizations suggest a global connection. As plate motion changes are likely controlled by changes in plate boundary forces (Gordon et al., 1978;Iaffaldano, 2014), we speculate that the effects of the India-Eurasia collision may have propagated throughout the global plate tectonic system. A comprehensive explanation of the connection between the 50-45 Ma events is beyond the scope of this paper; however, we stress the importance of a timescale

10.1029/2020JB020034
Journal of Geophysical Research: Solid Earth constrained by the global magnetic anomaly record to reliably reconstruct spreading rate changes and correlate major plate tectonic events.

Conclusions
We estimated here a new set of magnetic polarity BMDs spanning the chron C33-C13 interval in 154 ship tracks projected onto plate tectonic flow lines. The ship track data were assembled in summary BMDs over 13 ridge flank regions in the southern and northern Pacific, the southern Atlantic, and the Indian Ocean. This new set of BMDs extends substantially the South Atlantic-based distances originally compiled by Cande and Kent (1992). We used these BMDs to construct a Late Cretaceous-Eocene MQSD20 GPTS that minimizes the variability of spreading rates over all ridge flank regions and fits an up-to-date set of radioisotopic dates. At 50-45 Ma, MQSD20 shows a marked spreading rate decrease in the Indian Ocean and a contemporaneous increase in the South Atlantic and Northern Pacific. This spreading rate change coincides with the India-Eurasia collision and with the bend in the Hawaii-Emperor seamount chain.
The MQSD20 GPTS deliberately did not include astrochronology constraints in order to provide an independent source of information to check sediment cycle interpretations. The next step forward in GPTS construction will be to directly incorporate in the timescale information from astrochronology. Such an integration procedure will improve the usual approach, which is to build the timescale on the best data set that is typically taken to supersede other sources of information that are deemed less accurate. For example, future timescale development is often viewed as astronomical dating replacing GPTSs based on marine magnetic anomalies (e.g., Gradstein, 2012 p. 13;Hilgen et al., 2012 p. 947). A GPTS constructed on the basis of magnetic anomalies from multiple spreading centers, rather than from a single mid-ocean ridge, points to a better approach where diverse data sources are combined rather than selectively discarded. In this view, the GPTS is the result of an integration of astrochronology, radioisotopic dates, and magnetic anomaly data, where each piece of information is weighted by a measure of its uncertainty (e.g., Malinverno et al., 2012). The global set of BMDs and the Monte Carlo methods presented here provide the basis for this advance in timescale construction.

Data Availability Statement
Data supporting the conclusions of this study are in two open access data publications that list the 154 ship tracks used here, the geographic positions and distances to polarity block model boundaries (BMDs) in each of the original and projected ship tracks, and the summary BMDs for each of the 13 ridge flank regions (Malinverno et al., 2019a(Malinverno et al., , 2019b.

Acknowledgments
This study was supported by Award OCE-1535937 of the U.S. National Science Foundation. Thanks to two anonymous reviewers whose comments improved the paper. J. D. thanks G. C. Bhattacharya and his colleagues who, through Indo-French project CEFIPRA 3307-1, made some of the Indian Ocean data used in this study available in a common publication (Yatheesh et al., 2019). This is LDEO contribution number 8429 and IPGP contribution 4154.