Steady increase in geomagnetic reversal frequency after the Cretaceous Normal Superchron

The Earth's core is constantly and eciently cooled by mantle convection. The heat ux transferred from the core to the mantle through the core-mantle boundary (CMB) is critical for understanding the dynamics of solid Earth. Although it is dicult to estimate the CMB heat ux, its history could be reconstructed from geomagnetic reversal frequency. However, overlooked short geomagnetic reversals may exist in the geomagnetic polarity time scale (GPTS), which affects the estimation of the heat ux history. Here, we report four new high-precision 40 Ar/ 39 Ar ages of the Oligocene Ethiopian traps. The traps may contain undiscovered reversals in marine magnetic anomaly. Based on the ages, we identied new reversals in Chron C12n, which was not found in marine magnetic anomalies. Our non-parametric analysis of GPTS suggests four potential periods of missing geomagnetic reversals, which correspond to long polarity intervals in GPTS. We found that C12n correspond to one of the periods. This indicates that several undetected reversals may exist within or near the edge of long polarity intervals after the Cretaceous Normal Superchron (prolonged stable polarity period). Considering the undetected reversals, we conclude that the CMB heat ux increased more slowly and monotonically after the Superchron than that ever estimated.

We conducted four high-precision 40 Ar/ 39 Ar datings 12 (see Methods). We used the same samples as previous palaeomagnetic studies 6,13 . Before the dating experiments, we selected samples by inspecting thin sections (Fig. S1) to avoid contamination from secondary minerals. Samples A1-61 and A2-162 are from two reversely magnetized lava ows belonging to the R4 magnetozone. Samples A2-482 and A2-604 are from the normally magnetized lava ows belonging to N3 and N1 magnetozones, respectively. Four samples were dated by stepwise heating analysis (Table 1, S1; Fig. 2). Two samples (A1-61 and A2-162) yielded well-de ned age plateaus comprised of 82.7% and 82.1% of released gas, respectively (Table 1, S1; Fig. 2). The inverse isochron ages for the two age spectra are identical to the weighted average ages of plateau-forming steps within 2σ error (Table 1; Fig. S2). These data indicate that the two plateau ages (29.63 ± 0.14 Ma for sample A1-61 and 30.02 ± 0.22 Ma for sample A2-162) (2σ, respectively) are reliable eruption ages of the basalts. A2-482 showed a partially disturbed age spectrum comprised of 44.9% of released gas and did not have a plateau in a strict mean (Table 1, S1; Fig. 2). However, the inverse isochron age of sample A2-482 ( Fig. S2) is identical to the weighted average age of the middle-temperature steps and the 40 Ar/ 36 Ar intercept of the inverse isochron is consistent with the atmospheric ratio within 2σ error. Thus, the middle-temperature steps are regarded as plateau-forming steps, and the weighted average age (30.87 ± 0.22 Ma, 2σ) from these steps is interpreted as a reliable eruption age. Sample A2-604 gave an age spectrum comprised of 43.5% of the released gas and did not have a plateau in a strict mean either (Table 1, S1; Fig. 2). The age spectrum includes some apparent disturbance in low-temperature steps due to 39 Ar recoil in possible secondary phases. However, middle to high-temperature steps overlap with the plateau-forming steps of sample A2-162, which has a similar Ca/K plot with sample A2-604 (Fig. 2). These steps gave a weighted average age of 31.02 ± 0.24 Ma (2σ), which we interpret as plateau-forming steps and a reliable eruption age. Also, the four plateau ages are consistent with the altitudes in the Lima-Limo section (Fig. 3), which supports the reliability of the ages.
Age correlation with the geomagnetic polarity time scale The new ages enable a unique correlation of the magnetostratigraphy to the GPTS2020 3 . GPTS2020 is the latest version of the geomagnetic polarity time scale. This GPTS provides ages of geomagnetic reversal boundaries based on marine magnetic anomalies since middle Mesozoic, which is constrained by astronomical tuning. Our four 40 Ar/ 39 Ar ages are calibrated using the FCs standard with the age of 28.201 Ma 14 and the 40 K decay constant 15 , which are used in GPTS2020. Our ages are consistent with the previously reported high-precise ages from sanidine 9 (Fig. S3). The age of sample A1-61 from the R4 magnetozone (29.82 ± 0.14 Ma, 2σ) corresponds to late Chron C11n.2n of GPTS2020 (Fig. 3). The correlation is inconsistent with the measured reversed polarity of the samples. However, we interpret that the sample A1-61 belongs to Chron C11r because the age of sample A1-61 is closer to that of the top of Chron C11r than the bottom of Chron C11n.1r (Fig. 3). On the other hand, the 2σ error range of the age of sample A2-162 (30.21 ± 0.22 Ma, 2σ) is close to the interval of Chron C11r (Fig. 3).
The three putative excursions in the R4 magnetozone and the one in the R1 magnetozone have been previously reported 6 within Chrons C11r and C12r, respectively ( Fig. 3) based on the VGP cut-off method 16 . These intermediate palaeomagnetic directions may be correlative with Cryptochron C11r-1 and one of eight Cryptochrons (C12r-8 to C12r-1) in Chron C12r, respectively (Cryptochrons were de ned by Cande and Kent (17)). We tentatively assigned the age 29.97 Ma to the top of the magnetostratigraphy, which is the youngest age of Chron C11r, and 31.29 Ma to the bottom, which is the age of Cryptochron C12r-1 interpolated in GPTS2020 (Fig. 1b).
Based on our age constraints, four new reversals bounding the magnetozones from N1 to R3 are included in Chron C12n. We call these four brief reversals in Chron C12n "Lima-Limo reversals", which are not identi ed from the marine magnetic anomalies 17 . A low of the relative palaeomagnetic intensity in the early C12n observed in Oligocene marine sediments 18,19 may re ect the Lima-Limo reversals and suggest that this is a global event.
Geomagnetic reversal frequency estimation and missing reversals Geomagnetic reversal is a stochastic phenomenon called Poisson process 20 . Therefore, the Kernel Density Estimation (KDE), which is a sum of probability functions, is suitable to express the change in the geomagnetic reversal frequency. We used the xed and locally adaptive KDE 21 to estimate reversal frequency changes (see Methods). In general, the xed KDE is expressed as the sum of Gaussian functions and depends on the bandwidth, which determines the standard deviation of the averages of the Gaussian functions. The bandwidth is optimised across the entire data set. As a result, the bandwidth may partly be too wide to extract short-wavelength features, while apparent features may be generated when the bandwidth is too narrow. On the other hand, the locally adaptive KDE optimises bandwidth in accordance with data by minimising the difference between true unknown rates and estimated rates assuming the Poisson processes. Therefore, the locally adaptive KDE can reveal true features of the reversal frequency changes. The curve of the xed KDE has four reversal frequency minima after the Cretaceous Normal Superchron (CNS) (Fig. 4a, blue line, Table S2). However, a single peak and no minima (Fig. 4b, blue line, Table S2) is a true feature according to the locally adaptive KDE. This variation between 84 and 24 Ma is consistent with that estimated by a previous study 11 (Fig. 4b, gray line). Therefore, the short-wavelength changes in the reversal frequency expressed by the xed KDE could be artefacts. These apparent minima correspond to polarity chrons longer than 0.9 Myr in GPTS2020. However, true reversal frequency based on the locally adaptive KDE shows monotonical changes even during such periods. This discrepancy can be resolved if undiscovered geomagnetic reversals exist within or near the edges of these long chrons.
Next, we examined the change of reversal frequency curves when the Lima-Limo reversals were added to GPTS2020. For the xed KDE, after adding the Lima-Limo reversals, the decrease of the reversal frequency curve at around 30 Ma became smaller (Fig. 4a, red line, Table S2) and closer to that of the locally adaptive KDE mentioned bellow. These observations means that the Lima-Limo reversals compensate the part of reversals undetected in marine magnetic anomalies. For the locally adaptive KDE, the addition of the Lima-Limo reversals resulted in a slight decrease in the curvature of the reversal frequency curve from 50 to 25 Ma (Fig. 4b, red line, Table S2). The rate of the reversal frequency increase is 0.12 Myr -1 and 0.10 Myr -1 from 40 to 25 Ma before and after adding the Lima-Limo reversals, respectively, which implies that the increasing rate of the reversal frequency after CNS is gradual and close to constant. The bend of the reversal frequency curve around 45 Ma (Fig. 4b, red line) may be an artefact due to the undetected geomagnetic reversals.

CMB heat ux change after CNS
The slow and monotonous increase in the reversal frequency after CNS suggests continuous and steady changes in the CMB heat ux after CNS. This is because the reversal frequency re ects the variations of the CMB heat ux average according to the numerical geodynamo simulations 2 . A previous study linked reversal frequency to the growth and collapse of superplumes 22 . They interpreted the decrease in the reversal frequency as a collapse of a superplume. The collapse would trigger the production of hot plumes, and which leads the emplacement of the Large Igneous Provinces (LIPs) with a time lag of 30-60 Myr required for the hot plume to rise 22 . In other words, the timing of the emplacements of LIPs corresponds with the peaks of the reversal frequency with the time lag. However, over the past 160 Myr, our analysis shows that the reversal frequency had only three peaks, while LIPs have occurred nine times (Fig. 4b, c). Therefore, the change in the reversal frequency is alternatively interpreted to be linked with constant subduction and the fall of slabs onto the CMB, and not with sporadic superplume and activity of the LIPs. Based on this interpretation, the slow and monotonic increase in reversal frequency from 84 to 9.6 Ma suggests a gradual increase in the number of slabs reaching the CMB during this period. When further undetected geomagnetic reversals after CNS are discovered, the change in the reversal frequency during this period will become more gradual. Considering the observed time lag of 120 Myr between the subduction ux and the reversal frequency 23 , the reversal frequency after CNS would therefore be a constraint for plate reconstruction models from ~220 to ~120 Ma. Ages of the basalt samples from the Lima-Limo section were determined using the 40 Ar/ 39 Ar dating instrument at the Geological Survey of Japan/AIST. Details of the procedures are described in Ishizuka et al. (13). 20-25 mg of phenocryst-free groundmass, crushed and sieved to 250-500 μm in size, was analyzed using a step-wise heating procedure. The samples were treated in 6N HCl for 30 minutes and then 4N HNO 3 for 30 minutes at 95•C with stirring to remove any alteration products (clays and carbonates) present in interstitial spaces. After this treatment, samples were examined under a microscope, and microphenocrysts of clinopyroxene were removed. This procedure effectively separated and concentrated fresh plagioclase in the groundmass and microphenocrysts. Sample irradiation was done at the reactor of the Institute for Integrated Radiation and Nuclear Science, Kyoto University for four hours. Sanidine separated from the Fish Canyon Tuff (FC3) was used for the ux monitor and assigned an age of 27.5 Ma, which has been determined against our primary standard for our K-Ar laboratory 35 , Sori biotite, whose age is 91.2 Ma.
A CO 2 laser heating system (NEWWAVE MIR10-30) was used at continuous wave mode for sample heating. A faceted lens was used to obtain a 3.2 mm-diameter beam with homogenous energy distribution to ensure uniform heating of the samples during stepwise heating analysis. Argon isotopes were measured on a IsotopX NGX noble gas mass spectrometer tted with a Hamamatsu Photonics R4146 secondary electron multiplier in a peak-jumping mode.
Correction for interfering isotopes was achieved by analyses of CaF 2 and KFeSiO 4 glasses irradiated with the samples. The blank of the system including the mass spectrometer and the extraction line was 2.9 × 10−14 ml STP for 36 Ar, 1.4 × 10−13 ml STP for 37 Ar, 1.0 × 10−14 ml STP for 38 Ar, 1.2 × 10−14 ml STP for 39 Ar and 1.9 × 10−12 ml STP for 40 Ar. The blank analysis was done every two or three step analyses. All errors for 40 Ar/ 39 Ar results are reported at one standard deviation. Errors for ages include analytical uncertainties for Ar isotope analysis, correction for interfering isotopes and J value estimation. An error of 0.5% was assigned to J values as a pooled estimate during the course of this study. Results of Ar isotopic analyses and correction factors for interfering isotopes are presented in supplementary data (Table S1).
Plateau ages were calculated as weighted means of ages of plateau-forming steps, where each age was weighted by the inverse of its variance. The age plateaus were determined following the de nition 36 . Inverse isochrons were calculated using York's least-squares t 37 , which accommodates errors in both ratios and correlations of errors.

Geomagnetic reversal frequency calculation
The relative geomagnetic reversal frequency variations were computed in Python 3.7.3 using the sskernel.py function for the xed KDE and the ssvkernel.py function for the locally adaptive KDE. The methods are the Kernel Bandwidth optimization approach 18 with a Gaussian window function (https://pypi.org/project/adaptivekde/). The relative geomagnetic reversal frequency was converted to an absolute value using the average reversal frequency over the past 154.9 Ma calculated from GPTS2020. The only assumption of both methods is that the data are generated from a Poisson process. On the locally adaptive KDE, optimized bandwidths are determined locally based on minimization of the mean integrated square error between the unknown underlying rate and the estimated rate.    (3) for this study. When multiple age correlations were proposed in the previous studies, only the correlations preferred by the original authors in each study is presented here. We assumed that the Lima-Limo reversals are evenly distributed over the early half of Chron C12n.   (a) Reversal frequency curves since 154.9 Ma for GPTS2020 3 based on the xed kernel density estimation, before and after adding the newly found four reversals (the Lima-Limo reversals) to Chron C12n based on GPTS2020 (solid blue and red line, respectively). The bandwidth of the xed KDE before and after adding the new reversals are 2.551 and 2.619, respectively. Arrows indicate minima of the reversal frequency. The normal/reverse polarities of GPTS2020 are shown as black/white bars, which do not include the Lima-Limo reversals. The dashed black line shows the timing of the Lima-Limo reversals in Chron C12n. Gray area indicates intervals when the curve has changed signi cantly by adding the Lima-Limo reversals. We assumed that the Lima-Limo reversals are evenly distributed over the early half of Chron C12n. (b) Reversal frequency curves based on locally adaptive kernel density estimation, before and after adding the Lima-Limo reversals to Chron C12n (solid blue and