Research ArticleCLIMATOLOGY

Simulation of Eocene extreme warmth and high climate sensitivity through cloud feedbacks

See allHide authors and affiliations

Science Advances  18 Sep 2019:
Vol. 5, no. 9, eaax1874
DOI: 10.1126/sciadv.aax1874


The Early Eocene, a period of elevated atmospheric CO2 (>1000 ppmv), is considered an analog for future climate. Previous modeling attempts have been unable to reproduce major features of Eocene climate indicated by proxy data without substantial modification to the model physics. Here, we present simulations using a state-of-the-art climate model forced by proxy-estimated CO2 levels that capture the extreme surface warmth and reduced latitudinal temperature gradient of the Early Eocene and the warming of the Paleocene-Eocene Thermal Maximum. Our simulations exhibit increasing equilibrium climate sensitivity with warming and suggest an Eocene sensitivity of more than 6.6°C, much greater than the present-day value (4.2°C). This higher climate sensitivity is mainly attributable to the shortwave cloud feedback, which is linked primarily to cloud microphysical processes. Our findings highlight the role of small-scale cloud processes in determining large-scale climate changes and suggest a potential increase in climate sensitivity with future warming.


The Early Eocene [~56 to 48 million years (Ma) ago], the warmest interval in the Cenozoic (the past 66 Ma), was a period of elevated atmospheric CO2 concentrations [~1625 ± 760 parts per million by volume (ppmv)] (1, 2), substantially higher surface temperatures (3), and a latitudinal temperature gradient that was lower than today by at least 32% (4). The Eocene epoch began with the Paleocene-Eocene Thermal Maximum (PETM; ~56 Ma ago), the most severe of several short (105 to 106 years) hyperthermal events. During the PETM event, global surface temperature rose by 5° to 9°C (5, 6) in response to an estimated atmospheric CO2 increase of ~70 to 100%, suggesting an equilibrium climate sensitivity (ECS) of ~6°C (79). Simulating the extreme warmth of the Early Eocene and the large temperature increase in the PETM has been a challenge for climate models given their modest climate sensitivity (2.1° to 4.7°C) (10). The inability of climate models to match the warm conditions inferred from proxy evidence has been attributed to missing model components and physical processes, climate forcings, or misinterpretations and uncertainties in proxy reconstructions (1117).

Here, we conduct simulations of the Early Eocene using the Community Earth System Model version 1.2 (CESM1.2) with the Community Atmosphere Model version 5 (CAM5) (18) and compare them to proxy estimates of near-surface temperature. Following the Deep-Time Model Intercomparison Project (DeepMIP) protocol (2), four Eocene simulations are carried out with atmospheric CO2 concentrations of 1×, 3×, 6×, and 9× preindustrial levels (PILs; 285 ppmv). Other than CO2, the simulations use the same set of boundary conditions, including Early Eocene paleogeography, land-sea mask, vegetation distribution, and preindustrial (PI) non-CO2 greenhouse gas concentrations, soil properties, natural aerosol emissions, solar constant, and orbital parameters (Materials and Methods; figs. S1 and S2) (2). Our simulations capture major climatic features of the Early Eocene and the PETM in a state-of-the-art Earth system model forced by CO2 concentrations consistent with proxy reconstructions and without the addition of exotic forcings (16) or alteration of the model physics (19). That the CESM1.2 simulates the past Eocene extreme warmth within the uncertainty of proxy records provides additional confidence in this model’s simulation of high-CO2 future climates.


Model-data comparison of the Early Eocene temperature

Global mean surface temperature (GMST) in the Eocene 1× simulation is 4.5°C higher than in the PI simulation (Fig. 1A) due to the absence of ice sheets and anthropogenic aerosols, more widespread vegetation cover, and differences in the distribution of natural aerosols. Eocene GMST increases to 25.0°, 29.8°, and 35.5°C in the 3×, 6×, and 9× CO2 simulations, respectively. For comparison, we estimate the Early Eocene GMST from available terrestrial and marine proxy data to have been 29° ± 3°C (95% confidence interval; gray patch in Fig. 1A), 15°C greater than PI GMST and in close agreement with the latest estimate (Materials and Methods; table S1) (20). Interpolating our simulated GMSTs to the mean proxy estimate of CO2 for the Early Eocene (~1625 ppmv) yields a temperature of 29.4°C, a value in good agreement with proxy temperatures.

Fig. 1 Model-data comparison of the Early Eocene GMST and relative meridional SST gradient.

(A) GMST (in °C) in model simulations (markers) as a function of CO2 concentration compared with proxy estimates of the Early Eocene (gray patch) and the pre-PETM to PETM warming (black arrow). (B) Relative latitudinal gradient of SST in model simulations (markers) as a function of CO2 compared with estimates from proxies (gray patch). The CO2 ranges for the Early Eocene and PETM are from the latest estimates from proxies (1) and modeling (8), respectively. Eocene CESM1.2 simulations from this study are denoted as filled red circles. For reference, the PI simulation is marked as a black circle. Previous Eocene simulations (16, 19) including those from the Eocene Modeling Intercomparison Project (EoMIP) (3) are denoted by gray open markers. Note that methane concentrations in previous simulations that are different from the PI value have been converted to the equivalent CO2 concentration (Materials and Methods). In CCSM3KS simulations, the authors altered cloud properties substantially to account for presumed Eocene aerosol changes (16). Similarly, model physics are significantly altered in FAst Met Office/UK Universities Simulator (FAMOUS) (E17) (19).

Our simulations also capture the marked reduction in the meridional sea surface temperature (SST) gradient that is characteristic of the Early Eocene (4, 11). The 3×, 6×, and 9× simulations exhibit decreases in the SST gradient of 19, 24, and 32% relative to the PI case, all falling within the reconstructed range of 31 ± 13% for the Early Eocene (Fig. 1B; Materials and Methods). Simulated zonal mean surface temperatures over land and ocean compare well with proxy values at both high and low latitudes (Fig. 2). Our 3× and 9× simulations, which span the range of Early Eocene CO2 values (~3× to 8.5× PIL), bracket 80% of terrestrial and 50% of marine mean proxy temperatures. Site-by-site comparison of model and data paleotemperatures for the Early Eocene suggests a good model-data agreement with the root mean square error in simulations approaching the uncertainty in proxy records (fig. S3).

Fig. 2 Model simulated zonal mean temperature over land and ocean compared with proxy estimates.

(A) Zonal mean land surface temperature in the Eocene simulations compared with Early Eocene terrestrial proxy evidence (11). (B) Zonal mean SST in the Eocene simulations compared with published Early Eocene SST data, compiled in this study (table S1). Inferred temperatures using δ18O of planktic foraminifera, clumped isotopes, Mg/Ca of planktic foraminifera, and TEX86 are denoted as filled circles, upward-pointing triangles, squares, and downward-pointing triangles, respectively.

Our simulations also reproduce the PETM warming within proxy uncertainties of GMST and atmospheric CO2. From proxy records, we estimate GMST to have increased by ~5°C from ~27° (pre-PETM conditions) to ~32°C during the PETM (black arrow in Fig. 1A; table S2; Materials and Methods). PETM CO2 levels are not well constrained, and proxy and model estimates suggest an increase of approximately 70 to 100% over pre-PETM levels (79). Assuming a pre-PETM GMST of 27°C, we interpolate the corresponding CO2 change in our 3×, 6×, and 9× simulations and calculate a PETM warming of 4.6° to 6.8°C. Model simulations match the PETM warming in all 21 proxy records within their uncertainty, with a root mean square error of ~3°C (fig. S4).

The CESM1.2 with CAM5 shows marked improvement over earlier models (gray markers in Fig. 1) in its ability to simulate Eocene extreme warmth and the low meridional surface temperature gradient at CO2 levels that are within the range inferred from proxies (1). In Community Climate System Model version 3 (CCSM3), a predecessor of the CESM, atmospheric CO2 levels of more than 4560 ppmv were required to simulate an Eocene GMST of ~29°C (21). To address this limitation, past studies have simulated a warmer Eocene climate by introducing novel aerosol forcings (16) or by altering the model physics in ways that have not been validated against modern observational constraints (19). In contrast to these studies, the improvement in the CESM1.2 is due to the substantially larger climate sensitivity that is produced by self-consistent model physics under high CO2 levels. As we discuss below, this arises primarily from improved treatments of cloud microphysical processes in the atmosphere model CAM5 (22, 23).

Increased climate sensitivity with warming and cloud feedback

We estimate the ECS, defined as the GMST response to CO2 doubling, for the PI and the Eocene 1×, 3×, and 6× CO2 simulations using a slab ocean model (SOM) configuration (Materials and Methods). In the PI case, ECS is 4.2°C (Fig. 3A) and within the range of the latest estimates from the Intergovernmental Panel on Climate Change (IPCC) (10). ECS for the Eocene 1× case is 3.5°C, about 0.7°C lower than the PI ECS due primarily to a lower surface albedo feedback (Fig. 3B and table S3). Eocene ECS increases significantly with higher background CO2 levels to 6.6° and 9.7°C in the 3× and 6× simulations, respectively (Fig. 3A). In contrast, ECS in previous Eocene simulations (3) and in SOM simulations using CAM4 (Materials and Methods), the predecessor of CAM5, is 2.1° to 4.3°C and exhibits only minor increases (<1.5°C) with CO2 increase up to 8× PIL. Approximately 30% of the increase in ECS from CAM4 to CAM5 (from 3.2° to 4.2°C), under PI conditions, has been attributed to an increase in the efficacy of CO2 forcing (the radiative forcing per doubling CO2) as a result of the updated radiation scheme in CAM5 (24). Here, from offline radiation calculations, we estimate the efficacy of CO2 forcing to be 3.8 W m−2 for the PI simulation and 4.0, 4.7, and 5.2 W m−2 for the Eocene 1×, 3×, and 6× experiments using CAM5, respectively. A similar increase in efficacy of CO2 forcing has previously been found and attributed to the rapid nonfeedback cloud adjustments to radiative forcing and nonlogarithmic CO2 opacity (21). The increase in CO2 efficacy (by 0.7 and 1.2 W m−2) in the Eocene 3× and 6× runs accounts for ~20% (0.6° and 1.1°C) of the total ECS increase, assuming a linear relationship between radiative forcing and GMST. These results suggest that physical processes (climate feedbacks) other than the CO2 radiative forcing explain the high ECS in our Eocene simulations.

Fig. 3 ECS and climate feedback parameters in the Eocene simulations.

(A) The ECS (in °C) in the PI (black circle) and Eocene 1×, 3×, and 6× simulations (red filled circles) as a function of the atmospheric CO2 concentration. For comparison, ECSs in previous Eocene simulations (3) and in CAM4 slab ocean simulations are shown as gray open markers and blue filled circles, respectively. ECS for Eocene 6× is shown as a smaller marker than the other Eocene simulations, as it is estimated from the 6× and 9× experiments instead of using slab ocean simulations (Materials and Methods). (B) Climate feedback parameters (in W m−2 K−1) diagnosed from the two-way PRP method as a function of the atmospheric CO2 concentration in the Eocene simulations. Surface albedo (ALB), cloud (CLD), lapse rate (LPR), Planck (PLK), and water vapor (WVP) feedback parameters are shown. The cloud feedback parameter is further decomposed into shortwave and longwave components. For reference, feedback parameters in the PI simulation are denoted as open markers. Detailed numbers and uncertainty are shown in table S3. (C) Contributions to the shortwave cloud feedback parameter from cloud fraction, scattering, and absorption properties calculated from the approximated PRP method using the Eocene 1× and 9× simulations.

To quantify the strength of climate feedbacks responsible for the increased ECS with warming in our Eocene simulation, we computed the climate feedback parameter, defined as the increase in net downward radiation at top-of-the-atmosphere per unit of global warming, using a two-way partial radiative perturbation (PRP) calculation (25) (Materials and Methods; Fig. 3B and table S3). In the Eocene 1× simulation, surface albedo, cloud, lapse rate, Planck, and water vapor feedback parameters are 0.16, 0.50, −0.60, −3.03, and 1.97 W m−2 K−1, respectively, with a doubling of CO2. In the Eocene 6× simulation, the cloud feedback parameter increases by 110%, predominantly attributable to its shortwave component, and the water vapor feedback parameter rises by 67%. Meanwhile, albedo and lapse rate feedback parameters decrease, partly offsetting the increase in water vapor and cloud feedbacks. The intensification of the water vapor feedback results mainly from an increase in saturation vapor pressure and an increase in the height of the tropopause with warming (fig. S5) (26), responses that are consistent with previous Eocene simulations using CCSM3 (21). The doubling of the cloud shortwave feedback parameter from 1× to 6× PIL CO2 is, however, a fundamentally different response than in CCSM3 (21) and CAM4 SOM simulations, which exhibit insignificant or decreasing trends with warming up to 6× PIL CO2 (fig. S6).

To further assess the importance of water vapor and cloud feedbacks to ECS, we compare the feedback responses in CAM4 (with a low ECS) and CAM5 (with a high ECS) in an atmosphere-only configuration (Materials and Methods). Under identical Eocene boundary and sea surface conditions, both models exhibit similar increases in atmospheric water vapor content and tropopause height with warming (fig. S5), results that are consistent with previous studies showing a similar water vapor feedback between CAM4 and CAM5 (27, 28). In contrast, using an approximated PRP (APRP) method (Materials and Methods), we find that the magnitude of the shortwave cloud feedback parameter is much greater in CAM5 (0.40 W m−2 K−1) than in CAM4 (0.11 W m−2 K−1). Furthermore, in response to an identical GMST increase from 19.3° to 25.5°C, the shortwave cloud feedback parameter increases from 0.40 to 0.79 W m−2 K−1 in CAM5 but decreases slightly in CAM4 (fig. S6). These results are a strong confirmation that the shortwave cloud feedback is responsible for the high ECS and its increase with warming in CAM5 compared with CAM4.

Linkages of cloud feedback to cloud microphysical processes

Using the APRP method, we attribute the strong shortwave cloud feedback in our Eocene simulations, as compared to the overall weak cloud feedback in CAM4, to changes in cloud fraction and scattering (Fig. 3C and figs. S5C and S6B; Materials and Methods). With warming from our 1× to 6× simulations, low- and medium-level cloud cover decreases, especially at mid- and high latitudes (except over the Arctic; Fig. 4, A and B), a response that is different from the slight increase in cloud cover in CAM4 (Fig. 6) and consistent with the changes in the shortwave cloud forcing with warming. At the same time, cloud liquid droplet sizes increase from 1× to 6× simulations, with maximum values over mid-latitude regions rising from 12 to 15 μm (Fig. 4, E and F). In theory, when cloud liquid water content remains constant, clouds with larger droplets are less opaque and preferentially scatter in the forward direction, resulting in enhanced shortwave radiation to the surface (29). In our Eocene simulations, mid- and high-latitude in-cloud liquid water content is not fixed but increases (Fig. 4, C and D), partly offsetting the decrease in cloud opacity there. A decreased high-latitude cloud glaciation rate in a warmer climate could also decrease the cloud opacity (30). At lower latitudes, in-cloud liquid water content decreases with warming, which is another contributor to the decrease in cloud opacity. This thinning of clouds with warming is broadly consistent with a recent large eddy simulation of the subtropics (31). The combined reduction in cloud cover and cloud opacity increases surface shortwave radiation and surface warming.

Fig. 4 Zonal mean cloud cover, liquid water content, effective droplet size, and droplet number concentration in the Eocene 1× and 6× simulations.

(A) Zonal mean cloud cover (in %) in the Eocene 1× simulation. (B) As in (A) but for the Eocene 6× experiment. (C) Zonal mean in-cloud liquid water content (in g kg−1) in the Eocene 1× simulation. (D) As in (C) but for the Eocene 6× experiment. (E) Zonal mean in-cloud liquid droplet radius (in μm) in the Eocene 1× simulation. (F) As in (E) but for the Eocene 6× experiment. (G) Zonal mean in-cloud liquid droplet number concentration (in cm−3) in the Eocene 1× simulation. (H) As in (G) but for the Eocene 6× experiment. Zonal means of in-cloud variables are calculated for grid points with cloud liquid water greater than 1 parts per million.

The reduction in cloud coverage and the increase in cloud droplet size with warming in our simulations are ascribed primarily to the implementation of a new two-moment cloud microphysical scheme in CAM5 (22, 23). The new scheme predicts both cloud water mixing ratio and droplet number concentration, allowing a prognostic calculation of effective droplet radius. With the new scheme, cloud water content, particle size, and droplet number concentration in CAM5 are all in better agreement with recent satellite observations than previous models (22, 3234). In contrast, CAM4 includes a one-moment cloud microphysical scheme that does not allow cloud droplet size to change and, thus, does not capture the corresponding changes in cloud shortwave scattering (35).

The reduction in low-cloud cover in our simulations is hypothesized to arise from more efficient conversion of cloud water into precipitation under warming conditions via stronger accretion and autoconversion processes. As the two most important sinks of cloud condensates, autoconversion and accretion describe, respectively, the transformation of cloud particles to precipitation through coalescence and vapor diffusion and the growth of precipitation by droplet interception and accumulation of cloud water. In the Eocene 1× and 6× simulations, zonal mean in-cloud accretion and autoconversion rates increase everywhere in response to elevated CO2, with the maximum values increasing from ~1 to ~2 g kg−1 day−1 (Fig. 5). The increased accretion and autoconversion rates are associated with the increased mixing ratios of cloud water and rain water with warming and the slight reduction in droplet number concentration (Fig. 4, G and H) (22). We have confirmed the influence of the microphysical scheme on low-cloud cover through a series of CAM4 atmosphere-only simulations, in which we sequentially replaced the CAM4 physical parameterizations with the updated CAM5 ones (Materials and Methods). Replacing the CAM4 cloud microphysical parameterization with the new CAM5 scheme reproduces ~90% of the differences in low-cloud cover (Fig. 6). Other physical parameterizations play a minor role in decreasing low-cloud cover with warming.

Fig. 5 Increase in in-cloud accretion and autoconversion rates with warming in the Eocene simulations.

(A) Zonal mean in-cloud accretion rate (in g kg−1 day−1; conversion of cloud water to precipitation through rain droplets intercepting and collecting small cloud droplets) in the Eocene 1× simulation. (B) As in (A) but for the Eocene 6× experiment. (C) Zonal mean in-cloud autoconversion rate (in g kg−1 day−1; conversion of cloud water to precipitation by coalescence and vapor diffusion) in the Eocene 1× simulation. (D) As in (C) but for the Eocene 6× experiment. Zonal mean of in-cloud accretion and autoconversion rates are calculated for grid points with cloud cover greater than 1%.

Fig. 6 Role of CAM5 physical parameterizations in changing low-cloud cover.

Changes in low-cloud cover (in %) versus the GMST (in °C) in the Eocene atmosphere-only simulations using CAM4 with individual CAM5 physical parameterizations switched on sequentially (Materials and Methods). The physical parameterizations of radiation transfer (rad), cloud microphysics (micro), turbulence and shallow convection (uw), and cloud macrophysics (macro) are tested. For comparison, values of low-cloud cover have been realigned by subtracting the corresponding value in the 1× simulation. For reference, standard atmosphere-only simulations using CAM5 (red filled circles) and CAM4 (orange filled circles) are also shown.


The response of clouds to warming is responsible for much of the spread in model predictions of anthropogenically forced future climate change (10). Our study is one of many that have emphasized the role of clouds in amplifying CO2-induced warming (21, 30, 36, 37). These previous studies have attributed changes in cloud forcing to changes in SST patterns and the resulting decrease in tropospheric stability, thermodynamic effects from warmer SSTs, warming and drying by radiative forcing, increased convective mixing in the lower troposphere, and cloud interactions with aerosols (16, 17, 3638). Here, we identify cloud microphysical processes as an important factor in cloud forcing. The important role of cloud microphysics in determining climate sensitivity in the CESM1.2 justifies a greater effort to better understand and parameterize cloud microphysics through theory, observations, and simulations with multiple high-resolution models and super-parameterizations.

As the only period of high atmospheric CO2 in the Cenozoic, the Eocene is arguably the most critical target for benchmarking climate models that are tuned to present-day CO2 levels. Here, a state-of-the-art Earth system model simulates the extreme warmth and low meridional temperature gradient of the Eocene and PETM with age-appropriate estimates of past CO2 and without any alteration of the model physics. The good match of our Eocene simulations with proxy data is due to the fact that climate sensitivity in the CESM1.2 increases substantially with CO2-induced warming. The cloud processes responsible for the increased climate sensitivity in our Eocene simulations are also active under modern conditions. In a simulation with 4× CO2 based on the PI boundary conditions, the ECS increases by 33% to 4.2°C from 5.6°C in the PI simulation (fig. S7; Materials and Methods). These results suggest a higher climate sensitivity in a warmer future than typically estimated by the IPCC.


Model simulations

Fully coupled simulations. The CESM1.2 with the CAM5 (18) was used in this study. The released version of the CESM1.2 has been slightly modified to incorporate an upgrade in the later versions of the radiation code that corrects the missing diffusivity angle specifications for certain longwave bands, which is important for simulation of warm climates. All CESM model simulations were run with a horizontal resolution of 1.9° × 2.5° (latitude × longitude) for the atmosphere and land, and a nominal 1° for the ocean and sea ice. The atmosphere component CAM5 has 30 levels starting from ~60 m above surface to a model top at ∼2 hPa (~40 km). The Eocene simulations were initialized from an equilibrated PETM state using the CCSM3 (16) and integrated for more than 2000 years to ensure that the upper ocean reaches quasi-equilibrium (fig. S2). All Eocene simulations were run with the identical boundary conditions following the DeepMIP protocol (2) and differ only in atmospheric CO2 concentration. The Eocene boundary conditions included paleogeography, land-sea mask, vegetation distribution, and PI non-CO2 greenhouse gas concentrations, soil properties, natural aerosol emissions, solar constant, and orbital parameters (2, 39). The emissions of PI natural aerosol have been redistributed according to the Eocene paleogeography (40). For comparison, simulations based on the PI conditions with 1× and 4× PI CO2 were also conducted. All results presented here are based on climatological means of the past 100 years of each simulation.

SOM simulations. Following convention, SOM simulations (41) were used to calculate the ECS. In all SOM simulations, the CESM dynamic ocean model was replaced with a simple mixed-layer model, keeping all other model components unchanged. To reproduce the SST and sea ice conditions in the fully coupled simulation with dynamic ocean, surface heat flux convergence and mixed-layer depth were derived from ocean climatological state in a corresponding fully coupled run with identical boundary conditions and prescribed in the SOM simulation. SOM simulations typically reach equilibrium in fewer than 30 model years.

To calculate the ECS for different CO2 background states, two SOM simulations were conducted, one with the specified background CO2 (e.g., 1× or 3× for Eocene simulations and 1× or 4× for simulations with the PI conditions) and another with twice the background CO2 level. All runs were integrated for 60 years, with the past 30 years used for analysis. The ECS for a specific background state was estimated as the difference in 30-year averaged GMST between a SOM simulation and its corresponding doubled CO2 run. Eocene simulations experience a runaway greenhouse for CO2 at 12× PI levels, preventing us from directly calculating ECS for the 6× case using SOM simulations. Therefore, ECS for the 6× case was estimated from the fully coupled simulations. We first calculated the GMST difference to be 5.6°C between the 9× and 6× fully coupled simulations. We then used offline radiation code in CAM5 to calculate the adjusted radiative forcing within the 6× case by increasing CO2 to the 9× and 12× PI levels, which were 2.98 and 5.16 W m−2, respectively. ECS for the 6× case was then estimated to be 9.7°C by multiplying the GMST difference by the ratio of radiative forcing (5.16/2.98 = 1.73). We note that a similar runaway greenhouse has been reported in other models at much lower CO2 levels (e.g., 4× or 8×) (3). It is unclear whether the runaway greenhouse in the simulation is real. Further studies are needed to test its sensitivity to boundary conditions and model physics. Additional Eocene SOM simulations were carried out using CAM4 to compare its climate sensitivity and feedback strength to that of CAM5.

Atmosphere-only simulations. To investigate the differing cloud responses to warming in CAM4 and CAM5, Eocene atmosphere-only simulations were conducted for both models with prescribed climatological SST and sea ice derived from a corresponding fully coupled Eocene simulation with the same CO2 level. Any difference between these CAM5 and CAM4 atmosphere-only simulations is attributable to differences in the model physical parameterizations of radiation transfer, aerosol, boundary layer processes, shallow convection, and cloud microphysical and macrophysical processes. To further isolate the exact physical process responsible for the different cloud response between CAM5 and CAM4, additional atmosphere-only simulations were conducted using CAM4 with CAM5 physical parameterizations switched on sequentially. We first switched on the CAM5 radiation transfer scheme (rad), then the cloud microphysical parameterization (micro), then the boundary layer and shallow convection schemes (uw), and, lastly, the cloud macrophysical scheme along with a few other parameter changes (macro). The role of different aerosol schemes was diagnosed as the differences between standard CAM5 simulations and CAM4 simulations with CAM5 rad, micro, uw, and macro switched on. Note that when adding the uw scheme to CAM4, we also increased the vertical levels from 26 to 30, as is done in CAM5 (42). These atmosphere-only simulations enabled us to evaluate the impact of individual physical parameterizations on the cloud response to warming in CAM5.

Previous Eocene simulations. Eocene simulations using the CESM1.2 in this study were compared with published Eocene simulations, including those that participated in the Early Eocene Modeling Intercomparison Project (3, 11, 16, 19, 4346). Note that methane concentrations in previous simulations that differ from the PI value have been scaled to the equivalent CO2 concentration using published relationships (47).

Estimating GMST for the Early Eocene, pre-PETM, and PETM

To estimate the GMST for the Early Eocene, we made use of published compilation of terrestrial records (11) and compiled a new set of marine SST records spanning the Early Eocene (53.3 to 49.4 Ma), the pre-PETM period, and the PETM event in accordance with the DeepMIP protocol (2, 48). Specifically, the Eocene marine SST records consist of δ18O of planktic foraminifera from the Ocean Drilling Program (ODP) 690B (49), ODP 738 (50, 51), ODP 865 (15), Deep Sea Drilling Project (DSDP) 277 (52), Tanzania (53), Waipara (54), and Lodo Gulch (55); revised Mg/Ca temperatures (4) from ODP 865 (56), DSDP 277 (57), Hampden Beach (57), Tora NZ (57), Tawanui NZ (57), and Waipara (58, 59); TEX86 temperatures (60) from ODP 929 (61), ODP 959 (20), ODP 1172 (62), Arctic Coring Expedition (ACEX) (63), Hampden Beach (61), Hatchitigbee Bluff (64), South Dover Bridge (61), Tanzania (53), Waipara (54, 58, 61), Western Siberian Seaway (65), and Wilkes Land U1356A (66); and clumped isotope (Δ47) thermometry from Kutch India (4), Egem Belgium (4), and Hatchitigbee Bluff (64). The pre-PETM and PETM data consist of δ18O from Bass River (55), Wilson Lake (54, 67), DSDP 277 (52), ODP 865 (15), Tanzania (68), Nigeria SQ (6), Lodo Gulch (55), Tumey Gulch (55), Milville (69), ODP 689 (69), DSDP 401 (70), and DSDP 549 (71); revised Mg/Ca temperatures (4) from DSDP 401 (70), DSDP 277 (57), ODP 865 (56), DSDP 527 (72), ODP 1209 (73), Nigeria SQ (6), and Bass River (74); and TEX86 temperatures from Bass River (75), Fur Section and Store Bælt Section (North Sea) (76), Harrell Core (77), Nigeria 1B10A/B (6), Nigeria SQ (6), ODP 1172 (78), ODP 959 (20, 79), Wilson Lake (67), ACEX (63), Waipara (54), and Western Siberian Seaway (65). For accurate model-data comparison, paleolocations for all proxies were recalculated using the Herold14 reference frame (39), which is the paleogeography used for our simulations in accordance with DeepMIP recommendations (2). Proxies were also recalibrated in a consistent fashion. Δ47 and Mg/Ca data were converted to SST following Evans et al. (4). TEX86 data were calibrated using the BAYSPAR (linear bayesian spatially varying regression) approach (60). The TEXH86 calibration (80) was not used here due to its regression dilution, which will systematically underestimate warm Eocene temperatures (60). δ18O was recalibrated in a manner similar to the compilation in Lunt et al. (3). Using the Herold14 paleolocations, δ18O seawater estimates were drawn from the nearest grid cells from two Eocene isotope-enabled simulations conducted with Goddard Institute for Space Studies (GISS) model E-R (46, 81) and Hadley Centre Coupled Model version 3 (HadCM3) (82), respectively. These were converted to the Vienna Pee Dee Belemnite (VPDB) scale, and SSTs were calculated using the high-light and low-light calibrations of Bemis et al. (83). The upper and lower error bars for the δ18O therefore reflect both δ18O seawater uncertainty and calibration uncertainty. Note that for site DSDP 401, only the HadCM3 estimate of δ18O of seawater was used (the GISS estimate yielded an unrealistically low value of −3.13‰). The compiled proxy data average for each timeslice, along with paleolocations and corresponding source references, can be found in tables S1 and S2.

We estimated the Early Eocene GMST from proxy records using two methods. In the first method, we binned the proxy terrestrial and marine temperatures separately into 15° latitudinal bands and computed the arithmetic mean of land surface temperature and SST within each band. We next calculated the area-weighted global mean for land and ocean temperatures separately from available bands with records. Last, the GMST was obtained as the area-weighted average of land and ocean temperatures. In the second method, we first calculated proxy temperature anomalies from PI core-top temperatures. We then compiled a GMST anomaly following the same procedure as the first method. Absolute GMST for the Early Eocene was the sum of the complied GMST anomaly and PI GMST. The Early Eocene temperature records are spatially unevenly distributed with many fewer records in the tropics (fig. S3). As a result, the first method (using absolute temperature) underestimates the GMST, and the second method (using anomalous temperature) overestimates the GMST because of the high absolute and low anomalous Eocene temperature in the tropics. Our final GMST estimate is the average of results from the two methods, which removes part of the biases from spatially unevenly distributed proxies. We explored different latitudinal band sizes to estimate GMSTs and found similar results with a difference of <1°C. Years 1851–1900 of the Berkeley Earth surface temperature (84) were used to calculate the average climate state for the PI.

A simulation-aided Monte Carlo approach was developed to estimate the uncertainty of proxy GMST, including the propagation of calibration uncertainty from individual proxies and the sampling error from scarce records. To propagate uncertainty from individual proxies, we randomly drew the same number and type of surface temperature within a 5° box around each proxy (to incorporate uncertainty in paleolocations) in an Eocene simulation and calculated the GMST in the same way as we did for proxy records. We repeated this procedure for 10,000 iterations and calculated the 95% uncertainty interval to represent the uncertainty in proxy GMST propagated from individual records. To calculate sampling uncertainty from scarce records, we randomly drew the same number and type of surface temperatures from all model grids and compiled the GMST. We calculated the sampling uncertainty from 10,000 iterations. The total uncertainty was the sum of uncertainties from propagation and sampling. In this approach, we assumed that the uncertainty in proxy GMST due to sparse sampling and uncertain location is comparable to that from the CESM Eocene simulation. This assumption is reasonable, as our Eocene simulations closely capture the warming and reduced SST gradient in proxy reconstructions (Figs. 1 and 2, and figs. S3 and S4). Our Monte Carlo estimated uncertainty does not vary much (<0.5°C) whether the 3×, 6×, or 9× CO2 Eocene simulation was used in the procedure.

Calculating the latitudinal SST gradient

The latitudinal SST gradient in model simulations and PI observations was calculated as the SST difference between the averages in the low latitudes (|lat| < 30°) and the high latitudes (|lat| > 60°) (3). The SST gradient in Eocene proxy reconstructions was calculated following Evans et al. (4) We subtracted the mean of the benthic foraminifera Mg/Ca-derived deep ocean temperature record (85) from the mean of all available tropical SST data (i.e., records within 30°S to 30°N in our Early Eocene SST compilation; table S1) for the interval 49.4 to 53.3 Ma. The uncertainty was derived using a similar Monte Carlo approach as described above for compiling GMST. Latitudinal SST gradients calculated in this way may reflect maximum steepness endmembers (4).

The PRP and feedback analysis

The PRP approach (25) was adopted to diagnose the strength of cloud, water vapor, lapse rate, Plank, and surface albedo feedbacks in model simulations. PRP calculates top of the atmosphere radiative perturbations from an offline version of the model radiation code. First, we output high-frequency instantaneous radiation fields for both the control and the perturbed simulations. Then, offline radiation code was driven by the high-frequency radiation fields substituted one at a time. We adopted a two-way calculation; we first calculated the radiative perturbation by substituting fields from a lower CO2 run into a run with higher CO2 and then repeated the process substituting fields from a higher CO2 run into lower CO2 run. The final radiation perturbation was obtained by averaging perturbations from the two substitutions. This technique has the advantage of largely reducing the error from two correlated fields (25). Note that the stratospheric temperature was not substituted when calculating the lapse rate feedback. The tropopause was diagnosed from instantaneous model output following the method of Reichler et al. (86). The Planck feedback was obtained by simply perturbing all model levels by the surface temperature change. We used the Parallel Offline Radiative Transfer tool released together with the CESM to carry out the offline radiation calculation (87). The sampling frequency in this study is 73 model steps, which is found to be a good balance of sampling frequency, data size, and accuracy (87). In each case, the offline radiation code was integrated for five model years, with the past 4 years used for analysis.

A less-expansive APRP method (88) was used to further decompose the shortwave cloud feedback into contributions from changes in cloud fraction, scattering, and absorption and to compare clouds in CAM4 and CAM5, as the full PRP method demands large storage of high-frequency global fields and considerable computational resources. The APRP method has been found to be efficient and satisfactory for diagnosing shortwave cloud feedback parameters, with a difference from the full PRP method of less than 7% (88).


Supplementary material for this article is available at

Fig. S1. Topography and bathymetry at model resolution.

Fig. S2. Spin-up of the Eocene simulations.

Fig. S3. Model-data comparison of Early Eocene surface temperature.

Fig. S4. Model-data comparison of PETM warming in SST.

Fig. S5. Comparison of CAM5 and CAM4 Eocene atmosphere-only simulations.

Fig. S6. Comparison of temperature and shortwave cloud feedback parameters between CAM5 and CAM4 Eocene SOM simulations.

Fig. S7. Increase in climate sensitivity and cloud feedback parameter with warming under modern conditions.

Table S1. Compilation of SST proxies for the Early Eocene.

Table S2. Compilation of SST proxies for the pre-PETM and PETM.

Table S3. Feedback analysis from the PRP calculations.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank J. Kiehl, C. Shields, and M. Rothstein for providing the CESM code and boundary/initial condition files for our simulations. We thank C. Skinner, F. Ding, X. Huang, and A. Gettelman for the helpful discussions. We thank the editor D. Lea, D. Lunt, and three anonymous reviewers for their constructive comments leading to the improvement of the manuscript. We acknowledge the computational resources provided by the CESM Paleoclimate Working Group and the high-performance computing support from Cheyenne ( provided by NCAR’s CISL, sponsored by the NSF. Funding: This work was supported by Heising-Simons Foundation grant #2016-015 to C.J.P. and J.E.T. Author contributions: J.Z. and C.J.P. designed the study. J.Z. performed the numerical experiments and analyses with help from C.J.P. J.E.T. produced the SST synthesis. J.Z. and C.J.P. prepared the manuscript with input from J.E.T. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper are available in the Zenodo repository at doi:10.5281/zenodo.2642536. Further data related to this paper may be requested from J.Z. The CESM model code is available through the National Center for Atmospheric Research software development repository (

Stay Connected to Science Advances

Navigate This Article