Research ArticleOCEANOGRAPHY

Sea level and deep-sea temperature reconstructions suggest quasi-stable states and critical transitions over the past 40 million years

See allHide authors and affiliations

Science Advances  25 Jun 2021:
Vol. 7, no. 26, eabf5326
DOI: 10.1126/sciadv.abf5326


Sea level and deep-sea temperature variations are key indicators of global climate changes. For continuous records over millions of years, deep-sea carbonate microfossil–based δ18O (δc) records are indispensable because they reflect changes in both deep-sea temperature and seawater δ18O (δw); the latter are related to ice volume and, thus, to sea level changes. Deep-sea temperature is usually resolved using elemental ratios in the same benthic microfossil shells used for δc, with linear scaling of residual δw to sea level changes. Uncertainties are large and the linear-scaling assumption remains untested. Here, we present a new process-based approach to assess relationships between changes in sea level, mean ice sheet δ18O, and both deep-sea δw and temperature and find distinct nonlinearity between sea level and δw changes. Application to δc records over the past 40 million years suggests that Earth’s climate system has complex dynamical behavior, with threshold-like adjustments (critical transitions) that separate quasi-stable deep-sea temperature and ice-volume states.


Ice sheets are a pivotal component of the climate system during icehouse periods. The current icehouse period started at the Eocene-Oligocene Transition (EOT), 34 million years (Ma) ago, when the Antarctic Ice Sheet (AIS) developed, followed by fluctuations in its volume with northern hemisphere glacial variations (14). Modern warming drives climate in the opposite direction, and improved estimates are needed for the associated long-term equilibrium sea level rise. This requires continuous sea level records that extend at least through the Middle Pliocene (about 3 to 3.3 Ma ago), the last time atmospheric CO2 levels were similar to those of today (47). Records should preferably continue into more ancient times to ensure coherent results across the entire range of full glacial to completely ice-free climate states, which is essential for evaluating fundamental aspects such as climate-state dependencies of feedbacks, and climate tipping points.

Ice-volume changes have dominated sea level amplitudes since the EOT. Sea level per se, however, is only reasonably constrained in continuous records for the past 500 to 800 thousand years (ka), on the basis of records of marginal-sea isolation from the open ocean (811) and on statistical assessment across various methods (12). For more ancient times, deep-sea oxygen isotope records provide the dominant information source. We want to know seawater δ18O (δw), but the data come in the form of carbonate microfossil δ18O (δc) (13, 14), which represents combined influences of δw and Tw. Analytically, Tw influences may be corrected using Mg/Ca data (15, 16), but there are complications from oceanic Mg-concentration changes over time scales longer than a million years (1719). Moreover, Mg/Ca thermometry has 1° to 1.5°C uncertainties (1σ) in cold deep-sea environments (15), while Tw signals commonly amount to only a few degrees, causing unfavorable signal-to-noise ratios (16). Even if Tw could be resolved without uncertainties, the resolved δw record would still need to be translated into ice-volume/sea level changes. This is commonly done by assuming a linear ΔδwSL relationship, which is often “calibrated” to the Last Glacial Maximum (LGM)–to-modern δw gradient measured in marine sediment pore waters (2022). However, these assumptions remain untested.

The ubiquitous reliance on linear ΔδwSL approximations is unexpected given that, from first principles, the primary expectation should be that it is nonlinear. This is because ice sheets accumulate precipitation of increasingly negative δ18O composition as they grow because of increasingly intense Rayleigh distillation and because different ice sheets begin to form and develop their δ18O signature at different moments in time (and, hence, at different moments in global ΔSL history). Thus, the global mean ice δ18O will change with sea level, which imposes nonlinearity on the ΔδwSL relationship. Here, we evaluate this issue quantitatively using a new approach based on straightforward representations of underlying process relationships and interdependencies, and we apply our mutually consistent reconstruction method for key parameters in paleoceanographic reconstructions. We start with a simplified analytical assessment of the ΔδwSL relationship to illustrate its fundamental nature and sensitivity to key assumptions and uncertainties. Thereafter, we present applied assessments that use published sea level records to determine mutually consistent variations in ΔSL, ice sheet volumes (Vice), δice, Tw, and δw, culminating in a detailed assessment over the past 40 Ma. Sensitivity tests and a range of independent validation criteria are used to assess the robustness of results, and key avenues for further improvements are highlighted.

Our mutually consistent, process-based framework for ΔSL, Vice, δice, Tw, and δw permits complete system validation using multiple parameters rather than only one or two as in traditional approaches. We obtain new records for these parameters, which are of fundamental interest in (paleo)climate science, and that are independently cross-validated over at least the past 22 Ma. This approach provides a testable way to understand relationships among measurable parameters and underlying processes of interest in understanding climate change.


Relationship between precipitation δ18O and ice sheet volume

The volume (V) of ice sheets is a key determinant in the δ18O of their newly accumulating precipitation (δp). Past work has related low δp, reflected in the δ18O of accumulated ice (δice), to intense Rayleigh distillation of atmospheric vapor through successive vapor saturation/condensation cycles with decreasing temperature along vapor pathways from evaporative sources to condensation/precipitation sites [e.g., (23, 24), review in (25)]. Rayleigh distillation intensifies with cooling during atmospheric vapor transport from oceanic evaporation sites to high latitudes, from coastal regions to cold (winter) continental interiors, and from lower to higher altitudes. This is because vapor saturation drops nonlinearly with decreasing temperature (Clausius-Clapeyron relationship), while Rayleigh distillation follows a natural logarithmic relationship to more negative values with reduction in the fraction of remaining atmospheric vapor (25). As a result, relatively low-altitude seasonal snow over ground or small snow patches has much less negative δp than snow falling at the center of a large, high-altitude continental ice sheet. While this is an idealized representation and there are additional controlling factors (26), we address this both here with a sensitivity test (below) and in our final outlook for further method development.

Seasonal snow δp typically reaches −10 to −20‰ [e.g., (2729)]. Values become more negative with distance from the coast and with increasing altitude; for example, a trans-Alaskan continental gradient exists of −8.3‰/1000 km from the coast, with a superimposed orographic δ18O decrease of −6.8‰/km of altitude (28). As regional ice-albedo effects combine with high altitude, intense cooling in cold mountain regions with permanent ice (glaciers and ice caps) cause δp to typically reach −25‰ or even −30‰, with lowest values recorded in winter [e.g., (3033)]. Hence, a sharp initial δp drop is apparent between initial snow and sustained precipitation over nascent high-altitude, permanent ice bodies. Eventually, δp can reach −60‰ or lower in the coldest, most isolated locations of the intensely cold hearts of major, kilometers-high continental ice sheets, such as the Antarctic interior (34). Rayleigh distillation between evaporative sources and precipitation sites was less extreme for the North American Laurentide Ice Sheet (LIS) and the Eurasian Ice Sheet (EIS) during ice ages because of their positions in lower (“warmer”) latitudes closer to relatively warm oceanic moisture sources. For example, five isotope-enabled global climate circulation models indicate that the LIS reached a minimum LGM δp of about −38‰ (35).

We represent δp changes over ice sheets as a generic function of V, starting with an initial snow δp = −15‰ that drops rapidly with initial ice sheet volume buildup, followed by an exponential change with V toward very negative values over large ice sheets (Methods, Eq. 1, and fig. S1). We use two sets of equation constants to account for more intense Rayleigh distillation over high-latitude “cold” ice sheets than over lower-latitude warmer ice sheets. In our analytical assessment, sensitivity test 1 uses a different relationship that is simply linear (Methods and Eq. 2)—albeit physically implausible— to evaluate the importance of the selected function shape for the conclusions, given that our idealized representation may overlook additional controls on δp (26).

Analytical assessment

We use simple growth descriptions for the AIS as a single entity that represents both the East and West AIS, Greenland Ice Sheet (GrIS), EIS, and LIS, where the latter refers collectively to the LIS and Cascadian Ice Sheet. VAIS today is 57.8 m sea level equivalent (mseq) and VGrIS is 7.3 mseq (36) so that ΔSL = +65.1 m in an ice-free state. In this analytical assessment, we assume that AIS and GrIS grow proportionally between ΔSL = +65.1 and 0 m. Growth of continental ice sheets (i.e., glaciation) relative to the present caused negative ΔSL, mainly due to growth of LIS and EIS (which do not exist today) along with slight AIS and GrIS expansions. Common approximations for the LGM are that ΔSL was roughly −125 m because LIS had reached of the order of 70 mseq and EIS about half that, while AIS grew by roughly 15 mseq [see overview in (37)]. Here, we simply assume that GrIS grew another ~5 mseq, which may be an overestimate, but this does not appreciably affect results, especially because less GrIS growth would need to be compensated by more growth in other ice sheets to keep the sea level budget closed. For different glacial maxima, size dominance may have alternated between LIS and EIS [see overview in (37)], but this does not affect results either because both are similar types of “warm” ice sheet in our calculations. We assume initially that all expansions were linearly proportional to ΔSL. This gives individual ice sheet V variations relative to ΔSL (Fig. 1A) and their sum Vtot (Fig. 1C).

Fig. 1 Results of the analytical assessment.

(A) Volume of each ice sheet relative to total sea level change caused by all ice sheet changes. (B) Ice δ18O for each ice sheet. (C) Global volume- (or mass-) weighted mean ice δ18O for all accumulated ice (note that in terms of proportions, mass weighting and volume weighting are equivalent here). (D) Changes in seawater δ18O relative to sea level change. Solid lines indicate the main case (Methods and Eq. 1), dashed lines are for sensitivity test 1 (Methods and Eq. 2), and dotted lines are for sensitivity test 2. Shaded intervals in (D) indicate linear ΔδwSL relationships that would occur if global mean ice δ18O were constant at the indicated values. The overall steepness of lines in (D) is addressed in our applied assessments.

The schematic ice sheet V variations can now be combined with the chosen relationship between δp and ice sheet volume (Methods and Eq. 1 or 2) to determine δp developments for each ice sheet (Fig. 1B). We assume that all ice sheets are equilibrated with δp; that is, the mass-weighted mean δice of each ice sheet (indicated by δice¯) is assumed to equal δp (this simple instantaneous equilibration assumption is addressed in detail below in our applied assessments). Next, we calculate the mass-weighted mean global δice value for all accumulated ice (i.e., δice*; Methods and Eq. 3). Last, relative δw changes (i.e., Δδw) are determined using an ice-to-water density ratio of 0.9, a value for global ocean volume (Voc) of 3700 mseq, and ΔSL = − {(VAIS – 57.8) + (VGrIS – 7.3) + VLIS + VEIS} (Methods and Eq. 4). In sensitivity test 2, the AIS is set to grow nonlinearly during glaciation (ΔSL < 0 m), giving it a somewhat increased expansion rate during final glaciation stages. This test investigates the possible influences of a hypothetically enhanced AIS contribution to a marked final sea level drop (38) that coincided roughly with AIS expansion to its maximum LGM extent (39); the overall sea level budget is kept unchanged by setting LIS to grow nonlinearly in the opposite sense (Fig. 1A). We find (below) that this nonlinearity is too small to significantly affect the outcome of our analysis and the same holds if we invert the nonlinearity between AIS and LIS [the condition suggested in (38)].

The main case and both sensitivity tests reveal pronounced nonlinearity between ΔSL and δice* and, thus, between ΔSL and Δδw (Fig. 1, C and D). Similar ΔδwSL relationships for all three cases indicate that the pattern is robust with respect to the choice of V:δp relationship, details of ice sheet growth histories, or both (Fig. 1D). This is because the nonlinearity arises principally from addition of two much less negative ice sheets at ΔSL < 0 (which increases δice*) to the virtually fully formed and isotopically very negative AIS. Thus, the nonlinearity is a robust, unavoidable consequence of different glaciation thresholds for different ice sheets at different latitudes. Even if the same δp function were used for LIS and EIS as for AIS and GrIS, considerable nonlinearity would remain because LIS and EIS in their initial stages (up to 50 mseq) have more positive δp than AIS.

Applied assessments

Here, we use published sea level records to determine mutually consistent variations in ice sheet volumes, δice, Tw, and δw. A flow diagram for the entire method is presented in fig. S2. From this workflow, the initial step of regression-based conversion of δc changes into ΔSL is omitted when results are produced for published sea level records for comparison with δc-based results. Our method uses a straightforward model for growth of generic ice sheets with circular planoconvex lens shapes, axially symmetric parabolic profiles (40), and a constant height:radius aspect ratio ε = 2 × 10−3 (roughly the mean aspect ratio of the modern AIS). Unless specified otherwise, global ice-volume variations are determined in terms of sea level–equivalent changes, accounting for the density difference between ice and water. Model ice sheets retain their shape characteristics throughout, while sea level change between individual time steps is used to determine global net ice-volume growth or loss, which is partitioned over the different ice sheets as described in Methods (Eqs. 5 to 11). Thermosteric sea level changes are not considered.

The model contains one large and one smaller ice sheet that both disappear at the present-day sea level; these approximate LIS and EIS. It is irrelevant which of the two grows larger and which stays smaller, as all model ice sheets are generic, without geographic attribution, and both LIS and EIS approximations are taken to be warmer ice sheets for Rayleigh distillation (Methods and Eq. 1). The model approximates GrIS with VGrIS = 7.3 mseq at ΔSL = 0 m and AIS with VAIS = 57.8 mseq at ΔSL= 0 m (36). Ice-volume changes are driven in accordance with the applied sea level records and are attributed to the various ice sheets as discussed in Methods (Eqs. 5 to 11). δp changes over the ice sheets are calculated as in the main case of the analytical assessment above, using Methods and Eq. 1.

At this stage, we abandon the simplistic assumption that δice¯ for each ice sheet is instantaneously similar to δp and instead calculate its time evolution (Methods and Eqs. 12 to 14). The initial assumption needed to be abandoned because δp instantaneously tracks ice sheet dimension changes, whereas it takes time for δice¯ to equilibrate with changing conditions. Internal ice sheet properties reflect surface accumulation, downward and outward ice flow (ice dynamics), and ice loss through melting and calving. Thus, ice sheets contain contiguous ancient ice sequences within their interiors; in GrIS, this dates back to at least 125 ka ago (41) and in AIS to at least 800 ka ago (42, 43). AIS has been at (almost) its modern size long enough that even ancient ice formed from precipitation with δp close to present values, but in short-lived ice sheets, especially LIS and EIS, a noticeable lag is expected between ice-volume buildup and full δice¯ equilibration to the volume change. This lag is roughly 17 to 20 ka/mseq, based on modern maximum ice ages and the GrIS and AIS sizes, but may be ice sheet size and, thus, time dependent. Given this uncertainty, we refrain from applying a standard lag function and instead resolve δice¯ evolution directly per time step (Methods and Eq. 14). Here, δice¯j is calculated at time step j using ice sheet volume and δice¯j1 of preceding time step j − 1, and both gross mass accumulation with more negative δice (determined from the δp relationship above) and gross mass loss with δice that is (before achieving isotopic equilibration) less negative between time steps j − 1 and j. We use δice¯j1 for the latter term, along with a gross mass loss term that is determined by gross accumulation minus net accumulation/loss (where loss is negative). Net accumulation/loss is obtained from sea level variations, and gross accumulation (agrs) is calculated in an ice sheet surface-area–dependent manner directly proportional to modern global annual agrs of ~0.008 mseq [following rates in (44, 45)] for a global ice volume of 65.1 mseq.

Annual agrs is ~0.01% relative to ice volume and is set to change proportionally with ice sheet size, so the model ice sheets have isotopic equilibration time scales of the order of 104 years. This time scale is especially relevant for LIS and EIS, which existed only during glacial stages with life spans of order 104 to at most 105 years, and that typically reached maximum sizes only millennia before rapid glacial terminations. As a result, LIS and EIS were continually playing catch-up in terms of isotopic equilibration to their marked size changes. In the Late Pleistocene, the more permanent AIS (and, to a lesser extent, also GrIS) largely maintained conditions close to isotopic steady state. This did not hold during the Pliocene and older times, when considerable sea level changes well above 0 m imply major VGrIS and VAIS changes.

On the basis of ice-volume change, which is equated to sea level change (ΔSL), and the combination of all ice sheet V and δice¯ values, we can now determine Δδw (Methods and Eq. 15) and, thus, the ΔSL:Δδw relationship. Each ice sheet’s contribution to Δδw depends on its V and δice¯, so that the aforementioned time lag between volumetric changes and δice¯ equilibration holds considerable potential to influence ΔSL:Δδw nonlinearities.

Subtraction of the calculated Δδw record from δc [e.g., (13, 14)] resolves the δc component that relates to Tw changes by −0.25‰ °C−1 (15, 16). These Tw changes can then be compared with Mg/Ca-based Tw reconstructions to assess model-based Vice, δice, Tw, and δw variations. However, chronological offsets between model-input ΔSL and calculated δw on the one hand, and δc on the other hand, may cause large spurious Tw signals that obfuscate validation. Hence, closely synchronized records, or—better still—signals coregistered in a single sample set, are essential. We, therefore, only use Tw validation when this condition is met.

We apply our method to several input ΔSL records, each with different underpinning assumptions and uncertainties (Figs. 2 and 3A). We use a 550-ka Red Sea marginal-basin residence-time record multiplied by 1.1 to approximate global mean sea level changes (811), an 800-ka multimethod sea level stack (12), and a 5.3-Ma Mediterranean marginal-basin record that is known to give pre–1.2 Ma values that are biased too high (Fig. 2B) (46). These ΔSL reconstructions are not entirely independent of each other; notably, the Red Sea (811) and Mediterranean (46) records each contribute ~1/7 of the multimethod sea level stack (table S1) (12). This does not affect our approach because results from these ΔSL reconstructions are shown only for mutual comparison, to illustrate that our method is broadly applicable and not specific to any one realization.

Fig. 2 Sea level records used in this study.

(A and B) Comparison between the sea level solutions based on the Red Sea record [red; (811)], the statistical sea level record [blue; (12)], the Mediterranean Sea record that is known to be biased high before 1.2 Ma ago [magenta; (46)], the benthic δc stack-based sea level record [black; (13)], and western Mediterranean cave-based benchmarks [purple; (49, 50), where the youngest three points come from the more recent publication] [(A) is an expanded image of the last 800,000 years from (B)]. The same age scale is used in (B) and (C). (C) Comparison between the benthic δc stack-based sea level record [black; (13)] and the benthic δc megasplice-based sea level record [orange; (14)] (both based on the main case regression in fig. S3) and western Mediterranean cave-based benchmarks [purple; as in (B)]. The δc megasplice-based record has greater amplitudes, which may suggest that some amplitude reduction/smoothing has resulted from δc stack calculation. All chronologies are as presented in the original studies. In all three panels, the light blue rectangles indicate multiproxy best-fit sea level estimates for MIS 5e and 11; these rectangles are to scale only for sea level, not for age (48). Our two benthic δc-based sea level records reach the lower regions of these multiproxy best-fit interglacial sea level estimates, which likely results from failure of the benthic δc-based sea level records to fully resolve millennial-scale features, which are important for capturing range extremes (51).

Fig. 3 Comparison of our ΔSL, Tw, and δw with other observations.

(A) Sea level change from the sea level stack [(12), blue]; from the Red Sea record, × 1.1 to approximate global mean sea level variations [(11), orange]; from model-based δc deconvolution [(4), green]; and from lag-optimized, quadratic regression between (13) and (12) (fig. S3) forced to peak at ΔSL = 65.1 m (black), and from the regression’s upper 95% confidence limit (pink). Cyan boxes indicate maximum amplitude envelopes of New Zealand relative sea level fluctuations (52), vertically positioned for comparison. Green dots are western Mediterranean benchmarks (49, 50), with highlighted Middle Pliocene range (yellow box). (B) Tw changes derived using δc and δw, where δw is calculated from ΔSL [this study; black and pink using ΔSL in (A); orange using ΔSL from the unconstrained regression; fig. S3]. (C) to (E) all on the age scale of (E). (C) Last 1.6 Ma of (B) versus Antarctic (EDC) temperature on its independent chronology (54). Red dot represents a noble gas–based estimate of LGM global ocean cooling, plotted against the Tw axis (57). (D) As (C) but versus Mg/Ca-based Tw [(16), blue, individual data and 4-ka Gaussian smoothing]. (E) As (D), but for δw. We show ODP Site 1123 (16) alone at 0 to 0.35 and 1.45 to 1.6 Ma ago and a three-record δw stack including ODP 1123 for 0.35 to 1.45 Ma ago, with 1× and 2× bootstrap errors (58). Red dot indicates LGM sediment pore-water–based δw (2022).

To allow application of our method to more ancient times, we use two ΔSL records from deep-sea benthic foraminiferal δc records. These δc records are a 5.3-Ma stack (13) and the youngest 40 Ma of a 67-Ma global megasplice (Figs. 2, B and C, 3A, and 4A) (14). ΔSL fluctuations are inferred from these δc records using a lag-optimized quadratic regression between the 800-ka multimethod ΔSL stack and δc (12). This regression was originally constrained up to about ΔSL = +15 m and, here, needs extrapolation to ice-free conditions. Therefore, we reevaluate the quadratic regression while constraining it to peak at ΔSL = +65.1 m.

Fig. 4 Validations for our 40-Ma records.

(A) Sea level variations from δc (14) based on our main-case regression (black), upper 95% confidence bound (pink), and unconstrained quadratic regression (orange) (fig. S3); and previous model-based δc deconvolution [(4), green]. (B) Our results for AIS (black, pink, and orange as above) and northern ice volume (LIS + EIS + GrIS; red, only for the main regression for clarity). (C) Model-derived δw based on our ΔSL in (A) (black, pink, and orange) versus δc [(14), lilac]. (D) Tw changes based on δc residuals in(C) (black, pink, and orange). Regular italic numbers indicate positive validations: 1, western Mediterranean benchmarks (49, 50); 2, best- and maximum-amplitude sea level ranges (53); 3a and 3b, first major iceberg calving in Nordic Seas (59) and LIS calving onset (60), respectively; 4, AIS-volume variability range (61); 5, onset partial/ephemeral northern ice (1); 6, end of last intermittently ice-free period (3); 7, AIS glaciation onset (1, 3, 4, 64, 6670); 8, MMCT δw change of 0.35 ± 0.12‰ (62); 9a and 9b, loess fit through Mg/Ca Tw with 95% confidence bounds; at 34.5 to 37.5 Ma ago, data were questioned (67). Roman numerals indicate discrepancies: i, amplitudes of sea level change in North West Australia (63); ii, EOT sea level drop of 70 to 80 m (64) versus ~30 m in this study (Discussion); iii, 2.5°C EOT cooling (66) versus 3.5°C cooling here (Discussion); iv, MMCT cooling of 1.5° ± 0.5°C (62) versus 2.5° ± 0.5°C here.

Two other regressions are used as sensitivity tests. One represents the upper 95% confidence limit of the constrained quadratic regression. The other represents an unconstrained quadratic regression through the same dataset (see Methods and fig. S3). Results for the latter are shown in orange in Figs. 3A and 4A but are virtually indistinguishable from the main case (black lines) over the Plio-Pleistocene (Fig. 3) and top out at ~56 m in more ancient times so that this record suggests consistent ice presence before the EOT at 34 Ma ago (Fig. 4A). In view of reports that some ice may have existed in Antarctica before the EOT [e.g., (1, 47)], we cannot exclude the orange estimates. However, given that these estimates retain significant ice volume before the EOT, we consider them to be unlikely minimum sea level estimates. The main constrained regression (black) provides the lowest sea level estimates that are still consistent with an ice-free pre-EOT state and limited Middle Miocene ice conditions—these represent our best estimates. Conversion regressions from δc to ΔSL that are steeper than the upper 95% confidence limit of our forced regression (pink in Figs. 3A and 4A) lead to unrealistic, long-lasting, ice-free conditions through the Middle Miocene; we, thus, consider the pink line as an upper sea level estimate. Our approach differs from previous modeling that deconvolved δc into ΔSL, and Tw components, where the latter are assumed to be northern hemisphere–driven (4). We do not use that ΔSL record as input because its origin in δc would risk introducing circularity; we use it for comparison only.

Modeled time series for Vice and δice¯ based on the various ΔSL records are presented in figs. S4 to S6. Although differences between the input ΔSL records reverberate through the solutions in an absolute sense, the choice of ΔSL record, as long as it has a realistic amplitude and timing structure (relative to previous sea level records and ice volumes of the ice sheets concerned), does not affect relative relationships among Vice, δice¯, Δδw, and Tw because those relationships arise from the modeled processes rather than from input values. Therefore, plotting of parameters from individual model runs provides substantial information on time-dependent behavior among parameters, across the full range of sea level values covered by the input ΔSL records (Figs. 2 to 6). All variables of interest are mutually consistent so that validation of one typically implies validation of the others (Figs. 3, A and C to E, and 4, A to D). Before evaluating validation against observations, we make two initial observations.

Fig. 5 Diagrams of δice¯ relative to ice sheet volume for individual ice sheets for the five main model runs investigated in this study.

(A) Sea level record based on the Red Sea method (811), (B) sea level stack (12), (C) benthic δc stack-based sea level record (main regression in fig. S3) (13), (D) sea level record from the Mediterranean Sea method (46), and (E) youngest 40 Ma of the benthic δc megasplice-based record (main regression in fig. S3) (14). The less complete Vice:δice¯ representations are for shorter records that do not cover the full bipolar ice sheet response range; the more detailed (D) and (E) are more representative. The greater noise in (D) results from lower precision of the Pliocene Mediterranean ΔSL record. (F) Schematic illustration of the nature of the relationships in (A) to (E). In (F), 1 represents the trajectory associated with gradual ice-volume buildup, determined by continuous instantaneous ice volume–based adjustment of δp, and lagged adjustment of δice¯ according to the residence-time calculation discussed in Methods (Eqs. 12 to 14). 2 represents rapid ice loss during deglaciation that occurs at the δice¯ value attained just before deglaciation onset; and 3 represents adjustment at the end of deglaciation, when new ice starts to build up at the starting δp value. 4 is the trajectory associated with gradual partial glaciation (as 1), 5 is rapid partial deglaciation (as 2), and 6 represents more gradual δice¯ adjustment to conditions commensurate with the remaining ice volume after partial deglaciation.

Fig. 6 Regressions between ΔSL and Δδw.

(A) Regressions based on runs for three individual sea level approximations: light and dark blue based on the Red Sea record (811); orange and red based on a stack of different records (12); and gray and black based on a lag-optimized quadratic regression between benthic δc (13) and the sea level stack (12), forced to peak at ΔSL = 65.1 m (fig. S3). (B) All data from (A) combined, with regression in yellow. (C) As (B) with superimposed (magenta) results from a run based on the Mediterranean Sea record (46) and (light blue) results based on benthic δc (13) and the sea level stack of different records (12), constrained to peak at ΔSL = 65.1 m. (D) Black = green + magenta data from (C), and blue is all three datasets from (C) combined with regression in orange. The regression is a fifth-order polynomial for optimal visual fit through the inflection in the data cloud with least edge effects: Δδw = 9.6 × 10−11ΔSL5 + 1.9 × 10−8ΔSL4 + 2.5 × 10−7ΔSL3 − 1× 10−4ΔSL2 − 0.015ΔSL − 0.133.

First, slow isotopic equilibration to volumetric changes causes marked V versus δice¯ hysteresis for each ice sheet (Fig. 5). This arises from slow (lagging) δice¯ adjustments to growing ice sheets until the glacial maxima when the most negative δice¯ was achieved and subsequent rapid and large-amplitude deglacial net ice losses (Fig. 5). At each (partial) deglaciation, values shoot almost horizontally to the left in the hysteresis diagrams, leading into a jump that can be almost vertical—especially in the most marked cases of complete LIS and EIS deglaciation—as incipient glaciation starts again with new snow at δp = −15‰. The AIS only undergoes these fluctuations in the hysteresis diagrams when it experiences substantial volume reductions (i.e., at ΔSL > 0 m, during the Pliocene and more ancient times).

Second, the ΔδwSL relationships obtained from paleorecords (Fig. 6) clearly follow the same general ΔδwSL nonlinearity found in our analytical assessment (Fig. 1D). This substantially affects the way Δδw is reconstructed through time relative to traditional approaches that assume a linear ΔδwSL relationship (fig. S7). By resolving the mass-weighted mean global δice value for all accumulated ice (δice*) at every time step in accordance with evolving characteristics of the ice sheets, our reconstruction allows Δδw to shift seamlessly across lines obtained from ΔδwSL relationships based on different δice* values. This removes the need to assume which δice* value applies at specific times and how and when shifts between these values occurred. Our method thus provides a major step toward the necessity—specified in the most recent long-term sea level assessment—that “changes in the [… ΔδwSL …] calibration with evolution of ice-sheet size should be modeled” (3). We note that the shape of the ΔδwSL nonlinearity identified from our modeling remains robust even when combining data points from model runs for all ΔSL records, with all of their different underpinning assumptions and uncertainties (Fig. 6D). This indicates that the shape is largely insensitive to the distinct V:δice¯ hysteresis for individual ice sheets, which was omitted from our analytical assessment. There are impacts only on the overall slope of the lines and data distributions around the regression. However, such robustness does not imply that the model approximates reality well; it merely reflects internal model consistency. Validation against independent observations is needed to assess whether the calculated records are realistic (Discussion and Figs. 3, A, C, and D, and 4, A to D).


Validation of results

Clear linear progression in our method illustrates an absence of circularities; fundamental assumptions exist only at two stages (fig. S2). One assumption concerns regression-based derivation of ΔSL from δc (fig. S3). Yet, sensitivity tests with realistically different regressions (fig. S3) indicate that our δc-based ΔSL record is robust within about 10 m (total range; Figs. 3A and 4A). Diverse validations of the mutually consistent ΔSL, δice¯, δw, and Tw records support this robustness (below). The other assumption concerns the chosen δp versus ice-volume relationships. Sensitivity tests with alternative linear relationships indicate that our conclusions are robust regardless of the relationship shape (Fig. 1, C and D).

Below, we compare our reconstructions with independent observations. A first validation against observations is that even at the largest LIS volumes, LIS δice¯ barely reached −35‰ (Fig. 5), while values integrated over 100-ka glacial cycles within the past 500 ka are consistently between about −25 and −35‰ (fig. S5). Groundwaters in North America that reflect subglacial or preglacial recharge typically reach −26‰, with a possible −28‰ end member, while ground ice in Yukon and the Barnes Ice Cap, Nunavut (both in Canada) has values of −29 to −35‰ (35). These values agree well with our LIS δice¯ simulations.

Next, we observe reasonable agreement between our new δc regression–based ΔSL and independent multiproxy sea level reconstructions for interglacial Marine Isotope Stages (MIS) 5e and 11 (48), the Red Sea sea level record (11), and Middle Pliocene sea level benchmarks (Figs. 2, A and C, and 3A) (49, 50). The δc records that underpin our regression-based ΔSL do not fully resolve millennial-scale features, so our new (δc regression–based) ΔSL records are likely to miss short-term extremes of several-meters amplitude that can be important in interglacials like MIS 5e and 11 [e.g., (51)]. Compared with the Red Sea record (11), our new ΔSL records agree much better with the multiproxy assessment for MIS 11 (48), a further indication that the Red Sea record may underestimate sea level at that time. All Middle Pliocene sea level benchmarks shown (49, 50) agree with our ΔSL, although the youngest does so only at its extreme age uncertainty. Our ΔSL may offer a means to refine age estimates for benchmarks with relatively large dating uncertainties at times of large ΔSL variations.

Agreement with the 800-ka multiproxy sea level stack (Figs. 2B and 3A) (12) cannot be used for validation because it underpins the regression that provides our ΔSL (Methods and fig. S3). In addition, we note that a previous δc-based ΔSL deconvolution (4) also agrees well with our ΔSL between 0 and 3 Ma ago, albeit at somewhat muted amplitudes with glacials only reaching −100 m (Fig. 3A). Between 3 and 5.3 Ma ago (Fig. 3A), and extending to 13 Ma ago (Fig. 4A), however, that record (4) has much less ΔSL variability than our results, which requires comparison with additional independent evidence.

In the 3.3– to 3.0– and 2.9– to 2.65–Ma ago intervals, amplitude variations in both our ΔSL record and the previous δc-based record (4) agree with maximum relative sea level amplitude ranges of up to 27 and 35 m in New Zealand, respectively (Fig. 3A) (52). However, Middle Pliocene sea level benchmarks (Figs. 3A and 4A) (49, 50) and Miocene best and maximum sea level amplitude estimates from New Jersey and the Delaware Coastal Plain (Fig. 4A) (53) indicate that the previous δc-based ΔSL deconvolution (4) is anomalously flat and low between 3 and 13 Ma ago and that its amplitude variations between 13 and 22 Ma ago exceed independently established amplitude constraints (53). In contrast, our ΔSL record agrees well with both these independent lines of evidence (Figs. 3A and 4A). The difference may arise from the fact that our method operates in the same way across all time frames, based on generic warm and cold ice sheets. In contrast, that of de Boer et al. (4) was primarily set up to account for LIS and EIS variations, driven by mid-latitude to subpolar northern hemisphere temperatures inferred from δc, and used a tuning factor for AIS so that a large VAIS response was obtained across the EOT. Possibly, the tuning factor used by de Boer et al. (4) is too strong, resulting in overly sensitive VAIS responses through the Oligocene-Miocene, achieving a “full” AIS state that allowed no further expansions, which results in a flat ΔSL curve until northern hemisphere glaciation began. Alternatively, the assumption that mid-latitude to subpolar northern hemisphere temperatures could be inferred from δc (4) may be incorrect.

We favor the latter explanation, on the basis of comparison of our calculated Tw variations with independently determined Antarctic temperature variations at European Project for Ice Coring in Antarctica (EPICA) Dome C (EDC) (Fig. 3C) (54). This reveals excellent amplitude structure agreement between our Tw and TEDC, scaled as ΔTw ≈ 0.25 ΔTEDC. Offsets in this comparison only concern the chronologies between the two records, between 0.4 and 0.7 Ma ago. We infer that global mean deep-ocean temperature has consistently been set by temperature fluctuations in high southern latitudes rather than in the mid-latitude to subpolar northern hemisphere as was assumed previously (4). This is an important observation for understanding ocean-atmosphere feedbacks on these time scales. In addition, this observation may help in developing chronologies for old (partially deformed) ice-core sections and blue-ice sections, which are difficult to date because they lack stratigraphic continuity [e.g., (55, 56)], by adding independent age controls from Tw matching.

We also compare our Tw record with a 1.6-Ma Mg/Ca reconstruction for Ocean Drilling Program (ODP) Site 1123 from Chatham Rise, east of New Zealand (3290 m in water depth) (16). Reasonable agreement exists over the interval 0 to 1.2 Ma ago (Fig. 3C). Offsets occur at 1.3 to 1.45 Ma ago, where the Mg/Ca record suggests Tw drops to −4°C relative to present. Given a modern in situ bottom-water temperature of just above 1°C at Site 1123 (15) and adjustment to potential temperature of less than −0.2°C, a −4°C Tw drop at Site 1123 would suggest deep-water formation at surface temperatures just under −3°C. This is considerably below the seawater freezing temperature of about −1.8°C, which suggests that there is an issue with the Mg/Ca record in this interval. Additional temperature validation is provided by ice-core noble gas analyses, which indicate an LGM global mean ocean cooling of −2.57° ± 0.24°C (57), which agrees with our LGM global deep-sea cooling of −2.5° ± 0.3°C (Fig. 3, C and D). However, Mg/Ca temperatures at Site 1123 suggest LGM cooling of only −1° to −1.5°C (Fig. 3D).

The ODP Site 1123 Mg/Ca-based Tw has been used with coregistered δc to calculate δw (16). These data were later used as part of a three-record stack for the interval 0.35 to 1.45 Ma ago (58). These δw records compare reasonably with our δw results, especially in glacial intervals (Fig. 3E). This is a valid comparison to add to the Tw comparison because observational studies primarily resolve Tw and secondarily infer δw, while the opposite applies to the present study (fig. S2). Moreover, the studies use independent δc records. Additional δw validation comes from marine-sediment pore-water δw measurements. These reveal an LGM-modern gradient of about 1 ± 0.1‰ (2022), while we find 1.1 to 1.4‰ across the different runs; namely, 1.4‰ in the Red Sea–based example (811), 1.4‰ in the example for the statistical multiproxy sea level stack (12), 1.2‰ in the Mediterranean-based example (46), 1.1‰ based on our δc-to-ΔSL conversion using the benthic stack (13), and 1.2‰ based on our δc-to-ΔSL conversion using the benthic megasplice (14). The difference likely falls within assumptions and uncertainties in both methods, especially given that reevaluation of the pore-water method has indicated larger uncertainties (22) than considered originally (20, 21).

In Fig. 4B, additional validation criteria labeled 3a, 3b, and 5 indicate broad validation of key northern hemisphere (LIS + EIS + GrIS) glaciation stages (red) with ice-rafted debris (IRD) records from the Nordic Seas and wider North Atlantic (59, 60) and previous inferences of an onset of ephemeral or partial northern hemisphere ice development, respectively (1). Criterion 4 corroborates VAIS fluctuations (61), while 6 marks agreement about the end of the most recent intermittently ice-free period (3); we note that 4 and 6 are not fully independent because they partly rely on δc information (3, 61). AIS glaciation onset at the EOT (criterion 7) marks a well-studied greenhouse-icehouse transition that provides binary validation (ice versus no ice) although some suggest somewhat earlier major glaciation from ~36 Ma ago (47). In addition, a δw change of 0.35 ± 0.12‰ has been inferred across the Middle Miocene Climate Transition (MMCT) (62), which agrees well with our results (criterion 8). Last, a Mg/Ca-based Tw compilation over a broad Eocene-Oligocene interval indicates excellent agreement with our result between 37.5 and 40 Ma ago (criterion 9a) and reasonable agreement for the lowest Tw values around the EOT (criterion 9b; allowing for potential chronological offsets). In the Early and Middle Oligocene, values in the Mg/Ca compilation generally fall ~1°C above our estimates, albeit just within uncertainties. Between 31 and 28 Ma ago and after 26 Ma ago, this offset is reduced noticeably so that values agree well within uncertainties.

Use of multiparameter approaches to validate our reconstructions goes a long way to excluding potential diagenetic impacts on benthic δc. Following our method, diagenetic alteration of δc to more positive (negative) values would cause spuriously low (high) ΔSL anomalies that would imply larger (smaller) ice volume. This, in turn, would drive positive (negative) Δδw anomalies, whose subtraction from more positive (negative) δc would have a canceling effect, resulting in reconstruction of minor to negligible Tw changes. The only interval with limited (though still substantial) Tw variability is between 13 and 5 Ma ago (Fig. 4D), but ΔSL, AIS volume, and onset of northern hemisphere ice presence from our method remain well validated (criteria 1, 2, 4, and 5; Fig. 4, A and B), which indicates that low Tw variability in this interval does not arise from diagenetic alteration of δc.

There are also some discrepancies. Deep-sea cooling at the MMCT was estimated at 1.5° ± 0.5°C (62), which is ~1°C less than inferred here (discrepancy iv, Fig. 4D). In addition, at around that time, large-amplitude ΔSL fluctuations have been inferred from stratigraphic analyses of northeastern Australian sequences (Fig. 4A, discrepancy i) (63). However, these disagree not only with our results but also with the ΔSL record of de Boer (4) and with criterion 2 (53). Discrepancies ii and iii are now evaluated in more detail.

Previous cumulative EOT ΔSL estimates reach −70 to −80 m (3, 64), which greatly exceed our estimate of about −30 m (discrepancy ii). The ΔSL record of de Boer (4) indicates about −45 m of sea level change across the EOT (accounting for different pre-EOT values in Fig. 4A), similar to −30- to −50-m estimates of relative sea level fall from New Jersey Coastal Plain backstripping analysis (65); both are closer to our estimate than to the −70- to −80-m estimates. This EOT ΔSL discrepancy must be considered together with the difference between Mg/Ca-based ~2.5°C EOT cooling (66) and our inferred ~3.5°C cooling (Fig. 4D, discrepancy iii). The ~1°C difference scales to ~0.25‰ in δ18O. Given that the AIS at the time was only beginning to grow, its δice¯ would not have been very negative. This is evident in Fig. 4 where Δδw = 0.25‰ for ΔSL = −30 m at the EOT. If the ~2.5°C EOT cooling (66) is correct—which remains to be settled in view of agreement between our results and both high Eocene and lowest EOT Tw [criteria 9a and 9b (67)]—then the implied 0.25‰ of “missing” δ18O would correspond to another −30 m of ΔSL, giving a total of −60 m that is similar to the modern VAIS. Moreover, AIS may, at the time, have behaved more like a warm ice sheet (similar to the Pleistocene LIS and EIS), with early-growth δice¯ only around −25‰ and full-size δice¯ only around −40‰ (68). This would make the ΔSL of −60 m an underestimate and bring potential values closer to previous −70- to −80-m estimates (64). If this were the case, then a fundamental shift is implied in the global ΔδwSL relationship at some stage between 34 and 22 Ma ago. This might be related to a Rayleigh distillation change over the AIS from warm ice sheet behavior to cold ice sheet behavior so that the ΔδwSL relationship in Fig. 6 would apply only to later times (<22 Ma ago). It is also possible that Mg/Ca-based Tw estimates across the EOT are problematic; for example, because of ocean Mg concentration changes and carbonate saturation changes (69). Western Ross Sea sedimentary cycles suggest that sea level oscillations only reached ~20-m amplitudes across the EOT, similar to our ΔSL (70). Given literature disagreement about sea level and Tw changes across the EOT, it is not yet possible to either accept or reject our ΔSL result for that event.

Important evidence that is used to argue for large EOT VAIS buildup (or starting 36 Ma ago) concerns an onset of IRD deposition in marine sediments around Antarctica (47, 71), which indicates that ice large enough to support calving had (in places) made it to the coast. On the basis of the size of the continent and the general parabolic cross-sectional profile of ice sheets, this suggests, at first glance, that an ice sheet of roughly modern proportions had built up. However, similar arguments have been rejected for the Early Pleistocene LIS, on the basis of the fact that ice sheets on slippery regolith would remain much lower (because of reduced bed friction) and less voluminous than ice sheets positioned on bedrock after earlier glaciations removed most regolith (7274). We propose that consideration must be given to the possibility that the EOT AIS may have been less like the present-day AIS and more like an Early Pleistocene low-slung, slippery LIS (7275), which might explain IRD deposition even with VAIS only about half its present-day value. Our results suggest that a consistently large-volume/high AIS existed only since about 13 Ma ago (Figs. 4D and 7A).

Fig. 7 Variance analyses.

(A) Reconstructed ice volume in mseq for the southern hemisphere (SH) alone (AIS; light blue) and the global total (AIS plus GrIS, LIS, and EIS; dark blue). NH, northern hemisphere. (B) Sea level variations (gray), with trend line (red; see Methods), detrended sea level variations (green), and variance analysis [standard deviation (SD) per moving 500-ka window; blue]. (C) As (B) but for deep-seawater δ18O variations. (D) As (B) but for deep-sea temperature variations. Pale red, orange, and blue rectangles outline quasi-stable states, and yellow bars indicate transitions, as discussed in the text. Pale purple arrows indicate phases of Tw variance buildup, starting at around the times indicated with dashed purple lines. The question mark for the variance buildup from ~38.2 Ma ago indicates uncertainty in pinpointing the start of this process in records limited to the past 40 Ma. We used as guidance a study that identified an age of ~37 Ma ago for this event (76).

Wider implications

In our reconstructions, before deep northern hemisphere glaciations started 3 to 2.5 Ma ago, mean global deep-sea temperatures remained a few degrees above freezing, which indicates substantial deep-water supply from nonfreezing regions (Figs. 3, C and D, 4D, and 7D). Relative to the Late Pleistocene mean, Tw was 2° to 3°C higher at 3 to 13 Ma ago, 3° to 5°C higher at 15 to 33.5 Ma ago, and 7° to 9°C higher at 34.5 to 40 Ma ago, marking distinct quasi-stable states that remain evident when considering variability around those means (pale red, orange, and blue rectangles in Fig. 7D). These ranges and an absence of outlier patterns in our mutually consistent ΔSL, Δδw, and ΔTw reconstructions (Supplementary Materials) indicate that no major warm, saline deep-water contributions existed throughout the time scales investigated and that high-latitude deep-water formation processes were dominant.

The existence of quasi-stable Tw states is consistent with similar reconstructed ice volume behavior, which suggest virtually no ice before 34 Ma ago, small-to-mid size Antarctic ice at 15 to 33.5 Ma ago, and small-to-full size Antarctic ice with small northern hemisphere ice at 3 to 13 Ma ago (pale red, orange, and blue rectangles in Fig. 7A). These quasi-stable periods are separated by major transitions: the EOT, the end of the Middle Miocene Climatic Optimum (MMCO) and the Plio-Pleistocene transition to extensive bipolar glaciation (yellow bars in Fig. 7D). This resembles critical transition behavior in complex dynamical systems, where variance and autoregression increases occur before transitions (76). Variance analysis (Methods and Fig. 7D) reveals that ΔTw variance increased well before each transition—at ~38.2 Ma ago [previous work suggests ~37 Ma ago (76) and longer time series need to be considered to more precisely identify this onset in our method], 17.2 or even 24.8 Ma ago, and 6 Ma ago (dashed purple lines and magenta arrows in Fig. 7D)—and culminated during the transitions. Given that this is consistent with expectations for critical transitions in complex dynamical systems (76), we explore the potential drivers behind this behavior.

Earth’s climate system has diverse feedbacks that operate over different time scales, including short-term variability/“noise” (77, 78). Insolation, a key external climate forcing, is governed by astronomical cycles that lack multimillion-year trends (79) but that can precondition the system for major state changes at certain times (80, 81). Ice-albedo feedbacks operate over time scales that are too short to cause multimillion-year secular changes, but they will amplify any initial change as part of a feedback cascade, aiding transitions in regimes with multiple stable states. Plate tectonics is a suitable candidate for causing multimillion-year secular atmospheric CO2 changes and hence climate, through subtle but long-lasting modifications of the balance between volcanic CO2 outgassing and weathering-related CO2 sequestration (82). Note that this concerns long-term exogenic CO2 trends related to plate tectonics rather than shorter endogenic CO2 variations related to carbon-cycle feedbacks. CO2 levels changed from ~1000 parts per million (ppm) at around 40 Ma ago to ~180 ppm during Late Pleistocene glacial stages and ~280 ppm during interglacials (8385), which spans roughly two 2× CO2 changes. Given that radiative forcing change for each 2× CO2 change amounts to almost 4 Wm−2 (86), the CO2 reduction since 40 Ma ago amounts to roughly 8 Wm−2 reduction in radiative climate forcing. This represents 5° to 9°C cooling based on a 5 to 95% range for equilibrium climate sensitivity (86), in agreement with the 6°C cooling in our reconstructed Tw (Fig. 7D).

While major climate transitions separating relatively stable states are known from benthic δc compilations (13, 14), we have deconvolved these records into the contributing ice-volume and deep-sea temperature components. Adjusting more rapidly than ice volume, Tw represents a more direct reflection of global climate variability, while ice volume provides a longer time-integrated view of the changing climate state. We propose that long periods of increasing Tw variance (“flickering”) before the detected state transitions (Fig. 7) represent periods when gradually declining CO2 reached levels where multiple stable climate states exist (78, 87), the existing warmer state that supports less ice and a colder state that supports more ice. Stochastic shocks play an important role in triggering transitions before bifurcation points (87). We propose that cold- or warm-biased organization of climate system noise during cold or warm orbital extremes provided these shocks [akin to stochastic or nonlinear resonances; (78, 88)], similar to model-based findings for the EOT (81).

We emphasize that the entire past 40 Ma investigated here is dominated by cooling; apart from the MMCO onset, the interval contains no significant long-term warming. Our evidence, therefore, is largely indicative of transitions in the cooling direction, although inferences can also be made about potential warming transitions. Existence of long (2 Ma or more) flickering intervals preceding the detected transitions (Fig. 7) suggests that bistable regimes existed over long time periods. The eventual transition to a colder stable state with greater ice volume, therefore, likely occurred at a considerably lower CO2 threshold value than the return to a warmer stable state with smaller ice volume because of stabilizing feedbacks associated with large ice masses (89). Similar behavior has been inferred from Antarctic ice sheet modeling experiments, where the hysteresis relative to CO2 levels spans hundreds of parts per million (90). Using the bounds of flickering intervals in our reconstructions (Fig. 7), we infer that bistable regimes existed between at least 0 and 6 Ma ago, between 12.6 and 17.2 Ma ago, and—less well defined in our results—between 34 and 38.2 Ma ago. Precise CO2 reconstructions for those intervals may help to quantify climate state hysteresis relative to CO2, with implications for understanding potential future climate transitions.

Application and outlook

Our mutually consistent, process-based framework for reconstructing ΔSL, Vice, δice, Tw, and δw permits complete system validation using multiple parameters rather than only one or two as in traditional approaches. It also helps to identify unrealistic single-parameter fluctuations and/or to fill in information across data gaps so that more complete records can be obtained. Last, our results offer continuous, quantitative context to data from analytical methods such as sea level benchmarking or novel Tw assessment (e.g., clumped isotopes) with low temporal resolution and/or age uncertainties that preclude precise attribution to specific climate cycles.

Our mutually consistent ΔSL, δice, Tw, and δw records are well validated over at least the past 22 Ma, which comprise the full −130- to +65.1-m range of sea level variation. Validation is not (yet) as convincing for earlier times, so we infer that our ΔSL record (black in Fig. 4A) offers a well-validated sea level record only for the past 22 Ma; this record is obtained using ΔSL = −8.4 δc2 + 6.8 δc + 65.1, where δc is normalized to the most recent value of 3.23‰ in the Lisiecki and Raymo benthic δc stack (13). The two regression extremes (pink and orange in Fig. 4A) provide realistic total range bounds (for equations, see fig. S3). Before 22 Ma ago, our ΔSL results require further validation. Note that our assessment concerns orbital and longer time scale sea level variability and that no conclusions should be drawn with respect to millennial sea level variability without further extensive testing.

We find close amplitude structure agreement between our calculated Tw variations and independently determined Antarctic temperature variations at EDC (54), scaled as ΔTw ≈ 0.25 ΔTEDC (Fig. 3C). This stunning independent agreement provides strong validation of our approach. We infer that global mean deep-ocean temperature has consistently been set by high southern-latitude temperature fluctuations.

Diverse positive validations also indicate that our inferred ΔδwSL relationships offer reasonable approximations of reality. Hence, we propose that ice-volume/sea level reconstruction based on δw records from Tw-corrected δc data, or ice-volume corrections of δw records to identify regional hydrological influences, should no longer use linear transformations. Instead, we propose that improved results will be obtained using the nonlinear relationship found here through all instances (N = 54,644). The optimum regression fit to capture the shape of the data cloud is Δδw = 9.6 × 10−11ΔSL5 + 1.9 × 10−8ΔSL4 + 2.5 × 10−7ΔSL3 − 1 × 10−4ΔSL2 − 0.015ΔSL − 0.133. While polynomial orders higher than cubic do not statistically significantly improve the relationship, the fifth-order fit visually best captures the nonlinearity that is theoretically evident, with the least edge effects, and offers the best conversion equation.

Further research is needed to hone the approach presented here. A particularly promising route is to develop Vicep transforms using geographically specific, less idealized, δ18O-enabled ice sheet models that account for additional controls on δp (26), which can then be applied to the various ice sheet growth histories to determine more realistic δice¯ histories for each ice sheet. We advocate use of these transforms because integrating fully coupled δ18O-enabled ice sheet models over many millions of years will be computationally challenging. Making the suggested transform functions will provide a sound middle way between that and our use of idealized scenarios. In addition, there is a need for further paleoclimate and paleoceanographic benchmarks; for example, to resolve the origins of changes around the EOT and to provide further Oligocene-Miocene data for improved validation. Last, more detailed understanding of Mg/Ca changes over the past 40 Ma would improve deep-sea temperature reconstructions and facilitate more advanced validation of the approach presented here.

Our data-driven reconstructions suggest that climate-system changes during the descent from Eocene greenhouse to Late Pleistocene icehouse conditions were not gradual but comprised sharp threshold-like adjustments (critical transitions) between quasi-stable states. By analogy, it seems important for future climate change projections to allow for historically unprecedented adjustments, such as deep-water formation shifts, abrupt (partial) ice sheet collapse, or other cascading feedbacks.


ΔSL from δc records

Calculation of sea level change (ΔSL) from benthic δc records (13, 14) was undertaken using second-order polynomial regressions based on the lag-optimized comparison between δc and a statistical sea level assessment based on a variety of input data (12). The regressions and their equations are given in fig. S3. Uncertainty bounds account for both the regression uncertainty and additional data uncertainties in both X and Y directions. For the latter, we used normal distributions with 1σ = 0.1‰ for δc and 1σ = 2 m for ΔSL.

Isotope fractionation over the ice sheets

δp is determined in the main case of the analytical assessment using Eq. 1 and for sensitivity test 1 using Eq. 2. Constants were selected to achieve a reasonable visual fit through modern ice sheet values (fig. S1). For the main case of our analytical assessment and the subsequent applied assessments, we use a heuristically fitted function representative of the described trends that start with an initial snow δp = −15‰ that drops rapidly with initial ice sheet volume buildup, followed by an exponential change with V toward very negative values over large-scale ice sheets (fig. S1). To this end, we use an equation of the formδp=αe{(0.9V10)63}2(100V)(C+15)(1)

Here, C is a constant that is adjusted along with amplitude parameter α so that initial precipitation (for ice sheet volume V = 0 mseq) is precisely −15‰ in all cases. Extreme Rayleigh distillation at the isolated, high-latitude AIS causes very negative δp (42), which, in Eq. 1, requires α = 1 (with C = 102.55). GrIS is set to the same parameters. In contrast, a lower α is needed for LIS (and EIS) to account for less negative δp in response to more limited Rayleigh distillation between evaporative sources and precipitation sites, related to their positions in lower latitudes closer to relatively warm oceanic moisture sources. An α = 0.55 value (with C = 56.4) is used for LIS to reach a minimum LGM δp of roughly −35‰, as suggested by five isotope-enabled global climate circulation models (35). EIS is set to the same parameters.

Sensitivity test 1 instead uses linear δp changes (Eq. 2), although this is physically implausible. For AIS and GrIS, the slope-defining constants J and K are set to −60 and 66.3, respectively, and for LIS and EIS, these are set to −25 and 26, respectively, to get reasonable agreement with δp observations (fig. S1)δp=VJ15K15(2)

In sensitivity test 2, AIS glaciation is set to an arbitrary quadratic function so that AIS grows less rapidly during initial glaciation stages and more rapidly during later stages. The sea level budget is closed by setting an opposite LIS-volume anomaly. This is explored to evaluate impacts on the solution, which are found to be negligible.

Ice sheet and seawater δ18O

The global volume-weighted (here equivalent to mass-weighted) mean oxygen isotope composition of ice is given byδice*=δice_AISVAIS+δice_GrISVGrIS+δice_LISVLIS+δice_EISVEISVtot(3)

The imposed global seawater δ18O change is found usingΔδw=0.9 Vtotδice*Voc+ΔSL(4)

Here, 0.9 is the ice-to-water density ratio, Voc is global ocean volume (3700 mseq), and ΔSL = −{(VAIS − 57.8) + (VGrIS − 7.3) + VLIS + VEIS}. Note that, throughout, we consider ice-volume change to be inversely proportional to sea level change, accounting for density differences between water and ice. This is a simplification because some ice sheet portions may displace water so that they do not contribute to sea level change, whereas they do count toward ice-volume (and thus δw) change [e.g., (37)].

Individual ice sheet volume histories

In our applied assessments, ΔSL (in meters relative to present) is used to determine global net ice loss or gain per time step j (forward in time). Ice sheets are grown in relation to ΔSL using prescribed conditional functions. Note that we do not fully reduce the functions below, so it will be clearer what is calculated. AIS volume per time step is given byVAISj=|57.8+ΔSLj2 if 0<7.3+ΔSLj27.357.8ΔSLj if 7.3<ΔSLj57.80 if 57.8<ΔSLjVAISj1+zmin12515(ΔSLjzmin)2 otherwise(5)

Here, 57.8 stands for modern VAIS (in mseq), 7.3 for modern VGrIS (in mseq), 15 for the assumed LGM VAIS addition (in mseq), ΔSL for sea level position per time step (in meters relative to present), zmin for the maximum sea level drop in the ΔSL record used, and VAISj−1 at run initialization is simply set at modern VAIS (hence, all runs have an initialization period during which the solution stabilizes). The squared term in the last row provides a nonlinear AIS-excess buildup during further glaciation, relative to present, to approximate observations of less excess-ice growth at earlier stages and greater growth at later glacial stages (38).

GrIS volume is calculated similarly (with assumed LGM VGrIS addition of 5 mseq), asVGrISj=|7.3+ΔSLj2 if 0<7.3+ΔSLj27.30 if 7.3+ΔSLj20VGrISj1+(zmin1255ΔSLjzmin) otherwise(6)

For LIS and EIS combined, which exist only when ΔSL < 0 m, linear growth is determined in proportion to the ΔSL record. However, between them, we assume that one dominates growth in initial glaciation stages, with the other catching up at later stages, to roughly mimic moisture scavenging effects between ice sheets due to large-scale atmospheric circulation changes (91). This is done in a similar fashion to sensitivity test 2 of the analytical assessment with an arbitrary nonlinear gradient function, calculated for each time stepGj=1(0.5e)20ΔSLjzmin(7)

EIS volume is simply calculated taking into account this gradient function, and VLIS in addition compensates for the nonlinearity stipulated in VAIS at ΔSL < 0 m, to ensure that total ice-volume changes remain balanced with ΔSL (within rounding errors).

To calculate VLIS, there is a first stepVxj=|(zmin12570ΔSLjzmin)Gj(zmin12515ΔSLjzmin+zmin12515(ΔSLjzmin)2) ifΔSLj>00 otherwise(8)

Here, 125 mseq stands for the total ice-volume increase corresponding to the assumed LGM ΔSL, 70 mseq for the LIS component (35 mseq for EIS), and 15 mseq for the LGM expansion of VAIS, as before. Then, VLIS is given byVLISj=|Vxj if 0Vxjzmin125700 if Vxj<0zmin12570 otherwise(9)

To calculate VEIS, the first step isVsj=|ΔSLjzmin125(70+35)zmin(zmin12570ΔSLjzmin)Gj ifΔSLj>0ΔSLjzmin125(70+35)zmin otherwise(10)

In addition, VEIS is then given byVEISj=|Vsj if 0Vsjzmin125350 if Vsj<0zmin12535 otherwise(11)

Ice sheet weighted mean δ18O and seawater δ18O

For each ice sheet, the isotopic value of new precipitation (δp) is determined per time step using Eq. 1, for AIS and GrIS using α = 1 (C = 102.55), and for LIS and EIS using α = 55 (with C = 56.4). To calculate the volume-weighted mean δice (i.e., δice¯) for each ice sheet, we use a residence-time calculation. This starts with ice sheet gross accumulation rate (agrs) determination based on the area defined by its volume (with constant aspect ratio ε), simply scaled in proportion to a modern global annual agrs of ~0.008 mseq (44, 45) for a global ice volume of 65.1 mseq. Volume is obtained in relation to ΔSL according to Eqs. 5 to 11; volume is true volume (in m3) rather than scale volume (in mseq); this is indicated using V. The term ΔSL(LIS) is used to indicate the sea level change component attributed to the different ice sheets (here LIS). For a world ocean area Aoc of about 362 × 1012 m2, this givesVLISj=VLISjAoc0.9(12)

For the vertical area under a parabolic profile (with maximum height hmax, which is ε times radius r), the mean height is (2/3)hmax or (2/3r so that the volume of our standard planoconvex lens–shaped ice sheet with parabolic cross section is V = (2/3)ε π r3. Solving for r gives the information required to approximate the ice sheet area over which accumulation occurs (here, the example is given for LIS, but it is the same for all others)ALISj=π(32VLISjAoc0.9επ3)2(13)

Note that this concerns simply the surface area of the circle occupied by the ice sheet, ignoring the relevant (half) arc length of the parabola because the aspect ratio is so small that both numbers are virtually identical; hence, the circle area was chosen for computational efficiency. The rate agrs over that ice sheet at time step j then is agrsj = 0.008ALISj/Atot( mod ), for a modern global annual agrs = ~0.008 mseq when global ice volume is 65.1 mseq. Here, Atot(mod) is the surface area given by Eq. 13 for the modern total global ice volume of 65.1 mseq (10.7 × 1012 m2). The model calculates the isotopic contribution from that gross accumulation by its product with δp based on Eq. 1. At every time step, gross accumulation is balanced by a gross loss term, which is simply set to agrs minus the sea level change contribution from one time step to the next of the ice sheet considered (here, anetj = ΔSL(LIS)j − ΔSL(LIS)j−1). The model determines the isotopic impact of anetj using δice¯j1. The calculation is conditional in that periods with negative anet do not lead to changes in δice¯δice¯j=|δice¯j1 if anetj<0agrsjδpj(agrsjanetj)δice¯j1+VLISj1δice¯j1VLISj1+anetj otherwise(14)

Last, a constraint is used that when δice¯ from Eq. 14 drops below the model start value δp = −15‰ for ice growth, and then δice¯ is kept at −15‰. This is done to avoid spurious effects of rounding errors at extremely small ice volumes. We now have, per time step, ice volumes for each ice sheet and their weighted mean ice sheet δice¯. Thus, we solve Eq. 4 per ice sheet, per time step, to calculate each ice sheet contribution to a change in mean ocean δw using the world ocean volume Voc = 3700 mseq, giving (here for LIS)Δδw(LIS)j=0.9VLISjδice¯(LIS)jVoc+ΔSL(LIS)j(15)

Executed for each ice sheet, their cumulative effects can be determined per time step. The cumulative Δδw, normalized to the present-day value, can be subtracted from δc changes, giving a δc residual that indicates bottom-water temperature (Tw) changes based on a change of −0.25‰ per degree Celsius warming.

Thus, we have complete mutual consistency between all key parameters: ΔSL, ice sheet volumes, ice sheet weighted mean δice¯, global Δδw, and Tw. Mutual consistency means that if one changes, so must others. Hence, validation of model results can be undertaken using multiple parameters.

Variance analysis

For the variance analysis, we first detrended the records obtained using the benthic δc megasplice (14) with a 500-ka moving Gaussian smoothing window. We then calculated SDs of detrended records of residuals over 500-ka windows that were shifted in 100-ka steps across the record. For robustness, we consider only major signals revealed by this analysis. Autoregression analysis (76) is precluded because of the (partially) interpolated and/or smoothed nature of the time series considered.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Acknowledgments: Funding: This study was supported by U.K. Natural Environment Research Council grant NE/P011381/1 (to G.L.F. and E.J.R.). This study also contributes to Australian Research Council Discovery Project DP2000101157 (E.J.R., D.H., and G.L.F.). Author contributions: E.J.R. designed the project, developed the method, and led the write-up. All authors contributed to refining the method and validation, discussing results, and write-up. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are presented in the paper and/or the Supplementary Materials and will be archived at both and

Stay Connected to Science Advances

Navigate This Article