Global Mean and Relative Sea-Level Changes Over the Past 66 Myr: Implications for Early Eocene Ice Sheets

Highlights

• We present a new global mean geocentric sea level record integrating updated barystatic (ice volume), thermosteric(temperature), and long-term ocean basin volume estimates for the Cenozoic (past 66 Myr).

• We show two independent sea-level estimates (barystatic and backstripped relative) that have large (>25 m) Eocene sea-level falls, requiring significant ice sheets.

• Although the Early Eocene was ice-free at times, our two independent sea-level estimates show lowerings that refute current views of a completely ice-free Early Eocene.

Introduction

Earth’s climate and sea level have changed significantly from the end-Cretaceous mass extinction 66 million years ago (Ma) to the Anthropocene, with transitions from ice-free conditions during portions of the Eocene to bipolar ice sheets in the Quaternary. The Maastrichtian Age that ended the Cretaceous was relatively cool and may have had modest ice sheets (Miller et al., 2005) compared to the peak Cretaceous warmth (Huber et al., 2018). Following decadal to 100 kyr scale global temperature perturbations associated with asteroid impact (e.g., Vellekoop et al., 2019), Late Paleocene to Early Eocene (ca. 60–54 Ma) global climates warmed again and sea levels rose. Deep waters, reflecting high latitude surface temperatures, warmed to over 13°C compared to less than 2°C today (Savin et al., 1975; Kennett and Shackleton, 1976; Miller et al., 1987; Zachos et al., 2001). Peak Cenozoic warmth and sea level occurred in the Early Eocene (ca 54–48 Ma) with largely ice-free conditions (“Greenhouse World”; Miller et al., 1991, Miller et al., 2020; Zachos et al., 2001; “Hothouse World,” Westerhold et al., 2020). A secular cooling followed in the Middle to Late Eocene (ca. 48–34 Ma) with an increase in polar ice and general sea-level fall (Miller et al., 2020) during the Cool Greenhouse or Warmhouse World (Zachos et al., 2001; Miller et al., 1991; Miller et al., 2020; Westerhold et al., 2020). Our previous studies have shown that the generally warm Middle to Late Eocene interval has a record of modest-sized, ephemeral continental ice sheets (Figure 1; e.g., Browning et al., 1996; Miller et al., 2005, Miller et al., 2020). There was an abrupt shift in the earliest Oligocene (ca. 34 Ma) to a continental-scale Antarctica Ice Sheet (AIS) with peak ice volumes likely greater than today (Miller et al., 1987; Zachos et al., 2001; Miller et al., 2020). Sea-level and ice volume varied during the Oligocene to Early Miocene, with major deglaciation during the Miocene Climate Optimum (MCO; 17.0–14.9 Ma; Miller et al., 2020). A permanent East Antarctic ice sheet developed in the Middle Miocene (ca. 13.8 Ma). The Icehouse World of the past 34 Myr (Miller et al., 1991) can be divided into two states (The Coolhouse and Icehouse of Westerhold et al. (2020)), the former with large, but varying East Antarctic Ice Sheet (EAIS) in the Oligocene to Middle Miocene and the latter with a largely stable EAIS during the last 13.8 Ma. Northern hemisphere ice sheets existed back at least to the Middle Miocene (e.g., Wolf-Welling et al., 1996) and perhaps to the Middle Eocene (Stickley et al., 2009). Growth and decay of large Northern Hemisphere (Laurentide) ice sheets began in earnest during Quaternary Marine Isotope Stage 100 (ca. 2.6 Ma; aka “the ice ages,” summary in Miller et al., 2020).

FIGURE 1. Panel (A) Eocene to Oligocene Pacific benthic foraminiferal (Cibicidoides spp.) δ¹⁸O data (gray line) (Miller et al., 2020; this study) with an 800 kyr LOESS smoothed fit (solid black line); Mi1, Oi2a, Oi2, Oi1 are benthic foraminiferal δ¹⁸O maxima, LECO, Late Eocene Climate Optimum; MECO, Middle Eocene Climate Optimum; EEOC, Early Eocene Climate Optimum; Paleo., Paleocene; Mio., Miocene. Panel (B): timing of glacial events in Antarctica (compiled by Miller et al. (2008) and updated with Gulick et al., 2017; Barr et al., 2022); Panel (C): compilation of statistically derived estimates of Cenozoic pCO₂ (this study) from alkenone (yellow dots), boron (purple dots), paleosols (green dots), stomatal indices (blue dots). The best estimate of Cenozoic pCO₂ is represented by function shown with a heavy red line that runs through the boron isotope and is systematically offset from the other proxy measurements. The proxy-specific offsets were determined using a Gaussian process regression model and they are shown in Supplementary Figure S3.

Although this narrative of Cenozoic ice volume, climate, and sea-level history is largely accepted, several thorny problems and controversies remain. The record of ice sheets prior to the Oligocene is hotly disputed. The earliest Oligocene δ¹⁸O increase (Oi1; ca. 34 Ma; Miller et al., 1991) is often incorrectly referred to as the “initiation of Antarctic glaciation.” Instead, this was the first time since the Paleozoic that ice sheets covered most of the Antarctic continent. Browning et al. (1996) suggested that continental ice sheets developed in the Middle to Late Eocene based on continental margin and epicontinental sea RSL changes associated with major δ¹⁸O increases in deep sea benthic foraminifera. Miller et al. (2020) used deep Pacific benthic foraminiferal δ¹⁸O and Mg/Ca records to calculate seawater δ¹⁸O, hence sea-level changes, due to ice-volume variations. Using this method, they estimated over 25 m sea-level changes during parts of the Middle to Late Eocene and validated these changes with comparison to backstripped RSL changes; these independent estimates require significant ice sheets. The reconstructions of Miller et al. (2020) were deemed valid to the Middle Eocene (ca. 48 Ma), and the record older than 48 Ma was presented as less certain.

Most paleoclimate studies have assumed an ice-free Late Cretaceous to Eocene world (e.g., Zachos et al., 2001), though seismic and sequence stratigraphy studies suggested Myr scale changes of over 100 m (Vail et al., 1977; Haq et al., 1987). Such high amplitudes and implied rates (100 m/Myr; 10 cm/kyr) are difficult to explain in warm climate intervals such as the Early Eocene that are inferred to be ice-free since only ice-volume variations can explain such high rates of change (Pitman and Golovchenko, 1983; Miller et al., 2005; Miller et al., 2020). Here we address the question of whether there were ice sheets before 48 Ma (Early Eocene and older) and provide new insights and sensitivity studies on the size of ice sheets throughout the Cenozoic.

Sedimentological evidence (tills, ice-rafted detritus [IRD]) from Antarctica requires significant Middle to Late Eocene glaciation but records are ambiguous about Early Eocene glaciation (Figure 1). Early sedimentological evidence from on or near Antarctica pointed to the presence of at least mountain glaciers in Antarctica during the Eocene (Figure 1) (Margolis and Kennett, 1971; Birkenmajer, 1988; Barron et al., 1991; Birkenmajer et al., 2005). Recent mapping of Antarctica cirques suggest Paleocene to Eocene mountain glaciation (Figure 1) (Barr et al., 2022). Firm evidence of continental ice has been documented by grounded Middle Eocene tills in the Aurora Basin, offshore Antarctica that possibly extend into the Lower Eocene (Gulick et al., 2017), documenting significant continental glaciation (Figure 1). However, Early Eocene age estimates on the Antarctic tills (Gulick et al., 2017) and identification of Lower Eocene ice-rafted detritus in the Southern Ocean cores (e.g., Margolis and Kennett, 1971) are equivocal (Figure 1). The deep-sea records used as a proxy for ice volume and sea level have also been equivocal because estimates obtained from δ¹⁸O and Mg/Ca records have been poorly constrained due to uncertainties in the Mg/Ca record for the Early Eocene and older (Cramer et al., 2011; Miller et al., 2020). The null hypothesis that Antarctica was ice-free throughout the Early Eocene has been assumed to be true, and this perspective is justified by evidence for Early Eocene warmth (e.g., Inglis et al., 2020). We show here with new records of Early Eocene Mg/Ca data, that this hypothesis is incorrect.

Specifically, we address these questions by providing new Late Paleocene-Early Eocene (ca. 48–58 Ma) Mg/Ca data from Pacific Site 577, revising the existing Cenozoic Mg/Ca compilation of deep-water temperature by incorporating only Pacific and Southern Ocean data in the new Mg/Ca synthesis and providing a revised record of BSL changes due to changes in ice volume. We combine our BSL record with a new record of sea-level changes due to ocean volume (Wright et al., 2020) and thermosteric changes and provide a testable record of GMGSL. Our BSL and GMGSL estimates are similar to “backstripped” RSL estimates (accounting for compaction, loading, and thermal subsidence) from the US Mid-Atlantic Coastal Plain, and the BSL record we produce mirrors that of those RSL estimates accounting for the ocean basin volume and thermosteric effects, within errors of ±16 m since 42 Ma and ±27 m through the Cenozoic. Ours is not only the first BSL estimate for the Early Eocene, but it also provides a testable GMGSL record for the entire Cenozoic. We also compare with other published sea-level estimates (Haq et al., 1987; Snedden and Liu, 2010; Rohling et al., 2021; Rohling et al., 2022).

Background: Global, Relative, and Barystatic Sea Level

GMSL change is the global average variation in sea level caused by variation in ocean volume distributed across the area of the world oceans that is driven by variations in density that correspond with warmer or cooler global ocean temperatures (thermosteric sea level; TSL) and in the mass of water in the ocean due to changes in ice volume on land and in terrestrial water storage (BSL; Gregory et al., 2019). Today GMSL is rising at 3.3 mm/yr due to ocean warming and melting ice (glaciers and ice sheets; Nerem et al., 2010;¹) accelerating from a 20th century rise of ∼1.4 mm/yr (Hay et al., 2015). On long time scales (>1 Myr), changes in ocean basin volume also cause sea-level variations (OBVSL), which together with BSL and TSL, combine to cause global mean geocentric sea level (GMGSL) changes with respect to the terrestrial reference frame (Gregory et al., 2019).

Geologists have often used the term “eustasy” for global sea-level variation relative to a fixed datum, usually the center of the Earth (Posamentier et al., 1988). We avoid the term eustasy and adopt the more specific term GMGSL to refer to geocentric sea-level variation. We note that the temporal evolution of a geocentric reference frame cannot be readily constrained over geological timescales and avoid the terms eustasy and glacio-eustasy in favor of GMGSL and BSL. While the formal definition of GMGSL invokes a level relative to a geocentric reference frame, we treat changes in this height due to variations in the shape of the ocean basins as the effect attributable to the modification of both the volume of the ocean basins and the surface area in which the volume of ocean water resides. Therefore, we consider GMGSL change to reflect a value defined as the sum of volumetric changes to both the oceans and the ocean basins divided by the temporally varying surface area of the oceans. This conceptually supports the combination of estimates of OBVSL and GMSL change to derive an estimate of GMGSL change and separates the GMGSL metric from the height of the ocean surface relative to a position/datum that is unknown in the geologic past.

RSL is the difference in height between the sea surface and the solid Earth (Gregory et al., 2019), with solid Earth often referring to the Earth’s crust (Posamentier et al., 1988). RSL change at a particular place is controlled both by GMGSL 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), mantle dynamic topography (MDT), changes in the height of the geoid, and short-term (10–1000 yr scale) ocean dynamic sea-level changes. Changes in RSL and sediment supply produce transgressions (landward movement of the shore) and regressions (seaward movement of the shoreline due to RSL fall and/or increased sediment input). Transgressions, regressions, and unconformities are traditionally held to be the expected stratigraphic expression of sea-level change (Vail et al., 1977), though the phase relationship between sea-level and stratigraphic response is ambiguous considering other variables that affect sedimentation and accommodation (Neal and Abreu, 2009; Miller et al., 2018). Challenges in extracting a GMGSL signal from the stratigraphic record include uncertainties in sea-level proxies, chronologic limitations, and overprinting by regional and local effects. Nevertheless, the Cenozoic sea-level changes from seismic and sequence stratigraphy suggest large Myr scale changes (Vail et al., 1977; Haq et al., 1987) that are difficult to explain in warm, supposedly ice-free climate intervals such as the Early Eocene.

Long-term (>2–10 Myr) global sea-level changes have high amplitudes (>100 m) but slow rates of change (<10 m/Myr) reflecting the dominance of slow tectonic processes, especially changes in ocean crust production driven by spread rate and ridge length changes (Vail et al., 1977; Pitman and Golovchenko, 1983). Published estimates of long-term sea level have wide ranges (e.g., ∼100–350 m for the Late Cretaceous [ca. 80 Ma] peak; Miller et al., 2008), though estimates converge at low values for the past 40–50 Myr. In presenting a GMGSL estimate the error bars both within and between different proxies are quite large in the older part of our record (>50 Ma).

Previous studies have demonstrated that large continental ice sheets and attendant sea-level changes occurred not only in the Oligocene (e.g., Miller et al., 1991; Zachos et al., 1992) but also in the Middle to Late Eocene (Browning et al., 1996; Miller et al., 2005; Miller et al., 2020). However, the question of what drove sea-level variations during the purported ice-free Cretaceous to Early Eocene remains. We focus here on Eocene sea levels and climates, with new insights into the Early Eocene.

Materials and Methods

Time Scale and Benthic Foraminiferal δ¹⁸O_b Compilation

We use the Geological Time Scale (GTS) 2020. We used our previously published (Miller et al., 2020) compilation of deep sea benthic foraminiferal oxygen isotopes (δ¹⁸O_b; largely Cibicidoides spp. and thus ∼0.64‰ lower values than equilibrium) records from 0 to 66 Ma (where the Cretaceous/Paleogene boundary on the GTS is at 66 Ma). Our record splices kyr-sampled deep (>2 km) δ¹⁸O_b records exclusively from Pacific sites (Figures 1–3). All δ¹⁸O_b records have been astronomically tuned to the GTS by the publications cited in the synthesis. Our δ¹⁸O_b record is similar to δ¹⁸O_b splices of De Vleeschouwer et al. (2017) and Westerhold et al. (2020) who used a similar approach to assemble a synthesis of δ¹⁸O_b records to evaluate the response to Milankovitch forcing and to define climate states, although they used data from both the Pacific and Atlantic Oceans.

FIGURE 2. Global bathymetric map (Eckert IV Projection of data from the GEBCO Compilation Group, 2023 (doi:10.5285/f98b053b-0cbc-6c23-e053-6c86abc0af7b); showing sites used in the benthic foraminiferal δ¹⁸O splice (red), Mg/Ca compilation (blue), or both (red and blue). Scale is relative to Mean Sea Level. The map shows the elevation of the ice surface in Antarctica and Greenland.

FIGURE 3. Panel (A) ${\mathrm{\delta }}^{18}{O}_{Cibicidoides}$ splice (gray line) of Miller et al. (2020) with an 800 kyr LOESS smoothed fit (solid black line). Panel (B) Deep-sea Mg/Ca paleotemperature estimates (sites in legend), with a 1.5 Myr LOESS regression (red line) showing the long-term trend. The red line shows paleotemperature if 20% of the variability in ${\mathrm{\delta }}^{18}{\mathrm{O}}_{\text{sw}}$ calculated is due to temperature variation The black line is a LOESS regression of the ${pT}_{R_1}\left(t\right)$ record (red line) with an 800 kyr search window. Panel (C) The ${BSL}_{R_1}\left(t\right)$ record (light blue line) is sea level computed using a 0.13/10 m ${\mathrm{\delta }}^{18}{\mathrm{O}}_{\text{sw}}$ to sea level calibration, assuming that 80% of the variability in ${\mathrm{\delta }}^{18}{\mathrm{O}}_{\text{sw}}$ calculated using the paleotemperature equation (Eq. 1), ${\mathrm{\delta }}^{18}{O}_b$, and ${pT}_{R_0}\left(\mathrm{t}\right)$ is due to seawater variation in ${\mathrm{\delta }}^{18}\mathrm{O}$. The black line is an 800 kyr LOESS regression of the ${BSL}_{R_1}\left(t\right)$ record that can be compared with an 800 kyr LOESS regression of the Miller et al. (2020) BSL estimate (blue line) that was calculated identically to ${BSL}_{R_1}\left(t\right)$. Sea-level is referenced to 0 in the modern, though the LOESS regressions fit through the mean of last glacial-interglacial cycles. The red vertical line indicates ice-free at 71 m. Panel (D) Deep-sea Mg/Ca paleotemperature estimates (sites in legend). The bright red line shows paleotemperature variation if the amplitudes of Milankovitch-scale sea-level variations were modulated using Eq. 9 and Eq. 10 (${pT}_{R_2}\left(t\right)$), effectively reducing the variation in ${\mathrm{\delta }}^{18}{\mathrm{O}}_{\text{sw}}$ when the long-term estimate of continental ice volume (${BSL}_{R_{1,LT}}\left(t\right);$ Eq. 8) is low. The black line is a LOESS regression of the ${pT}_{R_2}\left(t\right)$ record (red line) with an 800 kyr search window. Panel (E) The ${BSL}_{R_2}\left(t\right)$ record, light blue line, is the sea-level variation that corresponds with the ${pT}_{R_2}\left(t\right)$ paleotemperature record in Panel (D) The black line is an 800 kyr LOESS regression of the ${BSL}_{R_2}\left(t\right)$ record that can be compared with an 800 kyr LOESS regression of the Miller et al. (2020) BSL estimate (blue line) that was calculated identically to ${BSL}_{R_1}\left(t\right)$. The red vertical line indicates ice-free conditions at 71 m. Pleist., Pleistocene; Plio., Pliocene.

Our strategy and that of Miller et al. (2020) differs from recent δ¹⁸O_b compilations (De Vleeschouwer et al., 2017; Westerhold et al., 2020) because we generally restricted our compilation to deep Pacific sites (>2 km paleodepth; the exception of Site 747 is discussed below) that lie in Pacific deep water (PDW). Sites located in PDW provide a monitor of mean deep water (MDW) changes because PDW comprises ∼60% of the modern ocean water, δ¹⁸O, and δ¹³C reservoirs and up to 80% of the Paleogene (Miller and Fairbanks, 1985). We view using sites from PDW as a prerequisite to using δ¹⁸O_b to reconstruct sea level and carbon cycle changes. Stable isotopic records from Atlantic locations are potentially overprinted by changes in deep ocean circulation. Using data from the smaller Atlantic basin and its proximity to the other known deep-water sources potentially records larger temperature changes and thus compromises isolating the ice volume δ¹⁸O_sw term. Because they have a larger temperature component, Atlantic δ¹⁸O records amplify orbital variations in δ¹⁸O_b that were the target of previous studies and are thus well-suited for these studies (De Vleeschouwer et al., 2017; Westerhold et al., 2020). However, by not using Atlantic data, our δ¹⁸O_b compilation is not as complete as previous splices in some intervals, though our target was primarily the Myr-scale sea-level and climate change, not the short-period Milankovitch changes.

New Mg/Ca Data

Previous benthic foraminiferal Mg/Ca syntheses (Cramer et al., 2011; Miller et al., 2020) relied heavily on determining deep sea paleotemperatures using Mg/Ca data from Shatsky Rise Site 1,209 for the Paleocene to Middle Eocene (Dutton et al., 2005; Dawber and Tripati, 2011), though data from Southern Ocean Kerguelen Plateau Site 757 (Billups and Schrag, 2003) and Southern Ocean Site 689 (Lear et al., 2000) were included. As noted by Cramer et al. (2011), data density was low in the Early Eocene (i.e., only 29 analyses from 58 to 51 Ma), and variability older than 48 Ma raised concerns about contamination by dolomitic overgrowth on the tests in the Early Eocene and older samples.

To address these concerns, we generated a new record of benthic foraminiferal Mg/Ca for the uppermost Paleocene to Lower Eocene at Deep Sea Drilling Project Site 577 located on the Shatsky Rise in the Pacific (Figure 2; Supplementary Table S1; Supplementary Figure S1) at a paleolatitude of ∼18°N and paleodepth of ∼2000 m (Miller et al., 1987). Site 577 is close geographically (75 km) to Site 1,209 and both are in PDW, though burial depths for coeval sections are less at Site 577 (65–85 m below seafloor). Approximately 5–30 specimens of Oridorsalis umbonatus (>150 μm) were picked from faunal census slides of Miller et al. (1987).

Foraminiferal tests were gently crushed between two glass plates to open the chambers to facilitate cleaning (Boyle and Keigwin, 1985; Rosenthal et al., 1997). Trace element analyses were measured on a Thermo Finnigan Element XR Sector Field Inductively Coupled Plasma Mass Spectrometer (SF-ICP-MS) operated in low resolution (m/Δm = 300) and medium resolution (m/Δm = 4,300) settings (Rosenthal et al., 1997). Samples with unreasonably low or high calcium elemental ratios (B/Ca, Mg/Ca, Al/Ca, Ti/Ca, Mn/Ca, Sr/Ca, Fe/Ca) were considered to be contaminated or possibly altered due to diagenesis and thus rejected from the reconstructions. Elemental ratios were monitored and corrected for matrix effects for each analytical run by analyzing a suite of standards using an internal spiked gravimetric standard (SGS) with the same elemental ratios but varying Ca concentrations (1.5 mmol/L to 8 mmol/L; Rosenthal et al., 1999). Ca concentrations in sample solutions were typically maintained in the range of 1.5–4 mmol/L to minimize matrix effects on Mg/Ca.

We generated new bottom water temperature estimates for 28 samples from Hole 577 from 66 to 86 m below seafloor and used magnetobiostratigraphic datum levels to correlate to the GPTS with ages of 57.1–47.6 Ma (Supplementary Figure S1; Supplementary Table S1). Bottom water temperatures were calculated from Mg/Ca ratios of benthic foraminifera using the calibration of Lear et al. (2002) and a Mg/Ca seawater of 2.5 mmol/mol (Farkaš et al., 2007). We added the Site 577 data into the Cramer et al. (2011) compilation but restricted that compilation to Pacific and Southern Ocean sites, still nearly doubling the previous data density for the Early Eocene.

Method of BSL Calculation

We generated an estimate of barystatic sea-level (BSL; Figures 3, 4) from: 1) the Pacific benthic foraminiferal δ¹⁸O splice (δ¹⁸O_b) of Miller et al. (2020); 2) bottom water paleotemperatures, $pT$, derived from Mg/Ca at the time, $t$, of each δ¹⁸O_b measurement, $i$, derived using an iterative process that leverages Mg/Ca measurements of Pacific benthic foraminifera and the δ¹⁸O_b record (Figure 3); 3) the paleotemperature equation for benthic foraminifera that produces an estimate in the variation in the δ¹⁸O values of seawater (δ¹⁸O_sw; Eq. 1) given $pT$ and ${\mathrm{\delta }}^{18}{O}_b$ of Cibicidoides ((Lynch-Stieglitz et al., 1999);

and 4) a calibration factor for the conversion of δ¹⁸O_sw change to sea-level change, taken as 0.13‰/10 m (Winnick and Caves, 2015). This calibration will vary with the δ¹⁸O composition of ice, a factor we address in the Discussion.

FIGURE 4. Comparison of inputs into our GMGSL estimate: 1) BSL from Figure 3B (black = smoothed line, blue = all points; 2) temperature component theta-a; and 3) ocean basin volume sea level (OBVSL; red curve) after Wright et al. (2020) showing pink error bars based on their estimate of age uncertainties associated with plate tectonic reconstructions. The red line indicates ice-free.

Paleotemperatures from Cenozoic Mg/Ca Records

The first step in deriving a record of bottom water paleotemperatures through time was to process each measurement, $j$, of benthic foraminiferal Mg/Ca ratio as described by Cramer et al. (2011). Specifically, using Eq. 2, the Mg/Ca data were corrected for: 1) carbonate ion saturation effect that varies with depth, $d$, of deposition; 2) the variation in carbonate saturation through time that is tied to variations in the Calcite Compensation Depth relative to modern (as defined by Pälike et al., 2012); and 3) the variation in the Mg/Ca ratio of seawater through time (as defined by Farkaš et al., 2007) that affects Mg/Ca ratios measured on benthic foraminifera tied to parameter $H$ (set to 0.7 as per Cramer et al., 2011). An estimate of paleotemperatures, ${pT}_{Mg/Ca}\left(t_j\right)$, was then generated using the corrected Mg/Ca ratios, ${Mg/Ca}_{corr}\left(t_j\right)$, and Eq. 3 (Equation 7b in Cramer et al., 2011).

A LOESS (locally estimated scatterplot smoothing) algorithm was applied to the record of ${pT}_{Mg/Ca}\left(t_j\right)$ using a 1.5 Myr window that truncates periods shorter than ∼500 kyr to derive an initial estimate of bottom water paleotemperature (Figure 3), ${pT}_{R_0}$ (t) (Eq. 4). Estimates of ${pT}_{Mg/Ca}\left(t_j\right)$ generated using ${Mg/Ca}_{bf}\left(t_j\right)$ data from Miocene Southern Ocean Site 747 were reduced by 2.5°C to account for the relatively shallow paleodepth of the site (∼1,500–1,700 m) as defined by differences in the mean between Site 747 and Pacific sites in the intervals of overlap.

BSL Calculation

The initial estimate of bottom water paleotemperature derived using the LOESS regression of the ${pT}_{Mg/Ca}\left(t_j\right)$ estimates, i.e., ${pT}_{R_0}\left(\mathrm{t}\right)$, was used within the paleotemperature equation (Eq. 3), interpolating values at times $t_i$ to derive an initial estimate of ${\mathrm{\delta }}^{18}{O}_{sw}\left(t_i\right)$. This initial estimate of ${\mathrm{\delta }}^{18}{O}_{sw}\left(t_i\right)$ was used to create an initial sea-level curve, ${BSL}_{R_0}\left(t_i\right)$, applying a ${\mathrm{\delta }}^{18}{O}_{sw}$ to sea-level calibration of 0.13‰/10 m (Eq. 5) (Figure 3C). In generating this initial curve, we follow Miller et al. (2020) in initially assuming 20% of the variation in ${\mathrm{\delta }}^{18}{O}_{sw}\left(t_i\right)$ to be the result of variation in bottom water temperature at Milankovitch periods of 19/23, 41, and quasi-100 kyr periods not accounted for in the long-term trend derived from the ${pT}_{Mg/Ca}\left(t_j\right)$ record. This assumption was based on the upper limit placed by Pleistocene Milankovitch changes of 33% of the temperature affect (Fairbanks, 1989), though the value used (20%–33%) is somewhat arbitrary and we test this assumption. We calculate the amount of temperature change required to reduce the variation in ${\mathrm{\delta }}^{18}{O}_{sw}\left(t_i\right)$ by 20% across the entire curve and account for this residual through a revision of the initial paleotemperature record, ${pT}_{R_0}\left(\mathrm{t}\right)$, to generate a paleotemperature curve that is supplemented by the record of ${\mathrm{\delta }}^{18}{O}_b$ variations (${pT}_{R_1}\left(\mathrm{t}\right)$; Eq. 6). The variation in sea level is then recalculated (Eq. 7), producing ${BSL}_{R_1}\left(t_i\right)$ (Figure 3C).

We test our initial assumption of 20% temperature contribution on Milankovitch scales through a recomputing our ${BSL}_{R_1}$ modeled as a function of paleotemperature and generated ${BSL}_{R_2}$. When sea level is higher, the fraction of the record of δ¹⁸O_b variation attributable to paleotemperature change is potentially higher than the 20% assumed in Figure 3C and by Miller et al. (2020) because of lower total ice volume changes. Therefore, we made a second correction to account for high-frequency variations in the estimate of paleotemperatures at higher paleotemperatures and sea-level values. The conversion of δ¹⁸O_b variation to paleotemperature is accommodated by 1) applying a LOESS algorithm with a 2 Myr window to the ${BSL}_{R_1}\left(t_i\right)$ sea-level curve (Eq. 8) to produce ${BSL}_{R_{1,LT}}\left(t_i\right)$; 2) converting the fraction of the residuals of the ${BSL}_{R_1}\left(t_i\right)$ (Figure 3C) relative to ${BSL}_{R_{1,LT}}\left(t_i\right)$ defined in Eq. 9 to variation in paleotemperature using Eq. 10 to produce a second revision of the paleotemperature record, ${pT}_{R_2}\left(t\right).$ (Figure 3D) In our model, the fraction of high-frequency sea-level variation converted to paleotemperature varies as a function of ${BSL}_{R_{1,LT}}\left(t_i\right)$, as defined in Eq. 9, and scales linearly from 0% at ${BSL}_{R_{1,LT}}\left(t_i\right)$ values below −50 m to 85% at ${BSL}_{R_{1,LT}}\left(t_i\right)$ values of 50 m, and again scales linearly from 85% at ${BSL}_{R_{1,LT}}\left(t_i\right)$ values of 50 m–100% at ${BSL}_{R_{1,LT}}\left(t_i\right)$ of 75 m. The variation in sea level is then recalculated (Eq. 11) producing ${BSL}_{R_2}\left(t_i\right)$, our best estimate of BSL change shown as BSL_R2 (Figure 3E).

We combine the sea-level variations generated by transfers of water mass from ice sheets to the oceans (BSL_R2 Figures 3E, 4), thermal expansion/contraction of the oceans (Figure 4), and the variations in the shape of the ocean basins (OBVSL: Figure 4) to generate an estimate of GMGSL (Figures 5, 6). GMGSL and all data used to calculate GMGSL including Age, BSL1, BSL2, Mg/Ca paleotemperature, thermosteric sea level, OBVSL of Wright et al. (2020) are given in the Supplementary Material and have been submitted online at NCEI².

FIGURE 5. Panel (A) GMGSL estimate (black line) compared with the Mid-Atlantic US Relative Sea Level (RSL; thick red line) versus the GTS 2020 time scale. Panel (B) Statistically derived estimates of Cenozoic pCO₂ from alkenone (yellow dots), boron (purple dots), paleosols (green dots), stomatal indices (blue dots), and ice core measurements (gray dots, only available for time 800 ka to present). The best estimate of Cenozoic pCO₂ is represented by function shown with a heavy red line that runs through the boron isotope and ice core data and is systematically offset from the other proxy measurements. From 66 Ma to 800 ka, the alkenone, paleosol, and stomatal proxy data supplement the variation of the red function, but the actual value of the pCO₂ estimate is weighted towards the values derived from the boron isotope proxy. From 800 ka to the present, the ice core measurements determine the value of pCO₂. The proxy-specific offsets were determined using a Gaussian process regression model and they are shown in Supplementary Figure S3.

FIGURE 6. Blow up of Eocene sea level showing our BSL ((A); black line smoothed and blue line data) and our GMGSL estimate from Figure 5 ((B); black line) with the Mid-Atlantic US Relative Sea Level (RSL; thick red line) versus the GTS 2020 time scale. RSL has been shifted to lower values by 50 m from 42 to 57 Ma in both panels.

Thermal Expansion and Tectonic Corrections

Here we calculate Cenozoic thermosteric sea level variation, $\mathrm{\theta }$, using the ${pT}_{R_2}\left(t\right)$ record of deep-sea paleotemperatures (Supplementary Figure S2). This is possible as deep-sea paleotemperatures are a reasonable representation of global mean surface temperatures for most of the Cenozoic (Hansen et al., 2013), with the possible exception of the Pleistocene when deep-sea temperatures approach physical constraints on a minimum temperature while global mean surface temperatures continue to decrease. To calculate $\mathrm{\theta }$, the variations in deep-sea paleotemperatures were first converted to ocean surface air temperatures using the 115 Myr regression of Valdes et al. (2021) that suggests 0.67°C of global mean surface ocean temperature change occurs for each degree C of deep-sea temperature change. The global mean surface air temperatures were converted to $\mathrm{\theta }$ using the regressions generated by Hieronymus (2019). Hieronymus (2019) suggested either 0.73 or 1.18 m of $\mathrm{\theta }$ for each degree C of global mean surface temperature warming. Applying these steps according to Eq. 12 or Eq. 13 lead to the calculation of ${\mathrm{\theta }}_a\left(t\right)$ or ${\mathrm{\theta }}_b\left(t\right)$, representative of the 1.18 m/°C or 0.73 m/°C regression constant for the conversion of surface temperature to $\mathrm{\theta }$, respectively. Alternatively, $\mathrm{\theta }$ can be calculated using the deep-sea ocean temperature variation as a representation of global mean ocean temperature variation. Every degree of global mean ocean temperature results in 0.7 m of sea-level rise due to thermal expansion (Eq. 14; Hieronymus, 2019) (Figure 4; Supplementary Figure S2).

The variation in sea level attributable to changes in the shape of the ocean basins through time, $\text{OBVSL}\left(\mathrm{t}\right)$, was obtained from Wright et al. (2020) (Figure 4). There are large errors in the OBVSL shown (Figure 5, errors shown are those of the tectonic model of Wright et al. (2020). The Wright et al. (2020) record of $\text{OBVSL}\left(\mathrm{t}\right)$ was added to the records of ${\mathrm{\theta }}_a\left(t_i\right)$ and ${BSL}_{R_2}\left(t_i\right)$ to generate a record of $GMSL\left(t_i\right)$ (Eq. 15; Figure 5).

Modeling Cenozoic CO₂ Variations

Early Eocene global warmth and high sea level are generally ascribed to high atmospheric CO₂ levels (e.g., Anagnostou et al., 2016) and we explore these relationships with an updated analysis of CO₂ records. We revise the Bayesian hierarchical model with Gaussian process priors (Rasmussen and Williams, 2006; Ashe et al., 2019) used by Miller et al. (2020) to estimate atmospheric CO₂ concentration over time. To model pCO₂ herein, we use the combination of paleosol and stomata proxy data compiled by Foster et al. (2017) and measurements from Antarctic ice cores (Bereiter et al., 2015) that were used by Miller et al. (2020), while replacing the boron and alkenone proxy data from Foster et al. (2017) with the updated datasets made available by Rae et al. (2021). The statistical model constructed by Miller et al. (2020) provided an estimate of atmospheric CO₂ concentrations through time based on the assumption that all proxy data within the Foster et al. (2017) compilation, i.e., measurements from alkenones, boron, liverworts, paleosols, and stomata, provided an equally accurate representation of atmospheric paleo-CO₂ concentrations. In this study, we model systematic differences between the less reliable alkenone, stomata, and paleosol measurements and the more accurate ice core and boron proxy measurements. In our revision, we treat the ice-core measurements and boron isotope proxy as representative of the true value of pCO₂ and allow for the model to account for temporally varying, statistically identifiable systematic differences between the other proxies that improve the overall fit of the model to all the data. Therefore, we can re-project the variations recorded by the variety of proxies accounting for systematic, proxy-specific offsets to probabilistically estimate the value of pCO₂ through geologic time as if it were being measured by ice-core measurements or the boron proxy.

To account for systematic differences between the pCO₂ measurement proxies, the CO₂ values, $y_i$ (Eq. 16), are modeled as the sum of the modeled CO₂ concentration, $f\left(t_i\right);$ the proxy-specific offset, $g_k\left(t_i\right)$; and the measurement value errors, ${\epsilon }_i^y$. We do not consider the temporal error associated with individual data points.

At the process level, the best estimate for the concentration of pCO₂ in the atmosphere through time is represented by a composite function, $f\left(t\right)$ (Eq. 17), comprising four subcomponents. This function, $f\left(t\right)$, provides an estimate of pCO₂ by leveraging variations observed within all proxy data within the Foster et al. (2017) database and ice core records (Bereiter et al., 2015). The shared variations observed within all data were used despite variations in absolute pCO₂ represented by individual proxies. These proxy-induced differences and other noise are modeled by composite function $g_k\left(t\right)$ (Eq. 18).

The component processes are each modeled with a Gaussian process prior. Three of the component functions that comprise $f\left(t\right)$ ($m\left(\mathrm{t}\right)$, $l_{0}t$, and $c_0$) are representative of a signal shared by all proxy data, and the final component (${\Delta }_{\text{ice}}\left(t\right)$) represents the precise fluctuations recorded by ice-core measurements over the past 800 kyr. This shared component, $m\left(\mathrm{t}\right)$, is a non-linear process that represents long-term variations in pCO₂. Additionally, a linear process, $l_{0}t$; and an offset, $c_0$, represents the long-term trend in pCO₂ as it is represented by the entirety of the pCO₂ data compilation, $y_i$. A second non-linear process, ${\Delta }_{\text{ice}}\left(t\right)$, captures the high-frequency changes that are only captured by the ice-core record, which only spans the past 800 kyr.

The function $g_k\left(t\right)$ accounts for the systematic, linear proxy-specific differences between the less reliable alkenone, stomata, paleosol, and liverwort proxies and the paleo-pCO₂ represented by the boron proxy data and ice-core measurements. Specifically, these differences are captured by a temporally linear term, $l_{k}t$, and an offset, $c_k$. The values of $k$, [1, 2, 3], denote the expression of $l_{k}t$ and $c_k$ specific to the alkenone, paleosol, and stomata data, respectively. Since there is no modeled component representative of a systematic difference between $f\left(t\right)$ and the boron proxy or ice-core measurements, the components of $f\left(t\right)$, contain values of pCO₂ through geological time that closely represent those boron proxy and ice-core measurements that are more reliable records of paleo pCO₂. Functions $f\left(t\right)$, $g_k\left(t\right)$, and their components are portrayed in Figures 1, 5. The Supplementary Material provides detailed information on the statistical modeling methods used to estimate Cenozoic atmospheric CO₂ concentrations, including information about model parameters (Supplementary Table S2).

Modeling Covariance of BSL and the NJ Margin Backstripping Record

We apply a Bayesian hierarchical model with Gaussian process priors to evaluate variation in sea level as recorded by ${BSL}_{R_2}\left(t_i\right)$ and backstripping studies of the NJ continental margin, represented by $a_i$ (Eq. 19) and $b_j$ (Eq. 20), respectively. The primary goal of this analysis was to evaluate Myr-scale covariance between the BSL record we produce and the sea-level record from the NJ continental margin accounting for variations likely attributable to OBVSL. The sea-level variations recorded by ${BSL}_{R_2}\left(t_i\right)$ and backstripping studies of the NJ continental margin were modeled over time, t, as function $g\left(t\right)$ (Eq. 21) and $h\left(t\right)$ (Eq. 22), respectively. These two modeled functions share both linear and non-linear components, $m\left(t\right)$ and $l\left(t\right)$, that are likely to represent BSL variations that are captured by both sets of measurements.

Both sets of measurements, $a_i$ (Eq. 19) and $b_j$ (Eq. 20), contained observational errors (${\epsilon }_i^a$ and ${\epsilon }_j^b$). Miller et al. (2020) previously calculated those errors to be normally distributed with a standard deviation of ±13 m for the ${BSL}_{R_2}\left(t_i\right)$ estimate and ±15 m for the backstripped NJ record. The component processes within $g\left(t\right)$ and $h\left(t\right)$ are each modeled using a mean-zero Gaussian process prior and are presumed to encapsulate discrete physical processes. The linear terms, $lt$ and $l_{h}t$, define the long-term variations within each dataset. The $l\left(t\right)$ component, shared between both records, signifies long-term ice-volume changes and the trend of variation in BSL. In contrast, $l_{h}t$ is unique to the backstripped record, $h\left(t\right)$, and accounts for long-term trends in regional vertical land motion, global tectonic influences on GMGSL, and the long-term trend in sea level due to thermal expansion. The non-linear $m\left(t\right)$ term shared by both $g\left(t\right)$ and $h\left(t\right)$ captures BSL variation attributable to changes in continental ice volume. The other non-linear term, $\Delta \left(\mathrm{t}\right)$, is unique to $h\left(t\right)$ and represents the differences between the sea-level changes recorded by ${BSL}_{R_2}\left(t_i\right)$ (i.e., $g\left(t\right)$) and those recorded via backstripping (i.e., $h\left(t\right)$) that are correlated over approximately a million-year scale. The white-noise terms (w_g(t) and w_h(t)) capture high-frequency variability and measurement error beyond that accounted for in the nominal error estimates. The constant terms ($c_1$ and $c_2$) capture differences in time-average sea-level estimates at each site.

To fit the model, we applied a Markov Chain Monte Carlo (MCMC) algorithm. The MCMC algorithm was used to estimate the values of hyperparameters that maximized the likelihood of $g\left(t\right)$ and $h\left(t\right)$ considering the underlying data, $a_i$ and $b_j$. Through an iterative process, the MCMC method samples a posterior distribution of optimized hyperparameter values (Supplementary Table S3) that could be used to generate posterior estimates of $g\left(t\right)$ and $h\left(t\right)$ and their constituent processes (see Discussion). The posterior estimates of $g\left(t\right)$ and $h\left(t\right)$ facilitated a robust comparison of the ${BSL}_{R_2}\left(t_i\right)$ estimate and the backstripped NJ margin sea-level record, shown as $h\left(t\right)-g\left(t\right)$. The Supplementary Material provides detailed information on the statistical modeling for the statistical modeling of ${BSL}_{R_2}$ and backstripped NJ margin sea-level record.

Results

We note that our methods are similar to and generally reproduce the bottom water temperature record of Cramer et al. (2011) for the last 34 Ma, albeit with better resolution in the Early Eocene. Our estimates in the Eocene are similar to but slightly cooler than Cramer et al. (2011), especially older than 48 Ma, which we attribute to the new Site 577 data. PDW (representing the bulk of the deep ocean reservoir) temperatures in the Site 577 record ranged from 5.0°C to 11.9°C (Supplementary Figure S1), with a mean of 9.1°C, that are ∼1°C cooler temperatures than the Early Eocene Cramer et al. (2011) compilation (Supplementary Figure S1), though both recorded spikes of ∼15°C representing hyperthermals. Earlier Paleocene (>57 Ma) Mg/Ca data are still lacking, and thermal history (hence sea level) is only constrained by a few points, though limited clumped isotope data from Meckler et al. (2022) (Figure 7) support our assumed earlier Paleocene (>57 Ma) bottom water temperatures of 9°C–10°C (see Discussion).

FIGURE 7. Sensitivity tests of BSL estimates. Panel (A) Comparison of three deep-sea Mg/Ca paleotemperature estimates. The black line plots the long-term ∼1.2 Myr trend of the Eq. 7b Mg/Ca calibration from Cramer et al. (2011), which was used by Miller et al. (2020) to calculate BSL. The Eq. 7b Mg/Ca calibration implies an influence of Mg/Ca_sw on benthic foraminiferal Mg/Ca ratios. The data points in the plot represent ${pT}_{Mg/Ca}\left(t\right)$, and the blue line shows the long-term trend of ${pT}_{Mg/Ca}\left(t\right)$ generated using a 1.5 Myr LOESS regression. This curve represents the Eq. 7b Mg/Ca calibration from Cramer et al. (2011), with an updated compilation of deep-sea Mg/Ca data and a less intricate trend algorithm that produces a closer fit to the Mg/Ca data points. The bright red line represents ${pT}_{R_2}\left(t\right)$, which is an estimate of the high-resolution temperature record derived using the initial temperature estimate, the ${\mathrm{\delta }}^{18}{O}_b$ splice, and a modulation of variation in ${\mathrm{\delta }}^{18}{O}_{sw}$ when the long-term estimate of continental ice volume (${BSL}_{R_{1,LT}}\left(t\right);$ Eq. 8) is low. The dark red line shows a 1.5 Myr LOESS regression of the Mg/Ca data compilation used to calculate ${pT}_{Mg/Ca}\left(t\right)$, but processed using the Lear et al. (2015) Mg/Ca-T calibration that implies a low sensitivity of benthic foraminiferal Mg/Ca to seawater Mg/Ca. Panel (B) The ${BSL}_{R_2}\left(t\right)$ record, light blue line, is the sea-level variation that corresponds with the ${pT}_{R_2}\left(t\right)$ paleotemperature record shown in the left panel. The black line is an 800 kyr LOESS regression of that ${BSL}_{R_2}\left(t\right)$ record. The black line corresponds to an 800 kyr LOESS regression of ${BSL}_{R_2}\left(t\right)$ if the Cramer et al. (2011) Eq. 7b Mg/Ca-T calibration was used to calculate paleotemperature. Similarly, the dark red line shows an 800 kyr LOESS regression of ${BSL}_{R_2}\left(t\right)$ if the Lear et al. (2015) Mg/Ca-T calibration was used to calculate paleotemperature. The red vertical line indicates ice-free conditions at 71 m. Panel (C) ${BSL}_{R_2}\left(t\right)$ with a ${\mathrm{\delta }}^{18}{\mathrm{O}}_{\text{sw}}$ to sea-level calibration of 0.13‰/10 m (blue line), 0.1‰/10 m (black line), and 0.05‰/10 m (dark red line) to demonstrate the sensitivity of ${BSL}_{R_2}\left(t\right)$ to end member calibrations. The red vertical line indicates ice-free conditions at 71 m.

Our update of the Cramer et al. (2011) synthesis still warrants further scrutiny, including new data and consideration of the importance of the carbonate ion effect. Lear et al. (2015) suggested the Mg/Ca_sw carbonate ion effect for Oridorsalis measurements is potentially insignificant (i.e., H = 0). The new Site 577 data provided are Oridorsalis, as are other data in the Cramer et al. (2011) compilation (including all of Site 1,218 and most of Site 1,209). To test the effect, we also computed bottom water temperatures and sea level with no correction for carbonate ion effect (i.e., H = 0) (see Discussion). The result yields bottom water temperatures that are slightly lower through the Cenozoic, though the results for the BSL changes are robust (see Discussion).

Sea level in a Purportedly Ice-free World

Sea level was higher in the absence of terrestrial ice sheets during at least parts of the Early Eocene, and thus modern storage volumes provide an upper BSL range. BSL estimates for ice-free conditions have ranged from 65 to 70 m (for discussion see Miller et al., 2005). More precise recent satellite-based estimates of BedMachine of 22.8 ± 0.5 million km³ of grounded ice sheets above sea level (above floatation) in Antarctica and of 2.99 ± 0.02 million km³ of ice in Greenland suggest a potential sea level equivalent of 58 ± 0.9 m and 7.4 ± 0.05 m in Antarctica and Greenland, respectively (Morlighem et al., 2020). Mountain glaciers contribute less than a meter to BSL change (0.32 ± 0.08 m; Farinotti et al., 2019). This total of 66 m modern sea-level equivalent does not include ice grounded below sea level (below floatation) of 6.6 million km³ (Morlighem et al., 2020) that would have raised sea level once isostatic compensation is realized and affected δ¹⁸O_sw. The 6.6 million km³ of ice stored below floatation in Antarctica (∼8 m) would yield ∼5.4 m of rise after isostatic compensation yielding a total sea level rise of ∼71 m if all modern ice were melted (Morlighem et al., 2020).

The 71 m ice-free estimate provides an important constraint on the δ¹⁸O-Mg/Ca based BSL estimates (Figure 4). The BSL sea-level record of Miller et al. (2020) and our uncorrected BSL₁ (Figure 3C) estimates show long-term (Myr scale) smoothed BSL estimates reaching this level, and numerous high frequency (20–100 kyr scale) events far exceeding it (up to 100 m above present) because these records have not been corrected for short term (<1 Myr) bottom water temperature variations. The times when these records exceed the 71 m upper limit are hyperthermal warming events (Figure 2C), as noted but uncorrected by Miller et al. (2020); for example, a 4°C hyperthermal warming would be expressed as greater than 60 m of artifactual sea level rise (assuming 0.2‰–0.25‰/°C and 0.13‰/10 m). Our corrected sea-level BSL₂ (Figure 3E) accounts for high-frequency variations of bottom water paleotemperatures at higher paleotemperatures and sea-level values and produces a BSL curve that does not exceed the ice-free limit and correctly shows the hyperthermal events not as ice-driven sea-level events but as thermal events (Figures 3E, 4).

It is appropriate to compare our BSL₂ estimate (which is essentially ice volume) with the 71 m ice-free constraints (Figures 3, 4). We see several time intervals when ice-free conditions were attained: Early Paleocene (ca. 64 Ma), Late Paleocene (ca. 57.5 Ma), several times during the Early Eocene (54, 51, 49.5, 48 Ma), the Middle Eocene Climate Optimum (40 Ma), and the Late Eocene (ca. 35 Ma; Figures 3, 4). Though the MCO was almost ice free, the BSL estimate still suggests ∼10–20 m of sea level equivalent ice volume (Figures 3, 4), supporting studies that suggest greater storage of ice in West Antarctica during the MCO and a smaller, but persistent ice sheet (Marschalek et al., 2021). Ice volume exceeded modern values numerous times in the Oligocene to Early Miocene (Figure 4). Wilson and Luyendyk (2009) have suggested that the area of the West Antarctic continent was larger in the later Eocene to Miocene than today due to rifting and subsequent thermal subsidence and erosion. They estimated that total Antarctic land area was 10% (e.g., capable of storing 7 m more SLE based on area alone) to 20% larger (capable of storing 12 m more SLE; Wilson and Luyendyk, 2009) than at present (Wilson and Luyendyk, 2009), consistent with previous estimates of an Oligocene ice sheet at times larger than modern (Zachos et al., 1996; Miller et al., 2005; Petersen and Schrag, 2015).

It is intriguing that the Early Eocene partly fits preconceived notions of an ice-free world, but even then, it appears that ice-free conditions existed for only about half of the 8 Myr of the Early Eocene (Figure 4). The Early Eocene peaks in BSL and ice-free conditions at 54, 51, 49.5, 48 Ma were followed by four BSL falls of 40, 30, 20, and 40 m, respectively (Figures 4, 6). Precise rates cannot be computed given the several 100 kyr smoothing used to construct this curve (peaking at 800 kyr, truncating periods shorter than ∼500 kyr) and our estimates have large error bars (±10 m). Nevertheless, BSL rates were at least 20–70 m/Myr (2–7 cm/kyr). Only ice-volume changes can explain these large, rapid changes, and we argue that ice-volume events punctuated the purportedly ice-free Early Eocene.

As previously shown, the Middle to Late Eocene was punctuated by several major sea-level falls (Figure 6; Browning et al., 1996; Miller et al., 2020) associated with glacial deposits on the East Antarctica continental shelf (Figure 1; Gulick et al., 2017). The Middle Eocene Climate Optimum (MECO; Bohaty et al., 2009; 40 Ma, GTS 2020) was preceded by several 15–25 m BSL variations and bracketed by a 40 m rise and a 40 m fall, with rates of BSL change 30–75 m/Myr (Figure 6). The ca. 37 Ma sea level fall after MECO was followed by a slow rise and another sharp BSL peak of ∼72 m at ca. 35 Ma that implies ice-free conditions, the Late Eocene Climatic Optimum (LECO; Figure 1). This short ice-free interval was followed by the Oi1 glaciation at ca. 33.8–34 Ma, greater than 60–80 m sea-level fall, and the beginning of the Icehouse World (Figure 6).

Relationship of δ¹⁸O_sw and Sea Level

We use a sea-level/δ¹⁸O_sw calibration of 01.3‰/10 m developed for the warm Pliocene (Winnick and Caves, 2015) that is supported by Mg/Ca bottom water paleotemperatures, by the relationships between modern δ¹⁸O_sw and ice volume, and by mass balance changes. Shackleton and Kennett (1975) estimated that melting modern ice sheets would contribute about −0.9 to −1.0‰ to δ¹⁸O_sw, with an implied mean ice sheet value (δ¹⁸O_ice) _of −50‰ and calibration of 0.14‰/10 m; more recent updates do not significantly change this estimate. The change in δ¹⁸O_b since the Early Eocene peak has been measured as ∼3.7‰ (Pacific Eocene values δ¹⁸O_Cibicidodies are −0.5‰; modern deep Pacific core top δ¹⁸O_Cibicidodies values are 3.2‰; Shackleton and Kennett, 1975; Savin et al., 1975; Miller et al., 1991; Miller et al., 2020; Zachos et al., 2001; Westerhold et al., 2020), which requires that ∼2.7‰ of the total Cenozoic δ¹⁸O_Cibicidodies increase was due to bottom water cooling. Given a slope of 0.21‰/°C from paleotemperature equation (Eq. 1), this predicts a ∼12°C–13°C cooling that agrees with that measured using Mg/Ca data (Figure 3). A δ¹⁸O_sw change of −0.9 to −1.0‰ for 71 m of sea level change (see above) empirically requires a mean calibration of 0.11‰–0.14‰/10 m, bracketing the 0.13‰/10 m value used here.

The relationship between δ¹⁸O_sw and sea level is not fixed (e.g., Mix and Ruddiman, 1984) and must have varied with contributions from different ice sheets. The different ice sheets have different δ¹⁸O_ice values: the Greenland Ice Sheet [GIS] is −35‰, the West Antarctica Ice Sheet [WAIS] is −41‰, and the East Antarctica Ice Sheet [EAIS] is −56.5‰ (Winnick and Caves, 2015). The Last Glacial Maximum calibration of 0.11‰/10 m (Fairbanks and Matthews, 1978) requires a δ¹⁸O_ice value of −40‰ for ice including the Laurentide ice sheet (Table 1). We note that the δ¹⁸O_sw and sea level relationship assumed implicitly predicts the δ¹⁸O_ice (Table 1). It might be assumed that on a warmer planet δ¹⁸O_ice would have higher values due to less meridional fractionation (Miller et al., 1991); higher values of δ¹⁸O_ice mean a larger sea-level change for a given δ¹⁸O_sw change. However, for the warm Pliocene, the higher contribution from WAIS and EAIS versus Laurentide (0.11‰/10 m) result in a higher calibration (0.13‰/10 m) and lower implied δ¹⁸O_ice even though the WAIS-EIAS δ¹⁸O_ice values were 1‰–4‰ higher than they are today (Winnick and Caves, 2015). The 0.13‰/10 m calibration is likely to be generally representative from 33 Ma and 2.55 Ma (Winnick and Caves, 2015). Previous studies have used 0.1‰/10 m for the Oligocene (Pekar et al., 2002) and we illustrate here that this yields higher amplitudes of sea level change.

TABLE 1. Relationship of δ¹⁸O_sw-sealevel calibration to δ¹⁸O_ice and sea-level changes for given δ¹⁸O_sw variations assuming mean depth of ocean of 3,700 m.

We acknowledge that further work modeling the δ¹⁸O_sw to sea-level calibration is warranted, but it will not alter our fundamental conclusion of significant Early Eocene ice sheets that are indicated by empirically derived δ¹⁸O_sw changes of ∼0.25‰ and independently derived backstripped RSL changes (Figure 6). For reference, the 0.13‰/10 m δ¹⁸O_sw to SL calibration generates ∼20 m of SL change for a 0.25‰ variation in δ¹⁸O_sw (Table 1). By using a 0.13/10 m δ¹⁸O_sw to BSL calibration, we likely underestimate the actual variation of Eocene sea level as a function of δ¹⁸O_sw change prior to the significant development of the continental scale Antarctic ice sheet (i.e., before 34 Ma). The calibration during this time likely could possibly range from 0.13 (our preferred) to 0.05‰ per 10 m (e.g., Rohling et al., 2022), generating a 20–54 m SL variation for a δ¹⁸O_sw variation of 0.25‰ (Table 1; Figure 7).

Uncertainties in Long-Term Tectonic Record

Our GMGSL estimate provides a working model for GMGSL change, though there are still considerably large errors in the long-term record (OBVSL; Figure 4) due to uncertainties in sea-floor production rates older than 40 Ma. As first shown by an error analysis by Kominz (1984), this is because about half of the ocean crust older than 55 Myr has been destroyed, and this is shown by increasing uncertainties in the OBVSL in the Eocene (Figure 4). As noted by Rowley (2002), any attempt to reconstruct ocean crust production rates for the Cretaceous to early Paleogene using plate tectonic models is fraught with uncertainty (errors shown in Figure 4 are those of the tectonic model of Wright et al., 2020). Van der Meer et al. (2022) obtained an estimate of long-term sea level by modeling changes in the Phanerozoic Sr-isotope record that suggest even higher tectonically driven sea levels (Figure 9). As shown here (Figures 5, 6), Early Eocene sea levels were certainly much higher than modern: the thermal expansion component was ∼6 ± 2 m higher (Supplementary Figure S2), the BSL component was 71 m, and the ocean basin component was ∼60 ± 80 m, for a peak of ∼143 ± 80 m (based on plate tectonic models).

Comparison With US Atlantic Continental Backstripped RSL

We have previously reported RSL estimates from backstripping (accounting for compaction, loading, and thermal subsidence) of the onshore US Atlantic continental margin Cretaceous to Miocene (New Jersey and Delaware; Miller et al., 2005) and offshore Miocene (Australian Marion Plateau, John et al. (2011); New Jersey continental shelf; Kominz et al., 2016). Backstripping of passive continental margins provides an RSL estimate that reflects GMSL, non-thermal tectonism, and other effects (e.g., oceanographic, rotational, gravitational) with errors of approximately ±10 m (Kominz et al., 2008). Miller et al. (2020) and Schmelz et al. (2021) applied a Bayesian hierarchical model with Gaussian process priors and established that their BSL estimates share a component that covaries on Myr scale with RSL estimates from backstripping of the US Mid-Atlantic Coastal Plain. This component is attributable to BSL sea level change but was only extended to 48 Ma due to uncertainties outlined above (Miller et al., 2020). Here, we update the comparison of BSL and RSL calculated by Miller et al. (2020); Figure 8.

FIGURE 8. Modeled estimates of sea-level change, recorded by the δ¹⁸O-Mg/Ca BSL proxy (Panel (A)) and by backstripping (RSL) studies from the US Mid-Atlantic margin (Panel (B)). The modeled estimates of BSL and the RSL record (i.e., $\mathrm{g}\left(\mathrm{t}\right)$ and $\mathrm{h}\left(\mathrm{t}\right)$, respectively) represents a combination of component processes that have been fitted using Gaussian process priors. A shared non-linear component, $m\left(t\right)$, represents the variation in sea level attributable to ice-volume changes that are captured by both the global and RSL margin records (Panel (C)). Another non-linear term, Δ(t), represents the million-year scale difference between the global BSL record and the US Mid Atlantic margin sea-level record (Panel (D)). The linear components $l\left(t\right)$ and $l_h\left(t\right)$, which represent the long-term (10 Myr-scale) trends of each sea-level record, demonstrate the long-term difference between the two records (Panel (E)). The difference between the BSL and RSL records is shown in full in Panel (F). The difference between records is minimal from ∼40 Ma to 10 Ma. The Myr-scale differences in the Miocene (Panels (D, F)) are attributable to mantle dynamic topography (MDT; Moucha et al., 2008; Kominz et al., 2016; Schmelz et al., 2021).

With our revised BSL (Figures 3, 4) and GMGSL (Figure 5) we here show comparisons of these independently derived sea-level estimates with backstripped RSL records for the entire Cenozoic (Figure 5), focusing on the Eocene (Figure 6). All show similar large Myr scale changes, even in the Early Eocene (Figure 6). The rises and falls of sea-level are faithfully mimicked in deep sea and continental margin records (Figure 6), even though both estimates (BSL and backstripping RSL) have large errors (±10–15 m in amplitude and ±0.5 Myr in age for margin studies). The BSL estimate is offset from the adjusted (50 m shifted) RSL estimate by ∼50 m (Figure 6A), reflecting primarily the OBVSL component as discussed below. Both estimates qualitatively show typical Myr scale variations of over 20–60+ m driven by ice-volume variations in the Icehouse. In the case of the Early Eocene, the changes of 20–40 m in much less than a million years in both records can only be explained by growth and decay of large (25%–60% of modern) ice sheets.

The U.S. mid-Atlantic margin RSL has been impacted by changes in MDT on the several-million-year scale. Miocene onshore U.S. mid-Atlantic backstripped records differ from offshore NJ backstripped records due to a non-thermal component attributable to MDT (Kominz et al., 2016). Schmelz et al. (2021) modeled these differences onshore-offshore and regionally and showed that these can be explained by regionally propagating changes in relative sea level during the Miocene due to MDT. We also quantify differences in the Eocene (Figure 8F) of ∼50 m as suggested by Miller et al. (2020).

We present a Bayesian analysis with Gaussian process priors that shows our new BSL estimate shares a component that covaries on Myr scale with RSL estimates from backstripping (Figure 8). We updated Bayesian comparisons of Miller et al. (2020) and model estimates of BSL (Figure 8A) and those from backstripping studies from the NJ margin (Figure 8B). The modeled estimates of BSL and the NJ SL record (i.e., $\mathrm{g}\left(\mathrm{t}\right)$ and $\mathrm{h}\left(\mathrm{t}\right)$, respectively) represent a combination of component processes that have been fitted using Gaussian process priors. A shared non-linear component, $m\left(t\right)$, represents the variation in sea level attributable to ice-volume changes that are captured by both the BSL and NJ margin records (Figure 8C). Another non-linear term, Δ(t), represents the million-year scale difference between the BSL record and the NJ margin sea-level record (Figure 8D). This component indicates that the BSL component shared by the two records are similar within ±16 m since 42 Ma, and ±26 m through the Cenozoic. The linear components $l\left(t\right)$ and $l_h\left(t\right)$, which represent the long-term (10 Myr-scale) trends of each sea-level record, demonstrate the long-term difference between the two records (Figure 8E). The long-term decrease in BSL is 1.01 ± 0.26 m/Myr. The difference between the two records is minimal from ∼40 Ma to 10 Ma. The Myr-scale differences in the Miocene (Figures 8D, F) have been attributed to mantle dynamic topography (MDT; Moucha et al., 2008; Kominz et al., 2016; Schmelz et al., 2021).

We quantify differences between the BSL and the NJ margin records earlier than ∼40 Ma of about 50–100 m (Figure 8D) that have been previously suggested. Miller et al. (2020) noted that the onshore RSL estimates overlay the δ¹⁸O-Mg/Ca based BSL estimate younger than 43 Ma, but the records were offset from each other by more than 50 m from 43 to 57 Ma (Figures 5, 6). They subtracted 50 m from the mid-Atlantic record from 57 to 43 Ma (as we did in Figures 5, 6) to register with the BSL estimate derived from δ¹⁸O-Mg/Ca and attributed the offset to one of three effects:

1) an elevated level of Eocene GMGSL attributable to OBVSL (Figures 5, 6). We account for by showing the shifted RSL overlays our GMGSL record (cf. Figures 6A, B), though note large uncertainties in OBVSL do not preclude its contribution.

2) overestimates of water depth in the deep water Eocene sections (e.g., assignments of paleodepth to 150 m could be shifted to 100 m); we reject this hypothesis based on noting that Lower to Middle Eocene sediments in New Jersey (Manasquan and Shark River Formations) are relatively deeper water than coeval along-basin strata in Virginia (Nanjemoy Formation): they are finer grained, contain greater percentages of planktonic foraminifera, and have a distinctly different deeper water fauna (Browning et al., 1997).

3) MDT effects have been hypothesized to have caused over 10 m of excess subsidence since the Eocene in the onshore New Jersey sites (Müller et al., 2008; Spasojević et al., 2008). As noted, our comparisons suggest about 50 m of excess subsidence on the RSL record in the Eocene (Figures 6B, 8D) in addition to the previously modeled changed in the Miocene (Schmelz et al., 2021).

Ultimately, discerning the source of the difference is a challenge because the errors on OBVSL and models of MDT older than 30 Ma are very large, which obscures the difference between the US Mid Atlantic margin record and BSL is the OBVSL signal or MDT. These residuals of BSL to the RSL of the US Mid Atlantic margin have modeled this in detail with similar data in the past (Schmelz et al., 2021), but without incorporating process-based models of MDT. Refining the model of these residuals is an avenue of future work that might better help discern the contributions of vertical land motion on the US Mid Atlantic margin attributable to MDT, and the record of OBVSL and GMGSL.

Comparison With Exxon Sea-Level Records

Comparisons with the 1987 Exxon Production Research Company (EPR) sea-level curve (Figure 7; Haq et al., 1987; Haq and Al-Qahtani, 2005; Snedden and Liu, 2010) reinforces the previously documented fact that their amplitudes are greatly exaggerated (Miller et al., 2005; Miller et al., 2020). Eocene sea-level falls on the Haq et al. (1987) EPR curve are over 100 m, and the mid-Oligocene fall is 160 m. Our records show typical Eocene falls of 20–40 m and the Oi1 fall of 60–80 m. The greatest mismatches with the Haq et al. (1987) are in the Early Oligocene, which is a high in EPR and a low in ours, and the general fact that their amplitudes are 2.5 times as great as ours. The EPR curves were generated with the assumption that all flooding surfaces match the long-term OBVSL of Pitman and Golovchenko (1983) and that only the mid-Oligocene fall reached amplitudes greater than 120 m (Miller et al., 2020). Snedden and Liu (2010) updated EPR curve better fits our GMGSL record due to its lower amplitudes and lower Oligocene sea levels. The Wright et al. (2020) OBVSL we use is similar to Pitman and Golovchenko (1983) in the Cenozoic and thus the long-term (10 Myr scale) comparisons of EPR and ours show some similarities. However, the EPR amplitudes are strikingly different and even the relative patterns do not entirely match, though as we have noted previously, the timing of the Myr-scale sea-level falls is similar. Thus, the EPR cycle charts can still be used to predict the timing of many sea-level falls in the Phanerozoic record, though the amplitudes of sea-level changes are not meaningful.

Comparison With Other δ¹⁸O_b Based BSL Estimates

Rohling et al. (2021) and Rohling et al. (2022) provide a synthesis of sea-level and temperature variations using δ¹⁸O records over the past 40 Myr that, at first glance, appears to be strikingly different from ours (Figure 9), at least for the Oligocene to Early Miocene. Their sea-level estimates are derived using a quadratic equation that relates δ¹⁸O_b to sea-level change (see Figure 6 in Rohling et al., 2022). The quadratic equation fits GMSL values from the Spratt and Lisiecki (2016) sea-level curve and δ¹⁸O_b values in Lisiecki and Raymo (2005), data that extend to 800 ka and span sea-level values from approximately −120 m to ∼15 m greater than present. This quadratic equation is extended beyond the bounds of those data to reach an ice-free value of 65.1 m at δ¹⁸O_b value of 0.5‰. Rohling et al. (2022) generated a Cenozoic sea-level curve by applying this quadratic equation to the Westerhold et al. (2020) splice of δ¹⁸O_b. The primary assumption underlying the validity of the curve is that the quadratic equation correctly identifies the ice-free value of δ¹⁸O_b.

FIGURE 9. Panel (A) Comparison of our BSL (blue lines) and smoothed record (black line) with the sea level record of Rohling et al. (2021, 2022) (red line). Panel (B) Comparison of our GMGSL (black line) and Haq et al. (1987) (red line recalibrated to GTS 2020), Snedden and Liu (2010) (green), and van der Meer et al. (2022) (blue). Note different horizontal scales between left and right panels. The red vertical line indicates ice-free conditions at 71 m.

We show that the modeled sea-level change (Rohling et al., 2021; 2022) is similar to our ${BSL}_{R_1}$, ${BSL}_{R_2}$, and the curve generated by Miller et al. (2020) in pattern, though their amplitudes are less than half of ours and there is a systematic offset (Figure 10, note top scale is twice the bottom). The similarity of the curves, both in timing and relative scaling of distinct sea-level events, is attributable to very high covariance between the Westerhold et al. (2020) and Miller et al. (2020) records of δ¹⁸O_b that underpin the calculation of the sea-level variation in both analyses. The difference in amplitude is attributable to the method of transformation of the δ¹⁸O_b data to sea-level values. The quadratic relationship applied by Rohling et al. (2022) implies that there is a greater proportion of temperature-driven fractionation effect at higher values of δ¹⁸O_b relative to δ¹⁸O_sw change. As discussed above, this is rational because higher values of δ¹⁸O_b correspond with higher temperatures that are less conducive to the growth of continental ice volume and thus less continental ice volume is available to melt, but the consequences of their model assumptions are contradicted by observations.

FIGURE 10. Blow up of Oligocene sea level. Our BSL estimates from Figure 3 (light blue) and the smoothed record (dark blue); compared to the Rohling et al. (2021, 2022) record (red); note the two are overlain using different scales, 25 m range for Rohling, 50 m range for this study.

Long-term bottom water paleotemperature records derived from Mg/Ca data and from clumped isotopes (Figure 7A) do not support the application of the quadratic regression applied in Rohling et al. (2022) to δ¹⁸O_b values or sea level values greater than those used to establish it (i.e., the span of GMSL values from the Spratt and Lisiecki (2016) sea-level curve [-120 m to ∼15 m] and the corresponding δ¹⁸O_b values in Lisiecki and Raymo 2005 [∼3‰–5‰]). Through the Oligocene, the bottom water temperature required to support a mean sea level of 45 m of sea level in Rohling et al. (2022) given a mean δ¹⁸O_b value of 2.3‰ (Oligocene average value from Westerhold et al., 2020; Miller et al., 2020) ranges from 0.2°C to 2.9°C. Thus, the average Oligocene temperature implied by values of δ¹⁸O_sw modeled by Rohling et al. (2022) is ∼2°C. The Mg/Ca proxy (Lear et al., 2000; Cramer et al., 2011; this study) suggest that bottom water temperatures in the Oligocene were an average of 7°C (5.5°C–8.5°C, 95% CI). Independent clumped isotopes methods of estimating bottom water temperature for this period are slightly warmer (∼7°C–9°C; Meckler et al., 2022), refuting the low temperatures required by the Rohling et al. (2022) calculation of sea level and their modeling of δ¹⁸O_sw.

Differences between the methods applied to derive sea level by Rohling et al. (2022) and those we apply have implications not only for sea-level magnitudes, but also for the values of δ¹⁸O_sw and the δ¹⁸O for ice sheets (δ¹⁸O_ice). The methods are quite different. Rohling et al. (2022) calculated δ¹⁸O_sw by budgeting sea-level change through a model of ice sheet growth and decay. This model generates a variable relationship between δ¹⁸O_sw and sea level because it incorporates functions that model the increasingly lower δ¹⁸O_ice values of continental ice sheets as they grow. As a result, δ¹⁸O_sw variation of a given magnitude corresponds with larger sea-level fluctuations at higher sea-level values within the Rohling et al. (2022) model than at lower sea-level values. Conversely, following Miller et al. (2020), we initially adopt a constant calibration of 0.13‰ variation in δ¹⁸O_sw for every 10 m of sea-level change. The Rohling et al. (2022) model assumes ice-free conditions at 0.5‰; in contrast, the calculations made by Miller et al. (2020) and herein applying Mg/Ca data as a constraint on bottom water temperatures produce ice-free conditions at δ¹⁸O_b values of −0.5‰ and sea-levels of 71 m.

The assumptions made in scaling δ¹⁸O_sw to sea level have implications for δ¹⁸O_ice values (Table 1). Our assumed calibration of 0.13‰/10 m requires ice sheets with a mean of −47‰ (Table 1). The Rohling et al. (2022) model produces a relationship between δ¹⁸O_sw and sea level that ranges from 0.08‰/10 m to 0.11‰/10 m (95%ile) with a mean of 0.10‰/10 m through the Oligocene; for the Late Eocene, their model yields a calibration of ∼0.04‰/10 m. The corresponding mean value of δ¹⁸O_ice required by their model through the Oligocene is −37‰ and for the Late Eocene is −15‰. This latter value is above the upper limit for precipitation of modern ice sheets (e.g., −17‰ values found in precipitation at the southern tip of Greenland or Patagonia; Darling et al., 2006).

The Rohling et al. (2022) model is an intriguing means of using δ¹⁸O_b to generate a sea-level record, though their assumptions result in lower amplitude variability in sea level relative to given variations in δ¹⁸O_b due to the low slope of their extrapolation. The Rohling et al. (2022) model implicitly does not allow sea level to vary more than about 10–20 m for δ¹⁸O_b values between 0.5‰ and 2‰, which is the main source of their low amplitudes. As noted above, the 0.13‰/10 m calibration is likely to be generally representative from 33 Ma and 2.55 Ma (Winnick and Caves, 2015). Coupling variations in sea-level to an ice-sheet model that characterizes the temporal variability of the δ¹⁸O_sw to sea-level calibration should be investigated further. However, the assumptions associated with Rohling et al. (2022) conversion of δ¹⁸O_b to sea-level (namely, the extension of the quadratic equation beyond the bounds of data) for the pre-Middle Miocene are contradicted by various empirical datasets. Specifically, datasets contradicting the Late Eocene to Early Miocene outputs of the Rohling et al. (2022) model include Mg/Ca and clumped isotopes paleothermometry and the excellent agreement between the Myr-year scale amplitude of our sea-level curve derived from Mg/Ca data and those from backstripped estimates of RSL from the US Mid Atlantic margin (Figures 5, 6).

Testing the Sensitivity of Sea Level to Mg/Ca Assumptions and δ¹⁸O_sw Calibrations

As noted throughout, scaling of δ¹⁸O_b to sea-level requires numerous assumptions, particularly paleotemperature, to obtain δ¹⁸O_sw and the relationship of δ¹⁸O_sw to sea level. We regard the agreement between two independent estimates, backstripping and δ¹⁸O_b-Mg/Ca approaches, as the best evidence that our estimates are reasonable and valid (Figures 5, 6). However, as noted by Raymo et al. (2018), errors on both methods are large (±10 m or worse), and we test our comparisons using end-member assumptions (Figure 7).

Here, we compute BSL estimates using various assumptions (Figure 7) to illustrate end-member scenarios for sea level using the δ¹⁸O_b-Mg/Ca method and differences due to assumed sea level to δ¹⁸O_sw calibrations. Mg/Ca thermometry has ±1° to ±2° uncertainty on a given measurement considering all uncertainties, which translates into ±0.21 to ±0.42‰ δ¹⁸O_sw and ±15 to ±30 m. The Lear et al. (2015) calibration implies a low sensitivity of benthic foraminiferal Mg/Ca to seawater Mg/Ca. Lear et al. (2015) demonstrate this finding to be applicable to intermediate depth Mg/Ca Oridorsalis umbonatus data. We extrapolate this implied relationship as an end-member sensitivity test for our assumptions in applying the Cramer et al. (2011) Eq. 7b calibration, which implies a universal effect of Mg/Ca_sw on benthic foraminiferal Mg/Ca. We show that applying the Lear et al. (2015) suggestion of no carbonate ion effect yields Mg/Ca bottom water temperatures that are 2°C–4°C cooler than ours (Figure 9, left panel, red line) and thus sea levels of 30–60 m higher, peaking over 100 m (Figure 9 center panel, red line), which are not reasonable. Independent bottom water temperature estimates from clumped isotopes provide additional constraints (Meckler et al., 2022) and show much warmer temperatures than obtained with the Lear et al. (2015) suggestion and are like or slightly warmer than our best estimate. One reason for the warmer temperatures from clumped isotopes is that they were measured on North Atlantic sites. Our bottom water temperature and sea level estimates are slightly cooler than Cramer et al. (2011); Figure 7, left panel), with sea levels slightly higher. Our proposed temperature record produces the most reasonable bottom water temperature estimates compared to clumped isotopes and sea-level estimates that stay below the ice-free limit of 71 m. They are most similar to those of Cramer et al. (2011), though with better resolution in the Early Eocene (Figure 7).

We also illustrate BSL estimates with three different calibrations (0.13‰/10 m; 0.1‰/10 m, and 0.05‰/10 m; Figure 7, right panel). The 0.05‰/10 m calibration (with its implied δ¹⁸O_ice of −17‰; Table 1) yields unreasonably high sea levels and amplitudes of sea-level change. The 0.1‰/10 m calibration (with its implied δ¹⁸O_ice of −37‰; Table 1) yields sea level above the ice-free line in the Eocene and higher amplitudes than backstripped estimates. The 0.13‰/10 m calibration (with its implied δ¹⁸O_ice of −47‰; Table 1) yields sea levels that do not greatly exceed the ice-free line (71 m; with smoothed estimates uniformly less than 72 m) and agree with backstripped estimates.

Peak Early Eocene Warmth, Cooling, and CO₂

The cause(s) of post-Early Eocene cooling and BSL fall is not well understood but, in part, must be related to decreasing atmospheric CO₂ (Figure 5; Foster et al., 2012; Rae et al., 2021). Eocene global climate warmth was associated with high atmospheric CO₂, with values estimated to be between about 800 and 1,000 ppm (Foster et al., 2012), though our revised analysis suggests an Early Eocene peak of 1,000–1,500 ppm (Figure 5). Here, we revise the Bayesian hierarchical model of Miller et al. (2020) that was used to estimate the atmospheric concentration of CO₂ through time and its uncertainty (Figures 1, 5). We use proxy data compiled by Foster et al. (2017) and measurements from Antarctic ice cores (Bereiter et al., 2015), but supplement the dataset with the recent MECO pCO₂ estimates modeled by Henehan et al. (2020). The statistical model constructed by Miller et al. (2020) provided an estimate of atmospheric CO₂ concentrations through time based on the assumption that all proxy data within the Foster et al. (2017) compilation provided an equally accurate representation of atmospheric paleo-CO₂ concentrations. In our revision, we treat the ice-core measurements and boron isotope proxy as representative of the true value of pCO₂ and allow for the model to account for temporally varying offsets to the other proxies that improve the overall fit of the model to all the data. Therefore, we reproject the variations recorded by the variety of proxies accounting for systematic, proxy-specific offsets to probabilistically estimate the value of pCO₂ through geologic time as if it were being measured by ice-core measurements or the boron proxy (Figure 5).

A recently published community synthesis of atmospheric CO₂ (The Cenozoic CO₂ Proxy Integration Project Consortium, 2023) compares well with ours, with similar patterns of change, maximum values exceeding 1,500 ppm at 50 Ma (EECO), subordinate peaks at ca. 40 Ma (MECO), 35 Ma (LECO), and 17 Ma (MCO), and decreasing values in the Early Oligocene. We suggest our Bayesian statistical method would be particularly appropriate to apply to their revised data and future datasets.

We reaffirm the previously established first-order relationship between sea level and CO₂ (Figure 5). BSL decreases by 1.0 ± 0.26 m/Myr through the Cenozoic, whereas CO₂ decreases at a rate of 15 ± 0.5 ppm/Myr. The high global BSL values in the Early Eocene are due to largely ice-free conditions associated with peak CO₂ values, causing sea-levels ∼70 m higher than present. The high values in the Early Eocene GMGSL record are also affected by the peak CO₂ values (1,000–1,500 ppm), with a contribution of about 9 m for thermal expansion, but the primary mechanism tying sea level and CO₂ is the growth and decay of continental ice sheets. The GMGSL also records ∼50 m of Cenozoic fall due to lesser ocean volume (Wright et al., 2020), that may or may not exist (considering the ±90 m error) or be directly tied to the CO₂ record.

The rate of production of ocean crust is the first-order control on the OBVSL record we used (Wright et al., 2020), and it may also be a factor that contributes to variation in atmospheric CO₂ (Berner et al., 1983). However, the record of OBVSL derived from crustal production estimates is highly uncertain (Wright et al., 2020) and does not agree well (Figure 9) with independent estimates (e.g., van der Meer et al., 2022). Some argue that the prior for variation in rate of ocean crust production should be that it does not change and that the data are not convincing enough to adjust that prior (Rowley, 2002), which would invalidate both the OBVSL contribution to our GMGSL curve and preclude a tectonic source for variation of CO₂. The relatively high values estimated in the Early Eocene CO₂ record (Figures 1, 7) also do not necessarily require high values of OBVSL contributing to high values of GMGSL, because the Cenozoic decrease in CO₂ can be attributed to a number of factors that can be reduced to a balance between drawdown of CO₂ (via mechanisms like silicate weathering and the deposition and sequestration of organic carbon) and flux of CO₂ (Berner et al., 1983). The variation in the rate of silicate weathering through the Cenozoic is likely the most important variable affecting the long term record of CO₂ (Kent and Muttoni, 2013) that results in long-term pCO₂ variations controlling the variation in sea-level rather than the opposite. Specifically, decreasing CO₂ values during the Middle Eocene to Oligocene contributed to cooling, glaciation, and BSL fall (Figure 5; Foster et al., 2012). The Oligocene and younger Icehouse world with continental-scale ice sheets was established at a threshold of 3–3.5 x (840–1,000 ppmv; Figure 5) pre-anthropogenic CO₂ levels (280 ppmv) as suggested by models (DeConto and Pollard, 2003). The subsequent, secular changes in sea level since the Early Eocene have been modulated at discrete, astronomically driven Milankovitch periods: 19/21 kyr precession, 41 and 1,200 kyr obliquity (tilt), and quasi-100, 405, and 2,400 kyr eccentricity that have resulted in glacial/interglacial cycles (Westerhold et al., 2020).

Conclusion

There is a consensus that global mean geocentric sea-level changes over the past 34 million years were controlled largely by changes in the size of continental ice sheets, though most studies assume an ice-free world prior to this time. We present an updated synthesis of global sea level changes that contradicts this assumption. We use two independent records from oxygen isotopes-trace element paleothermometry and continental margin sediments that require large (20–30 m) sea-level changes in the interval 48–34 million years ago (Middle to Late Eocene) that firmly indicate sea-level changes driven by growth and decay of significant ice sheets, in agreement with new evidence from the East Antarctica continental shelf (Gulick et al., 2017). More importantly, the interval from 56 to 48 million years ago (Early Eocene) was the warmest of the past 66 Myr, with very high (>1,000 ppm) CO₂, yet we show that this supposedly ice-free world was punctuated by four sea-level lowerings that suggest growth and decay of ice sheets driving 20–40 m BSL changes, challenging current views.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author. The data has been submitted to the NCEI site (https://www.ncei.noaa.gov/metadata/geoportal/#searchPanel).

Author Contributions

KM: lead author in charge of oversite of all aspects of this study. WS: lead for all statistical analyses. JB: contributed to design of project and assembly of all data sets and produced most figures. YR: oversaw Mg/Ca compilation. AH: generated Mg/Ca data. RK: contributed to evaluation of sea level record. JW: conceptual design of project and interpretation of oxygen isotopic record. All authors contributed to the article and approved the submitted version.

Funding

Supported by NSF grants OCE16-57013 (KM) and OCE14-63759 (KM and JB). WS was supported by a grant from DOE DE-FOA-002000 through Battelle.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

This paper was presented in hybrid form at the Sea level change in the past, present, and future conference held at the Geological Society of London, 06–07 February 2023. We thank the organizers of that meeting D. Lunt, E. Gasson, DvdM, and N. Barlow for the invitation to present and discuss. Reviews from JA, DvdM, and SG were particularly helpful. We thank the Office of Advanced Research Computing (OARC) at Rutgers University for providing access to the Amarel cluster and associated research computing resources. We thank D. Pak and M. Katz for Site 577 samples.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.escubed.org/articles/10.3389/esss.2023.10091/full#supplementary-material

Footnotes

¹https://sealevel.colorado.edu

²https://www.ncei.noaa.gov/metadata/geoportal/#searchPanel

References