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

See allHide authors and affiliations

Science Advances  15 May 2020:
Vol. 6, no. 20, eaaz1346
DOI: 10.1126/sciadv.aaz1346


Using Pacific benthic foraminiferal δ18O and Mg/Ca records, we derive a Cenozoic (66 Ma) global mean sea level (GMSL) estimate that records evolution from an ice-free Early Eocene to Quaternary bipolar ice sheets. These GMSL estimates are statistically similar to “backstripped” estimates from continental margins accounting for compaction, loading, and thermal subsidence. Peak warmth, elevated GMSL, high CO2, and ice-free “Hothouse” conditions (56 to 48 Ma) were followed by “Cool Greenhouse” (48 to 34 Ma) ice sheets (10 to 30 m changes). Continental-scale ice sheets (“Icehouse”) began ~34 Ma (>50 m changes), permanent East Antarctic ice sheets at 12.8 Ma, and bipolar glaciation at 2.5 Ma. The largest GMSL fall (27 to 20 ka; ~130 m) was followed by a >40 mm/yr rise (19 to 10 ka), a slowing (10 to 2 ka), and a stillstand until ~1900 CE, when rates began to rise. High long-term CO2 caused warm climates and high sea levels, with sea-level variability dominated by periodic Milankovitch cycles.


Cenozoic sea level

Sea-level change is one of the most obvious manifestations of climate change, both in the past as the Earth evolved from substantially ice-free hothouse/greenhouse worlds (e.g., Late Cretaceous to Eocene) to largely glaciated icehouse worlds [e.g., past 34 million years (Ma)] (1) and in the present-day world, where anthropogenically driven global mean sea level (GMSL) rise has accelerated from a rate of ~1.4 mm/year in the 20th century to an early 21st century rate of 3.1 ± 0.4 mm/year (2, 3). GMSL is the global mean of relative sea level (RSL) change and changes in response to changes in (i) the density of water in the ocean, due to changes in ocean temperature (primarily on decadal-centennial time scales); (ii) the mass of water in the ocean, due to changes in land ice volume and terrestrial water storage (operating on decadal-centennial and on astronomical times scales of ~20 ka to 2.4 Ma); and (iii) ocean basin volume [primarily on >1-Ma scale; summaries in (4)]. Geologists often use the term “eustasy” for GMSL with respect to some fixed datum level such as the center of the Earth (more properly termed global mean geocentric sea-level change) (5), but these datum levels are lacking or equivocal in the geologic record, and we avoid this term.

RSL has been defined in many ways, but we follow the definition of RSL as the difference in height between the sea surface and the solid Earth (5, 6). RSL change at a particular place is controlled both by GMSL change and by regional and local effects that include subsidence or uplift, including the effects of thermal subsidence, sediment loading, flexure, glacial isostatic adjustment (GIA) (6), mantle dynamic topography (MDT) (7), changes in the height of the geoid driven by the changing distribution of ice and ocean mass, and short-term (10- to 1000-year scale) ocean dynamic sea-level changes (6). RSL changes combine with the effects of sediment supply to produce regressions (seaward movement of the shoreline due to RSL fall and/or increased sediment input) and transgressions (landward movement of the shore, due to RSL rise and/or decreased sediment input). These, along with unconformities associated with RSL falls, are the predominant expression of sea-level change in the stratigraphic record (8). Extracting a GMSL signal from the stratigraphic record is challenging due to overprinting by regional and local effects, uncertainties in sea-level proxies, and limitations in chronological control.

The Cenozoic era (last ~66 Ma) exhibited changing CO2, climate, and sea levels, providing tests of interrelationships, especially during the last 50 Ma, when long-term global climate shifted from hothouse/greenhouse conditions, typified by high atmospheric CO2 and ice-free conditions, to fully glaciated icehouse conditions with low CO2 (1). Hothouse climates largely lacking ice sheets occurred during the Early Eocene (56 to 47.8 Ma), with atmospheric CO2 reaching levels higher than 800 parts per million (ppm) and maintaining a level more than two times higher than preanthropogenic levels of 280 ppm (Fig. 1) [e.g., (9)]. Previous studies suggested that a largely unipolar icehouse world began in the earliest Oligocene with continental-scale ice primarily in Antarctica paced by long obliquity (1.2 Ma) changes, although there is evidence for large-scale glaciation in the Eocene (1013), a topic we highlight here. A bipolar icehouse began ca. 2.7 to 2.55 Ma, at the beginning of the Quaternary, with the development of large, northern hemisphere (NH) ice sheets during the “ice ages” (14). This general view of cryospheric evolution is often debated or misunderstood; for example, the 2.55-Ma event has been dubbed the “initiation of NH glaciation,” despite evidence for NH glaciation extending to the Eocene (15). Here, we synthesize published benthic foraminiferal oxygen isotopic (δ18Obenthic) and Mg/Ca records and compare to previously published backstripped estimates to provide an updated overview of Cenozoic cryospheric and sea-level history. This paper is novel in providing the first δ18Obenthic splice for the entire Cenozoic, a detailed (Ma scale) sea-level record for the last 48 Ma based on the δ18Obenthic-Mg/Ca approach (Mg/Ca records older than 48 Ma are uncertain), and a comparison of two independent estimates of sea level, δ18Obenthic-Mg/Ca and backstripping.

Fig. 1 Summary of Cenozoic benthic foraminiferal δ18O, temperature, sea level, and CO2 records.

(A) Benthic foraminiferal δ18O splice reported to Cibicidoides spp. values. Modern is the core top value for δ18OCibicidoides; LGM is the LGM δ18OCibicidoides value; the icehouse line is placed at 1.8‰ in Cibicidoides, with values greater than 1.8‰ requiring major ice sheets (see text); vertical line at −0.5‰ separates cool greenhouse from hothouse worlds. Each site has a contrasting shade of blue, with the site number indicated on the left column in the corresponding color. (B) Temperature estimates and error bars from (33) based on Mg/Ca values in Pacific cores using their equation 7b. Blue is the smoothed record, removing periods shorter than 2 Ma. Black lines are 2σ error estimates. Red line is temperature assumed in Paleocene-Early Eocene. (C) Sea level obtained using our new benthic foraminiferal δ18O splice and δ 18Oseawater using temperatures from Mg/Ca (30); ice-free line (magenta) drawn 66 m above present, GIS-WAIS (light blue) drawn at 12 m above present, and Laurentide (dark blue) drawn at −50 m. Each site has a contrasting shade of blue, with the site number indicated on the left column in the corresponding color. Error bars as stated in the text and Materials and Methods. (D) Comparison of smoothed sea-level estimates from δ18O and Mg/Ca (blue; obtained by interpolating to 20-ka intervals and using a 49-point Gaussian convolution filter, removing periods shorter than 490 ka) with backstripped onshore NJ estimates [red; (4, 23)]. Magenta bar is the range of sea-level change on the Marion Plateau, East Australian margin (17). The NJ estimates for the Early to Middle Eocene were shifted by −50 m. Error bars as stated in text and Materials and Methods. (E) CO2 compilation and error bars of Foster et al. (47) and a compilation of measurements from Antarctic ice cores (98) fit with Bayesian hierarchical model that relates the raw CO2 data to geological time. The thick black line is the prediction of CO2 derived from the model, and the dashed and dotted lines demarcate the 1σ and 2σ confidence intervals for that prediction. Red, ice cores; purple, boron; blue, stomata; yellow, alkenone; gray, liverworts; green, paleosols. Inset map shows locations of Pacific sites and New Jersey for backstripped record.

Our view of Cenozoic sea-level changes has been influenced by seismic and sequence stratigraphic (including outcrops, wells, and well log) studies of Exxon Production Research Company (EPR) (8, 16). The EPR cycle charts provide a record of the timing of lowerings of sea level that have been tested and largely verified for portions of the Late Cretaceous to Miocene by academic drilling on the New Jersey, New Zealand, and Australian continental margins and the Bahamas and Maldives carbonate platforms (4, 1720). Although the timing of the Late Cretaceous to Miocene EPR curves appears to be approximately correct, their extremely high sea-level amplitudes (e.g., 160-m middle Oligocene event) are too high by a factor of 2 to 3 (see Discussion) (4, 17).

Extracting amplitudes of GMSL change from continental margin sequences is challenging, complicated by the effects listed above. “Backstripping” attempts to progressively account for the effects of sediment compaction, loading, and thermal subsidence on passive margin sections (21). After accounting for these effects, the residuals (“R2”) are an estimate of changes in GMSL, the gravitational and rotational effects of mass redistribution (primarily those due to GIA), and nonthermal subsidence or uplift, including both those associated with GIA and MDT (22, 23). Ocean dynamic sea-level changes are below the vertical and temporal resolution of backstripped records. Drilling of the middle Atlantic U.S. passive margin by International Ocean Discovery Program (IODP) and its predecessors (Legs 150, 150X, 174A, 174AX, and Expedition 313) has provided cores, downhole logs, and seismic profiles that were used to identify and backstrip Late Cretaceous to Miocene sequences and yield R2 estimates (4, 22, 23). The amplitudes of these R2 changes in the middle Atlantic U.S. Coastal Plain are on the order of 10 to 60 m; backstripping of the Miocene of the Marion Plateau, northeast Australia (17) (Fig. 1), verifies these amplitudes (Fig. 1) (17), contradicting the extremely high amplitudes for coeval events on the EPR curves (see Discussion). The middle Atlantic U.S. R2 record has also provided the first testable sea-level curve for the Late Cretaceous to Miocene because the data are in the public domain (23, 24).

Although the middle Atlantic U.S. R2 sea-level curve (Fig. 1C) provides a constraint on GMSL changes, it is overprinted by nonthermal tectonism. The middle Atlantic United States is a classic passive margin, dominated by thermal subsidence and loading; however, regional differences on this (23, 25) and other passive margins are caused by tectonism, although the mechanisms have been poorly known. The realization that MDT influenced surface topography on this margin [e.g., (7)] provides a plausible mechanism for the observations of differential uplift and subsidence [e.g., (23, 25)]. The role of mantle dynamic flow affecting vertical motion of the crust has been known for some time (26), but advances in high-resolution seismic tomography now allow modeling of these changes and reconstruction of past topography (27). Changes in MDT explain the vertical movement of continents that influence sea-level records but then question the significance of margin sea-level records. Another sea-level proxy is needed to sort out the effect of GMSL versus nonthermal tectonism, one afforded by benthic foraminiferal δ18O records.

δ18O and Mg/Ca syntheses

Foraminiferal δ18O records reflect variations in the temperature and δ18O of seawater (δ18Oseawater) during calcification. Changes in δ18Oseawater are due primarily to global changes resulting from growth or decay of continental ice sheets and regional and surface ocean changes due to Raleigh fractionation associated with evaporation and precipitation (“salinity effect”). These two unknowns (temperature and δ18Oseawater) are expressed in the paleotemperature equation (28)T=16.14.64 (δ18Ocalciteδ18Oseawater)+0.09 (δ18Ocalciteδ18Oseawater)2(1)

Deep-sea benthic foraminiferal oxygen isotopic records (δ18Obenthic) minimize the effects of evaporation, precipitation, and seasonal temperature variations. Deep-sea δ18Obenthic records reflect high-latitude temperatures (where deep and bottom waters originate) and ice-volume changes (29, 30). However, an independent temperature proxy is required to estimate δ18Oseawater changes caused by ice-volume variations.

Mg/Ca ratios in benthic foraminiferal tests (Mg/Cabenthic) provide an independent proxy for deep ocean temperature [e.g., (3133)]. We obtain an estimate of Cenozoic deep-water temperature changes using a published synthesis of benthic foraminifera Mg/Cabenthic records from the Pacific that was smoothed (>2 Ma periods; see Materials and Methods) (33). Application of this method is limited by the fact that the Mg/Ca ratio of seawater (Mg/Caseawater) has varied on the Ma scale [reviewed in (33)]. In addition, variations in seawater carbonate ion saturation state affect Mg/Cabenthic; for example, an abrupt deepening of the calcite compensation depth (CCD) during the Eocene-Oligocene transition (EOT) explains an apparent warming artifact in deep ocean coincident with glaciation of Antarctica (31, 33, 34). Cramer et al. (33) reviewed uncertainties in using Cenozoic Mg/Ca temperatures to obtain δ18Oseawater. They concluded that sea-level estimates derived from the δ18O-Mg/Ca method statistically agreed well with those from middle Atlantic U.S. backstripping, but their analysis was largely relegated to periods >2 Ma. Here, we apply their long-term paleotemperature estimates to ka scale sampled δ18Obenthic records to interrogate sea-level change primarily on time scales of Ma and shorter.

Previous studies have shown that astronomical (“Milankovitch”) forcing of ice sheet growth and decay causes beats in δ18Obenthic records with distinct periods (19/23-ka precession, 41-ka tilt, 95- and 125-ka short eccentricity, 405-ka long eccentricity, and 1.2-Ma long tilt) (13, 3538). The major Oi and Mi Ma-scale δ18Obenthic events of the Oligocene-Miocene (12) appear to be paced by the 1.2-Ma tilt cycle (13), implying a sea-level response at these periods. Here, we seek to constrain sea-level variations on the Ma scale using paleotemperature equation (Eq. 1) and the δ18O-Mg/Ca paleothermometer to estimate bottom water temperatures, solving the paleotemperature equation for δ18Oseawater, a proxy for ice-volume changes. Our strategy is to use δ18Obenthic from Pacific deep waters, which comprise 60% of the modern and 80% of the early Cenozoic global reservoir and that has apparently had one major southern source of deep water through the Cenozoic (38), minimizing the effects of temperature and “salinity” (evaporation-precipitation mixing) seen in Atlantic records due to deep circulation changes. No single site contains a complete Cenozoic record, so we splice together the most complete records of δ18Obenthic from Pacific sites.

We present a new splice of benthic foraminiferal records for the past 66 Ma similar to that in (37) for the past 35 Ma but composed entirely of Pacific records (enlarged map of sites fig. S1) and extending through the entire Cenozoic (Fig. 1A and fig. S2). Previous Cenozoic δ18Obenthic syntheses focused on stacks of records averaging multiple coeval records, but problems with interoceanic variations and inherent problems in stacking and smoothing limit their applicability to estimating sea-level amplitudes. Previous syntheses were either restricted to the Atlantic (12, 38), where δ18O variations may be overprinted by deep water circulation change, or used records from different ocean basins (39), with resulting artifacts in the stack. Cramer et al. (40) compiled data by ocean basin, and their stacks still provide the best overview of Cenozoic deep-water δ18O variations globally and by ocean basin, but like any record that is stacked and smoothed, it shows dampened amplitudes. Recent publication of δ18Obenthic values from the Pacific provides a remarkable record of mostly monogeneric Cibicidoides samples from a single site for given time intervals that allows us to compile the first astronomical-scale splice for the entire Cenozoic [Fig. 1A and fig. S1; Paleocene to Middle Eocene, site 1209 (41); late Middle Eocene to Early Miocene, site 1218 (36); Miocene, sites 1146, 1337, 1338 (4244); and Plio-Pleistocene, site 846 (45) and V19-30 (46)].

Here, we provide a Milankovitch-scale resolution δ18Obenthic spliced record back to ~66 Ma, using smoothed Mg/Ca paleotemperatures to derive an estimate of δ18Oseawater. We scale our δ18Oseawater to sea level (Fig. 1B). Our record has large uncertainties before 48 Ma due to uncertainties in Mg/Caseawater (33) and does not account for variations in temperature at Milankovitch scale. Nevertheless, it provides an excellent chronology of sea-level changes with an uncertainty of ±10 to ±20 m for the last 48 Ma, with larger errors before this (see Materials and Methods for discussion of errors). We compare the δ18O-Mg/Ca GMSL estimate with the middle Atlantic U.S. backstripped estimate (Fig. 1C), compare it to a compilation of Cenozoic CO2 estimates (Fig. 1D) (47), and explore the implications of our GMSL estimates for the Cenozoic evolution of the cryosphere.


The δ18O splice and Earth’s thermal states

It has been previously recognized that there was more than one climate state in greenhouse worlds, with warmer (“hothouse”) and cooler (“cool greenhouse”) modes (48). Benthic foraminiferal δ18O values in the spliced record can be used to differentiate hothouse (largely if not entirely ice free), cool greenhouse (with ephemeral ice sheets), and icehouse (those with a continental-scale ice sheet at one or both poles) climates. We use δ18OCibicidoides values to differentiate hothouse (values <−0.5‰), cool greenhouse (values −0.5 to 1.8‰), and icehouse (>1.8‰) worlds [see also (38, 48)]. We show ice sheets in the cool greenhouse of the Middle to Late Eocene (see Discussion), with largely ice-free conditions of the Hothouse Early Eocene.

Values of δ18OCibicidoides greater than ~1.8‰ beginning in the Oligocene indicate the presence of large ice sheets and development of an icehouse world (49). These high δ18O values indicate either bottom water temperatures colder than today (incompatible with an ice-free world) or the presence of an icehouse world (12). δ18OCibicidoides values did not reach modern values until the Miocene, and bottom waters were warmer than present (<1°C) until the late Neogene (e.g., 8°C in the Oligocene to early Middle Miocene cooling to 5°C in the Middle Miocene; fig. S1). Bottom water temperatures in the cool greenhouse (Middle to Late Eocene; 5 to 10°C) and the early icehouse (Oligocene–early Middle Miocene; 5 to 8°C) are unexpectedly warm (Fig. 1B). These warm deep-water temperatures (Fig. 1) occurred during the early icehouse when ice volumes exceeded modern (Fig. 1), refuting early expectations that warm deep waters indicating warm surface waters in the Antarctic source regions precluded large ice sheets in the Oligocene (29). Nevertheless, despite warm bottom waters, ice sheets likely extended to the margins of the Antarctic continent during the early icehouse, with ice volumes at times greater than modern (Fig. 1). The high δ18OCibicidoides values attained at the Last Glacial Maximum (LGM; ~4.7‰; ca. 20 to 27 ka) were not reached until after the middle Pleistocene transition (MPT; ca. 800 to 1100 ka) and were restricted to the largest glacials [marine isotope stages (MIS) 2, 4, 10, 12, and 16].

δ18O splice, Ma-scale climate, and Milankovitch beats

Our δ18Obenthic splice faithfully records the major known Ma-scale features of the Cenozoic record (Fig. 1A and fig. S1). A minor δ18O decrease (warming) occurred across the Cretaceous/Paleogene boundary (66.05 Ma). A period of stable δ18O values occurred in the Early Paleocene, and a δ18O increase (cooling) occurred across the Early/Late Paleocene boundary, reaching maximum values ca. 59 Ma (Fig. 1). Long-term warming from the Late Paleocene to Early Eocene (59 to 54 Ma) was punctuated by the Paleocene-Eocene thermal maximum (PETM) and numerous Early Eocene hyperthermals (Fig. 1). Minimum Ma-scale δ18O values during the Early Eocene climatic optimum (EECO; ca. 55 to 48 Ma) (39) were followed by a ~2‰ long-term δ18O increase in the Middle Eocene (ca. 48 to 40 Ma), punctuated by low δ18O values during the Middle Eocene Climate Optimum (MECO; ca. 40.5 to 40.0 Ma) (50). Late Eocene δ18O values were relatively stable (~1.0‰) except for a sharp 1‰ decrease and minimum at ca. 35 Ma, followed by the EOT and the beginning of the icehouse world with the Oi1 δ18O maximum (12).

Our icehouse (past 34 Ma) δ18O splice is similar to the “megasplice” in (37) because they are based on many of the same sites (Pacific sites 1146, 1337, 1338, and 1218) (Fig. 1A). However, Vleeschouwer et al. (37) also relied on Atlantic data (e.g., last 9 and 20 to 24 Ma), did not extend through the Paleogene, and focused on shorter (405 ka and shorter Milankovitch) periods. We see much less of a Late Oligocene δ18O decrease than that shown in (38), which we view as largely an artifact of combining Atlantic and Pacific data in their splice, similar to the Late Oligocene warming artifact in the stack in (39).

The EOT culminated in the earliest Oligocene δ18O maximum (Oi1) associated with the first continental-scale Antarctic ice sheets. It was followed by a general “middle” Oligocene δ18O maximum (Oi2 and its subdivisions) and the Mi1 δ18O maximum spanning the Oligocene/Miocene boundary (12). Early Miocene δ18O maxima (e.g., Mi1a, Mi1b, and Mi2) were paced by the 1.2-Ma tilt cycle (Fig. 2) (13) but are less prominent in Pacific records (Fig. 1A) than in Atlantic and Indian Ocean records. These events are real features of the temperature and ice-volume history but are amplified by deep-water temperature changes in the Atlantic. The Miocene climate optimum (MCO) is generally recognized in δ18Obenthic records by minimum Pacific values from 17.0 to 14.8 Ma (Fig. 1) (42).

Fig. 2 CWT spectrograms of forcing, δ18O, and sea level for the Cenozoic.

The lower white line is cone of influence above which data are comparable. The color scales of relative power are used to map the temporally localized amplitudes associated with the spectral components shown in the plot. The upper white solid line in (C) and (D) is the Nyquist frequency, and periods shorter than this are not discernible. Precession (19 and 23 ka), eccentricity (95, 125, 405, and 2400 ka), and tilt (41 and 1200 ka) periods are indicated. (A) CWT of La2004 full solution (97) for solar insulation at 65°N on June 21. (B) CWT of La2010 solution (86) for the eccentricity of the Earth-Moon barycenter. (C) CWT of the δ18O splice (Fig. 1A). (D) CWT of δ18O-Mg/Ca GMSL (Fig. 1B).

Even as the gradual Middle Eocene and sharp EOT δ18O increases marked changes in state, the Middle Miocene was also a change in state, as shown by spectral characteristics of the splice (Fig. 2). Previous authors (37) showed that, during the Oligocene to early Middle Miocene, δ18O variations were dominated by precession, tilt, and short eccentricity (95/125 ka) forced by Southern Hemisphere summer insolation, amplified by a dynamic Antarctic ice sheet. We note the same precession, tilt, and short eccentricity periods in our splice; we also looked at longer periods than (37) and show significant power in the >1-Ma range, with apparently strong long tilt (1.2 Ma) and very long eccentricity (2.4 Ma) forcing (Fig. 2). During the Early to Middle Miocene, strong quasi–100-ka periods recurred (Fig. 2) as noted previously (e.g., a 100 ka–dominated world has been identified in δ18O records from 19.2 to 20.8, 16.2 to 14.5, and 18.7 to 19.2 Ma). During the three stepwise δ18O increases of the middle Miocene climate transition [MMCT; steps at 14.8 (Mi2a), 13.8 (Mi3), and 12.8 Ma (Mi4)], physical evidence indicates that the East Antarctic ice sheet (EAIS) grew to near its present state and remained permanent (51) and the EAIS became too large to be strongly influenced by precession (37). As a result, late Middle Miocene to Early Pliocene δ18O variations were relatively small and paced by the short tilt (41 ka) cycle (Fig. 2) (37). During the Middle Miocene to Early Pliocene, large, Ma-scale Oi1 and Mi variations were muted, with minimal 2.4- and 1.2-Ma modulations, although 41-ka beat appears strong. The Early Pliocene was an interlude in the global cooling trend. Relatively high δ18O values occurred during the Pliocene climate optimum (PCO), followed by increasing δ18O values from 2.7 to 2.55 Ma, associated with the development of large NH ice sheets.

Cenozoic sea-level overview

Our sea-level record derived from δ18O-Mg/Ca has an “ice-free line” set at 66 m, the uncompensated change in water volume stored in continental ice sheets and glaciers (52, 53). The isostatic compensation due to water loading from melting all ice has often been estimated using Airy assumptions of a whole-ocean loading (a 66-m loading of water results in ~44 m of RSL rise if Airy compensation is assumed). The Airy assumption is not warranted for whole ocean compensation, and a whole Earth model is needed to compute the GMSL change with respect to a datum level. However, here, we are primarily concerned with measuring the ice-volume effect alone, but not its density (5), so placing an ice-free GMSL line of 66 m is a reasonable limit for ice-free conditions. Similarly, the modern Greenland ice sheet (GIS) and West Antarctic ice sheet (WAIS) contribute about 12 m of sea-level equivalent (52, 53), and sea levels above 12 m likely involve melting portions of the EAIS, yielding our “Greenland-WAIS” line. Last, we use the level achieved during MIS 100 at 2.4 Ma of −50 m as indicating significant Laurentide glaciation (“Laurentide” line; Fig. 1B).

Our overview (Fig. 1 and figs. S3 to S6) suggests large ice volumes (>20-m sea-level equivalent) during the Cenozoic except for much (but not all) of the EECO, MECO, the latest Eocene (ca. 35 Ma), and peak MCO interglacials. Laurentide ice volumes are not inferred until the Quaternary (<2.55 Ma), although continental ice storage larger than modern apparently occurred during glacial intervals of the Oligocene–Early Miocene and the late Middle to early Late Miocene (Fig. 1), and storage of NH ice sheets is implied at these times. Paleocene ice inventories as large or larger than modern are likely an artifact of poor temperature constraints due to uncertain Mg/Caseawater at that time (Fig. 1 and fig. S3; see Discussion) (33).

Although the δ18O-Mg/Ca sea-level curve (Fig. 1) has the same expected periodicities found in the δ18OCibicidoides record (Fig. 2), continental margin strata generally do not fully record higher frequency precession, short tilt, and short eccentricity sea-level changes, acting as a low-pass filter for sea-level changes. We smoothed the δ18O-Mg/Ca sea-level record using a Gaussian filter, removing periods shorter than 490 ka, and show that they are similar to R2 sea-level estimates obtained from the middle Atlantic U.S. margins (Fig. 1C; see statistical comparison in Materials and Methods). Both show a lowering of GMSL 40 to 50 m above present in the Eocene to near or below modern levels in the Oligocene, reflecting the development of a large EAIS. The timing of sea-level events is similar in both the δ18O-Mg/Ca and backstripped estimates, although the latter often do not record the lowstands; for example, the Middle Miocene backstripped onshore NJ record only captures the highstands (23), although these are captured in more complete offshore NJ records (fig. S5). Both sea-level estimates show well-developed, Ma-scale cycles in the Oligocene to Early Miocene (Fig. 1C and fig. S4). The agreement in timing differs by 1 Ma at Mi1 (Oligocene/Miocene boundary) but is otherwise excellent in the Middle Eocene to Miocene, at least on the Ma scale (see Materials and Methods, Fig. 1C, and fig. S7). Agreement between the two methods is poorer in the Paleocene to Early Eocene, and further studies of these intervals are needed (Fig. 1C and fig. S3).

We caution that Paleocene to Early Eocene δ18O-Mg/Ca records provide questionable temperature and thus sea-level estimates as noted in (33). For this interval, sea-level fluctuations greater than 70 m and ice volumes inferred greater than today give little confidence that we have properly constrained the thermal history, and thus, sea-level amplitudes older than ~48 Ma are suspect. The Paleocene to earliest Middle Eocene interval is also complicated by 405-year scale hyperthermal warmings and the PETM. Our curves do not account for these transient warmings, so the sea-level record reflects hyperthermals as artifacts with sea-level highstands above 75 m. Nevertheless, the pacing of sea-level change during this interval appears to have significant power at the 1- to 3-Ma period (Fig. 2), suggesting that we have captured the timing of astronomical sea-level modulations but can say little meaningful about the amplitudes.


Early Eocene and the ice-free conundrum

Peak warmth, peak sea levels, high CO2 [>800; Fig. 1; likely >1000 ppm (9, 54)], and mostly ice-free conditions (Fig. 1 and fig. S3) occurred in the hothouse Early Eocene (55.0 to 47.9 Ma). Although our Early Eocene δ18O-Mg/Ca variations are suspect and the NJ backstripped record influenced by changing baseline due to MDT at this time (7), there is ample geological evidence to strongly suggest ice-free conditions during the Early Eocene (1, 4, 10). Large (>25 m), rapid (<1 Ma) sea-level changes during the Early Eocene pose a conundrum, because the primary mechanism for explaining them is the growth and decay of continental ice sheets (4, 55). As in the Paleocene, the Early Eocene δ18O-Mg/Ca sea-level falls appear to be much higher in amplitude (50+ m) than the backstripped estimates of 15 to 30 m (Fig. 1 and fig. S3), and given the uncertainties in our temperature estimates before ca. 48 Ma, we acknowledge that we have not fully accounted for temperature variations. Nevertheless, major Ma-scale sea-level falls apparently occurred in both δ18O-Mg/Ca and backstripped estimates during the Early Eocene at 54.6, 52.8, 50.5, and 49 Ma, with early Middle Eocene falling at 47.5, 46.9, and 45.6 Ma (fig. S3). However, temperature estimates (33) are inversely proportional to δ18O-Mg/Ca sea-level variations (with high temperatures on the 2-Ma scale associated with low sea levels; Fig. 1), and we caution again that our sea-level amplitude estimates with this method may not be valid before 48 Ma.

We conclude that the Early Eocene was largely ice-free and that global sea-level falls of ~15 to 30 m were caused by ice-volume increases or by other mechanisms such as changes in groundwater storage. Although early studies (4, 56) assumed that changes in lakes and groundwater could only explain about 5 m of GMSL change, recent studies [e.g., (57)] have suggested that changes in groundwater storage (currently about the same volume as stored in ice sheets) can explain greenhouse sea-level changes of up to 10 to 40 m. We view this as an intriguing idea that can be tested by more reliable estimates of Early Eocene sea-level variations.

Middle Eocene to Late Eocene cool greenhouse

Cool greenhouse conditions of the Middle to Late Eocene were associated with sea-level changes of 15 to 40 m caused by growth and decay of small (25 to 35% of modern) ice sheets (Fig. 1 and fig. S4). The δ18O-Mg/Ca and backstripped estimates are similar despite the hiatuses in the latter (Fig. 1 and fig. S4), particularly before, during, and after MECO and during the late Late Eocene (ca. 35 Ma) deglaciation. Our result is supported by direct measurement of Mg/Ca and δ18O on backstripped Middle Eocene sequences at Bass River, NJ (58). Our δ18O-Mg/Ca–based estimates indicate that the Earth reached near ice-free conditions during MECO, which was bracketed by ~30 m GMSL falls, and during the late Late Eocene deglaciation (Fig. 1 and fig. S4). δ11B shows declining CO2 concentrations during the Middle to Late Eocene (59), supported by other proxies (Fig. 1D). During the Middle to Late Eocene, a drop in CO2 was associated with a ~1‰ δ18O increase (Fig. 1), ~3°C deep-water cooling (Fig. 1), and GMSL fall of ~25 m (inventory ~40% of modern), interrupted by warming associated with the MECO and the late Late Eocene events.

New observations from the Antarctic coast indicate that the cryospheric evolution of an Antarctic ice sheet began in the Middle Eocene or before (11), challenging previous assumptions. Marine-terminating glaciers existed at the Sabrina Coast adjacent to the Aurora subglacial basin by the Early to Middle Eocene, implying substantial ice volume before continental-scale ice sheets of the Oligocene (11). Age control on these strata is limited, and it cannot be conclusively shown that there were Early Eocene ice sheets, but certainly, marine terminating ice sheets existed by the Middle Eocene, indicating that Antarctica was partly glaciated intermittently during this cool greenhouse (11). It is also intriguing that ice-related debris appeared in the Arctic in the Middle Eocene (15).

Oligocene to Early Miocene large ice-volume variations

During the largely unipolar icehouse of the Oligocene to Early Miocene, the EAIS was not permanently developed, with periods of large-scale (~50 to 60 m sea-level equivalent) growth and decay. The O1 δ18O increase that heralded the onset of continental-scale glaciation (12, 60) was associated with a backstripped sea-level fall of 52 ± 10 m, although the lowstands were not fully sampled and this must be considered a minimum fall. Our δ18O-Mg/Ca estimate shows two EOT steps (EOT-1 and Oi1) as indicated by previous studies (61): ~25 m at EOT-1 (~33.9 Ma) and 50 m at Oi1 (33.65 Ma). Backstripping indicates a large deglaciation at ca. 32 Ma (Oi1a) with a 44 ± 10 m rise, suggesting loss of two-thirds of ice that was gained across the EOT glaciation; in contrast, δ18O-Mg/Ca estimates suggest approximately a 35-m rise (~50% deglaciation). Both proxies suggest that the EAIS was relatively unstable throughout the Oligocene, being paced by 1.2-Ma tilt cycle (Fig. 2) (13). Both Oi2 and its attendant subdivisions and Mi1 (ca. 23 Ma) had sea-level lowerings only slightly larger than at Oi1 (−30 to −40 m), suggesting an upper limit on ice storage on Antarctica. Similarly, the maximum sea levels attained (+30 to 40 m) are consistent but well below the ice-free line, suggesting the presence of a moderate to large (25 to 35 m) ice sheet in Antarctica. This pattern continued into the Early Miocene with 40 to 60 m variations on the long tilt (1.2 Ma) scale (13), when it was punctuated by the MCO beginning at ca. 17 Ma.

Miocene climate optimum

From an ice-volume perspective, the MCO was a period of reduced ice volume from 17.0 to 13.8 Ma, although from a δ18O perspective the MCO ended at 14.8 Ma (42). Near ice-free conditions were attained during peak 405 ka cycles (Fig. 1), which may have been the most recent ice-free period on Earth. The MCO was punctuated by the Mi2 (ca. 16.0 Ma) and Mi3a (14.8 Ma) δ18O increases and attendant sea-level falls of ~40 and ~30 m, ending with the Mi3 (13.8 Ma) δ18O increase (~50 m fall). Other than the Ma-scale Mi2 and Mi3a events, sea-level changes were small (<20 m) during the MCO.

MMCT and Late Miocene

From δ18O perspective (42), the MMCT consists of three major Ma-scale δ18Obenthic increases, coolings, and attendant sea-level falls: Mi3a (14.8 Ma; ~30-m fall and ~0.7°C cooling), Mi3 (13.8 Ma; ~50-m fall and ~1.2°C cooling), and Mi4 (12.8 Ma; 20- to 30-m sea-level fall and ~1.0°C cooling), culminating in the “permanent EAIS” (62) with GMSL less than 12 m above present (Fig. 1). The Mi3a and Mi3 δ18Obenthic increases reflect more ice growth than temperature (i.e., two-thirds of the δ18Obenthic change was due to ice), and the Mi4 increase had similar contributions by ice and temperature. Overall, the MMCT was a net cooling of ~3°C. Although the net increase in ice was 60-m sea-level equivalent across the MMCT, the inventory after the MCO was not much larger than before; however, after the MCO-MMCT, the subsequent EAIS was relatively stable, with smaller ~20- to 30-m Ma-scale sea-level variations. Thus, the ice sheet in the Oligocene to early Middle Miocene was very dynamic [not unexpected considering the warm deep waters (~6° to 8°C) and thus warm Antarctic surface waters] and much less dynamic after the MMCT. After the MMCT, deep-sea temperature remained low (<5°C), and the inferred cooling of high-latitude surface waters was associated with increased stability of the EAIS and its poor response to orbital forcing (37).

Antarctic continental records support the Miocene ice sheet history inferred from our sea-level records. Miocene warmth before 14 to 15 Ma is indicated in East Antarctica by temperate alpine glaciations, colluvium and soil development, and plant macrofossils (51, 63). There was a major change from wet-based glaciation 14.07 ± 0.05 Ma to the onset of cold-based glaciations at 13.85 ± 0.03 Ma (51), along with a 6 to 7°C cooling of sea surface temperatures near Antarctica between 14.2 and 13.8 Ma (64). By ca. 12.8 Ma, the development of permanent cold-based glaciers is indicated by limited ice expansion, absence of plant macrofossils, little soil development, and long-term ice preservation (65).

Sea level in the late Middle to Late Miocene was unexpectedly constant, with one ~30-m Ma scale fall at ca. 8.2 Ma and sea level seldom exceeding 20 m above present. During the Late Miocene, sea level was paced by the 41-ka tilt cycle (Fig. 2) (37), and there was little sea-level change during the Messinian salinity crisis initiated at 5.96 ± 02 Ma, nor with the isolation of the Mediterranean Sea from 5.59 to 5.33 Ma (although there is at a splice in the record at 5.32 Ma). The amplitude of Milankovitch-driven sea-level variations increased during the Early Pliocene (ca. 5.33 to 4 Ma), with progressively higher peak sea levels.

Pliocene climate optimum

Sea level rose in the Early Pliocene ca. 4.7 Ma, reaching highs that had not been consistently seen since the MCO. From a sea-level perspective, the Pliocene is marked by three intervals with sea level well (~10 to 20 m) above modern: 4.6 to 4.1, 3.9 to 3.3, and 3.3 to 2.85 Ma. The latter interval has been studied in detail by the PRISM project (66) and sometimes referred to as the PCO or mid-PCO (67), during which sea level peaked in Chron G17 (ca. 3 Ma) at ~20 m (68, 69). There are large uncertainties in the peak sea level at ca. 3 Ma (67, 68), and it is difficult to impossible to place firmer amplitude constraints [errors estimated as ±10 m (68) may be larger (67); see Materials and Methods]. However, benthic foraminiferal δ18O values suggest GMSL peaks of 20 m or less (68). GMSL higher than 12 m above modern requires loss of ice sheets in Greenland, West Antarctica, and sensitive areas of East Antarctica, the Wilkes, and Aurora Basins (70). This interval is of keen interest, because global temperatures were >2°C warmer than today (71) at times when atmospheric CO2 concentrations were on the order of those in 2020 CE (~400 ppm), and thus, the equilibrium sea-level state is relevant to ice sheet trajectories in coming centuries. The peaks between 3.9 and 3.3 Ma were even higher, reaching a peak of ~30 m during Gi13 (ca. 3.8 Ma; Fig. 1), and thus requiring some melting of the EAIS.

Quaternary ice ages

Large sea-level changes with lowerings up to 120 to 130 m below present were restricted to the past 2.7 Ma (Fig. 3), indicating development of full-scale continental-scale NH ice sheets during the Quaternary. A large (>1‰) δ18O increase associated with MIS100 (2.54 Ma) was first-order correlated with the apparent first appearance of dominant ice-rafted debris in the northern North Atlantic (14). Sometimes dubbed “the initiation of NH glaciation,” this change should be more accurately termed the “beginning of the NH ice ages” because of ample evidence for NH glaciation (GIS or larger scale) before the Pliocene (72). The interval from MIS G10 (ca. 2.8 Ma) to MIS100 represents a progressive increase in glacial-interglacial δ18O amplitudes, suggesting that this was a continuous process rather than an abrupt event. Our new sea-level estimates suggest that this event was marked by progressively larger glaciations with little trend in deglaciations. From 2.5 Ma to roughly 1 Ma, 41-ka tilt forcing continued to pace sea-level changes (Fig. 2), as shown in previous studies.

Fig. 3 Blowup of Quaternary GMSL and CO2 estimates.

Each site has a contrasting shade of blue, with the site number indicated on the left column in the corresponding color. Sea-level estimates of Waelbroeck et al. (94) (green) based on deep-sea δ18O and Rohling et al. (95) based on scaling of Red Sea data (magenta) are shown. The black line is the prediction of CO2 through time based on the data compiled by Foster et al. (47), and the dashed and dotted lines demarcate the 1σ and 2σ confidence intervals for that prediction, respectively. The Laurentide line is drawn at −50 m.

The largest sea-level changes (up to 130 m below present) of the Cenozoic were associated with quasi–100-ka terminations of the last 800 ka (Figs. 3 and 4). The major precessional (19 and 23 ka) and tilt (41 ka) scale sea-level lowerings of the past 800 ka were on the order of 10 to 60 m (Figs. 3 and 4) and were directly forced by Milankovitch pacing (Fig. 2). During the mid-Pleistocene transition (ca.1 to 0.8 Ma), large ice sheets began to display stronger quasi–100-ka cyclicity (Fig. 2). Note that our analysis shows that 95/125-ka power occurs in δ18O records before 1 Ma (e.g., last 3 Ma and older than 14 Ma; Fig. 2). Nevertheless, sea-level records of the last Ma first show the typical, dominant quasi–100-ka “saw tooth” cycles noted in the δ18O record (Figs. 2 to 4). The 405-ka cycle was prominent in the Miocene δ18O and sea-level records until ~12.8 Ma but is poorly developed to absent in the Quaternary (Fig. 2). The cause of the shift to dominant quasi–100-ka periods during the mid-Pleistocene transition and the absence of the 405-ka cycle in Quaternary δ18O records are unknown, although the 100-ka beat is likely due to amplification by CO2 (73). Long-term decreasing atmospheric CO2 during the Quaternary (from ~320 to 250 ppm) may have reached a threshold, resulting in the return of a 100-ka beat that had previously been important in the early history of the EAIS (ca. 34 to 12.8 Ma; Fig. 2).

Fig. 4 Late Pleistocene GMSL estimates.

Sea-level estimates of Fairbanks (99) from Barbados (green), Waelbroeck et al. (94) using deep-sea δ18Obenthic (light blue), and Rohling et al. (95) based on scaling of Red Sea δ18Oplanktonic data (red) are shown. Our sea-level estimates use a splice of site 845 (dark blue) and V19-30 (black).

On the basis of our observed responses of δ18O and sea level to forcing (Fig. 2), we note that large ice sheets (e.g., EAIS in the Oligocene-Middle Miocene; the Laurentide Ice Sheet (LIS) during the last 2.55 Ma) initially responded directly to Milankovitch changes due to precessional and tilt forcing as previously noted for the Quaternary (74). We suggest that these ice sheets were kicked into new states by the size of the large ice sheets themselves (Figs. 1 and 2) and, because of their large sizes, responded to different forcings. For example, the development of a permanent EAIS by 12.8 Ma resulted in a change from responding to precession, tilt, and eccentricity to subdued to absent response to eccentricity and precessional forcing that had previously been strong [see also (38)]; the 41-ka tilt cycle dominated ice sheet and sea-level response from ca. 12.8 to 1.0 Ma following the development of a permanent EAIS. During the mid-Pleistocene transition, very large, 100-ka paced LIS were amplified by 100-ka changes in CO2 (73) from ~180 (glacial) to 280 ppm (terminations).

The LGM (MIS2; ca. 27 to 20 ka) stands as the largest lowering of GMSL (~130 m) of the Mesozoic-Cenozoic (Fig. 4), although GMSL was higher than modern during MIS5e, MIS9, MIS11, and perhaps MIS 31. During the last deglaciation (ca. 19 to 10 ka), GMSL rise exceeded 40 to 45 mm/year (75), providing an upper limit on known rates of GMSL rise. Rates before radiocarbon ages are less certain, although the sea-level rises exceeded 10 mm/year during terminations. Sea-level rise progressively slowed during the Holocene until the late 19th to early 20th century when rates began to rise from near 0 to 1.2 mm/year in the early 20th century to a late 20th and 21st century rise of 3.1 ± 0.4 mm/year (2, 3).

CO2, climate, and sea level

Sea level follows long-term trends of atmospheric CO2, with high sea levels associated with high CO2 and warm climates. CO2 played an important role with high CO2 maintaining warmth in the Eocene (with values >800 to 1000 ppm; Fig. 1) associated with largely ice-free conditions and high sea levels. Generally, decreasing CO2 values during Middle Eocene to Oligocene led to cooling and glaciation (Fig. 1), although a secondary CO2 increase at ca. 35 to 36 Ma may be associated with the late Late Eocene warming. The cause of the CO2 decrease over the past 50 Ma has been widely discussed and debated but must be due to long-term (107-year) changes in CO2 sources (e.g., higher CO2 associated with inferred higher ocean crust production rates) or more likely the effectiveness of sinks CO2 (e.g., increased weathering associated with uplift of the Himalayas or exposure of basalts in tropical regions) (76, 77).

We note that sea-level variations mostly follow Milankovitch astronomical forcing, responding to orbital eccentricity, tilt (including the 1.2 Ma cycle), and precession as the dominant pacing (Fig. 2). Relationships among proxy records of CO2, climate, and sea level are not always clear (Fig. 1). For example, large changes in Miocene cryospheric evolution and deep-water temperature (Fig. 1; e.g., the MCO and MMCT) occurred with only what appear to be modest changes in CO2. The MCO is an interesting case in point where B isotopes indicate variations from 300 to 500 ppm during this interval of ~4°C warmer global temperatures above present (78). Our records that suggest nearly ice-free conditions occurred during the MCO and are thus intriguing if this is an equilibrium state for warming levels that will be attained in this century or the next century under sustained greenhouse gas emissions. Similarly, the long-known mismatch between Eocene to Oligocene CO2 proxies (with the main CO2 drop in the later Oligocene) and the earliest Oligocene glaciation has been perplexing, although we show that decreases in CO2 during the Middle to Late Eocene are associated with bottom water cooling and ice growth. In all cases, we note that errors on pre-ice core CO2 values are large ( hundreds of parts per million) and could mask closer relationships. We look forward to more refined estimates of Eocene to Miocene CO2 to test these relationships in detail.

Amplitudes of sea-level change

The estimates provided here from scaling of deep-sea δ18Obenthic and backstripping have large errors (±10 to ±20 m; see Materials and Methods). Given this, the agreement shown between the two estimates (Fig. 1D) is remarkable (see Materials and Methods for statistical comparisons of the two records). Despite the uncertainties, both estimates show typical Ma-scale variations of 20 to 50 m driven by ice-volume variations. Similar Ma-scale amplitudes have been found from backstripping the Australian margin for the Miocene (17) and the New Zealand margin in the Pliocene (79). Higher-frequency Milankovitch-scale sea-level changes vary from <10 to 20 m before the Quaternary up to 120 m in the last glacial cycle.

Convergence of Eocene to Miocene estimates from different methods (Fig. 1D) and from far-flung parts of the world calls into questions the very high (up to 160 m) estimates of the EPR sea-level curves (8082). A comprehensive review of Cretaceous sea-level records similarly casts doubt on the high EPR estimates (83). The EPR estimates are not reproducible because the data are not published, and methods used to produce their sea-level estimates do not rigorously account for subsidence and other processes. Their estimates were obtained by scaling major flooding surfaces observed in the stratigraphic records to long-term (107-year) sea-level estimates derived from sea-floor spreading (~250 m above present at 80 Ma ago), with lowstands chosen to not exceed Pleistocene sea-level falling except during the middle Oligocene (24). The relative amplitudes of events in the EPR curves also appear to be incorrect. For example, the largest EPR falls are inferred for the middle Oligocene and late Middle Miocene, yet the largest GMSL falls are at the EOT (~55 m) (4) and during peak Quaternary glacials when sea level was ~120 to 130 m lower (e.g., the best constrained lowstand at Barbados of 120 m below present yields a GIA-corrected GMSL of 127 m; Fig. 4) (84). We again note that the timing of the EPR curves still appears to be useful and that these “cycle charts” should be viewed as useful in their original intention as a chronology of sea-level changes, although they are largely devoid of amplitude information.

Conclusions and future work

We present a Cenozoic sea-level estimate derived using a new deep-sea Pacific benthic foraminiferal δ18O splice and Mg/Ca records to constrain long-term temperature, providing a record of ice-volume changes. Although our δ18O-Mg/Ca–based GMSL estimates have large errors (±10 to ±20 m; see Materials and Methods) and are probably not correct older than 48 Ma due to uncertainties in the Mg/Ca of seawater, major sea-level falls from the Middle Eocene to Miocene show similar amplitude and timing as those observed on passive continent margins, indicating a global cause, GMSL. Amplitudes of GMSL are 10 to 60 m on the Ma scale, agreeing with estimates obtained from backstripping middle Atlantic U.S. passive margin records. Our sea-level history constrains cryospheric evolution over the past 66 Ma, with ice-free conditions during most of the Early Eocene, MECO, latest Eocene, and possibly the MCO, with ice sheets (up to 40-m sea-level equivalent) in the Middle to Late Eocene greenhouse and with continental-scale Antarctic ice sheets beginning in the Oligocene. From 34 to 13.8 Ma, the EAIS varied from larger than today (sea-level ~35 m below present) to nearly ice-free (~50 m above present) but became permanent during the MMCT ca. 12.8 Ma. The PCO saw loss of ice sheets in Greenland, West Antarctica, and likely sensitive regions of East Antarctica. Quaternary sea-level variations were the largest of the Cenozoic, as large NH ice sheets grew and decayed.


Time scale and δ18Obenthic records

We use the geological time scale (GTS) 2012 (85) and its recent updates (86) throughout. We compiled published δ18Obenthic 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 δ18O records (Fig. 1 and fig. S2). With publication of recent δ18Obenthic 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 δ18Obenthic 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, 2325).

The efforts of the IODP paleoceanographic community are evident in the resolution of the δ18Obenthic 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 (δ18OCibicidoides) using published (89) calibrations. V19-30 and site 865 data were obtained from Uvigerina (Cibicidoides = δ18OUvigerina – 0.64‰). Data from the middle Middle Eocene at site 1209 are from both Nuttallides truempyi and Oridorsalis umbonatus (88) that were corrected [Cibicidoides = (δ18ONutallides + 0.1)/0.89‰ =δ18OOridorsalis – 0.28‰]. All other records were run using Cibicidoides-Planulina (mostly Cibicidoides mundulus/praemundulus or Planulina wuellerstorfi). The δ18Obenthic data used in Fig. 1A, temperature [equation 7A in (33)], the δ18Obenthic-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 δ18Obenthic 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)T=16.14.76 [δ18Obenthic(δ18Oseawater0.27)](2)to solve for δ18Oseawater, yielding a δ18Obenthic-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 δ18Oseawater 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 δ18Oseawater estimate was scaled to GMSL changes (Fig. 1) using a revised δ18Oseawater 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 δ18Oice, and the calibration of 0.13‰/10 m requires a δ18Oice of −50‰. The δ18Oice 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 (Δδ18Oseawater, 1.2‰) and 16 m on typical changes (Δδ18Oseawater, 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), δ18Obenthic values corrected to Cibicidoides, δ18Oseawater, 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 δ18Obenthic 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/Caseawater, assuming a high Mg/Cabenthic sensitivity, or inaccuracies in the Mg/Cabenthic 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 δ18O-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 δ18O-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 δ18O-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 δ18O-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 δ18O-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 δ18O-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 δ18O-Mg/Ca–based and backstripped sea-level estimates, ai (Eq. 3) and bj (Eq. 4), respectively. We modeled the recorded sea-level variations from the δ18O-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]ai=g(ti)+εia(3)bj=h(tj)+εjb(4)g(t)=m(t)+l(t)+wg(t)+cl(5)h(t)=m(t)+Δ(t)+l(t)+lh(t)+wh(t)+cl+clh(6)

The observational errors (εiaand εjb) were treated as having SDs of ±10 m for the δ18O-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 lh(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), lh(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 [wg(t) and wh(t)] capture high-frequency variability and measurement error beyond that accounted for in the nominal error estimates. The constant terms (c1 and c2) 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), lh(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, lh(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 δ18O-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 δ18O data (Fig. 2) and the δ18O-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 δ18O data and δ18O-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 δ18O and δ18O-Mg/Ca GMSL spectrograms.

To evaluate the shortest periods resolvable for a given time, we plotted the “local” Nyquist frequency of the δ18O data and δ18O-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 δ18O 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 δ18O record and δ18O-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.

CO2 records

We constructed a Bayesian hierarchical model to interpolate the atmospheric concentration of CO2 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 CO2 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 CO2 dataset (47). The statistical model allows us to generate the distribution of the atmospheric concentration of CO2 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 CO2 at measured points in time, as well as the uncertainty associated with both the CO2 concentration and temporal measurements; (ii) a process level that represents the distribution of CO2 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 CO2 data comprises the data level of our hierarchical model. The CO2 values, yi (Eq. 7), and the measurement ages, ti (Eq. 8), are modeled as the sum of their representative values (t¯i is the mean measurement age) and their errors, εiy and εit, respectivelyyi=f(ti)+εiy(7)ti=t¯i+εit(8)

At the process level, the concentration of CO2 through time f(t) is represented as a function of five component processes (Eq. 9).f(t)=P1(t)+P2(t)+C1(t)+C2(t)+C3(t)(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 CO2 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-CO2 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 CO2 concentrations, such as variations in outgassing or paleogeographic configuration, if present. The hyperparameter distributions that best represented the CO2 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 CO2 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 CO2 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 CO2 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.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank the IODP paleoceanographic community for their efforts on generating the exquisite records (most notably A. Holbourn and T. Westerhold), R. Fairbanks for his original (1984) suggestion to “scale δ18O records to sea level,” numerous colleagues who participated in generation of the “NJ sea-level curve” (especially M. Kominz, P. Sugarman, and S. Pekar), B. Cramer for his Mg/Ca synthesis, and the founders of pre-Quaternary paleoceanography and cryospheric communities (especially J. Kennett, N. Shackleton, S. Savin, C. Denton, and R. K. Matthews) who saw much of the broad picture, although we may now differ on critical details. We thank M. Simmons (Halliburton) for comments and three anonymous reviewers. We thank the Office of Advanced Research Computing (OARC) at Rutgers University for providing access to the Amarel cluster and associated research computing resources. Funding: This study was supported by NSF grants OCE16-57013 (K.G.M.) and OCE14-63759 (K.G.M. and J.V.B.). Author contributions: K.G.M. oversaw the studies and wrote the manuscript. J.V.B. contributed to assembly and scaling of stable isotopic data. W.J.S. conducted all time series analyses and computed statistical comparisons of records. R.E.K. contributed to discussions of sea-level change and statistical evaluations. G.S.M. contributed to the NJ sea-level records. J.D.W. contributed to discussion and understanding of stable isotopes and cryospheric history. 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 present in the paper and/or the Supplementary Materials. δ18O and sea-level curves posted on NGDC site Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article