A method is proposed for combining information from several emergent constraints into a probabilistic estimate for a climate sensitivity proxy Y such as equilibrium climate sensitivity (ECS). The method is based on fitting a multivariate Gaussian PDF for Y and the emergent constraints using an ensemble of global climate models (GCMs); it can be viewed as a form of multiple linear regression of Y on the constraints. The method accounts for uncertainties in sampling this multidimensional PDF with a small number of models, for observational uncertainties in the constraints, and for overconfidence about the correlation of the constraints with the climate sensitivity. Its general form (Method C) accounts for correlations between the constraints. Method C becomes less robust when some constraints are too strongly related to each other; this can be mitigated using regularization approaches such as ridge regression. An illuminating special case, Method U, neglects any correlations between constraints except through their mutual relationship to the climate proxy; it is more robust to small GCM sample size and is appealingly interpretable. These methods are applied to ECS and the climate feedback parameter using a previously published set of 11 possible emergent constraints derived from climate models in the Coupled Model Intercomparison Project (CMIP). The ±2σ posterior range of ECS for Method C with no overconfidence adjustment is 4.3 ± 0.7 K. For Method U with a large overconfidence adjustment, it is 4.0 ± 1.3 K. This study adds confidence to past findings that most constraints predict higher climate sensitivity than the CMIP mean.
Climate change is a defining problem of our time. It is hard to plan for future warming without knowing its magnitude, but our ±1σ “likely” confidence range for equilibrium climate sensitivity (ECS; the global-average surface warming due to doubling CO2 and letting the climate re-equilibrate) is currently 1.5–4.5 K (IPCC 2013)—which is disturbingly large. This uncertainty has persisted for decades despite large advances in our understanding of the climate system (Knutti et al. 2017).
Emergent constraints offer a possible path to narrowing this spread. An emergent constraint is a current-climate quantity that has skill at predicting future changes in climate. Such predictors may be valuable shortcuts to the complex and uncertain process of directly simulating climate change in a general circulation model (GCM) or inferring it from imperfect observational records. Because the physical processes governing climate change are generally the same ones that control present-day seasonal, weather-scale, and diurnal variations, it is likely that real emergent constraints exist. Hall and Qu (2006) was one of the first papers to identify such a constraint. They found that the seasonal cycle of snow albedo over Northern Hemisphere land is tightly correlated with snow albedo feedback over this region in 17 model simulations from phase 3 of the Coupled Model Intercomparison Project (CMIP3). This emergent constraint has an obvious motivation: surface warming reduces snow cover irrespective of whether that warming is due to seasonal changes in insolation or CO2-induced climate change. Nearly 40 other emergent constraints have been proposed since 2006 (Hall et al. 2019), although few have had such a satisfying physical explanation.
Several limitations and assumptions apply to the use of emergent constraints to predict ECS. First, the ECS simulated by a climate model is generally estimated from an integration of finite length (customarily 150 years in CMIP5) that is not fully equilibrated and keeps certain physics fixed (e.g., vegetation type and land ice). The response of that model (or the real climate system) to a time-varying radiative forcing is not determined purely by the ECS, but also depends upon natural and forced changes in the pattern of surface warming (e.g., Armour et al. 2013; Gregory and Andrews 2016). ECS is also a problematic target for emergent constraints because it arises from interaction between many processes. As a result, it is questionable whether any single current-climate variable would explain a large fraction of ECS variability. This is why Klein and Hall (2015) suggest that emergent constraints should be targeted toward a single climate feedback mechanism (e.g., snow cover) whenever possible. Nevertheless, many studies (including this one) focus on emergent constraints for global climate sensitivity proxies such as ECS or the climate feedback parameter λ (Cess et al. 1989) because of their importance. Last, emergent constraints in general derive from a blend of scientific reasoning and a posteriori optimization to maximize their correlation with ECS or λ over a modest set of GCMs, and only the most promising constraints are likely to be published. This suggests a risk of constraints being “overconfident”—that is, better correlated with ECS or λ over the GCMs on which they were first tested and optimized than in another independent set of GCMs. A statistical analysis by Caldwell et al. (2014) showed that a posteriori data mining of the limited set of GCMs in the CMIP5 archive could easily yield apparent constraints with a strong sample correlation to ECS despite there being no true correlation.
Emergent constraints have already been noted to predict larger climate sensitivity than expected from other lines of evidence (Tian 2015; Klein and Hall 2015). If agreement between constraints gives us confidence in their predictions, this is an alarming finding. A goal of this paper is to develop an approach for combining emergent constraints to provide a confidence range for ECS, or any other climate sensitivity proxy, while accounting for the issues just raised.
A key idealization that we make is that the climate proxy and the constraints have a multivariate Gaussian distribution. This idealization is not perfect. Both ECS and several of the emergent constraints that we use marginally fail statistical tests of normality, as we will describe in section 2b. This should make us wary about the tails of any estimated probability density function (PDF), which are most sensitive to the Gaussian assumption. We will use our theory to estimate Gaussian PDFs for ECS and also for its approximate inverse λ, the climate sensitivity parameter. Those PDFs cannot be fully consistent with each other because if λ is Gaussian then ECS is not. This is most problematic if these PDFs encompass a broad range of climate sensitivities, which ours do not. We have chosen to highlight ECS because it has been a popular target for emergent constraints. It might be more plausible to assume that λ, which can be composed as the sum of physical feedbacks, has a Gaussian distribution; that will yield a PDF of ECS with a longer positive tail. Readers should be aware of the inevitable trade-off between the simplicity of our approach and the limitations encompassed by the assumptions that go into it.
In mathematical terms, our approach is a straightforward application of the theory of multivariate Gaussian distributions developed by statisticians many decades ago. This theory may be unfamiliar to most climate scientists, but our method ultimately reduces to an application of multiple linear regression.
The emergent constraints used and the relevant data are described in section 2, and terminology is described in section 3. In section 4, we derive and apply our method to individual emergent constraints with observational uncertainty. Section 5 discusses correlations between the constraints. Section 6 uses Method C (described therein) to derive a Gaussian PDF of the climate proxy given multiple correlated constraints and shows its connection to multiple linear regression. It then applies Method C to sets of 4 and 11 constraints, including the use of robust regression to reduce artifacts in the 11-constraint case arising from strong sample correlations between constraints that must be estimated using an overly small set of GCMs. Section 7 presents and applies Method U, a special case of Method C that neglects any correlations between constraints outside their mutual correlation with the climate proxy. Method U can be more easily analyzed and interpreted, and it is used to explain why the full set of constraints confidently predicts a higher climate sensitivity than almost any single constraint. It is also more robust to our small sample of GCMs. Section 8 presents and applies an overconfidence adjustment to account for overfitting. Section 9 presents conclusions.
a. Choice of emergent constraints
For this study, we rely on 11 emergent constraints evaluated in Caldwell et al. (2018, hereinafter CZK18). These include the four constraints CZK18 judged to be “credible” (significantly correlated with ECS and supported by a physical mechanism that correctly identifies dominant physical processes and geographical regions that create this correlation), and seven constraints they judged to be “uncertain” or “unclear” (significantly correlated with ECS but not amenable to the above assessment of credibility). We will call these constraints “possible.” Two other “unclear” Klein constraints from CZK18 had to be excluded from our analysis for technical reasons described in section 2b. We also excluded six constraints assessed not to be credible in Table 4 of CZK18. Short explanations of each constraint that we used along with original citations and evaluations from CZK18 are provided in Table 1.
b. Do the constraints have Gaussian distributions across the GCMs?
With our small subset of GCMs, we cannot fully test the plausibility of a multivariate Gaussian distribution for describing the joint PDF of a climate proxy such as ECS and a large set of emergent constraints. However, we can test the Gaussianity of the univariate marginal distributions of ECS and the individual constraints using a test of Shapiro and Wilk (1965). The ECS and 4 of the 11 constraints that we will use (1, 2, 8, and 11) fail this test at the 95% significance level (p < 0.05), but none of these fail at the 99% significance level. Since the GCMs are not fully independent, these significance levels may be overstated. Overall, we conclude that the Gaussianity assumption is marginally acceptable for our dataset.
c. Observational estimates of constraints
CZK18 focused on the evaluation of emergent constraints using model data, while this study aims to use the observed values of those constraints to make climate sensitivity predictions. This requires observational estimates (including uncertainty) for the constraints. It would be ideal to obtain these estimates directly from the original data sources, but this is impractical given the number and diversity of constraints we use. Thus, we rely almost exclusively on values communicated by the papers originally proposing each constraint. These studies employed a variety of approaches and levels of detail in describing observational uncertainty. As a result, we are forced to make approximations to achieve uniformity of observed uncertainty estimates across constraints.
For simplicity, our analysis assumes observational uncertainty is normally distributed. The PDF of observed values for constraint i is specified by its mean μo,i and standard deviation σo,i. Although convenient, this assumption is not appropriate for two “Klein” constraints discussed by CZK18, which are based on positive semidefinite measures of model skill. Hence these two constraints were excluded from our analysis.
Our observed values and the information used to construct them are summarized in Table 2. Studies that provide mean and some multiple of the standard deviation were trivial to process. For studies which provide bounds for a given confidence level, we compute the number of standard deviations for that confidence level for a normal distribution, and we rescale the quoted range to estimate σo,i. Where several estimates of μo,j and σo,j were provided, we average the estimated means, and we increase σo,j such that μo,i ± 1σo,i just encompasses all of the individual estimated ±1σ ranges. For constraints that provide only minimum and maximum credible values (often taken from a pair of observations), we take the average of these values as the mean and 1/2 the distance between these values as the standard deviation. Because two samples provide a very poor sense of uncertainty, we occasionally use extra information from papers and/or expert judgement to modify these values, as noted in Table 2.
3. Terminology and covariance estimation
Our mathematical nomenclature is as follows. Capitalized Latin letters denote random variables, and lowercase versions of the same letter indicate particular values of these variables. Vectors are boldfaced. We define to be a climate sensitivity proxy such as the equilibrium climate sensitivity or a climate feedback strength, for which the constraints are derived. A single emergent constraint variable is denoted . A collection of n emergent constraints will be labeled . Versions of these random variables that have been normalized to have zero mean and variance of 1 are similarly denoted, but without the tilde. The PDF of any random variable U is p(u), and similarly for multivariate distributions.
The main mathematical formula that we use is the joint PDF of the components of a column vector U of m zero-mean Gaussian random variables that are known to have an m × m covariance matrix with determinant :
Calculating correlation and covariance from GCM samples
A key input to our analysis is the (n + 1) × (n + 1) covariance matrix between the climate proxy and the n constraints , derived from the available sample of GCMs. An important complication is that not all GCMs provide the data needed to compute all constraints. Using all available GCMs for calculating each needed covariance, rather than just the eight models that supplied data for all 11 constraints, is essential to obtaining an adequate sample size. The method used to do this must preserve the positive definiteness of the covariance matrix for the multivariate Gaussian method to provide stable results.
The approach that we settled on is to build a GCM covariance matrix based on the best possible estimates of the correlation coefficients. We compute the correlation coefficient between each pair (i, j) of constraints using all GCMs for which both constraints are available. We use a similar approach for GCM-based correlation coefficients between the climate proxy and the jth constraint, as well as for calculating the standard deviation of each constraint and of the climate proxy across all GCMs. The elements of the covariance matrix are computed as
Note that in general, a different set of GCMs is used for computing each of the three terms on the right-hand side.
Here and in the rest of the paper, rows and columns of the covariance matrix are indexed starting at 0, index 0 corresponds to the climate proxy, and indices 1 through n correspond to the n constraints.
4. Climate sensitivity PDF from a single constraint
In this section we describe our approach for computing a PDF of the climate sensitivity proxy from a single constraint . This is a useful first step toward treating multiple constraints. Like Bowman et al. (2018) and others, we first estimate a joint PDF of and from the GCMs, and then we apply our observational knowledge about to derive a constrained PDF of . A key assumption, questioned by Williamson and Sansom (2019), is that the GCM-derived joint PDF is applicable to the real climate (i.e., is a suitable prior for interpreting an observation of the constraint). Also like Bowman et al. (2018), we make the important simplifying assumption that the multivariate PDFs that we estimate are Gaussian, allowing them to be described in terms of a vector of means and a covariance matrix. See Cox et al. (2018), Brient and Schneider (2016), Bowman et al. (2018), and Williamson and Sansom (2019) for other approaches to deriving a PDF of a climate sensitivity proxy from a single emergent constraint, and for further discussions about issues with applying emergent constraints. Our approach captures all sources of uncertainty without sacrificing simplicity, and it is easily extensible to multiple constraints.
An emergent constraint is based on a GCM-based relationship between and . Such a relationship should not be trusted far outside the range of GCM values. Indeed, if the observed value of the constraint lies well outside the expected GCM range, we might interpret this as a physical shortcoming of the GCMs that requires further attention, rather than a solid basis for inferring that the climate sensitivity proxy lies outside its GCM range.
Philosophically, this frames our mathematical representation of emergent constraints. Unlike prior studies, we do not start by performing a GCM-based linear regression to determine from an observationally constrained . Instead we estimate a joint Gaussian PDF between and by substituting their 2 × 2 sample covariance matrix into (1). In contrast to linear regression, this retains information about the GCM-preferred range of . It tacitly assumes that relations between the climate proxy and the constraints are nearly linear. It also assumes that all GCMs have equal value in estimating how the emergent constraint is related to the climate proxy, whether or not they predict realistic values of the proxy. That assumption has been reasonably criticized (e.g., Brient 2020) but it is a fundamental premise of emergent constraints that the underlying relationship with the climate proxy should rely on a mechanism sufficiently robust as to be insensitive to details of the GCM physical formulation, even though those details are important to actually obtaining an observationally consistent value of the constraint.
If we could exactly observe that , we could substitute this observed value into the joint PDF to obtain the conditional PDF of . In the language of Bayesian analysis, this is the posterior probability for based on the GCM-only Gaussian prior and the observation of the constraint. However, there are two further practical complications to consider.
First, we cannot exactly observe the true value of . To handle observational uncertainty, we define a random variable , the estimated constraint, which is the sum of plus a normally distributed observational error with zero mean and the observed standard deviation σo. This observational error is assumed to be independent of and the physically determined (true) value of . Thus is a Gaussian random variable whose variance is the sum of its variance across GCMs and its observational variance. Furthermore, and have a Gaussian joint PDF determined by their covariance matrix,
where is computed using Eq. (2).
Second, the process of formulation and selection of emergent constraints may result in overconfidence, that is, constraints that are more highly correlated with the climate sensitivity proxy than would be obtained from a different independent random sample of GCMs, if such existed. We counteract overconfidence by artificially reducing the covariance between X and without changing the other elements of . We have not explicitly included this adjustment in our single-constraint analysis because a goal of that analysis is to evaluate overlap between PDFs for each constraint, and ignoring overconfidence provides a lower bound for that overlap. We consider the sensitivity of a multiple-constraint analysis to an overconfidence correction in section 8.
It is convenient to work with standardized variables with a mean of 0 and a standard deviation of 1:
The normalization of the standardized X accounts for both its variance across GCMs and its observational uncertainty. Overlines indicate averages over the GCM ensemble.
The covariance matrix of Y and X, derived from the sample of GCMs and adjusted for observational constraint uncertainty, is
where is the correlation coefficient between Y and X, or equivalently between and . Because of the observational uncertainty, r is smaller in magnitude than the GCM-estimated correlation coefficient between the climate proxy and the constraint .
We translate the observational estimate of the constraint, , into the standardized form
We condition the joint PDF of Y and X on this known value of X to obtain a Gaussian posterior for Y:
where, substituting Eq. (7) for , the logarithmic probability is
From this formula, we can read off the mean y(1) and standard deviation σ(1) of the Gaussian posterior distribution of Y:
The superscript (1) denotes that this is a one-constraint estimate. We will define the estimated posterior “range” of Y as lying within 2 standard deviations of the mean, that is, y(1) ± 2σ(1). Later, we will use a subscript i to denote an estimate that is based on constraint i.
Past studies of emergent constraints have typically used linear regression to quantify the relationship between the constraint and the climate sensitivity proxy. If we ignore observational uncertainty we could obtain the above result by regressing the climate sensitivity proxy Y on the constraint X. In this case, r would be the GCM-based correlation coefficient between and with no adjustment for the observational uncertainty. The best fit regression line y = rx matches the posterior mean y(1) when evaluated at xo. The residual in Y around that fit has a standard error (1 − r2)1/2 that matches σ(1).
Our approach naturally generalizes this regression method. It incorporates observational uncertainty in the constraint. For a single constraint, this was also done by Bowman et al. (2018) under similar assumptions but using a different mathematical approach. It is reassuring that, after accounting for our different notation and normalization, our formulas in Eqs. (10) and (11) are isomorphic to Eqs. (18) and (23) of Bowman et al. (2018). Unlike earlier work our approach extends naturally to many constraints, for which it becomes a form of multiple linear regression.
Single-constraint results for ECS
Table 3 gives the correlation coefficients ri between each constraint i and two choices of climate proxy Y (ECS and climate feedback parameter λ). To simplify the ensuing discussion, we henceforth “sign correct” all constraints so that ri > 0, by flipping the sign of those constraints that are negatively correlated with Y (indicated by a minus sign in the “Sign” column). For ECS, we also tabulate the GCM-based correlation coefficient before accounting for observational uncertainty. It will always be larger than ri since the observational uncertainty is assumed to be uncorrelated with ECS. This table also gives normalized constraint values xo,i (in units of standard deviation) for the sign-corrected constraints, such that xo,i > 0 favors y > 0 (climate proxy larger than the GCM mean). Eight of the 11 constraints have positive xo,i, with values up to 2.4 for constraint 1 (Sherwood D). Constraints 2 (Brient Shal), 9 (Lipat), and 11 (Cox) have modestly negative xo,i in the range from −0.4 to −0.8.
Last, Table 3 gives , the ratio of the observational uncertainty to the GCM-based standard deviation for each constraint [this is the inverse square root of the signal-to-noise ratio defined in Eq. (19) of Bowman et al. (2018)]. For the credible constraints 1–4 and constraint 10 (Siler), this ratio is less than 0.5 and observational uncertainty is relatively unimportant to the posterior range of Y. For the remaining 6 constraints, the ratio is larger than 0.5. For these constraints, observational uncertainty cannot be neglected—it substantially reduces ri, broadens the posterior range, and moves toward zero (i.e., it moves the posterior mean Y toward the GCM mean). This is most pronounced for Constraint 6 (Qu), with a ratio of 1.7; the ratio lies between 0.64 and 1.01 for the remaining five constraints.
The formulas in the previous section were nondimensional. To redimensionalize a constraint value x, we scale with its standard deviation (derived as the quotient of σo from Table 2 and from Table 3) and add its mean (from Table 2). We redimensionalize ECS by scaling with the GCM standard deviation K and adding the GCM mean K, which can be read off the “GCM Ensemble” line of Table 4; similarly for the climate sensitivity parameter λ.
Our approach is illustrated in Fig. 1, using the Su constraint as an example. For clarity, the figure is presented in terms of original rather than standardized variables. This constraint is constructed by computing the climatological-average latitude–height cross section of zonal average relative humidity (a proxy for the Hadley circulation) from model and observations, then computing for each model the slope of model versus observed values using each (latitude, height) value as a separate data point. Perfect agreement with observations occurs when the (nondimensional) slope is 1. Observations are inferred from satellite remote sensing; their uncertainties impart to the constraint a relatively large observational standard error of 0.25, as indicated by the PDF along the y axis. The cyan ellipse shows a contour of joint probability density between ECS and the Su constraint derived from our collection of CMIP GCMs (black dots). The maximum joint density is at the center of this ellipse, which is the centroid of the CMIP data points. The green ellipse accounts for the substantial observational uncertainty in the Su constraint, which vertically broadens the PDF. The posterior for ECS (red) is a cross section of the bivariate PDF at the best-guess observed constraint value of 1 (horizontal blue dashed line).
The posterior PDFs of ECS given each constraint separately are shown in Fig. 2, along with the PDF of ECS from CMIP3 + CMIP5 models, shown both as a histogram and a Gaussian fit. Credible constraints are shown in Fig. 2a and possible constraints are shown in Fig. 2b. Their means vary from 3.0 K (constraint 9 = Lipat) to 4.0 K (constraint 1 = Sherwood D), with a standard deviation of 0.5 K (constraint 3 = Zhai) to 0.7 K (constraint 6 = Qu). Eight of the 11 constraints have PDFs peaked at an ECS greater than the GCM mean. Constraint credibility has no systematic effect on peak probabilities or distribution widths.
These Gaussian PDFs are derived from the formulas for nondimensional posterior mean [Eq. (10)] and standard deviation [Eq. (11)] using the information in Table 3. These are redimensionalized as discussed above.
These PDFs do not account for possible overconfidence, which would further broaden their width. Even so, all of the constraints have substantial overlap with each other. This is reassuring because if two constraints had disjoint PDFs we would be forced to conclude that at least one of them must be wrong (which would be an disturbing finding worthy of further study). Results for the climate sensitivity parameter λ (not shown) are similar.
5. Dependence between constraints
Combining constraints is only useful if those constraints provide some independent information. Interdependence of constraints was investigated in CZK18 by computing correlations between each pair of constraints. CZK18 found that constraints were more correlated than expected by chance but that identifying pairs of significantly correlated constraints on the basis of related physical explanations was generally unsuccessful. They noted that some pairwise correlation is expected because all constraints are by construction correlated with ECS.
For our 11 constraints, there are 55 pairs of constraints. If we flip the sign of constraints that are negatively correlated with ECS, then 52 of the 55 constraint pairs are positively correlated across the GCMs, which is strong evidence for such mutual dependence; 14 of these have a positive correlation coefficient exceeding 0.4, and none have a correlation coefficient more negative than −0.25.
To correct for the mutual dependence between constraints that is associated with ECS, we compute the partial correlation coefficient (https://www.encyclopediaofmath.org/index.php/Partial_correlation_coefficient) between constraints Xi and Xj given ECS Y (denoted hereinafter by a subscript 0):
This is the correlation coefficient between the residuals after Xi and Xj have been linearly regressed on Y. Here rij is the sample correlation coefficient between Xi and Y, and similarly for ri, while rij is the sample correlation coefficient between Xi and Xj. We call two constraints with a partial correlation of zero “conditionally uncorrelated”; for the multivariate Gaussian distributions assumed in this paper, this is equivalent to conditional independence. The partial correlation coefficient is the correlation coefficient can vary between −1 and 1. Positive sign-corrected partial correlation indicates constraints that covary in the same sense as we would expect based on their correlation with Y.
Correcting for the mutual dependence with ECS removes some but not all of the covariation between our 11 constraints. Figure 3 shows the sign-corrected pairwise partial correlations. Well over one-half (38 of 55) are positive. This is suggestive, but are these partial correlations statistically significant? As in CZK18, choosing an appropriate number of degrees of freedom is difficult because there are complicated structural dependences between models (Masson and Knutti 2011; Knutti et al. 2013; Sanderson et al. 2015). Recall that some constraints could not be computed using outputs available from all GCMs, so each rij0 was computed on the basis of an (i, j)-dependent subset of the GCMs. Following CZK18, we handle this issue by using a fairly lax 90% two-sided test and by assuming each GCM that goes into the calculation of a particular partial correlation coefficient rij0 is independent. Both assumptions favor false positives, but partial correlations deemed not to be significant would almost certainly also be deemed insignificant with other reasonable assumptions. For a typical number of contributing GCMs (20), a partial correlation of magnitude 0.38 or larger is significant by this standard, based on a t test.
We found that only 6 of 55 constraint pairs have a significantly positive partial correlation; these are shown by orange shading in Fig. 3. Conversely, 1 of 55 constraint pairs has a significantly negative partial correlation, shown by purple shading in Fig. 3. Some expected relationships between constraints are borne out, like tight correlation between Siler and Volodin. Other constraints are significantly correlated even though their explanations seem to be unrelated, like Brient Shal and Brient Alb. This motivates our Method C (for “correlated”), which accounts for partial correlation between constraints. However, other constraint pairs, like Zhai and Brient Alb, are weakly correlated even though they share a physical explanation, and only 11% of the constraint pairs have a partial correlation above the 90% significance threshold. Thus, it is also a reasonable overall assumption to neglect partial correlations among our set of constraints, motivating the simpler Method U (for “uncorrelated”) discussed in section 7.
6. Method C: Climate sensitivity PDF from multiple correlated constraints
The method introduced in section 4 generalizes transparently to the case of multiple constraints with observational estimates having uncertainties σ0,1, …, σ0,n. We form a column vector U of length n + 1 whose components are the climate sensitivity proxy and the s. The components of its (n + 1) × (n + 1) covariance matrix modified for observational uncertainty are
Here δij is 1 if i = j and 0 otherwise. The matrix is guaranteed to be positive semidefinite if all covariances are computed with the same set of GCMs, but this is not guaranteed if different elements of are computed with different subsets of GCMs, as we are forced to do in this study.
We standardize and the constraints as before. The covariance matrix between the standardized ECS Y and the standardized constraints Xi is a correlation matrix with components
We also derive standardized values xo,i for each observational estimate .
We need the conditional distribution of Y given known values xo,i for all the random variables Xi. We will show that this is actually an application of multiple linear regression of the climate proxy Y on the constraints Xi, which is in turn a natural generalization of the one-constraint case.
Eaton (1983) provides convenient formulas for the means of the conditional distribution of a multivariate Gaussian distribution given known values of some of the random variables, which we specialize to the case at hand. To use Eaton’s formulas, we must break into blocks:
Here, CYY = 1 is 1 × 1; CYX = r is the column vector of ris, i = 1, …, n; and is the n × n covariance matrix of the standardized constraints Xi, which is with its first row and column removed. We must also know the (unconditional) mean and variance of Y. Since Y is standardized, these are μY = 0 and var(Y) = 1. Last, we need the mean μX = 0 of the column n-vector of standardized constraints.
Then Eaton’s formula for the mean of Y, given the observed estimates xo,i of the Xis, is
where, since is symmetric,
The corresponding formula for the conditional variance of Y is
These formulas show that Method C is equivalent to multiple linear regression of Y on the Xis using the GCMs [with no constant term in Eq. (17), since the predictors and predictand have mean of 0 by construction]. In particular, Eq. (18) is precisely the “normal equation” of multiple regression. Equation (20) is the expected mean squared residual of this prediction. A subtlety is that, unlike in standard multiple linear regression, is not just the empirical covariance matrix of all of the constraints sampled across all of the GCMs, because of our needs to handle missing models and to correct for observational uncertainty of the constraints. Thus, we must construct and solve the normal equations rather than using standard statistical analysis software.
b. Sampling uncertainty in the covariance matrix
In theory, Method C solves our problem for multiple correlated constraints, because it specifies the Gaussian posterior PDF of the climate proxy Y in terms of known quantities. In practice, if the covariance matrix is derived from a finite sample of GCMs it is sensitive to sampling uncertainty, especially if the number of constraints is comparable to the number of GCMs. For instance, if we have 10 constraints and 20 GCMs, we are using 20 samples to estimate an 11-dimensional correlation matrix with 11 × 12 ÷ 2 = 66 independent entries, which is a highly underconstrained problem. Thus, we anticipate that Method C may fail or give spurious results, especially if partial correlations between the constraints are important (multicollinearity). In section 6c, we apply Method C to ECS estimation, and we consider a penalized regression approach to improve its robustness.
c. Applying method C to ECS
1) Four credible constraints
The red curve in Fig. 4a shows the posterior PDF for ECS estimated using Method C based on the four sign-corrected credible constraints. It is compared with other PDFs, including the Gaussian CMIP prior (black dashed) and the average of the single-constraint PDFs (solid black). Methods U and U3 will be discussed in the following section.
Method C gives a strikingly higher and narrower PDF than the average of the single-constraint PDFs, which in turn is shifted slightly higher than the CMIP mean PDF. This may seem surprising but can be viewed as the result of strong evidence (3 of 4 constraints) agreeing that ECS is likely more positive than the GCM mean. Mathematically stated, of the four observed xo,i, only xo,2 < 0.
The nondimensional ECS standard deviation, derived from Eq. (20), is and does not depend on the constraint values. Redimensionalizing using the mean (3.24 K) and standard deviation (0.75 K) of the GCM ensemble yields a 4-constraint ±2σ ECS range of 4.14 ± 0.82 K.
Constraint 3 (Zhai) has the strongest weight a3 = 0.48 on the Method C ECS mean. Constraint 4 (Brient Alb) has only 60% of the weight of constraint 3 despite having a similar correlation r ≈ 0.7 with ECS; this is due to its covariances with the other constraints. All four constraints have positive weights ai > 0. This positivity property is physically reasonable and desirable. For instance, if constraint i has observed value xo,i > 0, its correlation with ECS is ri > 0, and all other constraints have observed values of zero, then we expect y > 0; this is only possible if ai > 0. However, weight positivity is not guaranteed by Method C.
Figure 4b assesses the sensitivity to the GCM sampling by redoing the 4-constraint analysis with each of the GCMs omitted in turn when computing the needed correlation matrix. The posterior ECS PDF is robust to this test. Removing any two GCMs provides very similar results, with the most likely value of ECS always being 3.9–4.2 K (not shown).
2) All constraints
Method C can also be applied to all 11 constraints, since the 12 × 12 estimated correlation matrix has all positive eigenvalues. It gives a ± 2σ ECS range 4.28 ± 0.68 K, shown as the dashed red curve in Fig. 4c. This PDF has a slightly higher mean and a narrower distribution of ECS than the 4-constraint estimate. Eight of the 11 weights ai are positive. Two are marginally negative, and one, a10 = −0.47, is strongly negative. This arises because of multicollinearity; constraint 10 (Siler) has very strong partial correlations with some of the other constraints, exceeding 0.5 with constraints 2 (Brient Shal) and 5 (Volodin). These partial correlations are also quite uncertain due to the small number of GCMs. Because the observed value xo,10 = 0.35 is small, the large negative weight does not have a strong direct impact on the most probable ECS. Nevertheless, this suggests the Method C ECS PDF may be less robust to sampling uncertainty in the correlation matrix with 11 constraints than with 4 constraints. This defeats the point of using more constraints.
Figure 4d assesses this by recomputing posterior PDFs after removing a single GCM, as in Fig. 4b. The posterior ECS PDFs are mostly robust to this test, but less so than with the four constraints, since there are now three clear outliers among the 11 ranges. The mean ECS ranges from 3.7 to 4.6 K across these cases and from 3.6 to 5 K when two GCMs are removed (not shown).
3) Ridge regression to improve robustness
Ridge regression (Hoerl 1962) is a technique for penalizing large constraint weights ai (specifically, the sum of squares of the weights) that can result from multicollinearity and a poorly conditioned covariance matrix. Because samples from all GCMs are not available for every predictor (constraint) and predictand (climate proxy), we can only use regression implementations based on a modified covariance matrix. This precludes most popular robust regression methods, e.g., “LASSO” (Santosa and Symes 1986), which are not trivial to adapt to this problem formulation. However, ridge regression can be cast as multiple regression with an additional diagonal term in the covariance matrix, that is,
A larger value of the user-chosen regularization parameter γ > 0 corresponds to stronger weight penalization. This comes at the expense of potential bias in the estimated predictand toward its a priori mean, in this case the GCM mean of the climate proxy.
For the 11-constraint problem, all the sign-corrected constraint weights became nonnegative for γ ≥ γc = 0.64. This is strong regularization; the sum of the squared weights is less than 20% as large as the unregularized solution. The zero weight at γc corresponds to constraint 10 (Siler), which had a strongly negative sign-corrected weight in the unregularized regression. The PDF for γc is shown as the solid red line in Fig. 4c and is labeled “C, RidgeRegr.” It has a smaller mean and broader width than the unregularized Method C but is qualitatively similar.
The condition number of the unregularized is 34. A value over 30 implies a substantial risk that the results will be sensitive to small uncertainties in the covariance matrix (see https://en.wikipedia.org/wiki/Multicollinearity and references cited therein). Ridge regression with γ = γc desirably reduces this down to 6. For the 4-constraint problem, the condition number of the unregularized is only 7, suggesting no need for regularization.
Ultimately, the choice of γ is subjective. A smaller value γ ≈ 0.25 reduces the ratio of the smallest negative weight to the largest positive weight by a factor of 3 relative to the unregularized solution, with much less reduction in the mean and increase in the width of the PDF relative to γc. This might provide a better trade-off between robustness and bias.
7. Method U: Estimation of Y from conditionally uncorrelated constraints
Overall, Fig. 3 suggests that most of the partial correlations between the constraint pairs are not that large and are statistically insignificant. Thus it is reasonable to consider the special case that all the constraints are mutually uncorrelated when conditioned on a given value of the climate sensitivity proxy Y. In that case, Method C reduces to a simpler Method U (for “uncorrelated”) that requires only estimates of the correlation coefficients ri between Y and the individual constraints (given in Table 3). These can be computed fairly reliably with 10 or more GCMs, eliminating the need for regularization approaches like ridge regression, even when there are many constraints. For a single constraint, Method U and Method C both reduce to the method already presented in section 4.
The Method U mean and variance of the climate proxy Y can be derived from the Method C Eqs. (17)–(20) using linear algebra. However, the same results are easier and more illuminating to derive directly from the joint multivariate Gaussian PDF
Since the constraints are assumed to have Gaussian PDFs, and they are uncorrelated when conditioned on the Gaussian variable Y, they are independent when conditioned on Y. Hence the joint PDF of the standardized constraints conditioned on Y can be written as a product:
Setting the constraints equal to their standardized observed values, we obtain the posterior PDF
where, with a slight rewrite of Eq. (25),
That is, the desired posterior is a product of Gaussian PDFs with means yo,i and standard deviation . For constraints with high correlation ri with the climate proxy, the weight wi is large and the standard deviation is small. The standardized GCM prior p(y) has been included in this notation as an i = 0 Gaussian PDF with mean 0 and standard deviation 1.
Equations (26) and (27) are the keys to understanding how multiple constraints combine to determine the posterior PDF of the climate proxy. The posterior will maximize for those y that are most consistent with all of the yo,is within their individual uncertainties .
For each constraint, the “proxy estimate” yo,i is the predicted value of the climate proxy we would obtain by linear regression applied to the single-constraint problem with the dependent and independent variables swapped, that is, regressing Xi on Y and using the linear fit to determine Y for the observed value Xi = xo,i. In that regression problem, since both Xi and Y are standard normal, the slope of the regression line is their correlation coefficient ri and the residual standard deviation of Xi is . Hence both the mean and standard error of the corresponding estimate of Y from this constraint are a factor 1/ri larger than for Xi. For constraints that are poorly correlated with Y, the proxy estimate yo,i can be quite large even if the observed constraint value xo,i is not large, but the uncertainty of this proxy estimate is even larger.
Continuing, we can write
Equation (30) has the form of weighted least squares estimation of the mean y of a Gaussian random variable Y based on samples yo,i, i = 0, …, n, with standard errors . Each “sample” corresponds to a proxy, except sample 0, which corresponds to the GCM prior. This is a standard problem in statistical estimation [e.g., chapter 3 of Strutz (2016)]. The resulting posterior PDF of Y is a Gaussian with mean and standard deviation
The normalized weights wi/ws sum to 1, so the posterior mean is a weighted average of estimates yo,i derived from individual constraints.
If several constraints all suggest large anomalies yo,i of the same sign, they can overwhelm the anchoring contribution of the GCM prior and combine to create a large , suggestive of a Y that is farther from the GCM mean than most or all of the single-constraint estimates, which are more strongly anchored to the GCM prior.
To connect to Method C, the Method U posterior mean of Y can also be written as a linear combination of the constraint values:
with constraint weights
The Method U constraint weights ai are appealingly interpretable. They are positive, and they increase with the correlation coefficients ri of the constraints with the climate proxy.
More conditionally uncorrelated constraints always increase ws, decrease , and narrow the posterior range of Y, especially those having a high correlation with Y.
a. Consistency with single-constraint posterior PDF
When applied to only one constraint, Methods C and U both reduce to the single-constraint posterior. Indeed, the Method U formulas (after dropping the index 1 to match the single-constraint notation)
match our single-constraint formulas in Eqs. (10) and (11). However, Method U provides an interpretation that generalizes to multiple constraints, namely that the posterior mean of Y is a weighted average of its GCM prior (which is zero), and the proxy estimate yo = x0/r from the constraint, multiplied by a normalized weight r2 < 1. In this interpretation, the contribution of the GCM prior “dilutes” the effect of the constraint by a factor of r2.
b. Method U results for ECS
The blue curve in Fig. 4a shows the Method U estimate of the ECS posterior given the four credible constraints. This Gaussian PDF has a ±2σ range of 3.3–4.9 K, similar to the 4-constraint estimate from Method C. That is, partial correlation among the four credible constraints has almost no impact on the posterior PDF. Indeed, according to Method U, the most likely ECS is
The coefficients are very similar to those in the analogous formula in Eq. (21) derived from Method C. The blue curve in Fig. 4c shows the Method U ECS posterior for all 11 constraints, which has a ±2σ range of 3.4–4.6 K. This is similar but slightly narrower than the 4-constraint posterior.
Method U is consistent with Method C, especially for the 4-constraint case, and has the attraction of interpretability. Figure 5a shows one way to visualize its results. Each constraint is represented by a colored vertical bar with height equal to the weight wi, located at an ECS redimensionalized from yo,i using the GCM mean and standard deviation. For any set of constraints, the posterior mean ECS is at the horizontal center of gravity (weighted average) of their bars, including the GCM prior (constraint 0). The 11-constraint posterior mean and uncertainty range are shown for comparison.
Recall that for each constraint, the weight wi is a strongly increasing function of the correlation coefficients ri with ECS. Most of the constraints have correlation coefficients with ECS of between 0.4 and 0.6 and modest weights of 0.1–0.5. Constraints 3 (Zhai) and 4 (Brient Alb) have larger weights near 1 because of their stronger correlation with ECS, so they are particularly influential to the ECS posterior range. “Constraint 0,” the GCM prior, has a weight of 1. With a single constraint, the GCM prior substantially moderates how far the posterior mean is likely to deviate from the GCM mean ECS, as we saw in Fig. 2. In combination, the four credible constraints significantly outweigh the GCM prior, reducing its influence. With 11 constraints, the GCM range of ECS only plays a modest role in shaping the posterior.
The observational proxy estimates yo,i span a wide range, but 8 of 11 exceed the GCM prior of 3.2 K and their 11-constraint median is nearly 4.5 K, so it is not surprising that the Method U all-constraint mean ECS estimate is 4 K. This is larger than the average of the single-constraint estimates, because those are more strongly anchored to the GCM prior. Constraint 1 (Sherwood D) has an anomalously large ECS estimate of nearly 8 K because its observational estimate is xo,1 = 2.4 standard deviations above the GCM mean, and because it has a relatively small correlation r1 ≈ 0.4 with ECS. However, that small r1 also causes it to have a small proxy weight w1.
Figure 5b shows the construction of the Method U posterior PDF of ECS as a product of independent Gaussian PDFs from individual constraints. This is easiest to visualize logarithmically as a summation of log-likelihood (the exponent of the Gaussian), which is a parabolic function of ECS for each constraint. The color shades show log-likelihood contributed by the GCM prior and each of the four credible constraints (adding the other seven constraints is trivial but would clutter the plot). Each constraint favors a different ECS range centered around its observational ECS estimate plotted in Fig. 5a, over which its log-likelihood is not too negative. The GCM prior weights the ECS toward the GCM-mean ECS of 3.2 K, but constraint 1 (Sherwood D) favors very high ECS, while constraint 2 (Brient Shal) favors low ECS. The resulting 4-constraint posterior log-likelihood is the black parabola bounding the lowest shaded region. For both this and the 11-constraint ECS posterior, no single constraint is responsible for the mean exceeding 4 K; that is a collective result of several constraints, lending more credence to the result.
8. Adjusting for overconfidence
We present a method to adjust Methods C and U for overconfidence by artificially reducing the correlation coefficients ri between the constraints and Y, while leaving the partial correlation coefficients between the constraints unchanged.
An obvious way to do this would be to scale all the correlation coefficients ri by the same factor. We use a related but slightly different approach. A constraint that is nearly perfectly correlated with ECS over 20+ GCMs should also be extremely highly correlated with ECS over a different set of GCMs. Thus we treat overconfidence as leading to a systematic underestimate by a user-specified factor α2 ≤ 1 of the variance in each of the constraints that is unexplained by regression onto Y. This amounts to scaling the ratio of the explained to unexplained variance by the specified factor α2:
for which the correspondingly reduced correlation coefficients are
For small values of ri, the reduced correlation coefficient is a factor of α as large as the original coefficient; for large ri, the fractional reduction is less.
To use Method C, we could adjust the correlation matrix to preserve the partial correlations rij0 between constraints while reducing the ris. However, if we are highly uncertain about ri, then our empirical estimates of the partial correlations are at least as uncertain and potentially meaningless.
Some overconfidence correction (α < 1) seems merited given our earlier arguments about a priori selection and optimization of constraints. One compelling basis for such a judgement is testing based on an independent set of GCMs. Some previous examples include some of the discredited constraints in CZK18, and the minimal correlation of the “credible” Sherwood D constraint with ECS across a perturbed parameter ensemble (Wagman and Jackson 2018). This latter work suggests a large uncertainty enhancement for constraint overconfidence is needed even for physically appealing emergent constraints. Further research on this issue is needed.
As a sensitivity test, we apply an extreme overconfidence adjustment α = ⅓ (a ninefold inflation of unexplained variance) to Method U. The green curves U3 in Figs. 4a and 4c show the resulting ECS PDFs for the 4-constraint and the 11-constraint cases. In both cases, the overconfidence correction widens the posterior range (Table 4) and moves the posterior mean ECS somewhat toward the GCM mean. The wider range is expected, because constraints weakly correlated with ECS are less informative.
We derive a new “Method C” of combining correlated emergent constraints to estimate a Gaussian PDF for ECS or any other global climate sensitivity proxy from a finite sample of GCMs. It accounts for observational uncertainty of the constraints and (optionally) for overconfidence in the estimated correlation between the constraints and the proxy. The method is based on approximating the joint PDF of the proxy and the constraints as multivariate Gaussian, and it can also be framed in terms of multiple linear regression. With our limited sample of 40 CMIP GCMs, most of which provided inadequate outputs to compute some constraints, this PDF is inadequately sampled, and therefore the method may not give robust results as the number of constraints becomes comparable to the number of GCMs. Ridge regression is shown to improve the robustness of the combined prediction but somewhat biases the estimated PDF of the climate proxy toward the GCM mean. We also present “Method U,” which neglects any partial correlations between emergent constraints that are not related to their joint correlation with the climate proxy; this is more robust when applied to a small sample of GCMs and is more interpretable.
We apply these methods to a set of four credible and seven other “possible” constraints from CZK18 and compare them with PDFs derived from single constraints. Taken singly, the 11 constraints imply ECS PDFs with ±2σ ranges having means between 3 and 4 K and ±2σ widths of ±0.8–1.2 K. Reassuringly, all of these constraints give overlapping PDFs for ECS, with eight of the 11 favoring ECS higher than the GCM mean. The 4-constraint ±2σ ECS range estimated by both Methods C and U is close to 4 ± 0.8 K, which lies exclusively in the upper half of the GCM range.
Using all 11 constraints, the Method C ECS range showed some sensitivity to omission of individual GCMs from our dataset but universally favored ECS well above the GCM mean. With a larger sample of GCMs, the robustness of Method C should further improve. With 11 constraints, Method U produced a similar but narrower ECS range when compared with 4 constraints. Relative to the single-constraint analysis, the Method U multiconstraint analysis generates PDFs that are less like the GCM prior. Since a large majority of the constraints individually favor ECS values above the GCM mean, combining them more strongly favors that result. This is an important and interesting result of this analysis.
We propose a user-chosen adjustment factor α to account for constraint overconfidence. This factor reduces the correlation coefficients of all the constraints with the climate sensitivity proxy while leaving their conditional correlations with each other unaltered. With a strong overconfidence adjustment α = ⅓, the ±2σ 11-constraint estimated ECS range of Method U doubles in width to 2.7–5.3 K, but its mean remains close to 4 K.
The same hierarchy of approaches was applied to the climate sensitivity parameter λ = −3.7 W m−2/ECS, with similar results shown in Table 4. Arguably, it is more plausible to fit λ with a Gaussian distribution than ECS. A Gaussian λ implies a skewed distribution for ECS. Indeed, the ±2σ range of λ for Method U using 11 constraints implies an ECS range of 3.4–6.1 K, that is, a substantially longer positive tail and a slightly truncated negative tail relative to the alternate assumption of Gaussian ECS.
We conclude that climate sensitivity estimated from combining the most reasonable current emergent constraints is very likely above the CMIP3/5 GCM mean of 3.2 K and has roughly even odds of exceeding 4 K. To better interpret and bolster this provocative result, we should continue to search for more physically motivated emergent constraints aimed at regime-specific cloud feedbacks (e.g., Qu et al. 2013). We will also examine how the 11 considered emergent constraints fare in predicting ECS and other measures of climate sensitivity across the GCMs in the ongoing CMIP6 intercomparison, several of which have ECS well in excess of 4 K.
We acknowledge the modeling groups, the Program for Climate Model Diagnosis and Intercomparison (PCMDI), and the World Climate Research Program’s Working Group on Coupled Modeling for making available the CMIP3 and CMIP5 multimodel datasets, supported by the U.S. Department of Energy (DOE) Office of Science. This work was supported by the Office of Science (BER) at Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344 and BER’s Regional and Global Climate Modeling (RGCM) Program. Steve Klein instigated this study and heavily influenced early drafts. Yuying Zhang provided estimates of the Volodin constraint from multiple satellites.