## Abstract

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

## INTRODUCTION

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 (*3*–*7*) 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 CO_{2} 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/m^{2} 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 (*3*–*8*, *10*), λ = (*F* − *H*)/*T*. *F* represents anomalies in downward radiative forcing, and the difference gives the upward radiative response, *R* = *F* − *H*. 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 CO_{2}, giving an estimated ECS of *F*_{2×}λ^{−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 (*11*–*18*). 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.

## RESULTS

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*, *22*–*24*), or any other linear system, using eigenmode decomposition (*25*)(1)where * represents the convolution operator. Each eigenmode of the temperature response, *T*_{n}, 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*_{2×}/*F*_{2×}. 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 (*25*–*28*) 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, *R*_{n} = λ_{n}*T*_{n}, 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*_{2×}. 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 CO_{2} 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 (*25*–*27*, *30*), although results differ slightly across GCMs (Table 1 and table S1).

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 CO_{2}. A broad range of 2.9 to 5.3 W/m^{2} [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*_{2×}, 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*_{2×} 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/m^{2}, 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*).

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 fit(2)Estimates of sensitivity based on the historical record are obtained by dividing *F*_{2×} 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*_{2×}/λ(*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/m^{2} 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).

The faster modes of response have larger magnitudes of λ, with best estimates of 1.6 W/m^{2} per°C (CI, 0.4 to 4.0 W/m^{2} per°C) for the ultrafast mode and 1.4 W/m^{2} per°C (CI, 0.8 to 2.8 W/m^{2} per°C) for the fast mode, as compared to 0.8 W/m^{2} per°C (CI, 0.3 to 1.6 W/m^{2} 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).

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*, *37*–*39*), 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.

## DISCUSSION

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.

## MATERIALS AND METHODS

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 CO_{2} 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 2*F*_{2×}, permits simplification of Eq. 1 (*25*), to(3)where the evolution toward an equilibrium warming of 2*T*_{2×} 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 eigenmodes(4)toward a radiative equilibrium wherein *R* equals the radiative forcing, 2*F*_{2×}. Fractional contributions of eigenmode *n* are given by β_{n}. It follows that the radiative feedback associated with mode *n* is(5)where *F*_{2×}/*T*_{2×} 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* = *F* − *R*, yielding(6)for the abrupt CO_{2} 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 σ^{2}_{T}/(1 − ρ^{2}_{T}) and σ^{2}_{H}/(1 − ρ^{2}_{H}). 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, and . 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*_{2×} and *F*_{2×} had exponential priors with means of 10°C and 10 W/m^{2}, 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/m^{2}, 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*_{2×}, *T*_{2×}, τ_{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 (*F* − *H*)/*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*_{2×}/λ(*t*). Note that values of *F*_{2×} associated with each posterior draw could vary from the 3.7 W/m^{2} assumed in the AR5 estimate of historical forcing. To propagate this uncertainty, we multiplied *F* by *F*_{2×}/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 MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/3/7/e1602821/DC1

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.

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.

- Copyright © 2017 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).