The Relationship Between the Global Mean Deep-Sea and Surface Temperature During the Early Eocene

Estimates of global mean near-surface air temperature (global SAT) for the Cenozoic era rely largely on paleo-proxy data of deep-sea temperature (DST), with the assumption that changes in global SAT covary with changes in the global mean deep-sea temperature (global DST) and global mean sea-surface temperature (global SST). We tested the validity of this assumption by analyzing the relationship between global SST, SAT, and DST using 25 different model simulations from the Deep-Time Model Intercomparison Project simulating the early Eocene Climatic Optimum (EECO) with varying CO 2 levels. Similar to the modern situation, we find limited spatial variability in DST, indicating that local DST estimates can be regarded as a first order representative of global DST. In line with previously assumed relationships, linear regression analysis indicates that both global DST and SAT respond stronger to changes in atmospheric CO 2 than global SST by a similar factor. Consequently, this model-based analysis validates the assumption that changes in global DST can be used to estimate changes in global SAT during the early Cenozoic. Paleo-proxy estimates of global DST, SST, and SAT during EECO show the best fit with model simulations with a 1,680 ppm atmospheric CO 2 level. This matches paleo-proxies of EECO atmospheric CO 2 , indicating a good fit between models and proxy-data.


Introduction
An important indicator of past, present, and future climate states is the global mean surface air temperature (global SAT): the mean near-surface air temperature above both land and sea.Current insights into the global SAT evolution during the Cenozoic era, as presented in IPCC AR6 (Gulev et al., 2021), rely on estimates of deep-sea temperature (DST) from oxygen isotope compositions of benthic foraminifera (Hansen et al., 2013;Westerhold et al., 2020; Figure 1c).The DST estimates are translated into global SAT estimates by assuming a 1-on-1 relationship between changes in global mean DST and SAT prior to the Pliocene using the following line of reasoning of Hansen et al. (2013): 1.The global DST is essentially determined by the sea surface temperature (SST) at high-latitude sites of deep-water formation, where surface water sinks to the bottom.Hence, global DST exhibits the same sensitivity to changes in radiative forcing (like the atmospheric CO 2 concentration) as the high-latitude SST. 2. Due to polar amplification, the high-latitude SST, and hence also global DST, responds stronger to changes in radiative forcing than the global mean SST. 3. The surface temperature change is greater over land than over the ocean.Hence, global SAT responds stronger to changes in radiative forcing than global SST. 4. Global DST and SAT are equally sensitive to changes in radiative forcing, when assuming that global SST underestimates the sensitivities of global DST and SAT by the same factor.
The validity of these underlying assumptions is essential to our understanding of past and future climate states.
Model-based analyses of the relationship between these global mean temperature indicators are sparse due to the long runtime needed to obtain energy equilibrium in the deep-sea: model simulations by Li et al. (2013) show that the full ocean requires more than 4,500 yr of climate simulation to equilibrate.Valdes et al. (2021) used an extensive set of long-run climate model simulations of different time-periods covering the whole Phanerozoic era to investigate the relationship between DST and global SAT.They found an overall linear relationship between the two over the Phanerozoic, yet with variable slope and coefficient of determination depending on the geologic time period under consideration.
To build on this analysis, we investigated the relationships between global DST, SST, and SAT for one specific time period during the Cenozoic: the early Eocene Climate Optimum (EECO,.The EECO is the warmest interval of the Cenozoic, with atmospheric CO 2 concentration reaching values between 1,150 and 2,000 ppm (Anagnostou et al., 2020;Rae et al., 2021).In case of continued high CO 2 emissions (Shared Socioeconomic Pathway 5-8.5), the atmospheric CO 2 concentration can reach comparable levels in the 22nd century (Chen et al., 2021).Reconstructions of the EECO climate can deepen our understanding of how the climate system operates under such high CO 2 conditions (Burke et al., 2018;Tierney et al., 2020), but the EECO does not represent a direct analog for the future climate due to differences in land cover and continental configuration (Figures 1a and 1b).In terms of ocean circulation, the North Atlantic Ocean was narrower than nowadays, limiting the connectivity with the Arctic Ocean (Herold et al., 2014;O'Regan et al., 2011).Furthermore, the oceanic gateway between Antarctica and Australia was more restrictive compared to today, inhibiting the circumpolar current (Bijl et al., 2011;Cramwinckel et al., 2020;Herold et al., 2014;Sauermilch et al., 2021).
We used an ensemble of model simulations of the EECO climate assuming different atmospheric CO 2 levels, performed as part of the Deep-Time Model Intercomparison Project (DeepMIP; Lunt et al., 2021Lunt et al., , 2017)).Hence, we analyzed the global temperature relationships under fixed non-CO 2 climatic factors (such as continental configuration and astronomical conditions) and variable CO 2 level.Most of the DeepMIP model simulations have a run-time of 4,500 yr or more, allowing us to study the long-term relationship between global DST, SST, and SAT.In addition to this model-based analysis, we performed a model-data comparison, building on the analysis by Lunt et al. (2021).They found that the DeepMIP simulations with CO 2 levels of 1,120 and 1,680 ppm (four and six times preindustrial) are consistent with paleo-proxies of both global SAT and CO 2 level during the EECO.In our analysis, we additionally consider paleo-proxies of global DST and SST and analyze the model-data alignment for the combination of the three global mean temperature indicators and the CO 2 levels. 10.1029/2022PA004532 3 of 18

Output of EECO Model Simulations
As part of the DeepMIP, model simulations for the EECO were performed with nine Earth System Models based on a standard set of boundary conditions (Lunt et al., 2021(Lunt et al., , 2017)).Each model performed a preindustrial (PI) control run, and several EECO simulations at various atmospheric CO 2 concentrations ranging from the PI CO 2 concentration of 280 ppm to nine times this baseline (2,520 ppm; Table 1).The model simulation output data sets contain monthly averages of variables over the last 100 simulation years on the model-specific grid (Tables 1 and 2), except for the NorESM ocean variables which are available as annual instead of monthly averages.Additionally, we used ocean temperature time series that were available for all EECO simulations except that of INMCM, with the underlying horizontal ocean cross-section differing per model (Table S1 in Supporting Information S1).

Data of Modern Ocean Temperature
Estimates of present-day global DST, SST, and SAT were added as a reference in the analysis (Table S2 in Supporting Information S1).Modern global DST and SST were estimated from three data sets of global current ocean temperature data (E.U.Copernicus Marine Service Information, 2021a, 2021b, 2021c).An estimate of modern global SAT was derived from the ERA5 reanalysis data set (Hersbach et al., 2019).Note.The resolution presented for the ocean grid is the nominal resolution, as most models use a non-rectangular grid for their ocean component (such as a bipolar, tripolar, or displaced pole grid).The depth levels of all models become thicker with increasing depth.Note that several models start with 1 run and branch off other runs after a spin-up period.The run-time presented is the total run-time, including possible spin-up.Detailed information on the simulations of each model is provided in Lunt et al. (2021).Additionally, Zhang et al. (2020) and Zhu et al. (2019) serve as references for the DeepMIP simulations of CESM and IPSL respectively.The models are referred to by the specified short name in the remainder of this article.

EECO Paleo-Proxies
Several existing proxy compilations were used for the paleo-proxy estimates of global DST, SST, SAT, and CO 2 level during EECO.As a basis for own estimates of global SAT and SST, we used the data set with 27 marine and 70 terrestrial surface temperature proxies compiled as part of DeepMIP (Hollis et al., 2019), and as given in the Supporting Information of Inglis et al. (2020) (Table S3 in Supporting Information S1).Additionally, we included the estimated global SAT ranges from Inglis et al. (2020) and Zhu et al. (2019).The paleo-proxy range for global DST during EECO used in our analysis is based on 5 DST estimates that fall in the early Eocene interval of elevated deep ocean temperatures from the Supporting Information of Meckler et al. (2022) (Table S4 in Supporting Information S1).Finally, the range for atmospheric CO 2 level during the EECO is based on the studies of Anagnostou et al. (2020) and Rae et al. (2021).Details of these paleo-proxies are given in Text S1.1 in Supporting Information S1.

Evaluation of Model Simulation Energy Equilibrium
We assessed the degree of energy equilibrium obtained by the model simulations to determine whether the model-based temperature estimates can be regarded as equilibrium values.We based our model simulation selection on an inspection of three indicators: • Model simulation run time.Based on the results of Li et al. (2013), we set the minimal run time of the simulations at 4,500 yr.• Global annual mean ocean temperature trend at the end of the simulation.We set the criterion for a sufficiently small trend at an absolute trend smaller than 0.5°C/1,000 yr.• Global annual mean top-of-atmosphere (TOA) energy imbalance over the last 100 simulation years.For the criterion regarding a sufficiently small TOA energy imbalance, we followed Lunt et al. (2021) who used a value less than 0.3 W/m 2 or a value similar to that of the PI control run.
The global annual TOA energy imbalance per model simulation was calculated as the temporally and spatially weighted average of the difference between incoming shortwave radiation and outgoing shortwave and longwave radiation (SW↓ -SW↑ -LW↑, Table 2).It primarily reflects the energy equilibrium at the surface, while the run time and ocean temperature trend determine the degree of energy equilibrium in the deep ocean.
The model simulations of CESM, INMCM, IPSL, and NorESM do not meet the minimal run time requirement (Table S5 in Supporting Information S1).Furthermore, the model simulations of CESM and NorESM 4xCO 2 have an absolute ocean temperature trend greater than 0.5°C/1,000 yr at the end of the simulations (Table S5 in Supporting Information S1).For CESM 6xCO 2 and 9xCO 2 and NorESM 4xCO 2 , the deep ocean is still taking up heat, while the simulations CESM 1xCO 2 and 3xCO 2 show a considerable cooling trend.Lastly, the TOA energy imbalance criterion does not hold for CESM 3xCO 2 , 6xCO 2 and 9xCO 2 , INMCM and both runs of IPSL (Table S5 in Supporting Information S1).Note that COSMOS and MIROC have a TOA energy imbalance above 0.3 W/m 2 , which is however similar in value to that of the PI control run.Combining these findings, we conclude that: • The simulations of COSMOS, GFDL, HadCM High and Low Resolution, and MIROC are close to equilibrium.
• The simulations of CESM, IPSL, and NorESM have not reached full equilibrium yet.
• There is insufficient information to correctly interpret the degree of equilibrium of INMCM.
Based on this assessment, we decided to exclude INMCM from the analysis and perform the analysis for two sets of model simulations: the simulations close to equilibrium only and the full set of simulations (except INMCM).
The results of INMCM are not used in the analysis but are included in the figures and tables, for reference and completeness.

Estimation of Model-Based Global Temperatures
The global DST per model simulation was estimated using the monthly mean potential sea water temperature data (Table 2).First, the annual mean temperature was calculated as the average over the 12 months for all models (except NorESM).We defined the deep-sea as the part of the ocean ranging from a depth of 2,000 m to the ocean floor.This demarcation corresponds with the study of deep ocean warming by Desbruyères et al. (2016) and the assessment of ocean heat uptake in IPCC AR6 in which the ocean below 2,000 m is considered separately from the upper layers (Gulev et al., 2021).As most model depth levels do not exactly coincide with this boundary, we used a model-specific deep-sea definition starting at the depth level closest to 2,000 m (Table S1 in Supporting Information S1).The global DST was calculated as the volume weighted average temperature of the 3-dimensional grid cells ranging from the model-specific deep-sea depth level to the bottom (Figure S1 in Supporting Information S1).The global SST and SAT estimations were calculated similarly to those of global DST, yet with area weighted averaging as the variables are given on a 2-dimensional spatial grid (Table 2).
The methods described are also used to estimate the modern values of global DST, SST, and SAT (Table S2 in Supporting Information S1).
For each model simulation, the annual mean local DST was plotted on a (longitude, latitude) grid on a Robinson projection to visually inspect the results and assess the usability of the data for the intended analysis.We made two changes to the data underlying the analysis based on this inspection.First, there are small outlier sections in some EECO model simulations, where the DST values were several °C higher or lower than those of the surrounding area.These sections span a couple of grid cells in the Arctic gateway east of Greenland for COSMOS, HadCM High Resolution, INMCM, and MIROC, and between the Australian and Antarctic continents for COSMOS and IPSL.The outliers were manually removed from the data.Second, we chose to exclude the Arctic Ocean from the calculations of global DST and SST.Inspection of the model ocean bathymetry showed that the Arctic was poorly connected with the rest of the ocean during EECO, which limits deep-water exchange.For example, in the CESM simulations the bulk of the gateway has a depth of about 50 m (Figure 1a).The Arctic was removed from the data by excluding the surface grid cells with latitudes above 65°N.We performed sensitivity analyses to assess the impact of our definition of the deep-sea starting at 2,000 m, and the exclusion of the Arctic from the analysis.
Note that global SAT is defined as the temperature of the air above the land and ocean.A related and commonly used global temperature indicator is the global mean surface temperature, defined as the mean temperature of the air above the land surface and the temperature of the water at the sea surface.Global mean surface temperature and global SAT are strongly related, but physically distinct temperature estimates and long-term differences between them are estimated to be at most 10% (Gulev et al., 2021).As the analysis of Hansen et al. (2013) does not specify whether the paleo proxy DST estimates are translated into global mean surface temperature or global SAT, we henceforth make no distinction and use the abbreviation global SAT to indicate either estimate.

Assessment of Spatial Variability of DST
The spatial variability of the DST per simulation was analyzed through the weighted distribution of absolute differences between local DST and global DST for all (longitude, latitude) locations defined by the surface grid cell, weighted by the volume of the deep-sea below the corresponding surface grid cell.We assessed this distribution by considering several percentiles.Due to the weighting, the interpretation of the xth percentile is that x% of the total deep-sea volume has an absolute difference between local DST and global DST less than the percentile value (for x between 0 and 100).Furthermore, we calculated the annual mean DST for the five main oceanic basins: North Atlantic Ocean, South Atlantic Ocean, Southern Ocean, Pacific Ocean, and Indian Ocean (Figure S2 in Supporting Information S1).The basin borders are defined such that coastal areas are excluded, in order to obtain stable open ocean basins.The basins are comprised of one or more rectangular (longitude, latitude) boxes following the demarcation in Table S6 in Supporting Information S1.This demarcation is applied to all models equally, except for NorESM.This model uses a slightly different paleogeography, and the basin demarcations are adjusted accordingly.The difference between the DST per basin and the global DST is considered both per model simulation and grouped per oceanic basin for all model simulations together.

Estimation of Global Temperature Relationships
We analyzed the relationships between the global temperature indicators by means of linear regression.Four temperature combinations were considered to test the assumptions set forth by Hansen et al. (2013): • Antarctic winter SST versus global DST for assumption 1, • Global SST versus global DST for assumption 2, • Global SAT versus global SST for assumption 3, • Global SAT versus global DST for assumption 4.
For each combination, two linear regressions were performed: one based on all EECO simulations (except INMCM), and one excluding the simulations of CESM, IPSL, and NorESM (and INMCM).The linear regression analyses result in expression of the first temperature (the dependent variable) as a linear function of the second temperature (the independent variable):  dependent = slope ⋅ independent + intercept.Estimates of the slope and intercept coefficients were obtained via the ordinary least squares method.The fit of the result was assessed by the coefficient of determination R 2 .

Estimation of Proxy-Based Global Temperatures
Proxy-based global SST and SAT estimates were obtained by applying the method D surf -2 of Inglis et al. (2020) using the data set of 27 marine and 70 terrestrial local temperature proxies (Table S3 in Supporting Information S1).This method requires two model simulations with different CO 2 levels (low and high CO 2 ).Assuming a linear relationship between the local and global temperature estimates in the two model simulations, the local proxy-based temperature estimate is scaled to an estimate of global mean temperature.The 95% confidence intervals (CI) were estimated following the approach of Farnsworth et al. (2019)

DeepMIP Model Estimates of Global Temperatures
As expected, the global mean DST, SST, and SAT increase with increasing CO 2 level in the EECO simulations (Figure 2 and Table S5 in Supporting Information S1).Furthermore, the temperatures in the 1xCO 2 EECO runs are higher than those of the PI runs for all models with a 1xCO 2 simulation.This difference is due to non-CO 2 -related effects such as differences in continental configuration, vegetation, and presence of ice sheets that are incorporated in the boundary conditions of the EECO and PI simulations (Lunt et al., 2021).The average difference between the 1xCO 2 run and the PI control run is 2.1°C for global DST, 3.5°C for global SST, and 3.8°C for global SAT.Hence, the response to non-CO 2 climatic factors seems stronger at the surface than in the deepsea.However, the difference in global DST between the EECO 1xCO 2 and PI runs varies considerable between models, from as low as 0.5°C for COSMOS and MIROC to almost 5°C for CESM.This could indicate that the models vary in their sensitivity of global DST to the non-CO 2 effects.We note however that CESM 1xCO 2 still shows a cooling trend in the ocean temperature (Table S5 in Supporting Information S1), which can also explain the relatively warm global DST estimate compared to PI.
We defined the deep-sea as the part of the ocean below 2,000 m depth.To understand the sensitivity of our results to the chosen demarcation, the difference between the global DST starting at 2,000 m and global DST starting at depth levels between 1,000 and 4,000 m was analyzed (Figure S3 in Supporting Information S1).On average, the global DST increases when starting at a shallower depth level, meaning that the ocean temperature decreases with depth level.For most model simulations, the difference between the global DST estimate used in the analyses and the other global DST estimates is less than 0.6°C.The model that shows the most variability is INMCM.Other model simulations with a relatively high sensitivity regarding deep-sea definition are the CESM 9xCO 2 and the NorESM simulations.A possible explanation is that the ocean has not reached thermal equilibrium yet in the model simulation and the temperature of the (deep-) ocean is still adjusting to the prescribed CO 2 scenario.The model simulations CESM 9xCO 2 and NorESM 4xCO 2 show a considerable warming trend for the full ocean at the end of the simulation (Table S5 in Supporting Information S1), indicating that the ocean is indeed still taking up heat.The NorESM 2xCO 2 simulation does not have such a warming trend for the full ocean, yet the different ocean layers could still be exchanging energy.The relatively large difference between global DST starting at 4,000 and 1,000 m for this model simulation could indicate that the deep ocean is still warming.

Spatial Variability in DST
A characteristic feature of the deep-sea is its spatial temperature homogeneity.Essentially this means that global DST is for an important part determined by the surface temperature at the location of deep-water formation, and there is limited change in DST during the subsequent deep-ocean circulation.The location of deep-water formation areas, continental configuration, and bathymetry can cause spatial variability in DST.Insight into the spatial variability of DST during paleo climates is essential, as spatial homogeneity is an important prerequisite for using local paleo proxies of DST as indicator of global DST.
In most of the DeepMIP EECO simulations, the spatial variability for DST is limited to ±2.5°C difference with the global mean (Figure 3).There is a warm tendency in the parts of the deep-sea that are relatively shallow, such as the Mid Atlantic Ridge and coastal zones around North America and Africa.This can be clearly seen in the simulations of CESM, MIROC, and NorESM.The Southern Ocean is generally the coldest ocean basin in the DeepMIP EECO simulations.However, several models, such as COSMOS, GFDL, HadCM, and IPSL, show a relatively large temperature difference between the deep-sea south and east of the Australian continent.The spatial DST distribution of these relatively warm and cold regions appears to be model-dependent and persists for different increasing CO 2 levels within the same model.This suggests that the spatial DST distribution is an inherent characteristic of the model and likely related to the bathymetry and nature of the Tasman Seaway between Antarctica and Australia limiting throughflow (Zhang et al., 2022).
Proxy data as well as other modeling studies support such a limited flow through the Tasman Seaway.Dinoflagellate cyst assemblages indicate the existence of two separate surface currents on either side of the Tasman Seaway during the early Eocene, with a proto-Antarctic Circumpolar Current starting to establish in the Middle Eocene (Bijl et al., 2011;Cramwinckel et al., 2020).Limited connection of the gateways in the Southern Ocean  The preindustrial control runs and the modern ocean all show a relatively warm deep-sea in the North Atlantic and Mediterranean, and a relatively cold Antarctic Ocean deep-sea (Figure S4 in Supporting Information S1).Overall, the PI simulations and the modern ocean agree in the spatial pattern of DST variability.However, several models, such as CESM, GFDL, INMCM, and NorESM, show a warmer North Atlantic than the modern data.Furthermore, several models overestimate the temperature of the Arctic Ocean, such as CESM, COSMOS, HadCM, IPSL, and MIROC (Figure S4 in Supporting Information S1).
To further analyze the spatial variability of DST, Figure 4   possible cause for this spatial variability is the weak meridional overturning circulation displayed in these simulations (Zhang et al., 2022).The other models show a DST variability of no more than 2°C when excluding the outliers beyond the 99th percentile.For all model simulations except NorESM 2xCO 2 (and INMCM), 50% of the deep-sea volume differs less than 0.4°C from the global DST, 75% less than 0.7°C, and 90% less than 1.5°C.When only considering the five main open ocean basins depicted in Figure S2 in Supporting Information S1, it holds that 50% of the deep-sea volume differs less than 0.4°C from the global DST, 75% less than 0.6°C, and 90% less than 0.8°C, for all model simulations except NorESM 2xCO 2 (and INMCM) (Figure S5 in Supporting Information S1).This indicates that the open ocean basins are more homogeneous in DST than the parts of the deep-sea that lie closer to the continents.Model resolution can also impact the model results in these relatively shallow parts of the ocean near the continents, as coastal regions are more difficult to resolve on a relatively coarse ocean model grid.
Based on the difference between the DST per oceanic basin and the global DST, the mean DST of the Pacific Ocean agree best with the global DST (grouped for all models except INMCM in Figure 5 and per model simulation in Figure S6 in Supporting Information S1).This can partly be explained by the relatively large area covered by the Pacific (Figure S2 in Supporting Information S1).In most simulations, the Southern Ocean is colder, and the North Atlantic Ocean is warmer than the global average.The results for the Indian Ocean and the South Atlantic Ocean show more variability with simulation and model.For 84% of the 125 combinations of 25 model simulations (excluding INMCM) and five basins, the absolute difference between basin and global mean DST is less than 0.5°C, and for 95% of the combinations it is less than 1°C.From this analysis, we conclude that the spatial variability for DST in basins is low and that local open ocean paleoceanographic proxies of DST can be A possible cause of inter-basin temperature variations can be diffusion of heat from the upper layers to the deep ocean.This would mean that the global DST is not fully determined by the SST in areas of deep-water formation, but also by the SST of other ocean basins.Particularly in basins that are further away from the deep-water formation area, deep-water temperatures could increase through heat transfer by mixing or diffusion.In the DeepMIP EECO simulations, this could for example be the case for the North Atlantic since there is no deep-water formation in this area.An interesting topic of further research would be to investigate the effect of changing vertical mixing coefficients in model simulations of the Eocene Ocean, as has been done for the Miocene and Pliocene (Lohmann et al., 2022).If heat transfer significantly contributes to DST in warm climates such as the EECO, the deep ocean can take longer to equilibrate.In this regard, we note that several model simulations still show a cooling or warming trend in the deep or full ocean temperature at the end of the simulation.Predominantly the simulations of CESM, IPSL, and NorESM have not reached full equilibrium yet (Table S5 in Supporting Information S1).In the EECO 3xCO 2 simulations, the strength of the overturning circulation is 18 Sv or more (Zhang et al., 2022).This suggests that a complete circulation of the deep-ocean water can be expected to take less than 2,400 yr in the case of a sufficiently deep global overturning cell, based on an ocean volume of 1.335 × 10 18 m 3 .If the DST is fully determined by the SST at deep-water formation areas at high latitudes and there is no contribution of heat transfer at lower latitudes, an energy equilibrium is expected to be reached after a complete overturning circulation.Hence, the presence of a cooling or warming trend in ocean temperature after multiple millennia can indicate that diffusion of heat from the upper layers to the deep ocean plays a role in determining DST as well.This aligns with the results of Li et al. (2013).In their model simulation based on modern continental configuration and bathymetry, the global SST has reached equilibrium after about 1,200 yr and the full ocean after 4,600 yr, displaying a near-uniform warming with water depth.Hence, in this experiment, it took the deep ocean 3,400 yr to reach equilibrium after the sea-surface had fully responded to a climate perturbation.
Looking at paleo-proxy based research into the inter-basin temperature variability during the EECO, Cramer et al. (2009) found that  18 O isotopic differences between the North Atlantic, South Atlantic, Pacific, and Southern Ocean were at most 0.5‰ in the early Eocene.The  18 O are most negative for the North Atlantic Ocean, aligning with the relatively warm deep ocean in the DeepMIP model simulations.The inter-basin difference of 0.5‰ corresponds to a DST variability of at most 2°C, assuming an approximate 4°C temperature change per 1‰ change in  18 O in an ice-free world (Hansen et al., 2013;Kim & O'Neil, 1997).This proxy-based estimate is in alignment with the finding that most DeepMIP model simulations exhibit a DST variability of no more than 2°C (Figure 4).

Temperature Relationships
We checked the four assumptions underlying the reasoning of Hansen et al. (2013) by means of linear regression (Figure 6).The first assumption is that global DST is essentially determined by the SST at areas of deep-water formation.In the DeepMIP simulations, the Antarctic region is the predominant area of deep-water formation, and the formation is strongest in the Antarctic winter months June-August (Zhang et al., 2022).The existence of a single source of deep-water in the Southern Ocean during the EECO aligns with proxy-compilations of stable carbon isotope ratios and fish debris neodymium isotope signatures (Zhang et al., 2022).Therefore, we analyzed the relationship between global DST and Antarctic winter SST.The NorESM simulations are excluded from the analysis of the Antarctic winter SST, as only annual mean SST estimates are available for NorESM.We defined the Antarctic region of the Southern Ocean as the part consisting of surface grid cells with latitudes south of 60°S.Note that this simplified definition does not include all deep-water formation areas, as the GFDL model simulations also show deep-water formation in the northern Pacific (Zhang et al., 2022).The linear regression fits the model data well, based on the R 2 value of at least 0.92 (Figure 6a).The slopes of the regressions (0.97 and 0.91) indicate that there is indeed a strong relationship between Antarctic winter SST and global DST, yet Antarctic winter SST is a little less sensitive to changes in CO 2 than global DST.
The linear regressions between the global temperature indicators also fit the model data well, with R 2 values of at least 0.95 (Figures 6b-6d).The slopes of the regressions (both 0.84) between global DST and global SST imply that global DST is more sensitive to CO 2 changes than global SST, validating the second assumption of Hansen et al. (2013).Furthermore, the slopes of the regressions (1.17 and 1.14) between global SST and global SAT indicate that global SAT is also more sensitive to CO 2 changes than global SST.This is in alignment with assumption three of Hansen et al. (2013).
The fourth assumption involves the relationship between global DST and SAT.This relationship is of special interest, as the global SAT evolution during the Cenozoic presented in IPCC AR6 is based on the translation of paleo-proxies of DST into global SAT estimates (Gulev et al., 2021).For the regression based on all models except INMCM, the slope of 0.99 implies that changes in global DST relate almost 1-on-1 to changes in global SAT.This means that the slopes of the linear regressions between global DST and SST (0.84) and between global SST and SAT (1.17) cancel out each other almost completely.Regarding only the models close to equilibrium, the global SAT is about 4% less sensitive to CO 2 changes than the global DST (slope 0.96).Overall, our linear regression results align with the reasoning set forth by Hansen et al. (2013) regarding the sensitivities of global DST, SST, and SAT to a change in atmospheric CO 2 .
We assessed the sensitivity of the linear regression results by duplicating the linear regression for definitions of the deep-sea starting at 1,500, 2,500, and 3,000 m and for global DST and SST estimates including the Arctic (Table 3).The slopes and fits of the linear regression results are comparable to the results of the main analysis (Figure 6).Hence, the relationships between changes in global DST, SST, and SAT are little affected by the deepsea definition and exclusion of the Arctic.Regarding the deep-sea definition, the intercepts of the regressions with global SST and SAT show a slightly increasing trend with the starting level of the deep-sea.This agrees with the decreasing value of global DST with starting depth (Figure S3 in Supporting Information S1): the absolute difference between global SST/SAT and global DST increases when the deep-sea threshold is deeper.The intercepts of the regression between global DST and Antarctic winter SST show a slightly decreasing trend with the starting level of the deep-sea.Because the Antarctic winter SST is lower than the global DST (Figure 6), this also agrees with the decrease of global DST with starting depth.
The linear regression results based on the DeepMIP EECO simulations validate the Hansen et al. (2013) assumption for the ice-free conditions of the early Eocene but cannot immediately be translated to other geologic periods with completely different boundary conditions.The research of Valdes et al. (2021) covers the whole Phanerozoic era and thus a wide range of continental configurations.Analyzing data from this 540 Myr long period, they find a relatively strong linear relationship between global DST and SAT yet with variable slope.Their overall slope for the whole Phanerozoic is 0.64 (R 2 = 0.74), which increases to 0.67 (R 2 = 0.90) when only considering the last 115 Ma.They conclude that the relationship between global DST and Antarctic winter SST depends on the location of deep-water formation, and the strength of the overturning circulation.Furthermore, the relationship between Antarctic winter SST and global SAT depends on the strength of ice albedo feedbacks.As a result, the combined relationship between global DST and SAT changes over the Phanerozoic due to differences in climate state.S2 in Supporting Information S1).The colored plus-signs mark combinations of proxy-based global DST, SST, and SAT estimates, and the vertical and horizontal lines extending from the plus-signs indicate the corresponding 95% confidence intervals (Inglis et al., 2020;Meckler et al., 2022;Zhu et al., 2019;Text S1.2 in Supporting Information S1).
Both our analysis and that of Valdes et al. (2021) are based on long model simulations.However, we note that some of the DeepMIP simulations have not reached full energy equilibrium yet.Model intercomparison projects like DeepMIP provide a solid basis for an array of analyses of the simulated climate.Striving for energy equilibrium in the model simulations performed within the standardized boundary conditions of intercomparison projects is necessary for the in-depth analysis of the role of the deep-sea in the climate system.Hence, millennia-long model simulations, ideally ran to full equilibrium, are indispensable for an accurate analysis of the relationship between deep-sea and surface temperatures and a good test of the assumptions of Hansen et al. (2013).

Comparison of Model-Based and Proxy-Based Global Temperatures
The proxy-based estimate of atmospheric CO 2 during EECO is 1,150-2,000 ppm with a best estimate of 1,500 ppm (Anagnostou et al., 2020;Rae et al., 2021).This translates to a range of 4.1 and 7.1 times the preindustrial CO 2 level, with a best estimate of 5.4.We derived proxy-based temperature estimates of 19.3°C (95% CI: 17.6-21.0)for global DST,) for global SST,) for global SAT based on marine and terrestrial proxies, and 33.1°C (95% CI: 31.0-35.0)for global SAT based on marine proxies (Text S1.2 in Supporting Information S1).Comparing the proxy-based and model-based temperature estimates, we see that the 6xCO 2 simulations of GFDL and CESM correspond best with the proxy reconstructions (Figure 6).This CO 2 level of 1,680 ppm corresponds well with the proxy-based range of CO 2 of 1,150-2,000 ppm.Hence, there is an agreement between models and proxies concerning the relationship between global DST, SST, SAT, and CO 2 level.These results are in line with the model-data comparison of global metrics by Lunt et al. (2021).They conclude that the 6xCO 2 level agrees best with proxies of the EECO atmospheric CO 2 level and that the global SAT from the 4xCO 2 and 6xCO 2 model simulations are consistent with paleo-proxies.Based on our combined analysis of global SAT, SST, and DST, we conclude that the 6xCO 2 simulations show the best overall fit.The proxy-based global DST estimate depicted in Figure 6 is derived from clumped isotopes (Δ 47 ) data (Meckler et al., 2022), because clumped isotopes analysis has the advantage over  18 O analysis that no assumption needs to be made regarding the  18 O of seawater (Eiler, 2007;Ghosh et al., 2006).The estimated global DST for EECO is several degrees higher than the range of 12°C-14°C based on previous  18 O estimates (e.g., Cramer et al., 2011;Hansen et al., 2013).Note that Meckler et al. (2022) estimated DST mainly for the North Atlantic, which has a slightly higher mean DST in the DeepMIP model simulations (Figure 5).However, the difference between basin mean and global mean DST is limited, which is also demonstrated by the similar temperature reconstruction from the South Atlantic (Table S4 in S4 in Supporting Information S1).However, the dependence of global DST value on the starting depth of the deep-sea is limited: the absolute difference between a global DST value starting at 1,500 or 3,000 m depth and that starting at 2,000 m is no more than 0.2°C in most model simulations (Figure S3 in Supporting Information S1).This indicates that the impact of the different paleo-depths of the proxies on the DST estimate is also limited.The correspondence between the regression-derived and proxy-based estimates of global SST indicates that there is model-data agreement regarding seawater temperatures.On the other hand, the global SAT estimates inferred from the regression relationship are on the higher end of the proxy-based estimates of global SAT.Note that there is still considerable uncertainty in the estimation of global SAT from proxies, the outcome varies both with the proxies taken into account and the methodology to deduce a global mean from the local proxies.The first is shown by the difference in our estimates of global SAT based on both the terrestrial and marine proxies and only the marine proxies.The analysis of Inglis et al. (2020) exemplifies the dependence on methodology: the best-estimate global SAT is based on six different methods, with underlying best estimates ranges from a low of 22.8°C to a high of 29.8°C.Since our regression-based estimated global SAT is inferred from recent Δ 47 -estimates of global DST, it is higher than previous estimates based on paleothermometers such as  18 O (Hansen et al., 2013;Westerhold et al., 2020;Figure 1c).A reconsideration of EECO global SAT to higher estimates would impact the EECO Equilibrium Climate Sensitivity (ECS), which is defined as the equilibrium change in global SAT that results from the change in radiative forcing following a doubling of the atmospheric CO 2 concentration.Hence, by definition, the ECS is highly dependent on the estimates of global SAT and CO 2 concentration.Assuming a fixed CO 2 concentration, a higher global SAT directly results in a higher ECS.Constraining temperature and CO 2 estimates of key geologic intervals is therefore a necessity for an accurate estimation of ECS.

Conclusion
Using millennia-long model runs simulating the early Eocene Climate Optimum, we have investigated the longterm responses of global mean DST, sea-surface temperature, and surface air temperature to different atmos- We estimated that the global DST during the EECO was 19.3°C (95% CI: 17.6-21.0)based on recent Δ 47 data from the Atlantic by Meckler et al. (2022).This DST estimate is consistent with our marine proxy-based estimates of global SST of 33.7°C (95% CI: 31.5-35.7)and global SAT of 33.1°C (95% CI: 31.0-35.0)and the obtained linear regressions based on the DeepMIP model simulations.The estimated global SAT is higher than previous compilations by Inglis et al. (2020) and Zhu et al. (2019).The best fit of the global DST, SST, and SAT paleo proxy estimates is found with the 6xCO 2 DeepMIP model simulations, corresponding to a CO 2 level of 1,680 ppm.This is furthermore consistent with paleo proxies of atmospheric CO 2 during EECO of 1,150-2,000 ppm (Anagnostou et al., 2020;Rae et al., 2021).Hence, the model simulations align well with paleo proxies in their combination of CO 2 level, global DST, SST, and SAT.Further research into other Cenozoic climate states is needed to test the robustness of these results in order to use the relationship in a broader context.

Figure 1 .
Figure 1.Top: Continental configuration during the early Eocene (a) and the preindustrial (b) with shared color bars for continent elevation and ocean depth.Continental ice sheets during the preindustrial are depicted in white.Bottom: Evolution of the global mean surface air temperature during the Cenozoic from Hansen et al. (2013) and Westerhold et al. (2020).The Eocene Climatic Optimum (EECO) is indicated with gray shading and the global surface air temperature during the preindustrial with a gray dashed line.
which entails testing the skill of the method by applying it to model-based instead of proxy-based local temperature estimates.The proxy-based estimate of global DST and 95% CI were derived from the results of Meckler et al. (2022) using the calibration of Meinicke et al. (2021) and the approach of De Winter et al. (2022) to obtain a weighted average and standard error.Details of these methods are given in Text S1.2 in Supporting Information S1.

Figure 2 .
Figure 2. Global DST (left), global SST (middle), and global SAT (right) relative to the PI value as a function of CO 2 multiplication factor for the DeepMIP EECO model simulations.The temperature values of each model are depicted with a different color and symbol and the values of different simulations from the same model are connected with a line.The models that are close to equilibrium are depicted with symbols with a colored center, while the models that have not reached full equilibrium yet are shown as symbols with a white center and thick colored outline.The CO 2 multiplication factors on the horizontal axis are plotted on a logarithmic scale.The color scheme of the figures in this article, except for Figure 3, is based on Wong (2011).

Figure 3 .
Figure 3. Difference between local and global DST for the DeepMIP EECO model simulations plotted on a Robinson projection.The plots show the continents in black and parts of the ocean shallower than 2,000 m are not colored.The color scale gradually changes from dark blue at −4°C to dark red at +4°C.Differences between −11.2°C and −4°C are depicted in darkest blue and only appear in the Arctic region.Differences between 4°C and 7.5°C are depicted in darkest red and only appear in the simulations of INMCM, MIROC, and NorESM.
presents boxplots of the absolute difference between the local DST and the global DST for all the simulations, weighted by volume of the deep-sea.The COSMOS model has the most homogeneous DST, with 98% of the deep-sea volume differing less than 0.5°C from the mean.The models INMCM, MIROC and NorESM show most spatial variability in DST.Regarding NorESM, a

Figure 4 .
Figure 4. Absolute difference between local and global deep-sea temperature (DST) in °C for all DeepMIP EECO model simulations for all (longitude, latitude) locations defined by the surface grid cell, weighted by the volume of the deep-sea below the corresponding surface grid cell.

Figure 5 .
Figure 5. Difference between basin and global deep-sea temperature (DST) in °C for all the DeepMIP EECO model simulations except INMCM and the open ocean basins depicted in Figure S2 in Supporting Information S1.The Arctic Ocean (defined as all latitudes north of 65°N) is plotted for reference with a different range for the horizontal axis.

Figure 6 .
Figure 6.Relationship between (a) global deep-sea temperature (DST) and Antarctic winter sea surface temperature (SST), (b) global DST and global SST, (c) global SST and global SAT, and (d) global DST and global SAT.The regression results based on all models except INMCM are depicted in blue, and those excluding CESM, INMCM, IPSL, and NorESM in pink, with 95% confidence bands in shading and standard error (SE) of the slope in parentheses.The results of the Eocene Climatic Optimum simulations except for INMCM are depicted with symbols indicating the model, and colors indicating the CO 2 multiplication factor.The simulations included in both regressions are shown as symbols with a colored center, while the simulations excluded from the pink regression line are shown as symbols with a white center and thick colored outline.The results of the PI runs are depicted with small gray symbols, and the results of INMCM are depicted with gray crosses.Reference values for modern temperatures are depicted with white stars (TableS2in Supporting Information S1).The colored plus-signs mark combinations of proxy-based global DST, SST, and SAT estimates, and the vertical and horizontal lines extending from the plus-signs indicate the corresponding 95% confidence intervals(Inglis et al., 2020;Meckler et al., 2022;Zhu et al., 2019; Text S1.2  in Supporting Information S1).
Supporting Information S1).This indicates that our calculated global DST estimate based on North and South Atlantic data from Meckler et al. (2022) could be regarded as representative of global DST during EECO.Additional clumped isotopes reconstructions of DST from other ocean basins will need to confirm this in the future.Furthermore, note that the Δ 47 data from Meckler et al. (2022) are derived from fossils from open ocean basins with paleo-depths of 1,500 and 3,050 m (Figure S7 in Supporting Information S1; Table pheric CO 2 levels.The model simulations show limited spatial variability in DST: for the open ocean it holds that 90% of the deep-sea volume differs less than 0.8°C from the global DST.This indicates that local DST estimates can be regarded as reasonably representative of global DST.We estimate strong linear relationships between global DST, SST, and SAT for the different EECO simulations using linear regression.The results indicate that global DST and global SAT both respond more strongly to changes in atmospheric CO 2 than global SST by similar factors.Hence, the responses of global DST and SAT to atmospheric CO 2 changes are comparable in magnitude.This indicates that changes in global DST can be used to estimate changes in global SAT, validating the assumptions made byHansen et al. (2013) for the DeepMIP EECO model simulations.

Table 2
DeepMIP Variables Used in the Analysis

Table 3
Regression Results Between Global DST, Antarctic Winter SST, Global SST, and Global SAT for Different Definitions of Deep-Sea and With the Arctic Region Included Meckler et al. (2022)ECO climate based on model-simulations shows a strong linear relationship between global DST, SST, and SAT during EECO (Figure6).The model-based relationships obtained through linear regression can be used together with a proxy-based temperature estimate to translate this temperature estimate into estimates of the other two global temperature indicators.We used our estimate of EECO global DST based on data ofMeckler et al. (2022)as input to derive estimated ranges of global SST and SAT.Combining the linear relationship between global DST and SST based on all the EECO simulations (except INMCM) together with the global DST proxy-based estimate, we obtain an estimate of global SST of 34.7°C (95% CI: 33.3-36.2).The linear regression excluding CESM, INMCM, IPSL, and NorESM leads to a global SST of 34.9°C (95% CI: 33.5-36.3).These ranges show good agreement with our derived proxy-based global SST estimate.This indicates that the model-based linear relationship is consistent with the combination of proxy-based global DST and SST estimates.Repeating this calculation for global SAT gives temperatures of 33.1°C (95% CI: 31.5-34.8)(only excluding INMCM) and 33.0°C (95% CI: 31.4-34.6)(excluding CESM, INMCM, IPSL, and NorESM).Both estimates are high compared to our proxy-based global SAT estimate based on marine and terrestrial proxies and the best-estimates of 27.0°C of Inglis et al. (2020) and 29°C of Zhu et al. (2019) but align well with our global SAT estimate based on marine proxies only.