Research ArticleCLIMATOLOGY

Slow climate mode reconciles historical and model-based estimates of climate sensitivity

See allHide authors and affiliations

Science Advances  05 Jul 2017:
Vol. 3, no. 7, e1602821
DOI: 10.1126/sciadv.1602821


The latest Intergovernmental Panel on Climate Change Assessment Report widened the equilibrium climate sensitivity (ECS) range from 2° to 4.5°C to an updated range of 1.5° to 4.5°C in order to account for the lack of consensus between estimates based on models and historical observations. The historical ECS estimates range from 1.5° to 3°C and are derived assuming a linear radiative response to warming. A Bayesian methodology applied to 24 models, however, documents curvature in the radiative response to warming from an evolving contribution of interannual to centennial modes of radiative response. Centennial modes display stronger amplifying feedbacks and ultimately contribute 28 to 68% (90% credible interval) of equilibrium warming, yet they comprise only 1 to 7% of current warming. Accounting for these unresolved centennial contributions brings historical records into agreement with model-derived ECS estimates.

  • Climate Sensitivity
  • Radiative Forcing
  • Climate Feedbacks
  • Bayesian Inference


Estimates of equilibrium climate sensitivity (ECS) from global climate models (GCMs) (1) and paleoclimate records (2) are generally consistent with a range of 2° to 4.5°C, but a number of studies based on historical instrumental records (37) yield a lower range of 1.5° to 3°C (8). These systematically lower observational estimates have been interpreted as demonstrating that GCMs are overly sensitive to CO2 forcing and that the ultimate amount of warming that Earth would experience at a given concentration of greenhouse gases is less than previously thought (5, 7).

A major challenge in inferring ECS from instrumental records is that the current climate system is not in energetic equilibrium. On average, Earth’s surface currently takes up between 0.1 and 0.9 W/m2 more heat than it loses (9), in which this rate of heating is denoted as H. To extrapolate to the temperature at which radiative equilibrium would be reestablished, a proportionality between changes in outgoing radiation, R, and global mean temperature, T, is generally assumed (38, 10), λ = (FH)/T. F represents anomalies in downward radiative forcing, and the difference gives the upward radiative response, R = FH. If λ is assumed constant over time, then ECS can be inferred by zeroing out H and assigning F equal to the forcing associated with a doubling of atmospheric CO2, giving an estimated ECS of Fλ−1.

The validity of this extrapolation is questionable, however, because of the well-documented time dependence of λ, or the net radiative feedback, in climate simulations (1118). This time dependence is a prime suspect for the systematic differences between GCMs and historical estimates of ECS (19). So far, efforts at resolving the discrepancy have been stymied by the lack of methodology to quantitatively account for changes in λ (8). We derive a generalized fit of the change in λ in 24 GCM simulations from phase 5 of the Climate Model Intercomparison Project (CMIP5) (20) and use the results to estimate the bias inherent in historical energy budget constraints on climate sensitivity.


Time dependence in the net radiative feedback of GCMs is related to changing patterns of warming (12, 17, 21) and interannual tropospheric adjustment (22). Various simple linear models have been postulated to describe this variability in λ. One set includes a fast-responding upper ocean that is coupled to a more slowly responding deep ocean (14, 23), sometimes including an efficacy term that modifies how much surface warming is associated with oceanic heat fluxes (13, 24). Another simple model represents the presence of distinct regions that have their own response time scales and feedbacks (16).

What model formulations are most useful for representing the heterogeneous evolution of the forced response of GCMs remains a subject of ongoing research, but it is possible to generically represent generically how global mean temperature evolves in existing simple models (14, 16, 2224), or any other linear system, using eigenmode decomposition (25)Embedded Image(1)where * represents the convolution operator. Each eigenmode of the temperature response, Tn, characterizes the rate and magnitude at which global temperature is excited by radiative forcing. The rate of response is governed by the time scale, τn, and the ultimate magnitude by a fractional contribution, αn, of the ECS, T/F. The amount of warming that is realized at time t is also a function of the structure of the forcing, F, where eigenmodes with small values of τ are more strongly excited by high-frequency variations than those with larger values.

We expand the classical eigenmode decomposition of the temperature response (2528) by jointly considering the decomposition of the radiative response. In any linear model of the climate system, the eigenmodes of the radiative response are directly proportional to those of the temperature response, Rn = λnTn, with each eigenmode potentially having a distinct net feedback parameter λn (29). In particular, eigenmodes are inferred for simulated global average temperature, T(t), and top-of-atmosphere energy flux, H(t), using a full Bayesian methodology that accounts for autocorrelated noise in T and H, uncertainties in extrapolating T and H to full equilibrium, and uncertainties in forcing, F. Performing the inference on both T and H allows for a joint posterior distribution of the eigenmode decomposition and the equilibrium sensitivity without a priori specifying an ECS, as is necessary with other methods (23, 25). Furthermore, the approach accounts for the fact that coefficients governing the distribution of allowable autocorrelations and uncertainties (the so-called hyperparameters) themselves need to be inferred from the data.

Despite the advantages of the eigenmode approach, there is also a limitation: Although the climate system is anticipated to evolve along a continuum of climate modes (30), the duration of the model simulations constrains us to fitting a small number of modes. Three eigenmodes are used on the basis of fewer modes leading to systematic discrepancies between our fit and the simulations during initial model adjustment (fig. S1), and the use of more eigenmodes failing to improve the fit, as judged by a Bayesian information criterion. This number of eigenmodes has been used elsewhere to describe the global temperature response of GCMs to both an abrupt increase in CO2 and a slow ramp-up of the forcing (25, 26). Our Bayesian fits generally parse the 140- to 300-year-long simulations into annual, decadal, and multicentennial modes, consistent with previous results (2527, 30), although results differ slightly across GCMs (Table 1 and table S1).

Table 1 Contributions of fast and slow modes to equilibrium and historical warming.

The posterior median and 5 to 95% CI values for the time scale, τ, and feedback factors, λ, of individual eigenmodes are provided, along with the fraction of warming contributed in equilibrium and from historical forcing as of 2011.

View this table:

Figure 1 illustrates the evolution of Earth’s energy imbalance with global warming in a representative GCM, and results from 23 additional GCMs are shown in figs. S2 to S5. Toward the limit of zero forced temperature response, our Bayesian inversion provides an estimate of the radiative forcing associated with a doubling of CO2. A broad range of 2.9 to 5.3 W/m2 [90% credible interval (CI)] across the ensemble reflects rapid initial temperature adjustments that limit the information that can be extracted from global mean values and stochastic radiative imbalances contributed by internal climate variability that obscure the forced response. Simulations using fixed sea surface temperatures (SSTs) can be used to independently constrain F, but these may give estimates that are biased low because land temperatures are allowed to evolve independently and tropospheric adjustment further opposes radiative forcing (31). Our maximum likelihood values of F are larger than the fixed SST estimates from the 10 models reported in Assessment Report 5 (AR5) of the Intergovernmental Panel on Climate Change (IPCC) (1) by a median of 0.1 W/m2, with only one estimate indicating a smaller value. Better determination of the forcing simulated within GCMs and the forcing historically experienced by the actual climate system would improve understanding of climate sensitivity (32).

Fig. 1 Evolution of top-of-atmosphere energy flux as a function of global temperature.

Annual average (brown dots) and 10-year running average values (black line) of top-of-atmosphere energy flux (TOA flux or H) and global-average temperature (T) are shown for a representative CMIP5 GCM (NorESM1-M) along with draws from the Bayesian fits (pink lines) and their median (red line). The distributions of radiative forcing (green line along y axis) and ECS (purple line along x axis) are estimated as part of the Bayesian fit. Curvature in the fit can be compared against a constant equilibrium feedback parameter that connects median radiative forcing to median warming (dashed gray line), where convexity indicates a decrease in the magnitude of λ(t) with time. Note that the line of constant λ used here connects median TOA and ECS values, whereas a regression line from the Gregory method (10) is used elsewhere. Also, note that both T and H are halved to facilitate comparison of these results using quadrupled CO2 forcing with the standard definition of ECS using doubled CO2. Ticks at top of graph indicate simulation year. NorESM1-M was chosen on the basis of having the minimum total deviation from the ensemble median in ECS, ICS, and F (see figs. S4 to S7 for the other 23 GCMs.)

The other limit indicated in Fig. 1 is toward full equilibration of top-of-the atmosphere radiation, but in which the approach of H toward zero typically involves an inflection. The radiative response, λ(t), or the negative of the slope of H with respect to T, typically decreases in magnitude, such that greater warming is required to achieve radiative equilibrium. For these simulations, equilibrium warming equates to ECS and has a median of 3.5°C across the ensemble, with a 5 to 95% CI of 2.2° to 6.1°C. The skewed distribution of ECS (33) can be seen in Fig. 1 as the geometric result of uncertainties in λ(t) having an asymmetric influence on the temperature at which top-of-atmosphere energy flux reaches zero. Individual GCM ECS estimates have maximum likelihood values similar to those reported by the IPCC (table S1). However, the 95th percentile of ECS increases from 4.5° to 6.1°C due to accounting for inflections in λ(t) more fully as compared to previous estimates (34, 35).

As simulations warm, the net feedback operating at any instant represents a weighted average across each mode in the Bayesian fitEmbedded Image(2)Estimates of sensitivity based on the historical record are obtained by dividing F by λ(t). Typically, λ is assumed constant, but we include t to emphasize that Eq. 2 indicates time dependence. We refer to these estimates as instantaneous climate sensitivity, ICS(t) = F/λ(t). ICS is sometimes referred to as effective ECS (8) when instrumental estimates are discussed but is distinct from other definitions of effective ECS used for GCMs (15). (Our point is that ICS may be an ineffective estimate of ECS.)

Equation 2 represents how different modes of warming cause variations in ICS over time. The slowest mode has a characteristic time scale of 350 years (CI, 180 to 960 years) and contributes a median of 44% (CI, 28 to 66%) of the total inferred equilibrium warming across the Bayesian fits. Historical forcing, however, only weakly excites this slow mode, such that it only contributes 3% (CI, 1 to 7%) of the historical warming (Fig. 2 and Table 1). This contribution is computed by applying the best estimate of effective historical forcing from 1750 to 2011 provided in IPCC AR5 (36) to the Bayesian fit associated with each GCM. This small response reflects that the slowest mode still remembers volcanic cooling episodes prevalent between 1750 and 1900 and that it is slow to realize warming from more recent increases in radiative forcing. In comparison, had the 2011 forcing of 2.2 W/m2 been continuously applied since 1750, the slow mode would instead contribute 29% (CI, 18 to 43%) of the warming. The structure of historical forcing is thus strongly biased toward sampling the faster modes of response with characteristic time scales of 0.8 years (CI, 0.2 to 2.6 years) and 9 years (CI, 4 to 37 years).

Fig. 2 Evolution of fast and slow modes under historical forcing.

Response of the CMIP5 GCMs to AR5 historical forcing from 1750 to 2011 (brown) and to an abrupt increase in 1750 to the 2011 forcing value of 2.2 W/m2 (blue). Responses are derived from the medians of our Bayesian fit, which permits for parsing ultrafast (mode 1, solid lines), fast (mode 2, dashed line), and slow (mode 3, dash-dot lines) contributions. Equilibrium contributions (green lines), referred to as committed warming in this context (43, 44), are also parsed according to mode for a net forcing of 2.2 W/m2. Mode 3 accounts for 44% of equilibrium warming but only 3% of present warming (see Table 1 for the median and 5 to 95% range for relevant parameters).

The faster modes of response have larger magnitudes of λ, with best estimates of 1.6 W/m2 per°C (CI, 0.4 to 4.0 W/m2 per°C) for the ultrafast mode and 1.4 W/m2 per°C (CI, 0.8 to 2.8 W/m2 per°C) for the fast mode, as compared to 0.8 W/m2 per°C (CI, 0.3 to 1.6 W/m2 per°C) for the slow mode. Larger magnitudes of λ indicate that less warming is required to achieve radiative equilibrium. Thus, sampling essentially only the fast modes in historical estimates biases ICS toward lower values than ECS. In particular, sampling ICS in a manner that is consistent with observational estimates leads to an ICS of 2.5°C (CI, 1.6° to 4.2°C), or 1°C cooler than ECS (Fig. 3).

Fig. 3 Equilibrium and instantaneous climate sensitivity distributions.

(A) Distribution of ECS from 5000 posterior draws of our Bayesian fit to each of 24 GCMs (indicated by colors). Aggregating across the posterior draws for all GCMs yields a median of 3.4°C and a 5 to 95% CI of 2.2° to 6.1°C. (B) Similar to (A) but for instantaneous climate sensitivity. ICS is obtained by applying AR5 historical forcing to our Bayesian fits and has a median of 2.5°C and a 5 to 95% CI of 1.6° to 4.2°C. The range of ICS values estimated in historical studies (vertical dashed black lines) (8) bracket the most likely GCM values, demonstrating consistency between observational and GCM results when they are appropriately compared.

It is also possible to explore the spatial patterns associated with each mode through regression onto the temperature variability expressed at individual grid boxes (29). As can be expected from differences in effective thermal inertia (14, 3739), the fast and intermediate modes project primarily onto continental regions, whereas the slow mode projects most strongly onto regions associated with oceanic upwelling (fig. S6). The implied pattern of transient warming is in good agreement with historical observations and model simulations (fig. S7). Furthermore, our inference of the slow eigenmode being relatively more powerful agrees with simulations indicating a strong net radiative feedback associated with warming in the Eastern Equatorial Pacific (16) and the Southern Ocean (17), as well as with simulations showing that prescribing historical SSTs yields a smaller net feedback relative to allowing upwelling regions to more fully equilibrate (18, 21). Our results thus indicate that the slow warming of the Eastern Equatorial Pacific and Southern Ocean are primarily responsible for the distinction between ECS and ICS.


A recent review of observationally based estimates of ICS shows a median of 2°C and an 80% range of 1.6° to 3°C, wherein the range is obtained by trimming the highest and lowest estimates from across the 10 available studies (8). This range of observationally inferred ICS values is contained within and centered on the maximum likelihood value of the distribution of ICS values that we estimate from applying historical forcing to the GCM eigenmode decomposition (Fig. 3B). Accounting for further issues involving forcing efficacy (40, 41) and temperature sampling methodology (42) might narrow the distributions and bring them into even closer agreement. Note that although historical GCM runs are consistent with instrumental temperature records (43), ICS cannot be directly determined from these runs because the radiative forcing is insufficiently known, generally having been diagnosed under the explicit assumption of constant radiative feedback (35).

Slow feedback contributions to warming only weakly manifest in the current climate system because the rise in greenhouse gas concentrations that occurred primarily in the past 50 years has forced the system for a short period relative to the multicentennial time scale of the dominant slow mode. This small manifestation makes it difficult, if not impossible, to accurately estimate ECS from historical observations and highlights the importance of paleoclimate records for the purpose of observationally constraining slow feedback processes (2). Looking forward, these results also highlight the fact that the warming current greenhouse that gas concentrations commit us to (44, 45) will evolve according to a slow mode of response that is distinct from present warming patterns.


Eigenmode decomposition was performed on global and annual average temperature, T, and top-of atmosphere energy flux, H, using a full Bayesian statistical inference methodology. Simulations involving an instantaneous quadrupling of CO2 were considered across 24 GCMs (table S1). These simulations have the highest signal-to-noise ratio of the available CMIP5 simulations, do not suffer from issues of collinearity between forcing and response, as in simulations with linear increases in forcing, and are not subject to uncertainties in the structure of the applied forcing, as in the historical GCM simulations; and the sustained forcing most strongly excite the slowest modes of response. In addition, values of ECS stated by the IPCC (1) are based on these abrupt quadrupling simulations (34, 35).

The constant forcing, denoted as 2F, permits simplification of Eq. 1 (25), toEmbedded Image(3)where the evolution toward an equilibrium warming of 2T involves three eigenmodes that have contributions of relative magnitudes, αn, and characteristic time scales, τn.

The radiative response to increased global average surface temperature is represented as evolving along the same eigenmodesEmbedded Image(4)toward a radiative equilibrium wherein R equals the radiative forcing, 2F. Fractional contributions of eigenmode n are given by βn. It follows that the radiative feedback associated with mode n isEmbedded Image(5)where F/T is the net equilibrium feedback. Previously described simple linear models for changes in λ(t) (14, 19, 23, 25) assumed three or fewer modes such that their temporal evolution can be fully described using Eqs. 3 to 5 (29).

The radiative response, R, is not directly available from GCM simulations. However, the top-of-atmosphere energy flux, H, could be related to the radiative response as H = FR, yieldingEmbedded Image(6)for the abrupt CO2 quadrupling experiment. Note that responses were corrected for GCM drift using the standard methodology (35), wherein drift values for T and H were computed in preindustrial control runs and subtracted from the quadrupling simulation.

Bayesian inference was performed by conditioning the parameters in Eqs. 3 and 6 on time series of T and H from the CMIP5 simulations. Variability that was not linearly related to the forced response in T and H, which can arise, for example, from internal modes of variability, such as El Niño, was represented as a normally distributed first-order autoregressive processes with autocorrelation parameters ρT and ρH and variances σ2T/(1 − ρ2T) and σ2H/(1 − ρ2H). This formulation was supported by the posterior residuals (figs. S8 and S9) that had a structure consistent with an autoregressive process and exhibited only minor deviations from a normal distribution when subjected to a Kolmogorov-Smirnov test that accounts for autocorrelation (46).

Uniform priors were prescribed for parameters defined on bounded intervals. Autocorrelations ρT and ρH were defined on the interval 0 to 1, and αn and βn were defined on the surfaces defined by the sum of points, Embedded Image and Embedded Image. Unbounded parameters were given weakly informative exponential priors that helped the sampling scheme converge more quickly by focusing the sampling away from implausible regions. T and F had exponential priors with means of 10°C and 10 W/m2, respectively; time scales τ1, τ2, and τ3 had exponential priors with means of 10, 100, and 1000 years, respectively; and SDs σT and σH had exponential priors with means of 10°C and 10 W/m2, respectively. The posterior maximum likelihood value for each parameter was consistent within 1% even if the means of all exponential priors were increased or decreased by a factor of 10, whereas the upper 95% value associated with each parameter was consistent within 10%. This low sensitivity to drastic changes in priors demonstrated that the simulation results from GCMs determined the posterior estimates.

Full Bayesian statistical inference was performed using the statistical modeling platform Stan (47). In particular, a joint posterior distribution was obtained for F, T, τn, αn, βn, σT, σH, ρT, and ρH using a Hamiltonian version of Markov chain Monte Carlo with a No-U-Turn Sampler (48). Posterior distributions for radiative feedbacks, λn, were computed from the joint posterior using Eq. 5. Results shown in the main text were from 5000 draws across five chains of 1000 realizations. Each chain was sampled after an initial burn-in period of 500 draws. To provide a complete representation of our results, we provided these draws from the joint posterior in data file S1.

ICS was inferred for the historical period by applying historical forcing to the estimated eigenmode response for each GCM. In particular, we applied the best estimate AR5 historical forcing (36), as used in previous studies (5, 7, 8), to draws from the posterior distribution of the eigenmode decomposition for each GCM to obtain the time series of T and H.

In the main text, λ(t) was computed as (FH)/T, where anomalies in forcing, heat uptake, and global temperature were computed as the departures in 2011 from preindustrial conditions, and ICS was obtained as F/λ(t). Note that values of F associated with each posterior draw could vary from the 3.7 W/m2 assumed in the AR5 estimate of historical forcing. To propagate this uncertainty, we multiplied F by F/3.7 for each draw before obtaining the values of H and T. We obtained a median ICS value across all draws and GCMs of 2.5°C (CI, 1.5° to 5.4°C). Inferences of “effective ECS” from the observational record, which we referred to as ICS, had been obtained using different intervals, but these alternate choices did not much influence our estimates. For example, computing λ(t) = (ΔF − ΔH)/ΔT between the intervals 1859 to 1883 and 1995 to 2011, following Otto et al. (5), we obtained an ICS of 2.6°C (CI, 1.6° to 4.2°C). Similarly, using intervals of 1860 to 1879 and 2000 to 2009, following Gregory et al. (3), we obtained an ICS of 2.6°C (CI, 1.6° to 4.1°C), and using the difference between values in 1955 and 2011, following Roe and Armour (4), we obtained an ICS of 2.3°C (CI, 1.5° to 4.2°C).


Supplementary material for this article is available at

Supplementary Text

fig. S1. Structure of residuals.

figs. S2 to S5. Figure 1 continued.

fig. S6. Spatial projection of the eigenmodes.

fig. S7. Historical temperature anomalies.

fig. S8. Autocovariance of temperature residuals.

fig. S9. Autocovariance of energy flux residuals.

table S1. Posterior parameter values for each GCM.

data file S1. Ensemble of posterior draws.

References (4952)

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: K. Armour, T. Laepple, K. McKinnon, A. Rhines, A. Stine, M. Tingley, C. Wunsch, and three anonymous reviewers provided helpful insights. CMIP5 data were obtained from the Earth System Grid Federation, the World Data Center for Climate in Hamburg, and the Geophysical Fluid Dynamics Laboratory. Funding: This work was funded by NSF grant 1304309. Author contributions: C.P. and P.J.H. designed the study and wrote the manuscript. C.P. performed the analysis. Competing interests: The authors declare that they have no competing financial 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 may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article