More recently, the community has started to move to “deep‐time” (here defined as time periods older than ∼5 million years ago) warm climates, during which CO 2 concentrations were very high (∼1,000 ppmv or more), in an attempt to find time periods where the climate signal and forcings were largest, with the expectation that signal‐to‐noise ratio will be reduced and climate sensitivity estimates will be more robust (e.g., Anagnostou et al., 2016 ; Cramwinckel et al., 2018 ; Shaffer et al., 2016 ; Zeebe et al., 2009 ). However, many such time periods were essentially free from continental ice sheets, and tectonic changes over millions of years mean that continental positions, mountain ranges, and ocean bathymetry were all very different compared to modern (e.g., Markwick, 2007 ). As such, it is possible that estimates of paleo climate sensitivity from these time periods may have limited relevance for assessing modern climate sensitivity, if climate sensitivity is paleogeography dependent. Furthermore, the upper end of CO 2 concentrations during many of these time periods are so high that nonlinearities in forcings and feedbacks (Bloch‐Johnson et al., 2015 ; Caballero & Huber, 2013 ) may further complicate their relationship to the short‐term future.

Given the societal importance of climate sensitivity, it is unsurprising that many studies have turned to observations of past climates in an attempt to further refine, and test, estimates of this key parameter (e.g., see reviews by Rohling et al., 2018 ; Schmidt et al., 2014 ). Any such paleofocused attempt to characterize climate sensitivity requires knowledge of (a) the magnitude of the climatic response to a particular forcing and (b) the magnitude of the forcing itself. As such, many attempts to characterize climate sensitivity from the geological record have focused on the relatively recent past, such as the Last Glacial Maximum (e.g., Crucifix, 2006 ; Hansen et al., 2008 ; Hargreaves et al., 2012 ), or the Pleistocene (e.g., Friedrich et al., 2016 ; Kohler et al., 2018 ) where data sets for temperature change have relatively good spatial coverage, and the forcings are well understood and quantifiable. However, such studies still have relatively large uncertainties due to difficulties associated with the interpretation of paleoproxy records and/or have often conflated feedbacks and forcings and so have either underestimated or overestimated sensitivity (PALAEOSENS Project Members, 2012 ). Furthermore, the forcing during the Last Glacial Maximum and much of the Pleistocene was negative relative to modern, and so the relevance for future climate change, where forcing is positive, may be limited. Other studies have focused on warm periods within the last 5 million years, such as the Pliocene (e.g., Hargreaves & Annan, 2016 ; Haywood et al., 2013 ; Lunt et al., 2010 ; Martínez‐Botí et al., 2015 ), when CO 2 concentrations were similar to that of today; however, even these studies do not explore CO 2 concentrations of ∼560 ppmv and above, that are more relevant to climate sensitivity and typical end‐of‐century scenarios (Moss et al., 2010 ).

Climate sensitivity, in its very broadest terms, is a measure of the response of the climate system to an external forcing. More specifically, it is defined as the global annual mean near‐surface air temperature equilibrium response to a defined forcing, usually either 1 W/m 2 (e.g., PALAEOSENS Project Members, 2012 ) or a doubling of CO 2 (∼3.7 W/m 2 , e.g., Charney et al., 1979 ). Despite (or perhaps due to) its simplicity, it is a widely used metric and is a key variable in many impact (e.g., Donner et al., 2005 ), policy (e.g., Rogelj et al., 2014 ), and economic (e.g., Hope, 2006 ) assessments. The Intergovernmental Panel on Climate Change ( 2013 ) assesses that climate sensitivity (defined in terms of a doubling of CO 2 ) is likely to be in the range 1.5 to 4.5 °C, very unlikely to be more than 6 °C, and extremely unlikely to be less than 1 °C.

At the end of the model integration, all the simulations are very well spun up, with only small residual trends in ocean temperature, even in the deep ocean (the mean absolute global mean temperature trend at 2,700‐m depth of all simulations in the last 1,000 years of simulation is 0.06 °C per millennium; see supporting information Figure S1 ). In addition, Gregory plots (Gregory et al., 2004 ) of the energy balance at the top of the atmosphere (TOA) against global mean surface temperature ( supporting information Figure S2 ) show that the energy balance is very close to zero for all simulations and that any residual energy balance is decreasing over time (the mean absolute TOA inbalance across all simulations is 0.036 W/m 2 , and the maximum TOA inbalance is +0.09 W/m 2 in the Cenomanian at ×4). However, there is a drift in the global mean salinity due primarily to the rigid‐lid approximation in the ocean model, which includes minimum and maximum salinity thresholds (0 and 45 psu, respectively). This, coupled with the presence of enclosed or restricted basins (such as the Arctic Ocean) with a net input of freshwater from the surrounding river catchments, means that the restricted basins become very fresh and reach the minimum salinity threshold at the surface while the rest of the ocean tends to increase in salinity. We do not expect this drift to affect our results because (i) the absolute salinity does not have a strong effect on the density gradients due to the linearity of the equation of state (Bryan & Cox, 1972 ) in the regime of the majority of the global ocean, (ii) the enclosed basins that reach the minimum salinity threshold have a relatively small area and do not affect the global ocean circulation, and (iii) the maximum salinity threshold is only reached in a few isolated grid cells in the tropics in the later Eocene simulations (supporting information Figure S3 ). Indeed, the mixed layer depths in our simulations (supporting information Figure S6 ) are very similar to those presented in Lunt et al. ( 2016 , their supplement Figure S10 ), indicating that the ocean circulation is relatively unchanged after 9,000 years, despite the gradual drift in salinity and gradually equilibrating temperatures.

The climate model used for the simulations is very similar to the HadCM3BL‐M2.1aD model described and evaluated in Valdes et al. ( 2017 ), except that we include a modification to the ozone profile to ensure that the model does not develop a runaway warming at ×4 CO 2 , as discussed in Lunt et al. ( 2016 ). The model includes representations of “Charney” feedbacks associated with changes in clouds, snow/seaice, lapse rate and water vapor, and the “Earth system” feedback associated with changes in vegetation (Lunt et al., 2010 ), but does not include many other Earth system feedbacks associated with changes to, for example, dust, methane, nitrogen dioxide, or ozone. The sign of many of these missing feedbacks is uncertain (Intergovernmental Panel on Climate Change, 2013 ), so it is not clear if our simulations will over or underestimate “true” long‐term Earth system sensitivity.

The climate model and the ×4 simulations are identical to those described in detail in Lunt et al. ( 2016 ), except that here they have been run for longer (10,422 years here compared with 1,422 years), in order to approach more closely to equilibrium, and they use the fully dynamic mode (as opposed to equilibrium mode) of vegetation model.. The ×2 simulations are branched off from the ×4 simulations after 422 years and run for 10,000 years. The only difference between the simulation of each stage relative to another stage is the prescribed paleogeography (see Figure S1 in the supplement of Lunt et al., 2016 ) and the solar constant (see Figure 2 of Lunt et al., 2016 ). All simulations have a modern astronomical configuration and are initialized from a homogeneous global ocean salinity and an idealized zonal mean ocean temperature structure, which is a cosine function of latitude (see section 2.5 of Lunt et al., 2016 ).

In order to explore the evolution of climate and climate sensitivity over geological timescales, we carry out an ensemble of 19 pairs of fully coupled atmosphere‐ocean‐vegetation climate model simulations with varying paleogeography and solar luminosity appropriate for each geological stage from the earliest Cretaceous (Berriasian stage, ∼143 Ma) to the latest Eocene (Priabonian stage, ∼36 Ma), at ×2 and ×4 preindustrial atmospheric CO 2 concentrations (i.e., 560 and 1,120 ppmv).

3 Results

3.1 Climate and Climate Sensitivity Figure 1a shows the global annual mean near‐surface air temperature for the high (×4) and low (×2) CO 2 simulations (and see Table S1 in the supporting information). Both sets of simulations show a similar trend (the r correlation coefficient between the two time series is 0.86; supporting information Figure S4). In particular, there is a cooling in the very earliest Cretaceous between the Berriasian and Valanginian at 140 Ma, a warming from ∼130 to 90 Ma peaking at the Turonian, and a cooling from 90 to 80 Ma to a minimum in the Campanian. In the Paleocene and Eocene the trends are less consistent across the two sets of simulations, but include a maximum in the Ypresian. Figure 1 Open in figure viewer PowerPoint (a) Global annual mean near‐surface air temperature in the ×2 (dark blue circles) and ×4 (red circles) CO 2 simulations, and the Gelasian ×1 simulation (green square). For each set of simulations the curved black line is a spline fit to the data, of resolution about 5 million years. The nonlinear model fit to the ×2, ×4 simulations is shown with light blue circles. (b) Global mean climate sensitivity (i.e., ×4 minus ×2) over time (black circles). Nonlinear model sensitivity is shown with light blue circles. Red vertical lines show those time periods that have a switch in ocean circulation between the ×2 and ×4 simulations. (c) Ocean area over time (as a fraction of the global surface area, calculated from the paleogeographies at the model resolution). (d) Solar constant over time (W/m2). The global mean temperature in our ×4 simulations has a greater total range across the geological stages than that in the shorter simulations of Lunt et al. (2016; ∼3 °C compared with ∼1 °C; compare our Figure 1 with their Figures 5 and 7), and the regional expression of this range also differs, but the main conclusion of Lunt et al. (2016), that changes in paleogeography have a substantial effect on regional climate that must be considered when interpreting long proxy records, is robust. Figure 1b shows the climate sensitivity of each stage (the difference in global mean temperature between the ×2 and ×4 simulations); regional plots are shown in supporting information Figure S5. The mean climate sensitivity over this period is 4.2 °C, the minimum is 3.7 °C (Berriasian stage, ∼145 Ma), and the maximum is 5.3 °C (Maastrichtian stage, ∼70 Ma). The change in global climate sensitivity over time includes a gradual increase in sensitivity from the earliest Cretaceous to the latest Cretaceous, a decrease to a minimum in the Paleocene, and then an increase followed by a decrease in the Eocene. Our modeled range of climate sensitivities in the Eocene (3.8 to 4.7 °C) is consistent with estimates of 3.5 to 8.9 °C (Cramwinckel et al., 2018) and 2.9 to 4.0 °C (Anagnostou et al., 2016), inferred from temperature and CO 2 proxies for this time period. A quantitative analysis of the driving mechanisms is carried out in section 3.2, but in qualitative terms, the global mean temperature in the simulations appears to be driven by a combination of CO 2 (explaining the offset in temperature between ×2 and ×4), increasing solar constant (explaining a general increase in temperature from the earliest Cretaceous to latest Eocene in both ×2 and ×4), plus a more variable signal that must be associated with changing paleogeography. This paleogeographic signal appears to be strongly related to the area of ocean in the prescribed paleogeographies (Figure 1c); the ocean area correlates with both the ×4 (r correlation coefficient = 0.63) and ×2 simulations (r correlation coefficient = 0.42). The mechanism for this paleogeographic forcing is discussed in section 3.2. However, there are some differences between ×2 and ×4 that do not fit this qualitative relationship and that can be seen in the evolution of climate sensitivity. In particular, there is a large swing in sensitivity around the Cretaceous‐Paleocene boundary, with relatively high sensitivity at the Maastrichtian and low sensitivity in the Selandian, due to a relatively warm Maastrichtian at ×4 and a relatively warm Selandian at ×2 (Figures 1a and 1b). These anomalies are associated with changes in ocean circulation that take place as a function of CO 2 at these two stages (supporting information Figure S6). The ocean circulation around this time has two possible states: (i) a relatively warm global state characterized by deep water formation in the south Pacific Ocean (Maastrichtian at ×4 and Selandian at ×2; supporting information Figures S6b and S6g) and (ii) a relatively cold global state characterized by deep water formation in the South Atlantic and Indian Oceans (Maastrichtian at ×2 and Selandian at ×4; supporting information Figures S6c and S6f). Although changes in ocean circulation may not be expected to affect the global mean temperature to a great extent, in this case increased poleward heat transport in the Pacific basin results in a greater increase in global mean temperature than that caused by increased poleward heat transport in the Atlantic/Indian basins. This enhanced Pacific poleward heat transport may be related to the greater areal extent of the Pacific basin and is enhanced by sea ice albedo feedbacks around Antarctica. As such, the Cretaceous‐Paleocene boundary is a particularly sensitive time period to relatively small changes in paleogeography and/or CO 2 , resulting in anomalously large or small climate sensitivities. For all other time periods, the deep water formation occurs in the same regions at ×4 as at ×2, although the region varies as a function of paleogeography (e.g., in the South Atlantic in the Priabonian, supporting information Figures S6a and S6e; in the North Pacific in the Berriasian, supporting information Figures S6d and S6h).

3.2 Quantitative Role of Nonlinear Feedbacks in Determining Climate Sensitivity 1984 F all (W/m2), where F all is the sum of three separate forcings: F co2 due to a change in CO 2 , F solar due to a change in solar constant, and F geog due to a change in paleogeography. Furthermore, the temperature difference relative to the reference state, ΔT, is given by (1) λ, the feedback parameter (W · m−2·K−1), is negative. If we choose the Turonian at ×2 as our reference simulation, then F co2 is 0 for all the ×2 simulations and equal to ΔR for all the ×4 simulations, where ΔR is the radiative forcing for a doubling of CO 2 . F solar can be approximated by (1−α p )ΔS 0 /4, where ΔS 0 is the change in solar constant relative to the Turonian and α p is the planetary albedo. For the paleogeographic forcing, F geog , the apparent importance of ocean area discussed in section 1980 2019 2006 2016 F geog , can be approximated by Δα s S 0 /4, where Δα s is the change in surface albedo relative to the Turonian. As well as the qualitative considerations above, a quantitative analysis of the simulations is also possible. Here, we focus on the global mean response through application of a radiative forcing and feedback framework. Traditional linear radiative forcing theory (Hansen et al.,) states that if we choose one simulation as a reference state, all other simulations can be considered as a perturbation to this reference state, caused by a forcing(W/m), whereis the sum of three separate forcings:due to a change in COdue to a change in solar constant, anddue to a change in paleogeography. Furthermore, the temperature difference relative to the reference state, Δ, is given bywhere, the feedback parameter (W · m·K), is negative. If we choose the Turonian at ×2 as our reference simulation, thenis 0 for all the ×2 simulations and equal to Δfor all the ×4 simulations, where Δis the radiative forcing for a doubling of COcan be approximated by (1−)Δ/4, where Δis the change in solar constant relative to the Turonian andis the planetary albedo. For the paleogeographic forcing,, the apparent importance of ocean area discussed in section 3.1 suggests a role for surface albedo radiative forcing. Using a simple energy balance model, Barron et al. () showed that the albedo forcing, due to changes in continental versus ocean area caused by changes in sea level, had a substantial effect on global mean temperature on geological timescales. The importance of ocean area through surface albedo forcing is also supported by more recent work (Pohl et al.,). However, others have suggested that the primary role of paleogeography on global mean temperature is through seasonality effects associated with the degree of fragmentation of continents (Donnadieu et al.,; Ladant & Donnadieu,). Nonetheless, assuming changes in surface albedo to be the primary forcing due to paleogeography,, can be approximated by Δ/4, where Δis the change in surface albedo relative to the Turonian. R, which is constant, giving a constant climate sensitivity of −ΔR/λ. This is clearly at odds with the variable climate sensitivity in the full climate model (Figure a (W·m−2·K−2), into equation 2015 (2) However, this linear framework predicts climate sensitivity to be constant over geological time, because for each pair of ×2,×4 simulations, the only difference in forcing is equal to Δ, which is constant, giving a constant climate sensitivity of −Δ. This is clearly at odds with the variable climate sensitivity in the full climate model (Figure 1 c), as also highlighted by the non‐unity slope of the best fit line between the ×2 and ×4 simulations ( supporting information Figure S4 ). As such, we introduce a nonlinear parameter,(W·m·K), into equation 1 , following Bloch‐Johnson et al. (), such that a represents a temperature dependence (i.e., a nonlinearity) of the strength of the feedbacks. Solving the resulting quadratic equation, and taking the physically realistic negative root, gives (3) e solar and e geog , which allow the efficacy (i.e., the temperature response due to a unit forcing, Hansen et al., 2005 2 (for which the efficacy is equal to 1). In order to calculate the CO 2 forcing (F CO2 ) and the solar forcing (F solar ) we take ΔR=3.7 W/m2 (Houghton et al., 2001 S 0 from the solar constants prescribed in each simulation (Figure F geog ), we simply assume that the surface albedo change, Δα s =(α l −α o )ΔA o where α l =0.14 and α l =0.06 are typical global mean albedos of land and ocean, respectively, in the simulations, and ΔA o is the relative change in area of ocean given in Figure λ, a, e solar , and e geog . A minimal subjective tuning shows that λ=−1 W·m−2·K−1, a=0.04 W·m−2·K−2, e solar =0.7, and e geog =1.25 gives a good qualitative match to our results for absolute temperature at ×2 and ×4, and climate sensitivity (Figures λ=−1 W·m−2·K−1 corresponds to a linear climate sensitivity of 3.7 °C per CO 2 doubling, which is within the range estimated by the IPCC. We also note that a=0.04 W·m−2·K−2 is within the range estimated by Bloch‐Johnson et al. ( 2015 a ≤ 0.06 W·m−2·K−2 from various climate model simulations carried out at multiple CO 2 concentrations. Finally, we note that a value of e geog close to unity also supports the assumption that surface albedo forcing associated with ocean area is a mechanism for variations in F geog that is consistent with simple radiative forcing considerations. In principle, these four unknowns (λ, a, e solar , and e geog ) could be objectively tuned to give a best fit to the data in Figure 2015 2015 λ and a, this occurs at a total forcing of +6.25 W/m2 relative to the Turonian at ×2, or about 1,800 ppmv of CO 2 . The parameterrepresents a temperature dependence (i.e., a nonlinearity) of the strength of the feedbacks. Solving the resulting quadratic equation, and taking the physically realistic negative root, giveswhere we have included factorsand, which allow the efficacy (i.e., the temperature response due to a unit forcing, Hansen et al.,) of the solar and paleogeography terms to differ from that of CO(for which the efficacy is equal to 1). In order to calculate the COforcing () and the solar forcing () we take Δ=3.7 W/m(Houghton et al.,) and Δfrom the solar constants prescribed in each simulation (Figure 1 d). For the paleogeographic forcing (), we simply assume that the surface albedo change, Δ=()Δwhere=0.14 and=0.06 are typical global mean albedos of land and ocean, respectively, in the simulations, and Δis the relative change in area of ocean given in Figure 1 c. In order to account for the greater impact on radiative forcing of albedo changes in the tropics, we calculate the paleogeographic forcing as a function of latitude, accounting for varying incoming solar radiation, and weighting by area. This leaves equation 3 with four unknowns;, and. A minimal subjective tuning shows that=−1 W·m·K=0.04 W·m·K=0.7, and=1.25 gives a good qualitative match to our results for absolute temperature at ×2 and ×4, and climate sensitivity (Figures 1 a and 1 b, light blue circles). Two exceptions are the Maastrichtian and Selandian stages, which we have already shown are associated with changes in the mode of ocean circulation, which are not captured in this simple forcings/feedbacks framework. We note that=−1 W·m·Kcorresponds to a linear climate sensitivity of 3.7 °C per COdoubling, which is within the range estimated by the IPCC. We also note that=0.04 W·m·Kis within the range estimated by Bloch‐Johnson et al. () of −0.04 ≤≤ 0.06 W·m·Kfrom various climate model simulations carried out at multiple COconcentrations. Finally, we note that a value ofclose to unity also supports the assumption that surface albedo forcing associated with ocean area is a mechanism for variations inthat is consistent with simple radiative forcing considerations. In principle, these four unknowns (, and) could be objectively tuned to give a best fit to the data in Figure 1 ; however, here we are most interested in the general principle, which is that the evolution in climate sensitivity can be explained by nonlinear feedbacks in the framework of Bloch‐Johnson et al. (). It is worth noting that, as discussed by Bloch‐Johnson et al. (), a negative term inside the square root in equation 3 can be interpreted as the system developing a runaway greenhouse. For our values ofand, this occurs at a total forcing of +6.25 W/mrelative to the Turonian at ×2, or about 1,800 ppmv of CO In order to explore the origin of these nonlinear feedbacks on a global scale, we carry out a zonal mean energy balance analysis of the simulations, following Heinemann et al. (2009). This shows that both albedo changes and emissivity changes contribute to the nonlinearity in response to the solar and paleogeographic forcing. Given the relative warmth of the deep‐time model simulations, it is likely that the nonlinearity in albedo response is driven by nonlinearities in short wave cloud feedbacks, rather than snow and sea ice feedbacks. Indeed, regionally it is the tropics that dominate the nonlinearity in climate sensitivity (supporting information Figures S7a–S7c). The changes in ocean area are dominated by changes in the Northern Hemisphere midlatitudes; the maximum in ocean area is during the Turonian (mid‐Cretaceous) when the Western Interior Seaway in North America is extensive (see Figure S1 in the supplement of Lunt et al., 2016). The ocean area in the paleogeographies is consistent with sea level reconstructions of the Phanerozoic, such as Müller et al. (2008), derived primarily from the age‐area and depth‐area distribution of the modern‐day ocean floor. These reconstructions also show a maximum in the mid‐Cretaceous; this agreement is unsurprising given that sea level reconstructions inform the production of the paleogeographies. However, the details of the coastlines in the paleogeographies are uncertain, in particular for greenhouse climates (Sømme et al., 2009). Furthermore, the nonlinearity found in our simulations supports previous studies that found that past climate sensitivity was a function of baseline CO 2 (e.g., Pohl et al., 2014, for the Ordovician) and suggests that further nonlinear behavior would occur if we carried out simulations at ×1 or ×8 CO 2 .

3.3 Implications for Assessment of Future Climate Sensitivity It has been proposed that observational evidence of the climate sensitivity of the past may be useful for constraining future climate sensitivity. To test this hypothesis, we explore the relationship between modern and past climate sensitivity. To this end, in addition to the Cretaceous, Paleocene, and Eocene simulations discussed above, we also carry out a pair of Gelasian (early Pleistocene, ∼2 Ma) simulations, spun up in an identical manner to the other simulations, but at CO 2 concentrations of ×1 and ×2 (Figure 1a). The Gelasian is the closest available paleogeography to modern that has been produced using the same procedures as the deep‐time paleogeographies. In addition to the solar constant, paleogeography and CO 2 , this pair differs from the others in that it has extensive Antarctic and Greenland ice sheets. The Gelasian climate sensitivity (3.7 °C) is at the lower end of the ensemble. Given the dependence of climate sensitivity on temperature (i.e., the nonlinearity of climate sensitivity) in the model simulations, this implies that for the modern, the relatively small ocean area (Figure 1c) combined with the presence of the Antarctic and Greenland ice sheets, more than offsets the relatively high solar constant (Figure 1d), resulting in a relatively low climate sensitivity (and relatively cold ×2 simulation in the context of the increasing trend through the Cretaceous‐Paleocene‐Eocene). The fact that CO 2 increases from ×1 to ×2 rather than ×2 to ×4 likely plays an additional role in the low modern sensitivity, through nonlinearities in both feedbacks and CO 2 forcing (Caballero & Huber, 2013). The regional expression of climate sensitivity in the Gelasian includes an area of cooling in the North Pacific (supporting information Figure S5a), but this does not appear to be associated with changes in deep water formation. Cooling in this region is not seen in CO 2 ‐doubling model simulations associated with CMIP5 (Intergovernmental Panel on Climate Change, 2013). This could either be because the Gelasian has some paleogeographic differences to modern (notably, land in the region of what is today the Barents Sea, absence of a Hudson Bay, and a slightly open Panama Seaway), or because these simulations are run much longer than CMIP5 simulations. If we take the Gelasian results as representative of the modern climate sensitivity, then the implication is that the use of observational geological data from the Cretaceous, Paleocene, or Eocene (e.g., Anagnostou et al., 2016; Cramwinckel et al., 2018) to constrain future climate sensitivity should be carried out with care and could in general lead to higher estimates than are appropriate for modern. Indeed, the dependence of feedbacks on baseline temperature highlighted in this work provide a plausible explanation for why Anagnostou et al. (2016) found a higher climate sensitivity (∼4 °C per CO 2 doubling) in the relatively warm early Eocene climatic optimum than in the relatively cool late Eocene (∼3 °C per CO 2 doubling). Our findings are also consistent with Cramwinckel et al. (2018, their Extended Data Figure 9), who, relative to the late Eocene, find a higher sensitivity for the warm EECO than for the cooler middle Eocene.