1 Introduction

The need to understand factors influencing climate feedback and sensitivity make simple models indispensable as means to investigate the climate response to changes in the Earth’s energy balance. In this paper we discuss the relationship between the energy balance of the climate system and its temperature response. We focus on the temporal behavior using the conceptual framework of four commonly used conceptual models. We reveal similarity and essential differences among these models, and we provide insight into the climate response in the Max Planck Institute Earth System Model 1.2 (MPI-ESM1.2) as a counterpart to our theoretical considerations. On mathematical grounds, we question the difference between the concept of ocean heat uptake efficacy and a pattern effect due to regional feedbacks weighted by a geographic pattern of surface warming. Further, we challenge the common assumption of these concepts that the characteristic timescales on which the climate system adjusts are independent of the strength of the external perturbation.

A popular framework to analyze the radiative response to external radiative forcing is the linear forcing-feedback framework (e.g. Gregory et al. 2004). It is based on a linear relationship between global net top-of-the-atmosphere (TOA) radiative flux N, radiative forcing F and surface temperature perturbation T which in turn actuates the radiative feedback parameter \(\lambda\);

$$\begin{aligned} N=F+\lambda T. \end{aligned}$$
(1)

The linear forcing-feedback framework Eq. (1) is an appropriate first-order description for the relationship between N and T in experiments of state-of-the-art atmosphere-ocean coupled general circulation models (AOGCMs) in which the preindustrial atmospheric CO\(_2\) concentration is abruptly doubled and then held constant. Details, such as the possibility of a nonlinear dependence on T are omitted from this model by construction. Numerical experiments with more comprehensive models suggest that these details may matter, and that Eq. (1) becomes an increasingly poor description as the forcing is increased to 4 or 8 times preindustrial CO\(_2\). In these cases, there is a tendency for the radiative response in several coupled climate models to decrease over time and with higher temperatures (e.g. Andrews et al. 2015; Bloch-Johnson et al. 2015). If one assumes that Eq. (1) holds at each instance of time, this can be interpreted as nonconstant global feedback. This is modeled by allowing \(\lambda\) to vary over time and rearrange Eq. (1) to obtain

$$\begin{aligned} \lambda (t)=\frac{N(t)-F(t)}{T(t)}, \end{aligned}$$
(2)

which is sometimes paraphrased as effective feedback (Senior and Mitchell 2000; Williams et al. 2008; Armour et al. 2013). To avoid ambiguities in terminology, we refer to \(\lambda (t)\) as the radiative response of the climate system relative to the global mean surface temperature perturbation T. The outcome of Eq. (2) may, however, also influenced by additional state-variables. It is itself a specific approach to analyze the climate response to external perturbations and provides a simple and vivid description of the radiative response. Unfortunately, to the extent that \(\lambda\) changes with t it does not alone provide a framework for understanding the origin of such changes which becomes evident during the course of this study.

Currently, there are two different branches of literature on conceptual models to explain the more complex climate response in AOGCM experiments. The first branch suggests that on decadal timescales the variation of \(\lambda (t)\) arises from different radiative responses associated with a fast and a slow component of the climate system (two-timescale approaches). The second branch accounts for the possibility that the strength of radiative feedbacks depends explicitly on surface temperature change (feedback temperature dependence).

Hasselmann (1979) already recognized the multiple timescale structure of the climate response, but at this time the global feedback was assumed to depend on the global mean surface temperature only. Recently, these characteristic adjustment timescales have been explicitly linked to nonconstant global feedback associating them with different state-variables each of which influences the radiative response (e.g. Held et al. 2010; Winton et al. 2010; Geoffroy et al. 2013a, b; Armour et al. 2013; Andrews et al. 2015; Hedemann et al. 2019). Held et al. (2010) and Geoffroy et al. (2013a, b) describe a model based on two distinct timescales to reproduce the transient global mean temperature response of a wide range of AOGCMs. They introduce a globally averaged two-layer ocean model with an efficacy parameter modifying the upper ocean response to heat uptake. Hereafter it is referred to as a two-layer model with efficacy, or just the two-layer model. Two or more adjustment timescales can also be explicitly linked to a one-layer model with two regions. In contrast to Held et al. (2010) and Geoffroy et al. (2013a), Armour et al. (2013) suggest a geographic approach in which constant regional feedbacks are weighted by an evolving pattern of surface warming. This approach is motivated by the idea that the adjustment over parts of the Earth’s surface which are weakly coupled to the state of the deep-ocean is more rapid than over regions which are more strongly coupled to the state of the deep-ocean. Consequently, this can be called as pattern effect (Stevens et al. 2016). In addition to time-varying sensitivity, two-region models have also been used to analyze the effect of perturbation heat transport on stability properties and equilibrium sensitivity (e.g. Bates 2012, 2016) or to explore polar amplification (e.g. Langen and Alexeev 2007). In all of these models, the time-variation of \(\lambda (t)\) does not depend on the strength of the forcing– in this sense these models are linear in forcing. The tendency of such models to have a value of \(\lambda (t)\) that varies with time is thus often paraphrased as a time-dependent feedback though it arises from a dependence of the radiative response on an additional state-variable with its own distinct temporal evolution. Using AOGCM output, e.g. Senior and Mitchell (2000) and Hedemann et al. (2019) discuss time-dependent feedback as a way to understand the radiative response in more comprehensive climate models.

A different approach to representing variation in the radiative response has been to relax the assumption that a temperature perturbation would be proportional to the associated radiative forcing F (e.g. Roe and Baker 2007; Zaliapin and Ghil 2010; Roe and Armour 2011; Meraner et al. 2013; Bloch-Johnson et al. 2015). Studies adopting this approach investigate the long-term response on centennial to millennial timescales as well as the nonlinear change in equilibrium sensitivity between different forcing strengths. Conceptual models on nonlinear feedbacks of the climate system can be traced back at least to Budyko (1969) and Sellers (1969). In a recent study, Bloch-Johnson et al. (2015) extend the linear forcing-feedback framework Eq. (1) by a quadratic coefficient a representing second-order feedback temperature dependence and solve for stationary solutions, \(-F=\lambda T + a T^2\). They estimate the range of likely feedback temperature dependencies for various GCM simulations and find \(-0.04 \le a \le 0.06\hbox { W m}^{-2}\hbox { K}^{-2}\). This range is in line with Roe and Armour (2011) who find \(-0.058 \le a \le 0.06\hbox { W m}^{-2}\hbox { K}^{-2}\). Positive feedback temperature dependence makes the feedback of the climate system continuously less effective as the surface temperature increases, while the opposite is true for negative feedback temperature dependence. In this connection, the radiative response depends explicitly on the climate state; i.e., feedback temperature dependence gives rise to state-dependent feedback such that the variation of \(\lambda (t)\) is changed as the strength of the forcing is altered.

We formulate the different conceptual models discussed above as simple ordinary differential equations and study their dynamics. To provide a basis, in Sect. 2 we describe the climate response found in an AOGCM. Then, in Sect. 3, we discuss dynamical systems of the two-timescale approaches and demonstrate that the two-layer model with efficacy can be projected onto a two-region model. In Sect. 4 we analyze dynamical systems which incorporate feedback temperature dependence to point out that the temporal behavior of the climate system, or more precisely the temporal behavior of the comprehensive model developed at our institute, is indeed forcing-dependent.

2 Noncostant global feedback in MPI-ESM1.2

To understand the radiative response of the climate system, state-of-the-art climate models such as the Max Planck Institute Earth System Model 1.2 (MPI-ESM1.2), thoroughly described in Mauritsen et al. (2018), can be forced by step forcing—usually an abrupt increase in atmospheric CO\(_2\) concentration. Such experiments reveal changes or common features in the relationship between N and T (e.g. Gregory et al. 2004). We briefly illustrate the MPI-ESM1.2 response to an abrupt increase in atmospheric CO\(_2\). A detailed analysis of radiative feedbacks in MPI-ESM1.2 or in its atmospheric component ECHAM6.3 is provided by Mauritsen et al. (2018) and Hedemann et al. (2019). MPI-ESM1.2 is forced by abrupt 2, 4, 8 and 16 times preindustrial CO\(_2\) (284.7 ppm) and is integrated forward in time to 1000 years. In this case, the radiative response is related to radiative feedbacks which emerge on decadal and centennial timescales such as the water vapor, lapse rate, clouds and snow/sea ice albedo feedbacks.

Fig. 1
figure 1

The relationship between the TOA imbalance N and surface temperature change T in MPI-ESM1.2 abrupt CO\(_2\) experiments using annual means. The black line represents a linear least-square regression from year 1–25. The dashed line is the extrapolation of this regression to \(N=0\)

Full size image

Figure 1 shows the MPI-ESM1.2 relationship between global and annual mean N and T. A common feature of these experiments is that there is an abrupt change in this relationship on decadal timescales. The model output for each CO\(_2\) doubling deviates significantly from the linear least-square regression (solid black line) which is computed using only output from year 1 to 25. Extrapolating the initial response (dashed lines), as would be appropriate if the radiative response to temperature were constant, leads to an underestimation of the equilibrium temperature response, increasingly so for larger forcing. In this connection, the initial variation of \(\lambda (t)\) can be explained by distinct timescales on which the climate system adjusts. This is at the heart of the theory of the two-timescale approaches (Sect. 3).

At the same time, the change in the long-term temperature response between the consecutive doublings substantially increases with warmer climates. State-of-the-art models reveal a nearly logarithmic dependence of the radiative forcing F on atmospheric CO\(_2\) concentration as predicted by theory which can be traced back to Arrhenius (1896). That is, the change in the associated radiative forcing F imposed by each doubling cannot explain the nonlinear rise in the equilibrium temperature response: extrapolating the relationship between N and T in MPI-ESM1.2 linearly from year 100 to 1000 gives a change in equilibrium sensitivity from \(2\times \text {CO}_2\) to \(4\times \text {CO}_2\) of about 3.6 K and from the \(4\times \text {CO}_2\) experiment to the \(8\times \text {CO}_2\) of about 4.9 K, while the change from the \(8\times \text {CO}_2\) experiment to the \(16\times \text {CO}_2\) is roughly 9.8 K. As discussed in the introduction, this type of nonlinearity has been attributed to feedback temperature dependence (e.g. Meraner et al. 2013; Bloch-Johnson et al. 2015).

Fig. 2
figure 2

The MPI-ESM1.2 temporal behavior in abrupt CO\(_2\) experiments using annual means. a The fraction of the equilibrium response realized at a given time. For a linear response the lines should lie atop one another. b The time difference of specific response times (\(\tau _{\mathrm {80}}\), \(\tau _{\mathrm {75}}\), \(\tau _{\mathrm {70}}\)) to the e-folding time \(\tau _{\mathrm {63}}\). Dashed lines are quadratic fits to the data points

Full size image

After 1000 years the \(8\times \text {CO}_2\) and \(16\times \text {CO}_2\) experiment are still far from steady state, whereas the \(2\times \text {CO}_2\) experiment is nearly equilibrated. After a forcing is applied, the temperature relaxes to a new equilibrium at specific timescales and the temperature increase actuates radiative feedbacks that amplify or attenuate surface warming. That the adjustment time-scale depends on the forcing is demonstrated by comparing at what time a given fraction of the equilibrium warming is realized as a function of forcing (Fig. 2a). The e-folding time \(\tau _{\mathrm {63}}\) and \(\tau _{\mathrm {80}}\), defined as the time when a fraction of 63 and 80 percent of the equilibrium response is reached, are compared. Evidently, the response times change as the radiative forcing is altered, and these changes occur in a nonlinear way. The change in response times acts on all timescales while the absolute change between different CO\(_2\) doublings increases with higher forcing and higher fractions of the equilibrium response. To highlight this behavior, we plot the difference between specific response times (\(\tau _{\mathrm {80}}\), \(\tau _{\mathrm {75}}\), \(\tau _{\mathrm {70}}\)) to the e-folding time \(\tau _{\mathrm {63}}\) (Fig. 2b). To compute specific response times, the MPI-ESM1.2 temperature output has been smoothed by running mean with a window of thirty years. The difference between these response times increases approximately quadratically as atmospheric CO\(_2\) is progressively doubled. Section 4 explains this phenomenon.

3 Two-timescale approaches

In the linear forcing-feedback framework Eq. (1) the imbalance is caused by the radiative forcing F which is equivalent to a radiative perturbation imposed by a radiative forcing agent such as CO\(_2\). In the following we explore the temporal dynamics of the energy budget and temperature response of the conceptual models subject to a step (Heaviside) radiative forcing, F. This is conceptually equivalent to the forcing imposed by an abrupt increase in atmospheric CO\(_2\) in a more comprehensive model. The time-dependent unit-area perturbation equation for a zero-dimensional model like Eq. (1) can be generalized by

$$\begin{aligned} \frac{\mathrm {d} E}{\mathrm {d} t}= N \end{aligned}$$
(3)

where N is the downward energy flux perturbation per unit area at TOA, and E is the system’s enthalpy per unit area. We assume that the atmosphere is always in energetic balance because its heat capacity is much smaller than the oceanic heat capacity. This way, we can parameterize the change in the climate system’s enthalpy by the surface temperature perturbation T multiplied by the effective oceanic heat capacity C such that

$$\begin{aligned} \frac{\mathrm {d} E}{\mathrm {d} t} \equiv C\frac{\mathrm {d} T}{\mathrm {d} t}=N. \end{aligned}$$
(4)

We obtain an ordinary differential equation in which surface temperature is the state-variable of an initial value problem. Different models in the literature propose different parametrizations for N. The ocean heat uptake efficacy model is a globally averaged two-layer model (Held et al. 2010; Geoffroy et al. 2013a, b) and based on two state-variables which represent an upper ocean- or mixed-layer temperature (T, comprising atmosphere, land and ocean) and a deep-ocean temperature (\(T_\mathrm {D}\)). As a counterpart, we consider a configuration of the pattern effect which comprises two distinct regions. We demonstrate that, on mathematical grounds, these models can describe the same global temperature evolution and radiative response.

3.1 Two-layer model

The two-layer model with ocean heat uptake efficacy is given by

$$\begin{aligned}&{C}\frac{\mathrm {d} T}{\mathrm {d} t}= F + \lambda _{\mathrm {eq}} T - \varepsilon \eta (T-T_{\mathrm {D}}) \end{aligned}$$
(5)
$$\begin{aligned}&{C_{\mathrm {D}}} \frac{\mathrm { d} T_{\mathrm {D}}}{\mathrm {d} t}= \eta (T-T_{\mathrm {D}}) \end{aligned}$$
(6)

where \(C \ll C_{\mathrm {D}}\) are the heat capacities of the upper- and deep-ocean respectively, \(\lambda _{\mathrm {eq}}\) is the equilibrium feedback parameter, \(\eta\) is the heat transport efficiency and \(\varepsilon\) the efficacy factor for ocean heat uptake. The concept of ocean heat uptake efficacy has been originally introduced by Winton et al. (2010) in analogy to the efficacy of radiative forcing agents (e.g. Hansen 2005). Following Stevens et al. (2016) we prefer to think instead of the radiative response depending on the disequilibrium of the deep-ocean because a certain surface temperature pattern may accompany this disequilibrium. This can be more easily seen by writing the TOA imbalanbce \(N= C \frac{\mathrm {d} T}{\mathrm {d} t} + {C_\mathrm {D}} \frac{\mathrm {d} T_\mathrm {D}}{\mathrm {d} t}\) as

$$\begin{aligned} N=F + \lambda _{\mathrm {eq}} T - (\varepsilon -1) \eta (T-T_\mathrm {D}). \end{aligned}$$
(7)

Hence, \((\varepsilon -1) \eta (T-T_{\mathrm {D}})\) is an additional radiative flux which is lost to space, while the heat transport into the deep-ocean corresponds to \(\eta (T-T_{\mathrm {D}})\).

To begin with, we briefly characterize basic features of the two-layer model. Subsequently, we project the two-layer model onto a two-region model. The initial change of N with respect to T is

$$\begin{aligned} \frac{\mathrm {d}N}{\mathrm {d}T} \vert _{t=0}= \lambda _{\mathrm {eq}} - (\varepsilon -1) \eta , \end{aligned}$$
(8)

while it approaches \(\lambda _{\mathrm {eq}}\) as the mixed-layer and deep-ocean temperatures converge. Thus, equilibrium sensitivity is given by \(-F/\lambda _{\mathrm {eq}}\) and independent of the coupling between the mixed-layer and deep-ocean, whereas the coupling determines the transient response. More precisely, nonunitary efficacy causes a temporal change in the radiative response that attenuates warming; i.e., \(\varepsilon>1\) leads to a decrease in the magnitude of the radiative response over time as the deep ocean heat uptake changes. Applying Eq. (2), the evolution of \(\lambda (t)\) is solely determined by the ratio between T and \(T_\mathrm {D}\) and it is a monotonically increasing function of time starting from a strong negative feedback. Strictly speaking, the time-variation of \(\lambda (t)\), or time-dependent feedback, does not depend explicitly on time but implictly because it emerges from an additional state-variable which also adjusts on forcing-invariant timescales. In the case of the two-layer ocean model, \(\lambda (t)\) increases or decreases gradually over time because it is relative to the surface temperature perturbation T. This gradual increase or decrease, however, leads to an apparant abrupt change in the relationship between N and T, and the strength of this change is determined by the efficacy factor. The abrupt change in MPI-ESM1.2 on decadal timescales implies an efficacy factor above unity.

Table 1 Parameters of the analytical solution for the two-layer model
Full size table

The analytical solution for each layer, thoroughly described in Geoffroy et al. (2013a, b) but presented here independently to highlight important features, consists of the sum of the equilibrium response and two modes which increase exponentially towards zero with their characteristic timescales. These modes are characterized by the fast and slow adjustment timescale \(\tau _\mathrm {f}\) and \(\tau _\mathrm {s}\), their amplitudes \(\psi _\mathrm {f}\) and \(\psi _\mathrm {s}\) which sum up to the initial ocean heat uptake temperature, and the mode parameters \(\zeta _\mathrm {f}\) and \(\zeta _\mathrm {s}\) which scale the adjustment of the deep-ocean layer (parameters are listed in Table 1);

$$\begin{aligned} T(t)= \frac{F}{-\lambda _{\mathrm {eq}}}+ \psi _\mathrm {f} e^{-t/\tau _\mathrm {f}} + \psi _\mathrm {s} e^{-t/\tau _\mathrm {s}} \end{aligned}$$
(9)
$$\begin{aligned} T_\mathrm {D}(t)= \frac{F}{-\lambda _{\mathrm {eq}}} + \zeta _\mathrm {f} \psi _\mathrm {f} e^{-t/\tau _\mathrm {f}} + \zeta _\mathrm {s} \psi _\mathrm {s} e^{-t/\tau _\mathrm {s}}. \end{aligned}$$
(10)

The amplitudes \(\psi _\mathrm {f}\), \(\psi _\mathrm {s}\) and \(\zeta _\mathrm {s} \psi _\mathrm {s}\) correspond to the negative heat uptake temperatures in the mixed-layer and the deep-ocean, and \(\zeta _\mathrm {f} \psi _\mathrm {f}\) is the positive contribution of the mixed-layer to the deep-ocean adjustment. The fast timescale \(\tau _\mathrm {f}\) is related to the small heat capacity of the mixed-layer and decreases as the coupling \(\varepsilon \eta\) or the strength of the feedback parameters \(\lambda _{\mathrm {eq}}\) increases. The slow timescale \(\tau _\mathrm {s}\) is related to the high heat capacity of the deep-ocean and likewise decreases as the strength of \(\eta\) or \(\lambda _{\mathrm {eq}}\) increases. However, the opposite is true for \(\varepsilon\), because the efficacy factor causes a temporary enhancement of the radiative response such that the heat exchange between the layers is prolonged.

3.2 Two-region model

Regional feedbacks can also be formulated as linearization about a regional energy budget. A common approach is a first-order linearization which gives

$$\begin{aligned} N(r,t)=F(r,t)+\lambda (r) T(r,t) - Q(r) \end{aligned}$$
(11)

where r denotes (not necessarily contiguous) regions and Q(r) denotes heat transport between these regions (e.g. Armour et al. 2013; Rose and Rayborn 2016). Even though regional feedbacks \(\lambda (r)\) are taken as constant this system may introduce nonconstant global feedback, because regional feedbacks are weighted by a pattern of surface warming. In the following we assume a two-region model in which the climate response in each region is characterized by a specific heat capacity C(r) and a feedback parameter \(\lambda (r)\). In analogy to the two-layer or two-region model we assume Heaviside step forcing input in both regions and five model parameters. The global mean TOA imbalance is given by

$$\begin{aligned} {\overline{N}}= F + \chi \lambda _\mathrm {1} T_\mathrm {1}+ (1-\chi ) \lambda _\mathrm {2} T_\mathrm {2} \end{aligned}$$
(12)

where the bar denotes the spatial mean and \(\chi\) and \((1-\chi )\) are the corresponding area-fractions. The evolution of \(\lambda (t)\) results from the weighting of the regional feedback parameters by the regional temperatures normalized by the global mean temperature perturbation (Armour et al. 2013). Further, we assume that the two regions are not coupled. Regions such as latitudinal bands can hardly be taken as decoupled, since oceanic and atmospheric energy transport alter the regional energy budget significantly. That is, dynamical aspects of climate sensitivity become important and the regional radiative response is related to remote forcing. Adding a parametrization for \(-Q(r)\) which scales linearly with the temperature gradient, the characteristic timescales \(\tau _\mathrm {1}\) and \(\tau _\mathrm {2}\) are relatively shortened while the nonlinearity between \({\overline{N}}\) and \({\overline{T}}\) is partially compensated for (“Appendix”). Langen and Alexeev (2007) and Bates (2012, 2016) focus on implications of perturbation heat transport for the temperature response in tropical-extratropical two-region models. Here, instead, we assume that an AOGCM can be aggregated on two different regions. The analytical solution for regional temperature change in the two-region model is characterized by the sum of the equilibrium response \(-F/\lambda (r)\) and a single mode which converges exponentially towards zero on its own (regional) timescale;

$$\begin{aligned} T(t,r)=\frac{-F}{\lambda (r)} (1-e^{-t/\tau (r)}). \end{aligned}$$
(13)

The previously discussed two-layer model with ocean heat uptake efficacy is a coupled model. However, it is a linear system, and this allows us to find analytical relationships between it and the uncoupled two-region model. These relationships are based on the condition that the globally averaged surface temperature T and TOA imbalance N of the two-layer model are equal to the global mean temperature perturbation \({\overline{T}}\) and the global mean TOA imbalance \({\overline{N}}\) of the two-region model. We assume that the radiative response of the climate system can be aggregated on these different regions and link the parameters for (regional) feedbacks to the fast component of the climate response and to the slow component of the climate response. The feedback parameters are then given by

$$\begin{aligned}&\lambda _\mathrm {1}=\lambda _{\mathrm {eq}} + (\varepsilon \eta - \eta ) (\zeta _\mathrm {f}-1) \end{aligned}$$
(14)
$$\begin{aligned}&\lambda _\mathrm {2}=\lambda _{\mathrm {eq}} + (\varepsilon \eta - \eta ) (\zeta _\mathrm {s}-1) \end{aligned}$$
(15)

where difference in the radiative response coefficients is set by \(\zeta _\mathrm {f}\) and \(\zeta _\mathrm {s}\). The term \(\zeta _\mathrm {f}\) is negative and scales the contribution of the mixed-layer to the deep-ocean adjustment. The term \(\zeta _\mathrm {s}\) is positive and scales the slow mode or heat uptake temperature of the deep-ocean component itself. Thus, the regional feedback parameters can be decomposed into a background value \(\lambda _{\mathrm {eq}} -\varepsilon \eta + \eta\) and a contribution from the fast or slow mode. The regional timescales must correspond to \(\tau _\mathrm {f}\) and \(\tau _\mathrm {s}\) and the regional heat capacity is thus given by

$$\begin{aligned} C(r)=-\lambda (r) \tau (r). \end{aligned}$$
(16)

Finally, the regional sensitivity or radiative response is weighted by the corresponding area-fraction \(\chi\) and \((1-\chi )\) to result in \({\overline{T}}\) and \({\overline{N}}\). This weighting is done by

$$\begin{aligned} \chi = \left( \frac{\lambda _\mathrm {2}}{\lambda _{\mathrm {eq}}}-1\right) / \left( \frac{\lambda _\mathrm {2}}{\lambda _\mathrm {1}}-1\right) . \end{aligned}$$
(17)

The projection of the two-layer model with efficacy onto two regions is not limited to Heaviside step forcing input as it is independent of the form of the forcing.

Fig. 3
figure 3

Ocean heat uptake efficacy and pattern effect. a The regional feedback parameters \(\lambda _\mathrm {1}\) and \(\lambda _\mathrm {2}\) as a function of the efficacy factor \(\varepsilon\) for \(\lambda _{\mathrm {eq}} =-1.3\hbox { W m}^{-2}\hbox { K}^{-1}\), \(\eta =0.7\hbox { W m}^{-2}\hbox { K}^{-1}\) while \(C=10\hbox { W year m}^{-2}\hbox { K}^{-1}\) and \(C_\mathrm {D}=100\hbox { W year m}^{-2}\hbox { K}^{-1}\). b The projection of CMIP5 model output onto the regional feedback parameters \(\lambda _\mathrm {1}\) and \(\lambda _\mathrm {2}\)

Full size image

Based on an idealized configuration, in Fig. 3 the feedback parameters \(\lambda _\mathrm {1}\) and \(\lambda _\mathrm {2}\) are plotted as a function of ocean heat uptake efficacy \(\varepsilon\). An efficacy factor below unity implies a feedback parameter \(\lambda _\mathrm {1}\) of smaller magnitude compared to a strong feedback parameter \(\lambda _\mathrm {2}\) in the second region. In this case, \(\lambda _\mathrm {1}\) is actuated by fast temperature adjustment, while \(\lambda _\mathrm {2}\) is actuated by slow temperature adjustment on a slow timescale which originally characterized the adjustment due to the deep-ocean component. An efficacy factor above unity, indicated by the MPI-ESM1.2 abrupt CO\(_2\) experiments, implies a strong feedback parameter \(\lambda _\mathrm {1}\) compared to a feedback parameter \(\lambda _\mathrm {2}\) of smaller magnitude in the second region. Likewise, \(\lambda _\mathrm {1}\) is actuated by fast temperature adjustment, while \(\lambda _\mathrm {2}\) is actuated by slow temperature adjustment. In the case of unitary efficacy, the analytical solution is still characterized by a fast timescale \(\tau _\mathrm {f}\) and a slow timescale \(\tau _\mathrm {s}\). However, the relationship between T and N in the two-layer model would be linear with slope \(\lambda _{\mathrm {eq}}\) and no pattern of surface warming and feedback could be resolved. That is, the fast and the slow adjustment timescales must be directly coupled to radiative feedbacks to cause time-variation of \(\lambda (t)\). In Fig. 3b we plot the range of \(\lambda _\mathrm {1}\) and \(\lambda _\mathrm {2}\) for output from the Coupled Model Intercomparison Project (CMIP5) (\(4\times \text {CO}_2\)) using parameter estimates from Geoffroy et al. (2013b) who fit the two-layer model to 16 CMIP5 AOGCMs. In line with the MPI-ESM1.2 experiments, most of these models exhibit an efficacy factor above unity (\(C_\mathrm {1} \ll C_\mathrm {2}\) and \(\lambda _\mathrm {1}<\lambda _2\)). However, these models differ in their temporal adjustment, which is discussed in the following.

Fig. 4
figure 4

The fast timescale \(\tau _{\mathrm {f}}\) and the slow timescale \(\tau _{\mathrm {s}}\) as a function of the magnitude of the regional feedback parameter using the projection of the two-layer model onto the two-region model as well as CMIP5 output. For the idealized case (curves) we set \(\lambda _{\mathrm {eq}} =-1.3\hbox { W m}^{-2}\hbox { K}^{-1}\), \(\eta =0.7\hbox { W m}^{-2}\hbox { K}^{-1}\) while \(C=10\hbox { W year m}^{-2}\hbox { K}^{-1}\) and \(C_\mathrm {D}=100\hbox { W year }\hbox {m}^{-2}\hbox { K}^{-1}\)

Full size image

According to the projection of the two-layer model onto the two-region model, we plot the fast timescale \(\tau _{\mathrm {f}}\) and the slow timescale \(\tau _s\) against the magnitude of \(\lambda _{\mathrm {1}}\) and \(\lambda _{\mathrm {2}}\) (Fig. 4). Again, we use an idealized configuration of the two-layer model (curves) and vary the efficacy factor to compute regional feedbacks. Furthermore, we plot the pairs of regional feedbacks and timescales for the output from CMIP5 models. The regional approach establishes a simple relationship between regional feedbacks and regional timescales, \(\tau (r)=-C(r)/\lambda (r)\). That is, the regional timescales are inversely proportional to the feedback parameters. The question arises whether CMIP5 models describe a unique relationship or wether the relationship is appropriate to characterize complex climate model behavior. For a given configuration, the relationship between \(\lambda (r)\) and \(\tau (r)\) is independent of the global mean equilibrium feedback parameter: using the idealized case, changing the equilibrium feedback parameter \(\lambda _{\mathrm {eq}}\) changes the regional feedback parameters and timescales in such a way that we move along the same curves to higher or lower values. Using CMIP5 output, the fast timescales are approximately linearly related to the regional feedback parameters in the region of fast adjustment. The effective heat capacity of this region is small and \(\tau _{\mathrm {f}}\) is determined by the magnitude of \(\lambda _{\mathrm {1}}\). This is commonly observed in regions which are weakly coupled to the state of the deep-ocean as found in low latitudes. As far as the slow timescale is concerned, the idealized configuration shows an exponential increase of \(\tau _{\mathrm {s}}\) as \(\lambda _{\mathrm {2}}\) gets more positive. The CMIP5 models deviate from the idealized configuration because they differ in C, \(C_{\mathrm {D}}\) and \(\eta\), and these model discrepancies induce changes in the temporal behavior of the climate response. The variation of these intertia parameters would change the relationship between \(\tau (r)\) and \(\lambda (r)\), and this way we would move along a different curve. The model parameters fitted to CMIP5 models, however, do support our general understanding that a more sensitive region adjusts on a longer timescale, and that this timescale changes in a nonlinear way as the magnitude of the radiative feedback is reduced.

How valid is such a regional approach or two-region model conceptually? In their experimental model study based on aquaplanet simulations Haugstad et al. (2017) show that local TOA radiative feedbacks depend on the pattern of the climate forcing modified by ocean heat uptake but the same radiative resonse arises when the surface temperature pattern induced by that forcing is prescribed. In that respect, ocean heat uptake induces a surface temperature pattern but it is secondary how this surface temperature pattern is induced because the same radiative feedbacks govern the relationship between N and T. We have shown theoretically that, from the perspective of global N and global T, there is no difference in the radiative response and temperature evolution between the two-layer model and the two-region model. Whereas the two-layer and two-region approaches are mathematically equivalent, the former may be attractive in some cases simply because it does not predicate a fixed spatial distribution of regional feedbacks.

4 Feedback temperature dependence

The concepts discussed in the previous section describe a temporal adjustment on forcing-invariant timescales. In this section, following Bloch-Johnson et al. (2015), we add an explicitly state-dependent component and analyze conceptually the implications of feedback temperature dependence for the transient response of the climate system. We first provide insight into the analytical solution of a single mixed-layer to demonstrate that feedback temperature dependence introduces a timescale that depends on the strength of the forcing. Finally, we analyze the interplay between the two-timescale approach and feedback temperature dependence in the case of the two-layer ocean model, since the relationship between N and T in MPI-ESM1.2 suggests both a time-dependent feedback due to the interaction of different state-variables and state-dependent feedback due to feedback temperature dependence.

4.1 Forcing-dependent timescale

Radiative feedbacks may change in strength as the surface temperature T increases, and this temperature dependence can cause a nonlinear temperature response. The first-order feedback framework \({\varDelta }N= \frac{\partial N}{\partial T} {\varDelta }T\) can be extended by retaining terms up to second-order in the Taylor Series expansion of N. This results in nonlinear relationships as

$$\begin{aligned} {\varDelta }N= \frac{\partial N}{\partial T} {\varDelta }T + \frac{1}{2} \frac{\partial ^2 N}{\partial T^2} {\varDelta }T^2 \end{aligned}$$
(18)

gives a quadratic relationship; i.e., the quadratic coefficient \(\frac{1}{2}\frac{\partial ^2 N}{\partial T^2} \vert _{T}\) relates the change in the radiative feedback to T. We call the quadratic coefficient feedback temperature dependence and denote it by a. Using the zero-dimensional model Eq. (1) extended by second-order feedback temperature dependence and solving for the MPI-ESM1.2 equilibrium response for \(2\times \text {CO}_2\), \(4\times \text {CO}_2\) and \(8\times \text {CO}_2\) gives a positive quadratic coefficient a of roughly 0.04 W m\(^{-2}\hbox { K}^{-2}\). In the following we assume positive feedback temperature dependence.

We add second-order feedback temperature dependence to a time-dependent mixed-layer approximation to explore associated changes in the temporal behavior. In general, a mixed-layer approximation can be obtained by combining Eqs. (1) and (4). Alternatively, one can assume that \(C_\mathrm {D} \rightarrow \infty\) such that the two-layer ocean model Eq. (5, 6) reduces to the globally-averaged one-layer model with deep ocean heat uptake.

$$\begin{aligned} C \frac{\mathrm {d}T}{\mathrm {d}t}= F+\lambda T -\varepsilon \eta T. \end{aligned}$$
(19)

In Eq. (19) heat is removed from the mixed layer at constant rate \(\varepsilon \eta\), and mathematically there is no difference to the time-dependent formulation of the linear forcing-feedback framework which decribes each zone in the two-region model (\(\varepsilon \eta\) acts as an additional radiative response parameter). To avoid confusion for later analysis, we aggregate first-order radiative response coefficients in \({\varLambda }= \lambda -\varepsilon \eta\). Again, we assume Heaviside step function input for the radiative forcing F. Adding the quadratic coefficient representing feedback temperature dependence, our problem is described by a first-order nonlinear differential equation with constant coefficients;

$$\begin{aligned} C \frac{\mathrm {d}T}{\mathrm {d}t}= F+ {\varLambda }T + a T^2. \end{aligned}$$
(20)

Substituting

$$\begin{aligned} T= - \frac{C}{a} \frac{\mathrm {d} S/\mathrm {d} t}{S} \end{aligned}$$
(21)

this equation can be transformed into the second-order linear differential equation

$$\begin{aligned} \frac{\mathrm {d}^2S}{\mathrm {d}t^2}= \frac{1}{C}\left( {\varLambda }\frac{\mathrm {d}S}{\mathrm {d}t} - \frac{aF}{C}S\right) . \end{aligned}$$
(22)

After re-transformation and applying the initial condition \(T(0)=0\) the analytical solution reads

$$\begin{aligned} \begin{aligned} T(t)= C\frac{\alpha _\mathrm {2}}{a} \frac{e^{\alpha _\mathrm {1} t}-e^{\alpha _\mathrm {2} t}}{\frac{-\alpha _\mathrm {2}}{\alpha _\mathrm {1}} e^{\alpha _\mathrm {1} t}+e^{\alpha _\mathrm {2} t}} \,\, \text {or} \,\, T(t)= C \frac{\alpha _\mathrm {2}}{a} \frac{1-e^{{t}/{\tau }}}{\frac{\alpha _\mathrm {2}}{\alpha _\mathrm {1}}e^{{t}/{\tau }}-1}\\ \text {where} \, \, \alpha _{\mathrm {1,2}}=\frac{{\varLambda }}{2C} \pm \frac{\sqrt{{\varLambda }^2-4 a F}}{2C} \, \, \text {and}\, \, \tau =\frac{C}{\sqrt{{\varLambda }^2-4 a F}}. \end{aligned} \end{aligned}$$
(23)

Stability depends additionally on the radiative forcing F and the quadratic coefficient a, and we get a quadratic runaway in the relationship between N and T in the case that \(F>\frac{{\varLambda }^2}{4a}\). Likewise, the characteristic timescale \(\tau\) on which temperature adjusts depends additionally on F and a, and changes in F or a lead to temporal changes in the adjustment of the system. In this connection, the temporal beavior of the climate response is determined by the strength of the forcing.

A stronger step forcing or a more positive a causes the equilibrium response to be reached far later in time compared to a configuration without feedback temperature dependence. In MPI-ESM1.2 the \(8\times \text {CO}_2\) and \(16\times \text {CO}_2\) experiments are still far from stationarity after 1000 years of integration, while the \(2\times \text {CO}_2\) experiment is nearly equilibrated. Figure 5 shows an example of forcing-induced changes in the normalized timescale of the mixed-layer model with feedback temperature dependence against \(a \rightarrow 0\); i.e.,

$$\begin{aligned} \tau /\tau _0=\frac{-{\varLambda }}{\sqrt{{\varLambda }^2-4 a F}}. \end{aligned}$$
(24)
Fig. 5
figure 5

The normalized change of the characteristic timescale \(\tau\) in the mixed-layer model as a function of F. The radiative response coefficient \({\varLambda }\) is set to \({\varLambda }=-1.3 \hbox { W m}^{-2}\hbox { K}^{-1}\)

Full size image

In this case, the choice of the effective heat capacity is arbitrary because the timescale \(\tau\) is proportional to C. The characteristic timescale increases in a nonlinear way as the forcing is increased, and this nonlinearity is reinforced by a more positive quadratic coefficient. The positive contribution from feedback temperature dependence to the temperature or radiative response is counteracted by the first-order radiative response coefficients comprised in \({\varLambda }\).

4.2 Two timescales and feedback temperature dependence

The one-layer model with feedback temperature dependence adjusts on a single timescale and the relationship between N and T is characterized by a continuous function with smoothly changing slope which can describe the long-term response in MPI-ESM1.2 on centennial timescales. The change in the relationship between N and T on a decadal timescale, however, cannot be reproduced due to the missing multiple timescale structure. We combine the two-timescale approach and state-dependent feedback by adding feedback temperature dependence to the mixed-layer energy budget of the two-layer ocean model Eqs. (5, 6). The TOA imbalance N is now given by

$$\begin{aligned} N_{\varepsilon \ne 1}= F+ \lambda _\mathrm {b} T + a T^2 - (\varepsilon -1) \eta (T-T_\mathrm {D}) \end{aligned}$$
(25)

with \(\lambda _\mathrm {b}\) the background or initial feedback parameter: in Eq. (5) we replace \(\lambda _\mathrm {eq}\) by \(\lambda _\mathrm {T}(T)= \lambda _\mathrm {b} + a T\). Solving for stationary solutions \(-F=\lambda _\mathrm {b} T + a T^2\) gives equilibrium sensitivity which is still independent of the coupling between the layers.

We split the surface adjustment T into a fast response \(T_\mathrm {f}\) when the deep-ocean is not yet significantly perturbed and a slow component \(T_\mathrm {s}\) when the deep-ocean responds; that is, \(T = T_\mathrm {f}+ T_\mathrm {s}\).

Fig. 6
figure 6

The fast (dashed) and the slow (solid) mode of the two-layer ocean model with feedback temperature dependence, for a different feedback temperature dependencies a and for b different step function input for the radiative forcing F. In (A) the forcing is set to \(F=12\hbox { W m}^{-2}\) (\(\sim 3\times F_{2 \times \text{ CO }_2}\)). In b the quadratic coefficient is set to \(a=0.03\hbox { W m}^{-2}\hbox { K}^{-2}\). The background feedback parameter is \(\lambda _\mathrm {b}=-1.3\hbox { W m}^{-2}\hbox { K}^{-1}\), \(\eta =0.7\hbox { W m}^{-2}\hbox { K}^{-1}\), \(\varepsilon =1.5\), \(C=10\hbox { W yr m}^{-2}\hbox { K}^{-1}\) and \(C_\mathrm {D}=100\hbox { W yr m}^{-2}\hbox { K}^{-1}\)

Full size image

This allows us to characterize the interplay between the two-timescale approach and feedback temperature dependence. The fast timescale, on the order of a few years for reasonable parameters, and the related fast response \(T_\mathrm {f}\) can be approximated by

$$\begin{aligned}&\tau _\mathrm {f} \approx {C}/{\sqrt{{\varLambda }^2-4aF}} \,\, \text {and} \nonumber \\&\quad T_{\mathrm {f}} \approx \frac{-{\varLambda }- \sqrt{{\varLambda }^2-4aF}}{2 a} \,\, \text {with} \,\, {\varLambda }={\lambda _\mathrm {b}-\varepsilon \eta }. \end{aligned}$$
(26)

We assume that \(T_\mathrm {f}\) is constant due to the fast adjustment of the mixed-layer; making the fast adjustment time-dependent leads to Eq. (21) with the analytical solution described by Eq. (23).  On timescales longer than \(\tau _\mathrm {f}\), the slow component is approximated by

$$\begin{aligned} T_\mathrm {s} (t) \approx \frac{\sqrt{{\varLambda }^2 - 4 a F}-{\sqrt{{\varLambda }^2- 4 a F - 4 a \varepsilon \eta T_\mathrm {D} (t)}}}{2a} \end{aligned}$$
(27)

where the term \(a \varepsilon \eta T_\mathrm {D} (t)\) is the time-dependent control. Equations (26) and (27) indicate that there is no linear superposition between the two-timescale approach and feedback temperature dependence. State-dependent feedback and time-dependent feedback are entangled, and thus indistinguishable in a single forcing realization. If time-dependent feedback and state-dependent feedback were not coupled, we could configure a purely time-dependent model and a purely state-dependent model to obtain the same output.

The choice of parameters, however, determines to what extent the nonlinear behavior of the climate response is related to the fast mode \(T_\mathrm {f}\) and to the slow mode \(T_\mathrm {s}\). We use an idealized case and plot \(T_\mathrm {f}\) and \(T_\mathrm {s}\) for different feedback temperature dependencies (Fig. 6a) and for different forcing input (Fig. 6b). For simplicity we choose \(F_{2\times \text {CO}_2}=4\hbox { W m}^{-2}\) and take multiples thereof. Furthermore, we choose a high ocean heat uptake efficacy \(\varepsilon = 1.5\). In Fig. 6a the effect of the different feedback temperature dependencies on the fast response is counteracted by \({\varLambda }= \lambda -\varepsilon \eta\). The first-order coefficient \({\varLambda }\) comprises the high ocean heat uptake efficacy \(\varepsilon\) and \(T_\mathrm {f}\) becomes thus approximately independent of a. The effect of feedback temperature dependence on the radiative response is reinforced as the slow component evolves. In contrast, both the fast component and the slow component change significantly as the radiative forcing F is linearly increased (Fig. 6b). The fast response depends sensitively on the magnitude of \({\varLambda }\), and time-dependent feedback dictates to what extent \(T_\mathrm {f}\) and \(T_\mathrm {s}\) behave nonlinearly as we change F. In the case of an efficacy factor below unity, the magnitude of the slow mode \(T_\mathrm {s}\) is reduced compared to an efficacy factor above unity. Likewise, the degree to which the slow mode \(T_\mathrm {s}\) behaves nonlinearly diminishes, because the radiative response is initially suppressed and subsequently enhanced.

Fig. 7
figure 7

Time-dependent feedback and state-dependent feedback in MPI-ESM1.2 based on the two-layer ocean model with feedback temperature dependence. The background feedback parameter is set to \(\lambda _\mathrm {b}=-1.65\hbox { W m}^{-2}\hbox { K}^{-1}\), \(a=0.04\hbox { W m}^{-2}\hbox { K}^{-2}\), \(\eta =0.8\hbox { W m}^{-2}\hbox { K}^{-1}\), \(\varepsilon =1.2\), \(C=5\hbox { W year m}^{-2}\hbox { K}^{-1}\) and \(C_\mathrm {D}=150\hbox { W year m}^{-2}\hbox { K}^{-1}\)

Full size image

Nonconstant global feedback in MPI-ESM1.2 can be explained by the combined effect of time-dependent feedback and state-dependent feedback. We can determine purely time-dependent feedback and purely state-dependent feedback based on parameter estimates for the two-layer ocean model fitted to AOGCM output. Parameter estimates are taken from Mauritsen et al. (2018) who fit the two-layer model with feedback temperature dependence to the different MPI-ESM1.2 abrupt CO\(_2\) experiments. In Fig. 7 we plot time-dependent feedback (\(a=0\) and \(\varepsilon \ne 1\)) and state-dependent feedback (\(a \ne 0\) and \(\varepsilon =1\)) in these experiments. Time-dependent feedback can explain the abrupt change in the relationship between N and T in the MPI-ESM1.2 on decadal timescales as well as the variation in the \(2\times \text {CO}_2\) experiment. For larger forcing such as \(4\times \text {CO}_2\) and \(8\times \text {CO}_2\), the effect of feedback temperature dependence on the radiative response becomes increasingly important and dictates the long-term response. We could quantify contributions from a and \(\varepsilon\) to the change of N or the variation of \(\lambda (t)\) in the combined run (\(a \ne 0\) and \(\varepsilon \ne 1\)). However, these contributions are relative to the temperature adjustments of the system while the characteristic timescales on which the system adjust themselves depend on F; i.e., these contributions are related to the combined effect of time-dependent and state-dependent feedback. Interpreting time-dependent feedback in AOGCM experiments when both time-dependent feedback and state-dependent feedback are present thus depends sensitively on the conceptual framework chosen and the assumptions made on the temperature-dependent feedback parameter.

Fig. 8
figure 8

The timescale \(\tau (t)\) for time-dependent feedback and state-dependent feedback in MPI-ESM1.2 based on the two-layer ocean model with feedback temperature dependence. The background feedback parameter is set to \(\lambda _\mathrm {b}=-1.65\hbox { W m}^{-2}\hbox { K}^{-1}\), \(a=0.04\hbox { W m}^{-2}\hbox { K}^{-2}\), \(\eta =0.8\hbox { W m}^{-2}\hbox { K}^{-1}\), \(\varepsilon =1.2\), \(C=5\hbox { W yr m}^{-2}\hbox { K}^{-1}\) and \(C_\mathrm {D}=150\hbox { W yr m}^{-2}\hbox { K}^{-1}\). In the case of time-dependent feedback, the lines simply lie atop one another

Full size image

To highlight the differences in the temporal adjustment between state-dependent and time-dependent models as well as their combination, we define a timescale \(\tau (t)\) which is allowed to vary over time. We define this timescale as \(\tau (t)= -t/\ln ( \frac{T(t)-T_{\mathrm {eq}}}{-T_{\mathrm {eq}}})\) where T(t) is the globally averaged surface temperature perturbation and \(T_{\mathrm {eq}}\) is the equilibrium response. That is, we assume a single mode. Figure 8 shows \(\tau (t)\) for purely time-dependent feedback, purely state-dependent feedback and their fitted combination. In general, \(\tau (t)\) increases monotonically from relatively small values which are related to the fast mode of the climate response to relatively high values which are related to amplified warming due to the slow component on a centennial timescale. In the case of purely time-dependent feedback, the variation of \(\tau (t)\) over time is independent of the strength of the radiative focing F, since the characteristic timescales of the two-timescale approaches are constant and the amplitudes of the temperature adjustment are linear in forcing. In the case of purely state-dependent feedback due to feedback temperature dependence and unitary efficacy, the variation of \(\tau (t)\) depends additionally on the radiative forcing F. Due to the interaction of time-dependent feedback and state-dependent feedback the temporal adjustment is even prolonged and the effective timescales deviate from purely state-dependent feedback despite the same equilibrium response. Understanding the forcing-dependency of the characteristic adjustment and its interaction with time-dependent feedbacks is therefore important to narrow the timing of transient warming on a centennial timescale.

5 Conclusion

The two-timescale approaches give rise to time-dependent feedback but we paraphrase time-dependence as implicit because it emerges from the interplay of at least two state-variables and relates to the evolution of the radiative response over time rather than to an explicit time-dependence. The recent conceptual frameworks of the two-timescale approaches, the two-layer model with efficacy and the two-region model, can describe the same global temperature evolution and radiative response. On mathematical grounds, the two-layer ocean model with nonunitary efficacy can be projected onto a two-region model, and we can use the regional model feedbacks and timescales to characterize complex climate model behavior. The mathematical equivalence raises the question which specific mechanism causes time-dependent feedback. We cannot reject one or the other framework given only global mean output for the TOA imbalance N and surface temperature change T. We suggest that the radiative response of the climate system can be aggregated on different regions such as regions directly coupled to the state of the deep-ocean and not directly coupled to the state of the deep-ocean.

In constrast to time-dependent feedback, state-dependent feedback due to feedback temperature dependence introduces a nonlinearity which makes the timescale on which the climate system adjusts forcing-dependent. The finding that the timescales on which the climate system adjusts depend on the strength of the forcing has not been addressed to date. As we alter the forcing, we change the system’s response times and a certain fraction of the equilibrium response or temperature level is reached earlier or later in time. Further, this forcing-dependency and the combined effect or entanglement of time-dependent feedback and state-dependent feedback poses challenges for determining contributions from different sources such as feedback temperature dependence or a pattern effect. Single forcing realizations may not constrain the model parameters, and this complicates attempts to understand temperature change especially in more strongly forced climates on a centennial timescale but also the details of the instrumental record.