A mathematical approach to understanding emergent constraints

One of the approaches to constrain uncertainty in climate models is the identification of emergent constraints. These are physically explainable empirical relationships between a particular simulated characteristic of the current climate and future climate change from an ensemble of climate models, which can be exploited using current observations. In this paper, we develop a theory to understand the appearance of such emergent constraints. Based on this theory, we also propose a classification for emergent constraints, and applications are shown for several idealized climate models.

1 Introduction

Improving the accuracy of climate projections is one of the most important challenges in climate modelling. The uncertainty can be reduced by the development of more and more sophisticated global climate models, capturing more processes and scales. However, the societal importance of climate projections calls for a faster pace of improvement and alternative approaches that aim to better determine the accuracy of existing models. One of the proposed methods to accomplish this has been the use of so-called emergent constraints, where current observations are used to constrain future projections (Collins et al., 2012).

In multi-model ensembles of complex climate models, an apparent linear relation can be found between short-term and long-term changing variables. More credibility is attached to models that match the observed variability or trend well over the recent period. In this way, current observations provide a constraint to long-term trends. The observed variable is called the predictor, while the variable that is to be constrained is called the predictand (Klein and Hall, 2015). In recent years, emergent constraints have been found for Arctic warming, snow-albedo feedback, tropical carbon, the global precipitation among other variables (Allen and Ingram, 2002; Bracegirdle and Stephenson, 2013; Hall and Qu, 2006; Wenzel et al., 2014) and more recently, climate sensitivity (Cox et al., 2018).

A prominent example is the emergent constraint found in Hall and Qu (2006) where an emergent relationship was found between the strength of the snow-albedo feedback (SAF) on a seasonal timescale and the SAF under global warming in a Coupled Model Intercomparison Project phase 3 (CMIP3) ensemble. They also elucidated the key physical process behind the emergent constraint. Models where the maximum albedo of snow is highest have the largest SAF on both timescales because the contrast between snow-covered and snow-free areas is high (Qu and Hall, 2007).

However, a more general dynamical picture on how emergent constraints occur in multi-model ensembles or even in a parameter ensemble of a single model is still lacking. Under which circumstances are these constraints expected to arise? Some emergent constraints may be spurious and could arise because of shared errors in a particular multi-model ensemble (Bracegirdle and Stephenson, 2013). A mathematical framework is desired to identify spurious constraints and to give an indication as to where new emergent constraints might arise.

Here, we investigate how and under what conditions emergent constraints appear and what can be learned about the physics of the climate system. We will use linear response theory (LRT) to address the problem of forcing-response relations on different timescales (Risken, 1996). Ruelle demonstrated that LRT can be extended to study the response of non-equilibrium systems to external forcing. As with the fluctuation–dissipation theorem, Ruelle's LRT uses the statistical properties of the unforced (equilibrium) state only, but it does not assume (quasi)-equilibrium. Recently, LRT has been proposed as a rigorous framework for computing the response of the climate system and its applicability has been tested on the Lorenz-96 model and on the idealized global climate model PlaSim (Lucarini and Sarno, 2011; Ragone et al., 2016).

The paper is organized as follows. To obtain an understanding of emergent constraints, we start by formulating a mathematical framework in terms of susceptibilities by making use of LRT (Sect. 2). This results in explicit expressions for the appearance of emergent constraints in terms of susceptibility functions. In Sect. 3, a classification scheme for emergent constraints is proposed. Then, in Sect. 4, applications are presented for conceptual climate models, such as Ornstein–Uhlenbeck processes in one and two dimensions, an energy balance model and the PlaSim model. The results are summarized and discussed in Sect. 5.

2 Response functions

In this section, explicit expressions are given for response functions of the state of a dynamical system which depends on a single parameter and which is subjected to a non-stationary forcing. Such response functions are used in the following section to classify the different emergent constraints. Rigorous results for linear response properties of large class of general stochastic systems was obtained by Hairer and Majda (2010) (linear response theory for non-equilibrium systems was developed by Ruelle, 1998, 2009).

We illustrate the approach using the general one-dimensional forced stochastic differential equation (SDE):

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }\mathrm{d}{X}_t=\left(-V^{\prime }\left(X_t\right)+F\left(t\right)\right)\mathrm{d}t+\sqrt{\mathit{\sigma }}\mathrm{d}{W}_t.\end{array}$$

Here, V(x) is a smooth confining potential, meaning that a equilibrium solution exists for the unforced system (Pavliotis, 2014), and F(t) is a prescribed forcing. Furthermore, σ is the noise amplitude and the associated Wiener process is indicated by W_t. Usually, the potential depends on a parameter.

The probability density function of the unforced (F(t)=0) system, say $\stackrel{\mathrm{‾}}{p}$, satisfies the Fokker–Planck equation:

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{p}}{\partial t}}={\displaystyle \frac{\partial \left(V^{\prime }\right(x\left)\stackrel{\mathrm{‾}}{p}\right)}{\partial x}}+{\displaystyle \frac{\mathit{\sigma }}{\mathrm{2}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{‾}}{p}}{\partial x^{\mathrm{2}}}}={\mathcal{L}}^{*}\left(\stackrel{\mathrm{‾}}{p}\right),\end{array}$$

which defines the Fokker–Planck operator ℒ^*. The equilibrium distribution of the unforced system, here indicated by ${\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}$, is given by

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}\left(x\right)={\displaystyle \frac{\mathrm{1}}{Z}}{e}^{\frac{-\mathrm{2}V\left(x\right)}{\mathit{\sigma }}};Z=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{\frac{-\mathrm{2}V\left(x\right)}{\mathit{\sigma }}}\mathrm{d}x.\end{array}$$

Linear response theory (Ragone et al., 2016) provides an expression for the change in the expectation value of the change in an observable O (e.g. the temperature, ice extent or the standard deviation of either), say ΔO(t) when the system is forced, compared to the unforced case, i.e.

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }O\left(t\right)=E\left[O\left(X_t\right)\right]-E\left[O_{\mathrm{e}}\left(X_t\right)\right],\end{array}$$

where again the subscript “e” indicates the equilibrium of the unforced system. It follows that

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{\Delta }O\left(t\right)=\underset{\mathrm{0}}{\overset{t}{\int }}{R}_{\mathrm{O}}(t-s)F\left(s\right)\mathrm{d}s;\\ \text{(5)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{O}}\left(t\right)=\mathcal{H}\left(t\right)\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}O\left(x\right)e^{{\mathcal{L}}^{*}t}\left(-{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}{\partial x}}\right)\mathrm{d}x,\end{array}$$

where R_O(t) is the response function, which is extended to be zero for t<0 to ensure causality with a Heaviside function ℋ(t). When Eq. (5) is Fourier transformed, we find, using the convolution theorem,

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }\mathcal{F}\left(\mathrm{\Delta }O\right(t\left)\right)\left(\mathit{\omega }\right)=\mathit{\chi }\left(\mathit{\omega }\right)\widehat{F}\left(\mathit{\omega }\right),\end{array}$$

where the Fourier transform χ(ω) of the response function R_O(t) is the susceptibility. If we take a cosine forcing, i.e. F(t)=F₀cosω₀t, then $\widehat{F}\left(\mathit{\omega }\right)=F_{\mathrm{0}}\mathit{\pi }\left(\mathit{\delta }\right(\mathit{\omega }-{\mathit{\omega }}_{\mathrm{0}})+\mathit{\delta }(\mathit{\omega }+{\mathit{\omega }}_{\mathrm{0}}\left)\right)$, so once we know χ(ω), we can determine the response ΔO(t).

In the Appendix, it is shown that when we take the identity operator O=x as the observable, thus taking the mean value of this variable, the response function and its corresponding susceptibility can be written as

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }{R}_{\mathrm{O}}\left(t\right)={\displaystyle \frac{\mathrm{2}}{\mathit{\sigma }}}\sum _{l=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\beta }}_l{e}^{-{\mathit{\lambda }}_{l}t},\mathit{\chi }\left(\mathit{\omega }\right)={\displaystyle \frac{\mathrm{2}}{\mathit{\sigma }}}\sum _{l=\mathrm{1}}^{\mathrm{\infty }}{\displaystyle \frac{{\mathit{\beta }}_l}{{\mathit{\lambda }}_l+i\mathit{\omega }}},\end{array}$$

where λ_l indicates the eigenvalues of the so-called generator ℒ,

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }\mathcal{L}u=V^{\prime }\left(x\right){\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\mathit{\sigma }}{\mathrm{2}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial x^{\mathrm{2}}}},\end{array}$$

and the β_l represents projection coefficients that indicate how strongly the system responds to the forcing; see the Appendix for a more detailed description.

The amplitude A of the response to a periodic forcing F(t)=F₀cosω₀t is determined by the absolute value of the susceptibility:

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }A\left(\mathrm{\Delta }X\right(t\left)\right)\left({\mathit{\omega }}_{\mathrm{0}}\right)={\displaystyle \frac{\mathrm{2}{F}_{\mathrm{0}}}{\mathit{\sigma }}}\sum _{l=\mathrm{1}}^{\mathrm{\infty }}{\displaystyle \frac{{\mathit{\beta }}_l}{\sqrt{{\mathit{\lambda }}_l^{\mathrm{2}}+{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}}}}.\end{array}$$

When the predictor is a response to a uniform frequency forcing, it can be expressed as in Eq. (9). The predictand will generally also have this form, as a forcing with a long timescale can be approximated by a low-frequency forcing. The previous analysis can be generalized to more dimensions. In two dimensions, for example, with a state vector ${\mathit{Y}}_t=(Y_{\mathrm{1}t},Y_{\mathrm{2}t}{)}^T$, the SDE becomes

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }\mathrm{d}{\mathit{Y}}_t=\left(-\mathrm{\nabla }V\left({\mathit{Y}}_t\right)+F\left(t\right)\widehat{i}\right)\mathrm{d}t+\sqrt{\mathit{\sigma }}{\mathbf{I}}_{\mathrm{2}}\mathrm{d}{W}_t,\end{array}$$

where the term $F\left(t\right)\widehat{i}$ denotes a forcing in the direction of the first variable and I₂ the identity matrix. As shown in Pavliotis (2014), the derivation of the response function follows the one-dimensional case closely, resulting in

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle }{R}_{Y_{\mathrm{1}}}\left(t\right)={\displaystyle \frac{\mathrm{2}}{\mathit{\sigma }}}\sum _{l=\mathrm{1}}^{\mathrm{\infty }}{g}_l{e}^{-{\mathit{\lambda }}_{l}t};\text{(12)}& {\displaystyle }& {\displaystyle }{R}_{Y_{\mathrm{2}}}\left(t\right)={\displaystyle \frac{\mathrm{2}}{\mathit{\sigma }}}\sum _{l=\mathrm{1}}^{\mathrm{\infty }}{h}_l{e}^{-{\mathit{\lambda }}_{l}t},\end{array}$$

where g_l and h_l are again projection coefficients; g_l and h_l contain a term describing strength of the response in Y₁ and Y₂, respectively. The derivation of the exact terms is given in the Appendix. Calculating the response is analogous to the one-dimensional case, so that the Fourier transforms of the response functions are given by

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle }A\left(\mathrm{\Delta }{Y}_{\mathrm{1}}\left(t\right)\right)\left({\mathit{\omega }}_{\mathrm{0}}\right)={\displaystyle \frac{\mathrm{2}}{\mathit{\sigma }}}\sum _{l=\mathrm{1}}^{\mathrm{\infty }}{\displaystyle \frac{g_l}{\sqrt{{\mathit{\lambda }}_l^{\mathrm{2}}+{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}}}};\text{(14)}& {\displaystyle }& {\displaystyle }A\left(\mathrm{\Delta }{Y}_{\mathrm{2}}\left(t\right)\right)\left({\mathit{\omega }}_{\mathrm{0}}\right)={\displaystyle \frac{\mathrm{2}}{\mathit{\sigma }}}\sum _{l=\mathrm{1}}^{\mathrm{\infty }}{\displaystyle \frac{h_l}{\sqrt{{\mathit{\lambda }}_l^{\mathrm{2}}+{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}}}}.\end{array}$$

Note that, generalizing to uncoupled multidimensional systems, the eigenfunctions are found to be the tensor products of the eigenfunctions in the one-dimensional case, while the corresponding eigenvalues are the sum of the eigenvalues in the one-dimensional case.

3 Classification of emergent constraints

Although a wide set of different emergent constraints has been found, no attempts have been made to classify them so far using dynamical criteria. Here, a classification is proposed based on the time characteristics of the predictor and on the relationship between the predictor and the predictand. Using this classification, assessment of their applicability becomes easier. Furthermore, a classification is a prerequisite for a dynamical description of emergent constraints.

Figure 1 The emergent constraint on snow-albedo feedback $\frac{\mathrm{\Delta }{\mathit{\alpha }}_{\mathrm{s}}}{\mathrm{\Delta }{T}_{\mathrm{s}}}$ (from Hall and Qu, 2006, α_s given in units of percent). This is an example of a direct emergent constraint (it links the SAF in both past and future time) and a dynamical emergent constraint (it uses a response to a seasonal forcing as its predictor).

Firstly, an emergent constraint can be either direct or indirect. In the direct case, the predictor and predictand are the same observable, while in the indirect case they are not. In the latter case, the predictor variable and predictand variable have to be closely linked, for instance, via a physical process. We make a further distinction between static and dynamic emergent constraints. In a dynamic emergent constraint, a response to a known, or sometimes even unknown, forcing in the (present-day) predictor is linked to the response of the (future) predictand under the same (or a similar) forcing. For example, the forcing can be the annual cycle of solar radiation but can also be caused by ENSO or historical climate change. In a static emergent constraint, a relationship between the time-independent quantity of the unforced system in the present-day (predictor) is linked to the response in a quantity under climate change.

Knutti et al. (2006)Hall and Qu (2006); Qu and Hall (2014)Boe et al. (2009)Clement et al. (2009)Trenberth and Fasullo (2010)Fasullo and Trenberth (2012)Bracegirdle and Stephenson (2013)Gordon and Klein (2014)Qu et al. (2014)Sherwood et al. (2014)Su et al. (2014)Wenzel et al. (2014)Tian (2015)Kwiatkowski et al. (2017)Cox et al. (2018)

Table 1 Application of our classification of emergent constraints to a selection of examples found in literature. DD is a direct dynamical constraint, DS is a direct static constraint, and IS is an indirect static constraint, while ID denotes indirect dynamical emergent constraints. Abbreviations are as follows: RH: relative humidity, ITCZ: Intertropical Convergence Zone, TOA: top of atmosphere, SH: Southern Hemisphere, ECS: equilibrium climate sensitivity, LLC: low-level cloud, SAF: snow-albedo feedback, SAT: surface air temperature; GMST: global mean surface temperature. The emergent constraint found by Trenberth and Fasullo (2010) seems to be spurious: no physical mechanism was proposed and it did not appear in different ensembles, such as CMIP5 (Grise et al., 2015).

As an illustration, we apply our classification to examples of emergent constraints found in the literature in Table 1. Although this is not a complete overview, examples are found of the four types of emergent constraints. There are many examples of direct dynamical constraints, such as the one involving the snow-albedo feedback shown in Fig. 1 (Hall and Qu, 2006). Dynamic direct emergent constraints are the most intuitive. As long as the variations in the predictor are of a sufficient amplitude compared to those of the predictand, a correlation between the predictor and predictand automatically points towards a common physical basis, for example, a common dynamical response to an external forcing. The direct static emergent constraint found by Bracegirdle and Stephenson (2013) makes use of spatial patterns. All of the indirect constraints involve equilibrium climate sensitivity as the predictand. Often the mean of some variable with some known bias in the model ensemble is linked to ECS. For instance, in Tian (2015), the asymmetry bias in the Intertropical Convergence Zone (ITCZ) is linked to climate sensitivity. An example of a dynamical indirect emergent constraint is provided by Cox et al. (2018), who relate a function of autocorrelation of global surface temperature to ECS. In this case, the short-time forcing is assumed to be caused by internal variability.

Based on the response function theory in Sect. 2, we further elaborate on the classification and also discuss conditions for each type of constraint for a dynamical system with varying parameters (which defines the ensemble of models).

For a direct dynamical emergent constraint, in the standard case of a linear relationship, the relation has the following form: Predictand = C^st× Predictor, where C^st is a constant that is independent of the parameter used to generate the ensemble of models. Rewriting this, the ratio of the responses to forcing of frequencies ω₁ and ω₂ should be constant over the (parameter) ensemble members e_i. For the simple case of two forcings that only differ in frequency, we find the condition from the ratio of the susceptibilities SR as

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle }\mathrm{SR}\left(e\right)={\displaystyle \frac{A\left(\mathrm{\Delta }O\right(t\left)\right)\left({\mathit{\omega }}_{\mathrm{2}}\right)}{A\left(\mathrm{\Delta }O\right(t\left)\right)\left({\mathit{\omega }}_{\mathrm{1}}\right)}}={\displaystyle \frac{\sum _{l=\mathrm{1}}^{\mathrm{\infty }}\frac{{\mathit{\beta }}_l}{\sqrt{{\mathit{\lambda }}_l^{\mathrm{2}}+{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}}}}}{\sum _{l=\mathrm{1}}^{\mathrm{\infty }}\frac{{\mathit{\beta }}_l}{\sqrt{{\mathit{\lambda }}_l^{\mathrm{2}}+{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}}}}}=C^{st}.\end{array}$$

One variable (B) can act as forcing to a second variable (O), while being itself forced externally (F). The predictor and predictand are then given by the quotient of the response functions of O and B. A further complication is that often the forcing patterns are not exactly the same for the short (F₂) and long (F₁) periodic forcing. In this case, (Eq. 15) has to be adjusted to

$$\begin{array}{}\text{(16)}& {\displaystyle }{\displaystyle }\mathrm{SR}\left(e\right)={\displaystyle \frac{A\left(\mathrm{\Delta }O\left(t\right){|}_{F_{\mathrm{2}}}\right)\left({\mathit{\omega }}_{\mathrm{2}}\right)}{A\left(\mathrm{\Delta }O\left(t\right){|}_{F_{\mathrm{1}}}\right)\left({\mathit{\omega }}_{\mathrm{1}}\right)}}{\displaystyle \frac{A\left(\mathrm{\Delta }B\left(t\right){|}_{F_{\mathrm{1}}}\right)\left({\mathit{\omega }}_{\mathrm{1}}\right)}{A\left(\mathrm{\Delta }B\left(t\right){|}_{F_{\mathrm{2}}}\right)\left({\mathit{\omega }}_{\mathrm{2}}\right)}}=C^{st}.\end{array}$$

This is further discussed in the example of the idealized energy balance model.

Physically, we expect that the same mechanism is responsible for the response at a short and long timescale to obtain this type of emergent constraint. The system should have response times smaller than the timescale of the forcing or equivalently: the generator should have eigenvalues λ larger than the frequency of the forcing. Naturally, the response times $\frac{\mathrm{1}}{\mathit{\lambda }}$ of the dominant processes are expected to be at least smaller than the timescale of the slow forcing $\frac{\mathrm{1}}{{\mathit{\omega }}_{\mathrm{1}}}$.

Mathematically, the ratio in Eq. (15) becomes 1 in the case that all eigenvalues λ_l are much larger than the forcing frequencies. Interestingly, the linear relation breaks down in the case that the fast forcing has the same order of magnitude as the eigenvalues of the dominant terms in the susceptibility. Under the assumption of a single dominant term in the susceptibility and a slow forcing with frequency ω₂→0, the first correction term to the slope-one linear relation between predictand y and predictor x, is cubic in x.

In the case of indirect dynamical emergent constraints, a relationship between a predictor Y₁ and a predictand Y₂ is found. Assuming the predictor Y₁ is again a response to some forcing, we can repeat the analysis above for direct constraints for a system of two dimensions, where a forcing is added in one direction. Mutatis mutandis, a condition very similar to Eq. (15) is found, as

$$\begin{array}{}\text{(17)}& {\displaystyle \frac{A\left(\mathrm{\Delta }{Y}_{\mathrm{1}}\left(t\right)\right)\left({\mathit{\omega }}_{\mathrm{2}}\right)}{A\left(\mathrm{\Delta }{Y}_{\mathrm{2}}\left(t\right)\right)\left({\mathit{\omega }}_{\mathrm{1}}\right)}}={\displaystyle \frac{\sum _{l=\mathrm{1}}^{\mathrm{\infty }}\frac{g_l}{\sqrt{{\mathit{\lambda }}_l^{\mathrm{2}}+{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}}}}}{\sum _{l=\mathrm{1}}^{\mathrm{\infty }}\frac{h_l}{\sqrt{{\mathit{\lambda }}_l^{\mathrm{2}}+{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}}}}}=C^{st},\end{array}$$

where g_l and h_l are defined as in Eq. (12). For an emerging constraint to exist, the projection terms of the different observables should thus change in a similar fashion under the change in parameter.

Static direct constraints link the mean of an observable (predictor) to a change in the system under a specific forcing (predictand). Note that the susceptibility only contains information about the response to such forcing. Even in the limit of ω→0, it denotes the linear response of the system, without any information on the mean state (Lucarini and Sarno, 2011). So, to derive the condition for a linear relationship, the mean $E\left[O_{\mathrm{e}}\right(X_t\left)\right]=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}O\left(x\right)\mathrm{d}x$ and the susceptibility at frequency ω₁ are used.

For static emergent constraints, the linear relationship between the predictand and the predictor is not expected to pass through the origin, since the predictor will in general be non-zero. Therefore, an additional term I is added to the ratio, denoting the intercept of the line between the predictor's mean state and the predictand. Instead, the susceptibility is compared to the mean state and the following condition is derived, where C^st should again be a constant that is independent of parameter(s), which are used to generate the ensemble:

$$\begin{array}{}\text{(18)}& {\displaystyle }{\displaystyle \frac{E\left[O_{\mathrm{1}t}\right]-I}{A\left(\mathrm{\Delta }{O}_{\mathrm{2}}\left(t\right)\right)\left({\mathit{\omega }}_{\mathrm{1}}\right)}}={\displaystyle \frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}{O}_{\mathrm{1}}\left(x\right)\mathrm{d}x-I}{\sum _{l=\mathrm{1}}^{\mathrm{\infty }}\frac{\mathrm{2}}{\mathit{\sigma }}\frac{h_l}{\sqrt{{\mathit{\lambda }}_l^{\mathrm{2}}+{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}}}}}=C^{st}.\end{array}$$

Again, C^st can either be positive or negative, depending on the physics under consideration. This equation is both valid for direct and indirect static emergent relationships; in the case of a direct constraint, O₁=O₂ and the term h_l contains information about the response of O₁ to a forcing, while in the indirect case O₁≠O₂ and h_l contains information about O₂.

As an illustration of the theory from Sect. 2 and a direct dynamical emergent constraint, we take the Ornstein–Uhlenbeck process (OU process). Here, $V^{\prime }\left(x\right)=\mathit{\gamma }x$, where γ is a parameter that indicates the steepness of the potential. The eigenvalues and eigenfunctions of the generator are given by (Pavliotis, 2014)

$$\begin{array}{}\text{(19)}& {\displaystyle }{\displaystyle }{\mathit{\lambda }}_l=\mathit{\gamma }l;{\mathit{\varphi }}_l\left(x\right)={\displaystyle \frac{\mathrm{1}}{\sqrt{l\mathrm{!}}}}{H}_n\left(\sqrt{{\displaystyle \frac{\mathrm{2}\mathit{\gamma }}{\mathit{\sigma }}}}x\right),\end{array}$$

where H_n are the Hermite polynomials. For the Ornstein–Uhlenbeck case, the ratio of response amplitudes reduces to

$$\begin{array}{}\text{(20)}& {\displaystyle }{\displaystyle }\mathrm{SR}\left(\mathit{\gamma }\right)={\displaystyle \frac{{\mathit{\beta }}_{\mathrm{1}}/\sqrt{{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}}}}{{\mathit{\beta }}_{\mathrm{1}}/\sqrt{{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}}}}={\displaystyle \frac{\sqrt{\mathrm{1}+{\left({\mathit{\omega }}_{\mathrm{1}}/\mathit{\gamma }\right)}^{\mathrm{2}}}}{\sqrt{\mathrm{1}+{\left({\mathit{\omega }}_{\mathrm{2}}/\mathit{\gamma }\right)}^{\mathrm{2}}}}},\end{array}$$

since both the observable and the derivative of the potential are orthogonal to all eigenfunctions other than ϕ₁. This ratio is dependent on γ. In the case of γ≫ω_i for $i\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$, this ratio is nearly 1 and an emergent relationship is present for a model ensemble generated by varying γ.

4 Application to idealized climate models

From the previous sections, it appears that the computation of the eigensolutions of the generator of the dynamical system are central to determine whether an emergent constraint will appear or not. In this section, we will provide examples using idealized climate models.

The eigenvalues and eigenfunctions of the generator were numerically determined using the fact that the eigenvalues of the Fokker–Planck operator ℒ^* are equal to those of the generator and that the eigenfunctions can be computed from the transformation ${\mathit{\varphi }}_l={\mathit{\varphi }}_l^{*}/{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}$. The Fokker–Planck operator was discretized with use of Chang–Cooper algorithm (Chang and Cooper, 1970). Eigenvalues and eigenvectors were determined using an implicitly restarted Arnoldi method (Lehoucq et al., 1998). Explicit simulations of the SDEs were performed using a stochastic Runge–Kutta method (Kloeden and Platen, 1992).

Figure 2 (a) Response to forcings at two different frequencies of the one-dimensional Ornstein–Uhlenbeck process. Shown is the average of a 500-member simulation of trajectories. (b) The susceptibility at these frequencies, whose ratio is given in the inset figure. This is an example of a direct dynamical emergent relationship.

Figure 3 Eigenvalue spectrum for (a) δ=0.2 and (b) δ=0.5. The dashed line corresponds to the frequency ω₂ of the fast forcing. (c, d) Corresponding susceptibilities, with their ratio in the inset figures. This is an example of an indirect dynamical emergent relationship. Note that for reasons of numerical stability, the range of γ₁ is different than that of γ in Fig. 2.

4.1 Ornstein–Uhlenbeck cases

First, the one-dimensional Ornstein–Uhlenbeck process is considered with SDE

$$\begin{array}{}\text{(21)}& {\displaystyle }{\displaystyle }\mathrm{d}{X}_t=\left(-\mathit{\gamma }{X}_t+F_i\left(t\right)\right)\mathrm{d}t+\sqrt{\mathit{\sigma }}\mathrm{d}{W}_t,\end{array}$$

forcing F_i(t)=sin2πtω_i and frequencies ω₁=0.001 and ω₂=0.1. A parameter ensemble is created by varying γ. In this case, analytic solutions exist for the eigenvalues and eigenvectors of the generator. Eigenvectors and eigenvalues were determined using the Chang–Cooper scheme on a domain [−25, 25] with Δx=0.25. The numerically computed susceptibilities, as shown in Fig. 2b, are in agreement with the analytic ones and capture the response (Fig. 2a) well, as expected in this linear case.

In the two-dimensional Ornstein–Uhlenbeck case, the same forcing F_i(t) is added but only in the first dimension. The governing SDE is given by

$$\begin{array}{}\text{(22)}& {\displaystyle }{\displaystyle }\mathrm{d}{\mathit{X}}_t=\left[\left(\begin{array}{cc}-{\mathit{\gamma }}_{\mathrm{1}}& \mathit{\delta }\\ \mathit{\delta }& -{\mathit{\gamma }}_{\mathrm{2}}\end{array}\right){\mathit{X}}_t+\left(\begin{array}{c}{F}_i\left(t\right)\\ \mathrm{0}\end{array}\right)\right]\mathrm{d}t+\sqrt{\mathit{\sigma }}\left(\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}\end{array}\right)\mathrm{d}{\mathit{W}}_t,\end{array}$$

and a parameter ensemble is generated by changing the damping rate γ₁. Two ensembles are compared with δ=0.2 in the first ensemble and δ=0.5 in the second. The damping term γ₂ is held constant at γ₂=0.6.

In Fig. 3, the eigenvalues and susceptibility ratios are plotted. In the case of a relatively weak coupling (δ=0.2), all non-zero eigenvalues are larger than the fast forcing frequency ω₂, so the system response time is smaller than the forcing timescales. On the other hand, the strong coupling (δ=0.5) leads to a slowdown of the system, so that some eigenvalues now become smaller than ω₂. In these cases (γ₁<0.5), the system does not have time to portray the full response to a forcing, while for others (γ₁>0.5) it does. Consequently, the strength of the response actually decreases for γ₁<0.5. Directly calculating the expectation value as the mean of 500 stochastic trajectories confirms this result (not shown).

Figure 4 (a, b) Projections g_l (of predictor variable) and h_l (of predictand variable) for a weakly coupled two-dimensional OU system with δ=0.2 and (c, d) the same for δ=0.5.

The results in Fig. 4 show a large variation over the ensemble in the projection term of the predictor on the eigenfunctions (g_l; see the Appendix). In contrast, the product of the two projection terms in the predictand (h_l) changes relatively little over the ensemble for both coupling strengths. Even though the projection terms now play a significant role in determining the response, the eigenvalues still determine whether the relation is linear (fast compared to forcing) or non-linear (similar size to forcing frequency). In the weak-coupling system, the susceptibility ratio is almost constant, and an emergent linear relationship is found. The strong-coupling system only portrays an emergent relationship for certain regimes (low or high γ₁). A case can be made though that the highly coupled system is the system for which finding an emergent constraint is more likely, because the strength of the response is substantially higher and a better signal-to-noise ratio can be obtained.

4.2 Energy balance model

In this section, a specific emergent constraint is examined in more detail, namely the one pertaining to the SAF first described by Hall and Qu (2006). They found a correlation between SAF on a seasonal scale and SAF as a result of climate change. In models with a high snow albedo, the contrast between snow-covered and bare surfaces was largest and consequently the sensitivity to changes in temperature was largest (Qu and Hall, 2007). To study this emergent constraint, we modify a simple energy balance model based on the seminal work by Budyko (1969) and Sellers (1969). The albedo is made temperature dependent following Fraedrich (1979), and a stochastic term is added following Sutera (1981). A parameter in the albedo function will be used to define a parameter ensemble.

Table 2 Constants for the energy balance model.

With constant albedo, the energy balance model reads

$$\begin{array}{}\text{(23)}& {\displaystyle }{\displaystyle }\mathrm{d}T={\displaystyle \frac{\mathrm{1}}{c_T}}\left(Q(\mathrm{1}-\mathit{\alpha })+A\ln {\displaystyle \frac{C}{C_{\mathrm{ref}}}}+G-\mathit{ϵ}{\mathit{\sigma }}_{\mathrm{B}}{T}^{\mathrm{4}}\right)\mathrm{d}t+\sqrt{{\mathit{\sigma }}_T}\mathrm{d}{W}_t,\end{array}$$

where dT is the temperature change, c_T the atmospheric heat capacity, Q the solar insolation, α the albedo, C the concentration of greenhouse gases, C_ref a reference concentration, G the radiative forcing due to the reference greenhouse gas concentration, σ_B the Stefan–Boltzmann constant and ϵ the emissivity of the Earth. The standard parameter values for this model can be found in Table 2. The parameters of the albedo function are chosen to ensure that no bistability is present in the model, in which case LRT would break down.

Before examining the snow-albedo feedback, note that, for some variables, notably the climate sensitivity, a simple energy balance model (EBM) can react differently to forcing from solar insolation or greenhouse gases. This can be determined from, with $H=G+A\ln \frac{C}{C_{\mathrm{ref}}}$ and for a value of ϵ=1,

$$\begin{array}{}\text{(24)}& {\displaystyle }& {\displaystyle \frac{\partial }{\partial \mathit{\alpha }}}{\displaystyle \frac{\partial T}{\partial Q}}={\displaystyle \frac{{\mathit{\sigma }}^{\mathrm{1}/\mathrm{4}}}{\mathrm{4}\left(Q\right(\mathrm{1}-\mathit{\alpha })+H{)}^{\mathrm{3}/\mathrm{4}}}}\left({\displaystyle \frac{\mathrm{3}Q(\mathrm{1}-\mathit{\alpha })}{Q(\mathrm{1}-\mathit{\alpha })+H}}-\mathrm{1}\right)<\mathrm{0};\text{(25)}& {\displaystyle }& {\displaystyle \frac{\partial }{\partial \mathit{\alpha }}}{\displaystyle \frac{\partial T}{\partial H}}={\displaystyle \frac{\mathrm{3}Q{\mathit{\sigma }}^{\mathrm{1}/\mathrm{4}}}{\mathrm{16}\left(Q\right(\mathrm{1}-\mathit{\alpha })+H){)}^{\mathrm{7}/\mathrm{4}}}}>\mathrm{0}.\end{array}$$

Sensitivity to greenhouse forcing decreases when albedo decreases, while sensitivity to solar insolation (seasonal sensitivity) increases for an increasing albedo, using typical values for Q and H.

To mimic the physical mechanism behind the emergent constraint, the albedo is taken to be temperature dependent; i.e. for low (high) temperatures, the albedo is high (low). A logistic function is used to model this effect:

$$\begin{array}{}\text{(26)}& {\displaystyle }{\displaystyle }\mathit{\alpha }\left(T\right)={\mathit{\alpha }}_{\min }+{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{amp}}}{\mathrm{1}+\exp \left(k\left(T-T_{\mathrm{h}}\right)\right)}},\end{array}$$

where α_min is the minimum albedo, α_amp is the amplitude, k is a steepness factor, and T_h is the temperature at which half of the amplitude is reached. The amplitude α_amp is the parameter that is varied over the ensemble.

In the first case, the insolation forcing is given by $Q=Q_{\mathrm{0}}(\mathrm{1}+Q_{\mathrm{s}}\sin \mathrm{2}\mathit{\pi }t/\mathit{\tau })$, where τ corresponds to 1 year and Q_s is a seasonal modulation amplitude, with parameter values shown in Table 2. The snow-albedo feedback term is then computed by dividing the amplitude of the albedo cycle by the amplitude of the temperature cycle. A second case is considered in which the greenhouse gas concentration C is increased 0.3 % per year from 295 ppmv over a period of 300 years. Here, the snow-albedo feedback is computed by dividing the total albedo response by the total temperature response. In each case, the variance of the noise σ_T in Eq. (23) was chosen as 10⁻⁷ K² s⁻¹. Changing this parameter does not influence the eigenvalues as expected from the theory (Pavliotis, 2014). While the projections of the eigenvalues and eigenfunctions did change slightly, the susceptibility ratio was not influenced significantly by a variation of the σ_T (halving and doubling of σ, not shown). In the computation of the solution of the Fokker–Planck equation using the Chang–Cooper scheme, we used a resolution of 1 K which is sufficient to accurately determine the eigenvalues and eigenfunctions of the generator.

As mentioned above, application of Eq. (15) is not self-evident. Considering temperature to be a forcing ignores the fact that temperature responds differently to seasonal and greenhouse gas forcing, as shown in Eq. (25). Secondly, using dα∕dT as the observable directly does not work either. Linear response theory does not give the expectation value of the observable, but the expectation value of the deviation due to the forcing, while we are interested in the change due to a parameter change.

Instead, the SAF can be described by two observables: SAF is determined by taking the ratio of the susceptibilities of albedo to temperature. Therefore, we use the modified Eq. (16):

$$\begin{array}{ll}{\displaystyle }\mathrm{RFS}\left({\mathit{\alpha }}_{\mathrm{amp}}\right)& {\displaystyle }={\displaystyle \frac{A\left(\mathrm{\Delta }\mathit{\alpha }\left(t\right){|}_Q\right)\left({\mathit{\omega }}_{\mathrm{2}}\right)}{A\left(\mathrm{\Delta }\mathit{\alpha }\left(t\right){|}_C\right)\left({\mathit{\omega }}_{\mathrm{1}}\right)}}:{\displaystyle \frac{A\left(\mathrm{\Delta }T\left(t\right){|}_C\right)\left({\mathit{\omega }}_{\mathrm{1}}\right)}{A\left(\mathrm{\Delta }T\left(t\right){|}_Q\right)\left({\mathit{\omega }}_{\mathrm{2}}\right)}}\\ \text{(27)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\sum _{l=\mathrm{1}}^{\mathrm{\infty }}\frac{{\mathit{\alpha }}_l}{\sqrt{{\mathit{\lambda }}_l^{\mathrm{2}}+{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}}}}}{\sum _{l=\mathrm{1}}^{\mathrm{\infty }}\frac{{\mathit{\gamma }}_l}{\sqrt{{\mathit{\lambda }}_l^{\mathrm{2}}+{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}}}}}:{\displaystyle \frac{\sum _{l=\mathrm{1}}^{\mathrm{\infty }}\frac{{\mathit{\delta }}_l}{\sqrt{{\mathit{\lambda }}_l^{\mathrm{2}}+{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}}}}{\sum _{l=\mathrm{1}}^{\mathrm{\infty }}\frac{{\mathit{\beta }}_l}{\sqrt{{\mathit{\lambda }}_l^{\mathrm{2}}+{\mathit{\omega }}_{\mathrm{2}}^{\mathrm{2}}}}}}=C^{st},\end{array}$$

where

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\alpha }}_l=\langle \mathit{\alpha },{\mathit{\varphi }}_l{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}\langle (\mathrm{1}-\mathit{\alpha }(T\left)\right)V^{\prime }\left(T\right),{\mathit{\varphi }}_l{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}},\\ \text{(28)}& {\displaystyle }& {\displaystyle }{\mathit{\gamma }}_l=\langle \mathit{\alpha },{\mathit{\varphi }}_l{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}\langle V^{\prime }\left(T\right),{\mathit{\varphi }}_l{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}\end{array}$$

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\beta }}_l=\langle T,{\mathit{\varphi }}_l{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}\langle (\mathrm{1}-\mathit{\alpha }(T\left)\right)V^{\prime }\left(T\right),{\mathit{\varphi }}_l{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}},\\ \text{(29)}& {\displaystyle }& {\displaystyle }{\mathit{\delta }}_l=\langle T,{\mathit{\varphi }}_l{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}\langle V^{\prime }\left(T\right),{\mathit{\varphi }}_l{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}.\end{array}$$

In the case that the susceptibilities are all dominated by one term with index l, this reduces to $C^{st}=\left({\mathit{\alpha }}_l{\mathit{\delta }}_l\right)/\left({\mathit{\beta }}_l{\mathit{\gamma }}_l\right)=\mathrm{1}$.

In Fig. 5, the sensitivity of temperature to varying amplitude of the albedo function is shown, as well as the sensitivity of the snow-albedo feedback and condition for the existence of an emergent constraint. As shown in Fig. 5a, no emergent relationship is found for climate sensitivity, a feature that was analytically found in the case of constant albedo. In Fig. 5b, the emergent constraint on SAF is shown. In the warm regime (low albedo, lower line in the figure), the SAF becomes larger for larger α_amp. The larger the maximum albedo, the steeper the logistic albedo function. A second effect also takes place, the higher the maximum albedo, the warmer it gets. Consequently, sensitivity of the albedo function is smaller. This decrease in sensitivity also takes place in the cold regime; the colder it gets, the less sensitive the albedo gets. In the cold regime, it is clear that this second mechanism dominates. The results can be reproduced by use of LRT, as shown in Fig. 5c and d. The discrepancies disappear when forcing is small; the climate change forcing in particular is causing most of the differences.

One can extend the energy balance model by representing the response of snow and ice explicitly as a relaxation towards the logistic reference albedo function α(T) given in Eq. (26). This gives the extended model

$$\begin{array}{}\text{(30)}& {\displaystyle }& {\displaystyle }\mathrm{d}T={\displaystyle \frac{\mathrm{1}}{c_T}}\left(Q(\mathrm{1}-\stackrel{\mathrm{‾}}{\mathit{\alpha }})+H-\mathit{ϵ}\mathit{\sigma }{T}^{\mathrm{4}}\right)\mathrm{d}t+\sqrt{{\mathit{\sigma }}_T}\mathrm{d}{W}_t\text{(31)}& {\displaystyle }& {\displaystyle }\mathrm{d}\stackrel{\mathrm{‾}}{\mathit{\alpha }}=-{\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{s}}}}\left(\stackrel{\mathrm{‾}}{\mathit{\alpha }}-\mathit{\alpha }\left(T\right)\right)\mathrm{d}t+\sqrt{{\mathit{\sigma }}_{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}}\mathrm{d}{W}_t,\end{array}$$

where ${\mathit{\tau }}_{\mathrm{s}}=\mathrm{4}\times {\mathrm{10}}^{\mathrm{6}}$ s is the response time of the albedo. The drift term in the Fokker–Planck equation corresponding to Eq. (31) is not the gradient of a potential but the eigensolutions of the generator can of course still be computed numerically.

Figure 5 (a) The relation between temperature response to the seasonal cycle and the temperature response to greenhouse gas forcing. (b) The strength of the snow-albedo feedback to solar and greenhouse gas forcing on different timescales. In the inset: their ratio as a function of α_max. For clearness, panels (a, b) are shown without noise. (c) The susceptibilities for temperature as the observable. (d) The ratio of albedo and temperature susceptibilities and their ratio (RFS).

Extending the model with an explicit albedo function does not change the dynamics of the system significantly, nor the eigenvalues and eigenvectors. Figure 6b shows the eigenvalues of the extended EBM to be almost exactly equal to the eigenvalues of the original model, the imaginary parts continuing to be zero. The projection coefficients are very similar as well (not shown). Thus, the inclusion of a smaller timescale does not improve the response.

4.3 PlaSim

To bridge the gap between parameter ensembles in simple dynamical systems and Earth system models, the SAF emergent constraint is further examined in PlaSim. PlaSim is a numerical model of intermediate complexity, developed at the University of Hamburg to provide a fairly realistic present climate which can still be simulated on a personal computer (Fraedrich et al., 2005). The atmospheric dynamics are modelled using the primitive equations formulated for temperature, vorticity, divergence and surface pressure. Moisture is included by transport of water vapour. The equations are solved using the spectral method. A full set of parameterizations is used for unresolved processes such as long- and shortwave radiation with interactive clouds, boundary layer fluxes of latent and sensible heat and diffusion.

In this climate model, snow albedo is a function of surface temperature T_s, snow depth and vegetation cover. The bare soil snow albedo in PlaSim is described by

$$\begin{array}{}\text{(32)}& {\displaystyle }{\displaystyle }{A}_{\mathrm{snow}}=\left\{\begin{array}{ll}{A}_{\max },& \mathrm{if}{T}_{\mathrm{s}}\le \mathrm{10}{}^{\circ }\mathrm{C}\\ A_{\min },& \mathrm{if}{T}_{\mathrm{s}}>\mathrm{0}{}^{\circ }\mathrm{C}\\ A_{\min }-\left(A_{\max }-A_{\min }\right)\frac{T_{\mathrm{s}}}{\mathrm{10}}& \mathrm{otherwise}\end{array}\right..\end{array}$$

This equation is modified in the presence of vegetation and in the case of shallow snow depth; see Lunkeit et al. (2011) for more details. A set of simulations was performed with A_max varying between 0.650 and 0.900. The historical forcing in PlaSim was approximated by a CO₂ increase from 295 ppm at a rate of 0.3 % per year in the 20th century and 1 % per year in the 21st century before it stabilized at 720 ppm; a 50-year spin-up corresponding to the period 1850–1900 was used.

Figure 6 (a) Eigenvalues of the EBM depending on the amplitude of the albedo function for the simple EBM. The zero eigenvalues correspond to the invariant measure. (b) The extended EBM. (c) Albedo projection terms for solar forcing (α_l) as defined in Eq. (29) where the markers denote l. (d) Same for temperature β_l (e, f) projections terms for greenhouse gas forcing γ_l and δ_l, respectively.

In Fig. 7, the PlaSim results are shown which can be compared to the results from Hall and Qu (2006) in Fig. 1. Note that the variation in CMIP3 is significantly larger than the variation found in PlaSim but that the PlaSim results fit on the relation found by Hall and Qu (2006). Variations in other parameterizations, such as the maximum snow albedo over forested regions, increase the spread in PlaSim SAF further (not shown). This simulation shows that the constraint that emerges in a multi-model ensemble with structurally different formulations of the snow response can to some extent also be reproduced using variations in one parameter. This provides the justification for simplifying further to energy balance models to examine the SAF emergent constraint.

Figure 7 Same as Fig. 1 but now with results from PlaSim.

5 Summary, discussion and conclusions

In this paper, we have presented a dynamical framework behind the occurrence of emergent constraints in parameter-dependent stochastic dynamical systems. In these systems, emergent constraints are related to ratios of response functions which can be determined using linear response theory. It was shown that for a large class of systems, these ratios could be expressed in terms of eigenvalues and projections on eigenvectors of the generator of the system.

A classification of emergent constraints was given and several types could be distinguished depending on whether similar (direct) or different (indirect) observables are considered and whether a response in present-day climate (dynamical) or the time-independent part of present-day climate (static) is linked to a response of the future climate system. For a linear dynamical emergent constraint, the ratio of susceptibilities at the two frequencies under consideration should be a positive constant over the ensemble. When the response is computed with respect to an internal variable (in contrast to an external forcing), a condition is posed on the susceptibilities of the two observables in the system. Static constraints are encountered when a linear relationship is found between the expectation value of the observable and the susceptibility at the frequency of the forcing.

Examples were given using several idealized climate models. In particular, the emergent constraints involving the snow-albedo feedback was considered in detail. We found that linear dynamical emergent relationships can occur when the timescale of the system, indicated by the eigenvalues, changes with the parameter and is smaller than the forcing timescales. This is of particular interest because differences in response size between climate models is often determined by feedback strength in climate systems. Larger feedbacks give rise to larger timescales (Roe, 2009), which is reflected in the eigenvalues of the generator. For an emergent constraint on a feedback quantity, a more complicated constraint mechanism occurs, where one has to take into account the response to two different observables, which typically have different timescales. When the condition of the predictor's timescale being smaller than the forcing timescale is not met, deviations from linearity occur. When the linearity of the relation is exploited in further analysis, such as in the interpretation of emergent constraints by Wenzel et al. (2014), this might lead to a bias in the estimate of the predictand.

Modelling emergent constraints with conceptual models is justified when different Earth system models (ESMs) are closely related and structural differences can be parameterized. This can, for instance, be tested using an intermediate complexity model with full parameterization of the process under consideration.

The classification of emergent constraints provided gives a hint to which kind of emergent constraints one can look out for in an ensemble of high-dimensional global climate models (GCMs). To find an emergent constraint for climate sensitivity by data mining in a CMIP5 ensemble proved fruitless (Caldwell et al., 2014). Using the susceptibilities to find new emergent constraints does not seem to have a direct advantage above directly looking for plausible correlations, but susceptibilities might provide additional information. For example, when a susceptibility shows a resonance at a certain frequency over the ensemble of models, this could suggest that the same feedback is present in all simulations.

In a high-dimensional dynamical system, eigenfunctions and eigenvalues can be accessed with the help of transfer operators, associated with the propagation of probability densities associated with the Fokker–Planck operator. The eigenfunctions that lie on the invariant measure are then computed by making use of the ergodic properties of the climate system. To overcome the burden of high dimensionality, a reduced transfer operator can be computed from a very long simulation, from which the eigenfunctions on the attractor are approximated (Tantet, 2016). However, a forcing on the system does not generally lie only on the attractor and should be split into a part parallel and perpendicular to the attractor. Consequently, the eigenvectors off the attractor cannot a priori be ignored (Lucarini and Sarno, 2011). Gritsun and Lucarini (2017) showed that indeed for some geophysical systems, specifically quasi-geostrophic flow with orographic forcing, the response to the forcing may have no resemblance to the unforced variability in the same range of spatial and temporal scales.

In conclusion, while the current theoretical framework provides an understanding on how emergent constraints may arise in low-dimensional stochastic dynamical systems, its application to output from GCMs, in particular in finding novel and useful emergent constraints, is a challenging issue for future work.

Appendix A

For A=x, we find from Eq. (5) that

$$\begin{array}{}\text{(A1)}& {\displaystyle }{\displaystyle }{R}_{\mathrm{A}}\left(t\right)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}x{e}^{{\mathcal{L}}^{*}t}\left(-{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}{\partial x}}\right)\mathrm{d}x.\end{array}$$

Using the expression for the equilibrium solution ${\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}$ from Eq. (3), we find

$$\begin{array}{}\text{(A2)}& {\displaystyle }{\displaystyle }-{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}{\partial x}}={\displaystyle \frac{\mathrm{2}}{\mathit{\sigma }}}{V}^{\prime }\left(x\right){\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}},\end{array}$$

and hence Eq. (33) becomes

$$\begin{array}{}\text{(A3)}& {\displaystyle }{\displaystyle }{R}_{\mathrm{A}}\left(t\right)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}x{e}^{{\mathcal{L}}^{*}t}\left({\displaystyle \frac{\mathrm{2}}{\mathit{\sigma }}}{V}^{\prime }\left(x\right){\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}\right)\mathrm{d}x.\end{array}$$

With the standard L² inner product, the adjoint of ℒ determined as $\langle {\mathcal{L}}^{*}g,h\rangle =\langle g,\mathcal{L}h\rangle$, where ℒ is the generator of the OU process, is given by

$$\begin{array}{}\text{(A4)}& {\displaystyle }{\displaystyle }\mathcal{L}u=V^{\prime }\left(x\right){\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\mathit{\sigma }}{\mathrm{2}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial x^{\mathrm{2}}}}.\end{array}$$

Using this property in Eq. (35), we find

$$\begin{array}{}\text{(A5)}& {\displaystyle }{\displaystyle }\langle x,e^{{\mathcal{L}}^{*}t}\left(V^{\prime }\left(x\right){\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}\right)\rangle =\langle e^{\mathcal{L}t}x,V^{\prime }\left(x\right){\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}\rangle ,\end{array}$$

and hence

$$\begin{array}{}\text{(A6)}& {\displaystyle }{\displaystyle }{R}_{\mathrm{A}}\left(t\right)={\displaystyle \frac{\mathrm{2}}{\mathit{\sigma }}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{\mathcal{L}t}\left(x\right)V^{\prime }\left(x\right){\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}\mathrm{d}x.\end{array}$$

Next, an inner product $\langle g,h{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}$ is defined as

$$\begin{array}{}\text{(A7)}& {\displaystyle }{\displaystyle }\langle g,h{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}gh{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}\mathrm{d}x.\end{array}$$

As a next step, let λ_l and ϕ_l be the eigenvalues of the generator, i.e. solutions v of

$$\begin{array}{}\text{(A8)}& {\displaystyle }{\displaystyle }\mathcal{L}\mathit{\varphi }=-\mathit{\lambda }\mathit{\varphi }.\end{array}$$

For reversible processes, these eigenvalues are real, positive and discrete under the inner product $\langle ,{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}$. The eigenfunctions form a complete orthonormal basis, such that $\langle {\mathit{\varphi }}_n,{\mathit{\varphi }}_m{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}={\mathit{\delta }}_{nm}$ (Pavliotis, 2014). Now, e^ℒt(x) represents solutions u(x, t) of the problem

$$\begin{array}{}\text{(A9)}& {\displaystyle }{\displaystyle \frac{\partial u}{\partial t}}=\mathcal{L}u,\end{array}$$

with initial condition $u(x,\mathrm{0})=x$. We can expand u into eigenfunctions as

$$\begin{array}{}\text{(A10)}& {\displaystyle }{\displaystyle }u(x,t)=\sum _{l=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\alpha }}_l{\mathit{\varphi }}_l\left(x\right)e^{-{\mathit{\lambda }}_{l}t}.\end{array}$$

From the initial condition, we find

$$\begin{array}{}\text{(A11)}& {\displaystyle }{\displaystyle }\sum _{l=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\alpha }}_l{\mathit{\varphi }}_l\left(x\right)=x,\end{array}$$

and using the orthogonality of the ϕ_l under the inner product $\langle ,{\rangle }_{{\overline{p}}_{\mathrm{e}}}$, we find

$$\begin{array}{}\text{(A12)}& {\displaystyle }{\displaystyle }{\mathit{\alpha }}_l=\langle x,{\mathit{\varphi }}_l{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}.\end{array}$$

On the other hand, substituting the expression for u into Eq. (38) gives

$$\begin{array}{}\text{(A13)}& {\displaystyle }{\displaystyle }\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sum _{l=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\alpha }}_l{\mathit{\varphi }}_l\left(x\right)e^{-{\mathit{\lambda }}_{l}t}{V}^{\prime }\left(x\right)\stackrel{\mathrm{‾}}{p}\mathrm{ed}x=\sum _{l=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\beta }}_l{e}^{-{\mathit{\lambda }}_{l}t},\end{array}$$

where

$$\begin{array}{}\text{(A14)}& {\displaystyle }{\displaystyle }{\mathit{\beta }}_l={\mathit{\alpha }}_l\langle V^{\prime }\left(x\right),{\mathit{\varphi }}_l{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}=\langle x,{\mathit{\varphi }}_l{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}\langle V^{\prime }\left(x\right),{\mathit{\varphi }}_l{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}.\end{array}$$

Repeating the derivation with a general observable A=f(x) gives $\langle f\left(\mathit{x}\right),{\mathit{\varphi }}_l{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}\langle V^{\prime }\left(\mathit{x}\right),{\mathit{\varphi }}_l{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}$. The first term in β_l denotes the projection of the observable on the eigenfunctions and could intuitively be interpreted (for l>0) as the amenability of the observable to change. The second projection term in β_l can be understood to be the amenability of the whole system to change under the influence of the forcing field. In Eq. (12), those observables are Y₁ and Y₂, so that $g_l=\langle Y_{\mathrm{1}},{\mathit{\varphi }}_l{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}\langle V^{\prime }\left(\mathit{x}\right),{\mathit{\varphi }}_l{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}$ and $h_l=\langle Y_{\mathrm{2}},{\mathit{\varphi }}_l{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}\langle V^{\prime }\left(\mathit{x}\right),{\mathit{\varphi }}_l{\rangle }_{{\stackrel{\mathrm{‾}}{p}}_{\mathrm{e}}}$.