The sensitivity of the climate to CO2 forcing depends on spatially varying radiative feedbacks that act both locally and nonlocally. We assess whether a method employing multiple regression can be used to estimate local and nonlocal radiative feedbacks from internal variability. We test this method on millennial-length simulations performed with six coupled atmosphere–ocean general circulation models (AOGCMs). Given the spatial pattern of warming, the method does quite well at recreating the top-of-atmosphere flux response for most regions of Earth, except over the Southern Ocean where it consistently overestimates the change, leading to an overestimate of the sensitivity. For five of the six models, the method finds that local feedbacks are positive due to cloud processes, balanced by negative nonlocal shortwave cloud feedbacks associated with regions of tropical convection. For four of these models, the magnitudes of both are comparable to the Planck feedback, so that changes in the ratio between them could lead to large changes in climate sensitivity. The positive local feedback explains why observational studies that estimate spatial feedbacks using only local regressions predict an unstable climate. The method implies that sensitivity in these AOGCMs increases over time due to a reduction in the share of warming occurring in tropical convecting regions and the resulting weakening of associated shortwave cloud and longwave clear-sky feedbacks. Our results provide a step toward an observational estimate of time-varying climate sensitivity by demonstrating that many aspects of spatial feedbacks appear to be the same between internal variability and the forced response.
Forecasting global warming is one of climate science’s key challenges. As the atmospheric carbon dioxide concentration increases, the planet’s radiation of energy to space becomes less than its absorption of sunlight (Arrhenius 1896). This energy imbalance, the radiative forcing, warms the surface, setting off processes (radiative feedbacks) that close the imbalance, restoring the system to a new steady state. We call the global average of the radiative feedbacks the climate feedback [also called the climate feedback parameter (National Research Council 1979) or the thermal damping rate (Dessler 2013)]. The total warming in response to a given increase in CO2 is thus determined by the resulting radiative forcing and the climate feedback (National Research Council 1979). The rate of warming also involves the thermal inertia of the surface, mostly due to oceanic heat uptake (Gregory et al. 2002). Uncertainty in the climate feedback contributes the most to uncertainty in future warming (Otto et al. 2013; Lewis and Curry 2015; Lutsko and Popp 2019), in part because of the inverse relationship between feedback and sensitivity (Roe and Baker 2007).
Directly simulating radiative feedbacks is difficult primarily because cloud feedbacks depend on small-scale processes (Wetherald and Manabe 1988). Alternatively, the climate feedback can be inferred from observations, either by solving for it using the observed warming, observed deep ocean heat uptake, and simulated radiative forcing (Gregory et al. 2002; Otto et al. 2013) or by analyzing how the planet’s energy imbalance changes as the surface temperatures varies month-to-month or year-to-year (Forster and Gregory 2006; Murphy et al. 2009; Dessler 2010; Cox et al. 2018; Lutsko and Takahashi 2018; Jiménez-de-la-Cuesta and Mauritsen 2019; Libardoni et al. 2019). These observational methods often assume that the climate feedback is constant, but many studies have shown that it typically changes with time in simulations (e.g., Murphy 1995; Watterson 2000; Senior and Mitchell 2000; Armour et al. 2013; Jonko et al. 2013; Andrews et al. 2015). While the temperature dependence of feedbacks can cause this to occur under sufficient (and likely strong) warming (Meraner et al. 2013; Bloch-Johnson et al. 2015), the change occurs even after relatively small amounts of warming (e.g., Armour et al. 2013; Andrews et al. 2015; Rugenstein et al. 2016). Since warming in different regions sets off radiative feedbacks of different strengths, the inconstancy of the climate feedback is likely caused by the change in the spatial pattern of warming with time (Winton et al. 2010; Armour et al. 2013). Since the temperature pattern associated with internal variability differs from the forced response, we should expect the climate feedback associated with each to differ (Dessler 2013; Colman and Hanson 2017), and in fact the climate feedback appears to vary across the historical record (Gregory and Andrews 2016; Fueglistaler 2019). The climate feedback may vary between historical and future warming (Zhou et al. 2016; Armour 2017; Proistosescu and Huybers 2017; Andrews et al. 2018), although the importance of this effect may be modest (Lewis and Curry 2018).
Recent modeling work has explored a new framework in which the climate feedback is a linear combination of radiative feedbacks associated with different regions of the surface, weighted by the temperature change in each region (Zhou et al. 2017; Dong et al. 2019). This assumes that the spatial radiative feedbacks themselves are constant, with only the map of surface temperature change evolving. This paper explores a corollary: since internal variability creates an ever-changing pattern of surface temperature and top-of-atmosphere radiative imbalance, a sufficiently long record of this variability should exhibit the behavior of these spatial radiative feedbacks. In this paper, we propose and evaluate a multiple regression (MR) method to estimate the spatial radiative feedbacks of six atmosphere–ocean general circulation models from control simulations, which we compare to existing methods for estimating feedbacks from internal variability (section 2). We do so in spite of the known bias in regression methods related to stochastic variation in top-of-atmosphere fluxes (Spencer and Braswell 2008, 2011; Choi et al. 2014; Proistosescu et al. 2018). We test the method by convolving the estimated spatial feedbacks with warming patterns from forced simulations performed with the respective models (section 3), assessing the method’s accuracy in recreating aspects of the forced response. We discuss insights the MR method provides into climate dynamics, such as the competing nature of local and nonlocal cloud feedbacks (section 4), and summarize our findings (section 5).
2. Illustrating the MR method with a conceptual model
In this section, we present a method for predicting spatial feedbacks from records of unforced variability using multiple regression. We first set up a conceptual climate model designed to illustrate the method and capture some features of the complex climate models discussed in section 3. This conceptual model has two regions of equal area. In each, the change in surface temperature (Ti) is proportional to the net energy gain of that region, which is the sum of the net downward top-of-atmosphere (TOA) radiative flux (Ni), the net gain from horizontal energy transport from the atmosphere and ocean combined (−H in region 1, H in region 2), and additional random forcing (Fsurf,i):
where ci is the surface thermal inertia associated with region i. This model can be re-expressed in terms of anomalies relative to an initial equilibrium state, so that we consider instead of Ti, Ni, H, and Fsurf,i. We assume that heat transport is proportional to the temperature gradient between the two regions:
Changes in a region’s top-of-atmosphere radiative fluxes are caused by radiative feedbacks (λi,j, which represents the influence of surface temperature in region j on the net TOA flux in region i), radiative forcing due to changes in a forcing agent such as an increase in CO2, and radiative forcing due to random atmospheric fluctuations that occur independently of surface temperature (FTOA,i):
The terms λ1,1 and λ2,2 are local radiative feedbacks, while λ1,2 and λ2,1 are nonlocal radiative feedbacks (where our sign convention ensures that a negative λ implies a negative, stabilizing feedback).
Nonlocal radiative feedbacks (Rugenstein et al. 2016; Zhou et al. 2017; Po-Chedley et al. 2018; Dong et al. 2019) are changes in a region’s top-of-atmosphere flux that occur due to changes in surface temperature elsewhere, independent of local surface temperature changes. For example, in Fig. 1, regions 1 and 2 represent the convecting and subsiding branches of an overturning cell respectively. Surface warming in region 1 propagates vertically, warming region 1’s free troposphere, and then horizontally into the free troposphere of region 2, increasing H′. Region 2 now has a warmer troposphere, which radiates more, decreasing . The resulting horizontal advection may also increase the humidity of region 2’s free troposphere, increasing . Assuming region 2 has a subsidence-induced boundary layer inversion, its low cloud cover could also increase, causing a further decrease in . All of these changes in occur independently of any changes in , and conspire to make λ2,1 positive or negative.
We note that an increase in H′ will also increase directly [Eq. (1); Feldl and Roe 2013b]. While this latter effect is connected to nonlocal radiative feedbacks in that both occur due to horizontal fluxes of heat and moisture, the two effects are different, and can disagree in the sign of the resulting surface warming, as demonstrated by the above example. While the influence of H′ on surface temperature is important for understanding the evolution of the spatial pattern of warming, in this paper we are focused only on the influence of surface temperature on TOA radiative fluxes, and so we focus on nonlocal radiative feedbacks.
Suppose that region 1 has a weak positive local feedback, λ1,1 = 0.5 W m−2 K−1, and a stronger negative nonlocal feedback, λ2,1 = −2 W m−2 K−1, so that the combined feedback from region 1 is λ1 = −1.5 Wm−2 K−1 (red solid line, Fig. 2b). We also assume that the surface temperature of the subsiding region 2 has no net effect on TOA fluxes, so that λ1,2 = λ2,2 = λ2 = 0 W m−2 K−1 (gray solid line in Fig. 2b). We assume that region 2’s thermal inertia is much larger than region 1’s, representing more ocean heat uptake in this region (see the appendix for details).
We define the global climate feedback λ to be the dependence of the globally averaged net TOA flux on the globally averaged surface temperature, that is
where , and a bar over a quantity indicates the global average of that quantity. We do not have to use an anomaly for because is 0 in equilibrium. Note that even though the spatial feedbacks Λ are constant, the global feedback λ can change with time because of the evolving spatial pattern of warming (dT/dt)(t).
We perform two 5000-yr experiments: a “control” experiment, where all variations in and are due to random forcing at the surface and TOA , and an “abrupt4×” experiment in which the time series also respond to an initial step forcing akin to a quadrupling of CO2 concentration .
For the abrupt4× simulation, the climate feedback changes significantly around year 20. We therefore define two forced feedbacks, λ4×,early and λ4×,late, which are the slopes of the linear regressions of against taken over years 1 to 20 and years 21 to 5000 respectively (Fig. 2c). Before these regressions are taken, we average each annual time series (gray dots) over roughly exponentially increasing time periods (colored dots). The change in feedback between the periods is Δλ4× ≡ λ4×,late − λ4×,early.
We seek a method to predict λ4×,early, λ4×,late, and Δλ4× given and (internal variability) and (the spatial pattern of warming). The simplest method would be to regress annual averages of against to get the resulting regression slope λcontrol (the slope of the blue line in Fig. 2a) and to assume that λ4×,early = λ4×,late = λcontrol (Forster and Gregory 2006; Murphy et al. 2009; Dessler 2010). We call this the “global” method because it uses information about changes in global surface temperature only.
The radiative feedbacks associated with temperature change induced by random forcing (i.e., Fsurf and FTOA) differ from those induced by uniform greenhouse forcing (Dessler 2013; Colman and Hanson 2017; Proistosescu et al. 2018). Our conceptual model illustrates how this can arise from spatial variation. Since the thermal inertia in region 2 is larger, most of the temperature variability occurs in region 1, so that λcontrol is weighted toward the feedbacks associated with this region (λcontrol ≈ λ1,1 + λ2,1). The spatial pattern of warming in the forced response is initially dominated by region 1 as well, once more because it has the lowest thermal inertia. As a result, the global method predicts λ4×,early well (see Figs. 2c,d). However, the global method always predicts Δλ4x = 0, as it assumes a constant λ. Since warming moves to region 2 over time and λ1,2 + λ2,2 > λ1,1 + λ2,1, Δλ4× is positive. As a result, the global method underpredicts the warming of the abrupt4× simulation by about 1.5 K (Fig. 2c). To address this shortcoming, we need a method that accounts for the spatial variation of feedbacks.
The “local” method is a commonly used method [Boer and Yu 2003b; Crook et al. (2011); see the local method in Feldl and Roe (2013a), Brown et al. (2016), and Trenberth et al. (2015)] for estimating spatial feedbacks. In this method, we construct where λi,local is the result of regressing against . Taking the dot product of λlocal with then provides an estimate of that we can use to estimate λ4×,early, λ4×,late, and Δλ4×.
This method assumes all radiative feedbacks are local, while allowing for the nonlocal effects of heat transport (Feldl and Roe 2013b). However, if there are nonlocal radiative feedbacks, then the local method can miss or conflate their effects. In region 1, estimates of λ1,local tend toward λ1,1 = 0.5 W m−2 K−1 (dotted red line, Fig. 2b), missing the negative nonlocal feedback λ2,1. Since the early period is dominated by warming in region 1, the local method overestimates λ4×,early (where “overestimating” implies that the estimate of λ4×,early is more positive than the true value, even if both are negative, resulting in an overestimate of the sensitivity). On the other hand, tends to be positively correlated with , due to heat transport, while tends to be anticorrelated with because λ2,1 is negative. As a result, the local method predicts that λ2,local is negative (dotted gray line, Fig. 2b), even though has no net influence on N. Since contributes more to warming over time, the local method incorrectly predicts a more negative feedback (Figs. 2c,d). Similar discrepancies can occur when local feedbacks are used to diagnose feedbacks in GCMs, which may explain instances when the local method fails to predict feedback changes properly (Rose et al. 2014). We need a method that includes nonlocal feedbacks while accounting for correlation between temperature in different regions.
We propose a multiple regression method (MR method), which estimates the local and nonlocal feedbacks associated with (i.e., the influence of and on ) by regressing against both regions simultaneously:
In least squares multiple regression, λi,j,MR is the same as the slope of the regression of against , where the star indicates that each time series is the residual after regressing against the surface temperatures in all non-j regions (see the appendix). This removes the effect of correlations between surface temperature in different regions giving spurious feedbacks, as with λ2,local above. Multiple regression has been used to estimate other surface temperature-dependent feedbacks from internal variability, although not radiative feedbacks (Liu et al. 2008; Li et al. 2012; Li and Forest 2014; Liu et al. 2018). The dashed lines in Fig. 2b show that, given sufficient time, the MR method predicts the local and nonlocal feedbacks in each region, so that when we multiply the full matrix of estimated spatial feedbacks by to estimate Nabrupt4× (t), the resulting estimates λ4×,early, λ4×,late, and Δλ4× are accurate (Figs. 2c,d). Therefore, for this example, the MR method is able to account for the difference in climate feedback between internal variability and the forced response.
Random fluctuations in N influence T via planetary energy gain at the same time that T influences N via radiative feedbacks. As a result, T will tend to lag N with a positive correlation, while N will lag T with a negative correlation, so that regressions taken without a lag will be biased toward 0 (Spencer and Braswell 2008, 2011; Choi et al. 2014; Proistosescu et al. 2018). This issue does not occur for random forcing at the surface, which only affects N indirectly through radiative feedbacks. Therefore, the more stochastic forcing that occurs at TOA (FTOA) as opposed to the surface (Fsurf), the more the regression of N versus T will overestimate the true radiative feedback. For the example in Fig. 2, Fsurf,1 and Fsurf,2 are white noise with variance 20 W2 m−4, while FTOA,1 and FTOA,2 are white noise with variance 5 W2 m−4. Figure S1 in the online supplemental material shows a case where these variances are 10 and 15 W2 m−4 respectively, with the result that all three regression methods overestimate λ4×,early and λ4×,late, while underestimating Δλ4×. In other words, given sufficient random TOA forcing, regression estimates of spatial feedbacks will be biased. We consider this bias in discussing our results in the next section.
It should be mentioned that Proistosescu et al. (2018) model ENSO variability as a distinct additional mechanism by which N and T mutually influence each other, which similarly leads to overestimates of λ from regression-based methods. As part of their model, they assume that T influences N with a lag of about three months. Since this is beyond the time scale of most atmospheric processes, we assume that this feedback propagates in part through the ocean, so that the atmospheric component may still operate through the same spatial feedbacks that operate under other forms of variability and under the forced response (e.g., it could occur due to a “tropical atmospheric bridge” mechanism; Klein et al. 1999).
3. Using the MR method on AOGCMs
To test the methods discussed above on atmosphere–ocean general circulation models (AOGCMs), we use simulations from LongRunMIP, an archive of fully coupled millennial-length simulations of complex climate models (Rugenstein et al. 2019). We chose the six models with millennial-length control and abrupt4× simulations for which we have monthly output. Details of these models and simulations are given in Table S1 in the online supplemental material.
We alter the three methods from section 2 to reflect the more complex nature of AOGCMs:
CO2 forcing can lead to atmospheric changes that are independent of surface warming. These “adjustments” to forcing occur mostly within the first year (e.g., Gregory and Webb 2008). We remove this year from our analysis, redefining our early period to be years 2 to 20.
For AOGCMs, there are more than two regions with distinct behaviors. Dividing our models into n regions, Eq. (7) becomes
giving a system of n equations
where Λ is a matrix of feedbacks λi,j. Each equation in this system has n − 1 degrees of freedom, so n must be smaller than the length of the control simulation, and preferably much smaller given the significant spatial correlation of surface temperature. For simplicity, we divide the surface equally in latitude and longitude, although this may miss features of the climate system. Since our control simulations last at least 1000 years (Table S1), we use a 15° by 15° grid, giving 288 regions (Fig. 3).
Circulations, and therefore radiative feedbacks, change with season. Thus, we compute feedbacks for each season individually, first by averaging all monthly time series into seasonal time series (where the seasons are DJF, MAM, JJA, and SON), and then performing a separate regression for each season [e.g., all DJF values of against all DFJ values of ], creating a set of four feedbacks. We multiply each month of by the relevant seasonal feedback, and take the annual average to estimate . We compare seasonal averages to other approaches in Tables S2 and S3. While seasonal averaging tends to reduce the error in the MR method, the qualitative behavior of the different methods is not affected by the choice of time averaging.
Figures 3 and 4 show versus of the control and abrupt4× simulations of the six models respectively. Figure 4 also shows estimated using the three methods, assuming that each estimate starts with the true value of at year 2. The solid lines in Fig. 4 are local regressions of against performed using LOESS (locally estimated scatterplot smoothing; Cleveland and Devlin 1988; see the appendix for more detail). We can use the slopes of these lines mapped against the time series of to estimate feedbacks as a function of time (lines in Fig. 5).
Although there is a range of feedback values between models, all six forced simulations have a feedback that gets less negative with time (black lines), consistent with past results for similar models (Andrews et al. 2015). The MR method (green lines) matches or overestimates the feedback value, with this error tending to decrease with time. This error can range from ~1 W m−2 K−1 for the early years of CESM1.0.4 and GISS-E2-R (i.e., at least half of the feedback strength itself) to roughly 0 for HadCM3L. The MR method correctly predicts that the feedback gets less negative with time, although for some of the models it underestimates the magnitude of the change.
The global method (blue) overestimates the early feedback. Since the global method is agnostic about the pattern of surface warming, the predicted feedback is mostly constant except for small differences due to changes in the seasonal distribution of warming and in seasonal feedbacks (e.g., the early years of HadCM3L). As a result, as the true feedback increases with time, it becomes more positive than the global estimate for half the models. For some models, this allows the global method to more accurately forecast the equilibrium warming than the other methods, albeit due to compensating errors in the early and later periods (i.e., CESM1.0.4 and MPI-ESM1.2 in Fig. 4).
The local method (orange) predicts a positive feedback for all models except GISS-E2-R, implying a climate unstable to external forcing, and does not predict the increase in feedback with time seen in all models.
The dots in Fig. 5 represent estimates of λ4×,early and λ4×,late (feedbacks before and after year 20; see the appendix for details). We visualize the estimates of these feedbacks and their difference using a scatterplot (black dots in Fig. 6), as in Fig. 2d. The global and MR methods perform similarly for λ4×,early and λ4×,late, while the MR method gets closer to accurately predicting Δλ4×, consistent with the discussion around Fig. 4 and reflected by the root-mean-square errors in Table 1 (for feedback values for all models and components, see Tables S7 and S8).
The terms N′ and λ can be expressed as the sum of shortwave (SW) and longwave (LW) terms, which can be separated in turn into clear-sky (fluxes recalculated as if no clouds were present) and cloud terms (the residual of total and clear-sky terms; cloud feedbacks defined this way may include changes in cloud masking rather than in clouds themselves; Soden et al. 2004).
Examining these components individually shows that the error in λ4×,early in the MR and global methods is due primarily to SW cloud feedbacks (red markers in Figs. 6a and 6b). Both the MR and global methods have smaller errors in λ4×,late (Figs. 6d,e), but for the MR method this is caused by a reduction in the error in SW cloud, while for the global method this is due to offsetting errors in the SW and LW cloud feedbacks (see also Table 1). Cloud feedbacks are similarly the cause of the local method’s large overestimation, while the local method outperforms the other methods at predicting the primarily local SW clear feedback (Table 1). Note that the global method has a relatively small error for the LW clear feedback, consistent with Lutsko and Takahashi (2018). The increase in feedback with time (Δλ4×) and the variation in this increase between models is driven by the SW cloud feedback (Figs. 6g–i). The MR method has the smallest error in estimating Δλ4×, with this error tending to be an underestimate. Figures S2–S5 show feedback time series plots for all component fluxes.
All methods examined contain some degree of error. We can find the geographic source of these errors by looking at the true and estimated normalized change in (multimodel mean in Fig. 7; errors in the multimodel mean and for individual models in Figs. S6–S8), calculated by taking the finite difference in between the first and last part of the indicated time period, where each part contains similar amounts of warming (see the appendix). The difference is normalized by the global temperature change, allowing intermodel comparison. For the global method, we make this estimate by regressing against [the “global” method in Feldl and Roe (2013a) and the “local contribution” in Boer and Yu (2003a,b), Crook et al. (2011), Zelinka et al. (2012), and Andrews et al. (2015)] and using this as the predicted normalized change in .
The MR method does quite well at recreating the multimodel spatial pattern of TOA flux change, both for net and component fluxes (Figs. S9–S12), with the exception of regions south of 30°S and the North Atlantic. The MR method also overestimates the change in these regions in individual models (Figs. S6–S8). The error in these regions has contributions from all component fluxes, foremost the SW cloud feedback (for multimodel mean component flux errors, see Figs. S13–S17). For all periods, models, and fluxes except for SW clear sky (which is primarily a local feedback), the MR method outperforms the other two methods when scored by the area-weighted root-mean-square error (Table 2; for comparison with annual or monthly approaches, see Table S3; for values for individual models, see Table S4; for details on the error metric, see the appendix). Specifically, the global method has large compensating errors, especially in the tropics, and the local method overestimates the change almost everywhere (Figs. S6–S8).
There are several potential explanations for the MR method’s overestimate for TOA fluxes south of 30°S and over the North Atlantic. These may be regions where there is significantly more stochastic forcing at TOA than at the surface, resulting in a similar overestimation to that discussed in section 2 and shown in Fig. S1. Alternatively, the spatial feedbacks that influence N′ in these regions may be nonlinear, either in that they change in value as the world warms (e.g., a reduction in the strength of the SW clear feedback once sea ice melts), or the effect of warming in different regions combines nonlinearly, as might occur in response to circulation changes such as a shift in the midlatitude jet; or surface fluxes may influence N′ there independently of surface warming. Further research is needed to diagnose this error.
In spite of this overestimate, the MR method can be used to explain the multimodel forced TOA flux response for roughly three-quarters of Earth using feedbacks estimated from internal variability (see Table S5 and S6, which show the same error metrics as Tables 1 and 2, using only TOA fluxes north of 30°S). We now discuss the spatial feedbacks estimated by the MR method, as well as some of their implications.
We first test if the spatial feedbacks estimated using the MR method exhibit behavior broadly consistent with physically modeled feedbacks. The ith column of Λ represents the change in N′ from warming in region i. Zhou et al. (2017) performed fixed-SST experiments with the CAM5 model where the temperature in region i was perturbed. The top row of Fig. 8 shows spatial cloud feedbacks for three representative regions calculated using this approach. The bottom row shows the multimodel and multiseason mean response for warming in similar regions estimated by the MR method. For both approaches, warming in the extratropics or in regions of tropical subsidence produces cloud feedbacks that are mostly local and positive, while warming in tropical convecting regions has significant nonlocal feedbacks that are mostly negative. Since the models, region sizes, and degree of perturbation differ, the details and magnitudes of the feedbacks differ. Further, the fixed-SST method allows land temperatures to evolve freely, so that regions that have significant nonlocal effects, like tropical convecting regions, can cause large changes in TOA fluxes over land (Fig. 8b). The MR method is able to estimate land feedbacks directly, so that TOA flux changes due to land warming are not included in these tropical convecting feedbacks (Fig. 8e). See also Fig. 4 in Dong et al. (2019).
The top left panel of Fig. 9 shows a map of the multimodel and multimonth mean spatial feedbacks estimated by the MR method: the change in caused by warming in each region divided by that region’s fractional area (so that smaller polar regions do not have artificially smaller feedbacks). Spatial feedbacks are strongly negative in regions of tropical convection (e.g., Indonesia and Central America) and are mostly positive over the tropical oceans in regions of atmospheric subsidence as well as much of the extratropical oceans, in keeping with the examples from Fig. 8. These strongly negative feedbacks are robust when feedbacks are recalculated using just the first or second half of the control simulations (Figs. S18–S22), although outside these regions there is some noise, with the sign of roughly a third of the net feedback cells differing between the first and second halves. The variation in the spatial pattern is largely determined by the SW cloud feedback (bottom left panel, Fig. 9; for all flux components, see Figs. S19–S22).
a. Local and nonlocal feedbacks
The MR method allows us to split spatial feedbacks into local (the diagonal elements of Λ, giving the influence of warming on TOA fluxes directly overhead) and nonlocal components (the off-diagonal elements of Λ), and to calculate the local and nonlocal components of the map of spatial feedbacks (middle and right columns of Fig. 9 respectively). We note that the devision between local and nonlocal feedbacks depends on grid resolution, with local feedbacks in coarser grids incorporating more nonlocal processes. For the grid considered in this paper, the local feedback is positive almost everywhere, due to cloud feedbacks (Figs. S21 and S22): in the tropics and in subtropical subsiding regions, local warming reduces lower tropospheric stability, leading to a loss of low clouds and a positive SW cloud feedback (Klein and Hartmann 1993; Wood and Bretherton 2006; Zhou et al. 2017; Dong et al. 2019). This result holds for each AOGCM except for GISS-E2-R, which lacks a positive local SW cloud feedback (see Fig. S24 and Table S8). For most models, there is a partially compensating negative local LW cloud feedback in tropical convecting regions, possibly due to an iris effect (Lindzen et al. 2001; Mauritsen and Stevens 2015). Outside of the tropics, there is a positive local LW cloud feedback, possibly associated with an increase in middle and high cloudiness as convection increases (Zelinka et al. 2012).
Positive local feedbacks provide an explanation for observational studies that use the local method to predict spatial feedbacks, finding that they are positive over much of Earth and in the global mean (Brown et al. 2016; Trenberth et al. 2015). For example, the multimodel mean feedbacks estimated using the local method (top middle panel, Fig. S23) resemble the feedbacks in the upper right panel of Fig. 10 in Trenberth et al. (2015). While local method feedbacks can differ from the local component of MR method feedbacks due to correlation between temperature in different regions as discussed in section 2, the observational studies provide evidence that real world local feedbacks are substantially positive. If we use the MR method to estimate the local components of λ4×,early and λ4×,late (Table S8), we get positive values for all models except GISS-E2-R. For these models, the mean estimated local feedback is 3.37 W m−2 K−1 for the early period and 3.13 W m−2 K−1 for the late period (Table S8).
The MR method implies that in the absence of negative nonlocal feedbacks, five out of six of these AOGCMs would be unstable to radiative forcing, even accounting for the dominant stabilizing Planck feedback. The MR method predicts that there are strongly negative nonlocal feedbacks coming from regions of tropical convection (upper right panel, Fig. 9), largely due to the SW cloud feedback (lower right panel). This is consistent with tropical convecting regions behaving similarly to region 1 of the conceptual model from section 2: surface warming in the convecting tropics propagates throughout the tropical free troposphere, increasing the temperature aloft while leaving surface temperatures alone. This increases the lower tropospheric stability, and thus low cloud cover (a negative SW cloud feedback), as well as the troposphere’s outgoing longwave radiation (a negative LW clear feedback) (Rose and Rayborn 2016; Andrews and Webb 2017; Ceppi and Gregory 2017; Klein et al. 2017; Zhou et al. 2017; Dong et al. 2019). Note that incorporating these nonlocal interactions changes both local and total values of the LW clear feedback, giving different values than studies that analyze this feedback purely locally (e.g., Koll and Cronin 2018).
For the five models with positive local components, the average nonlocal component of the abrupt4x feedbacks is −4.21 W m−2 K−1 for the early period and −3.69 W m−2 K−1 for the late period (Table S8), so that the net forced climate feedback is a small residual between competing local and nonlocal feedbacks, with local and nonlocal feedbacks strongly anticorrelated between different models (Table S8; the correlation coefficient for early period non-GISS-E2-R local versus nonlocal feedbacks is −0.96, and for late is −0.98). A modest shift in the relative strength of these feedbacks (e.g., due to a shift in circulation) could lead to large changes in climate sensitivity; an increase in the local feedback of only a third would be enough to make these AOGCMs unstable (local and nonlocal feedbacks differ by ~1 W m−2 K−1, which is on average roughly a third of the magnitude of the local feedback for the non-GISS-E2-R models). Additional research is needed to understand what mechanisms cause the anticorrelation between local and nonlocal feedback strength, and whether we expect this cancellation to hold in different climate states. Given that the local/nonlocal cancellation does not hold in all contexts—for example, the nonlocal feedback’s seasonal cycle has a larger amplitude and is more latitudinally constrained than the local feedback’s seasonal cycle (Fig. S25)—it is unlikely that this cancellation is purely a statistical artifact. Our findings have a bearing on exoplanet research, as they suggest that it may be harder to have a cloudy atmosphere with a stable climate than previously thought (Leconte et al. 2013), potentially reducing the chance of finding habitable worlds.
b. The cause of the increase in climate feedback over time
For all six models, the change in feedback with time (Δλ4×) is positive, primarily because of the SW cloud feedback, and secondarily the LW clear feedback (Fig. 4 and Table S7). The MR method gets the correct sign of Δλ4× but underestimates this increase for each model, once more primarily due to the SW cloud feedback (Table S8).
We can estimate how much the change in the spatial pattern of warming with time (Fig. 10a) contributes to Δλ4× by multiplying this change by the MR estimate of the spatial pattern of feedbacks for each flux component (Fig. 9; see also Figs. S18–S22). The resulting maps show the contribution of the change in warming pattern to the change in feedback (Figs. 10b–f).
The MR method identifies two main latitude bands that contribute to the increase in feedback with time: the tropics, whose convecting regions increase the SW cloud and LW clear feedbacks [less warming in these regions reduces the role of the strongly negative nonlocal feedbacks discussed above, consistent with Andrews and Webb (2017), Ceppi and Gregory (2017), Dong et al. (2019), and Fueglistaler (2019)]; and the Southern Ocean, which increases the SW clear feedback (due to the delayed warming in this region leading to the delayed melting of sea ice). The MR method estimates that the LW clear sky and SW cloud feedback have offsetting negative contributions in the Southern Ocean. While the LW clear sky offset is consistent with the total change in the LW clear feedback being small, and with the LW clear TOA flux change getting more negative in the Southern Ocean due to a more strongly negative local feedback (zonal figures in the top row of Fig. S19), the change in the SW cloud TOA flux is too negative in this region (lower left panel of Fig. S17), suggesting that the SW cloud negative contribution is an error, and is likely the reason for the MR method’s underestimate of Δλ4x.
While the exact evolution of temperature patterns in the tropics in AOGCMs may be incorrect due to cold-tongue biases (Seager et al. 2019), our findings match with Dong et al. (2019) in that as long as the feedbacks in tropical convecting regions are far more negative than anywhere else, the delayed warming in regions of ocean heat uptake will ensure an increase in sensitivity over time. Observational evidence suggests that depends on tropical midtropospheric temperatures (Dessler et al. 2018; Ceppi and Gregory 2019; Fueglistaler 2019), supporting our argument that a reduction in the share of surface warming occurring in the tropical convecting regions that set these temperatures likely influences Earth’s sensitivity.
The global climate feedback, one of the key parameters in determining future climate change, is inconstant in part because radiative feedbacks vary spatially. The MR method estimates these spatial feedbacks from records of internal variability, and improves upon existing methods for doing so by incorporating both local and nonlocal radiative responses to surface warming. For the six atmosphere–ocean general circulation models studied, the spatial feedbacks estimated by the MR method applied to the pattern of surface warming recreate the spatial pattern of top-of-atmosphere flux response to forcing more accurately than existing methods, as well as providing better estimates of the change in feedback with time. The method consistently overestimates the change in TOA flux over the Southern Ocean and North Atlantic, and so overestimates the sensitivity. The method finds that that there are significant negative nonlocal feedbacks associated with regions of tropical convection, and that the reduction in the share of warming that occurs in these regions over time contributes to an increase in the global feedback with time in these models, consistent with recent studies (Andrews and Webb 2017; Ceppi and Gregory 2017; Dong et al. 2019; Fueglistaler 2019).
The MR method finds that five of the six AOGCMs have strongly positive local cloud feedbacks countered by strongly negative nonlocal cloud feedbacks. These positive local feedbacks may explain why studies that use local regressions to estimate spatial feedbacks from observed internal variability find that they are on average positive (Brown et al. 2016; Trenberth et al. 2015). While the AOGCMs exhibit an anticorrelation between local and nonlocal feedbacks, a small relative shift in the balance between these feedbacks could cause large changes in sensitivity, and such shifts may be relevant for paleoclimate or future warming. Given the large magnitudes associated with these local and nonlocal cloud feedbacks, it may be harder for cloudy exoplanets to have stable atmospheres, reducing the chances of finding habitable worlds.
Spatial feedbacks estimated from observations could potentially improve warming forecasts and serve as emerging constraints on AOGCMs. The success of the MR method for most fluxes and regions of Earth (with the important exception of Southern Ocean cloud feedbacks) suggests that many of the spatial feedbacks at work under global warming are observable under internal variability. Challenges remain to applying the MR method to observations. We would need to reduce the information necessary to fit our statistical model to be less than the length of the satellite record; to remove changes in forcing from records of top-of-atmosphere fluxes; and to account for systematic biases in the observations themselves. We would also need to account for regions of Earth and states of the climate where the MR method is biased, such as for Southern Ocean cloud feedbacks. Furthermore, since spatial feedbacks are just one link in the coupled energy balance of the climate, we would need complementary theory to complete the forecast of future warming, particularly its spatial pattern. Still, our results suggest that the processes that will determine the sensitivity in both the near and far future may be observable today.
This paper was inspired by conversations with Kyle Armour, Cristian Proistosescu, Daniel Koll, and Elizabeth Moyer and was improved further by discussions with B. B. Cael, Yue Dong, Jonathan Gregory, Malte Jansen, Thorsten Mauritsen, Brian Rose, Hansi Singh, M Wu, and Chen Zhou. The research would not have been possible without the efforts of the contributors to the LongRunMIP project. The authors thank Chen Zhou for permission to use the top row in Fig. 8. We acknowledge support from the National Science Foundation under NSF Award 1623064. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant Agreement 786427, project “Couplet”). Our work was completed with resources provided by the University of Chicago Research Computing Center, with special thanks to Hossein Pourreza. We also thank three anonymous reviewers for their insightful comments and our editor, Timothy Delsole, for his guidance throughout the publication process.
Data and Methods
a. Data/code access
b. Matrix and vector notation
Note that in the main body of the text, time is treated as continuous, so that time series are written as functions [e.g., T(t) is the evolving spatial pattern of warming]. Since the appendix documents the calculations we have employed, it treats time as discrete, and so time is instead treated as an additional dimension (e.g., T is the evolving spatial pattern of warming). Therefore, a vector in the main body of the text refers to a spatial pattern, while a vector in the appendix refers to a time series of a scalar value (such as a global average).
c. Conceptual model
The conceptual model is a system of stochastic differential equations:
The thermal inertia ci is defined as miρcp, where ρ and cp are the density and specific heat of ocean water respectively, and mi is an equivalent mixed layer depth; m1 is 50 m, and m2 is 1000 m. are both 0 W m−2 (8 W m−2) for the control (abrupt4×) simulation. λ1,1 = 0.5 W m−2 K−1, λ2,1 = −2 W m−2 K−1, λ1,2 = λ2,2 = 0 W m−2 K−1, and γ = 2 W m−2 K−1. The terms Fsurf and FTOA are white noise processes. In the example shown in Fig. 2, the variance of Fsurf,1 and Fsurf,2 is 40 W m−2 and the variance of FTOA,1 and FTOA,2 is 5 W m−2, while for the example in Fig. S1, the variance of Fsurf,1 and Fsurf,2 is 10 W m−2 and the variance of FTOA,1 and FTOA,2 is 15 W m−2.
d. The multiple regression method
Suppose that we have a time series of surface temperatures and TOA radiative fluxes of the Earth, real or simulated, where the surface of Earth is regridded into ngrid (288) regions, and where we have ntime years of monthly observations. For each season s (1 ≤ s ≤ 4), we can define an ntime × ngrid matrix , where the element in row i and column j, Ti,j,s, is the surface temperature in region j during season s of year i. We can also define a matrix of anomalies, , where
To estimate the spatial feedbacks associated with a TOA radiative flux of type f (where f is one of net, LW clear, SW clear, LW cloud, or SW cloud) and season s, we first define an ntime × ngrid matrix of anomalies , which is analogous to above (N from the main body of the text is Rnet). We can fit the statistical model defined in Eq. (9) using least squares to solve for seasonal spatial feedbacks (Λf,s):
Seasonal feedbacks are used in section 3, but section 2 uses an annual version, in which case instead of a set of four seasonal feedback matrices, only one feedback matrix estimated using the above Eq. (A3), with the difference that the time series are annual averages. The “monthly” approach in section 1.2.1 of the SI is the same as the seasonal approach in Eq. (A3), except instead of 4 regressions, 12 are performed, with all time series being monthly averages sampled every 12 months. The “all months” approach instead performs only one regression, just like the annual approach, except that monthly average time series are used instead of annual averages (the logic being that even though months may have different properties, there may be an advantage in maximizing the data available to fit a regression).
e. Estimating the forced response
1) Forced feedbacks
Suppose that we have a ntime,abrupt4×-yr-long abrupt4× simulation of a GCM for which we have spatial feedbacks estimated from a control run. We then define an early period (years 2 to 20) and a late period (years 21 to ntime,abrupt4×). The true feedbacks λf,p for the abrupt4× simulation during each period p (where p is early or late) are defined as the slope of the least squares fit of the linear regression of the time series of globally averaged TOA flux anomalies of type f from the abrupt4× simulation , against the globally averaged surface temperature anomalies from the abrupt4x simulation :
where the curly brackets denote that the time series are averaged over exponentially longer periods, with annual averages for the first decade increasing to centennial averages by the simulation’s end, and the p subscript denotes whether values from before or after year 20 are used. Note that and are vectors with as many entries as years in the abrupt4× simulation (1000 yr).
We can make estimates of these feedbacks using the MR method by first estimating the abrupt4× simulation’s TOA radiative flux of type f for each month of the year m by multiplying the surface temperature time series of that abrupt4× simulation for that month, by the spatial feedbacks for that month’s season:
We use months instead of seasonal averages because our seasons do not start in January, and this approach allows us to have annual averages that start in January. These monthly time series can then be turned into annual averages , and then global averages , allowing us to estimate the feedbacks for period p by performing the same least squares fit as above:
2) Spatial patterns of TOA flux change
We quantify the normalized spatial pattern of TOA radiative flux change of flux type f across a period p by taking a finite difference approach, taking the mean value of during two parts of the period and subtracting the first part from the second (where the divisions for the early period are years 2–6 and 7–20, and the divisions for the late period are 21–170 and 171–ntime,abrupt4x, with both divisions chosen to allow for substantial warming in each period), and then dividing this by the average change in the globally averaged surface temperature between these two periods:
where tstart,p and tend,p are the first and last years in period p, respectively, where tmid,p is 6 for the early period and 170 for late period, where is the element in the ith row and jth column of , and where Tabrupt4×,i is the ith element in Tabrupt4×. Finite difference is used instead of regressing values against a global average because the presence of local and nonlocal feedbacks causes nonlinear relationships between and [or ], which would lead to biased estimates of change from a linear regression.
We calculate two types of errors: feedback errors (Table 1 and Table S2), and spatial errors (Table 2 and Table S3). We add a subscript g to our feedbacks and spatial patterns of TOA flux change to signify that they belong to the GCM g, where g is one of CCSM3, CESM1.0.4, GISS-E2-R, HadCM3L, IPSL-CM5A, and MPI-ESM1.2. The feedback error is given by the root-mean-square error:
where nGCMs is 6, the number of AOGCMs. The spatial error is measured by taking the area-weighted root-mean-square error of the spatial estimate
where ai is the area of the ith grid cell. For the spatial errors in the main body of the paper, this is taken on multimodel mean values of and . For the same calculation for individual models (Table S4 and Figs. S6–S8 in the supplemental materials), values for each model are used instead.
g. Other methods to calculate feedbacks
We consider two other methods for deriving spatial feedbacks, estimating abrupt4x feedbacks, and estimating spatial patterns of TOA flux change.
1) The global method
The seasonal version of the “global” method used in the main body of the paper is estimated using the least squares fit on this regression:
where Ts and Rf,s are globally and seasonally averaged time series of control simulation surface temperature and TOA flux f respectively, sampled every fourth seasonal value so that all elements of the time series are from season s. The four seasonal feedbacks are used to recreate estimates of the global averaged time series Rf,abrupt4×, which in turn is used, as above, to estimate abrupt4x feedbacks. Once more, different averaging of the control time series and groupings of regression equations can be used to make the annual, monthly, and all months versions of this method featured in Tables S3 and S4.
The normalized spatial pattern of TOA flux change can be found by first estimating the “local contribution” (Boer and Yu 2003a,b; Crook et al. 2011; Zelinka et al. 2012; Andrews et al. 2015), using Eq. (1), but replacing the time series vector with the spatial time series matrix from above, and replacing the single feedback λglobal,f,s with the spatial vector of feedbacks, λglobal,f.
2) The local method
The “local” method assumes the statistical model
Spatial feedbacks are estimated using least squares:
h. Local regression
We use LOESS (i.e., locally estimated scatterplot smoothing; Cleveland and Devlin 1988) to take local regression of scatterplots of versus . LOESS uses a weighted regression of a certain number of nearest neighbors, in our case 30. Full details can be found in the code for this paper listed above and in the LocallyWeightedRegression.jl Julia package (https://github.com/juliohm/LocallyWeightedRegression.jl).