Spatial Radiative Feedbacks from Internal Variability Using Multiple Regression (original) (raw)

1. Introduction

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 (T i) is proportional to the net energy gain of that region, which is the sum of the net downward top-of-atmosphere (TOA) radiative flux (N i), 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 (_F_surf,i):

c1dT1dt=N1−H+Fsurf,1,(1)

c2dT2dt=N2+H+Fsurf,2,(2)

where c i 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

Ti′,Ni′,H′,and Fsurf,i′

instead of T i, N i, H, and _F_surf,i. We assume that heat transport is proportional to the temperature gradient between the two regions:

H′=γ⁡(T1′−T2′).(3)

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

⁡(FCO2,i)

, and radiative forcing due to random atmospheric fluctuations that occur independently of surface temperature (_F_TOA,i):

N1′=λ1,1T1′+λ1,2T2′+FCO2,1+FTOA,1,(4)

N2′=λ2,1T1′+λ2,2T2′+FCO2,2+FTOA,2.(5)

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 N2′. The resulting horizontal advection may also increase the humidity of region 2’s free troposphere, increasing N2′. Assuming region 2 has a subsidence-induced boundary layer inversion, its low cloud cover could also increase, causing a further decrease in N2′. All of these changes in N2′ occur independently of any changes in T2′, and conspire to make _λ_2,1 positive or negative.

Fig. 1.

Fig. 1.

Fig. 1.

A schematic representation of the conceptual model used in section 2, consisting of an overturning cell with a convecting region (region 1) and a subsiding region (region 2). Warming of the surface temperature _T_1 has nonlocal effects: it increases the horizontal heat transport H, and it changes properties of the atmosphere aloft in region 2 that affect its net top-of-atmosphere radiative flux, _N_2, for instance by warming its free troposphere, increasing its lower tropospheric stability, and therefore increasing its low cloud cover. The dependence of _N_2 on _T_1 (holding _T_2 fixed) is an example of a nonlocal radiative feedback.

Citation: Journal of Climate 33, 10; 10.1175/JCLI-D-19-0396.1

We note that an increase in _H_′ will also increase T2′ 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).

Fig. 2.

Fig. 2.

Fig. 2.

Two experiments are performed with the conceptual model in Eq. (1): (a),(b) an unforced “control” simulation and (c),(d) a forced “abrupt4×” simulation. Values of N¯ vs T¯′ from each experiment are given by the black dots in (a) and (c), representing annual averages for the control simulation and exponentially increasing averages for the abrupt4× simulation. The global method assumes that the slope of the regression in (a) (blue line) gives the slope of the black dots in the lower left panel, underestimating the increase in this slope over time [blue lines and markers in (c),(d)]. The local method regresses Ni′ against Ti′ to estimate λ i for both regions [dotted lines in (b)], which leads to an overestimate of the combined feedback associated with region 1 [_λ_1 = _λ_1,1 + _λ_2,1, dotted red line in (b)], and therefore an overestimate of the feedback early on [orange lines and markers in (c),(d)]. The MR method, given sufficient years to regress over, correctly estimates all spatial feedbacks [dashed lines in (b)], accurately predicting the feedbacks and its change with time [green lines and markers in (c),(d)].

Citation: Journal of Climate 33, 10; 10.1175/JCLI-D-19-0396.1

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

λ⁡(t)=∂N¯∂T¯⁡(t)=∑ [ΛdTdt⁡(t)]dT¯dt⁡(t),(6)

where

T=[T1T2], Λ=[λ1,1λ2,1λ1,2λ2,2]

, and a bar over a quantity indicates the global average of that quantity. We do not have to use an anomaly for

N¯

because

N¯

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 Tcontrol′⁡(t) and Ncontrol′⁡(t) are due to random forcing at the surface {Fsurf′⁡(t)=[Fsurf,1′⁡(t)Fsurf,2′⁡(t)]} and TOA {FTOA⁡(t)=[FTOA,1′⁡(t)FTOA,2′⁡(t)]}, and an “abrupt4×” experiment in which the time series Tabrupt4×′⁡(t) and Nabrupt4×′⁡(t) also respond to an initial step forcing akin to a quadrupling of CO2 concentration ⁡(FCO2,1=FCO2,2=8 W m−2).

For the abrupt4× simulation, the climate feedback λ=∂N¯/∂T¯′ 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 N¯abrupt4×⁡(t) against T¯abrupt4×′⁡(t) 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 Tcontrol′⁡(t) and Ncontrol′⁡(t) (internal variability) and Tabrupt4×′⁡(t) (the spatial pattern of warming). The simplest method would be to regress annual averages of N¯control⁡(t) against T¯control⁡(t) 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 ⁡(FCO2) (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 Δ_λ_4_x = 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 λlocal=[λ1,localλ2,local] where λ i,local is the result of regressing Ni,control′⁡(t) against Ti,control′⁡(t). Taking the dot product of **λ**local with Tabrupt4×′⁡(t) then provides an estimate of Nabrupt4×′⁡(t) 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, T2′ tends to be positively correlated with T1′, due to heat transport, while T1′ tends to be anticorrelated with N2′ because _λ_2,1 is negative. As a result, the local method predicts that _λ_2,local is negative (dotted gray line, Fig. 2b), even though T2′ has no net influence on N. Since T2′ 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

Ni′

(i.e., the influence of

T1′

and

T2′

on

Ni′

) by regressing

Ni,control′⁡(t)

against both regions simultaneously:

Ni,control′⁡(t)=λi,1,MRT1,control′⁡(t)+λi,2,MRT2,control′⁡(t)+FTOA,i.(7)

In least squares multiple regression, λ i,j,MR is the same as the slope of the regression of

Ni,control′⁡(t)*

against

Tj,control′⁡(t)*

, 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

ΛMR=(λ1,1,MRλ1,2,MRλ2,1,MRλ2,2,MR)

by

Tabrupt,4×′⁡(t)

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, _F_surf,1 and _F_surf,2 are white noise with variance 20 W2 m−4, while _F_TOA,1 and _F_TOA,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:

Ni′⁡(t)=λi,1,MRT1,control′⁡(t)+λi,2,MRT2,control′⁡(t)+…+λi,n,MRTn,control′⁡(t)+FTOA,i,(8)

N′⁡(t)=ΛT′⁡(t)+FTOA,(9)

Fig. 3.

Fig. 3.

Fig. 3.

Plots of N¯ vs T¯′ for control simulations of six coupled atmosphere–ocean general circulation models (see Table S1 for details). We use the simulations to estimate spatial feedbacks using the global, local, and MR methods. We regrid simulations to 15° × 15° grids, giving 288 regions.

Citation: Journal of Climate 33, 10; 10.1175/JCLI-D-19-0396.1

Figures 3 and 4 show N¯ versus T¯′ of the control and abrupt4× simulations of the six models respectively. Figure 4 also shows N¯ estimated using the three methods, assuming that each estimate starts with the true value of N¯ at year 2. The solid lines in Fig. 4 are local regressions of N¯ against T¯′ 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 T¯ to estimate feedbacks as a function of time (lines in Fig. 5).

Fig. 4.

Fig. 4.

Fig. 4.

N¯ vs T¯′ for abrupt4× simulations of the same six GCMs from Fig. 3 (black dots). Colored dots show estimates of N¯abrupt4×⁡(t) made using the spatial feedbacks inferred from each model’s control simulation and its spatial pattern of warming [Tabrupt4×′⁡(t)] using the three methods described in the text; year one is not included in any method. Larger dots represent averages taken over exponentially increasing periods. Smaller gray dots show all years. Solid lines show local regressions using LOESS. Global estimates for GISS-E2-R do not appear because it is nearly identical with MR estimates.

Citation: Journal of Climate 33, 10; 10.1175/JCLI-D-19-0396.1

Fig. 5.

Fig. 5.

Fig. 5.

True and estimated abrupt4× feedbacks as a function of time calculated using slopes of the local regression from Fig. 4 (solid lines). Vertical dotted lines show the division between the early (2–20 yr) and late (21 yr–end) periods. Dots show true and estimated values of _λ_4×,early and _λ_4×,late. Feedbacks get more positive over time for all models. The MR and global methods initially overestimate feedbacks. The MR estimate increases with time as well, while the global method predicts a roughly constant feedback. The local method greatly overestimates the true feedback for all models except GISS-E2-R. Figures S2–S5 give the same plot for component fluxes.

Citation: Journal of Climate 33, 10; 10.1175/JCLI-D-19-0396.1

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).

Fig. 6.

Fig. 6.

Fig. 6.

True vs estimated feedbacks for the (a)–(c) early and (d)–(f) late periods, and (g)–(i) the change between them. Black dots give values for the net feedback, while colored markers give values of the component feedbacks, which sum to the net feedback. The MR and global methods overestimate the early feedback due to SW cloud (red) feedbacks. The MR estimate of the late period in (d) has a small error across all components, while the global estimate in (e) has a smaller net error due to offsetting errors between LW and SW cloud feedbacks. As in Fig. 5, the MR method is able to capture some of the change in feedback, while the global method does not. The local method greatly overestimates the net feedback, primarily due to cloud feedbacks. Numerical values of the feedbacks are given in Table S7 and S8.

Citation: Journal of Climate 33, 10; 10.1175/JCLI-D-19-0396.1

Table 1.

Feedback errors. Root-mean-square errors of estimates of abrupt4× feedbacks (_λ_4×,early, _λ_4×,late) and their change with time (Δ_λ_4×), for net TOA fluxes and each component flux (in W m−2 K−1) and for the seasonal versions of the three methods presented in section 2 (see the appendix for details). For annual and monthly values, see Table S2, and for fluxes north of 30°S, see Table S5.

Table 1.

Table 1.

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 N4×′ (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 N4×′⁡(t) 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 Ncontrol′⁡(t) against T¯control′⁡(t) [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 N4×′.

Fig. 7.

Fig. 7.

Fig. 7.

Multimodel mean spatial pattern of net TOA flux change associated with the (top) early and (middle) late periods, and (bottom) the change between them, calculated by taking the finite difference across each period. Changes are normalized by the total warming in each period, giving units of W m−2 K−1. The MR method is close to the true pattern except for overestimates south of 30°S and during the early period in the North Atlantic. This holds for individual flux components as well (Figs. S9–S17). The global and local methods both have substantial errors over most of the globe. Figures S6–S8 show errors (estimates − true values) for the multimodel mean and individual models.

Citation: Journal of Climate 33, 10; 10.1175/JCLI-D-19-0396.1

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).

Table 2.

Spatial errors. The model-mean area-weighted root-mean-square error of estimates of the warming-normalized change in TOA fluxes during the early and late periods of the abrupt4× simulations, and the change in pattern between these period (see the appendix for details). All values have units of W m−2 K−1. For annual and monthly versions in addition to seasonal, see Table S2, for individual models see Table S4, and for fluxes north of 30°S, see Table S6.

Table 2.

Table 2.

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.

4. Discussion

We first test if the spatial feedbacks estimated using the MR method exhibit behavior broadly consistent with physically modeled feedbacks. The _i_th 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).

Fig. 8.

Fig. 8.

Fig. 8.

Net cloud feedbacks associated with warming in regions circled in green estimated for CAM5 (a)–(c) by Zhou et al. (2017) using fixed-SST experiments or (d)–(f) as a multimodel and multiseason mean using the MR method. For perturbations outside of tropical convecting regions [in (a), (c), (d), and (f)], the effects are mostly local and positive, while perturbations in tropical convecting regions have significant negative nonlocal effects in many regions of Earth [in (b) and (e)]. Note that the fixed-SST experiments in (b) allow some land warming in response to these perturbations, while the MR method in (e) is agnostic about whether the surface is land or ocean, and so does not include resulting land warming.

Citation: Journal of Climate 33, 10; 10.1175/JCLI-D-19-0396.1

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 N¯ 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).

Fig. 9.

Fig. 9.

Fig. 9.

Multimodel and multiseason mean spatial feedbacks estimated by the MR method. (a) The estimated change in N¯ caused by warming a degree in each cell as weighted by the cell’s area. This is the sum of (b) local changes in N, which are almost uniformly positive, and (c) nonlocal changes, which are usually negative, especially in regions of tropical convection. (d)–(f) The competing positive local and negative nonlocal components are primarily due to the SW cloud feedback. For maps of all flux components and assessments of uncertainty, see Figs. S18–S22. For spatial feedbacks of all methods, see Fig. S23. Compare with estimates of spatial feedbacks for CAM4 in Fig. 5c of Dong et al. (2019).

Citation: Journal of Climate 33, 10; 10.1175/JCLI-D-19-0396.1

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).

Fig. 10.

Fig. 10.

Fig. 10.

(a) The multimodel mean change in the pattern of warming between the abrupt4× early and late period, showing a shift toward regions of deep ocean heat uptake. (b)–(f) Multiplying this pattern by MR-estimated spatial feedbacks gives an estimate of each grid cell’s contribution to the change in feedback with time, Δ_λ_4×. Although the resulting patterns are patchy, there are positive contributions from tropical convecting regions via the SW cloud and LW clear feedbacks, and from regions of Southern Ocean sea ice in the SW clear feedback, as shown by the accompanying zonal averages. The LW clear feedback has a compensating negative term from the Southern Ocean, so that its total estimated contribution to Δ_λ_4× is smaller than the SW cloud feedback’s (e.g., Fig. S2 vs Fig S5).

Citation: Journal of Climate 33, 10; 10.1175/JCLI-D-19-0396.1

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 Δ_λ_4_x_.

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 N¯ 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.

5. Conclusions

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.

Acknowledgments

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.

APPENDIX

Data and Methods

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:

c1dT1′dt=N1′−H′+Fsurf,1,

c2dT2′dt=N2′+H′+Fsurf,2,

where

H′=γ⁡(T1′−T2′)

and

N1′=λ1,1T1′+λ1,2T2′+FCO2,1+FTOA,1,(A1)

N2′=λ2,1T1′+λ2,2T2′+FCO2,2+FTOA,2.(A2)

The thermal inertia c i is defined as m i ρc p, where ρ and c p are the density and specific heat of ocean water respectively, and m i is an equivalent mixed layer depth; _m_1 is 50 m, and _m_2 is 1000 m.

FCO2,1=FCO2,2

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 _F_surf,1 and _F_surf,2 is 40 W m−2 and the variance of _F_TOA,1 and _F_TOA,2 is 5 W m−2, while for the example in Fig. S1, the variance of _F_surf,1 and _F_surf,2 is 10 W m−2 and the variance of _F_TOA,1 and _F_TOA,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 _n_grid (288) regions, and where we have _n_time years of monthly observations. For each season s (1 ≤ s ≤ 4), we can define an _n_time × _n_grid matrix

Tm

, where the element in row i and column j, T i,j,s, is the surface temperature in region j during season s of year i. We can also define a matrix of anomalies,

Ts′

, where

Ts′=[T1,1,sT1,2,s…T1,ngrid,sT2,1,sT2,2,s…T2,ngrid,s⋮⋮⋱⋮Tntime,1,sTntime,2,s…Tntime,ngrid,s] −1ntime[∑i=1ntimeTi,1,s∑i=1ntimeTi,2,s…∑i=1ntimeTi,ngrid,s∑i=1ntimeTi,1,s∑i=1ntimeTi,2,s…∑i=1ntimeTi,ngrid,s⋮⋮⋱⋮∑i=1ntimeTi,1,s∑i=1ntimeTi,2,s…∑i=1ntimeTi,ngrid,s]

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 _n_time × _n_grid matrix of anomalies

Rf,s′

, which is analogous to

Ts′

above (N from the main body of the text is R_net). We can fit the statistical model defined in Eq. (9) using least squares to solve for seasonal spatial feedbacks (Λ_f,s):

Λf,s=[λf,1,1λf,1,2…λf,1,ngridλf,2,1λf,2,2…λf,2,ngrid⋮⋮⋱⋮λf,ngrid,1λf,ngrid,2…λf,ngrid,ngrid]=⁡(Ts′TTs′)−1Ts′TRf,s′.(A3)

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 _n_time,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 _n_time,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

⁡(Rf,abrupt4x′)

, against the globally averaged surface temperature anomalies from the abrupt4x simulation

Tabrupt4×′

:

λabrupt4×,f,p={Tabrupt4×′}p⋅{Rf,abrupt4×′}p‖{Tabrupt4×′}p‖2,(A4)

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

Rf,abrupt4×′

and

Tabrupt4×′

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,

Tm,abrupt4×′⁡(an ntime,abrupt4××ngrid matrix)

by the spatial feedbacks for that month’s season:

R^f,m,abrupt4×′=Tm,abrupt4×′Λf,s⁡(m)(A5)

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

R^f,m,abrupt4×′

can then be turned into annual averages

R^f,abrupt4×′

, and then global averages

R^f,abrupt4×′

, allowing us to estimate the feedbacks for period p by performing the same least squares fit as above:

λ^abrupt4×,f,p={Tabrupt4×,p′}⋅{R^f,abrupt4×,p′}‖{Tabrupt4×,p′}‖2.(A6)

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

Rf,abrupt4×′

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–_n_time,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:

ΔRf,abrupt4×,p′=[∑i=tmid,p+1tend,p⁡(Rf,abrupt4×,i,1′Rf,abrupt4×,i,2′⋮Rf,abrupt4×,ngrid′)−∑i=tstart,ptmid,p⁡(Rf,abrupt4×,i,1′Rf,abrupt4×,i,2′⋮Rf,abrupt4×,ngrid′)]⁡(∑i=tmid,p+1tend,pTabrupt4×,i−∑i=tstart,ptmid,pTabrupt4×,i),(A7)

where _t_start,p and _t_end,p are the first and last years in period p, respectively, where _t_mid,p is 6 for the early period and 170 for late period, where

Rf,abrupt4×,i,j′

is the element in the _i_th row and _j_th column of

Rf,abrupt4×′

, and where _T_abrupt4×,i is the _i_th 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

Ni′⁡(t)

and

Ti′⁡(t)

[or

T¯′⁡(t)

], which would lead to biased estimates of change from a linear regression.

f. Errors

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:

εfeedback,f,p=1nGCMs∑g∈GCMs⁡(λ^f,abrupt4×,p,g−λf,abrupt4×,p,g)2,(A8)

where _n_GCMs is 6, the number of AOGCMs. The spatial error is measured by taking the area-weighted root-mean-square error of the spatial estimate

εspatial,f,p=∑i=1ngrid⁡(ΔR˜^f,abrupt4×,p,i′−ΔR˜f,abrupt4×,p,i′)2ai∑i=1ngridai,(A9)

where a i is the area of the _i_th grid cell. For the spatial errors in the main body of the paper, this is taken on multimodel mean values of

ΔR˜^f,abrupt4×,p,i′

and

ΔR˜f,abrupt4×,p,i′

. 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:

λglobal,f,s=Ts′⋅Rf,s′‖Ts′‖2,(A10)

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 Rf,s′ with the spatial time series matrix Rf,s′ 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

Ri′⁡(t)=λlocal,iTi′⁡(t)+ε⁡(t) for each region i.(A11)

Spatial feedbacks are estimated using least squares:

λlocal,f=[λlocal,f,1λlocal,f,2⋮λlocal,f,ngrid]=[T1′⋅Rf,1′‖T1′‖2T2′⋅Rf,2′‖T2′‖2⋮Tngrid′⋅Rf,ngrid′||Tngrid′||2],(A12)

where

Ti′

and

Rf,i′

are the _i_th rows of

T′

and

Rf′

respectively. We can then generate estimates of

Rf,abrupt4×′

as above. We apply these estimates to Eqs. (A6) and (A7) to estimate forced global feedbacks and spatial patterns of TOA flux change.

h. Local regression

We use LOESS (i.e., locally estimated scatterplot smoothing; Cleveland and Devlin 1988) to take local regression of scatterplots of N¯ versus T¯′. 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).

REFERENCES