Introduction

The early Eocene (56–48 Ma) was the warmest period of the Cenozoic, characterized by very high atmospheric CO2 concentrations (≥1000 ppm)1,2, reduced equator to pole temperature gradients, and absence of continental ice sheets3,4,5,6. As such, the early Eocene hothouse is considered a potential analogue for our future climate7. Accurate temperature reconstructions are essential to obtain a better understanding of the functioning of the climate system under high CO2 levels. Temperature proxy data can be used to estimate global mean temperatures and to validate the performance of model simulations of the Eocene climate2,8,9,10. In comparison to the sea surface, the deep ocean is spatially uniform in temperature and is considered a relatively stable component in the climate system because of its high heat capacity11,12,13. For these reasons, deep-sea temperature is commonly used to infer the state of the global climate7,9,14,15,16.

Over the past three decades, much of our knowledge of the early Eocene climate on both long (>105 yr) and short (<105 yr) time scales has been based on deep-sea temperature estimates derived from oxygen isotope measurements of microfossil carbonate skeletons of benthic foraminifera (δ18Ob). For instance, negative excursions in high-resolution stable oxygen and carbon isotope records from deep marine sediments have revealed the periodic occurrence of multiple transient (10–100 kyr) episodes of global warming and ocean acidification (hyperthermal events; e.g. PETM, ETM2, and ETM3) (refs. 12,16,17,18,19,20,21), generally linked to massive release of isotopically light carbon into the ocean-atmosphere system4,18,22,23.

However, the reliability of using δ18Ob for temperature reconstructions is hampered by uncertainties in non-thermal factors. Besides temperature, this proxy reflects the δ18O composition of the seawater (δ18Osw), which is poorly constrained for the Earth’s past24,25. The traditional view is that with the absence of ice sheets in the Eocene, a fixed ice-free δ18Osw value can be applied to convert the δ18Ob signal directly into deep-sea temperature12,16,26. Yet, in this approach, other factors, such as salinity of deep water masses, are generally ignored as potential sources that could also have caused changes in δ18Osw (refs. 24,25). Furthermore, it is found that foraminiferal oxygen isotopes are sensitive to seawater pH and species-specific physiological factors (vital effects)24,25,27,28,29,30.

In an attempt to overcome these uncertainties, estimates of early Eocene deep-sea temperatures have additionally been obtained from benthic foraminiferal magnesium/calcium (Mg/Ca) ratios11,31,32. However, while this temperature proxy has the advantage of not being influenced by δ18Osw, it depends on the Mg/Ca composition of the seawater (Mg/Casw) (refs. 32,33). Although Mg and Ca have a conservative behaviour and are therefore equally distributed in the ocean, this Mg/Casw ratio appears to have varied on long time scales (>1 Myr) in the past34,35. Uncertainties in Mg/Casw reconstructions and lack of sufficient understanding of the effect of these changes on foraminiferal Mg/Ca can bias temperature estimates32,33,36,37,38,39. In addition, the foraminiferal incorporation of Mg and Ca is affected by the carbonate chemistry of the ocean and by vital effects that require species-specific calibrations36,40,41,42,43,44. The dependency on other factors complicates the use of oxygen isotopes and Mg/Ca in reconstructing ocean temperature. Available deep-sea temperature reconstructions are not as robust as commonly appreciated and estimates from new proxies are hence desired.

The carbonate clumped isotope thermometer (Δ47) shows potential to improve temperature reconstructions as it is insensitive to the isotope composition of the seawater45,46 and is largely unaffected by pH47,48,49,50 and physiological factors51,52,53. This proxy uses the temperature dependence of the degree of bonding of two heavy rare isotopes (13C and 18O) within the carbonate ion, with more clumping favoured at cooler temperatures45,46. Recent analytical developments now ensure inter-laboratory consistency of measurements54,55 and allow analysis of relatively small sample sizes (e.g. foraminifera)56,57,58. These advances in the application of the clumped isotope proxy have enabled reconstruction of marine temperatures independent from non-thermal uncertainties59,60,61,62,63. So far, Δ47-based reconstructions of the relatively cooler Miocene and middle Eocene have revealed a warmer deep ocean than traditionally accepted61,62,63,64, raising the question of whether other periods in the Cenozoic also experienced warmer deep-sea temperatures than assumed.

Here we report Δ47-based deep-sea temperature estimates across two early Eocene hyperthermal events, namely Eocene Thermal Maximum 2 (ETM2 or H1) and H2 that occurred ~2 Myr (~54 Ma) after the Paleocene-Eocene Thermal Maximum (PETM) (refs. 65,66). High-resolution (~1 kyr) benthic foraminiferal carbon (δ13Cb) and oxygen isotope records were previously generated from four ODP Sites (1262, 1263, 1265, and 1267) at the Walvis Ridge, Southeastern Atlantic and from one ODP Site (690) at Maud Rise in the Weddell Sea, Southern Ocean66 (Supplementary Fig. S1) with paleowater depths ranging from ~1500 to ~3600 m (refs. 17,67,68,69). The δ18Ob records of all these sites show almost identical values, suggesting that they reflect similar deep-sea conditions66. We re-examined these sites and applied paired stable and clumped isotope analyses on the benthic foraminiferal species Nuttallides truempyi and Oridorsalis umbonatus for accurate reconstruction of deep-sea temperatures in the South Atlantic across these two hyperthermals. Overall, we find warmer deep-sea temperatures from Δ47 than estimates based on δ18Ob, which has implications for the conceptual understanding of the δ18Osw composition and/or pH effects that govern the oxygen isotope proxy. Based on our findings, we propose a re-evaluation of the interpretation of δ18Ob records in the geological past.

Results and discussion

Stable and clumped isotope analysis

Our carbon and oxygen isotope measurements (Fig. 1) are consistent with the previously published δ13Cb and δ18Ob records66. Distinct negative excursions in δ13Cb (−1.4‰ and −0.8‰, respectively) and δ18Ob (−0.8‰ and −0.5‰, respectively) characterize ETM2 and H2, showing that ETM2 is the larger event of the two. These negative δ13Cb and δ18Ob excursions suggest the injection of large amounts of 13C-depleted carbon into the climate system and elevated deep-sea temperatures, respectively65,66. Limited isotope data are obtained for the peak ETM2 interval at the Walvis Ridge sites, as benthic foraminifera are rare and small-sized in the associated red clay layer, known as the Elmo horizon65,66.

Fig. 1: South Atlantic benthic foraminiferal δ13C and δ18O records across ETM2 and H2.
figure 1

The figure shows the δ13C and δ18O measurements of benthic foraminifera N. truempyi (a and b, respectively) and O. umbonatus (c and d, respectively). The offsets in the stable isotopes between the two species suggest species-specific vital effects30. The sites from the Walvis Ridge used for benthic foraminiferal stable and clumped isotope analysis include ODP Sites 1262, 1263, 1265, and 1267. Paleowater depth of these sites ranges from ~1500 to ~3600 m in the order 1263, 1265, 1267, 1262 (shallow to deep)17,67. Maud Rise is represented by ODP Site 690 with a ~2100 m paleowater depth68,69. The stable isotope data are plotted against age following the age models of Stap et al.66 and Westerhold et al.126.

Full size image

Due to the sporadic natural abundance of 13C–18O bonds within carbonate ions45,46, clumped isotope analysis requires large samples to achieve high analytical precision54,55,56,57,58, which is needed in paleoclimate reconstructions. Averaging of numerous replicate Δ47 measurements (>30 replicates of ~100 µg each) is required to mitigate large analytical uncertainty54,55,56,57,58. The strategy of averaging Δ47 measurements for precise reconstruction of small temperature variability becomes especially important when sample size is limited70. In our study, the abundance of N. truempyi and O. umbonatus is too low in the samples at all sites to take sufficient measurements for the precise reconstruction of deep-sea temperature change across ETM2 and H2. As a solution to acquire the desired precision (<3 °C at 95% Confidence Interval), we sorted the Δ47 measurements from all sites based on their corresponding δ18Ob values and compiled three average Δ47-based temperatures that represent the average background state (high δ18Ob), average slope (middle δ18Ob), and average hyperthermal peak (low δ18Ob) conditions across the entire studied period70. Thus instead of showing deep-sea temperatures in a time series, we plotted three temperature bins that contain measurements of both the ETM2 and H2 intervals, which are representative of three different climate states. Unequal variance (Welch’s) t-tests were performed between the ‘high δ18Ob’ and ‘low δ18Ob’ Δ47 measurements to determine the most optimal bin size for the background and hyperthermal temperatures70 (Supplementary Table S1 and Fig. S2).

Deep-sea temperatures

Our Δ47-based reconstruction indicates very warm deep waters with average background temperatures of 13.5 ± 1.9 °C (95% Confidence Interval) and average hyperthermal peak temperatures of 16.9 ± 2.3 °C (95% CI) (Fig. 2b; Supplementary Fig. S3 for temperatures per site). Statistically, these background and average hyperthermal peak temperature bins are significantly different from each other at 98% CI (p-value of 0.02). The peak temperature of ETM2 is expected to be higher than our reconstructed average of 16.9 ± 2.3 °C (95% CI). This is because the bin of the hyperthermal peak temperature combines Δ47 measurements of both the ETM2 and H2 events, in which ETM2 has almost twice the magnitude in δ18Ob of H2 (ref. 66) (Figs. 1 and 2a). As such, this approach creates a smoothing effect in which the peak temperature of ETM2 is underestimated.

Fig. 2: Deep-sea temperatures and δ18O seawater composition across ETM2–H2.
figure 2

a δ18Ob values derived from N. truempyi and O. umbonatus that are corrected towards Cibicidoides (assumed seawater equilibrium). b Deep-sea temperatures based on Δ47 and δ18Ob including the number of Δ47 (and δ18Ob) measurements that are contained in the background, slope, and hyperthermal peak temperature bins. The yellow (background), red (slope), and purple (hyperthermal peak) bars indicate the range in δ18Ob values corresponding to the Δ47 measurements that were used for compiling the three bins. The dark grey error bars on the Δ47-based temperatures represent fully propagated analytical and calibration uncertainties as 68% (solid) and 95% (dashed) CI (Methods). The dashed red envelope on the δ18Ob-based deep-sea temperatures displays the uncertainty (95% CI) of the δ18Ob-temperature calibration of Marchitto et al.77. Analytical uncertainty on the δ18Ob-based temperatures of the three bins is very small (±0.1 °C 95% CI) and is therefore contained within the data symbols (Methods). c δ18Osw (in ‰ VSMOW) of the background, slope, and hyperthermal peak bins including uncertainties as 68% (solid) and 95% (dashed) CI. Values of δ18Osw were calculated using the mean Δ47-based temperatures in combination with the average δ18Ob values of the three bins in the δ18Ob-temperature relationship of Marchitto et al.77 (Methods). The horizontal red and blue dashed lines show the assumed δ18Osw value based on ice-free conditions and the present-day δ18Osw value, respectively.

Full size image

In terms of absolute values, our deep-sea temperatures from clumped isotopes are in broad agreement with benthic Mg/Ca-based records available from the Pacific, which indicate mean temperatures at 54 Ma of 15.0 or 16.3 °C with an uncertainty of ±3.0 °C (90% CI) based on two calibrations32,39,40,43 (Supplementary Fig. S4). However, Mg/Ca records from the Paleocene to the early Eocene provide questionable trends in deep-sea temperature, as the Paleocene stands out to be warmer than the early Eocene, which does not become apparent from δ18Ob records12,16 (Supplementary Fig. S4). These unexpected temperatures suggest that these records may be biased before the middle Eocene (>48 Ma) due to changing Mg/Casw and different foraminiferal Mg/Ca-temperature sensitivity32,39, which makes them less reliable.

To compare the Δ47-based deep-sea temperatures with those derived from δ18Ob estimates, we used the same approach as applied in the revised Cenozoic δ18Ob-based deep-sea temperature compilation by Cramer et al.32. Traditionally, it is common practice12 to apply the classic δ18O-temperature relationship of Shackleton71 based on few core top data, together with the assumption that the δ18O composition of modern Oridorsalis is in equilibrium with the seawater composition72,73,74. Cibicidoides is the foraminiferal genus mostly used for deep ocean reconstructions12 and is usually corrected for seawater equilibrium since it was found that the modern genus has a 0.64‰ offset in δ18O with Oridorsalis72,73,74. However, more recent studies show evidence that modern Cibicidoides are actually in close isotopic equilibrium with seawater and derive a δ18O-temperature relationship from this genus based on a large dataset75,76,77. To calculate early Cenozoic deep-sea temperatures from other foraminiferal species with reference to Cibicidoides, well-documented correction factors are used to account for interspecies δ18Ob offsets30,32. We adjusted the N. truempyi and O. umbonatus δ18O values towards Cibicidoides and subsequently calculated deep-sea temperatures from these values following Marchitto et al.77 (Methods and Fig. 2b). For the mean δ18Osw of the global ocean used in δ18Ob-based temperature reconstructions, we adopted a fixed value of −1‰ VSMOW based on the assumption of absence of ice sheets in the early Eocene26.

In comparison to the δ18Ob-based deep-sea temperatures of the three bins, our reconstructions from clumped isotopes are on average three degrees warmer, i.e. 13.5 ± 1.9 °C versus 10.8 ± 0.1 °C for the background, 15.0 ± 2.0 °C versus 12.0 ± 0.1 °C for the slope, and 16.9 ± 2.3 °C versus 13.6 ± 0.2 °C (all errors as 95% CI) for the average hyperthermal peak temperatures (Fig. 2b). Note that if we would instead apply the traditional δ18O-temperature equation71 and assume Oridorsalis δ18O representing equilibrium, the derived temperature differences between clumped isotopes and oxygen isotopes are larger as the δ18Ob-based temperatures become ~1 °C cooler using that approach (Supplementary Fig. S4). The significant temperature discrepancy between Δ47 and δ18Ob means that there is an apparent bias in the deep-sea temperatures from oxygen isotopes towards cooler values. This underestimation of temperature by δ18Ob suggests that the used assumptions on non-thermal factors are incorrect.

Constraints on early Eocene δ18Osw composition

An important consequence of the warmer Δ47-derived deep-sea temperatures is that they imply a mean value for deep water δ18Osw of approximately −0.35 ± 0.4‰ VSMOW (95% CI), which is significantly higher than the usually assumed value (Fig. 2c). At first glance, this enriched δ18Osw value would suggest the presence of substantial continental ice sheets in the early Eocene. However, this hypothesis is unlikely given that subtropical temperatures and vegetation prevailed at both poles3,5,78,79,80. Although small ephemeral ice caps might have existed on Antarctica and in the Arctic region, records of ice-rafted debris in marine sediments indicate that polar ice formation was not initiated before the middle Eocene81,82. Hence, we turn to other factors that may explain a relatively high δ18Osw composition.

The distribution of different water masses in the early Eocene could have caused spatial differences in the δ18Osw value of the ocean. The early Eocene deep Atlantic basin might have been filled with more saline waters, characterized by high δ18Osw. Assigning exact numbers on the salinity is complicated by the fact that the salinity–δ18Osw relationship is not constant across the global ocean, as observed in modern water masses76. Deep water formation is thought to have predominantly occurred in the high-latitude Southern Ocean83,84. Although, saline deep waters might have formed more easily at the lower latitudes, such as the Tethys Ocean, where high evaporation would have increased salinity and hence density of the surface water85,86. However, a direct implication of this situation would be that other regions of the deep ocean must have had much fresher waters (with lower δ18Osw). Large regional salinity differences are required to meet the global seawater isotope budget and thus to balance the global ocean δ18Osw towards ice-free conditions. At first sight, such large heterogeneity in the ocean, i.e. large spatial differences in δ18Osw, is difficult to reconcile with relatively small gradients in δ13Cb and δ18Ob that have been found between the ocean basins in the early Eocene20,21,87. However, we acknowledge that similar δ18Ob values across different ocean basins do not necessarily exclude the possibility that water masses would have had different properties in terms of temperature and salinity. Although speculative, it might be that warmer saltier water masses and cooler fresher water masses resulted in the same δ18Ob values being recorded by the benthic foraminifera. In this situation, an intensified hydrological cycle in the Eocene hothouse may have caused increased regional salinity differences in the ocean78,80,88,89,90, which might have enabled the formation of deep water masses with different properties.

Changes in substantial groundwater storage on land and isotope interaction between seawater and oceanic crust are factors that may have affected the δ18Osw composition on a global scale. Much like ice sheets, large aquifers filled with isotopically depleted water would change the δ18Osw towards more positive values91. It has been hypothesized that the charging and discharging of these groundwater reservoirs could explain apparent large sea level fluctuations in the hothouse climate of the Cretaceous91. However, geological evidence for large aquifer existence is not available and climate models have not yet been able to simulate increased groundwater storage in warm climates92. On million-year time scales, isotopic exchange of ocean water cycling through the oceanic crust at mid-ocean ridges may change through faster and slower seafloor spreading93,94. The classic idea has been that seafloor spreading rates and associated volcanic activity were higher in the Cretaceous and early Eocene times, which is often used as one explanation for the high atmospheric CO2 levels4,95,96. Studies propose that an intensification of high-temperature basalt-seawater interactions (e.g. circulation of hydrothermal fluids) during faster seafloor spreading would have released more 18O from the ocean crust, enriching the seawater with heavy oxygen93,94. Afterwards, low-temperature processes (e.g. submarine weathering) would have become dominant at mid-ocean ridges when spreading rates decreased. This low-temperature basalt alteration would then have preferentially removed 18O from the seawater, decreasing the ice-free δ18Osw value towards more recent times93,94. However, there is no consensus if seafloor spreading rates were truly higher in these hothouse climates, as tectonic reconstructions and model simulations find conflicting results97,98,99. In summary, elevated δ18Osw in the early Eocene deep ocean can be ascribed to several factors. Our current understanding does not point towards one explanation as the most probable and a high δ18Osw could result from a combination of these factors as well.

Potential pH effect on δ18Ob and species-specific isotope fractionation

In addition to δ18Osw, other non-thermal factors have been considered to influence the foraminiferal δ18O composition. An underestimation of δ18Osw may not be necessary to explain the apparent cool temperature bias of δ18Ob. In fact, the assumed δ18Osw value of −1‰ VSMOW may still be supported by our warm Δ47-based deep-sea temperatures if a pH effect on the incorporation of the oxygen isotopes in the foraminiferal shell is taken into account. Boron isotope analyses on planktic foraminifera indicate that pH at the sea surface of the Eocene, which is directly linked to the high atmospheric CO2 concentrations, was approximately 0.5 units lower than today (~7.7 versus ~8.2, respectively)1,2,100. At the moment, such absolute pH reconstructions are not available for the deep ocean. The pH at the seafloor does not necessarily follow the sea surface pH in exactly the same way, because degradation of organic matter, carbonate dissolution, and ocean circulation affect the distribution of dissolved inorganic carbon species in the deep ocean101,102,103. Despite this, combined proxy-model and model studies suggest lower pH in the Eocene Atlantic deep ocean as well104,105,106. Theory predicts that pH influences the oxygen isotope composition of foraminiferal calcite27,28, which is confirmed in culture experiments29. Decreasing pH generally causes foraminiferal δ18O to increase and hence biases temperature to cooler values, but the magnitude of this effect has been found to vary between species29. In theory, clumped isotopes are also affected by pH to some extent (although much less than oxygen isotopes), which would work in the same direction on the reconstructed temperatures47,50. However, this pH effect on Δ47 has not been observed yet in experiments48,49. Important to mention is that a pH effect on foraminiferal δ18O has been successfully demonstrated for planktic foraminifera29, yet it is uncertain whether this effect also exists in benthic foraminifera77. Nonetheless, the underestimated early Eocene deep-sea temperatures shown by δ18Ob are consistent with a scenario of highly reduced seawater pH.

Changes and uncertainties in species-specific isotope fractionation on geological time scales are also important to consider as a possible explanation for the relatively cool δ18Ob-based deep-sea temperatures. While core top and deep time studies suggest that species-specific effects do not influence the Δ47 signal of foraminifera51,52,53,59,62, offsets in δ18O (and δ13C) are evident among different benthic foraminiferal species30,74. These interspecies offsets in δ18O also appear to change over time. For instance, a δ18O offset of 0.64‰ is recorded between Cibicidoides and Oridorsalis in the late Cenozoic74, while this offset is different in the early Cenozoic, i.e. 0.28‰ (ref. 30). These changing isotopic offsets suggest that vital effects may not have been constant over the Cenozoic because of evolutionary changes in the foraminifera30. These offsets are however difficult to quantify, considering that changes in pH and its potential effect on δ18Ob (that may be species-specific as well29) may also have played a role through geological time30. Hence, it is complicated to establish an assumed species oxygen isotope equilibrium value in deep time with certainty. Another problem is that there seems to be no commonly accepted view on which modern benthic foraminiferal species precipitate their skeletons in equilibrium with seawater, in part because the physiological processes in the fractionation of oxygen isotopes during calcification are not well understood77. As a consequence, conflicting methods are used in the literature to reconstruct deep-sea temperatures from Cenozoic δ18Ob compilations12,16,32.

Deep-sea warming during hyperthermal events

Our reconstructed background and average hyperthermal peak temperatures (13.5 ± 1.9 °C and 16.9 ± 2.3 °C, respectively) from clumped isotopes imply average deep-sea warming of 3.4 ± 2.9 °C (95% CI) during ETM2–H2. Oxygen isotope thermometry may underestimate deep-sea temperature increase during hyperthermals when considering a pH effect on δ18Ob. Boron isotope measurements on planktic foraminifera from the PETM and ETM2 indicate further acidification of the sea surface of approximately −0.3 and −0.2 pH units, respectively107,108,109. In absence of deep-sea pH proxy reconstructions, numerical models show that this drop in surface pH is propagated into the deep ocean105. The effect of an even lower deep-sea pH would potentially dampen the negative δ18Ob excursion of these events to about three-quarters of the original magnitude, resulting in a significant underestimation of the deep-sea warming29,105. While the increase in temperature from Δ47 and δ18Ob (2.8 ± 0.2 °C 95% CI based on the latter proxy) in our reconstructions seems to be fairly similar (Fig. 2b), we note that the temperature uncertainty in the Δ47-based temperatures does allow the possibility that the true deep-sea warming is larger than indicated by δ18Ob. In that case, the oxygen isotopes would underestimate the temperature rise by one degree or more.

Alternatively, potential changes in deep ocean circulation during these events may also have influenced the magnitude of the δ18Ob excursion, as these circulation changes would have affected δ18Osw of the deep Atlantic. Changes in grain size and benthic faunal assemblages from Atlantic sites are interpreted to indicate a switch in the Atlantic circulation during ETM2 (refs. 110,111). If such circulation changes were accompanied by an increase in salinity (and hence higher δ18Osw) of deep water masses, this would represent an alternative mechanism that causes an apparent δ18Ob underestimation of the hyperthermal deep-sea warming. Future efforts in reconstructing pH of the early Eocene deep ocean, investigating pH effects on benthic foraminiferal δ18O, development of salinity proxies112,113, as well as more Eocene deep-sea temperature reconstructions from clumped isotopes are needed to provide better constraints on the effects of pH and salinity on δ18Ob records.

Hyperthermal deep-sea temperature reconstructions based on Mg/Ca are currently only available for the PETM and indicate 4–5 °C warming31. A direct comparison of the magnitude of negative δ18Ob excursions shows that the average ETM2–H2 excursion (0.6‰) is about one-quarter the size of that of the PETM (~2.5‰) (refs. 65,66,114). In this respect, the average deep-sea warming of 3.4 ± 2.9 °C (95% CI) that we estimate for ETM2–H2 from clumped isotopes is larger than the warming one would expect when taking into account the recorded PETM temperature rise from Mg/Ca. This discrepancy implies that either the current Mg/Ca-based temperatures largely underestimate the deep-sea warming during the PETM or that the relative magnitude of the δ18Ob excursions between ETM2–H2 and PETM cannot be simply extrapolated in terms of temperature.

Secular changes in Mg/Casw are likely to act on much longer time scales than hyperthermal variability33. An underestimation of hyperthermal deep-sea temperature warming by Mg/Ca may result from a lowering of the Mg/Ca-temperature sensitivity at low Mg/Casw values in the early Eocene37,44. However, effects of much lower seawater pH on foraminiferal Mg/Ca may act in the opposite direction on the reconstructed magnitude of temperature change and might therefore cancel out the effect of a lower Mg/Ca-temperature sensitivity44. The relative importance of these two factors is however uncertain, making it difficult to constrain to which extent they balance each other out44.

Alternatively, Mg/Ca-based deep-sea warming of the PETM may be accurate if the ETM2–H2 and PETM were characterized by similar magnitudes of deep-sea warming despite their different δ18Ob excursions. The large negative δ18Ob excursion of the PETM might be a result of a large contribution of a decrease in δ18Osw, caused by a circulation switch to a much fresher water mass during the event31. However, numerical models contradict such a scenario and instead find a reduction in deep water formation in the southern high latitudes and development of source regions in the saline surface waters of the Tethys Ocean during the PETM (refs. 85,86). Combined clumped isotope and Mg/Ca reconstructions across the PETM and/or ETM2 events are required to better assess the performance of Mg/Ca in estimating deep-sea temperature variability in the Eocene.

Conclusions and further implications

In short, our study shows that robust marine temperature reconstructions from foraminiferal oxygen isotopes, and Mg/Ca as well, can only be achieved with independent constraints on non-thermal effects. Without these constraints, the reliability of the classic reconstructions of the early Eocene12,16,32 should be questioned. Our findings corroborate similar findings in other periods of the Cenozoic, such as the Miocene62,63,64. These findings warrant a revision of the traditional Cenozoic history of deep-sea temperatures that has provided us with a presumed benchmark of global climate evolution over the past 66 Myr (refs. 12,16,32).

Finally, our warmer temperatures may have implications for estimates of global mean temperatures and climate sensitivity in the early Eocene if we assume that our reconstructed South Atlantic deep-sea temperatures are representative of the global deep ocean. Available middle Eocene (~43 Ma) Atlantic Δ47-based reconstructions from multiple sites indicate almost identical deep-sea temperatures between the North and South Atlantic61. Reconstructions from the other ocean basins (e.g. Pacific) are needed to further test this assumption of globally uniform deep-sea temperatures in the Eocene. Yet, if the early Eocene deep ocean was warmer on a global scale than previously thought, these warmer temperatures would result in higher estimates of global mean temperatures and climate sensitivity than previously obtained from δ18Ob (refs. 7,9,15,16). Such higher estimates would be more consistent with constraints from sea surface temperature reconstructions2,9. Through providing independent and accurate temperature reconstructions, clumped isotope thermometry opens up new opportunities to evaluate existing views on the homogeneity of the deep ocean over the Cenozoic20,21,87.

Methods

Foraminiferal preservation and sampling

We measured stable and clumped isotopes on monospecific multi-specimen samples of Nuttallides truempyi and Oridorsalis umbonatus, which were picked from the >150 µm sediment size fraction of the samples of Stap et al.66. Grouping Δ47 measurements from different benthic species to produce one temperature estimate requires the assumption of absence of foraminiferal species-specific effects on Δ47, which is supported by both core top and deep time studies51,52,53,59,62. Scanning electron microscope (SEM) images were previously generated to examine the preservation in the studied interval and show that the foraminiferal tests are moderately altered by recrystallization (formation of secondary calcite)66. Nevertheless, the effect of recrystallization of primary calcite on the Δ47 signal in benthic foraminifera is considered negligible61,62,115, since this diagenesis usually occurs during early burial116,117, and therefore records similar temperatures as the seafloor temperatures. Due to limiting numbers of picked foraminifera per sample, we needed to pool foraminifers from neighboring samples to obtain sufficient foraminiferal material (~25 specimens per measurement), which resulted in averaging in the depth domain. Before isotope analysis, foraminiferal samples were cleaned by removing adherent contaminants, for instance, nannofossils, organics, and clay particles. The foraminifers were gently cracked between two glass plates to enable any chamber infill to escape and subsequently ultrasonicated in deionized water at least two times for 30 seconds. The test fragments were rinsed until the suspended particles were removed (at least two times after one ultrasonication) and oven-dried at 40 °C for one night.

Clumped isotope analysis

In total, 333 successful stable and clumped isotope sample measurements were generated across the ETM2 and H2 intervals. Our analytical methods generally follow the methodology reported in previous work54,55,56,57,58, and a more detailed description can be found in the Supplementary Methods. We measured samples of each ~70–100 µg foraminiferal material on a Thermo Scientific Kiel IV (Thermo Fisher Scientific, Waltham, MA, USA) carbonate preparation device coupled to a Thermo Scientific MAT 253 Plus mass spectrometer at Utrecht University. In addition, we measured an approximately equal number of carbonate standards. We employed three main standards (ETH-1, 2, and 3 in an approximate 1:1:5 ratio, respectively118), which differ in δ13C, δ18O, and Δ47 composition, to correct the sample measurements. Two additional check standards (IAEA-C2 and Merck) were measured in each run to monitor the long-term reproducibility of the instrument. External reproducibility (1 standard deviation) in Δ47 of IAEA-C2 after correction was 0.033‰. The δ13C and δ18O values (reported relative to the VPDB scale) of IAEA-C2 showed an external reproducibility (1 standard deviation) of 0.05‰ and 0.12‰, respectively119.

Temperature and δ18Osw calculations

As the relationship between Δ47 and temperature is nonlinear (Eq. (1)), we first determined three average values from the final corrected Δ47 measurements of the background, slope, and hyperthermal peak bins. To calculate deep-sea temperatures from these three mean Δ47 values, we employed a recent foraminifer-based calibration53, recalculated to the InterCarb-Carbon Dioxide Equilibrium Scale (I-CDES) in Meinicke et al.120 using the recently updated values of the carbonate standards:55

$${\Delta }_{47}({{{{{\rm{I}}}}}}{\mbox{-}}{{{{{{\rm{CDES}}}}}}}_{90^\circ {{{{{\rm{C}}}}}}})=(0.0397\pm 0.0011)\times {10}^{6}/{T}^{2}+(0.1518\pm 0.0128)$$
(1)

in which temperature is given in K. This Δ47-temperature relationship represents a composite of previously published foraminifer-based calibrations51,52 with new foraminiferal data53 and has a larger data density in the range of ocean temperatures compared to other calibrations. Within uncertainty, this regression is indistinguishable from calibrations based on inorganic calcite54,121,122.

Calculation of the δ18Ob-based deep-sea temperatures follows the approach that is used in the reconstruction of Cramer et al.32 (see discussion in Cramer et al.32 and Marchitto et al.77). First, we converted δ18O of N. truempyi and O. umbonatus towards Cibicidoides values for assumed seawater equilibrium using the correction factors from Katz et al.30 [δ18OCib = (δ18ONutt + 0.10)/0.89 and δ18OCib = δ18OOrid − 0.28]. Subsequently, deep-sea temperatures were determined following the δ18O-temperature relationship of Cibicidoides. In our study, we applied the improved calibration for this genus of Marchitto et al.77 (Eq. (2)), which is based on a large dataset existing of both previous and new core top measurements:

$${\delta }^{18}{{{{{{\rm{O}}}}}}}_{{{\rm{b}}}} (\textperthousand \,{{{{{\rm{VPDB}}}}}})-{\delta }^{18}{{{{{{\rm{O}}}}}}}_{{{{{{\rm{sw}}}}}}}(\textperthousand \,{{{{{\rm{VSMOW}}}}}})+0.27\\ =(-0.245\pm 0.005)\times T +{\,}(0.0011\pm 0.0002)\times {T}^{2}+(3.58\pm 0.02)$$
(2)

whereby temperature is represented in °C. We used a fixed value of −1‰ VSMOW for the mean δ18Osw representing ice-free conditions26. This quadratic regression shows good agreement with the regression of Kim & O’Neil123, which is based on inorganic calcite precipitated between 10 and 40 °C (refs. 75,123).

Finally, we calculated δ18Osw for the background, slope, and hyperthermal bins, which allows for testing the traditional assumptions on δ18Osw in an ice-free world. In Eq. (2), we used the average δ18Ob (corrected to Cibicidoides) values of the three climate states in combination with the mean Δ47-based temperatures to calculate the three corresponding δ18Osw values.

Uncertainties and error propagation

For our data representation and visualization, we present uncertainties on the three average Δ47-based temperatures and calculated δ18Osw values as 68% and 95% CI. For the Δ47-based temperatures, these error bars represent combined analytical (spread of the group of Δ47 measurements) and calibration uncertainties. To propagate uncertainties in the Δ47-based temperatures, one must first know the covariance of the slope and intercept of the Δ47-temperature calibration of Meinicke et al.53 (recalculated in Meinicke et al.120), as the slope and intercept of the regression line are dependent on each other. Niklas Meinicke calculated 100,000 regression slope-intercept pairs from the recalculated calibration dataset using a York bootstrap fit, similarly as performed in Meinicke et al.53 for the original calibration53,124. Using these values, we calculated the covariance between the calibration slope and intercept and subsequently propagated the analytical and calibration uncertainties into Δ47-based temperature using provided Matlab scripts (Supplementary Data 3 and acknowledgements). The temperature error propagation follows the approach of Huntington et al.125, which is based on standard error propagation formulas for correlated Gaussian errors (for a detailed explanation of the math see supporting information of Huntington et al.125). Most of the uncertainty in the Δ47-based temperatures is due to the analytical precision of the measurements and only a very small part comes from calibration uncertainty (~6% in our error bars)51,63,125. To determine uncertainties in the calculated δ18Osw, we used the Δ47-based temperatures at 68% and 95% CI in Eq. (2) to focus on the uncertainty that is derived from the precision of the Δ47 measurements, as this is the main contributor to the uncertainty in δ18Osw.

Regarding the temperatures from δ18Ob, analytical uncertainties on the temperatures of the three bins are very small (±0.1 °C 95% CI) in comparison to Δ47 and hence fall within the plot symbols (Fig. 2b). The dashed red envelope in Fig. 2b represents the uncertainty (95% CI) of the δ18O-temperature calibration of Marchitto et al.77. In the main text, we mention uncertainties on δ18Ob-based temperatures of the three bins as combined analytical and calibration uncertainties (in 95% CI).