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

We estimate ice-volume driven (barystatic; BSL) sea-level changes for the Cenozoic using new Mg/Ca data from 58 to 48 Ma and a revised analysis of Mg/Ca trends over the past 66 Myr. We combine records of BSL, temperature-driven sea level, and long-term ocean basin volume variations to derive a new global mean geocentric sea level (GMGSL; “ eustatic ” ) estimate. Bayesian analysis with Gaussian process priors shows that our BSL estimate shares a component that covaries on the Myr scale with “ backstripped ” relative sea-level (RSL) estimates (accounting for compaction, loading, and thermal subsidence) from the US Mid-Atlantic Coastal Plain, validating our method and estimates with errors of ±10 m. Peak warmth, elevated GMGSL and BSL, high CO 2 , and ice-free conditions occurred at times in the Paleocene to Eocene (ca. 64, 57.5, 35 Ma) and in much of the Early Eocene (55 – 48 Ma). However, our new results show that the Early Eocene was punctuated at speci ﬁ c times by several Myr-scale sea level lowerings (~20 – 40 m) that require growth and decay of signi ﬁ cant continental ice sheets even in the supposedly “ ice-free ” world. Continental-scale ice sheets waxed and waned beginning ca. 34 Ma ( > 50 m BSL changes), with near complete collapse during the Miocene Climate Optimum (17 – 14.8 Ma). Both the BSL and RSL estimates have markedly higher Oligocene to Early Miocene Myr-scale amplitudes (20 – 60 m) than recently published δ 18 O-based estimates ( < 20 m) and much lower estimates than those of Exxon Production Research ( > 100 m), leading us to reject those estimates. The US Mid-Atlantic margin RSL was dominated by GMGSL but was overprinted by changes in mantle dynamic topography on the several Myr scale, showing approximately 50 m higher Eocene estimates and regionally propagating Miocene RSL changes.


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).Sealevel and ice volume varied during the Oligocene to Early Miocene, with major deglaciation during the Miocene Climate Optimum (MCO; 17.0-14.9Ma; 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  (Miller et al., 2020; this study) with an 800 kyr LOESS smoothed fit (solid black line); Mi1, Oi2a, Oi2, Oi1 are benthic foraminiferal δ 18 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 2 (this study) from alkenone (yellow dots), boron (purple dots), paleosols (green dots), stomatal indices (blue dots).The best estimate of Cenozoic pCO 2 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.
(EAIS) in the Oligocene to Middle Miocene and the latter with a largely stable EAIS during the last 13.8Ma.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).
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 δ 18 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 δ 18 O increases in deep sea benthic foraminifera.Miller et al. (2020) used deep Pacific benthic foraminiferal δ 18 O and Mg/Ca records to calculate seawater δ 18 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 δ 18 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; 1 ) 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.

Time Scale and Benthic Foraminiferal δ 18 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 (δ 18 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) δ 18 O b records exclusively from Pacific sites (Figures 1-3).All δ 18 O b records have been astronomically tuned to the GTS by the publications cited in the synthesis.Our δ 18 O b record is similar to δ 18 O b splices of De Vleeschouwer et al. (2017) and Westerhold et al. (2020) who used a similar approach to assemble a synthesis of δ 18 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.
Our strategy and that of Miller et al. (2020) differs from recent δ 18 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, δ 18 O, and δ 13 C reservoirs and up to 80% of the Paleogene (Miller and Fairbanks, 1985).We view using sites from PDW as a prerequisite to using δ 18 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 δ 18 O sw term.Because they have a larger temperature component, Atlantic δ 18 O records amplify orbital variations in δ 18 O b that were the target of previous studies and are thus wellsuited for these studies (De Vleeschouwer et al., 2017;  The BSL R1 (t) record (light blue line) is sea level computed using a 0.13/ 10 m δ 18 O sw to sea level calibration, assuming that 80% of the variability in δ 18 O sw calculated using the paleotemperature equation (Eq.1), δ 18 O b , and pT R0 (t) is due to seawater variation in δ 18 O.The black line is an 800 kyr LOESS regression of the BSL R1 (t) 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 R1 (t).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 R2 (t)), effectively reducing the variation in δ 18 O sw when the long-term estimate of continental ice volume (BSL R1,LT (t); Eq. 8) is low.The black line is a LOESS regression of the pT R2 (t) record (red line) with an 800 kyr search window.Panel (E) The BSL R2 (t) record, light blue line, is the sea-level variation that corresponds with the pT R2 (t) paleotemperature record in Panel (D) The black line is an 800 kyr LOESS regression of the BSL R2 (t) 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 R1 (t).The red vertical line indicates icefree conditions at 71 m.Pleist., Pleistocene; Plio., Pliocene.Westerhold et al., 2020).However, by not using Atlantic data, our δ 18 O b compilation is not as complete as previous splices in some intervals, though our target was primarily the Myrscale 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.6Ma (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 δ 18 O splice (δ 18 O b ) of Miller et al. (2020); 2) bottom water paleotemperatures, pT, derived from Mg/Ca at the time, t, of each δ 18 O b measurement, i, derived using an iterative process that leverages Mg/Ca measurements of Pacific benthic foraminifera and the δ 18 O b record (Figure 3); 3) the paleotemperature equation for benthic foraminifera that produces an estimate in the variation in the δ 18 O values of seawater (δ 18 O sw ; Eq. 1) given pT and δ 18 O b of Cibicidoides ((Lynch-Stieglitz et al., 1999); and 4) a calibration factor for the conversion of δ 18 O sw change to sea-level change, taken as 0.13‰/10 m (Winnick and Caves, 2015).This calibration will vary with the δ 18 O composition of ice, a factor we address in the Discussion.

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 (t j ), was then generated using the corrected Mg/Ca ratios, Mg/Ca corr (t j ), 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 (t j ) 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 (t j ) generated using Mg/Ca bf (t j ) 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.

Mg Ca corr t j
Mg Ca bf t j − 0.009

BSL Calculation
The initial estimate of bottom water paleotemperature derived using the LOESS regression of the pT Mg/Ca (t j ) estimates, i.e., pT R 0 (t), was used within the paleotemperature equation (Eq.3), interpolating values at times t i to derive an initial estimate of δ 18 O sw (t i ).This initial estimate of δ 18 O sw (t i ) was used to create an initial sea-level curve, BSL R 0 (t i ), applying a δ 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 δ 18 O sw (t i ) 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 (t j ) 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 δ 18 O sw (t i ) by 20% across the entire curve and account for this residual through a revision of the initial paleotemperature record, pT R 0 (t), to generate a paleotemperature curve that is supplemented by the record of δ 18 O b variations (pT R1 (t); Eq. 6).
BSL R 0 t i ( ) 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 δ 18 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 δ 18 O b variation to paleotemperature is accommodated by 1) applying a LOESS algorithm with a 2 Myr window to the BSL R 1 (t i ) sea-level curve (Eq.8) to produce BSL R 1,LT (t i ); 2) converting the fraction of the residuals of the BSL R 1 (t i ) (Figure 3C) relative to BSL R 1,LT (t i ) defined in Eq. 9 to variation in paleotemperature using Eq. 10 to produce a second revision of the paleotemperature record, pT R 2 (t).(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 (t i ), as defined in Eq. 9, and scales linearly from 0% at BSL R 1,LT (t i ) values below −50 m to 85% at BSL R 1,LT (t i ) values of 50 m, and again scales linearly from 85% at BSL R 1,LT (t i ) values of 50 m-100% at BSL R 1,LT (t i ) of 75 m.The variation in sea level is then recalculated (Eq.11) producing BSL R 2 (t i ), 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 NCEI2 .

Thermal Expansion and Tectonic Corrections
Here we calculate Cenozoic thermosteric sea level variation, θ, using the pT R 2 (t) 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 θ, 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 θ using the regressions generated by Hieronymus (2019).Hieronymus (2019) suggested either 0.73 or 1.18 m of θ 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 θ a (t) or θ b (t), representative of the 1.18 m/ °C or 0.73 m/ °C regression constant for the conversion of surface temperature to θ, respectively.Alternatively, θ 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 sealevel 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, OBVSL(t), was , and ice core measurements (gray dots, only available for time 800 ka to present).The best estimate of Cenozoic pCO 2 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 2 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 2 .The proxy-specific offsets were determined using a Gaussian process regression model and they are shown in Supplementary Figure S3.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 OBVSL(t) was added to the records of θ a (t i ) and BSL R 2 (t i ) to generate a record of GMSL(t i ) (Eq. 15; Figure 5).

GMGSL t
Modeling Cenozoic CO 2 Variations Early Eocene global warmth and high sea level are generally ascribed to high atmospheric CO 2 levels (e.g., Anagnostou et al., 2016) and we explore these relationships with an updated analysis of CO 2 records.We revise the Bayesian hierarchical model with Gaussian process priors (Rasmussen and Williams, 2006;Ashe et al., 2019)  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 2 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, proxyspecific offsets to probabilistically estimate the value of pCO 2 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 2 measurement proxies, the CO 2 values, y i (Eq.16), are modeled as the sum of the modeled CO 2 concentration, f(t i ); the proxyspecific offset, g k (t i ); and the measurement value errors, ε y i .We do not consider the temporal error associated with individual data points.
At the process level, the best estimate for the concentration of pCO 2 in the atmosphere through time is represented by a composite function, f(t) (Eq.17), comprising four subcomponents.This function, f(t), provides an estimate of pCO 2 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 2 represented by individual proxies.These proxy-induced differences and other noise are modeled by composite function g k (t) (Eq.18).
The component processes are each modeled with a Gaussian process prior.Three of the component functions that comprise f(t) (m(t), l 0 t, and c 0 ) are representative of a signal shared by all proxy data, and the final component (Δ ice (t)) represents the precise fluctuations recorded by icecore measurements over the past 800 kyr.This shared component, m(t), is a non-linear process that represents long-term variations in pCO 2 .Additionally, a linear process, l 0 t; and an offset, c 0 , represents the long-term trend in pCO 2 as it is represented by the entirety of the pCO 2 data compilation, y i .A second non-linear process, Δ ice (t), captures the highfrequency changes that are only captured by the ice-core record, which only spans the past 800 kyr.
The function g k (t) accounts for the systematic, linear proxyspecific differences between the less reliable alkenone, stomata, paleosol, and liverwort proxies and the paleo-pCO 2 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(t) and the boron proxy or ice-core measurements, the components of f(t), contain values of pCO 2 through geological time that closely represent those boron proxy and ice-core measurements that are more reliable records of paleo pCO 2 .Functions f(t), g k (t), 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 2 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 (t i ) 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 (t i ) and backstripping studies of the NJ continental margin were modeled over time, t, as function g(t) (Eq.21) and h(t) (Eq.22), respectively.These two modeled functions share both linear and non-linear components, m(t) and l(t), 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 (ε a i and ε b j ).Miller et al. ( 2020) previously calculated those errors to be normally distributed with a standard deviation of ±13 m for the BSL R 2 (t i ) estimate and ±15 m for the backstripped NJ record.The component processes within g(t) and h(t) 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(t) component, shared between both records, signifies longterm ice-volume changes and the trend of variation in BSL.In contrast, l h t is unique to the backstripped record, h(t), 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(t) term shared by both g(t) and h(t) captures BSL variation attributable to changes in continental ice volume.The other non-linear term, Δ(t), is unique to h(t) and represents the differences between the sea-level changes recorded by BSL R 2 (t i ) (i.e., g(t)) and those recorded via backstripping (i.e., h(t)) 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 timeaverage 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(t) and h(t) 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(t) and h(t) and their constituent processes (see Discussion).The posterior estimates of g(t) and h(t) facilitated a robust comparison of the BSL R 2 (t i ) estimate and the backstripped NJ margin sea-level record, shown as h(t) − g(t).
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).
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 3 of grounded ice sheets above sea level (above floatation) in Antarctica and of 2.99 ± 0.02 million km 3 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 sealevel equivalent does not include ice grounded below sea level (below floatation) of 6.6 million km 3 (Morlighem et al., 2020) that would have raised sea level once isostatic compensation is realized and affected δ 18 O sw .The 6.6 million km 3 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 δ 18 O-Mg/Ca based BSL estimates (Figure 4).The BSL sea-level record of Miller et al. (2020) and our uncorrected BSL 1 (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 2 (Figure 3E) accounts for high-frequency variations of bottom water paleotemperatures at higher paleotemperatures and sealevel 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 2 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).

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.
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 onshoreoffshore 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., g(t) and h(t), respectively) represent a combination of component processes that have been fitted using Gaussian process priors.A shared non-linear component, m(t), 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(t) and l h (t), 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 δ 18 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 δ 18 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  Moucha et al., 2008;Kominz et al., 2016;Schmelz et al., 2021).
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 longterm 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 δ 18 O b Based BSL Estimates Rohling et al. (2021) and Rohling et al. (2022) provide a synthesis of sea-level and temperature variations using δ 18 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 δ 18 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 δ 18 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 δ 18 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 δ 18 O b .The primary assumption underlying the validity of the curve is that the quadratic equation correctly identifies the ice-free value of δ 18 O b .
We show that the modeled sea-level change (Rohling et al., 2021;2022) is similar to our BSL R1 , 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 δ 18 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 δ 18 O b data to sea-level values.The quadratic relationship applied by Rohling et al. (2022) implies that there is a greater proportion of temperaturedriven fractionation effect at higher values of δ 18 O b relative to δ 18 O sw change.As discussed above, this is rational because higher values of δ 18 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.
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)  ).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 δ 18 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 δ 18 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 δ 18 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 δ 18 O sw and the δ 18 O for ice sheets (δ 18 O ice ).The methods are quite different.Rohling et al. (2022)  The assumptions made in scaling δ 18 O sw to sea level have implications for δ 18 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 δ 18 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 δ 18 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 δ 18 O b to generate a sea-level record, though their assumptions result in lower amplitude variability in sea level relative to given variations in δ 18 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 δ 18 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 δ 18 O sw to sea-level calibration should be investigated further.However, the assumptions associated with Rohling et al. (2022)   Testing the Sensitivity of Sea Level to Mg/Ca Assumptions and δ 18 O sw Calibrations As noted throughout, scaling of δ 18 O b to sea-level requires numerous assumptions, particularly paleotemperature, to obtain δ 18 O sw and the relationship of δ 18 O sw to sea level.We regard the agreement between two independent estimates, backstripping and δ 18 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 δ 18 O b -Mg/Ca method and differences due to assumed sea level to δ 18 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‰ δ 18 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 δ 18 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 δ 18 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 δ 18 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 2
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 2 (Figure 5; Foster et al., 2012;Rae et al., 2021).Eocene global climate warmth was associated with high atmospheric CO 2 , 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 2 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 2 estimates modeled by Henehan et al. (2020).The statistical model constructed by Miller et al. (2020) provided an estimate of atmospheric CO 2 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 2 concentrations.In our revision, we treat the ice-core measurements and boron isotope proxy as representative of the true value of pCO 2 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 2 through geologic time as if it were being

FIGURE 1 |
FIGURE 1 | Panel (A) Eocene to Oligocene Pacific benthic foraminiferal (Cibicidoides spp.) δ 18 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 δ 18 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 byMiller et al. (2008) and updated withGulick et al., 2017;Barr et al., 2022); Panel (C): compilation of statistically derived estimates of Cenozoic pCO 2 (this study) from alkenone (yellow dots), boron (purple dots), paleosols (green dots), stomatal indices (blue dots).The best estimate of Cenozoic pCO 2 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 FigureS3.

FIGURE 2 |
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 δ 18 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 |
FIGURE 3 | Panel (A) δ 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 δ 18 O sw calculated is due to temperature variation The black line is a LOESS regression of the pT R1 (t) record (red line) with an 800 kyr search window.Panel (C)The BSL R1 (t) record (light blue line) is sea level computed using a 0.13/ 10 m δ 18 O sw to sea level calibration, assuming that 80% of the variability in δ 18 O sw calculated using the paleotemperature equation (Eq.1), δ 18 O b , and pT R0 (t) is due to seawater variation in δ 18 O.The black line is an 800 kyr LOESS regression of the BSL R1 (t) record that can be compared with an 800 kyr LOESS regression of theMiller et al. (2020) BSL estimate (blue line) that was calculated identically to BSL R1 (t).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 R2 (t)), effectively reducing the variation in δ 18 O sw when the long-term estimate of continental ice volume (BSL R1,LT (t); Eq. 8) is low.The black line is a LOESS regression of the pT R2 (t) record (red line) with an 800 kyr search window.Panel (E) The BSL R2 (t) record, light blue line, is the sea-level variation that corresponds with the pT R2 (t) paleotemperature record in Panel (D) The black line is an 800 kyr LOESS regression of the BSL R2 (t) record that can be compared with an 800 kyr LOESS regression of theMiller et al. (2020) BSL estimate (blue line) that was calculated identically to BSL R1 (t).The red vertical line indicates icefree conditions at 71 m.Pleist., Pleistocene; Plio., Pliocene.

FIGURE 4 |
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.

FIGURE 5 |
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 2 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 2 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 2 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 2 .The proxy-specific offsets were determined using a Gaussian process regression model and they are shown in Supplementary FigureS3.

FIGURE 6 |
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.

FIGURE 7 |
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 (t), and the blue line shows the long-term trend of pT Mg/Ca (t) 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 R2 (t), which is an estimate of the high-resolution temperature record derived using the initial temperature estimate, the δ 18 O b splice, and a modulation of variation in δ 18 O sw when the long-term estimate of continental ice volume (BSL R1,LT (t); 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 (t), 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 R2 (t) record, light blue line, is the sea-level variation that corresponds with the pT R2 (t) paleotemperature record shown in the left panel.The black line is an 800 kyr LOESS regression of that BSL R2 (t) record.The black line corresponds to an 800 kyr LOESS regression of BSL R2 (t) 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 R2 (t) 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 R2 (t) with a δ 18 O 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 R2 (t) to end member calibrations.The red vertical line indicates ice-free conditions at 71 m.

FIGURE 8 |
FIGURE 8 | Modeled estimates of sea-level change, recorded by the δ 18 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., g(t) and h(t), respectively) represents a combination of component processes that have been fitted using Gaussian process priors.A shared non-linear component, m(t), 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(t) and l h (t), 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).
to δ 18 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 δ 18 O b values in Lisiecki and Raymo 2005 [~3‰-5‰] calculated δ 18 O sw by budgeting sea-level change through a model of ice sheet growth and decay.This model generates a variable relationship between δ 18 O sw and sea level because it incorporates functions that model the increasingly lower δ 18 O ice values of continental ice sheets as they grow.As a result, δ 18 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 δ 18 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 δ 18 O b values of −0.5‰ and sealevels of 71 m.
conversion of δ 18 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 ).

FIGURE 9 |
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.

FIGURE 10 |
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.
Miller et al. (2020)15)020)to estimate atmospheric CO 2 concentration over time.To model pCO 2 herein, we use the combination of paleosol and stomata proxy data compiled byFoster et al. (2017)and measurements from Antarctic ice cores(Bereiter et al., 2015)that were used byMiller et al. (2020), while replacing the boron and alkenone proxy data fromFoster et al. (2017)with the updated datasets made available byRae et al. (2021).The statistical model constructed byMiller et al. (2020)provided an estimate of atmospheric CO 2 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 2

TABLE 1 |
Relationship of δ 18 O sw -sealevel calibration to δ 18 O ice and sea-level changes for given δ 18 O sw variations assuming mean depth of ocean of 3,700 m.