Despite major advances in climate science over the last 30 years, persistent uncertainties in projections of future climate change remain. Climate projections are produced with increasingly complex models which attempt to represent key processes in the Earth system, including atmospheric and oceanic circulations, convection, clouds, snow, sea-ice, vegetation and interactions with the carbon cycle. Uncertainties in the representation of these processes feed through into a range of projections from the many state-of-theart climate models now being developed and used worldwide. For example, despite major improvements in climate models, the range of equilibrium global warming due to doubling carbon dioxide still spans a range of more than three. Here we review a promising way to make use of the ensemble of climate models to reduce the uncertainties in the sensitivities of the real climate system. The emergent constraint approach uses the model ensemble to identify a relationship between an uncertain aspect of the future climate and an observable variation or trend in the contemporary climate. This review summarises previous published work on emergent constraints, and discusses the huge promise and potential dangers of the approach. Most importantly, it argues that emergent constraints should be based on well-founded physical principles such as the fluctuation-dissipation theorem. It is hoped that this review will stimulate physicists to contribute to the rapidly-developing field of emergent constraints on climate projections, bringing to it much needed rigour and physical insights. Numerical methods have become a standard technique to simulate complex systems. The equations governing components of such systems may be well-known. But their solutions cannot be solved analytically, creating a need for numerical approaches. Because of the discretization of time and space inherent in numerical techniques, modelling complex systems must involve ways to include effects of unresolved processes. Often there is no 'first principles' approach to do this. Typically, the effects of unresolved processes are included by resorting to quasi-empirical relationships between them and explicitly-resolved variables, otherwise known as 'parameterization'. There are usually multiple defensible ways to parameterize unresolved processes. Thus, independently developed models of the same complex system might incorporate different parameterization choices. The more models that are independently developed, the greater the diversity of approaches for modelling the same natural system. A classic example of this model diversity is the use of numerical models of the atmosphere to predict hurricane development. Initial conditions are imposed on a model at some time, and it is integrated forward in time to produce simulations of critical hurricane features such as intensity and track. This approach might be replicated for multiple models incorporating different parameterization choices, producing an 'ensemble' of hurricane forecasts. The spread in the forecasts is a measure of the uncertainty in the future hurricane behaviour, given the range of plausible approaches to atmospheric modelling. Over time, with enough hurricanes and associated predictions, the various models can be evaluated for their prediction skill. Certain parameterization choices may emerge as producing systematically better predictions. The models can then be rebuilt or re-calibrated with the better choices. Over time the ensemble will become more skilful, with less spread. In fact, this is approximately the process that has resulted in dramatic improvements not only in hurricane prediction, but in weather forecasting generally over the past seven decades. Earth's climate is another example of a complex system whose governing equations can only be solved through numerical methods. (In fact, the dynamical equations for the atmospheric component of a climate model are typically almost identical to those in the hurricane models referred to above.) As expected, there are a variety of plausible approaches to parameterization in the components of Earth system models (ESMs) 1 . Thus, modelling groups throughout the world have built a few dozen ESMs with different approaches to parameterization. Because of these differences, these models produce different future climate states, even when the same scenario of radiative forcing is imposed (associated, for example, with an increase in greenhouse gases). A classic climate change experiment is to double CO 2 concentrations in the atmospheric component of an ESM, and measure the surface warming that occurs after the simulation has equilibrated (Manabe and Wetherald, 1975), an important number in climate science referred to as equilibrium climate sensitivity (ECS). Reaching a true equilibrium with a full complexity ESM requires long, computationally expensive simulations of thousands of years so ECS is usually estimated from shorter duration CO 2 doubling experiments. When these experiments are done with contemporary ESMs, the spread in values across the ensemble of ESMs is large, between 1.5 and 6 degrees Celsius (Forster et al., 2019), although true equilibrium values are higher (median 17% higher, Rugenstein et al. (2020)). The international climate science community has organised itself to generate scenarios of greenhouse gas (GHG) emissions that result in more realistic future radiative forcing than the CO 2 doubling experiment. These correspond to scenarios of controls (or lack thereof) on future GHG emissions (Moss et al., 2010;Riahi et al., 2017). The scenarios are designed so that the same radiative forcing time series is imposed on each ESM. This allows for inter-comparison of ESM response at some future specified time, say the end of the 21st century. The ESM responses vary significantly across the ensemble, an indicator of deep uncertainty in Earths' climate future, simply due to the variety of plausible ways to handle parameterization within ESM components. If each ESM is treated as being an equally plausible analogue of the real climate system, known as 'model democracy', an undesirably large factor of three difference in future response to man-made GHG response currently results. Apart from the subject of this review, there are several techniques used to reduce this uncertainty and much research effort has been dedicated to this important task. One approach is to give ESMs that simulate the real world 'better' (according to some chosen metric) more weight in the future projection. There are also methods that do not use ESMs at all, by constraining future response from past observations and those that use a combination of all of these. Many more details can be found in e.g. Eyring et al. (2019); Sherwood et al.
(2020) and chapters 9 and 11-14 in IPCC (2013). The approach taken to improve hurricane forecasting models -to evaluate their performance and adjust parameterizations over the course of multiple prediction cycles -is unfortunately impractical in the case of projecting future climate with ESMs. With future climate, we have only one realisation of the real system's trajectory, and our mandate is to predict it as best we can now. To compound the problem, the climate system involves many more components than the atmosphere alone, including the ocean, the land surface, glaciers and ice sheets, and the marine and terrestrial biospheres. The engineering required to model each component and allow them to interact in a simulation is impressive. But the components can interact with one another in information flows are so complex that it is not easy to predict how behavior in one component might affect behavior in another. For example, parameterization reformulation or adjustment in one component can significantly affect the simulated state of another component for reasons that are not always clear (e.g. Donner et al. (2011)). To circumvent the impossibility of directly evaluating ESMs for their ability to simulate a future climate state which has not yet been observed, a new technique has emerged over the past decade and a half, known as the emergent constraint (EC) approach . The basic idea is to identify an element (X) of the observable climate that varies significantly across the ESM ensemble, and which exhibits a statisticallysignificant relationship with variations in some important variable (Y ) describing the ESM's future simulated state. If we call this relationship f , then Y = f (X) + , where  is a relatively small departure from f . Since X is an element of the observable climate, it is a quantity that can be measured. The relationship f may then place a E m e r g e n t R e l a t i o n s h i p inter-model uncertainty 1 Schematic showing the most common procedure used to derive emergent constraints on Earth system sensitivities. An ensemble of ESMs (each red dot is an individual ESM) running the same experiment (the pdf on the right-hand yaxis represents the spread in the ensemble) is used to identify an emergent relationship (black dashed line with grey uncertainty range), between an uncertain Earth system sensitivity Y (y-axis) and an observed trend or variation X (x-axis). An observation of the trend or variation (blue pdf on the x-axis) can then be combined with the model-based emergent relationship to derive an emergent constraint on the Earth system sensitivity (green pdf on the left-hand y-axis). useful constraint on Y , provided the measurement uncertainty in X is small compared to the range of simulated values. This constraint is 'emergent' because the emergent relationship f cannot be diagnosed from a single ESM. It only becomes apparent when the full ensemble is analysed. If the relationship f arises from model physics or dynamics common to the ESMs, then reducing the spread in X through reformulation or adjustment of parameterizations ought to result in spread reduction in Y . Corresponding reduction in the spread in other future climate variables affected by Y should also occur. If the process is repeated for enough variables X and Y , then we can imagine that overall simulated spread in multiple aspects of future climate would gradually be reduced. We also note that ECs can be used directly to ascertain the most likely values for a particular Y . To further illustrate the concept, the subject of the first ). In this case, the future climate variable of interest (Y ) is the snow albedo feedback (SAF), a climate mechanism characterized by the retreat of highly reflective snow cover under climate warming and the associated reduction in surface albedo, which amplifies warming and pro-motes further snow melt (Bony et al., 2006;Hall, 2004;Thackeray and Fletcher, 2016). However, modern ESMs disagree on the strength of this feedback, exhibiting a nearly threefold spread (Qu and Hall, 2014;Thackeray et al., 2018). SAF also occurs each spring in the current climate when Northern Hemisphere snow cover recedes from its winter peak to its summer minimum, enhancing seasonal warming in the process (X). A comparison between simulated SAF under climate change (Y ) and in the seasonal cycle context (X), uncovers a strong linear relationship in three different generations of ESM ensembles . The strength of this relationship tells us that seasonal SAF is likely highly predictive of climate change SAF, while the consistency across model generations illustrates robustness to out-of-sample testing. Furthermore, we can quantify the seasonal SAF using satellite-derived observations, thus allowing for model bias to be properly assessed. For example, an ESM that underestimates observed SAF in the seasonal context is likely to underestimate SAF in future climate. Producing a more accurate SAF in ESMs is more than just an academic exercise as variability in SAF can account for a significant portion of uncertainty in 21st century projections of warming across Northern Hemisphere extratropical land (Qu and Hall, 2014). Thus, by reducing variability in SAF we can expect similar reductions in projections of regional warming. Figure 3 highlights the first application of this and several other ECs from the literature (pertaining to the carbon cycle and climate sensitivity), all of which will be further discussed in Section V. In this review article, we synthesise what is now understood about ECs. We discuss how relationships f might emerge in ESM ensembles in section II, as well as the theory that underlies ECs in section III. In section IV we note pitfalls associated with ECs and how to guard against them. In section V we survey the ECs that have been found in the climate system, and describe the statistics that can be used to apply them rigorously in section VI.
Finally, we discuss how techniques similar to ECs have been or might be used in other fields, and consider the outlook for development of the technique going forward in section VII before concluding in section VIII. ECs are possible because of 'emergent relationships' appearing in an ensemble of ESMs. An emergent relationship, Y = f (X) + , between an element of the observable climate X, something that can be measured now in the real world, and the unknown future response Y one wants to know, can then be used to place a constraint on the real world value of Y via observations of X to give an EC. For example, in Cox et al. (2018a) the observable was the variability in global annual mean sur- face air temperature during roughly the last one hundred years and the unknown future response was equilibrium climate sensitivity (ECS, defined as the magnitude of the Earth's warming to doubled CO 2 levels), one of the oldest and most important numbers in climate science. Despite decades of research the latest ESMs still vary widely in their predictions of this number by a factor of three (Forster et al., 2019;Zelinka et al., 2020). Williamson et al. (2019)). This mechanism from high to low dimensional relationship emergence and others are given in the following subsections. First we discuss the differences in the types of climate model ensembles commonly used to find emergent relationships. Model simulations are not perfect reproductions of the system they are designed to emulate. Ensembles of models are used to get a handle on (i) initial condition uncertainty, (ii) parametric uncertainty and (iii) structural uncertainty. Uncertainty in the initial state is particularly important for numerical weather prediction. As the equations governing the weather are chaotic, forecast solutions with very small differences in their initial states, equally likely to be the 'real' initial state, can diverge strongly after just a few days. This sensitivity to initial conditions means running the weather model just once, even if the model is a perfect reproduction of the real world, may produce a very different forecast from the one actually experienced. By running the same weather model many times with different, but equally plausible initial states chosen to sample the most unstable and divergent regions of the model phase space, the likelihood of experiencing a particular weather forecast can be constructed. Ensemble forecasting is now a standard tool in numerical weather prediction (Epstein, 1969;Leith, 1974;Molteni et al., 1996;Toth and Kalnay, 1993). Initial value ensembles in climate applications are usually used to assess natural climate variability (Kay et al., 2015;Maher et al., 2019). Climatic variables are (generally speaking) the long term statistics of weather and their prediction is less impacted by initial state uncertainty relative to parametric and structural uncertainties (Hawkins andSutton, 2009, 2011). Parametric uncertainty arises from uncertainty in the values of constants ('parameters') in quasi-empirical relations ('parameterizations') used to model the effects of unresolved, subgrid scale processes. Parametric uncertainty can be sampled by running the same model with different values of these parameters. ESM model ensembles that do this are known as perturbed physics or perturbed parameter ensembles (PPEs, Murphy et al. (2004)). Examples of ESM PPEs are the climateprediction.net experiment (Stainforth et al., 2005) and the Met Office's Quantifying Uncertainties in Model Projections (QUMP) ensemble (Collins et al., 2011) that both used the HadCM3 ESM. Structural uncertainty arises from uncertainty in the functional form of equations. Although a lot of the model equations in ESMs are well known, there are multiple, equally defensible parameterization schemes for unresolved processes (e.g. convective precipitation and cloud radiative properties) and so different equations can appear in different ESMs. ESM ensembles that sample structural uncertainty are called multi-model ensembles (MMEs) and include the coupled model intercomparison project (CMIP) ensembles (Eyring et al., 2016;Meehl et al., 2007a;Taylor et al., 2011) used to inform the Inter-governmental Panel on Climate Change (IPCC) reports (Collins et al., 2013;Cubasch et al., 2001;Meehl et al., 2007b) and used extensively in much of the EC literature. The CMIP MME ensembles consist of a number of different ESMs developed by different international groups thus sampling structural uncertainty. Some CMIP groups also submit multiple initial value runs and occasionally multiple physics parameter runs performed with the same model, sampling initial condition and parametric uncertainty to some degree. Numerical weather prediction ensemble forecasting has also incorporated more parametric and structural uncertainty into their forecasting as time has gone on (Palmer, 2019). Most ECs have derived from MME ensembles to date although there are some examples from PPEs e.g. Knutti et al. (2006). ECs derived from MMEs are, in our opinion, the most believable. Each model in a MME could be thought of as that particular group's best guess of the Earth system. PPEs on the other hand, generally vary only a handful of parameters and thus sample just a few dimensions of the model's phase space. Due to the few varied dimensions, relationships hard coded in the fixed equations between model variables become easy to find (based on the authors experience with the QUMP PPE (Lambert et al., 2013)). These ECs need to be backed up with physical arguments and ideally tested in MMEs. Collins et al. (2011) found that long wave cloud feedback was highly correlated to climate sensitivity in the QUMP HadCM3 PPE ensemble, however this correlation was absent in (smaller sample size) MMEs. Yokohata et al. (2010) also found that strong correlations in a PPE with HadSM3 were not always present in the equivalent MIROC3.2 PPE. Further, Yokohata et al. (2010) found the mechanism for variations in climate sensitivity was different. Values of parameters chosen to correctly give the proper balance between opposing processes in MMEs could be skewed to one of the processes by simultaneously varying them in a PPE thus reducing the realism of each ensemble member solution. However, PPEs can be useful for testing ECs (e.g. Kamae et al. (2016); Wagman and Jackson (2018) and section IV.I) and quantifying parametric uncertainty (section VII.F).
A starting null hypothesis is that emergent relationships occur by chance and are not indicative of a deeper mechanistic relationship (Hall et al., 2019). Data mining an ESM ensemble for high correlations between pairs of variables may fall into this category (see also section IV.A). If an ESM ensemble is reasonably small and the number of variables output by each ESM is high, the expected number of variable pairs with high correlation is high purely by chance (Caldwell et al., 2014). More precisely, an ESM labelled by the index i in an ensemble of n ESMs i = {1, 2, . . . , n} calculates a large number m of different variables X l , each variable labelled with index l and l = {1, 2, . . . , m}. In this data mining scenario, a pair of variables (X l i , X k i ) is chosen for the prospective emergent constraint, one for the observable X l and the other, the response X k one wishes to know in the real world (previously written Y ). In an indiscriminate data mining approach, every possible pair of variables are used to create 1 2 m(m  1) data sets, each labelled S kl . Each dataset has n elements, with each element corresponding to one of the n ESMs. If the pair of variables making up S kl is a good candidate for an emergent relationship and we assume this relationship is linear, then the correlation r kl in S kl should be high. Correlation here is defined as r kl = cov(X l , X k )/ X l  X k and  X is the standard deviation of X. The number of ESMs in a model ensemble, n, is typically small being around 10 to 40 in the state-of-the-art CMIP ensembles although numbers increase with each successive generation (Eyring et al., 2016;Meehl et al., 2007a;Taylor et al., 2011). If we calculate the correlation r kl within every one of the 1 2 m(m  1) possible datasets, some will have high correlations purely by chance. The likelihood increases as the number of models in the ensemble n gets smaller giving less data points in each S kl . While the fraction of datasets with high correlations |r| > |r high |, does not increase as the number of variables m gets larger, the total number of datasets with correlation above |r high | will increase simply because there are more possible variable pairs to correlate. Thus one must be careful that any prospective emergent constraint makes physical sense, particularly in current, small ESM ensembles to avoid this pitfall. To illustrate this point in a worst case type scenario where every ESM output variable is uncorrelated to any other output variable or any of the other ESMs, we can calculate the correlations expected between all possible pairs of variables as a function of ensemble size n. Figure 4 shows the fraction of all possible datasets with a particular correlation for a 'toy' ESM ensemble of n models where each of the toy ESMs is a modelled as a random number generator. This is a rather extreme and unrealistic example, each ESM producing m random output variables, the lth one labelled X l i , drawn from a normal distribution with unit standard deviation, X l i  N (0, 1). Even though there are no mechanistic relationships in this example by construction, a fraction of the datasets have strong correlations, particularly in smaller ensemble sizes (small n). This is purely by chance. Each ESM is not a real ESM; it is represented by a random variable drawn from a normal distribution of unit standard deviation and each of the n ESMs produces m = 100 output variables giving 4950 different possible variable pairs to correlate with each other. This example is designed to be a worst case scenario yet still shows calculated correlation even between uncorrelated objects can be large in small enough ensembles. In the left panel, histograms of fraction of datasets with correlation falling in a particular interval are plotted. In the right hand panel, the total fraction of datasets with magnitude of correlation greater than |r| is plotted. Ensemble size is varied for toy ESM ensemble sizes of n = 5, n = 30 and n = 100, C. Low dimensional relationships 'emerge' from high dimensional ESMs In the last example even though all n x m 'toy' ESM variables in the ensemble were uncorrelated from each other by construction, high correlations were still possible. This is more likely in small ensembles. In reality there are mechanistic relationships between X and Y in any particular ESM. Real ESMs are not random number generators, the output of a ESM is the numerical solution of a large set of coupled, nonlinear equations codifying relationships between variables representing the ocean and atmosphere dynamics and thermodynamics and biogeochemical feedbacks 2 . However the exact relationship between any pair of X and Y is not solely a function of just these two variables, it is generally also a function of many others that prescribe the exact state of the ESM at every point in space and time on the ESM grid. The minimal number of (non-unique) variables, d, required to describe the exact state of the ESM at time t can be thought of as a point in d dimensional phase space with coordinates (x 1 , x 2 , . . . , x d ) and in the case of a ESM, the dimension of this space is extremely high 3 . As the ESM state changes with time, a trajectory is drawn between these phase space coordinates connecting the past and present states. Variables describing the precise state of the ESM such as temperature or precipitation, the X and the Y , may be functions of all the d phase space coordinates as well as time (d+1) i.e.
X = f X (t, x 1 , x 2 , . . . , x d ) and Y = f Y (t, x 1 , x 2 , . . . , x d ). Although the mechanistic relationship between X and Y is generally a function of all d + 1 coordinates, climatic variables of interest are usually, although not exclusively, long term temporal and/or spatial averages that can be thought of as (time invariant) attractors of the phase space reducing the information relative to a precise configuration at any one time. In many circumstances the mechanistic relation between X and Y can be approximated to a good degree by a much smaller dimensional subset of the full phase space reducing the effective dimension. The success of statistical mechanics and science in general 4 relies on this effective dimension reduction, that certain properties of systems with many degrees of freedom can be well approximated by fewer state parameters. For example, a container of gas is, as far as we know at present, most completely described by the equations of quantum mechanics. In principle, one could solve these equations for each of the N gas molecules, at each position and time in the container. The dimension of this space increases exponentially with N , requiring a huge number of degrees of freedom to specify the state exactly even for just a handful of molecules. This could be useful if the question one wanted to answer required detailed knowledge of the full quantum state. However, if we 'only' want to know the configuration of positions and momenta of the gas molecules, one can use an approximation of the quantum theory. Treating the gas as a bunch of identical spheres each with a mass and a different position and momentum, the state can be described by a point in 6N dimensional phase space, a reduction in the phase space dimension yet still a massive number of dimensions for a realistic gas. Generally detailed knowledge of a precise system state is not desired or required. Bulk properties describing 3 To give an idea of a lower bound on the dimension of the phase space required to describe a typical CMIP6 era ESM, the atmosphere is typically cut into a total number of  200 x 200 x 100 = 4 x 10 6 spatial cubes. A single atmospheric prognostic variable will need at least this many phase space dimensions at each point in time to specify it exactly. A lower estimate for the number of prognostic variables required just to specify the physical state of the atmosphere is O(10). This already gives d  10 7 even before including the ocean, land and biogeochemisty. 4 Low dimensional models can often mimic the responses of more complex ones to a good approximation and this is essentially why science is successful i.e. we can comprehend and predict the high dimension real world using (relatively) simple models. the mean overall state are a useful way of reducing information and making sense of complex systems. Going back to the container of gas, kinetic theory provides a way of mapping the many degrees of freedom from each of the molecules to the few bulk parameters of the ideal gas law, under reasonable assumptions (see for example Schekochihin (2015)). That is, the 6N dimensions effectively reduces to the three encoded by the relationship pV = nRT . The bulk variable pressure p felt as a force per unit area on the container wall and is related to how often and how vigorously on average molecules with mean square velocity v 2 hit the sides of the container (p = n 3V v 2 , n is the number of moles in the container of volume V ), whereas the temperature T is a measure of the most probable kinetic energy of a molecule in the container. R is the ideal gas constant. This approach works well when N is large as one finds more of the 6N states correspond to the ideal gas law. Effective dimension reduction in ESMs can also happen when we want answers to questions about bulk mean states of the climate. For instance, relations between bulk variables such as the annual global mean surface air temperature and the net incoming radiation are well modelled by the few degrees of freedom in a simple energy balance model (see for example Gregory (2000); Hasselmann (1976); and Wigley and Raper (1990)) even though the exact state of a ESM is given by all prognostic variables at every point and multiple times on its spatial grid. If one wanted to ask detailed questions about the state of a ESM such as 'How many days in March does it rain more than x mm in Madagascar?' then the full ESM theory is best placed, however if the question was 'How much does the global mean temperature increase from a spatially uniform doubling of atmospheric CO 2 ', a lot of that extra detail seems unnecessary. The analogy we attempt to present here is that ESMs are analogous to the solutions to the laws of quantum mechanics and the effective, reduced dimension relationship for a particular ESM i, Y = f i (X) + is analogous to the ideal gas law. The analogy is far from perfect, ESMs are of course not as simple as containers of large numbers of identical molecules and there is, as yet, no such elegant route between the two solutions provided by the statistical mechanical or kinetic theory recipes. The best candidates for ESM ideal gas laws to serve as theoretical bases for emergent relationships are simplified, analytically soluble models of the climate or its subsystems such as the energy balance models previously mentioned (Gregory, 2000;Hasselmann, 1976;Wigley and Raper, 1990). These act more as testable, plausible hypotheses of the bulk behaviour of ESMs rather than derivations from the basic laws that constitute them.
So far we have discussed how a physically plausible, mechanistic, low dimensional relationship f i between X and Y may appear from a particular ESM i, Y = f i (X) + . We have not yet discussed how low di-mensional relationships, the emergent relationship, could appear across an ensemble of ESMs. Although each ESM is different, they should be equally plausible models of the (X, Y ) relationship in the real world and it should not be a surprise that ESM solutions across a model ensemble should be also be mechanistically related. This is because many of equations are based on well established and tested physics (such as the Navier-Stokes equations) and are common amongst ESMs. Indeed, if the all the equations were well known and could be numerically integrated exactly, the ESMs should produce identical outputs (and the problem of climate modelling would be solved and emergent relationships would disappear). However, some equations can not be exactly numerically integrated due to limitations on temporal and spatial resolution. One then has to parameterize subgrid scale processes and this may be a source of difference between ESM solutions, each model giving a different (X, Y ) pair (see section II.A). Although they are our best attempts at understanding the climate, some ESM responses diverge on important questions. Such a scenario is where the emergent constraint approach may be helpful. To have a chance of finding an emergent relationship, the set of n (X i , Y i ) data points should be different (i.e. (Y 1 , X 1 ) = (Y 2 , X 2 )) and span a wide enough range for the emergent relationship to become apparent. The model relationship f i also needs to be shared amongst the model ensemble (f i (X)  f j (X)i, j = {1, 2, . . . , n}). Provided such a relationship exists, the emergent relationship f (X) can be determined. If this model relation is also shared by the real world, the real world's response can also be determined (along with its uncertainty) to give an emergent constraint (see section II.E). A related mechanism for a low dimensional relationship to emerge between X and Y occurs when the magnitude of the same physical process, correlated to the observable X, differs appreciatively across the ESM ensemble and the magnitude of this process heavily determines the size of the response Y . This is also a case of effective dimension reduction -the range in response is dominated by the dimensions of phase space of just that physical process. An example of this happening was provided by Caldwell et al. (2018) who evaluated ECs on equilibrium climate sensitivity (ECS) in the CMIP3 and CMIP5 ESM ensembles. That is, the targeted unknown response, Y , amongst these ECs was the same, ECS, although different observables X were found to be well correlated to it. Examples of a few of the X used in the different ECs were the strength of resolved-scale mixing between the boundary layer and lower troposphere in the trop-ical east Pacific and Atlantic (Sherwood et al., 2014), error in the distribution of cloud-top pressure and optical thickness for regions between 60 o N/S (Klein et al., 2013), fraction of tropical clouds with tops below 850 mb whose tops are also below 950 mb (Brient et al., 2016) and variability in global mean air temperature (Cox et al., 2018a). Caldwell et al. (2018) showed all of these observables were also highly correlated to the short wave (SW) cloud feedback, a strong and uncertain feedback on the resulting global temperature. The SW cloud feedback in turn was shown to be highly correlated to intermodel variations in ECS within the CMIP5 ensemble. To paraphrase Caldwell et al. (2018), '...intermodel variations in cloud feedback were so big that they left a strong imprint on intermodel variations in ECS. This means that fields that are strongly correlated with ECS are probably correlated with the SW cloud feedback (and vice versa).' E. What is needed for an EC? In addition to having an emergent relationship in a model ensemble, an EC also ideally requires a few other things to be useful. The range of X in the ESM ensemble should be large relative to the uncertainty in the value of the real world observable. Uncertainty in each of the model X should also be small relative to the range in model X. Ideally the real world value of the observable X obs should also fall well within the range of values in the ESM ensemble to avoid extrapolation issues. One would like a large range in the values of the ESM ensemble responses relative to the uncertainty in each of the individual ESM responses. This reduces the uncertainty in the real world value of Y derived from the EC. A particular ESM labelled i should have a physically plausible relationship f i between the observable and response Y = f i (X) + . It should be based on a simplified theory that predicts a functional response and has clear, testable assumptions that can be falsified independently. In addition, all models in the ensemble should share roughly the same relation between observable and response i.e. f i = f j . The stronger this relationship is, the more useful the emergent constraint due to the po-tential uncertainty reduction in Y . And this relationship should also be shared by the real world. Ideally you would have a large number of independent ESMs in your ensemble (larger sample, lower error in estimates, less chance for high correlations by chance). See section IV.E for further discussion of ESM ensemble independence. Many ;Kidston and Gerber, 2010;Lin et al., 2017;Lipat et al., 2017;Massonnet et al., 2012;Selten et al., 2020;Siler et al., 2018;Simpson and Polvani, 2016;Tian, 2015;Volodin, 2008).
However, a growing number of papers are also now relating Earth System sensitivities to observable temporal variations, such as trends (Boe et al., 2009;Jimenez-de-la Cuesta and Mauritsen, 2019;Nijsse et al., 2020;Tokarska et al., 2020), interannual variability (Clement et al., 2009;Cox et al., 2018aCox et al., , 2013;;Kwiatkowski et al., 2017;O'Gorman, 2012;Qu et al., 2015;Wenzel et al., 2014) and seasonal cycles (Hall and Qu, 2006;Knutti et al., 2006;Qu and Hall, 2014;Thackeray and Hall, 2019;Wenzel et al., 2016;Zhai et al., 2015). This section discusses the theoretical basis for emergent constraints based on such temporal variations. In the case of changes in relatively fast variables, such as seasonal snow-cover (Qu and Hall, 2007), or marine phytoplankton concentration (Kwiatkowski et al., 2017), there may be a fairly straightforward near one-to-one relationship between the short-term variability and the longer-term sensitivity, because the fast variable will be in a quasi-equilibrium state even with short-term climate variations. For slower variables (such as the forest carbon storage), short-term variations are more likely to measure fluxes (or equivalently the rate of change of the store). In this case, finding a constraint on future changes in the store requires multiplying the flux sensitivity to short-term variations by a characteristic timescale for each model. In some cases, the characteristic timescale may be similar across the model spectrum, leading to a simple emergent relationship between the short-term flux sensitivity and the long-term sensitivity of the store (Cox et al., 2013;Wenzel et al., 2014). In general though, converting a flux sensitivity to a store sensitivity requires an independent estimate of the characteristic timescale of the store, which itself requires a theoretical basis (Williamson et al., 2019). The Fluctuation-Dissipation Theorem (FDT) provides one such theoretical framework (Kubo, 1966). The FDT relates the sensitivity of a system to a small external forcing to the response of the same system to its own internally-generated fluctuations. FDT is therefore a potential theoretical basis for ECs, as it links the natural variability of a system to its sensitivity. The size of the forcing is important because the theorem only strictly applies to near-equilibrium linear systems. Nevertheless, FDT-based approaches have had a huge impact in statistical physics, including Einstein's work on Brownian motion (Einstein, 1905), and the understanding of Johnson-Nyquist noise in electrical circuits (Johnson, 1928;Nyquist, 1928) . The first proposal to apply FDT to the climate system came almost 40 years ago (Leith, 1975), but there has been a recent resurgence of interest in this area fuelled by methodological advances and detailed comparison of FDT-derived estimates of climate change to climate model simulations (Majda et al., 2010). In principle, it may be possible to use FDT to get good estimates of the response of the real climate system to small forcing (e.g. that due to doubling CO 2 ) purely from accurate long-term climate observations that reveal the full spectrum of natural fluctuations of the climate (Bell, 1980;Schwartz, 2007). Unfortunately though, the length of the detailed climate records required to achieve this are typically much longer than those available (Kirk-Davidoff, 2009), and include the contemporary period when the mean climate is changing. In the EC technique model projections are instead used to define emergent relationships between observable variations and future climate (Hall et al., 2019), and then specific observations provide a selection principle to constrain the range of future climate projections from these model-defined relationships (Eyring et al., 2019). The EC technique has yielded many proposed constraints on aspects of the future climate (Hall et al., 2019) and the carbon cycle (Cox, 2019). However, this power also comes with a danger -that blind data-mining of multidimensional model outputs could lead to spurious and misleading constraints (Caldwell et al., 2014). Protection from this risk can come from tests of robustness across different model ensembles, and from basing the search for ECs on firm theoretical and mathematical foundations. Theory-led ECs can also be considered as hypotheses that can be tested against the ensemble of complex ESMs. The remainder of this section aims to describe the reasons why relationships between Earth System variability and sensitivity are ubiquitous, and to provides some examples of the emergent relationships to be expected under different types of time-varying forcing (such as seasonal cycles and longer-term trends). We can write the time evolution of a dynamical system variable V in the general form: where F is the net flux into the system, and z i are the environmental variables that affect this net flux. Equation 1represents a local conservation law as it implies that the variable V can only change if the net flux is nonzero. When the environmental variables are in a steady state, z i = z i (0), the variable V has an equilibrium state V = V (0), defined by: In climate science we are typically interested in how a climate system variable, such as global mean temperature or land carbon storage in the tropics, varies with environmental factors, such as the atmospheric carbon dioxide concentration or the climate in the tropics. Where the associated perturbations z i are small compared to z i (0), equation 1 can be linearly expanded about the equilibrium state: where V is the resulting perturbation to the state variable, and the partial derivatives are calculated around the initial equilibrium {V (0), z i (0)}. The initial equilibrium is stable if F V < 0, and we can then define an effective timescale to perturbations of V as  = ( F V ) 1 . Equation 3 can therefore be written as: The righthand-side of equation 4 can be viewed as the 'external' forcing factors that produce changes in V (e.g. the radiative forcing due to increasing atmospheric CO 2 ). For a permanent time-invariant change in the environmental variable z i , V will be changed by an amount V eq where: Defining the sensitivity of V to the environmental variable z i as  i = V eq /z i , we can rewrite this equation as: where  i = F/z i .
Equation 6implies that the sensitivity of V to a step change in z i will be proportional to the sensitivity of the net flux F , with a constant of proportionality which is the characteristic lifetime of perturbations to V . In climate change research, the equilibrium sensitivities ( i ) are often the things we would most like to constrain (e.g. Equilibrium Climate Sensitivity, ECS). For emergent constraints we therefore need to find relationships between sensitivities  i and observable variations in the climate system. In this subsection we derive candidate emergent relationships for different timevariations in the environmental variables z i . First we rewrite equation 4 in the form: As we limit ourselves here to this linear model with a single timescale  , the emergent relationships presented below are intended to be illustrative, or as simple hypotheses to be tested against the outputs from complex models. As an idealised representation of the response of our system to diurnal and seasonal forcing, we first consider sinusoidal environmental variations of angular frequency  and amplitude a z , for which z i = a z e it . The solution to equation 7 under these circumstances (after initial transients have died-down) is: where  = arctan ( ). Therefore V also has a sinusoidal variation but with a phase-lag relative to the forcing of  which asymptotes to /2 as   . Substituting z i = A z e it . The amplitude of the sinusoidal variation in V , a v , is proportional to a v but also depends on the frequency of the sinusoidal forcing: Using equation 6 we can therefore write the sensitivity  i in terms of the ratio of the sinusoidal amplitudes: In the limit of slowly-varying forcing compared to the system timescale (  0) this formula reduces to : This is the limit of very fast variables that are in a quasiequilibrium with the sinusoidal forcing, such as seasonal snow-cover (Qu and Hall, 2007) or marine primary production (Kwiatkowski et al., 2017). However, in the opposite high-frequency limit ( 1) the emergent relationship also depends linearly on the frequency of the forcing: 2. Linearly-increasing forcing Contemporary climate change is largely being driven by an approximately exponential rate of increase in atmospheric carbon dioxide above the pre-industrial level, which yields a radiative forcing that is approximately linear in time. This has motivated a long-running series of idealised climate model experiments which prescribe an exponential increase in carbon dioxide, including the 1% per year runs which are used to define the concept of the Transient Climate Response (TCR) . In addition, many other Earth System changes can be approximated by linear trends. In this sub-section we therefore consider emergent relationships under a linear increase in the environmental variable z i : where  is the linear rate of increase z i which begins at time t = 0. Under this idealised forcing the solution to equation 7 is From equation 5 we note that the sensitivity  i =   i , so equation 14 can be rewritten as: The exponential decay-term here is a transient response to the sudden switch-on of the linear trend at t = 0. Once this term has died-out the transient solution lags the quasi-equilibrium solution by  years such that: This equation represents a potential emergent relationship between the sensitivity  i and a transient change V (t). Emergent relationships have also been proposed between interannual variability and sensitivities (Cox et al., 2013;Nijsse et al., 2019;Schwartz, 2007), assuming that the environmental variable z i is approximately gaussian white-noise. Under those circumstances equation 7 becomes the widely-used Ornstein-Uhlenbeck equation of statistical physics (Uhlenbeck and Ornstein, 1930). Standard solutions relating the variance of V to the variance of z i can be derived by integrating the sinusoidal solution given by equation 8 over all frequencies . For interannual variability this yields the following relationship for the ratio of standard deviations: and an even simpler relationship for the lag-1 autocorrelation of Y : These equations can be combined to yield an equation for the sensitivity  i without needing to know the timescale  (Cox et al., 2018a) : Similar formulae can be derived for more sophisticated representations, such as two-box and semi-infinite heat diffusion models of ocean heat uptake (Williamson et al., 2019). Uncertainty in future projections makes adaptation planning difficult, and so there is tremendous pressure on climate researchers to provide much more refined predictions of expected large-scale environmental change as atmospheric greenhouse gases (GHGs) rise. However, in the absence of a full knowledge of all climate processes and their parameterizations, such deficiencies will continue to cause substantial ESM differences. For this reason, the method of emergent constraints has attracted substantial attention. Such interest is because ECs offer a method to potentially 'short cut' current deficiencies in process understanding, by providing better estimates of change. Or, as a minimum, ECs offer more reliable estimates of bulk aggregated parameters of the climate system of interest to policy. We believe ECs do provide a route to reduce uncertainty, and consider it is appropriate for the technique to substantially underpin societal decisions both regarding adaptation planning and parallel mitigation programs to reduce GHG emissions. However this brings tremendous responsibility, and so the method must be reliable and robust. It is therefore prudent to discuss the multiple circumstances in which the method may fail, thus altering to where caution is needed. We now describe such situations. One approach with ECs is to simply 'data mine' within climate model ensembles, such as CMIP5 and CMIP6, searching for two quantities in each model that when analyzed form a statistically significant inter-ESM regression. In these circumstances, the researcher is not led by process understanding, a hypothesis to be tested, or indeed intuition (see section II.B for an extreme example). However, there is a logic that in some circumstances, this approach might be valid.
In the climate system, there may be links between different parts of the Earth system, and including between a contemporary measurable quantity and an attribute change in a future climate state, that are not immediately obvious as being connected. One set of examples could be similar to the on-going discovery of teleconnections in the climate system, where strong correlations are noted between variations in two parts that are a substantial geographical distance apart (Nigam and Baxter, 2015). Besides spatial distances, discovered teleconnections can contain lags, although unlike ECs that are designed to project decades ahead, the delays are months. These delays are often used to estimate expected meteorological conditions in the months ahead. For example, the status of El Nino, can strongly weight the probabilities of particular weather features for land regions and across the tropics and for the season ahead. Although El Nino is an obvious system fluctuation to investigate for its impacts elsewhere, over the last couple of decades, researchers have been investigating other candidates e.g. Feldstein (2000). Given the sometimes diverse form of teleconnections, their discovery can require a data mining approach, as comparison between two different parts of the climate system, additionally offset in space and time, are unlikely to be found through intuition-led inspection of ESMs. If novel ECs exist, between diverse parts of the Earth system, then such mining methods may be needed to aid their discovery. If ECs are discovered by statistical methods, then subsequent process analysis may reveal the connection to have a strong physical basis. Hall et al. (2019) refers to this as 'EC confirmation', although this may take many years to achieve, and especially if additional ESM diagnostics are needed to confirm the size of any transport equation terms between different spatial locations. However, whilst such an approach is expected to reveal new ECs, caution is needed. The concern is that statistically-significant regressions are likely to be found in a small number of instances, simply by chance, and the resultant EC is therefore not a robust indicator of future change (i.e. section II.B). That is, and presented loosely, if statistical mining finds a large number of new ECs, all with non-zero regressions forming them and at a 90% confidence level, then roughly one-in-ten will be invalid. This highlights the need to undertake process-based confirmation of ECs found through statistical means. The term 'p-hacking' was first introduced by Nuzzo (2014), and alerts to the risk of self-selecting only scientific findings that are statistically significant. In the context here, that would be ECs having their regressions at small p values, suggestive of a low probability of an X  Y relation occurring merely by chance. The concept has strong similarities to the dangers noted in data mining, where scanning across a range of potential ECs could lead investigators to only concentrate on a few with low p-values, with the attendant risk they occur by chance. However, p-hacking goes further and cautions over other decisions that researchers may make, in a rush to find a significant result. Behaviours could include a deliberate selection of one-sided tests and the stopping of sampling upon finding a low p-value. It could also involve adjusting the parameters, relationships between data, or time frequencies considered until lower p-values are discovered. In the case of the climate system, and as an example, the assessment of attributes of the global hydrological cycle offers many time scales of interest. Precipitation statistics with potential to form different X and Y quantities of ECs, range from short, intense rainfall events through to seasons or even decades. To account for overconfidence Bretherton and Caldwell (2020) suggests that ECs be corrected by scaling up the unexplained variance by a user-defined factor across the board. ECs most frequently link a fluctuating quantity that is measurable for the contemporary period, to either a future change (e.g. extent of polar ice sheet cover), or an invariant system attribute that describes change (e.g. climate sensitivity). A potential concern for the validity of an EC is if there is a missing process in every model that affects the regression forming the basis for the EC. If, for instance, that error affected the X-axis fluctuating quantity, but the measured quantity was accurate, then this would introduce a bias in to the projection of the future Y -axis quantity. Such a missing process could itself impact ESM performance at all modelled levels of atmospheric GHG concentrations. One potential example is that many ESMs underestimate interception loss, by vegetation, of rainfall e.g. Lian et al. (2018) and Yang et al. (2018). Interception loss is the return of water from rainfall to the atmosphere, and that has not passed through soils but instead is temporarily held on leaves and branches instead. If an emergent constraint on, for instance, a feature of the atmospheric part of the hydrological cycle that depended on the overall land-atmosphere exchange of water, then the EC could have an overly strong dependence on plant evaporation (i.e. to compensate for the low modelled values of interception loss). If, then, the true X value was taken from data of actual plant evaporation, this would introduce a bias in the estimate of the Y quantity. A straightforward corollary to this is if an EC depends on an X quantity that is modelled well, but measurements of it contain a bias, then this too would lead to inaccurate estimates of the value of the Y quantity. Alternatively, other process omissions could become more important as GHG levels rise. Until recently, an example of a component frequently not modelled in ESMs is the terrestrial nitrogen cycle (Thomas et al., 2015). Nitrogen limitation could become an increasingly important factor as CO 2 concentrations rise, by suppressing levels of vegetation fertilisation from the increased CO 2 in the atmosphere.
The lack of nitrogen cycle inclusion in the land components of ESMs could therefore affect ECs linking current variations of the carbon cycle to future carbon stores, leading to potential overestimates of the capability of terrestrial ecosystems to offset future anthropogenic CO 2 emissions. That said, in Wenzel et al. (2014), two ESMs that did have the nitrogen cycle included were close to an EC regression line in an ensemble of carbon cycle only ESMs. That analysis linked contemporary fluctuations in temperature and atmospheric CO 2 concentrations to the long-term stability of terrestrial carbon stores under global warming. A corollary to single model errors is the case where two contemporary processes are defined, but together they allow a model to perform well on some metrics as a consequence of compensating errors. There is quite substantial anecdotal evidence that many ESMs have these issues, as model developers often discover that improvement of one component, by removing the compensation, actually results in less accurate contemporary climate simulations. Poor performance is only then resolved when the process representation of the second feature is also upgraded. Such compensations present a risk for ECs. Specifically, if the X-axis quantity of an EC relates to just one of the two process that compensate, then that will introduce a bias. Returning to the interception example, if we consider the deficiency recognised but accounted for by parameterizing overly strong plant evaporation amounts, any EC dependent on plant evaporation (X) would introduce a bias in prediction of Y if actual measurements of X are used. Ultimately, climate science requires a set of ESMs that make similar projections and are highly accurate. Hence the on-going attempt to include all relevant processes in ESMs, and to remove compensating errors, remains an appropriate activity to advance such predictive tools. However, a co-benefit is that such ESM enhancements will likely raise the accuracy of any ECs that are used to constrain understanding where uncertainty remains. It is arguable that moderate human adjustment to the climate system is relatively linear, and so the expected meteorological changes increase proportionally with changes to the combined radiative forcing of atmospheric GHG concentrations. If this same linearity extends to higher frequency system responses such as interannual variability, then observations of contemporary fluctuations may give a strong indication of the system's response to a more permanent forcing from increased GHGs. One example is Cox et al. (2013), where the change to tropical atmospheric CO 2 concentration in response to annual temperature variation may project any expected loss of rainforest carbon stores under sustained global warming. Hence, there is an assumption of linearity present in each individual ESM, or at least a relationship that is monotonic in increasing temperature. When plotting the simultaneous variations in temperature and CO 2 (X) against response of future carbon stores in response to warming (Y ), then a further linear inter-ESM relationship is found, and this is the EC itself. A legitimate question, therefore, is what happens should the planetary system move beyond linear responses, and further, start to pass through tipping points (Lenton et al., 2008)? Many examples of such potential nonlinear behaviours are conjectured, where a relatively small increase in radiative forcing could cause major system changes. Potential climate tipping points such as the collapse of the Atlantic meridional circulation (Stocker and Wright, 1991), Amazon dieback (Cox et al., 2000), Greenland ice sheet melt (Toniazzo et al., 2004) and others (Drijfhout et al., 2015) are presently thought to be low probability, high risk events although there is evidence subsystems of the climate have 'tipped' many times in the past (Bond et al., 1992). Most of these 'tips' have resulted in regional rather than global climate changes and have been emulated through simpler equation systems, amenable to bifurcation analysis to characterise the nonlinearity (e.g. Dijkstra (2013)). Once past a large-scale tipping point, the Earth System might have radically different responses and feedbacks, and so behave very differently. This is likely to impact on the EC approach is two ways. First, for each individual ESM, for sufficiently changed climate, the contemporary fluctuating quantity may no longer give accurate information in that model on future store size. Second, if there are system tipping points, but an EC is only validated for moderate changes to the climate system (so, for instance, the Y -axis quantity is only derived from ESMs operated under lower emission scenarios), then using it to project to more severe climate alteration might be inappropriate. A more abstract but relevant question is whether information from pre-tipping point conditions has any validity at constraining climatic conditions beyond one? Do tipping points invalidate 'information exchange' between different GHG levels when passing through one? Two ESMs may have very similar responses up to a tipping point, and thereby appear as points near each other on a standard EC X  Y plot, but marked different responses beyond a tipping point, thus weakening an ECbased inter-ESM regression. Also in the abstract is a view that ECs are capitalising on the 'hard-wiring' implicit in ESMs, as based on model parameterizations constrained from contemporary measurement datasets. Beyond any major tipping point represents a world that has not yet been witnessed and so not measured for equation parameterization purposes. We hope that climate researchers can be encouraged to investigate further whether the planet passing through a tipping point invalidates ECs, and possibly also the period just beforehand. This could be in the form of a conceptual model, or by analysis of ESMs in the CMIP5 ensemble. ESM simulations with a tipping point could be split in to a future period pre-tipping point, and one post-tipping point for a high emissions scenario, testing the EC on both modelled time frames. One obstacle to undertaking such analysis is that although ESMs do exhibit tipping points, they vary markedly between models regarding their location, component of the Earth system affected, and level of global warming (Drijfhout et al., 2015).
Notable is that the strong regional features of tipping points implies that any ECs that may be affected by their presence also have to be localised in application. Possibly of more concern is that if tipping points affect the predictive capability of ECs, there will be little inter-ESM consensus on when and how this might occur. This lack of understanding is due to major model differences in predictions of tipping point occurrence, or even existence (Drijfhout et al., 2015), and so an 'EC-type' modulation factor to the original EC and common across models is unlikely to be discovered. ECs rely on the statistics of regressions, which in turn assume an independence of data points. However, individual ESMs may not be completely independent. For parts of the climate system known to respond to wellestablished and well-understood physical processes, commonality in models is to be expected. The concern is for the model parts that are suspected as characterising actual processes less well, and where ECs seek to constrain this uncertainty. The lack of independence may take the form of individual research centres offering multiple model versions but at different resolutions. Alternatively, some components of ESMs are shared between research centres, or are coded in similar ways (Knutti et al., 2013). The lack of independence could appear in the range of effective bulk parameters sampled, such as equilibrium climate sensitivity. Another possibility is there may be a common component existing between models, and that an EC seeks to reveal, but the ESM differences to reveal this are in the fluctuating forcings. These fluctuating quantities could have relatively low sampling due to similarities between how ESMs calculate their values. As part of developing comprehensive process understanding of discovered ECs, Hall et al. (2019) encourage 'out-of-sample' testing. Such testing is where the EC relationship is checked to be valid in additional ESMs that are not part of the original set used to initially find the EC. New simulations from individual modelling centres can be analysed as they become available. A more comprehensive test is to search within new ensembles, and so for instance an EC found in the CMIP5 climate model ensemble can be checked for its presence in the newly-released CMIP6 set of ESM simulations. However this too may not be a completely independent assessment, if modelling centres retain substantial amounts of model code and parameterizations between ensemble contributions. We note that the out-of-sampling testing performed by Schlund et al. (2020), that assesses if the ECs for climate sensitivity, ECS, found in the CMIP5 ensemble remain valid for the CMIP6 models. As commented elsewhere, they find that for CMIP6, the CMIP5-based ECs have less predictive capability, and also give generally higher ECS values. Schlund et al. (2020) note that the majority of emergent constraints for ECS are related to some extent to cloud feedbacks, and this is a major on-going area of climate research to create numerical cloud schemes with strong predictive skill. As many cloud schemes are currently under development, then their newness implies that this is an example of little carrying of common code from CMIP5 to CMIP6. From that perspective, CMIP6 models have independent features from CMIP5, validating using this more recent ensemble for out-of-sample testing. However, if it is the new cloud schemes that cause the CMIP6derived ECs to have less predictive capability, then this suggests that previously estimated EC-based uncertainties on ECS from the CMIP5 models are overly narrow. This form of error differs from that of Section IV.C that considers how a missing process in all models may create a systematic bias in an EC projection. Here, instead, the suspected missing processes associated with cloud dynamics and feedbacks are such that there likely remains substantial uncertainty in how to model them accurately. Alternative cloud schemes may weaken ECs by causing a larger spread around their regression lines, expanding the bounds on predictions of quantity Y for quantity X and for which contemporary data also exists. F. What to do when different ECs are found for the same quantity, but differ in value, or differ between ensembles? In some instances, multiple and different contemporary measurements have been suggested, via ECs, to be able to predict the same quantity. This is particularly the situation where a broad range of climate attributes have been used to evaluate equilibrium climate sensitivity (ECS). However, and notably, cloud feedback (Brient and Schneider, 2016;Brient et al., 2016;Sherwood et al., 2014;Zhai et al., 2015) or top-of-atmosphere (TOA) radiation flux (Brown and Caldeira, 2017) based emergent constraints tend to project higher ECS values than surface temperature based constraints (Cox et al., 2018a). Investigations need to continue to ascertain why these ECs disagree. Some disagreements may be due to spurious correlations, affecting Y values. However, Caldwell et al. (2018) show that these emergent relationships are highly correlated. With high correlation, but different ECS estimates, this implies that differences may instead be due to measurement biases in some X value observations. Alternatively, the real world may not be sharing the same responses as the models, which could suggest a persistent error across ESMs for at least one of the ECs. As ECs come under increasing scrutiny due to their growing widespread use, more needs to be understood about how they operate and any limitations. In that context, we suggest understanding the differences between projections of the same quantity will be a highly rewarding research path. However in order to generate a single EC-based projection, Bretherton and Caldwell (2020) present statistical methods to merge multiple ECs of the same quantity in to a single range of uncertainty for the quantity being estimated. That analysis focuses on combining the multiple ECs used by other researchers to estimate ECS. As noted in Section IV.E, new ensembles provide an 'out-of-sample' test for existing ECs. The question then is how to use that new value, and especially if it is substantially different to the value for an earlier ensemble.
One possibility is to use the methods of Bretherton and Caldwell (2020) to merge the findings from two ensembles. A further possibility is to introduce a weighting to the Bretherton and Caldwell (2020) approach, but with a high influence for the more recent ensemble, corresponding to the hope that newer ESMs are better models. In the most general terms, a reduction of spread of climate models is to be highly welcomed, and especially if their convergence is on to projections that are accurate estimates of future change. Such convergence is good if what is learnt from ECs applied to previous ensembles is accurate, and this has caused climate modellers to make their new simulations achieve the constrained Y values (see section VII.B). Convergence may also occur simply because previously uncertain parts of the Earth system have become better understood. However, the rapid development of ECs could cause a reduction in model spread that is potentially unwelcome, precluding further progress with the methodology. As Despite the aim of ECs being to reduce the effects of ESM uncertainty, we suggest that climate science remains alert to these two concerns, where the existence of an EC may overly reduce the model spread in any new ensemble (although this is not something that we have encountered to date). Such a reduction in a new ensemble, where model developers may wish to replicate the earlier EC, risks, in particular, the propagation of errors that might be prevalent in the earlier ensemble. In Section IV.A we heed the dangers of ECs based on data mining, noting that Hall et al. (2019) encourage that ECs discovered statistically should first move to having potential, i.e. contain an element of intuition as their confirmation. Detailed process understanding, maybe even as far as analytical assessment of related key climate-based differential equations, ideally then move an EC from having 'potential' to being 'verified'. At the same time, we suggest a potential danger whereby it is not possible to follow this cascade of increased certainly in validity of an EC. For example, there may be insufficient saved diagnostics from ESMs to build and study offline the process interactions that need understanding to achieve EC verification. The risk here is that important and informative potential ECs are dismissed as their understanding cannot be developed. Although this may ini-tially appear opposite to the suggestions in Section IV.A, we advise such ECs are not cast aside, thereby with the danger of losing valuable insights into the climate system. Instead such potential ECs from data mining can generate requests to climate research centres to provide related additional ESM outputs, to allow additional testing. Any additional ESM diagnostics may also reveal the extent to which a potential EC is dependent on one major physical process, or is caused by an amalgamation of multiple interacting effects. If more than one major process is present, then this may make it more difficult to transition from 'plausible' to 'verified'. The operation of ESMs is hugely computationally expensive and often there is only capacity in research centres to perform simulations for standardised scenarios of changes in future atmospheric greenhouse gas concentrations (e.g. the Representative Concentration Pathways, RCPs; Meinshausen et al. (2011)). Unfortunately, this precludes performing simulations with varied parameters, and especially for those where there is substantial uncertainty surrounding their true value (or effective value, if grid-dependent). We suggest that this lack of model investigation restricts further testing of ECs, and in a way that may either preclude gaining extra confidence in their reliability or alternatively fails to alert to one that may not be as robust as initially believed. We suggest that PPEs may require interpretation in two ways, and that are converse to each other. If a parameter value is changed, and if it relates to a part of the Earth system believed to be independent of the underlying processes associated with a particular emergent constraint, then for that ESM the Y value should not change in the EC. Such invariance both confirms EC robustness and helps illustrates that the EC aligns to the parts of the model expected from intuition or process understanding. The alternative situation is where a parameter believed to relate to an EC is perturbed, and potentially adjusted to the outer bounds of its expected value or beyond. The expectation then may be that for that particular ESM, the point moves to the bounds of the EC regression. Such an anomaly confirms that the EC is affected by the part of the climate system the parameter pertains to, and may indicate that outlier parameter values are not present in the other ESMs. If spare computational cycles become available, we suggest perturbed physics experiments with ESMs will be beneficial to the development and testing of ECs. Numerous ECs have been identified across physical and biogeochemical components of the Earth System (a handful of which have been discussed thus far), with a substantial increase in their number over the past decade. In order to capture the breadth of its application in the geosciences, we discuss previously documented ECs based on their broad geophysical classification here. A list of EC studies is provided in Table I. Note that this list is not exhaustive and may be easily superseded. Moreover, some of the references in Table I offer more than one potential current climate quantity to explain Y , in which case we list only one example. For example, Sherwood et al. (2014) provide three metrics that describe various aspects of lower-tropospheric mixing to explain equilibrium climate sensitivity. Rather than list all three metrics, we only discuss the one that appears most robust (Caldwell et al., 2018;Schlund et al., 2020). An important criterion for establishing credibility of an EC is to show robustness to the choice of ensemble. Therefore, Table I also documents the model ensemble(s) for which each EC appears to have value (i.e., a statistically significant emergent relationship is present).
In several cases this is only the ensemble from which the EC was derived as many of the more recent examples have not undergone out-of-sample verification at the time of this paper. We expect several upcoming publications to assess the validity of previously ) and Schlund et al. (2020)). A vast majority of ECs pertain to one of the following general topics: climate sensitivity, cloud feedbacks, cryospheric feedbacks and change, carbon cycle feedbacks, and the hydrologic cycle, but in theory the methodology can be used for any number of applications. There are also additional examples relating to topics such as radiative forcing (Bowman et al., 2013), regional air and sea surface temperature change (Lin et al., 2017;Selten et al., 2020;Sgubin et al., 2017), frequency of temperature extremes (Donat et al., 2018), and atmospheric circulation changes (Kidston and Gerber, 2010;Simpson and Polvani, 2016). For the purpose of this review, we will focus primarily on the main applications to date. Along with the differing applications, ECs can also be sorted by the type of constraint . One way to think of this is to group ECs by the time-scale of information that defines their X (e.g., from multi-decadal to extreme events). For example, a rather simple class of ECs relates climatological biases or multi-decadal trends in some quantity (X) to the future change in Y . Similarly, a collection of ECs use seasonal or interannual variations in X to constrain future changes in Y . In many of these instances, variability in some quantity to seasonal or interannual temperature variability is related to sensitivi- 1. ECs are sorted by the timescale of their current climate quantity and the type of future constrained quantity. ECs are also sorted by whether the future quantity relates to an impact or a physical or geochemical quantity. ties of that same quantity to future warming. Constraints using climatological biases or short-term sensitivities are also common for attempting to constrain invariant properties such as equilibrium climate sensitivity. Applications to future changes in higher-order climate statistics (e.g., extreme events) are more rare. This highlights a relatively unexplored area of EC research. By far the most popular application of the EC approach has been in relation to equilibrium climate sensitivity (ECS) (Table I). ECS is seemingly an ideal candidate for this technique because of its importance to predicting future warming rates and its consistently large spread across various GCMs, which has been a persistent feature of all model generations. (The likely range of ECS according to the IPCC (1.5 to 4.5K) has been largely unchanged for decades). A series of contemporary quantities (X) have been identified as being closely linked to ECS, many of which relate to simulated cloud feedbacks or properties. In general, these cloud-based constraints tend to suggest relatively high ECS values (3.5-4.1K) (Brient, 2020;Knutti et al., 2017). Relevant examples include the height of tropical low clouds (Brient et al., 2016), the sensitivity of the reflection by subtropical low clouds to sea-surface temperature (Brient and Schneider, 2016), seasonal sensitivity of low cloud between 20-40 deg N and 20-40 deg S to sea-surface temperature (Zhai et al., 2015), variability in relative humidity and cloud extent (Fasullo and Trenberth, 2012), vertically resolved zonally-average relative humidity and clouds between 45S-40N (Su et al., 2014), and the climatological difference between tropical and Southern Hemisphere mid-latitude cloud fraction (Volodin, 2008). Another suite of studies relate ECS to intermodel differences in the characteristics of historical radiative fluxes, such as climatological TOA energy balance in the Southern Hemisphere (Trenberth and Fasullo, 2010), features of TOA radiation fluxes (Brown and Caldeira, 2017;Huber et al., 2011;Tett et al., 2013), and cloud-sky radiative flux sensitivity to temperature (Lutsko and Takahashi, 2018). Alternatively, studies have also used temperature characteristics defined over a variety of timescales to constrain ECS. This approach tends to suggest relatively weaker ECS values than the others (Knutti et al., 2017). Examples of X include the seasonal cycle of temperature (Knutti et al., 2006), 20th century warming (Annan and Hargreaves, 2006), and statistics of interannual variability (Cox et al., 2018a). Similar to this, there is a collection of studies that have evaluated temperature changes over much longer timescales (e.g., global-mean cooling during the Last Glacial Maximum: 19-23 ka before present) as derived from paleoclimate data (Hargreaves et al., 2012;Schmidt et al., 2014). These paleoclimate-based constraints suggest ECS of 2.3 and 3.1 K, respectively. However, there is a large amount of uncertainty surrounding the model simulations of the LGM, the proxies used to constrain the intermodel relationship, and the validity of using data from these time periods as proxies for anthropogenic warming (Brient, 2020;Harrison et al., 2015). A variety of other metrics have also been proposed as being strongly related to ECS with large variability in their constrained predictands. The highest value (4.5K) is derived from a constraint using vertical mixing between the boundary layer and lower troposphere over tropical oceans (Sherwood et al., 2014). Tian (2015) use climatological precipitation in the Intertropical Convergence Zone (ITCZ) region to suggest a slightly weaker ECS (4.1K), while the weakest non-paleo estimate uses the climatological latitude of SH Hadley Cell edge in Dec-Feb to suggest an ECS of only 2.5K (Lipat et al., 2017). As discussed in Section IV.F, having such a large number of ECs for the same quantity can complicate their interpretation when the constraints suggest differing Y values. In this instance, rigorous testing of various proposed ECs is necessary to better understanding their validity. Along these lines, ECs on ECS have been heavily scrutinized over the past few years with assessments tied to their plausibility and robustness across both mul-timodel ensembles (Caldwell et al., 2018;Grise et al., 2015;Schlund et al., 2020) and, to a lesser extent, perturbed physics ensembles (Kamae et al., 2016;Wagman and Jackson, 2018;Zhao et al., 2016). Caldwell et al. (2018) provide the most robust assessment to date on the credibility of various constraints on ECS.
In doing so, they find that several of the aforementioned constraints lack a physically plausible mechanism connecting the X quantity with ECS (Siler et al., 2018;Su et al., 2014;Tian, 2015;Volodin, 2008). This limits how much faith can be put into the EC itself as the statistical relationship could be fortuitous. Furthermore, many of the constraints are found to be closely related (Caldwell et al., 2018). For example, one strongly correlated group relates aspects of present-day Southern Hemisphere cloud cover with ECS, suggesting they may all be capturing a single SH mechanism (Lipat et al., 2017;Siler et al., 2018;Trenberth and Fasullo, 2010;Volodin, 2008). Additional testing of EC robustness comes in the form of out-of-sample verification (i.e., evaluating the emergent relationship in a different model ensemble to the one from which the EC was originally derived). Although there are some known limitations to this approach related to the existence of common model code across generations (Section IV.E), it is still seen as a valuable exercise for evaluating ECs. For example, Grise et al. (2015) show that the Trenberth and Fasullo (2010) EC derived from CMIP3 is not valuable in CMIP5 because the emergent relationship only exists in a subset of models with unrealistic cloud properties in the Southern Hemisphere subtropics. Caldwell et al. (2018) extend this type of analysis to show that Fasullo and Trenberth (2012) is also non-existent in CMIP5. On the other hand, they find the Volodin (2008) EC to be present in both CMIP3 and CMIP5, despite lacking a physically plausible explanation. Schlund et al. (2020) recently assessed a number of the constraints on ECS documented here using the CMIP6 ensemble. They find that most of these constraints appear less skilful at predicting ECS in CMIP6, largely tied to differing representation of cloud processes. Three of the ECs discussed here (Brient and Schneider, 2016;Sherwood et al., 2014;Volodin, 2008) still exhibit statistically significant skill, while two others no longer have value (Lipat et al., 2017;Su et al., 2014). Lastly, two ECs are somewhat in between, what the authors define as "indeterminate", including Cox et al. (2018a) and Tian (2015) (Schlund et al., 2020). Further testing of these ECs has utilised perturbed physics ensembles to better sample the parametric uncertainty space associated with various emergent relationships (section II.A). Of particular interest to these studies is the EC proposed by Sherwood et al. (2014). Results from three different PPEs show that the physical mechanism proposed by this study may not be valid, while the emergent relationship between lower tropospheric mixing and ECS only exists when certain convection schemes are used (Kamae et al., 2016;Wagman and Jackson, 2018;Zhao et al., 2016). Similar doubts have also been reported for the Fasullo and Trenberth (2012) EC (Wagman and Jackson, 2018). This collection of analyses implies that the most robust constraints on ECS (determined by successful outof-sample testing and plausible physical mechanisms) are tied to the present-day response of subtropical low clouds to SST variability (Brient and Schneider, 2016;Zhai et al., 2015). However, one factor that may limit confidence in these constraints is the short observational record of low cloud characteristics (Brient, 2020). ECs from Sherwood et al. (2014) and Volodin (2008) are also statistically significant in three ensembles (Table I), but they face credibility issues tied to their physical mechanism (or lack thereof) (Caldwell et al., 2018;Wagman and Jackson, 2018). Despite these signs of robustness, the much weaker correlation values associated with these ECs in CMIP6 (Schlund et al., 2020) supports the belief that it is unlikely that a single predictor will be able to physically explain a large amount of ECS variability (Caldwell et al., 2018;Hall et al., 2019). Rather, several constraints on various feedback components will likely be needed to make a difference for such complex processes. In this fashion, a recent comprehensive assessment of ECS gave little weight to ECs on ECS when attempting to determine its likely value instead favouring observations of actual temperature change and climate feedbacks (Sherwood et al., 2020). This is also believed to be true for the global cloud feedback (described in more detail below) (Klein and Hall, 2015) given that the radiative effect of clouds in response to warming is expected to vary by cloud type and regime (Gettelman and Sherwood, 2016). Similar research efforts have recently suggested that historical decadal warming trends are strongly predictive of the transient climate response (TCR; defined as the amount of warming that occurs at the time of atmospheric carbon dioxide doubling, having increased by 1% each year) (Jimenez-de-la Cuesta and Mauritsen, 2019;Nijsse et al., 2020;Tokarska et al., 2020). TCR is of interest because it more closely resembles the way carbon dioxide concentrations have changed in the past. These constraints suggest that the best estimate of TCR is 1.67K in CMIP5 (Jimenez-de-la Cuesta and Mauritsen, 2019) and between 1.60K and 1.68K in CMIP6 (Nijsse et al., 2020;Tokarska et al., 2020). The difference between the latter two studies largely stems from the choice of historical time period evaluated (Nijsse et al., 2020). All three of these estimates are lower than the raw ensemble statistics (median of 1.95K in CMIP6). Similar constraints on future warming have also been developed for regional temperature changes. Both Lin et al. (2017) and Selten et al. (2020) find that present-day summer temperatures across central USA and Europe, respectively, are strongly tied to future continental summer warming. The authors hypothesise that the same physical mecha-nisms operating in the current climate (tied to precipitation biases and soil hydrology) also govern the future climate response. Incorporating observations allows for the future spread in European summer warming to be reduced by nearly 50% (Selten et al., 2020). One of the primary contributors to intermodel spread in climate sensitivity is the cloud feedback, which encompasses changes in how clouds modulate the radiation budget in response to warming. (For a recent synthesis of cloud feedbacks see Section 3.3 of Sherwood et al. (2020)).
The global mean cloud feedback is the most uncertain component of the total climate feedback (Bony et al., 2006;Ceppi et al., 2017;Sherwood et al., 2020), varying from -0.36 to 1.17 W m 2 K 1 in the latest generation of climate models (Zelinka et al., 2020). Strengthening of ensemble mean cloud feedback strength in CMIP6 is also thought to be largely responsible for a recent ECS increase in many models (Zelinka et al., 2020). Thus, many of the aforementioned constraints for ECS include cloud effects either directly or indirectly (Brient and Schneider, 2016;Brient et al., 2016;Fasullo and Trenberth, 2012;Qu et al., 2018;Su et al., 2014;Volodin, 2008;Zhai et al., 2015). This also means that there is great interest in constraining cloud feedbacks as well (Brient and Bony, 2013;Gordon and Klein, 2014;Qu et al., 2014;Siler et al., 2018;Zhou et al., 2015). For example, Siler et al. (2018) propose a constraint on future global mean cloud feedback using the climatological latitudinal gradient in cloud reflectivity (a quantity that can be derived from satellite observations). In general, this EC suggests that ESMs with lower cloud albedo in warm-SST regions and higher cloud albedo in cool-SST regions will exhibit greater cloud feedback (and to a lesser extent ECS). It suggests a global cloud feedback in the upper range of ESM estimates (0.58 +/-0.31 W m 2 K 1 ). However, this constraint lacks a testable physical explanation (Caldwell et al., 2018). Another study that seeks to constrain the global mean cloud feedback is Zhou et al. (2015). This paper highlights strong similarities between the global mean cloud feedback and a present-day analogue cloud feedback derived from interannual climate variability (Zhou et al., 2015). This follows the assumption of several ECs that processes operating on seasonal to interannual timescales are also likely to be occurring on longer timescales. However, the usefulness of this EC is limited because the observational uncertainty in the present-day quantity is large relative to inter-model spread. As previously mentioned, it is seen as unlikely that a single current climate quantity can account fully explain variability in the global cloud feedback (Klein and Hall, 2015). Thus, a number of studies have focused on iso-lating the individual components of the cloud feedback (high-cloud altitude, tropical marine low-cloud, tropical anvil cloud-area, land cloud amount, mid-latitude marine low cloud amount, and high-latitude low-cloud optical depth). For example, Po-Chedley et al. (2019) narrow in on the changes to upper troposphere clouds and relative humidity across the tropics. They show that future changes in upper troposphere cloud fraction under climate warming are strongly tied to the vertical gradient in climatological mean cloud fraction and relative humidity. Elsewhere, much of this research concerns tropical low clouds as ESMs struggle to represent them, they account for nearly half of the variance in the global mean cloud feedback (Klein et al., 2017;Zelinka et al., 2016). Moreover, intermodel spread in the low cloud feedback is closely tied to variability in ECS (Sherwood et al., 2020). An early EC along these lines showed that the future low cloud feedback sign (Y ) is related to current low cloud sensitivity as a result of natural variability (X) in CMIP3 (Clement et al., 2009). Unlike most other ECs, this example simply suggests the sign of a feedback rather than specific values. Thus, it is not always mentioned in reviews of this nature (Brient, 2020;Klein and Hall, 2015). More traditional ECs relating to this topic include the link between low cloud optical depth changes with warming and the optical depth response to temperature anomalies associated with natural variability (Gordon and Klein, 2014). This feedback is of limited importance globally, but makes an important contribution at high latitudes (Klein and Hall, 2015). The Gordon and Klein (2014) EC also has a simple thermodynamic explanation, which makes it a promising example. Low cloud amount changes under future warming play a much larger role in the global cloud feedback so this has naturally been an area of focus (Brient and Schneider, 2016;Qu et al., 2015Qu et al., , 2014)). Qu et al. (2014) identified an EC on this low cloud response through the low cloud response to temperature and stability anomalies derived from interannual variability. Although weaker in CMIP5 than CMIP3 (Klein and Hall, 2015), this EC has value in both ensembles and features a strong physical basis explaining the relationship (Qu et al., 2015(Qu et al., , 2014)). This type of process based evidence can then be used in combination to assemble a best estimate of the more complex global mean cloud feedback. Figure 7 of Sherwood et al. (2020) illustrates this concept to get a global cloud feedback of 0.45 +/-0.33 W m 2 K 1 . Going forward, it is imperative that we continue to improve our process-based understanding of individual cloud components, where potential physical mechanisms are more likely to be uncovered. Elsewhere, the EC technique has been used extensively to constrain elements of climate-carbon cycle feed-backs (Table I #28-36) (Cox et al., 2013;Hoffman et al., 2014;Mystakidis et al., 2017;Wenzel et al., 2014Wenzel et al., , 2016)). Climate-carbon cycle feedbacks are characterised by changes in terrestrial and ocean carbon storage in response to climate change, which modify the atmospheric CO 2 concentration, thus enhancing warming. Intermodel differences in feedback strengths lead to highly variable projections of future atmospheric carbon dioxide (more than 300 ppmv of intermodel spread at the end of the 21st century for a high-emissions scenario (Cox, 2019)). A rather simple relationship was uncovered between the mid-century atmospheric CO 2 concentration in ESMs and the simulated present-day concentration, but it was shown to weaken as the influence of varying climatecarbon feedbacks becomes more apparent later in the 21st century (Hoffman et al., 2014). This EC suggests a much narrower likely range for atmospheric CO 2 concentrations at the end of the century (947 +/-35 ppm) than the full ensemble (980 +/-161 ppm).
Much of the model uncertainty in future atmospheric CO 2 has been shown to stem from differences in the land carbon sink (Cox et al., 2013;Wenzel et al., 2014). Therefore, significant effort has gone into constraining this feature of the carbon cycle. A prominent example uses the historical sensitivity of annual atmospheric carbon dioxide growth-rate to temperature variability (which is strongly impacted by tropical land carbon storage fluctuations) to constrain the future loss of terrestrial tropical carbon (Cox et al., 2013;Wenzel et al., 2014). As with prior examples, the CO 2 growth rate sensitivity is an observable quantity, thus allowing for the EC to suggest future tropical land carbon stability. This EC has been proven robust across both C4MIP and CMIP5 ensembles (Wenzel et al., 2014) 5 , but has yet to undergo testing pertaining to CMIP6. This is an important EC because there was a fourfold simulated spread in tropical land carbon reductions per degree of warming in C 4 MIP (29-133 GtC K 1 ) (Jones et al., 2016), with the higher-end models suggesting potentially catastrophic dieback of the Amazon rainforest (Cox et al., 2004(Cox et al., , 2000)). However, the Cox et al. (2013) EC suggests tropical land carbon reductions slightly weaker than the ensemble mean (53+/-17 GtC K 1 ). Another carbon cycle EC uses the change in seasonal atmospheric CO 2 amplitude to constrain CO 2 fertilisation of photosynthesis on the extratropics (Wenzel et al., 2016). On a more regional basis, Winkler et al. (2019a) highlight a strong relationship between future increases in terrestrial Arctic gross primary productivity and historical increases in leaf area index (greening). These two ECs both suggest that most models are underestimating future changes in gross primary productivity across the high latitudes (Wenzel et al., 2016;Winkler et al., 2019a). While these ECs have yet to undergo out-of-sample verification, Winkler et al. (2019b) provides detailed evaluation of the factors contributing to uncertainty in the (Winkler et al., 2019a) constraint. They point to choices relating to the temporal period of the predictor variable, choice of observational dataset, and the rate of CO2 forcing as particularly large sources of uncertainty (Winkler et al., 2019b). Lastly, additional research identified constraints on terrestrial carbon cycle feedbacks using interannual variability in evapotranspiration, net biome productivity, and gross primary productivity. This EC suggests a 40% reduction in the climate-carbon feedback and a 30% reduction in the concentration-carbon feedback (Mystakidis et al., 2017). When it comes to understanding how the ocean carbon cycle will respond to climate change, one of the main challenges is to reduce uncertainty in estimates of tropical ocean primary productivity. Satellite observations of the interannual variability in ocean productivity resulting from ENSO-driven SST anomalies have been used to constrain highly uncertain projected changes in tropical marine primary productivity with warming (Kwiatkowski et al., 2017). There is a strong inverse relationship between NPP and SST anomalies in the observational record, but the number of ENSO events in the satellite record is limited. This EC suggests a substantial reduction in the long-term tropical NPP sensitivity (-3.4% K 1 to -2.4 % K 1 ) from the unconstrained ensemble average (-4.0 +/-2.2 % K 1 ). Another recently ) relates future global ocean carbon uptake to the contemporary carbon uptake by the Southern Ocean. In this example, models with anomalously low uptake in the current climate project low global uptake over the course of the 21st century. Uncertainty in simulated ocean uptake also translates to uncertainty in projections of future ocean acidification, the impacts of which are expected to be greatest in the Arctic ecosystem (Terhaar et al., 2020). Recent work uses CMIP5 models to identify an EC on future Arctic ocean acidification using the simulated density of Arctic ocean surface waters in the current climate (Terhaar et al., 2020). Observations of sea surface density when coupled with this emergent relationship imply that future Arctic Ocean acidification will be greater than previously expected. The Earth's high-latitude regions are rapidly warming and there is a great deal of uncertainty in how various components of the cryosphere (Earth's frozen surfaces) will evolve in the future (Meredith et al., 2019;Mudryk et al., 2020;Notz et al., 2020). Thus, ECs are an intriguing option for reducing intermodel spread here. One of the earliest documented examples of an emergent constraint is for the snow albedo feedback (Hall and Qu, 2006), which was introduced in Section I. This is a leading example of the EC technique because of the simplicity behind its physical mechanism and its robustness across several generations of models . An interesting point on the SAF EC is that despite substantial time since its initial publication, we have yet to see a meaningful reduction in its current climate quantity (X). Following on from its application to SAF, recent work has identified that the emergent constraint approach can also be applied to a similar process over Arctic sea ice (Thackeray and Hall, 2019). The sea ice albedo feedback (SIAF) is characterized by the enhancement of future warming through the reduction of surface albedo as a result of decreasing sea ice (Hall, 2004;Holland and Bitz, 2003). An analogue of this process operates in the seasonal cycle during the transition from maximum to minimum ice extent (X). These two processes are found to be closely linked across the CMIP5 ensemble, forming the basis of an EC. However, given the projected rapid loss of sea ice in the 21st century, this constraint exhibits a regime-dependence whereby the relationship between seasonal cycle and climate change SIAF begins to break down in the latter half of the 21st century, when most models exhibit an ice-free Arctic (Thackeray and Hall, 2019). This is similar to the concept of tipping points discussed in Section IV.D. In addition to the aforementioned radiative feedbacks, there is a large amount of uncertainty when it comes to the fate of a number of high-latitude Earth system components, notably sea ice and permafrost (Slater and Lawrence, 2013;Stroeve et al., 2012).
Given this uncertainty and the potential widespread implications associated with these changes, ECs can provide valuable information here. Notably, a number of studies have proposed constraints on projections of future Arctic sea ice (Boe et al., 2009;Liu et al., 2013;Massonnet et al., 2012). Boe et al. (2009) used the observed historical trend in September Arctic sea ice extent over the satellite era to constrain the time in the 21st century when the Arctic is likely to become ice-free during summer. Similar research by Massonnet et al. (2012) used a series of historical ice characteristics and trends in an attempt to constrain the time period when the Arctic is likely to become seasonally ice-free. Their results suggest ice-free conditions for September to begin somewhere between 2041-2068 under a high-emissions scenario, a significant reduction from the full CMIP5 spread that spans nearly 100 years. The presence of constraints between historical ice properties and the future change in both CMIP3 and CMIP5 is a promising sign, but we are not aware of CMIP6 testing for this EC thus far. Similar historical sea-ice metrics have also been used to inform future changes in high-latitude temperature variability (Borodina et al., 2017a). Permafrost underlies a significant portion of the Northern Hemisphere high-latitudes and is highly sensitive to warming (Lawrence et al., 2012), with its degradation expected to have detrimental climate and developmental effects (Teufel and Sushama, 2019). However, modern ESMs struggle to agree upon future changes in permafrost area (Slater and Lawrence, 2013). Chadburn et al. (2017) use a relationship between mean annual air temperature and permafrost area during the historical period (X) to constrain projections of future permafrost thaw (Y ). This constraint suggests that roughly 20% more permafrost will be lost per degree of warming than previously expected (Chadburn et al., 2017). Future changes to precipitation are of great interest to the climate science community because of wide-ranging impacts on natural and human systems. However, precipitation is highly variable in space and time, thus making it difficult to both observe and predict. These factors contribute to large intermodel differences in future projections of precipitation. Reducing this uncertainty is therefore vital for those making policy and infrastructure decisions based on model projections. ECs have been applied to multiple aspects of the hydrologic cycle across various scales (DeAngelis et al., 2015;Li et al., 2017;O'Gorman, 2012;Rowell, 2019;Watanabe et al., 2018). First, we consider the globally-averaged changes to the hydrologic cycle. In response to rising temperatures and increasing atmospheric water vapour, globalmean precipitation is projected to increase. This is because the atmosphere radiatively cools to space and this radiative cooling must be balanced by latent heat release from precipitation, thus setting the radiative-convective balance (DeAngelis et al., 2015;Pendergrass and Hartmann, 2014). Under a high-emission scenario, models project a future increase of between 2 to 10 % in globalmean precipitation (Kharin et al., 2013). This spread can be partially explained by differing rates of future warming (Flaschner et al., 2016), although a threefold spread still exists when the precipitation change is normalised by warming (Kharin et al., 2013). In an attempt to minimise this spread, DeAngelis et al. (2015) used the sensitivity of clear-sky shortwave radiation absorption to changes in column water vapour to constrain the global mean precipitation increase with warming (termed hydrologic sensitivity). This constraint suggests a 35 % reduction in the ensemble spread and a 40 % reduction in the ensemble mean hydrologic sensitivity. This work was followed by Watanabe et al. (2018), who use the climatological global mean surface longwave cloud radiative effect (CRE) to constrain the future surface longwave CRE feedback. In combination with the most-likely sensitivity of clear-sky shortwave radiative absorption to water vapour (DeAngelis et al., 2015), this EC suggests a constrained estimate of 1.8 % K 1 for global hydrologic sensitivity, down from an ensemble mean of 2.6 % K 1 (Watanabe et al., 2018). However, a recently ) EC is still present, its strength is weakened (Pendergrass, 2020). Moreover, the Watanabe et al. (2018) EC on HS is found to be nonexistent in the new ensemble. This highlights the importance of applying out-of-sample testing to the remainder of the ECs on hydrologic cycle components discussed below. Another key component of the global hydrologic cycle is terrestrial evapotranspiration (ET), which is made up of evaporation and biological transpiration (T). The ratio of T to ET is used to estimate land water fluxes, while because similar processes influence T and land-atmosphere carbon exchanges, better estimates of T/ET help to reduce uncertainty in global carbon cycle projections (Lian et al., 2018). Climate models exhibit large variability in the global strength of T/ET, spanning a factor of 3. Lian et al. (2018) found a strong relationship (r 2 = 0.93) between simulated local T/ET (averaged over grid-cells corresponding to in situ measurements) and global mean T/ET. This EC reveals that the constrained global T/ET ratio is significantly higher than the CMIP5 ensemble mean. Alternatively, several studies have attempted to use this technique to constrain the change in aspects of the hydrologic cycle over specific regions. Regional constraints have the potential for more policy and impacts relevance than globally averaged metrics. As a first attempt, O'Gorman (2012) proposed an EC on the change in future daily tropical precipitation extremes using the observed sensitivity of tropical precipitation extremes to temperature variability. They found a strong relationship between the sensitivity of very heavy precipitation (99.9th percentile) to both seasonal and future temperature changes. The constrained estimate based on observations (6-14% K 1 ) was substantially narrower than the intermodel spread across the CMIP3 ensemble (2-23% K 1 ). Similarly, Borodina et al. (2017b) found that the scaling of annual maximum one day precipitation (Rx1day) with warming global land temperatures can be used to constrain the future intensification of heavy rainfall across extratropical regions with sufficient data records. This EC suggests that the CMIP5 models are likely underestimating the future change in extreme rainfall where climatological rainfall intensity is high.
Another example links biases in simulated historical climate and future Indian Summer Monsoon (ISM) rainfall (Li et al., 2017). In this case, the authors find that models overestimating historical precipitation across the tropical western Pacific also exhibit larger increases in Monsoon rains in the future (Li et al., 2017). Using observations of western Pacific rainfall reduces the projected increase in ISM rainfall by nearly 50%. CMIP5 models also differ greatly in their projections of future total (March-April-May mean spanning a region from northwest Tanzania to southwest Ethiopia) -from a 20% decrease to a 120% increase (Rowell and Chadwick, 2018). Rowell (2019) use an EC between present day interannual SST-low level cloud sensitivity and future SST change over the Indian Ocean to assess the credibility of projections of East African precipitation. They find that one outlier model, which projects a doubling of seasonal precipitation is likely unreliable because of unrealistic SST/low cloud processes. This reduction in intermodel spread is another example of the value that an EC approach can provide. Along these lines, Lehner et al. (2019) recently suggested that the sensitivity of historical runoff to temperature and precipitation change across three watersheds in the Western US is closely tied to projections of future runoff across a series of ESMs. It is expected that regional applications of ECs to highly uncertain quantities like precipitation change will be a key area of future research. TABLE I: Collection of existing ECs. Note that some of these ECs involve correlations that are lower than those portrayed in Figure 1, with correspondingly less potential for uncertainty reduction. The ensemble(s) for which the EC appears to have value is also listed (note that many have only been tested on the ensemble that they were developed on). Accurately constraining the unknown future value of Y in the real world requires one to include and quan-tify all the possible sources of uncertainty in each step of the EC procedure. This section will give an overview of four types of uncertainty that have been incorporated into ECs so far: those stemming from uncertainty in the real world observation of X, uncertainty in X from internal variability, uncertainty in the functional form of the emergent relationship and uncertainty from ESMs being imperfect replications of the real world. The section concludes with a discussion on how to combine those four types of uncertainty in the resulting in EC. In addition, we discuss how to combine multiple emergent constraints of the same quantity, derived from alternative features of historical climate. There are multiple sources of observational error in the real world value of X. First, a lack of spatial and/or temporal coverage can be present and this may lead to biases if not taken into account (Cowtan and Way, 2014). There are two ways to handle missing data: it can be interpolated and extrapolated from existing data, or alternatively, model output can be filtered to reflect only locations and times for which observational data is present (e.g. AchutaRao et al. (2006); Cox et al. (2018a); Durack et al. (2014)). Secondly, observational records are of finite length, introducing additional uncertainty from low sample size. Standard errors quantifying finite size effects can be computed, but care should be taken when timeseries are autocorrelated as this increases the standard error by effectively reducing the sample size (Trenberth, 1984). For stationary processes, standard equations for autocorrelation errors can be found in Zhang (2006). Estimates of errors in instrumentation and data gathering are often available from literature e.g. Hennermann (last edited July 2018). If multiple observational data sets are available, these can be used to infer uncertainty (Kwiatkowski et al., 2017;Trenberth and Fasullo, 2010). The errors in these data sets might not be independent; different satellite products might for instance have the same biases. Independent sources of error  1 and  2 can simply be computed by  2 total =  2 1 +  2 2 . If observational uncertainty makes up a large percentage of overall uncertainty, care should be taken to assess whether to use a normal distribution to describe the probability density function (pdf). It may sometimes be possible to estimate a full pdf from measurements. Alternatively, stochastic reduced-form modelling of the system can be used to estimate the shape (Nijsse and Dijkstra, 2018;Williamson et al., 2019). Whatever method is used, it is important that this step of capturing observational error is not neglected (Hall et al., 2019). Like the real world, climate models have internal variability. Because of the finite length of the simulation or observed climate record, internal variability can have a significant impact on the estimation of the predictor. One possibility is to use very long model control simulations to estimate the size of internal variability if it is believed to be independent of forcing (Nijsse et al., 2019). Variability may however be dependent on forcing, and consequently, estimating it from a forced initial value ensemble (see section II.A) may be preferable (Tokarska et al., 2020). For instance, global inter-annual variability is expected to decrease in the future (Huntingford et al., 2013). Jimenez-de-la Cuesta and Mauritsen (2019) used the 100-member historical ensemble of MPI-ESM1.1 to quantify the effect of internal variability, whereas Nijsse et al. (2020) used all available historical initial value members from each CMIP6 model to estimate the mean model variability. Both used model estimates as a proxy for real internal variability. Reducing a high-dimensional climate model to a lower dimension brings some uncertainty. Clearly, not all the variance will be explained with only two variables, and performing a regression is a tool to quantify this. While most emergent constraints so far assume linear relationships between X and Y and use linear regression to infer the emergent relationship, the regression does not necessarily have to be linear (Bracegirdle and Stephenson, 2012;Nijsse and Dijkstra, 2018). If a linear relationship is imposed when, in reality, the relationship is nonlinear, additional errors will occur.
A nonlinear emergent relationship leads mostly to a non-normal pdf for Y with a standard deviation that is potentially significantly larger or smaller than for a linear fit . Of course, more data or clear prior information on the parameters is needed when fitting additional parameters. Regression dilution takes place as a consequence of errors in estimating model predictors (e.g. finite simulation length): when there is an error in the modelled explanatory variables, the slope of a linear regression fit will be smaller than without error and the intercept regresses towards the mean (Frost and Thompson, 2000). Multiple strategies to reduce this can be employed: taking the mean of a set of initial value simulations (Jimenez-de-la Cuesta and Mauritsen, 2019), using orthogonal distance regression, which takes into account both errors in the dependent and independent variable (Jimenez-de-la Cuesta and Mauritsen, 2019) or using a hierarchical Bayesian method that assume the 'true' independent variable is unknown (latent) and using the realizations to infer this true value simultaneously while performing the regression Illustration with synthetic data. (Nijsse et al., 2020;Sansom, 2014). Regression confidence intervals also depend on sample size. As climate models often have shared computer code, effective sample size is probably lower than the total number of models (Herger et al., 2018;Masson and Knutti, 2012;Pennell and Reichler, 2010), effectively reducing significance (Knutti et al., 2013). This commonality can be partially addressed by selecting only one model per modelling centre (Cox et al., 2018a;Sansom et al., 2019), where centres frequently offer multiple versions of alternative spatial resolutions or other small differences in model physics. However, even then effective sample size may be overestimated. As strongly related models in terms of code may have very different values for X and Y , it is not always clear how important corrections are (Nijsse et al., 2020). Almost paradoxically, the technique of emergent constraints relies on errors in models. Hence some type of errors are useful, while others simply contribute to widening the confidence interval. Model error comes from parameteric and from structural errors (see also section II.A). Parameter error is often examined by perturbed physics ensembles (PPEs), an ensemble where a single climate model is run repeatedly with systematically varied parameters designed to span a large range of model responses. Care should be taken to only include the physically adequate parameter choices, but not restrict the parameter space too much (Wagman and Jackson, 2018). Usually, due to the expense of running climate models, an ad hoc ensemble of different models is exploited to establish emergent constraints, the so called 'ensemble of opportunity'. In this case, models should be cho-sen in such a way that they are still comparable: If half the models contain a nitrogen cycle, while the other half does not, they are unlikely to fit on a single regression (Wenzel et al., 2014). Williamson and Sansom (2019) argue that due to structural uncertainty, the derived regression between X and Y should not be regarded as the real relationship, but instead as only informing the real relationship. E. Combining sources of uncertainty in an EC Williamson and Sansom (2019) argue on theoretical grounds that emergent constraints should be performed using a Bayesian framework, instead of the more common frequentists' framework using ordinary least square fitting (e.g. Wenzel et al. (2014)), in which priors are set for the regression parameters but not for the predictand Y . Bowman et al. (2018) formulated a Bayesian hierarchical statistical framework that combines uncertainty in observations with uncertainty related to the regression itself. An explicit formula for the probability distribution can be approximated incorporating these uncertainties: assuming normality and performing a linear regression with intercept a, slope b and confidence around the regression of  f , the probability density function for the predictand p(y) is computed by integrating the conditional probability function p(y|x) with the observational probability density p(x). This formula for p(y), rewritten in the notation of other statistical frameworks and under Gaussian assumptions, turns into Here X is the predictor, Y the predictand, and X obs is the best estimate of the observation. The observational uncertainty is denoted by  O . To better estimate the regression and incorporate internal climate variability explicitly, a second hierarchical model was developed independently by Sansom (2014) and Nijsse et al. (2020). Internal climate variability is incorporated both as part of the regression, and as an additional term in the uncertainty around observations, by using all initial value simulations X m,j of each model m and each initial value member j. Here formulated using Gaussian distributions: The internal variability is denoted by  I and is assumed to be independent of Y and X (and therefore model independent). The probability density function is sampled from estimates of the observed X obs and the parameters of the emergent relationship, giving a similar equation to Eq. 20. p(y) = N a + bN X obs ,  2 I +  2 O ,  f Hall et al. (2019) queries the extent to which we should trust in only partially confirmed ECs. This can be addressed formally, by inclusion of a researcher's assessment of how reliable a certain EC is. Williamson and Sansom (2019) describe a method explicitly incorporating the trust that is put in any particular EC, including for instance by how much the EC changes in different ensembles. They introduce an additional parameter  R that represents the uncertainty even after having taken all model information into account. This parameter is set subjectively, judged on the degree of physical trust the researcher has in the EC and features as p(y) = N (x,  2 f +  2 R ), where  is a covariance matrix of the regression parameters, for which priors are also provided with information about possible biases or errors in the regression parameters. A very simple method to combine multiple constraints was used in Brient (2020). They use a Gaussian kernel density estimation of a histogram of the best values of previous constraints.
To account for the variance, the posterior variance of the different ECs was included in a weighted Gaussian kernel density estimation. This method suffers from multiple drawbacks. Most prominently it cannot capture increasing confidence from having multiple independent constraints. Adding a new constraint with the same mean does not automatically lead to a narrower combined constraint, with the choice of kernel bandwidth remaining subjective. It further does not take into account to what extent the different emergent constraints are related to each other. In Bretherton and Caldwell (2020) multiple emergent constraints were combined using a multivariate Gaussian PDF, which can be viewed as a form of multilinear regression. Their 'method C' (for Correlated) includes information about correlations between different emergent relationships. Regularization was applied to deal with the strong collinearity between the emergent relationships. A second 'method U' (for Uncorrelated) used a smaller subset of ECs, those regarded as confirmed constraints. As collinearity is less important with fewer ECs this method simplifies C by dropping the covariance ECs. Both variants were extended with a transparent method to account for overconfidence in the EC: they scale the ratio of the explained to unexplained variance with a factor  2  1, reducing all correlation coefficients. where gaps exist, thus highlighting where the potential for EC growth is greatest going forward. For the most part, ECs have been traditionally focused on globallyaggregated quantities related to the mean state of climate (e.g., climate feedbacks, hydrologic sensitivity). However, recent applications to more regional features of the climate system are particularly encouraging (e.g. informing future monsoonal precipitation change over India; Li et al. (2017)). Seeing as these regional applications have the potential for greater policy impact than most globally-averaged metrics we anticipate much more focus in this area going forwards. Climate extremes are also of great societal importance, but from Figure 5 we can see that very few studies have attempted to constrain the uncertainty in their future changes. We anticipate a greater focus on these and other higher order moments of climate statistics going forward. This advance is spurred on by an ever-lengthening satellite record that is now long enough to robustly sample extreme precipitation events. Similarly, improving paleoclimate records offer great potential for constraining emergent relationships pertaining to aspects of the climate system that vary on long timescales. In terms of their Earth system components, proposed constraints are heavily skewed towards atmospheric and biogeochemical fields, while constraints pertaining to oceanography and ice sheets are lacking. The latter is a rather recent addition to ESMs, but as observational records continue to improve there could be great potential for reducing uncertainty in critical metrics like sea level rise. Lastly, most ECs pertain to local relationships, but there are likely many undiscovered remote relationships in the climate system, where the current variability in one region is strongly tied to the future change in another region through atmospheric teleconnections (Rowell, 2019). It is also believed that biases in the position of various climate features may be systematically tied to how the features respond to future forcing (Hall et al., 2019). Utilising these spatial characteristics of climate for ECs has begun (Kidston and Gerber, 2010;Simpson and Polvani, 2016), but better detection of these remote constraints likely requires improved analysis techniques e.g. Barnes et al. (2019). An appealing yet underused aspect of ECs is targeted model development for bias reduction in a particular aspect of climate change (Y ), although there are dangers to be aware of too (see section IV.G). If an EC has strong physical underpinnings, then we can expect that taking steps to reduce the spread in X will result in a corresponding spread reduction in Y . First, we must assess in detail how the structural and parametric makeup of the ESMs influence the spread in X. Through this analysis the best parameterizations or parameters for simulating X may be uncovered, thus providing guidelines for modelling centres on how they can efficiently reduce bias in X going forward. This type of analysis is encouraged for all well-established ECs, but has only been completed for ECs on hydrologic cycle intensification (DeAngelis et al., 2015), the sensitivity of extratropical cloud reflectivity to temperature (Gordon and Klein, 2014), and snow albedo feedback (Thackeray et al., 2018). It is also possible that these targeted development activities will lead to a reduction in spread for related attributes of climate change affected by Y . C. Use of conceptual models as the basis of emergent relationships and understanding of more complex ESMs In the earlier years of ECs there may have been a tendency to data mine for variables with high correlations within ESM ensembles, potentially resulting in bad ECs (see section II.B). However, data mining is not necessarily bad if it guides the search for independently testable mechanisms that may not have been obvious otherwise (section IV.A). One way of guarding against fortuitous correlations and hypothesizing new mechanisms is using conceptual, analytically soluble models as the basis of the emergent relationship. Building and solving conceptual models requires one to make assumptions about the real world. These assumptions are often clear, simple to understand and ideally, testable against real world and ESM ensemble data (Williamson et al., 2018). Provided the assumptions survive these tests, the solution to the conceptual model can form the basis of the emergent relationship. This testing of simplified theory not only aids understanding and intuition of what might (or might not) be going on in the complex ESMs but also adds a greater degree of rigour to ECs. There are potential dangers one should be aware of however; (i) p-hacking might occur (see section IV.B) and (ii) it might be that a process in the ESM one is trying to constraint is closer to the conceptual model in the spectrum of complexity than it is to reality. The use of conceptual models could therefore be used to justify the relationships arising from slightly less conceptual parameterizations which are implemented in the ESMs.
State-of-the-art ESM spatial resolution has continually increased and more physical and biogeochemical processes have been incorporated making these models ever more comprehensive representations of the Earth as time has gone on. They have become the main theoretical tool in climate science. However, due to their ever increasing complexity, their behaviour is difficult to understand and predict without running the models themselves. They are also frequently treated as one-stop shops and oracles for any question about contemporary and future climate. In the past, numerical simulation was more limited by computing power and researchers had to think carefully about which bits of the Earth system were important to include in their models to answer a particular question. This led to wider use of conceptual models specific to that question, analytical techniques and a better intuition due to the reduction of model complexity (Budyko, 1969;Dijkstra, 2013;Hasselmann, 1976;Saltzman, 2002;Wigley and Raper, 1990). Today, there is even more of an opportunity for the top down insights of specific conceptual models to meet and complement the comprehensive, bottom up approach of from state-of-the-art ESMs; there are many more high quality observations from satellites (Yang et al., 2013), ocean floats (Roemmich et al., 2019), and continuous temperature records; the global warming signal also has become clearer through the passage of time; there is also a large archive of past and present ESM simulations. With all the new data, it would be worth revisiting these simple, understandable models. We think a promising way forward to incorporate all the extra information from ESMs, conceptual models and observations is the EC approach. Much of the work on ECs has used just two scalar variables, a predictor X, and a response or predictand Y , related linearly. This could be extended to two or more predictors or predictands (Renoult et al., 2020). The predictand space to constrain will become larger however, and this might be difficult with small ESM ensembles. Most of the work to date has also assumed a linear emergent relationships, exceptions being the constraints on ECS using global mean temperature change in Jimenez-de-la Cuesta and Mauritsen (2019) and Nijsse et al. (2020). This is also not necessary if theoretical reasoning suggests a different relation is more suitable. Equivalently to the use of multiple predictors/predictands, one could imagine using vectors for X and Y , rather than scalars if theory suggested such a relationship. These vectors could for example encode spatial fields of a climate variable (Brown and Caldeira, 2017), spatial normal modes of variance, referred to a empirical orthogonal functions (EOFs) in meteorology (von Storch and Zwiers, 1999) or the leading linear dynamical modes, also known as principal oscillation patterns (POPs) in meteorology (Hasselmann, 1988;Williamson and Lenton, 2015). One prediction we make is that the discovery of new ECs will be led more by those asking questions as to how climate impacts of concern may adjust in the future. Many altered meteorological features expected to have the most detrimental effect on societal safety, livelihoods and well-being are the magnitude and frequency of extreme weather events. A focus on extremes may result in ECs that are different in composition to those discovered to date. Rather than a fluctuating and observable contemporary quantity (X) projecting a 'bulk' climate property (Y ) of importance to the future, instead a different format is required. That is, summary climate properties of large-scale features (or trends in them), as observable quantity X, require investigation to determine if they can estimate the future extreme frequency of quantity Y . In addition, a search for such ECs will focus attention on three further challenges. First, attributes of extremes and their changes might show substantial geographical variation. Second, by definition, extremes are rare, and so there are few data points available for investigation in ESM diagnostics where meteorological values are above thresholds that may be unsafe. Third, the future connection between observable X and extreme feature Y may involve a teleconnection between locations, and possibly with seasonal inertia. Presently, the most commonly used ESM ensembles in climate as well as EC research are multi-model ensembles (MMEs, see section II.A) allowing structural uncertainty to be explored and in more recent MMEs such as the CMIP datasets, initial condition uncertainty. Parametric uncertainty is also captured to some degree although not in a systematic way unlike in PPEs. It is therefore unclear whether all plausible model parameteric configurations are explored. See also sections IV.I and VI.D. The strengths of MMEs and PPEs could be combined in a 'super-ensemble' of PPEs generated from structurally different models. Apart from better assessment of confidence intervals, questions about the necessity of plausible ESM simulations for good ECs could be answered i.e. is it necessary to filter out unrealistic simulations to strengthen the EC, based on an observation unrelated to the predictor? Statistical emulators (Sacks et al., 1989) of ESMs may help to give enough effective simulations after filtering unrealistic simulations from a PPE. Use of machine learning (ML) techniques in ECs may prove fruitful in future research. Machine Learning (ML) (Webb and Copsey, 2011) is a wide definition but encompasses ECs, particularly supervised learning type al-gorithms. In a supervised ML approach, the algorithm builds a mathematical model from a set of training data that contains both inputs and desired outputs. In the case of ECs, this training data is from an ESM ensemble, the inputs being the observable X and the outputs being the future projection Y . The mathematical model is f i.e. Y = f (X) + . The mathematical model that has been a popular choice in ECs so far has been f (X) = aX + b where a and b are constants fitted via linear regression. Phrasing linear regression as an example of supervised ML could be seen to be an over complicated generalization. However, viewing ECs within this framework potentially allows one to take advantage of the existing tools within ML.
For example, the mathematical models fitted by the supervised ML algorithms could be more elaborate with more physical or statistically based parameters that give better explanatory power or more realism. The inputs and outputs equally could be more elaborate having extra spatial dimensions and/or extra variables (section VII.D). This could of course get as complicated as reproducing the original models that the inputs and outputs derive from. However, one ideally would like to walk the line between physical plausibility/understandability of the model and explanatory power to derive a minimal description. One also has to be aware of the number of data available to fit to a particular mathematical model. Presently, MMEs consist of around 30 ESMs limiting defensible fits of mathematical models to just a few parameters. However, one may be able to leverage the extra information in the ESMs spatial and temporal degrees of freedom to increase data size and therefore mathematical model complexity. Though our expertise is in the area of climate and Earth system dynamics, we have found a few examples of EC-like approaches in other fields. In the closely related field of weather prediction, an approach known as 'Model Output Statistics' is used to enhance the quality and relevance of weather forecasts produced by numerical models (e.g. Veenhuis (2013)). As models are used again and again to make forecasts, an archive of past simulated data accumulates. This archive can be analysed to diagnose systematic statistical relationships between model output and observed quantities of interest. (An example of a quantity of interest might be temperature at a particular weather station.) When the models are then used to make real-time forecasts, these relationships can be applied to the forecast output, generating values for the quantities of interest. This is similar to the EC approach in that the structure of model biases is ascertained and exploited to produce future predictions that are likely more realistic. We find another example in the field of cosmology, where models have been created to simulate galaxy development. Like climate models, these models have tunable parameters. When the model parameters are varied to produce an ensemble of simulations, statistical relationships among galaxy components across the ensemble can be analysed. In the study of Terrazas et al. (2020), their galaxy model produces strong statistical relationship between the black hole mass and the associated stellar mass (c.f. Figure 7 of that paper). Under certain conditions, observations of black hole and stellar mass mirror the simulated relationship. This suggests that given an observation of black hole mass for a galaxy without a corresponding stellar mass observation, the model could be used to make a meaningful prediction of stellar mass. This example is similar to ECs in the sense that the intrinsic model (and real world) physics connecting key quantities of interest can only be ascertained through analysis of a multi-model ensemble. In the field of economics, past forecasts of economic performance using multiple economic models have been analysed and compared to the actual ensuing economic performance. Felix et al. (2018) found that the skewness of the ensemble forecast's distribution is a strong predictor of economic surprises. The implication is that if a current ensemble of economic models produces a skewed forecast distribution, the likelihood of an economic surprise is elevated. Such an outcome might not be detectable from more conventional predictions, such as the ensemble mean forecast. This example has some similarities to the Model Output Statistics example from weather forecasting: Model biases are ascertained from multiple realisations of past forecast performance, and used to improve the current forecast. This brief survey of a few close cousins to the EC approach reveals that ECs and their variants are most likely to be useful where we are trying to understand complex systems where experimentation with the real system is not possible or is cumbersome. In such systems we often need to simulate dynamical evolution, initialising from a snapshot in time or at best a relatively short period of observations. Multiple models of the system are possible, and which model is the correct one is unknown. Emergent relationships across model variants are a way to reveal the deeper connections between observable elements of simulations and model outcomes of interest. There are likely many more examples of EC-like research being done on complex systems. These disparate communities appear to be unaware of one another, and could benefit from organised efforts to convene in workshops and meetings. The spread in climate change projections has not reduced substantially despite major advances in model resolution and process understanding and inclusion over recent decades. The lack of progress represents a disappointment for climate science, and hinders society's ability to plan for future impacts. These challenges cannot be overcome by just increasing the resolution of ESMs and including more and more physical and biogeochemical processes. A balanced interplay of observations as well as top down, conceptual and bottom up, comprehensive modelling and theory is required to achieve the required progress and scientific understanding. In this review, we have argued that the EC approach offers a promising way to incorporate all of these elements and ultimately reduce key uncertainties in future climate. The development of ECs has arisen from the requirement to reduce often substantial inter-ESM differences in projections of climate. Large model differences makes adaption planning difficult and risks spending funds set aside to help societies cope with climate change. For those tasked with formulating mitigation plans that will reduce greenhouse gas emissions in order to contain global warming, uncertainty prevents the formation of accurate trajectories of 'allowable' fossil fuel burning. In this review paper, we have shown existing applications of ECs, discussed theoretical issues surrounding them that likely require further investigation, and characterised a few circumstances where they may start to fail. Nevertheless, although this seemingly points to more research required, ECs are a key methodology to distil the often conflicting information from climate modelling centres across the world.
For this reason, we anticipate that the application of ECs will continue to grow, and results from their application will have an increasing role in projecting future change, and with strong representation in future UN reports by the IPCC. This review has looked in depth at how ECs have become a standard methodology, used by climate researchers to synthesize substantial inter-ESM prediction differences, into projections of change with lower levels of uncertainty. We make a critical assessment of the EC technique, and many of the open questions raised here may lead to exciting avenues of research over the years ahead. It may also be true that other research disciplines will find a role for the EC approach. However, it is always worth recalling that the main scientific issue remains pressing, as without substantial reductions in emissions, the climate is expected to change substantially and to a dangerous state for large fractions of society. For this reason, intense scrutiny of the EC method is highly welcome and appropriate. The relentless need for climate predictions is also likely to trigger the discovery of new ECs, caused instead by impacts led requirement for information taking precedence over more curiosity driven investigation. This paper benefited from discussions at the (AGCI) This paper benefited from discussions at the (AGCI) In a discussion paper by Renoult et al. (2020), a simple method was proposed to combine independent ECs to create a tighter estimate for Y . Where the regression is normally given as p(y|x), they propose to instead formulate the statistical model as p(x|y), allowing for a prior on (y) to be integrated into the emergent constraint as p(y) here is the posterior distribution of a previous EC. To make sure the two ECs are indeed independent, the authors state that observations need to be independent, and that insofar possible, the errors in models should also stem from different sources. Their example involved a warm and cold climate state for which temperature change was reconstructed. Temperature change is dominated by different processes in this case, so that model error can be considered independent to first order. This method is not consistent with other methods described above. Linear regression is typically not symmetric; regression where X predicts Y , p(y|x), describes a different function than regression where Y predicts X, p(x|y) (Smith, 2009), as illustrated in In this section we give some directions which we think are promising or exciting for future EC research. The collection of ECs described in Section V highlight key areas of focus to date, but this also demonstrates In a discussion paper by Renoult et al. (2020), a simple method was proposed to combine independent ECs to create a tighter estimate for Y . Where the regression is normally given as p(y|x), they propose to instead formulate the statistical model as p(x|y), allowing for a prior on (y) to be integrated into the emergent constraint as p(y) here is the posterior distribution of a previous EC. To make sure the two ECs are indeed independent, the authors state that observations need to be independent, and that insofar possible, the errors in models should also stem from different sources. Their example involved a warm and cold climate state for which temperature change was reconstructed. Temperature change is dominated by different processes in this case, so that model error can be considered independent to first order. This method is not consistent with other methods described above. Linear regression is typically not symmetric; regression where X predicts Y , p(y|x), describes a different function than regression where Y predicts X, p(x|y) (Smith, 2009), as illustrated in In this section we give some directions which we think are promising or exciting for future EC research. The collection of ECs described in Section V highlight key areas of focus to date, but this also demonstrates