Articles | Volume 16, issue 5
Research article
26 Oct 2020
Research article |  | 26 Oct 2020

Global mean surface temperature and climate sensitivity of the early Eocene Climatic Optimum (EECO), Paleocene–Eocene Thermal Maximum (PETM), and latest Paleocene

Gordon N. Inglis, Fran Bragg, Natalie J. Burls, Margot J. Cramwinckel, David Evans, Gavin L. Foster, Matthew Huber, Daniel J. Lunt, Nicholas Siler, Sebastian Steinig, Jessica E. Tierney, Richard Wilkinson, Eleni Anagnostou, Agatha M. de Boer, Tom Dunkley Jones, Kirsty M. Edgar, Christopher J. Hollis, David K. Hutchinson, and Richard D. Pancost

Accurate estimates of past global mean surface temperature (GMST) help to contextualise future climate change and are required to estimate the sensitivity of the climate system to CO2 forcing through Earth's history. Previous GMST estimates for the latest Paleocene and early Eocene (∼57 to 48 million years ago) span a wide range (∼9 to 23 C higher than pre-industrial) and prevent an accurate assessment of climate sensitivity during this extreme greenhouse climate interval. Using the most recent data compilations, we employ a multi-method experimental framework to calculate GMST during the three DeepMIP target intervals: (1) the latest Paleocene (∼57 Ma), (2) the Paleocene–Eocene Thermal Maximum (PETM; 56 Ma), and (3) the early Eocene Climatic Optimum (EECO; 53.3 to 49.1 Ma). Using six different methodologies, we find that the average GMST estimate (66 % confidence) during the latest Paleocene, PETM, and EECO was 26.3 C (22.3 to 28.3 C), 31.6 C (27.2 to 34.5 C), and 27.0 C (23.2 to 29.7 C), respectively. GMST estimates from the EECO are ∼10 to 16 C warmer than pre-industrial, higher than the estimate given by the Intergovernmental Panel on Climate Change (IPCC) 5th Assessment Report (9 to 14 C higher than pre-industrial). Leveraging the large “signal” associated with these extreme warm climates, we combine estimates of GMST and CO2 from the latest Paleocene, PETM, and EECO to calculate gross estimates of the average climate sensitivity between the early Paleogene and today. We demonstrate that “bulk” equilibrium climate sensitivity (ECS; 66 % confidence) during the latest Paleocene, PETM, and EECO is 4.5 C (2.4 to 6.8 C), 3.6 C (2.3 to 4.7 C), and 3.1 C (1.8 to 4.4 C) per doubling of CO2. These values are generally similar to those assessed by the IPCC (1.5 to 4.5 C per doubling CO2) but appear incompatible with low ECS values (<1.5 per doubling CO2).

1 Introduction

Under high growth and low mitigation scenarios, atmospheric carbon dioxide (CO2) could exceed 1000 parts per million (ppm) by the year 2100 (Stocker et al., 2013). The long-term response of the Earth system under such elevated CO2 concentrations remains uncertain (Stevens et al., 2016; Knutti et al., 2017; Hegerl et al., 2007). One way to better constrain these climate predictions is to examine intervals in the geological past during which greenhouse gas levels were similar to those predicted under future scenarios. This is the rationale behind the Deep-time Model Intercomparison Project (DeepMIP;, last access: 21 October 2020) which aims to investigate the behaviour of the Earth system in three high-CO2 climate states in the latest Paleocene and early Eocene (∼57–48 Ma) (Lunt et al., 2017; Hollis et al., 2019).

Sea surface temperature (SST) and land air temperature (LAT) proxies indicate that the latest Paleocene and early Eocene were characterised by global mean surface temperatures (GMSTs) much warmer than those of today (Cramwinckel et al., 2018; Farnsworth et al., 2019; Hansen et al., 2013; Zhu et al., 2019; Caballero and Huber, 2013). Having a robust quantitative estimate of the magnitude of warming at these times relative to modern is useful for two primary reasons: (1) it allows us to contextualise future climate change predictions by comparing the magnitude of future anthropogenic warming with the magnitude of past natural warming; (2) combined with knowledge of the climate forcing, it allows us to estimate climate sensitivity, a key metric for understanding how the climate system responds to CO2 forcing. Using different proxy data compilations (Hollis et al., 2012; Lunt et al., 2012), the 5th Intergovernmental Panel on Climate Change (IPCC) Assessment Report (AR5) stated that GMST was 9 to 14 C higher than for pre-industrial conditions (medium confidence) during the early Eocene (∼52 to 50 Ma) (Masson-Delmotte et al., 2014). However, subsequent studies indicate a wider range of estimates, from 9 to 23 C warmer than pre-industrial (Caballero and Huber, 2013; Cramwinckel et al., 2018; Farnsworth et al., 2019; Zhu et al., 2019; Fig. 1 and Table 1). It is an open question whether this range arises from inconsistencies between the methods used to estimate GMST, such as selection of proxy datasets, treatment of uncertainty, and/or analysis of different time intervals. This methodological variability has thwarted robust comparisons between GMST methodologies for key intervals through the latest Paleocene to early Eocene.

Figure 1Published GMST estimates during the early Paleogene (57 to 48 Ma). Dots represent average values. The horizontal limits on the individual dots represent the reported error, and y-axis labels refer to previous estimates (see Table 1).


Table 1Previous studies that have determined GMST for the early Eocene (EE), EECO, PETM, or latest Paleocene (LP); n/a indicates that no error bars were reported in the original publications.

Download Print Version | Download XLSX

Here we calculate GMST estimates within a consistent experimental framework for the target intervals outlined by DeepMIP: (i) the Early Eocene Climatic Optimum (EECO; 53.3 to 49.1 Ma), (ii) the Paleocene–Eocene Thermal Maximum (PETM; ca. 56 Ma), and (iii) the latest Paleocene (LP; ca. 57–56 Ma). We use six different methods to obtain new GMST estimates for these three time intervals by employing previously compiled SST and LAT estimates (Hollis et al., 2019) as well as bottom water temperature (BWT) estimates (Dunkley Jones et al., 2013; Cramer et al., 2009; Sexton et al., 2011; Littler et al., 2014; Laurentano et al., 2015; Westerhold et al., 2018; Barnet et al., 2019). We also undertake a suite of additional sensitivity studies to explore the influence of particular proxies on each GMST estimate. We then compile GMST estimates from all six methods to generate a “combined” GMST estimate for each time slice and use these, with existing estimates of CO2 (Gutjahr et al., 2017; Anagnostou et al., 2016), to develop new estimates of “bulk” equilibrium climate sensitivity (ECS) during the latest Paleocene, PETM, and EECO.

2 Methods and materials

Three different input datasets are used to calculate GMST: (1) dataset Dsurf, which consists of surface temperature estimates, both marine (sea surface temperature) and terrestrial, (2) dataset Ddeep, which consists of deep-water temperature estimates, and (3) dataset Dcomb, which consists of a combination of surface- and deep-water temperature estimates. Here we make use of six different methodologies, which are described in detail below, to estimate GMST from these datasets.

2.1 Dataset Dsurf

Dataset Dsurf is version 0.1 of the DeepMIP database, as described in Hollis et al. (2019) (Supplement). It consists of SSTs and LATs for the latest Paleocene, PETM, and EECO. The SSTs are derived from foraminiferal δ18O values, foraminiferal Mg∕Ca ratios, clumped isotopes (Δ47), and isoprenoid glycerol dialkyl glycerol tetraethers (GDGTs) (TEX86). Foraminiferal δ18O values and Mg∕Ca ratios are calibrated to SST following Hollis et al. (2019) and Evans et al. (2018), respectively. TEX86 values are calibrated to SST using BAYSPAR (Tierney and Tingley, 2014). Δ47 values are reported using the parameters and calibrations of the original publications (Evans et al., 2018; Keating-Bitonti et al., 2011). LATs are derived from leaf fossils, pollen assemblages, mammal δ18O values, paleosol δ18O values, paleosol climofunctions, and branched GDGTs. LAT estimates are calculated using the parameters and calibrations of the original publications (see Hollis et al., 2019, and references therein). The locations of the proxy datasets are shown in Fig. S1 in the Supplement using the paleomagnetic-based reference frame (Hollis et al., 2019). For each dataset, we utilise the uncertainty range of temperature estimates reported in Hollis et al. (2019).

Four methods (Dsurf-1, Dsurf-2, Dsurf-3, and Dsurf-4) are employed to calculate GMST from dataset Dsurf. These methods employ parametric (Dsurf-1, Dsurf-2, Dsurf-4) or non-parametric (Dsurf-3) functions to estimate temperature. We calculate GMST on the mantle-based reference frame and employ the rotations provided in Hollis et al. (2019). These differ very slightly from those utilised in the DeepMIP model simulations (Lunt et al., 2020). Each method conducts a “baseline” calculation that uses the SST and LAT data compiled in accordance with the DeepMIP protocols (i.e. Hollis et al., 2019). Our baseline calculation (Dsurf-baseline; Table 2) excludes δ18O values from recrystallised planktonic foraminifera because the resulting temperature estimates are biased by diagenesis toward significantly cooler temperatures than those derived from (i) the δ18O value of similarly aged and similarly located well-preserved foraminifera, (ii) foraminiferal Mg∕Ca ratios, and (iii) Δ47 values from larger benthic foraminifera (Pearson et al., 2001; Hollis et al., 2019, and references therein). For each method, we also conduct a series of illustrative subsampling calculations relative to Dsurf-baseline based on varying assumptions about the robustness of different proxies (Table 2). The first sensitivity experiment (Dsurf-Frosty; Table 2) includes δ18O values from recrystallised planktonic foraminifera. The second sensitivity experiment (Dsurf-NoTEX; Table 2) removes TEX86 values as these give slightly higher SSTs than other proxies, especially in the middle to high latitudes (Bijl et al., 2009; Hollis et al., 2012; Inglis et al., 2015). The third sensitivity experiment (Dsurf-NoMBT; Table 2) removes MBT(')  CBT values derived from marine sediment archives as they may suffer from a cool bias (Inglis et al., 2017; Hollis et al., 2019). The fourth sensitivity experiment (Dsurf-NoPaleosol; Table 2) removes mammal and paleosol δ18O values as well as paleosol climofunctions as these proxies may suffer from a cool bias (Hyland and Sheldon, 2013; Hollis et al., 2019). For each method, GMST is calculated for (i) the Early Eocene Climatic Optimum (EECO; 53.3 to 49.1 Ma), (ii) the Paleocene–Eocene Thermal Maximum (ca. 56 Ma), and (iii) the latest Paleocene (LP; ca. 57–56 Ma).

Table 2Baseline and optional subsampling experiments applied to Dsurf.

Download Print Version | Download XLSX


Method Dsurf-1 was first employed by Caballero and Huber (2013) to estimate GMST from early Eocene surface temperature proxies after it was recognised that pervasive recrystallisation of foraminiferal δ18O could overprint the original SST signal (e.g. Pearson et al., 2001, 2007). That study used data compilations (Huber and Caballero, 2011; Hollis et al., 2012) which were the predecessors to the DeepMIP compilation (Hollis et al., 2019).

Here, the anomalies of individual proxy temperature data points with respect to modern values at the corresponding paleolocation are first calculated. The time period used is between 1979 and 2018, and we used a climatology of the full ERA-Interim period (Dee et al., 2011). The calculation involves binning into low, middle, and high latitudes (30 N to 30 S, 30 to 60 N–S, and 60 to 90 N–S) and calculating the unweighted mean anomaly within these bins between the median reconstructed value at a given locality and the temperature in the modern system (from reanalysis). The geographically binned means are then weighted according to relative spherical area to calculate a globally weighted mean temperature anomaly between the paleo-time slice and modern. All samples are treated equally and considered independent. The associated errors are added in quadrature with the inter-sample standard deviation. These two sources of error were combined and normalised by the square root of the number of samples. This method is intended as an unsophisticated, brute-force approach to estimating GMST when dealing with many localities with poorly characterised errors in which there is a large difference between the reconstructed temperature at a given location and the modern equivalent. It is not intended to identify small changes in GMST; nor is it expected to work well under conditions in which temperature gradients are stronger than today, continents are far removed from their current configuration, or systematic errors are not readily mitigated by large sample size (i.e. when there are correlations in systematic errors between proxies). It is designed to be relatively straightforward to interpret and simple to reproduce without overly relying on climate models or sophisticated statistical models.

Various sanity checks have been performed to determine if the method is likely to produce useful results for a given sampling distribution and what corrections should be applied to optimise it. For example, if the modern temperature field is sampled using a geographic sampling distribution for a given time interval, what would the reconstructed modern temperature be? Sampling the modern global annual average surface temperature field in the reanalysis product ERA-5 yields a mean value of 15.1 C, but when resampled at the equivalent geographic distribution of our samples from the latest Paleocene, PETM, and EECO yields mean values for the modern of 16.9 C (±1.8C), 14.2 C (±1.7C), and 15.2 C (±1.1C), respectively. Thus, for the sampling densities and spatial structure of the early Paleogene, this method can approach the true value within ∼1.5C and the error propagation adequately characterises the error in this “perfect knowledge” scenario. Seeking precision beyond that range is unwarranted and, as indicated above, systematic biases are a serious concern. However, estimating the latest Paleocene and early Eocene GMST may be somewhat easier than estimating the modern GMST because temperature gradients were greatly reduced compared to modern. Huber and Caballero (2011) estimate a reduction to less than half the modern temperature gradient, whilst Evans et al. (2018) constrain the low- to high-latitude SST gradient to at least ∼30 % (±10 %) weaker than modern (Evans et al., 2018).

Alongside modern observations, we can also use paleoclimate model results to characterise how well the existing paleogeographic sampling network will impact results (Fig. 2). Here we utilise two CESM1 simulations, as described in Cramwinckel et al. (2018; EO3 and EO4). The two cases are chosen to minimise the magnitude of the correction to GMST, and the final result is not sensitive to the choice of reference simulation between these two (Supplement). For each interval, the difference between reconstructed global temperatures and the true paleoclimate model mean is <1 to 3 C. These comparisons demonstrate that this method produces estimates that are within random error given otherwise perfect knowledge. The errors introduced by limited paleogeographic sampling can be alleviated by incorporating the offset in mean values between the true paleoclimate model GMST and the sampled paleoclimate model GMST outlined above (Fig. 2). We utilise this offset to correct for systematic errors, but this is the only component in which paleoclimate model information is included in this GMST estimation methodology. This approach is best applied within the context of studying the random and systematic error structure as described above, and caution should be taken in using systematic corrections that are significantly bigger than the estimated random error. The underlying assumption is that the bias in the global mean estimate that exists due to uneven sampling is the same in the “proxy” Eocene world as in the “model” Eocene world, i.e. that the zonal and meridional gradients are well characterised by the model, even if the absolute temperatures are not.

Figure 2An illustration of method Dsurf-1 during the EECO. (a) Modelled early Eocene temperatures utilising CESM1.2 at pre-industrial CO2, (b) interpolated absolute SST reconstructions, and (c) data–model difference between (a) and (b).

We note that the magnitude of the global correction could be sensitive to different models and/or boundary conditions. To explore this further, we performed the same analysis using Community Earth System Model version 1.2 (CESM1.2) at CO2. This model simulation offers a major improvement over earlier models (Zhu et al., 2019) due to the improved treatment of cloud microphysics and is able to reproduce key features of the early Paleogene (e.g. the meridional SST gradient; Zhu et al., 2019; Lunt et al., 2020). We find that CESM1 ( and 16×CO2) and CESM1.2 (CO2) yield similar GMST estimates during the PETM, EECO, and latest Paleocene. For example, GMST values (obtained using Dsurf-baseline) during the EECO average 24.5, 24.6, and 25.2 C for CESM1 (×8 CO2), CESM1 (×16 CO2), and CESM1.2 (CO2), respectively. This indicates that the final result is not overly sensitive to the choice of reference simulation, at least within the CESM family. In the following sections, we only discuss CESM1 simulations to avoid circularity if the results from this paper are used to evaluate more recent simulations (e.g. CESM1.2; Lunt et al., 2020).


GMST estimates are calculated using the method described in Farnsworth et al. (2019), in which a transfer function is used to calculate global mean temperature from local proxy temperatures. The transfer function is generated from a pair of early Eocene climate model simulations carried out at two CO2 concentrations. The first simulations are the same CO2 and CO2 HadCM3L Eocene simulations from Farnsworth et al. (2019). The second simulations are the CO2 and CO2 CCSM3 simulations of Huber and Caballero (2011), also discussed in Lunt et al. (2012). The two models are configured for the Eocene with different paleogeographies (Table S1 in the Supplement). We provide a final estimate based on the mean of our two models.

The principal assumption of this approach is that global temperatures scale linearly with local temperatures and that a climate model can represent this scaling correctly (see below). The resulting GMST estimate is therefore independent of the climate sensitivity of the model but dependent on the modelled spatial distribution of temperature. For a single given proxy location with a local temperature estimate (Tproxy), Farnsworth et al. (2019) estimate global GMST (<T>inferred) as

(1) < T > inferred = < T low > + ( T proxy T low ) < T high > - < T low > T high - T low ,

where <Tlow> and <Thigh> are the global means of a low- and high-CO2 model simulation, respectively, and Tlow and Thigh are the local temperatures (same location as the proxy) from the same simulations. Tlowand Thigh represent local modelled SSTs or local modelled near-surface LATs (in contrast to Farnsworth et al. 2019, who only used local modelled near-surface LATs to calculate Tlow and Thigh, even if Tproxy was SST). If the proxy temperature is greater than Thigh or cooler than Tlow, then the inferred global mean is found by extrapolation rather than by interpolation and is therefore more uncertain (Fig. 3). This will be sensitive to the choice of model simulation; models that simulate less polar amplification (e.g. HadCM3L) are more likely to obtain <T>inferred (i.e. GMST) via extrapolation. We repeat this process for each proxy data location (Fig. 4) and take an average over all proxy locations as our best estimate of global mean temperature.

Figure 3An illustration of method Dsurf-2 for two sites: (a) Big Bend LAT in the EECO as diagnosed using HadCM3L and (b) DSDP Site 401 SST in the PETM as diagnosed using CCSM3. The vertical dashed line shows <T>inferred and the horizontal dashed line shows Tproxy, which intercept at the orange dot. The dark blue dots show the intercept of Tlow with <Tlow>, and the red dots show the intercept of Thigh with <Thigh>.


Figure 4Inferred global mean temperature (<T>inferred) using Dsurf-2 for (a) each EECO-aged LAT proxy as diagnosed using HadCM3L and (b) each PETM-aged SST proxy as diagnosed using CCSM3. For (a) and (b), the final estimate of global mean temperature is the average of all the individual sites. The solid line shows the continental outline in each model, and the dashed line shows the continental outline.

Recent work has demonstrated that CESM1.2 and GFDL model simulations offer a major improvement over earlier models (Zhu et al., 2019; Lunt et al., 2020). As such, we also calculated GMST using CESM1.2 ( and CO2; Zhu et al., 2019; Table S1) and GFDL ( and CO2; Hutchinson et al., 2018; Lunt et al., 2020; Table S1). We find that all four simulations (i.e. HadCM3L, CCSM3, CESM1.2, and GFDL) yield similar GMST estimates. For example, GMST during the PETM ranges between 32.3 and 34.5 C (Supplement). This demonstrates that Dsurf-2 is not overly sensitive to the climate model simulation. However, as CESM1.2 and GFDL have greater polar amplification than other models (e.g. HadCM3L), GMST is more likely to be found by interpolation (compare to extrapolation). To explore whether GMST scales linearly with local temperatures, we used CESM1.2 to recalculate GMST using the same method as above but using the CO2 simulation in place of the CO2 simulation. We find that GMST estimates are very similar (±0.4C). This is because, although the relationship between GMST and CO2 is non-linear (Zhu et al., 2019), the relationship between local and global temperature is relatively constant. In the following sections, we employ CCSM3 and HadCM3 simulations to avoid circularity if the results from this paper are used to evaluate more recent simulations (e.g. CESM1.2, GFDL; Lunt et al., 2020).


For Dsurf-3, GMST estimates are calculated using Gaussian process regression (Fig. 5; Bragg et al., 2020). In this method, temperature is treated as an unknown function of location, f(x). Many possible functions can fit the available proxy dataset. By using a Gaussian process model of the unknown function, we assume that temperature is a continuous and smoothly varying function of location, and once fitted to the data, the posterior mean of the model gives the most likely function form for the temperature. We use a Gaussian process prior and update it using the proxy data to obtain the posterior model, which we can then use to predict the surface temperatures on a global grid. Prior specification of the model is via a mean function E(f(x))=m(x) and a covariance function Cov(f(x),f(x))=k(x,x) (which tells us how correlated f(x) is with f(x)). We also specify the standard deviation of the observation uncertainty about each data point (σi2). If f=fx1,fxnT is a vector of temperature observations at each location xi, then the model is

(2) f N ( μ , Σ ) ,

where μi=m(xi) and Σij=kxi,xj+Ii=jσi2. The proxy temperatures are expressed as anomalies to either the marine or terrestrial present-day zonal mean temperature at the respective paleolatitude. We subtract the mean temperature anomaly (weighted by the paleolatitude) for each time period and core experiment prior to the analysis and therefore fit the model to the residuals. This means the predicted field will relax towards the mean surface warming in areas of no data coverage. The covariance function – which considers the clustering of proxy locations – describes the correlation between f(xi) and f(xj) in relation to the distance of xi and xj. We use a squared-exponential covariance function with Haversine distances replacing Euclidean distances so that correlation is a function of distance on the sphere.

Figure 5Predicted surface warming by Gaussian process regression using Dsurf-3 for the (a) latest Paleocene, (b) PETM, and (c) EECO. Anomalies are relative to the present-day zonal mean surface temperature. Circles (triangles) indicate all available SST (LAT) proxy data for the respective time slice used to train the model. Symbols for locations where multiple proxy reconstructions are available are slightly shifted in latitude for improved visibility.

A heteroscedastic noise model is used to weight the influence of individual proxy data by their associated uncertainty; i.e. the model will better fit reconstructions with a smaller reported error. Proxy uncertainties are taken from Hollis et al. (2019). Standard deviations for TEX86, Mg∕Ca, and δ18O records are derived from the reported 90 % confidence intervals (Hollis et al., 2019). A minimum value of 2.5 C for the standard deviation is assumed for all other methods. The output variances and length scale of the covariance function are estimated using their maximum likelihood values, obtained with the GPy Python package (GPy, 2012). We apply the method to the marine and terrestrial data separately and combine the masked fields afterwards to prevent mutual interference. We further constrain the lower bound of the length scale parameter to 2000 km to always fit a reasonably smooth surface, even in some continental areas with noisy proxy data (e.g. western North America). We note that our choice of the minimum length scale and the separation of land and ocean temperatures influence the predicted regional surface temperature patterns but do not significantly change our GMST estimates.

The Gaussian process approach provides probabilistic predictions of temperature values, i.e. uncertainty estimates of the predicted field. The uncertainty reported for an individual GMST estimate is calculated via random sampling. We generate 10 000 surfaces from a multivariate normal distribution based on the predicted mean and full covariance matrix and calculate the GMST for each sample. Uncertainty of the mean estimate is then defined as the standard deviation of the 10 000 random samples. Regional model uncertainty (expressed as standard deviation fields) is typically highest in areas with sparse data coverage (e.g. the Pacific Ocean and Southern Hemisphere landmasses; Fig. S2). The lower uncertainty for the latest Paleocene relative to the PETM and EECO is related to the smaller reported uncertainties in the proxy dataset rather than enhanced data coverage. The large spread in reconstructed terrestrial temperatures for North America during the PETM and EECO (Fig. S2) propagates through into relatively large uncertainties in the GMST estimates for these intervals.


For Dsurf-4, GMST estimates are calculated using a simple function of latitude (θ) tuned to best fit the proxy data:

(3) T θ a + b θ + c cos θ ,

where T(θ) is the Eocene zonal mean temperature, and the coefficients a, b, and c are chosen to minimise the sum of the squared residuals relative to Dsurf (i.e. the SST and LAT data from Hollis et al., 2019). This new model represents T(θ) well in the modern climate (Fig. S3) when supplied with a similar number of data points as in the Hollis et al. (2019) dataset, and it ensures a global solution that is consistent with the physical expectation that temperature should decrease – and the meridional gradient in temperature should increase – from the tropics toward the poles (Fig. S3).

For each data point, we account for three types of uncertainty (i.e. temperature, elevation, latitude). For temperature, we assume a skew-normal probability distribution based on the stated 90 % confidence intervals. Where uncertainty estimates are not given, we assume a (symmetric) normal distribution with a 90 % confidence interval of ±5 K. For elevation, we assume a skew-normal distribution with a 90 % confidence interval equal to the lowest and highest elevations of adjacent grid points in the paleotopography dataset of Herold et al. (2014), with a lower bound of zero.

T(θ) was estimated by sampling temperature, elevation, and latitude from their respective distributions at each location (Fig. S4), and a lapse-rate adjustment of 6 K km−1 was applied. Then, using a standard Monte Carlo bootstrapping method, the same number of data points was resampled via replacement, and the coefficients in Eq. (3) were found that best fit the subsampled data. This procedure was repeated 10 000 times to find a probability distribution of T(θ). The uncertainty associated with an individual GMST estimate is the standard deviation.

2.2 Dataset Ddeep

Dataset Ddeep consists of benthic foraminiferal δ18O-derived bottom water temperatures (BWTs) for the latest Paleocene, PETM, and EECO. The benthic foraminiferal δ18O dataset is based on previous compilations (Dunkley Jones et al., 2013; Cramer et al., 2009), updated to include more recently published datasets (Sexton et al., 2011; Littler et al., 2014; Laurentano et al., 2015; Westerhold et al., 2018; Barnet et al., 2019). The EECO dataset is sourced from 11 sites, providing spatial coverage of the Pacific, Atlantic, and Indian Ocean (DSDP/ODP Sites 401, 550, 577, 690, 702, 738, 865, 1209, 1258, 1262, and 1263). The PETM and latest Paleocene datasets are sourced from a compilation of nine and seven sites, respectively, differing from Dunkley-Jones et al. (2013) in that (i) more recent datasets were added, and (ii) PETM sites with a muted carbon isotope excursion (CIE) magnitude (<1.5 ‰) were excluded as these datasets may be missing the core PETM interval (Table S2). Benthic foraminifera δ18O values are adjusted to Cibicidoides following established methods (Cramer et al., 2009), allowing temperature to be calculated using Eq. (9) of Marchitto et al. (2014):

(4) ( δ cp - δ sw + 0.27 ) = - 0.245 ± 0.005 t + 0.0011 ± 0.0002 t 2 + 3.58 ± 0.02 ,

where t is bottom water temperature in Celsius, δcp is δ18O of CaCO3 on the Vienna Pee Dee Belemnite (VPDB) scale, and δsw is δ18O of seawater on the Standard Mean Ocean Water (SMOW). δsw is defined in accordance with the DeepMIP protocols (−1.00 ‰; see Hollis et al., 2019).


For Ddeep-1, GMST estimates are calculated following the method of Hansen et al. (2013), which utilises only the deep-ocean benthic foraminifera δ18O dataset, and we refer the reader to that study for a detailed justification of the approach. Briefly, for time periods prior to the Pliocene, GMST is scaled directly to deep-ocean temperature. Specifically, ΔGMST = ΔBWT prior to ∼5.3 Ma, with early Pliocene BWT and GMST calculated following Eqs. (3.5), (3.6), and (4.2) of Hansen et al. (2013). As such, the calculations presented here differ from those of Hansen et al. (2013) only in that (i) we use the revised benthic δ18O compilation described above rather than that of Zachos et al. (2008) and (ii) a different equation (Eq. 4) to convert δ18O to temperature.

2.3 Dataset Dcomb

Dataset Dcomb uses a combination of (tropical) surface- and deep-water temperature estimates. The deep-ocean dataset (Ddeep) is identical to that described in Sect. 2.2. The tropical SST dataset utilises all relevant surface-ocean proxy data from the DeepMIP database, i.e. those with a paleolatitude in the magnetic reference frame within 30 of the Equator. An expanded (relative to modern) definition of the tropics is used because tropical SST reconstructions are relatively sparse; 30 was chosen because it retains tropical SST data from several proxies for all three intervals, whilst SST seasonality remains relatively low within these latitudinal bounds.


For Dcomb-1, GMST estimates are calculated for each time interval based on the difference between tropical SSTs and deep-ocean BWTs (Evans et al., 2018) such that

(5) GMST = 0.5 ( tropical SST + BWT ) .

The fundamental assumptions of this approach are that (1) GMST can be approximated by global mean SST, (2) global mean SST is equivalent to the mean of the tropical and high-latitude regions, (3) benthic temperatures are representative of high-latitude surface temperatures, and (4) the temperature gradient between the abyss and high-latitude SST is fixed through time (see Sijp et al., 2011). To test these assumptions from a theoretical perspective, we modelled the shape of the latitudinal temperature gradient using a simple algebraic function (Fig. S5). These results suggest that Dcomb-1 may underestimate GMST by 0.75 to 1.25 C in the modern. We also compared GMST from the EO3 and EO4 model simulations of Cramwinckel et al. (2018) to that calculated using Dcomb-1 (Fig. S5) and find a similar cold bias during the Eocene (∼1 to 3 C). However, we note that these findings depend on the accuracy of the modelled deep-ocean temperatures.

Probability distributions for each time interval were computed as follows. In the case of the tropical SST data, 1000 subsamples were taken, following which a random normally distributed error was added to each data point in the DeepMIP compilation, including both calibration uncertainty and variance in the data where multiple reconstructions are available for a given site and time interval. Mean tropical SST was calculated for each of these subsamples. To provide a BWT dataset of the same size as the subsampled tropical SST data, 1000 normally distributed values were calculated for each time interval based on the mean ±1 SD variation of the pooled benthic δ18O data from all sites including calibration uncertainty.

3 Results and discussion

3.1 Comparison of surface and bottom water temperature-derived GMST estimates

The following section discusses our baseline GMST estimates calculated on the mantle-based reference frame only. During the latest Paleocene and PETM, GMST estimates derived from Dsurf-baseline average ∼27 and 33 C, respectively (Table 3; Fig. 6). These values are consistent with previous studies analysing the latest Paleocene (∼27C; Zhu et al., 2019) and PETM (∼32C; Zhu et al., 2019). During the EECO, GMST estimates calculated using Dsurf average ∼27C (Fig. 6). These values are up to 3 C lower compared to previous estimates from similar time intervals (ca. 29 to 30 C; Huber and Caballero, 2011; Caballero and Huber, 2013; Zhu et al., 2019). This is likely because we use an expanded LAT dataset (n=80) compared to previous studies (n=51; Huber and Caballero, 2011). Several of these proxies saturate between ∼25 and 29 C (e.g. leaf fossils, pollen assemblages, and brGDGTs; see Hollis et al., 2019, and references therein) and/or are impacted by non-temperature controls (e.g. paleosol climofunctions; see below) and could skew GMST estimates towards lower values. To confirm this, we calculated GMST values using LAT proxies only (Supplement). We show that LAT-only GMST estimates are up to 6 C lower than our baseline (SST + LAT) calculations, suggesting that EECO GMST estimates (Dsurf-baseline) may represent a minimum temperature constraint.

Table 3Individual GMST estimates for the latest Paleocene (LP), PETM, and EECO. Reported GMST estimates utilise baseline experiments except Dsurf-1 during the EECO, which uses Dsurf-NoPaleosol. GMST estimates are based on the mantle-based reference frame. Error bars on each individual method are the standard deviation (1σ), except Dsurf-1 and Dsurf-2, which use the standard error (1σx).

Download Print Version | Download XLSX

Figure 6GMST estimates during the (a) PETM, (b) EECO, and (c) latest Paleocene for each methodology. GMST estimates utilise baseline experiments except Dsurf-1 during the EECO, which uses Dsurf-NoPaleosol. GMST estimates are based on the mantle-based reference frame. Error bars on each individual method are the standard deviation (1σ), except Dsurf-1 and Dsurf-2, which use the standard error (1σx).


GMST estimates for the latest Paleocene, PETM, and EECO, calculated using Ddeep, are 25.8 C (±1.4C), 31.1 (±2.9C), and 28.0 C (±1.3C), respectively (Table 3; Fig. 6). These estimates are comparable to those derived from surface temperature proxies alone (Table 3). GMST estimates from the EECO are also comparable to previous estimates based on globally distributed benthic foraminifera data (∼28C; Hansen et al., 2013). As benthic foraminifera are less susceptible to diagenetic alteration than planktonic foraminifera (e.g. Edgar et al., 2013), this implies that benthic foraminiferal δ18O values could be used to provide the “fine temporal structure” of Cenozoic temperature change (e.g. Lunt et al., 2016; Hansen et al., 2013). However, we also urge caution as this approach scales GMST directly to BWT prior to the Pliocene and assumes that the characteristics of polar amplification are constant through time (see Evans et al., 2018; Cramwinckel et al., 2018). Changes in ice volume may also influence the benthic foraminiferal δ18O signal (see Hansen et al., 2013), and additional corrections are required before applying this method to other time intervals (e.g. the Eocene–Oligocene transition). Ddeep also implies that vertical ocean stratification is fixed, even though vertical ocean stratification has been proposed to change dramatically in the past (e.g. Sijp et al., 2013; Goldner et al., 2014) and may shift the slope and/or intercept of the relationship between BWT and GMST.

GMST estimates for the latest Paleocene, PETM, and EECO, calculated using Dcomb, are 21.6 C (±1.2C), 26.6 (±2.1C), and 22.8 C (±1.0C), respectively (Fig. 6). These estimates are consistently lower (up to 5 C) than GMST estimates derived using Dsurf and Ddeep. Although Dcomb-1 can estimate modern GMST within ∼1 to 2  C of measured values, whether this approach can be applied in greenhouse climates remains to be confirmed. As described above, we used CESM1 simulations (EO3 and EO4 from Cramwinckel et al., 2018) to compare the “true” model simulation GMST to that calculated using Dcomb-1 (Supplement). We find that Dcomb-1 underestimates GMST by 1 C during the Eocene when the model high-latitude SST is used as a proxy for the deep ocean and 2–3 C when the model deep-ocean temperature is used. As such, we suggest that Dcomb-1 may reflect a minimum GMST constraint. We suggest that variable weighting of the deep ocean and tropics could improve the Dcomb method in future studies (Eq. 5 gives an equal weighting to each).

3.2 Influence of different proxy datasets upon Dsurf-derived GMST estimates

To explore the importance of the proxies themselves for Dsurf-derived GMST estimates, we conducted a series of illustrative subsampling experiments relative to Dsurf-baseline (Table 2). This was performed for methods Dsurf-1, Dsurf-2, Dsurf-3, and Dsurf-4. In the first subsampling experiment (Dsurf-Frosty; Table 2), we include δ18O SST estimates from recrystallised planktonic foraminifera. This yields lower GMST estimates (<1 to 4 C; e.g. Figs. S6–S8) and is consistent amongst all four methods. This agrees with previous studies which indicate that δ18O values from recrystallised planktonic foraminifera are significantly colder than estimates derived from the δ18O value of well-preserved foraminifera (Pearson et al., 2001; Sexton et al., 2006; Edgar et al., 2015), foraminiferal Mg∕Ca ratios (Creech et al., 2010; Hollis et al., 2012), and clumped isotope values from larger benthic foraminifera (Evans et al., 2018).

The removal of TEX86 results in lower GMST estimates (∼1 to 4 C; e.g. Figs. S6–S8) across all methodologies (Dsurf-NoTEX; Table 2). This is consistent with previous studies which indicate that TEX86 gives slightly higher SSTs than other proxies, especially in the middle to high latitudes (e.g. Hollis et al., 2012; Inglis et al., 2015). The functional response of TEX86 at higher-than-modern SSTs remains relatively uncertain, which may explain why TEX86 gives slightly higher SSTs than other proxies (see discussion in Hollis et al., 2019). New indices or calibrations could help to reduce the uncertainty associated with TEX86-derived SST estimates beyond the modern calibration range. TEX86 values can also be complicated by the input of isoGDGTs from other sources (see discussion in Hollis et al., 2019). The DeepMIP database excludes samples with anomalous GDGT distributions (Hollis et al., 2019). However, a Gaussian process regression (GPR) model may help to better identify anomalous GDGT distributions in the sedimentary record using a nearest-neighbour distance metric (Eley et al., 2019). This methodology could be employed in future studies to further refine GDGT-based SST datasets, but this methodology is currently under review and is not considered here. Despite the caveats and concerns raised in previous work, the exclusion of TEX86 data shifts GMST by a relatively small amount.

The input of brGDGTs from archives other than mineral soils or peat can bias LAT estimates towards lower values (Inglis et al., 2017; Hollis et al., 2019), and the exclusion of MBT(')  CBT-derived LAT estimates could yield higher GMST values. Excluding MBT(')  CBT in marine sediments does yield slightly warmer GMST estimates (0.5 to 1.0 C). However, the impact of excluding MBT(')  CBT values is relatively minor because there are other proxies (e.g. pollen assemblages, leaf floral) which yield comparable LAT estimates in the regions where MBT(')  CBT values are removed (e.g. the SW Pacific).

The removal of δ18O values from paleosols and mammals as well as paleosol climofunctions (Dsurf-NoPaleosol; Table 2) also leads to slightly warmer GMST estimates (∼0.5C). This may be related to additional controls on paleosol and mammal δ18O values. This includes variations in the isotopic composition of rainfall (i.e. meteoric δ18O; Hyland and Sheldon, 2013), variations in soil water δ18O values (Hyland and Sheldon, 2013), and/or δ18O heterogeneity within nodules (e.g. Dworkin et al., 2005). Temperature estimates from paleosol climofunctions may also be prone to underestimation (e.g. Sheldon, 2009), and Hyland and Sheldon (2013) suggest that paleosol climofunctions are only applied as an indicator of relative temperature change. Intriguingly, the Dsurf-1 method yields much higher GMST estimates during the EECO when δ18O values from paleosols and mammals as well as paleosol climofunctions are excluded (∼3C higher than Dsurf-baseline). This is attributed to the inclusion of two “cold” LAT estimates from the Salta Basin, NW Argentina (Hyland et al., 2017), which overly influence GMST (e.g. Fig. 2). For Dsurf-1, a direct comparison of new and prior estimates (Caballero and Huber, 2013) can be made in which the only change has been the use of a newer data compilation. For our new estimate, the EECO is ∼4.5C colder than previous estimates (29.75 C; Caballero and Huber, 2013). Given that the floristic LAT estimates are identical between the DeepMIP compilation and the older compilation, the lower GMST estimates are largely due to the incorporation of additional LAT datasets (e.g. paleosol climofunctions).

3.3 A combined estimate of GMST during the DeepMIP target intervals

To derive a combined estimate of GMST during the latest Paleocene, PETM, and EECO, we employ a probabilistic approach, using Monte Carlo resampling with full propagation of errors. Our combined estimate employs GMST estimates from each baseline experiment (except Dsurf-1 for the EECO for which we use Dsurf-NoPaleosol; see discussion above). We generated 1 000 000 iterations for each of the six methods for each time interval (latest Paleocene, PETM, and EECO). In these iterations, the GMST estimates were randomly sampled with replacement within their full uncertainty envelopes, assuming a Gaussian distribution of errors. As the different GMST estimates ultimately derive from the same proxy dataset, we do not consider them to be independent. The resulting 6 000 000 GMST iterations for each time period are thus simply added into a single probability density function in order to fully represent uncertainty (Fig. 7). This is equivalent to a linear pooling approach with equal weights (Genest and Zidek, 1986). From this probability distribution, the median value and the upper and lower limits corresponding to 66 % and 90 % confidence limits were identified (Table 4).

Figure 7Probability density function of combined GMST during the DeepMIP intervals with full propagation of errors. GMST estimates are calculated on the mantle-based reference frame.


Table 4Combined GMST estimates (66 % and 90 % confidence intervals) during the (i) latest Paleocene (LP), (ii) PETM, and (iii) EECO.

Download Print Version | Download XLSX

Sequential removal of one GMST method at a time (jackknife resampling) was performed to examine the influence of a single method upon the average GMST estimate. Jackknifing reveals that no single method overly influences the mean GMST or 66 % confidence intervals during the latest Paleocene, PETM, or EECO (±1.5C; Supplement and Fig. S9). However, the removal of Dsurf-2 (which has relatively large error bars; Fig. 6) reduces the 90 % confidence interval (Supplement). We also show that removing Dcomb-1 removes the bimodality of the temperature distribution (Fig. S9). This is because Dcomb-1 is associated with consistently lower GMST estimates compared to other methods (see Sect. 3.1).

During the latest Paleocene, the average GMST estimate is 26.3 C and ranges between 22.3 and 28.3 C (66 % confidence interval; Table 4; Fig. 7). During the PETM, the average GMST is higher (31.6 C) and ranges between 27.2 and 34.5 C (66 % confidence interval; Table 4; Fig. 7). Assuming a pre-industrial GMST of 14 C, our average GMST estimates indicate that the latest Paleocene and PETM are 12.3 and 17.6 C warmer than pre-industrial, respectively. Our results indicate that GMST likely increased by ∼4 to 6 C between the latest Paleocene and PETM (66 % confidence), in keeping with previous estimates (Frieling et al., 2019; Dunkley Jones, 2013). During the EECO, the average GMST estimate is 27.0 C and likely ranges between 23.2 and 29.7 C (66 % confidence interval; Table 4; Fig. 7). Assuming a pre-industrial GMST of 14 C, our average GMST estimate indicates that the EECO is 13.0 C warmer than pre-industrial. The GMST anomaly for the EECO is ∼2C lower than previous studies (∼15C warmer than pre-industrial; Caballero and Huber, 2013; Zhu et al., 2019), but the median falls within the range quoted previously in the IPCC AR5 (9 to 14 C warmer than pre-industrial). The EECO is approximately 4 to 5 C colder than the PETM (66 % confidence). This is larger than previously suggested (∼3C; Zhu et al., 2019) and may be related to a cold bias in EECO GMST estimates (see Sect. 3.1).

3.4 Equilibrium climate sensitivity during the latest Paleocene, PETM, and EECO

Equilibrium climate sensitivity (ECS) can be defined as the equilibrium change in global near-surface air temperature resulting from a doubling in atmospheric CO2. Various “flavours” of ECS exist, some of which specifically exclude various feedback processes not always included in climate models, such as those associated with ice sheets, vegetation, or aerosols (Rohling et al., 2012). ECS may also be state-dependent (Caballero and Huber, 2013), and there is no reason to expect that it has not changed with time or as a function of climate state (Farnsworth et al., 2019; Zhu et al., 2020). Therefore, direct comparison of ECS in the past to modern conditions is a fraught enterprise. For our purposes we define a bulk ECS (ECSbulk) as being a gross estimate of ECS between our three intervals and pre-industrial, i.e.

(6) ECS bulk = ( Δ T CO 2 - vs - PI ) / ( Δ F CO 2 - vs - PI ) ,

where ΔTCO2-vs-PI is the temperature difference between pre-industrial and the time period of interest that can be attributed to CO2 forcing, and ΔFCO2-vs-PI is the CO2 forcing relative to pre-industrial. The result is then normalised to a CO2 forcing equal to a doubling of CO2. Such calculations have been performed previously (e.g. Anagnostou et al., 2016) and they provide some constraint on the range of climate sensitivity values that are relevant for near-modern prediction (Rohling et al., 2012). For example, Anagnostou et al. (2016) indicated that early Eocene ECS (excluding ice sheet feedbacks) falls within the range 2.1–4.6 C per CO2 doubling, with maximum probability for the EECO of 3.8 C. These values (2.1–4.6 C per CO2 doubling) are similar to the IPCC ECS range (1.5–4.5 C at 66 % confidence). Here we calculate bulk ECS estimates using the change in GMST and CO2 in the latest Paleocene, PETM, and EECO intervals with reference to the pre-industrial. Following the approach of Anagnostou et al. (2016) and using the forcing equation of Byrne and Goldblatt (2014), we first determine the relative change in climate forcing relative to pre-industrial (ΔFCO2-vs-PI):

(7) Δ F CO 2 - vs - PI = 5.32 ln ( C t / C PI ) + 0.39 [ ln ( C t / C PI ) ] 2 ,

where CPI is the atmospheric CO2 concentration during pre-industrial (278 ppm) and Ct refers to the CO2 reconstruction at a particular time in the Eocene. The mean proxy estimate of CO2 for the PETM is ∼2200 ppmv (+1904/-699 ppmv; 95 % confidence) (Gutjahr et al., 2017). The mean proxy estimate of CO2 for the LP is ∼870 ppmv (Gutjahr et al., 2017). The uncertainty of latest Paleocene CO2 represents 2 standard deviations of pre-PETM CO2 (Gutjahr et al., 2017), equal to ±400 ppm. The mean proxy estimate of CO2 for the EECO is ∼1625 ppmv (±750 ppmv; 95 % confidence) (Anagnostou et al., 2016; Hollis et al., 2019). To calculate bulk ECS, we then use radiative forcing from a doubling of CO2 from Byrne and Goldblatt (2014) to translate CO2 into forcing relative to pre-industrial (ΔFCO2):

(8) ECS = ( Δ T CO 2 - vs - PI ) / Δ F CO 2 - vs - PI × 3.875 ,

where GMST (ΔT) distributions are based on output generated via our Monte Carlo simulations (see Sect. 3.3). Some of the temperature anomaly of the latest Paleocene, PETM, and EECO is caused not by CO2 but by the different paleotopography, paleobathymetry, and solar constant compared with pre-industrial. Furthermore, we choose here to calculate an ECS that explicitly excludes feedbacks associated with vegetation, ice sheets, and aerosols, i.e. S[CO2,LI,VG,AE] in the nomenclature of Rohling et al. (2012). To account for these effects, we subtract a value of 4.5 C (Caballero and Huber, 2013; Zhu et al. 2019) from GMST; i.e.

(9) Δ T CO 2 - vs - PI = Δ GMST - 4.5 C .

Following Anagnostou et al. (2016), the uncertainty of the slow-feedback correction on ΔGMST follows a uniform “flat” probability (±1.5C). This value of 4.5 C is based upon a comparison of pre-industrial and Eocene simulations (both CO2) conducted with CESM1.2 (Zhu et al., 2019), which incorporates the paleogeographic, solar constant, ice sheet, vegetation, aerosol, and ice sheet changes from pre-industrial to Eocene. Our value is similar to previous studies which attribute ∼4 to 6 C to non-CO2 and non-aerosol forcings and feedbacks (Anagnostou et al., 2016; Caballero and Huber, 2013, Lunt et al., 2012). However, the sensitivity to these Eocene boundary conditions is likely model-dependent and this value may differ between model simulations. The uncertainties in our estimated ECS are the products of 10 000 realisations of the latest Paleocene, PETM, and EECO CO2 values as well as the respective ΔGMST estimate (the mean estimate and propagated uncertainty) based on randomly sampling each variable within its 66 % and 90 % confidence interval uncertainty envelope.

Table 5Estimates of ECS (66 % and 90 % confidence) during the (i) latest Paleocene (LP), (ii) PETM, and (iii) EECO.

Download Print Version | Download XLSX

S[CO2,LI,VG,AE] values (66% confidence) for the EECO and PETM average 0.80 (0.46 to 1.15) and 0.92 (0.60 to 1.20), respectively. This yields ECS estimates (66 % confidence) for the EECO and PETM compared to modern which average 3.1 C (1.8 to 4.4 C) and 3.6 C (2.3 to 4.7 C), respectively (Table 5; Fig. 8). These are broadly comparable to previous estimates from the early Eocene which account for paleogeography and other feedbacks (∼2.1 to 4.6 C; Anagnostou et al., 2016). They are also similar to those predicted by the IPCC (1.5 to 4.5 C per doubling CO2). S[CO2,LI,VG,AE] values (66 % confidence) during the latest Paleocene average 1.16 (0.61 to 1.75), which is somewhat higher than the other DeepMIP intervals. This yields ECS estimates (66 % confidence) for the latest Paleocene which average 4.5 C (2.4 to 6.8 C) (Table 5; Fig. 8). Higher ECS values are attributed to relatively high GMST estimates (∼26C) and relatively low CO2 values (∼870 ppm) during the latest Paleocene. As latest Paleocene CO2 estimates remain highly uncertain (Gutjahr et al., 2017; see above), new high-fidelity CO2 records are required to accurately constrain ECS during this time.

Figure 8Probability density function of bulk ECS during the latest Paleocene, PETM, and EECO that explicitly accounts for non-CO2 forcings of paleogeography and the solar constant, as well as feedbacks associated with land ice, vegetation, and aerosols (Zhu et al., 2019), i.e. S[CO2,LI,VG,AE] in the nomenclature of Rohling et al. (2012).

ECS may be strongly state-dependent, and model simulations indicate a non-linear increase in ECS at higher temperatures (Caballero and Huber, 2013; Zhu et al., 2019) due to changes in cloud feedbacks (Abbot et al., 2009; Caballero and Huber, 2010; Arnold et al., 2012; Zhu et al., 2019). This implies caution when relating geological estimates to modern climate predictions (e.g. Rohling et al., 2012; Zhu et al., 2020) and it may be more appropriate to calculate ECS between different time intervals (e.g. latest Paleocene to PETM; Shaffer et al., 2016). To this end, we also calculate ECS between the transition from the latest Paleocene to the PETM, assuming that non-CO2 forcings and feedbacks are negligible. This yields an ECS estimate of 3.6 C. However, early Paleogene CO2 estimates remain uncertain (Gutjahr et al., 2017), and well-synchronised, continuous, and high-resolution CO2 records are required to accurately constrain ECS during the DeepMIP intervals.

4 Conclusions

Using six different methods, we have quantified global mean surface temperatures (GMSTs) during the latest Paleocene, PETM, and EECO. GMST was calculated within a coordinated, experimental framework and utilised six methodologies including three different input datasets. After evaluating the impact of different proxy datasets upon GMST estimates, we combined all six methodologies to derive an average GMST value during the latest Paleocene, PETM, and EECO. We show that the “average” GMST estimate (66 % confidence) during the latest Paleocene, PETM, and EECO is 26.3 C (22.3 to 28.3 C), 31.6 C (27.2 to 34.5 C), and 27.0 C (23.2 to 29.7 C), respectively. Assuming a pre-industrial GMST of 14 C, the latest Paleocene, PETM, and EECO are 12.3 C, 17.6, and 13.0 C warmer than modern, respectively. Using our “combined” GMST estimate, we demonstrate that “bulk” ECS (66 % confidence) during the latest Paleocene, PETM, and EECO is 4.5 C (2.4 to 6.8 C), 3.6 C (2.3 to 4.7 C), and 3.1 C (1.8 to 4.4 C) per doubling of CO2. Taken together, this study improves our characterisation of the global mean temperature of these key time intervals, allowing future climate change to be put into the context of past changes and allowing us to provide a refined estimate of ECS.

Data availability

Data can be accessed via the Supplement or from the authors (contact email:


The supplement related to this article is available online at:

Author contributions

Authorship of this paper is organised into three tiers according to the contributions of each individual author. GNI (Tier I) organised the structure and writing of the paper, contributed to all sections of the text, and designed the figures. FB, NJB, MJC, DE, GLF, MH, DJL, NS, SS, JET, and RW (Tier II authors) assumed a leading role by contributing the methodologies used in the text. EA, AMdB, TDJ, KE, CJH, DKH, and RDP (Tier III authors) contributed intellectually to the text and figure design.

Competing interests

The authors declare that they have no conflict of interest.


We thank two anonymous reviewers whose thoughtful comments significantly improved the paper. This research was funded by NERC through NE/P01903X/1 and NE/N006828/1, both of which supported Gordon N. Inglis, Daniel J. Lunt, Sebastian Steinig, and Richard D. Pancost. Gordon N. Inglis was also supported by a GCRF Royal Society Dorothy Hodgkin Fellowship. Natalie J. Burls is supported by NSF AGS-1844380 and the Alfred P. Sloan Foundation as a Research Fellow. Fran Bragg, Daniel J. Lunt, and Richard Wilkinson were funded by the EPSRC Past Earth Network (EP/M008363/1). Matthew Huber was funded by NSF OPP 1842059. Tom Dunkley Jones, Kirsty M. Edgar, and Gavin L. Foster were supported by NERC grant NE/P013112/1. Agatha De Boer and David Hutchinson acknowledge support from the Swedish Research Council under project 2016-03912. GFDL numerical simulations were performed using resources provided by the Swedish National Infrastructure for Computing (SNIC) at NSC, Linköping. David Hutchinson was also supported by FORMAS project 2018-01621. The authors also thank Chris Poulsen and Jiang Zhu for assistance with the CESM1.2 model simulations.

Financial support

This research has been supported by the Natural Environment Research Council (grant nos. NE/P01903X/1 and NE/N006828/1), the National Science Foundation (grant nos. AGS-1844380 and OPP 1842059), the EPSRC Past Earth Network (EP/M008363/1), NERC (NE/P013112/1), the Swedish Research Council (grant nos. 2016-03912 and 2018-01621), and the Royal Society (grant no. DHF\R1\191178).

Review statement

This paper was edited by Yannick Donnadieu and reviewed by two anonymous referees.


Abbot, D. S., Huber, M., Bousquet, G., and Walker, C. C.: High-CO2 cloud radiative forcing feedback over both land and ocean in a global climate model, Geophys. Res. Lett., 36, L05702,, 2009. 

Anagnostou, E., John, E. H., Edgar, K. M., Foster, G. L., Ridgwell, A., Inglis, G. N., Pancost, R. D., Lunt, D. J., and Pearson, P. N.: Changing atmospheric CO2 concentration was the primary driver of early Cenozoic climate, Nature, 533, 380–384,, 2016. 

Arnold, N. P., Tziperman, E., and Farrell, B.: Abrupt transition to strong superrotation driven by equatorial wave resonance in an idealized GCM, J. Atmos. Sci., 69, 626–640, 2012. 

Barnet, J. S., Littler, K., Westerhold, T., Kroon, D., Leng, M. J., Bailey, I., Röhl, U., and Zachos, J. C.: A high-Fidelity benthic stable isotope record of late Cretaceous–early Eocene climate change and carbon-cycling, Paleoceanogr. Paleoclim., 34, 672–691, 2019. 

Bijl, P. K., Schouten, S., Sluijs, A., Reichart, G.-J., Zachos, J. C., and Brinkhuis, H.: Early Palaeogene temperature evolution of the southwest Pacific Ocean, Nature, 461, 776–779, 2009. 

Bragg, F. J., Paine, P., Saul, A., Lunt, D. J., Wilkinson, R., and Zammit-Mangion, A.: A Statistical Algorithm for Evaluating Palaeoclimate Simulations Against Geological Observations, Geoscientific Model Development, in preperation, 2020. 

Byrne, B. and Goldblatt, C.: Radiative forcing at high concentrations of well-mixed greenhouse gases, Geophys. Res. Lett., 41, 152–160, 2014. 

Caballero, R. and Huber, M.: Spontaneous transition to superrotation in warm climates simulated by CAM3, Geophys. Res. Lett., 37, L11701,, 2010. 

Caballero, R. and Huber, M.: State-dependent climate sensitivity in past warm climates and its implications for future climate projections, P. Natl. Acad. Sci. USA, 110, 14162–14167, 2013. 

Cramer, B. S., Toggweiler, J. R., Wright, J. D., Katz, M. E., and Miller, K. G.: Ocean overturning since the Late Cretaceous: Inferences from a new benthic foraminiferal isotope compilation, Paleoceanogr. Paleoclim., 24, PA4216,, 2009. 

Cramwinckel, M. J., Huber, M., Kocken, I. J., Agnini, C., Bijl, P. K., Bohaty, S. M., Frieling, J., Goldner, A., Hilgen, F. J., Kip, E. L., Peterse, F., van der Ploeg, R., Rohl, U., Schouten, S., and Sluijs, A.: Synchronous tropical and polar temperature evolution in the Eocene, Nature, 559, 382–386, 2018. 

Creech, J. B., Baker, J. A., Hollis, C. J., Morgans, H. E., and Smith, E. G.: Eocene sea temperatures for the mid-latitude southwest Pacific from Mg/Ca ratios in planktonic and benthic foraminifera, Earth Planet. Sci. Lett., 299, 483–495, 2010. 

Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Hólm, E. V., Isaksen, L., Kållberg, P., Köhler, M., Matricardi, M., McNally, A. P., Monge-Sanz, B. M., Morcrette, J.-J., Park, B.-K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J.-N., and Vitart, F.: The ERA-Interim reanalysis: configuration and performance of the data assimilation system, Q. J. Roy. Meteorol. Soc., 137, 553–597, 2011. 

Dunkley Jones, T., Lunt, D. J., Schmidt, D. N., Ridgwell, A., Sluijs, A., Valdes, P. J., and Maslin, M.: Climate model and proxy data constraints on ocean warming across the Paleocene–Eocene Thermal Maximum, Earth Sci. Rev., 125, 123–145, 2013. 

Dworkin, S. I., Nordt, L., and Atchley, S.: Determining terrestrial paleotemperatures using the oxygen isotopic composition of pedogenic carbonate, Earth Planet. Sci. Lett., 237, 56–68, 2005. 

Edgar, K. M., Anagnostou, E., Pearson, P. N., and Foster, G. L: Assessing the impact of diagenesis on δ11B, δ13C, δ18O, Sr/Ca and B/Ca values in fossil planktic foraminiferal calcite, Geochim. Cosmochim. Acta, 166, 89–209, 2015. 

Edgar, K. M., Pälike, H., and Wilson, P. A.: Testing the impact of diagenesis on the δ18O and δ13C of benthic foraminiferal calcite from a sediment burial depth transect in the equatorial Pacific, Paleoceanography, 28, 468–480, 2013. 

Eley, Y. L., Thompson, W., Greene, S. E., Mandel, I., Edgar, K., Bendle, J. A., and Dunkley Jones, T.: OPTiMAL: A new machine learning approach for GDGT-based palaeothermometry, Clim. Past Discuss.,, in review, 2019. 

Evans, D., Sagoo, N., Renema, W., Cotton, L. J., Müller, W., Todd, J. A., Saraswati, P. K., Stassen, P., Ziegler, M., Pearson, P. N., Valdes, P. J., and Affek, H. P.: Eocene greenhouse climate revealed by coupled clumped isotope-Mg/Ca thermometry, P. Natl. Acad. Sci. USA, 115, 1174–1179,, 2018. 

Farnsworth, A., Lunt, D., O'Brien, C., Foster, G., Inglis, G., Markwick, P., Pancost, R., and Robinson, S.: Climate sensitivity on geological timescales controlled by non-linear feedbacks and ocean circulation, Geophys. Res. Lett., 46, 9880–9889,, 2019. 

Frieling, J., Gebhardt, H., Huber, M., Adekeye, O. A., Akande, S. O., Reichart, G. J., Middelburg, J. J., Schouten, S., and Sluijs, A.: Extreme warmth and heat-stressed plankton in the tropics during the Paleocene-Eocene Thermal Maximum, Sci. Adv., 3, p.e1600891,, 2017. 

Genest, C. and Zidek, J. V.: Combining probability distributions: a critique and an annotated bibliography, Stat. Sci., 1, 147–148, 1986. 

Goldner, A., Herold, N., and Huber, M.: Antarctic glaciation caused ocean circulation changes at the Eocene – Oligocene transition, Nature, 511, 574–577, 2014. 

Gutjahr, M., Ridgwell, A., Sexton, P. F., Anagnostou, E., Pearson, P. N., Pälike, H., Norris, R. D., Thomas, E., and Foster, G. L.: Very large release of mostly volcanic carbon during the Palaeocene–Eocene Thermal Maximum, Nature, 548, 573–577,, 2017. 

Hansen, J., Sato, M., Russell, G., and Kharecha, P.: Climate sensitivity, sea level and atmospheric carbon dioxide, Philos. Trans. Royal Soc. A, 371, 20120294, 2013. 

Hegerl, G. C., Zwiers, F. W., Braconnot, P., Gillett, N. P., Luo, Y., Marengo Orsini, J., Nicholls, N., Penner, J. E., and Stott, P. A.: Understanding and attributing climate change, IPCC, 2007: Climate Change 2007: the physical science basis. contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change, 2007. 

Herold, N., Buzan, J., Seton, M., Goldner, A., Green, J. A. M., Müller, R. D., Markwick, P., and Huber, M.: A suite of early Eocene (∼55 Ma) climate model boundary conditions, Geosci. Model Dev., 7, 2077–2090,, 2014. 

Hollis, C. J., Taylor, K. W. R., Handley, L., Pancost, R. D., Huber, M., Creech, J. B., Hines, B. R., Crouch, E. M., Morgans, H. E. G., Crampton, J. S., Gibbs, S., Pearson, P. N., and Zachos, J. C.: Early Paleogene temperature history of the Southwest Pacific Ocean: Reconciling proxies and models, Earth Planet. Sci. Lett., 349–350, 53–66,, 2012. 

Hollis, C. J., Dunkley Jones, T., Anagnostou, E., Bijl, P. K., Cramwinckel, M. J., Cui, Y., Dickens, G. R., Edgar, K. M., Eley, Y., Evans, D., Foster, G. L., Frieling, J., Inglis, G. N., Kennedy, E. M., Kozdon, R., Lauretano, V., Lear, C. H., Littler, K., Lourens, L., Meckler, A. N., Naafs, B. D. A., Pälike, H., Pancost, R. D., Pearson, P. N., Röhl, U., Royer, D. L., Salzmann, U., Schubert, B. A., Seebeck, H., Sluijs, A., Speijer, R. P., Stassen, P., Tierney, J., Tripati, A., Wade, B., Westerhold, T., Witkowski, C., Zachos, J. C., Zhang, Y. G., Huber, M., and Lunt, D. J.: The DeepMIP contribution to PMIP4: methodologies for selection, compilation and analysis of latest Paleocene and early Eocene climate proxy data, incorporating version 0.1 of the DeepMIP database, Geosci. Model Dev., 12, 3149–3206,, 2019. 

Huber, M. and Caballero, R.: The early Eocene equable climate problem revisited, Clim. Past, 7, 603–633,, 2011. 

Hutchinson, D. K., de Boer, A. M., Coxall, H. K., Caballero, R., Nilsson, J., and Baatsen, M.: Climate sensitivity and meridional overturning circulation in the late Eocene using GFDL CM2.1, Clim. Past, 14, 789–810,, 2018. 

Hyland, E. G. and Sheldon, N. D.: Coupled CO2-climate response during the early Eocene climatic optimum, Palaeogeogr. Palaeoclim. Palaeoecol., 369, 125–135, 2013. 

Hyland, E. G., Sheldon, N. D., and Cotton, J. M.: Constraining the early Eocene climatic optimum: A terrestrial interhemispheric comparison, GSA Bull., 129, 244–252,, 2017. 

Inglis, G. N., Farnsworth, A., Lunt, D., Foster, G. L., Hollis, C. J., Pagani, M., Jardine, P. E., Pearson, P. N., Markwick, P., Galsworthy, A. M. J., Raynham, L., Taylor, K. W. R., and Pancost, R. D.: Descent toward the Icehouse: Eocene sea surface cooling inferred from GDGT distributions, Paleoceanography, 30, 1000–1020,, 2015. 

Inglis, G. N., Collinson, M. E., Riegel, W., Wilde, V., Farnsworth, A., Lunt, D. J., Valdes, P., Robson, B. E., Scott, A. C., Lenz, O. K., Naafs, B. D. A., and Pancost, R. D.: Mid-latitude continental temperatures through the early Eocene in western Europe, Earth Planet. Sci. Lett., 460, 86–96,, 2017. 

IPCC: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P. M., Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 1535 pp., 2013. 

Keating-Bitonti, C. R., Ivany, L. C., Affek, H. P., Douglas, P., and Samson, S. D.: Warm, not super-hot, temperatures in the early Eocene subtropics, Geology, 39, 771–774, 2011. 

Knutti, R., Rugenstein, M. A., and Hegerl, G. C.: Beyond equilibrium climate sensitivity, Nat. Geosci., 10, 727–736, 2017. 

Lauretano, V., Littler, K., Polling, M., Zachos, J. C., and Lourens, L. J.: Frequency, magnitude and character of hyperthermal events at the onset of the Early Eocene Climatic Optimum, Clim. Past, 11, 1313–1324,, 2015. 

Littler, K., Röhl, U., Westerhold, T., and Zachos, J. C.: A high-resolution benthic stable-isotope record for the South Atlantic: Implications for orbital-scale changes in Late Paleocene–Early Eocene climate and carbon cycling, Earth Planet. Sci. Lett., 401, 18–30, 2014. 

Lunt, D. J., Dunkley Jones, T., Heinemann, M., Huber, M., LeGrande, A., Winguth, A., Loptson, C., Marotzke, J., Roberts, C. D., Tindall, J., Valdes, P., and Winguth, C.: A model–data comparison for a multi-model ensemble of early Eocene atmosphere–ocean simulations: EoMIP, Clim. Past, 8, 1717–1736,, 2012. 

Lunt, D. J., Farnsworth, A., Loptson, C., Foster, G. L., Markwick, P., O'Brien, C. L., Pancost, R. D., Robinson, S. A., and Wrobel, N.: Palaeogeographic controls on climate and proxy interpretation, Clim. Past, 12, 1181–1198,, 2016. 

Lunt, D. J., Huber, M., Anagnostou, E., Baatsen, M. L. J., Caballero, R., DeConto, R., Dijkstra, H. A., Donnadieu, Y., Evans, D., Feng, R., Foster, G. L., Gasson, E., von der Heydt, A. S., Hollis, C. J., Inglis, G. N., Jones, S. M., Kiehl, J., Kirtland Turner, S., Korty, R. L., Kozdon, R., Krishnan, S., Ladant, J.-B., Langebroek, P., Lear, C. H., LeGrande, A. N., Littler, K., Markwick, P., Otto-Bliesner, B., Pearson, P., Poulsen, C. J., Salzmann, U., Shields, C., Snell, K., Stärz, M., Super, J., Tabor, C., Tierney, J. E., Tourte, G. J. L., Tripati, A., Upchurch, G. R., Wade, B. S., Wing, S. L., Winguth, A. M. E., Wright, N. M., Zachos, J. C., and Zeebe, R. E.: The DeepMIP contribution to PMIP4: experimental design for model simulations of the EECO, PETM, and pre-PETM (version 1.0), Geosci. Model Dev., 10, 889–901,, 2017. 

Lunt, D. J., Bragg, F., Chan, W.-L., Hutchinson, D. K., Ladant, J.-B., Niezgodzki, I., Steinig, S., Zhang, Z., Zhu, J., Abe-Ouchi, A., de Boer, A. M., Coxall, H. K., Donnadieu, Y., Knorr, G., Langebroek, P. M., Lohmann, G., Poulsen, C. J., Sepulchre, P., Tierney, J., Valdes, P. J., Dunkley Jones, T., Hollis, C. J., Huber, M., and Otto-Bliesner, B. L.: DeepMIP: Model intercomparison of early Eocene climatic optimum (EECO) large-scale climate features and comparison with proxy data, Clim. Past Discuss.,, in review, 2020. 

Marchitto, T., Curry, W., Lynch-Stieglitz, J., Bryan, S., Cobb, K., and Lund, D.: Improved oxygen isotope temperature calibrations for cosmopolitan benthic foraminifera, Geochim. Cosmochim. Acta, 130, 1–11, 2014. 

Masson-Delmotte, V., Schulz, M., Abe-Ouchi, A., Beer, J., Ganopolski, A., Gonzalez Rouco, J. F., Jansen, E., Lambeck, K., Luterbacher, J., Naish, T., Osborn, T., Otto-Bliesner, B., Quinn, T., Ramesh, R., Rojas, M., Shao, X., and Timmermann, A.: Information from Paleoclimate Archives, in: Climate Change 2013 – The Physical Science Basis: Working Group I Contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, 383–464, 2014. 

Pearson, P. N., Ditchfield, P. W., Singano, J., Harcourt-Brown, K. G., Nicholas, C. J., Olsson, R. K., Shackleton, N. J., and Hall, M. A.: Warm tropical sea surface temperatures in the Late Cretaceous and Eocene epochs, Nature, 413, 481–487, 2001. 

Pearson, P. N., van Dongen, B. E., Nicholas, C. J., Pancost, R. D., Schouten, S., Singano, J. M., and Wade, B. S.: Stable warm tropical climate through the Eocene Epoch, Geology, 35, 211–214,, 2007. 

Rohling, E. J., Sluijs, A., Dijkstra, H. A., Köhler, P., van de Wal, R. S. W., von der Heydt, A. S., Beerling, D. J., Berger, A., Bijl, P. K., Crucifix, M., DeConto, R., Drijfhout, S. S., Fedorov, A., Foster, G. L., Ganopolski, A., Hansen, J., Hönisch, B., Hooghiemstra, H., Huber, M., Huybers, P., Knutti, R., Lea, D. W., Lourens, L. J., Lunt, D., Masson-Delmotte, V., Medina-Elizalde, M., Otto-Bliesner, B., Pagani, M., Pälike, H., Renssen, H., Royer, D. L., Siddall, M., Valdes, P., Zachos, J. C., Zeebe, R. E., and Members, P. P.: Making sense of palaeoclimate sensitivity, Nature, 491, 683–691,, 2012.  

Sexton, P. F., Wilson, P. A., and Pearson, P. N.: Microstructural and geochemical perspectives on planktic foraminiferal preservation: “Glassy” versus “Frosty”, Geochem. Geophys. Geosyst., 7, Q12P19,, 2006. 

Sexton, P. F., Norris, R. D., Wilson, P. A., Pälike, H., Westerhold, T., Röhl, U., Bolton, C. T., and Gibbs, S.: Eocene global warming events driven by ventilation of oceanic dissolved organic carbon, Nature, 471, 349–352, 2011. 

Shaffer, G., Huber, M., Rondanelli, R., and Pedersen, J. O. P. J. G. R. L.: Deep time evidence for climate sensitivity increase with warming, Geophys. Res. Lett., 43, 6538–6545, 2016. 

Sheldon, N. D.: Non-marine records of climatic change across the Eocene-Oligocene transition, The Late Eocene Earth – Hothouse, Icehouse, and Impacts, Geol. Soc. Am. Special Paper, 452, 241–248, 2009. 

Sijp, W. P., England, M. H., and Huber, M.: Effect of the deepening of the Tasman Gateway on the global ocean, Paleoceanography, 26, PA4207,, 2011. 

Stevens, B., Sherwood, S. C., Bony, S., and Webb, M. J.: Prospects for narrowing bounds on Earth's equilibrium climate sensitivity, Earth's Future, 4, 512–522, 2016. 

Tierney, J. E. and Tingley, M. P.: A Bayesian, spatially-varying calibration model for the TEX86 proxy, Geochim. Cosmochim. Acta, 127, 83–106, 2014. 

Westerhold, T., Röhl, U., Donner, B., and Zachos, J. C.: Global extent of early Eocene hyperthermal events: A new Pacific benthic foraminiferal isotope record from Shatsky Rise (ODP Site 1209), Paleoceanogr. Paleoclim., 33, 626–642, 2018. 

Zachos, J. C., Dickens, G. R., and Zeebe, R. E.: An early Cenozoic perspective on greenhouse warming and carbon-cycle dynamics, Nature, 451, 279–283, 2008. 

Zhu, J., Poulsen, C. J., and Tierney, J. E.: Simulation of Eocene extreme warmth and high climate sensitivity through cloud feedbacks, Sci. Adv., 5, EAAX1874, 2019. 

Zhu, J., Poulsen, C. J., and Otto-Bliesner, B. L.: High climate sensitivity in CMIP6 model not supported by paleoclimate, Nat. Clim. Change, 10, 378–379, 2020. 

Short summary
This paper presents estimates of global mean surface temperatures and climate sensitivity during the early Paleogene (∼57–48 Ma). We employ a multi-method experimental approach and show that i) global mean surface temperatures range between 27 and 32°C and that ii) estimates of bulk equilibrium climate sensitivity (∼3 to 4.5°C) fall within the range predicted by the IPCC AR5 Report. This work improves our understanding of two key climate metrics during the early Paleogene.