Using climate models in climate assessments: a Bayesian perspective

1 Introduction

In the absence of alternative habitable planets, complex computer climate or Earth system models (ESMs) provide the only means of doing planet-wide experiments [e.g., (Jones et al., 2016; Gillett et al., 2016)]. They are also useful in isolating processes (Reed et al., 2025), projecting different futures (O'Neill et al., 2016), and understanding Earth's past climate states (Kageyama et al., 2021). There are now 49 separate modeling groups worldwide contributing nearly 100 models to the Coupled Model Intercomparison Project (Eyring et al., 2016), whose current (CMIP6) and next (CMIP7) generations will inform the next IPCC report. These ESMs differ in spatial resolution, in the number and type of processes they explicitly simulate, and in the parameterizations, they used to characterize sub-grid-scale processes. However, these ESMs are not independent of each other (Knutti, 2010; Knutti et al., 2010; Sanderson et al., 2015; Masson and Knutti, 2011). Developers read a common literature, discuss freely with each other, and often even share or inherit code (Kuma et al., 2023). Major advances in modeling often occur simultaneously across modeling centers as computing improves [e.g., the move to higher resolution (Gosling and Arnell, 2013) or the inclusion of complex processes such as aerosol – cloud interactions (Cherian and Quaas, 2020; Wang et al., 2021)]. Nor is it immediately obvious how the CMIP ensemble relates to the real-world values of quantities we wish to estimate: CMIP models are models, not samples of the real world, and as such will exclude processes that matter in reality (Lloyd and Winsberg, 2018). As a result, it is widely acknowledged that CMIP results are not an ensemble generated by an underlying random process, but an “ensemble of opportunity” (Tebaldi and Knutti, 2007; Knutti et al., 2010; Touzé-Peiffer et al., 2020) that provides a definitely incomplete, likely biased estimate of the reasonable spread of possible values.

How, then, should scientists analyze and report the results of different CMIP models? Here, we argue that a Bayesian framework provides the optimal lens for viewing the multi-model archive because it allows for the use of prior knowledge, is easily updated with new evidence and results, and – perhaps most importantly — forces us to make explicit the assumed connection between the models and the real world. We present several ways to translate CMIP output into an inference about the real world and show how these lead to different probability distributions for climate variables. We then show how these may be updated with new observations or modeling experiments, resulting in a transparent and iterable framework for climate assessment.

This article is intended to be both a guide for scientists wishing to use Bayesian methods to interpret CMIP evidence and an introduction to a set of analysis models available at https://github.com/netzeroasap/CMIP_BAYES. We hope it will be of use to scientists conducting assessments of quantities such as equilibrium climate sensitivity (Sherwood et al., 2020; Marvel and Webb, 2025), transient climate response (Tokarska et al., 2020), or transient climate response to emissions (MacDougall, 2015). To that end, we adopt an approach in which we present a series of interpretive models with gradually increasing complexity (Harte, 1985).

2 Description vs. inference

The CMIP6 ensemble mean equilibrium climate sensitivity (ECS) is 3.7 °C, and the sample standard deviation is 1.0 °C (Zelinka et al., 2020). The ensemble mean value of β_*L*, the terrestrial CO₂-fertilization effect (Gregory et al., 2009; Friedlingstein, 2015), is 0.85 PgC ppm⁻¹ with a sample standard deviation of 0.36 PgC ppm⁻¹ (Arora et al., 2020). There are 13 million dogs in the United Kingdom (McMillan et al., 2024). All three of these statements are true, but in isolation none of them tells us anything about the real-world value of ECS, β_*L*, or any other climate-related quantity. These are useful descriptive statements about CMIP6 (or dogs). But what we generally want to do is not description but inference: we want to know what these pieces of evidence teach us about the real world (McElreath, 2020).

Intuitively, the British dog population probably provides little information about the real-world equilibrium climate sensitivity. Equally intuitively, we might believe that the physics-based Earth system models participating in CMIP do say something about ECS. Our assumptions about the relationship between CMIP and the real world are what allows us to make inferential statements about Earth's climate.

Unfortunately, descriptive ways of presenting CMIP evidence can be misleading. We often present statistics such as the CMIP multimodel mean and standard deviation in a way that implies a particular connection between the ESMs and the real world. In practice, however, this choice corresponds to a particular interpretive model whose assumptions are rarely stated explicitly or subject to scrutiny (hereafter, we use “model” to mean such an interpretive model and “ESM” for any climate model participating in CMIP, whether or not it contains an interactive carbon cycle; this includes all complex general circulation models). If we aim to systematically include observational data and update our understanding as data and ESMs evolve, we must be explicit about how we are linking CMIP to the real world. Inference requires an analyst to make choices. The Bayesian framework makes these clear and subject to scientific criticism.

3 Inference with CMIP

CMIP simulations don't tell us everything: they are not perfect representations of reality. But they don't tell us nothing. Few scientists would dispute that these constantly refined and physics-based ESMs provide at least some information about reality. By using these ESMs, we implicitly treat their results as evidence of the real world.

The focus of this article is on the assessment of some scalar quantity X using CMIP. X may be a quantity such as the global mean temperature change that is readily observable in the real world. Or it may be a parameter (such as ECS) that is accessible only in ESMs and may be only indirectly constrained by other lines of evidence.

Bayes' theorem states

P(X|CMIP)=P(CMIP|X)P(X)P(CMIP)(1)

The term on the left-hand side of Equation 1 is generally what we want: the posterior, or our new estimate of the climate quantity X updated by the evidence of the CMIP ensemble. The term P(X) is the prior: our knowledge of X before using CMIP evidence, perhaps based on fundamental physics or intuition about the Earth system. The denominator P(CMIP) is a normalization constant ensuring that the sum of all resulting probabilities is 1. The likelihood P(CMIP|X) is where the meaningful analysis happens: we cannot update our understanding of X without it. Choosing a likelihood function means asking “for a given value of the “true” X, how likely would it result in the CMIP ensemble we have?” Answering this question requires us to make choices about CMIP and the real world. These choices, unavoidable in any analysis, reflect our knowledge and beliefs about the relationship between climate model results and the real world. The posterior in Equation 1 is not only conditioned on the CMIP6 ensemble, but also on the underlying model we use to link CMIP results to the quantity X we want to infer. Instead of merely describing the CMIP distribution, the resulting inference answers the question: what possible X are compatible with CMIP?

We show below how different interpretive models inform assessments of a climate variable X. If this quantity is not directly observable, Earth system model evidence becomes informative about the real world only through explicit interpretive models, which encode different assumptions and yield different inferences. The uncertainty we report is largely a consequence of how we decide to view the evidence.

3.1 A trivial model: a single ESM is “true”

We begin with a trivial example. Suppose we believe a single ESM—for example, UKESM1-0-LL (Sellar et al., 2020)—to be a perfect digital twin of the real world and all other ESMs to be irrelevant. In this case, our likelihood function would assign a value of 1 to X_{*UKESM*1 − 0−*LL*} and a value of zero to all other X_*m*. This interpretive model, manifested in the design of the likelihood function, yields a posterior P(X|CMIP) that is a delta function centered on the value of X in UKESM1-0-LL. If we start our analysis believing the ESM is perfect, no updating with the evidence of other CMIP ESMs is possible.

Using this interpretive model, we would be 100% confident that the equilibrium climate sensitivity is 5.36 K, that the terrestrial CO₂ fertilization effect β_*L* is 1.03 PgC ppm⁻¹, and so on. But this vastly increased (and misplaced) confidence does not mean we have meaningfully reduced uncertainties in these quantities. The narrow posteriors are simply a consequence of the way we interpret CMIP; the uncertainty estimate arises from the interpretive model we use to construct the likelihood. This particular model is highly unrealistic and not explicitly in use anywhere in the literature (although one could argue that climate studies or experiments using a single model implicitly assume it). Still, this simple example highlights something important: all posteriors arise from a specific way of interpreting CMIP.

3.2 The fixed effects model: ESMs differ by sampling error

Now, suppose we treat the value of X in the jth ensemble member of each ESM m as

Xmj=X+ϵmj(2)

Where ϵmj quantifies the “sampling error” that causes a simulated value to deviate from the real world. This is a Bayesian “fixed-effects” model, which, while far too simple for most uses, is nevertheless an important building block in the statistical model hierarchy we discuss here.

To construct the likelihood, we require an explicit model for the noise term ϵmj. In a situation where the sampling error results from numerical or implementation noise, we might model ϵmj for all ESMs and ensemble members as a random variable drawn from the same normal distribution denoted (N) with underlying standard deviation σ:

ϵmj~N(0,σ).(3)

If noise results from internal variability, a non-normal distribution with fatter tails, skew, or autoregressive behavior may be more appropriate (Marvel et al., 2025; Sura and Hannachi, 2015). The salient point is that the distribution we assign to ϵmj (as in Equation 3) encodes our assumptions about the source of the noise. It specifies why we think the values of X's simulated by each ensemble member may differ from one another (and from the real world).

Equation 2 along with the noise model Equation 3 specifies the underlying generative model for the CMIP ensemble. This means the likelihood of observing Xmj given the true underlying value of X is simply the familiar Gaussian probability density

P(Xmj|X)=1σ2πe-(Xmj-X)22σ2(4)

and the likelihood of the entire CMIP ensemble given X is

P(CMIP|X)=(12πσ2)Ntotexp(-12σ2∑m=1nESM∑j=1jm(Xmj-X)2)(5)

where n_*ESM* is the number of ESMs in the CMIP ensemble, j_*m* is the number of ensemble members for ESM m, and N_*tot* is the total number of simulations. Using the fixed-effect model, Equation 5 is the likelihood term used in the numerator of the right-hand side of Equation 1.

The fixed-effects model is often productively applied in analyses in which the variation between ensemble members is attributed to (for example) numerical noise, parameter choices, or internal variability-often in the context of a single ESM or multiple ESMs assumed to center on the same value and to have similar noise characteristics. Here, we show an example of how it may be applied to the multimodel ESM–not to recommend its use, but as an illustrative example. The fixed-effects model encodes the assumption that all ensemble members of all ESMs are otherwise identical in their representation of the world, statistically independent, and would therefore be expected to deliver the same result in the absence of sampling noise. The more ensemble members of a given ESM m, the more confident we can be in the “true” latent (i.e., a quantity that appears in the model but is never directly observed) value of X in that ESM. Given this model, the more confident we are in a single ESM's value of X, the more confident we can be in the true X. But this estimate is conditional on our assumption that all ensemble members differ from the real-world value of X only by a common sampling error. This is clearly an unrealistic assumption, but it is instructive to see where it goes wrong.

As an example, consider applying a fixed-effects model in which X is the forced component of the annual mean temperature anomaly (relative to the 1951–1980 base period) in the year 2000, denoted T. We want to calculate the posterior P(T|CMIP) of T (how did external forcing contribute to the temperature anomaly that year?) conditioned on the CMIP output. Here, we use results from biogeochemically coupled CMIP ESMs driven by emissions (Arora et al., 2020). If we assume a normal distribution for the noise, which represents both numerical implementation noise and internal variability in the annual mean temperature, the model has two parameters: T, the “true” real-world value, and σ, the width of the underlying noise distribution. This means we must specify priors on each of these variables.

3.2.1 Priors

The width of the underlying noise distribution cannot take negative values. The prior distribution we select reflects this: following standard advice for priors on variance (Gelman, 2006), we use a half-normal prior for σ:

σ~N+(10).(6)

Equation 6 is simply a normal distribution truncated so that the probability density for negative values is zero. We use a width of 10- an extremely conservative choice that allows a very wide variation in temperature due to internal variability.

We must also specify a prior on T. Here, it encodes our previous knowledge of the forced component of the year-2000 annual temperature anomaly. We test the sensitivity of this analysis to two choices: a wide prior in which the forced component of the temperature anomaly is considered almost completely unknown:

Wide prior: T ~ N(0,100)(7)

and a narrow prior reflecting a confident view that the global temperature has not been forced

Narrow prior: T ~ N(0,.01).(8)

Regardless of the prior on T, the fixed-effects model results in similar, extremely narrow posteriors (Figure 1, blue shaded distribution and blue dashed distribution). This posterior reflects the underlying assumption that each ESM simulation constitutes an independent sample. This means the more simulations we have as evidence, the more strongly we update our prior. On the other hand, when we have little evidence, we have to lean more on our prior knowledge.

3.3 The random-effects model: ESMs differ in structure

A more reasonable assumption is that different CMIP ESMs would be expected to differ due to structural factors: parameterizations, the inclusion or exclusion of various processes, resolution, etc. In that case, we would model a quantity X in the j^*th* ensemble member of ESM m as

Xmj=Xm+ϵmj(9)

where ϵmj is the noise. As with the fixed effects model, the more ensemble members we have, the more confident we can be in X_*m*. However, in this case X_*m* is the (latent) value of X in that ESM, not the true real-world X.

To link the ESM-simulated values X_*m* to the real world, we might allow for each ESM to have a specific bias b_*m*:

Xm=X+bm.(10)

If we have no reason to believe any given ESM to be more or less biased than the others, we might assume these to be distributed around zero:

bm~N(0,σstruct).(11)

This is a canonical random-effects (RE) model, commonly used in Bayesian meta-analysis (Gronau et al., 2021). The RE model contains a parameter σ_*m* for each ESM and also adds a parameter σ_*struct* that reflects ESM structural differences. This RE framework resembles the analysis choices that many climate scientists already make intuitively: when looking at a multimodel ensemble (where some ESMs contributed many more ensemble members than others), we know to take the mean over ensemble members first, then take the mean across models. The RE intuition is already used in the climate field, albeit seldom within an explicit Bayesian framework.

To illustrate this model, we return to the previous example in which X = T, the forced component of the year-2000 temperature anomaly. The orange shaded distribution in Figure 1a shows the resulting posterior P(T|CMIP) assuming half-normal priors on the structural spread σm~N+(1)∘C, half-normal priors on each ESM's internal variability σm~N+(1)∘C, and the wide prior (Equation 7) on T. These posteriors are far narrower than the spread in CMIP6 (magenta ticks along the x-axis in Figure 1a). While it may seem counterintuitive that the process of linking CMIP to the real world results in a distribution narrower than the standard deviation across CMIP, we stress that this posterior is a consequence of the assumptions about CMIP ESMs encoded in the random-effects model. This interpretive model assumes that CMIP ESMs have structural errors and that no ESM is more or less likely to be biased than any other. The resulting P(T|CMIP) is thus much narrower than the full CMIP distribution, which includes many random manifestations of internal variability.

Given the same prior on T, the posterior conditioned on an RE model is wider than that for the fixed-effects model, reflecting that in addition to internal variability, CMIP spread also results from structural differences between models. It is also shifted to the left relative to the fixed effects model, because the assumption that each ensemble member informs T_*m* rather than T directly no longer gives undue weight to higher-sensitivity ESMs with larger ensembles. The posterior reflects uncertainty in each ESM's implementation noise/internal variability as well as in the structural spread. As expected, the ESMs with larger ensembles have more constrained posteriors on σ_*m*, while uncertainty in those with fewer ensemble members is dominated by the prior on σ_*m* (Figure 1b).

The random-effects model is sensitive to prior knowledge about T and σ_*struct* by design. If the prior on the climate quantity X is tightly constrained but the CMIP structural spread σ_*struct* is unknown and weakly identified, the model can favor the values of T with large prior weight by inflating the structural spread (dashed orange line, Figure 1a). In plain language, if we're confident in T and uncertain about the structural spread in CMIP, then the ESM values far from the T prior will be explained away as structurally biased. In the case of the annual mean temperature, if we begin with a confident belief that there has been no warming, the fact that the CMIP models all say otherwise will be attributed to a problem with the ESMs, and we will not update our beliefs at all. If we're already sure what T is before looking at CMIP, CMIP data will not change our minds.

The random-effects model will give overconfident posteriors if CMIP ESMs consistently over- or underestimate a quantity due to incomplete representation, parameterization, or exclusion of physical processes in the ESMs. If a subset of models is biased, then the truth does not necessarily lie at the center of the model distribution as assumed by the RE model. But the CMIP6 archive, in isolation, does not tell us its relationship to the real world. If we think the ESMs may be biased, we need to use a model that lets us take this into account.

3.4 The correlated random-effects model: ESMs may all be biased in similar ways

The random-effects model above assumes that biases exist but are uncorrelated across the ESMs, and therefore CMIP results are distributed about the real-world value of the quantity of interest. This is equivalent to the paradigm of a statistically indistinguishable ensemble discussed in Annan and Hargreaves (2010), in which the true value of X is assumed to be drawn from the same distribution as the CMIP-simulated values. However, CMIP is not a collection of independently developed ESMs. We know that CMIP ESMs share code, schemes, parameterizations, and, in some cases, developers, consult a shared literature, and use the same observational datasets for tuning. CMIP might therefore be expected to have a common bias due to errors or parameterization choices shared across ESMs. To take this into account, we can write (for any quantity X):

Xmj=Xm+ϵm(12)

Xm=X+bshared+bm(13)

Equation 12 encodes the assumption that is: the j^*th* ensemble member of ESM m varies about the central value X_*m* in that ESM due to internal variability described by ϵ_*m*. In addition, Equation 13 encodes the assumption that the central value X_*m* in each ESM differs from the true value X by a bias b_*m* (specific to that ESM) and a common bias b_*shared* (shared by the entire CMIP ensemble).

Rather than use different unknown distributions for b_*m* and b_*shared*, we let ρ represent the fraction of the total bias shared across all ESMs. In this case, we model

bshared~N(0,ρσstruct2)(14)

bm~N(0,(1-ρ)σstruct2).(15)

In the correlated random-effects (CRE) model described by Equations 14, 15, the parameter ρ measures the correlation between the biases in any two models. If ρ = 0, there is no common bias: all of the ESMs are distributed about the “true” value of X, the CRE model reduces to the RE model, and unless we have strong prior reason to believe otherwise, the CMIP mean is our best estimate of X and confidence in that estimate will grow with the CMIP ensemble size. If, on the other hand, ρ = 1, all of the bias is shared: the CMIP ESMs effectively move collectively as one, and we have no real way of knowing how far away from the true X they might be. A nonzero ρ breaks the assumption of a statistically indistinguishable ensemble, because when we allow for a common bias among the ESMs we are admitting the possibility that the truth may be drawn from a different distribution than that of CMIP.

One limitation of the CRE model is that it assumes no ESM is likelier to be more biased than any other. It also ignores known ESM genealogy, assuming that no two ESMs are more related than any others. In this formulation, we don't know which models are related, only that they likely are. The advantage of the model is that it contains few parameters: rather than inferring an entire covariance matrix for CMIP ESMs, we only need to consider the overall bias and the fraction of it that is shared.

As an example, consider ECS estimates from CMIP6 as calculated in Zelinka et al. (2020).

If we assume a standard random-effects model with ρ = 0 and small internal variability σ_*m* = 0.1K across all models, the posterior P(ECS|CMIP) obtained with the uncorrelated random-effects model is extremely tightly constrained (blue line, Figure 2). But if we instead assume that exactly half of the bias is common across models (for example, due to a real-world process missing or poorly resolved in CMIP6), the resulting posterior is much wider (orange line, Figure 2). If the ESMs have a common bias, it becomes much more difficult to dismiss a high- or low-sensitivity ESM as an outlier — the correlation forces the CMIP6 ESMs to move together, and thus inflates the posterior. This is a crucial point: if all the CMIP ESMs might be biased in their estimation of ECS, we will be less certain of its true value as a result.

In general, we do not know ρ (although we can generally be confident based on the model genealogy literature that it is nonzero). To model this uncertainty, we might place an informative prior on ρ, allowing for the possibility that the common bias may be large or small. The green line in Figure 2 shows the result when we assume a beta distribution, chosen because it is defined on the interval [0,1] and has mean 0.5 [ρ~Beta(5, 5)]. The posterior for ρ in this case is extremely similar to the prior (Figure 2b). This means that our previous knowledge or belief about the structural biases in CMIP heavily informs our resulting conclusions. This is not particularly surprising: the CMIP archive cannot, on its own, tell us about the relationships between the participating ESMs. We need additional information about the origin or nature of the shared biases.

3.5 Accounting for process biases

What if we explicitly know the source of some common or specific biases in the ESMs? In the case of ECS, this might be related to cloud parameterizations or the absence of a particular feedback process. For other quantities, however, we know exactly which process is missing and have information about its likely effect.

As an example, we will consider another quantity simulated in CMIP Earth system models. Following Friedlingstein et al. (2003), we assume that the simulated changes in land carbon storage ΔC_*L* can be approximated as a linear combination of effects resulting from the change in atmospheric CO₂ concentrations ΔC_*A* and effects resulting from the change in global temperature ΔT:

ΔCL=βLΔCA+γLΔT.(16)

Here, we will use CMIP data to constrain the coefficient β_*L* that quantifies the land carbon fertilization effect: how terrestrial carbon stores change in response to increasing atmospheric CO₂. For increased clarity in this article, we will drop the L subscript that denotes land in Equation 16 and refer to this quantity as β.

We know of at least one clear source of bias in some of the Earth system models participating in CMIP. Both process-based modeling studies and experiments suggest that the ability of terrestrial ecosystems to respond to increases in atmospheric CO₂ may be constrained by nitrogen availability (Zaehle et al., 2015; Meyerholt et al., 2016; Davies-Barnard et al., 2020). However, roughly half of the CMIP6 ESMs do not explicitly simulate the consequences of terrestrial nitrogen (N) limitation for land carbon uptake (Gier et al., 2024). This missing process might therefore be expected to bias the β_*m* simulated by each of these ESMs, and thus CMIP6-derived estimates of β. What interpretive model might capture this expectation?

First, we must explicitly model how we expect the inclusion of this specific process to affect each ESM's β_*m*. It seems reasonable to assume that nitrogen limitation can scale the CO₂ fertilization effect up or down. Thus, in the real world, the CO₂ fertilization effect is the product of its unscaled value U and a scaling factor η_*N* that reflects the effects of nitrogen limitation:

β=ηNU(17)

The ESMs that include the process give us information about the scaled quantity β in Equation 17 while those that exclude it give us information only about U. The differences between these groups of ESMs provide information about η_*N*.

Thus, we model each ESM's calculated CO₂ fertilization effect β_*m* as

βm=ηmLmUm(18)

Um=U+ϵu(19)

ηm=η+ϵη(20)

where

Lm={1,if model m simulates N limitations0otherwise

In this model, an ESM's β_*m* gives us information about the unscaled effect and, if the process is included, the scaling factor as well (Equation 18). Each model's unscaled effect is assumed to be the real-world U plus a bias term (Equation 19), and similarly for the scaling factor (Equation 20). Crucially, the posterior estimates of β are now conditioned on two sources of evidence: CMIP6 ESMs, and a binary vector L that encodes information about which models contain the process. In the absence of the information in L and an interpretive model for the effect of including this process, a simple average over ESMs encodes the implicit assumption that nitrogen limitation has no effect (η = 1).

As a result, the posterior conditioned on both the CMIP6 archive and our knowledge about which models have nitrogen limitation, denoted P(β|CMIP, L), is shifted lower compared to the posterior P(β|CMIP (Figure 3a). This is because the posterior for the scaling factor (Figure 3b) indicates that η is very likely to be less than unity.

This hierarchical framework informs our estimates of η through the difference between ESMs with and without nitrogen limitation. It propagates the uncertainty in η (how does the inclusion of N limitations scale the CO₂ fertilization effect?) and the uncertainty in U (what is the intrinsic structural spread in the unscaled effect U) through to the final posterior. The major advantage of this framework is that it allows us to incorporate new evidence- not just from standard ESM runs, but from targeted process-based runs and potentially relevant observations, which we illustrate in Section 4.3.

3.6 Summary

There are many kinds of uncertainty. The parameters defined by the interpretive models we have introduced here relate to three fundamentally different sources: simulation spread due to noise or internal variability from running a model, inter-model spread due to differences among models, and the bias of the whole CMIP ensemble compared to reality. For any given ESM, the outcome of a single simulation is a combination of the forced signal and internal variability. Simply running the same ESM many times allows a much more confident diagnosis of its forced signal, as is common in the detection and attribution literature (Hegerl et al., 1996; Allen and Stott, 2003; Santer et al., 2011). The fixed-effects model would be appropriate in the single-ESM case, and as the sample size grows, our estimate of X becomes more confident. But this does not reduce the uncertainty or error of the ESM itself. Similarly, if we view CMIP as an ensemble that spans many ways of attempting to represent the real world, then the more modeling groups participate the more confidently we can identify the best estimate that CMIP can give us. Finally, these remaining errors or biases in CMIP are represented in the correlated random-effects model and can only be addressed if we incorporate specific information about the CMIP ensemble such as knowledge of processes known to be missing or treated incompletely in the ESMs.

Every level of complexity introduces new parameters into the model, which may or may not reflect important aspects of the real world. The goals of model selection are parsimony and accuracy: we want to interpret data in the simplest way possible, but no simpler.

Other models to interpret CMIP data are possible: nontrivial correlation matrices that reflect known model genealogies, weighting of “good” ESMs relative to “bad” ones, etc. Leo Tolstoy famously wrote that “happy families are all alike; every unhappy family is unhappy in its own way”. Similarly, all perfect Earth system models are alike, but all biased ESMs are biased in their own way. How are we to decide the right criteria for model weighting? Ultimately, this is a question of how (and whether) we can update our CMIP-based understanding with other lines of evidence.

4 Updating with new evidence

Here, we show how to incorporate new sources of evidence into our Bayesian framework. To parallel the structure of the previous sections, we first show how to update T with observations, ECS with emergent constraints, and β with process-based modeling.

4.1 Direct observations

If we have sufficiently confident observational evidence, the way we choose to interpret CMIP matters much less. Unless we use an extremely restrictive model (e.g., assuming a single ESM result must be the truth), the weight of the observational evidence will shift most reasonable posteriors regardless of how we interpret the ESM results.

Consider, for example, the realized annual mean temperature anomaly in 2000: how warm, exactly, was that year? This is a well-observed quantity. For illustrative purposes, we will update CMIP-derived posteriors using two observational ensembles: HadCRUT (Morice et al., 2021) and GISTEMP (Lenssen et al., 2019). We assume the uncertainty ensembles provided by both observational groups represent draws from appropriate posteriors and sequentially update our estimates using a hierarchical model:

T2000=T+ϵ(21)

THadCRUT=T2000+ϵHadCRUT(22)

TGISTEMP=T2000+ϵGISTEMP.(23)

Here, T in Equation 21 is the forced component whose posterior we obtained in Section 3.2. ϵ is the amplitude of internal variability in CMIP, which here we will take to be normally distributed with width 0.1K. ϵ_*HadCRUT* in Equation 22 and ϵ_*GISTEMP* in Equation 23 are the uncertainty terms in the published observational ensembles.

Because the HadCRUT and GISTEMP observations are highly constrained and in good agreement with each other, the updated posterior is nearly identical no matter which interpretive model or which priors used to interpret the CMIP archive (Figure 4). Plainly, if we have plenty of observational evidence, it doesn't matter what (or how) we think of the CMIP ESMs.

4.2 Emergent constraints

For a quantity that is unobservable and accessible only in ESMs, such tight observational constraints are impossible. However, it is possible to narrow the posteriors indirectly via so-called “emergent constraints” (ECs), in which the relationship (in CMIP) between an unobservable X and an observable quantity is assumed to hold in the real world, and thus to provide information about the true value of X (Hall and Qu, 2006; Williamson et al., 2021). Traditionally, such constraints are derived by using ordinary least-squares (OLS) regression to fit a line to CMIP data, treating the observable in CMIP as the predictor and the CMIP-simulated X as the dependent variable (Caldwell et al., 2018; Schlund et al., 2020). The estimated slope, intercept, and observed quantity are then used to predict X.

However, we do not want to predict the unobservable X from observations: we want to determine what values of X are compatible with the observations. We also want to smoothly propagate uncertainty in the slope and intercept through to the posterior on X, allow for the use of any prior knowledge of X, and treat the CMIP ensemble as noisy values informing the regression, not as a fixed design matrix — all goals that are not easy or possible to accomplish in the standard frequentist OLS framework. The difference between prediction and inference will be most important when the regression is noisy (as in Figure 5a) and the observations lie at the edge of the CMIP distribution (as is common in emergent constraint problems).

Thus, following Renoult et al. (2020), we instead adopt a fully Bayesian hierarchical formulation for emergent constraints, in which the observable Y_*m* in CMIP and Y_*obs* in the real world are modeled as

Ym=aXm+b+ϵreg(24)

Yobs=aX+b+ϵobs.(25)

Here, CMIP simulations of the unobservable variable(s) X_*m* and observable variable(s) Y_*m* are used to learn the slope a and intercept b, along with the regression uncertainty ϵreg~N(σreg) (Equation 24). The observations and their reported uncertainty ϵobs~N(σobs) are then used to infer the “true” X (Equation 25). Note that in this framework the role of predictor and predicted is reversed: a CMIP ESM's value of X_*m*, or rather the collection of processes and parameterizations that lead to it, is treated as predictive of the observable Y_*m*. We are then inferring a posterior for the latent “true” predictor P(X|data) conditioned on the real-world values of the observable. This has the advantage of correctly propagating uncertainty in the slope, intercept, regression noise, and observations through to the resulting posterior for X.

For example, consider two emergent constraints on ECS derived from CMIP6 values provided in Table A2 of Schlund et al. (2020). The EC proposed by Brient and Schneider (2016), is referred to here, as in (Schlund et al., 2020), as BRI. It links ECS to the albedo of low clouds over tropical oceans, measured by the covariance of cloud fraction with SSTs on interannual timescales. Zhai et al. (2015), hereafter ZHA, as in (Schlund et al., 2020), proposed an EC based on the response of the marine boundary layer cloud fraction to SST changes in subtropical ocean regions. In CMIP6, the correlation between BRI and ECS is −0.34 (although the correlation was found to be much stronger in CMIP5) and between ZHA and ECS is −0.79 (scatterplots, Figure 5a).

Assuming a wide prior on the climate sensitivity ECS~N(0,100), priors N(0, 5) for the slope a and y-intercept b, and a half-normal prior with the scale parameter 10 on the regression uncertainty σ_*reg*, the BRI emergent constraint provides very little information to update our understanding of ECS. The relationship between BRI and CMIP ECS is simply very weak in CMIP6, and many lines with many slopes can be drawn through the data. If we have no prior information about reasonable values of ECS, then BRI provides almost no constraint on it. The resulting 95% highest posterior density is (–57.4 K, 88.0 K). Conversely, the relationship between ZHA and ECS in CMIP6 is stronger, leading to a much narrower posterior with a 95% HDI (3.3, 4.6 K) with the same priors.

Emergent constraints depend on the assumption that the relationship between the observable and the unobservable parameter is correctly simulated across the CMIP ensemble. A strong correlation between the observable and the quantity of interest in CMIP leads to overconfident results if shared model bias inflates the apparent precision of the fitted relationship; if the ESMs are not independent draws from a distribution, we have fewer effective data points for fitting the regression. Consequently, our structural interpretation of the CMIP ensemble remains a critical design choice when using ESMs as indirect evidence for the real world.

By explicitly accounting for ESM-specific and shared biases using a correlated random-effects model, we can avoid becoming overconfident in what CMIP can reveal, prevent overreliance on particular emergent constraints, and account for the possibility that emergent relationships may reflect shared ESM biases rather than real-world links. Consequently, we extend the methodology presented in Renoult et al. (2020) by combining a multiple emergent constraints model with a generative model for CMIP. The random-effects model protects against over-interpreting a single emergent constraint by pulling the posterior range back toward the region consistent with CMIP simulations. However, multiple strong emergent constraints all acting in the same direction will combine to shift the posterior away from the CMIP distribution, leading to the conclusion that the CMIP ESMs are likely to share a common bias, that structural uncertainty in the CMIP models is larger than previously thought, or both.

Figure 5c shows the posterior P(ECS|CMIP) assuming a correlated random-effects model with unknown ρ~Beta(5, 5) (green line, identical to the green line in Figure 2a). Updating this posterior with the new evidence from the BRI emergent constraint results in a narrower posterior (purple distribution, Figure 5c) than using BRI alone (purple distribution, Figure 5b). In this case, our correlated random-effects model for CMIP restricts the available parameter space of ECS, allowing even a weak negative relationship between ECS and BRI to provide meaningful constraint. We can sequentially update this distribution with multiple emergent constraints, seamlessly propagating uncertainty in the simulated ECS, simulated observables, and observations. The brown distribution in Figure 5d is the posterior for ECS when ZHA is used to sequentially update the BRI-constrained ECS posterior. Because the relationship between ZHA and ECS is much stronger, the resulting distribution is pulled toward values compatible with the ZHA constraint.

The posterior distribution shown in pink in Figure 5d results from sequentially updating:

P(ECS∣CMIP,BRI,ZHA)∝P(ZHA∣ECS,CMIP)×P(BRI∣ECS,CMIP)×P(ECS∣CMIP).(26)

Equation 26 applies if the emergent constraints BRI and ZHA are conditionally independent of one another. However, this is not the case: both of these observables depend on similar marine low-cloud processes and can be expected to provide redundant information. Indeed, across CMIP6 models, BRI and ZHA are moderately correlated (r = 0.4). To take account for this covariance, we calculate a posterior in which we do not assume conditional independence:

P(ECS∣CMIP,BRI,ZHA)∝P(BRI,ZHA∣ECS,CMIP)×P(ECS∣CMIP).(27)

To calculate the joint likelihood P(BRI, ZHA|ECS, CMIP) we stack the emergent constraints into a vector **X_*m*** = (BRI_*m*, ZHA_*m*) and account for the correlation by modeling

Ym=aTXm+b+ϵ(28)

with a = (a_*BRI*, a_*ZHA*), b = (b_*BRI*, b_*ZHA*). In Equation 28 the covariance ϵ~MN(Σ) is modeled by a multivariate normal distribution that captures the relationship between observables. Σ here is a two-dimensional covariance matrix rather than a scalar. The diagonal entries represent the variance of each constraint, with the covariance between the constraints on the off-diagonal. We can then assume that the true climate observables follow the same emergent relationship:

Yobs=aTX+b+ϵ.(29)

The inclusion of the same observational evidence using Equations 28, 29 widens the posterior (pink distribution, Figure 5d) because the interpretive model no longer assumes that each constraint provides independent information (Bretherton and Caldwell, 2020). In short, we no longer double-count the portion of the constraint provided by both processes.

Emergent constraints, while useful, are frequently criticized in the scientific community. Various emergent constraints do not persist across different generations of CMIP (Schlund et al., 2020), some lack a clear underlying physical mechanism, and publication bias means that apparently strong correlations receive more attention than other possible processes (Caldwell et al., 2018). The methods presented in this section update the posterior distribution with new evidence and account for relationships among multiple processes. However, the most robust way to incorporate emergent constraint information is to explicitly model the physical mechanism connecting the simulated observable process to the quantity of interest in CMIP and to specify assumptions governing the relationship between the observable and the real-world value of that quantity. In the next section we demonstrate an example of such explicit process-based modeling.

4.3 Process-based modeling

The Bayesian framework allows us to incorporate new model-based evidence – not just from standard ESM runs, but from targeted process-based experiments (Stein and Alpert, 1993) and potentially relevant observations. Suppose, for example, we have direct estimates η^*i* of the nitrogen limitation scaling factor from a targeted process-based simulation i. The new evidence ηNi can then be used to update the posterior for η, which propagates through to the final posterior for β. We model this as

ηi=η+ϵproci.(30)

Figure 6 shows the impact of updating the posteriors using Equation 30 with data from two process-based model simulations: ACCESS-ESM 1.1 (Ziehn et al., 2021) (η = 0.83) and UKESM1-0-LL (Liddicoat et al., 2026) (η = 0.74). The existence of these simulations means we are faced with yet another choice: how representative of the “real” η_{_*N*} do we expect them to be? Once more, we must determine – as with the CMIP6 results – whether these estimates are likely to be biased, whether those biases are ESM-specific or structural, and whether we can attribute them to specific aspects of each ESM – for example, the presence or absence of submodels. Because we have only two simulations, we choose here to model the spread in process-based simulation scaling factors as

ϵproci~N(0,0.1).(31)

That is, in Equation 31 we assume a known spread σ = 0.1, no common structural bias (despite the fact that these process-based results arise from two ESMs also participating in CMIP), and that more process-based simulations will invariably better constrain the real-world scaling factor. (If we had more simulations, we might put a prior on σ_η and learn the structural spread from the data.) The red distribution in Figure 6 shows the resulting updated estimate for η when combining information from CMIP and L (i.e., the difference between CMIP ESMs with and without nitrogen limitation) the two process-based runs. This updated estimate then propagates through to a new estimate of X, in this case to slightly higher values due to the updated η.

If the process-based runs are assumed to be highly certain and close to one another, they will strongly constrain the scaling factor η, which propagates through to the posterior for β. The purple distribution in Figure 6 shows the result of adding hypothetical simulation data: 1000 synthetic process-based runs in which the scaling factor η^(*m*) = 1 exactly. This certainty translates into a wider posterior for β. If nitrogen limitation does not modify β, the CMIP spread must be attributed to alternative mechanisms.

This framework simplifies the fact that ESMs may face structural difficulties in simply turning on/off processes such as nitrogen limitation, depending on how they are implemented. It also neglects the nuances of process definition. For example, an ESM may incorporate nitrogen interactions, but no limitation, and this may be inappropriate to simply consider a process “on” or “off”. Moreover, ESMs that do not contain the process may partially account for it when tuning to historical data. These nuances may be addressed through the careful design of factorial experiments or accommodated by modifying the interpretive model to handle additional complexity.

This approach lends itself to more granular decomposition. For example, we have explicitly modeled the effect of nitrogen limitation on CMIP simulations of β. However, other processes that can strongly affect land carbon uptake and storage are also present in selected CMIP ESMs: fire, for example, or dynamic vegetation. To account for these mechanisms, we require a new model, perhaps one in which multiple processes can scale the carbon effect up or down:

βm=∏pηpLmpUm+bm+bshared(32)

where

Lmp={1,if model m simulates process p0otherwise.

This model lets us update our posterior for β by incorporating information about multiple processes — from targeted simulations, experiments, or physical intuition about their effects. If we instead use the simpler model described in Section 3.5, we can only use information about nitrogen limitation, with no clear way to incorporate our knowledge about other processes.

5 Discussion and conclusions

Because we hope to update our understanding of a quantity X with observational constraints or direct evidence continually, why should we not simply set the prior on X to be N(μ, σ) with μ and σ the CMIP ensemble mean and standard deviation, respectively? At first glance, this appears to be a straightforward way to reconcile standard IPCC practice with the proposed Bayesian framework. However, this is description, not inference, and provides us with no coherent means of updating the posterior (for μ or σ). In practice, we want to use CMIP to generate priors on the quantity of interest X and then update them to the best of our ability with non-CMIP evidence: observations, paleoclimate reconstructions, emergent constraints, or simulations performed with other kinds of models. Why does the prior matter?

Sometimes it doesn't. The posterior for the global mean temperature anomaly ΔT relative to preindustrial levels is extremely well-constrained. This would hold for a wide range of priors, even extremely uninformative ones, and even ad hoc choices such as the frequently used CMIP multi-model mean and standard deviation N(μ, σ). This is because we have excellent observations of global mean temperature, including several independently developed datasets that agree strongly with one another over the historical period.

Conversely, in cases where the evidence is indirect or weak, the choice of prior can matter greatly. This is expected from Bayesian principles: the less observational evidence we have, the more we must rely on our priors and/or the CMIP archive for our conclusions. In such settings, our assumptions about the link between CMIP and the real world can be of immense scientific importance.

By clearly stating our assumptions, we can evaluate them and replace them with more defensible approaches when necessary. For example, the sequential updating of ECS with different emergent constraints relies on the assumption that the constraints themselves are independent. Because they are not, it is more scientifically defensible to explicitly model their correlations, resulting in a wider posterior. Clearly specifying the interpretive approach makes it easier to refine or replace it as methodological requirements dictate.

The framework is flexible enough to accommodate other forms of statistical dependence between ESMs: perturbed parameter ensembles using the same ESM, higher- and lower-resolution versions, different initial or boundary conditions. Knowledge about these experimental setups constitutes evidence, just as the lookup table L specifying which models contain certain processes does. We already implicitly use this knowledge in analyses – the Bayesian methods make it explicit.

In short, we argue that a Bayesian approach to assessment is principled because it makes all unavoidable subjective choices transparent and open to scrutiny. It is also the appropriate analytical framework, because posterior distributions can be systematically updated as more evidence becomes available.

In climate assessments, the interpretive framework chosen matters as much as the evidence itself. The evidence is generated by the wider scientific community, which produces satellite and in situ observations, process-based simulations, theoretical understanding, and CMIP (and other MIP) simulations. The task of the assessors is twofold: first, to determine which evidence should be included, and second, to design and specify the models used to interpret it. In doing so, assessors should decompose the model to the level where evidence exists or may exist in the future given planned experiments, missions, or data collection. It is important not to decompose further than necessary: extra granularity will add parameters without adding value. The goal is to balance flexibility and simplicity, and to present the analysis framework clearly to the research community. Once a model is chosen, expert assessment can then be solicited by requesting priors, allowing individual participants to articulate, specify, and defend their subjective assessments of the data and relevant variables. In some cases, prior specifications will exert minimal influence, as the likelihood dominates the posterior. But in others, where we have little evidence, expert judgment carries more weight.

Subjective decisions and expert judgment are present in all scientific analyses. In the Bayesian framework presented here, they enter in three locations:

• What constitutes evidence?
• What model should we choose to interpret the evidence?
• What priors should we use?

This article has focused on the role of interpretive models in CMIP. The choice of a model determines the parameters on which we must place priors. In addition, it determines how we incorporate evidence into our analyses. All of these decisions are subjective, and all are unavoidable in any analysis. Bayesian methods surface these choices by making them explicit.

Generative AI statement

The author(s) declared that Generative AI was used in the creation of this manuscript. Generative AI was used to fix capitalization errors in the BibTeX reference file and to add a list of DOIs to existing citations. All citations were human-checked to verify accuracy. Generative AI was also used to package human-written code (checking import paths and making all scripts pip-installable). No manuscript text was written by AI.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Andre L. Belem, Fluminense Federal University, Brazil

Reviewed by: Nathan Gillett, Canadian Centre for Climate Modelling and Analysis (CCCma), Canada

Jeremy Fyke, Environment and Climate Change Canada British Columbia, Canada