Cenozoic sea-level and cryospheric evolution from deep-sea geochemical and continental margin records

Time scale and δ¹⁸O_benthic records

We use the geological time scale (GTS) 2012 (85) and its recent updates (86) throughout. We compiled published δ¹⁸O_benthic records from 0 to 65.86 Ma (where the Cretaceous/Paleogene boundary on the GTS is at 66.05 Ma) and constructed a spliced record based exclusively on deep Pacific benthic foraminiferal δ¹⁸O records (Fig. 1 and fig. S2). With publication of recent δ¹⁸O_benthic records (especially Shatsky Rise site 1209; Fig. 1 and fig. S2) (41), we were able to restrict our compilation to records from the deep (>2 km paleodepth) Pacific Ocean to minimize local or regional effects, despite minor δ¹⁸O_benthic variations that can affect records. All records have been astronomically tuned to the GTS by the publications cited. Site 846 (eastern equatorial Pacific, record 0.135 to 5.32 Ma) was tuned to a 65°N Milankovitch target with an uncertainty of less than a precession cycle (<10 ka) (45); older records were tuned using the 405 long eccentricity cycle, which is robust not only for the last 50 Ma (87) but also for much of the Phanerozoic. However, the time resolution could be worse than 405 ka if a cycle is missed or misidentified as may be possible for the Paleocene. Backstripped records have larger age uncertainties of ±0.5 Ma (4, 23–25).

The efforts of the IODP paleoceanographic community are evident in the resolution of the δ¹⁸O_benthic data that is excellent in the Neogene, nominal in the Oligocene, and nascent in the Paleocene-Eocene. Average sample resolution is 0.5 ka in the Late Pleistocene V19-30 record extending back to 133 ka (46); 2.5 ka in the site 846 record extending back to the base of the Pliocene [5.3 Ma; resolving precession, tilt, and eccentricity; data compiled and tuned in (45)]; 2 to 5 ka in the Miocene (resolving precession, tilt, and eccentricity) at sites 1146 (42), 1337 (43), and 1338 (44); and <6 ka in the Oligocene to Early Miocene (resolving tilt and eccentricity) at site 1218 (36). Sample resolution is lower in the late Middle to Late Eocene (58 ka, resolving just eccentricity) at site 1218 (36) and middle Middle Eocene at site 1209 (19 ka, resolving just eccentricity) (88). Paleocene to early Middle Eocene records are better resolved (9.9 ka, potentially resolving tilt and eccentricity) at site 1209 (41). There is a ~300-ka gap in the splice at 8.03 to 8.32 Ma and a 143-ka gap from 20.03 to 20.18 Ma.

All records were adjusted to the Cibicidoides–Planulina values (δ¹⁸O_Cibicidoides) using published (89) calibrations. V19-30 and site 865 data were obtained from Uvigerina (Cibicidoides = δ¹⁸O_Uvigerina – 0.64‰). Data from the middle Middle Eocene at site 1209 are from both Nuttallides truempyi and Oridorsalis umbonatus (88) that were corrected [Cibicidoides = (δ¹⁸O_Nutallides + 0.1)/0.89‰ =δ¹⁸O_Oridorsalis – 0.28‰]. All other records were run using Cibicidoides–Planulina (mostly Cibicidoides mundulus/praemundulus or Planulina wuellerstorfi). The δ¹⁸O_benthic data used in Fig. 1A, temperature [equation 7A in (33)], the δ¹⁸O_benthic-Mg/Ca–based sea-level estimate, and other sea-level curves are given in supplementary tables.

Sea-level estimate

We provide a sea-level estimate using our spliced δ¹⁸O_benthic record (Fig. 1A). We used the 2-Ma smoothed paleotemperatures (Fig. 1) of Cramer et al. (33), who used a low-pass filter that passes >80% of the amplitude for frequencies <0.5/Ma (wavelength, >2 Ma), ramping down to <20% of the amplitude for frequencies >1.25/Ma (wavelength, <0.8 Ma). We used equation 7b in (33) and an updated, simplified paleotemperature equation for benthic foraminifera (33, 90)

(2)to solve for δ¹⁸O_seawater, yielding a δ¹⁸O_benthic-temperature calibration of −0.21‰/°C appropriate for deep-sea temperatures (<13°C). We assume that shorter-term (<2 Ma) temperature changes comprise ~20% of the δ¹⁸O_seawater changes. Our GMSL estimate accounts for long-term temperature changes, although it does not remove glacial-interglacial temperature variability that is significant at least during peak eccentricity-scale interglacials (67). The resultant δ¹⁸O_seawater estimate was scaled to GMSL changes (Fig. 1) using a revised δ¹⁸O_seawater sea-level calibration of 0.13‰/10 m (91) [see (67) for discussion of uncertainties in this calibration]. Using this calibration does not account for changes in δ¹⁸O_ice, and the calibration of 0.13‰/10 m requires a δ¹⁸O_ice of −50‰. The δ¹⁸O_ice composition likely varied from −56.5‰ in the modern EAIS, −35‰ in the modern GIS, and −40‰ inferred for LIS (91). Our scaling using 0.13‰/10 m is necessarily an approximation given the end members of ice sheets and leads to errors as high as 30 m on the largest changes (Δδ¹⁸O_seawater, 1.2‰) and 16 m on typical changes (Δδ¹⁸O_seawater, 0.5‰). Changes in the calibration with evolution of ice sheet size should be modeled, but the calibration provides a useful first approximation considering the overall large errors in the method (±10 to ±20 m) (Fig. 1) (67, 68). Table S1 provides site, age (GTS2012), δ¹⁸O_benthic values corrected to Cibicidoides, δ¹⁸O_seawater, and sea level relative to modern and published backstripped sea-level estimates in (22) adjusted to GTS2012.

Our initial sea-level curve largely yielded realistic variations except for three intervals: Pleistocene interglacials (MIS 5e, 7, 9, 11, etc.), where sea levels are too high (approaching ice-free estimates) due to unrealistically cold input temperatures; the EOT, which initially yielded over 100-m sea-level fall; and the Paleocene, where sea levels appear to be too low (< modern) due to warm input temperatures. The long-term Pleistocene temperatures in (30) are reasonable (~0.4°C), but it has been long known that the Pleistocene deep sea follows a temperature hysteresis curve with a dominant cold mode (<1°C) and shifts to a warm mode (temperature ~1 to 2°C) during peak 100-ka scale interglacials (92). We iteratively fit the last interglacial cycle to known sea level during MIS5e (7 to 9 m) (93) and applied these temperatures (1.8°C) to major Middle to Late Pleistocene peak interglacials, tapering the temperature from the long-term estimates for the peak interglacials using a Gaussian filter. Our corrected sea-level curve for the Pleistocene agrees remarkably well with other sea-level estimates including those of Waelbroeck et al. (94), who used bottom-water temperature-corrected δ¹⁸O_benthic records to scale sea level and that derived from Red Sea (Figs. 3 and 4) (95).

It has been known since the pioneering studies of Lear et al. (31) that carbonate ion changes associated with the EOT and the drop in the CCD overprint the temperature record, yielding apparent warming across the event and unrealistic sea-level amplitudes of >100 m. Pusz et al. (34) empirically corrected for the carbonate ion change, and we used their calibration to remove an apparent warming effect of ~1.5°C; we applied their empirical correction to the sea-level curve, reducing the amplitude by 28 m from 34.17 to 34.30 Ma.

Cramer et al. (33) also noted that Paleocene deep-water temperatures obtained from Mg/Ca (Fig. 1) are biased and that this discrepancy can only be reasonably explained by changing Mg/Ca_seawater, assuming a high Mg/Ca_benthic sensitivity, or inaccuracies in the Mg/Ca_benthic ratio due to species or diagenetic effects. In the absence of independent temperature constraints, we assumed Paleocene bottom water temperatures older than 58 Ma of 9°C, linearly interpolated temperatures between 55.73 (15°C) and 58 Ma (9°C) to the bottom of the record. This yields maximum sea levels of ~60 to 70 m above present for the interval 58 to 66 Ma, consistent with an ice-free world. We note that the temperature, hence our sea-level record, older than ~48 Ma (early Middle Eocene) is suspect and requires additional evaluation.

To compare with continental margin sequences and focus on Ma-scale variability, we interpolated the records to constant sampling interval and filtered the scaled δ¹⁸O-Mg/Ca sea-level record with a Gaussian convolution filter, removing periods shorter than 490 ka (Fig. 1C). Continental margin sedimentation acts as a low-pass filter, generally preserving the Ma-scale cyclicity and aliasing higher-order variability. During intervals of high sedimentation rates (e.g., >10 cm/ka), higher-order (405 ka, quasi–100-ka periods) sea-level cycles can be preserved, but in general, comparison of Ma-scale variability (Fig. 1C) is appropriate.

The continental margin record is derived from backstripping onshore New Jersey and Delaware Coastal Plain records (4, 23, 25) as recently updated by (23) to GTS 2012/2016. Backstripping provides an estimate of GMSL and nonthermal tectonism with errors of ~±10 m (4, 17, 23, 24). However, Kominz et al. (23) showed that Miocene onshore NJ backstripped record presented here deviates from offshore NJ backstripped records due to nonthermal component attributable to MDT (fig. S5). These changes in MDT introduce significant uncertainties, with effects on the order of 10 to 50 m at rates of ~10 m/Ma (27). In our analysis, we note that the onshore backstripped record overlays the δ¹⁸O-Mg/Ca–based estimate younger than 43 Ma (Fig. 1C), but the records appear offset from each other by ~50 m from 43 to 66 Ma (fig. S7B). This may be due to several effects: (i) MDT effects caused excess Eocene subsidence in the backstripped sites, leading to higher Eocene sea levels; this amplitude of offset of ~50 m is predicted by models of MDT (7); (ii) overestimates of water depth in the deep water Eocene sections (e.g., assignments of paleodepth to 150 m could be shifted to 100 m to explain the discrepancy that is within the ±50-m error bar); or (iii) long-term GMSL fall of 50 m during the Eocene ascribed to various non–ice-related processes. Because the focus of this paper is to understand cryospheric evolution and particularly Ma-scale sea-level changes, we subtracted 50 m from the (23) record from 43 to 66 M, as predicted by MDT models (7), to register with the GMSL estimate derived from δ¹⁸O-Mg/Ca.

Stated estimates of uncertainty of sea-level estimates vary among published studies but vary from ±10 to ±20 m depending on method (backstripping or δ¹⁸O-Mg/Ca) and author [e.g., (67, 68)]. Backstripping analyses generally quote ±10 m RSL uncertainties, although uncertainties can be greater considering MDT and other effects (4, 17, 23, 24). Two recent studies have considered errors on comparisons of Pliocene sea-level backstripping and δ¹⁸O-Mg/Ca. Miller et al. (68) propagated the stated uncertainty for the individual measurements for both methods and obtained error estimates of ±10 m or better (1σ) and ± 20 m or better (2σ). Raymo et al. (67) considered all possible errors including diagenesis and concluded that biases of ±20 m are possible. Here, we consider that errors (1σ) are at least ±10 m but less than ±20 m on an individual Ma-scale sea-level variation.

Visual comparisons between δ¹⁸O-Mg/Ca–based and backstripped estimates are compelling (Fig. 1), and here, we provide additional statistical comparisons. We applied a Bayesian analysis with Gaussian process priors (96) to estimate sea-level variations recorded by both the δ¹⁸O-Mg/Ca–based and backstripped sea-level estimates, a_i (Eq. 3) and b_j (Eq. 4), respectively. We modeled the recorded sea-level variations from the δ¹⁸O-Mg/Ca–based method, g(t), and the backstripped records from the NJ margin, h(t), separately through time, t, as the combination of component processes with a nonlinear term and a linear term shared between them [m(t) and l(t); Eqs. 5 and 6]

(3)

(4)

(5)

(6)

The observational errors (

and

) were treated as having SDs of ±10 m for the δ¹⁸O-Mg/Ca–based GMSL estimate and ±15 m for the backstripped NJ record. The component processes are each modeled with a distinct mean-zero Gaussian process prior, and they are interpreted as capturing discrete physical processes. In this case, the linear terms, l(t) and l_h(t), capture the long-term variations in recorded sea level within each record. The sea-level variation captured by the shared l(t) component is interpreted as reflecting long-term ice-volume changes, the largest source of sea-level change shared between the two records. The linear term unique to h(t), l_h(t), is interpreted to include tectonic influences, as well as any other sources of long-term regional trends. The modeled process represented by the shared nonlinear term, m(t), captures variations largely attributable to the ice-volume signal. The nonlinear process, Δ(t), is unique to h(t), and it represents the differences between g(t) and h(t) correlated at an approximately million-year scale. The white-noise terms [w_g(t) and w_h(t)] capture high-frequency variability and measurement error beyond that accounted for in the nominal error estimates. The constant terms (c₁ and c₂) capture differences in time-average sea-level estimates at each site.

The hyperparameters defining prior expectations for the amplitude and scale of each of these component functions were estimated by sampling distributions of hyperparameter values that maximize the likelihood of the joint model of g(t) and h(t) given the two datasets (a and b) using a Markov chain Monte Carlo (MCMC) method. The posterior distributions of the optimized parameters (table S2) were then used to generate posterior estimates of g(t) and h(t) (fig. S7, A and B) and their constituent processes (fig. S7, C to E). The objective was to assess the differences between g(t) and h(t). A useful representation of the difference was generated by calculating an estimate for and the uncertainties of h(t) − g(t) (fig. S7F), which represents the difference between the two curves, and the associated uncertainties, after accounting for the correlations between them.

The analysis indicates that sea level declined from an average of 28 ± 6.8 m between 42 and 40 Ma to −0.5 ± 2.7 m between 12 and 10 Ma, with the shared linear trend component, l(t), indicating a long-term rate of decrease of 1.0 ± 0.5 m/Ma. The two records differ significantly for portions of the record (fig. S7F). Most notably, the backstripped record indicates higher sea-level values for long stretches of time from the middle Eocene through most of the Oligocene. This offset is demonstrated through the linear trend exclusive to h(t), l_h(t). The offset is 3.1 ± 4.8 m at 10 Ma and is 18 ± 20 m at 42 Ma, increasing backwards through time at a rate of 0.5 ± 0.6 m/Ma (fig. S7E). The between-record noise associated with errors correlated on an approximately million-year scale, represented by Δ(t) (fig. S7D), has an SD of approximately 14 m. This noise is the difference between the shared record, or g(t), and h(t) with the offset, l_h(t), accounted for. Although it is not possible to identify which of the two datasets is the source of this noise, the 14-m variance-based metric of difference derived from Δ(t) falls within the range of uncertainty estimates for the δ¹⁸O-Mg/Ca–based made by Miller et al. (68) and Raymo et al. (67).

Time series analyses

We generated two continuous wavelet transform (CWT) spectrograms using the δ¹⁸O data (Fig. 2) and the δ¹⁸O-Mg/Ca GMSL estimate (Fig. 2) data to examine the variation in spectral power associated with certain frequencies through time within these records. The two datasets were preprocessed by first removing the mean and long-term trends. A single linear function does not closely fit the entirety of either dataset and cannot be used to adequately remove the means and trends from them. Therefore, we used a LOESS regression algorithm with a 20-Ma window to establish the long-term trends through each series and subtracted the regression result from the original data. As a result, all relevant low frequency variations were maintained (conservatively, all frequencies higher than or equal to one that corresponds to a 5-Ma periodicity), and the entirety of each record had a consistent local mean at a 10- to 20-Ma scale. We then interpolated these data to generate records with a consistent 3-ka time step. Using the preprocessed δ¹⁸O data and δ¹⁸O-Mg/Ca sea-level estimates, we applied the CWT, which can be defined as the convolution of the time series data with a series of wavelets scaled to respond with amplitude variations that correspond with given frequencies in a Fourier transform. We used the MATLAB wavelet toolbox for wavelet analysis. The resulting CWT amplitudes at given coordinates in frequency and time were plotted using a color gradient. As a baseline for comparison, we also applied the CWT method to both the La2004 solution (97) for solar insulation at 65°N on June 21 and the La2010 solution (86) for the eccentricity of the Earth-Moon barycenter and plotted the outputs alongside the δ¹⁸O and δ¹⁸O-Mg/Ca GMSL spectrograms.

To evaluate the shortest periods resolvable for a given time, we plotted the “local” Nyquist frequency of the δ¹⁸O data and δ¹⁸O-Mg/Ca sea-level estimates over the CWT spectrograms (white line, Fig. 2), connecting to the “cone of influence” that demarcates the time window of usable amplitude information at each frequency. The local Nyquist frequency corresponds to a time-variable boundary marking the highest frequency that amplitude information should be possible to recover considering the original data density as calculated by (i) counting the number of δ¹⁸O data points within 500 ka of a given point in time, (ii) dividing range of time those points cover by the number of points counted to obtain the average data spacing for that window, and (iii) taking the inverse of two times the data spacing in time. This process was repeated in 1000-year steps through the entire δ¹⁸O record and δ¹⁸O-Mg/Ca sea-level estimates. Last, horizontal black lines were plotted at levels on the frequency axis (y axis) that correspond to the 19-, 23-, 41-, ~95-, ~125-, 405-, 1200-, and 2400-ka Milankovitch forcing periodicities.

CO₂ records

We constructed a Bayesian hierarchical model to interpolate the atmospheric concentration of CO₂ through time and estimate its uncertainty. The primary data are a combination of a compilation of proxy data arranged in (47) and a compilation of measurements from Antarctic ice cores (91). These data contain errors associated with the measurement of both the atmospheric CO₂ concentration and the age of individual samples. The variation in measurement uncertainty across this dataset is large, spanning from an average of ±1.3 ppm (1σ) and reliable age estimates for the ice core data (98) to hundreds of ppm and Ma age uncertainties for some of the proxy data that comprise the CO₂ dataset (47). The statistical model allows us to generate the distribution of the atmospheric concentration of CO₂ that represents the best estimate of its value through time, despite the noise associated with the cloud of individual data points.

The hierarchical model has three levels: (i) a data level that represents the measured concentration of CO₂ at measured points in time, as well as the uncertainty associated with both the CO₂ concentration and temporal measurements; (ii) a process level that represents the distribution of CO₂ concentrations through time as a Gaussian process prior; and (iii) a hyperparameter level that captures the distribution of parameter values controlling the structure of the statistical model.

The aforementioned combination of CO₂ data comprises the data level of our hierarchical model. The CO₂ values, y_i (Eq. 7), and the measurement ages, t_i (Eq. 8), are modeled as the sum of their representative values (

is the mean measurement age) and their errors,

and

, respectively

(7)

(8)

At the process level, the concentration of CO₂ through time f(t) is represented as a function of five component processes (Eq. 9).

(9)

These processes are each modeled with a distinct zero mean Gaussian process prior. The component processes likely capture discrete physical processes and measurement noise operating on different time scales. For example, the processes represented by P1 and P2 capture small-scale, high-frequency variations. While these processes could persist throughout geological time, they can only be captured in our data using the ice core samples that have a mean spacing of ~400 years. However, these ice core data only extend 800 ka into the past. The Foster et al. (47) record, in contrast, spans the past 400 Ma, and correlation at longer time scales is represented by processes C1, C2, and C3. The incorporation of the higher and lower frequency processes provides the model flexibility to represent recent (<800-ka-old) variations in CO₂ concentrations that we have a dense and reliable record for, alongside the longer-term record, which is sparse and less precise/accurate.

The aforementioned structure that controls the individual processes in Eq. 9 is defined within the third layer of our model, the hyperparameter level. These hyperparameters, free parameters within the functions used to model the paleo-CO₂ concentrations, define the distribution of f(t). The priors on the time scale parameters of functions P1, P2, C2, and C3 (Eq. 9) have been assigned with the intention of capturing the effect of known orbital frequencies, and those of function C1 (Eq. 9) were set to capture longer-term variations that may be associated with tectonic processes. Function C1 (Eq. 9) was also structured to be twice differentiable to capture large-scale and gradually varying tectonic processes affecting CO₂ concentrations, such as variations in outgassing or paleogeographic configuration, if present. The hyperparameter distributions that best represented the CO₂ measurements and proxy data were ultimately determined using an MCMC routine.

Last, to estimate the effect of the distribution of potential parameter values and temporal uncertainties on f(t), we apply a random sampling/Monte Carlo method that alters the hyperparameter values and ages of the input data drawing from their probability distributions. For each iteration in the Monte Carlo application, a sample of the modeled CO₂ concentrations spanning 66 Ma to present were drawn from the estimate of f(t), incorporating the uncertainty of that estimate. The resulting distribution derived from this Monte Carlo application allows us to generate a best estimate of atmospheric CO₂ concentrations and represent its uncertainty for the past 66 Ma. Our model structure also allows us to view the expression of the individual component processes or of combinations of processes. For the purposes of this study, the processes represented by functions C1 and C2 provide the most useful representation of the variations in atmospheric CO₂ through time (Fig. 1D), because they capture the relevant long-term trends, but do not include the high-amplitude, short-term fluctuations and the measurement noise captured by functions P1, P2, and C1.