On the role of moderator and mediator variables in conditional event attribution of heatwaves over Europe

Introduction

In the aftermath of a heatwave, we often ask whether the event would have unfolded differently in a world without anthropogenic climate change. Starting with [29], event attribution studies have answered this question based on the changing probability distribution of, for example, extreme temperatures in some region of interest. In establishing such statistics, very different events will often be grouped together and treated as samples from one distribution. As an example we consider the summers of 2006 and 2007 in central Europe, both of which were approximately $1.5^\circ\mathrm{C}$ warmer than the climatological mean in the reference period 1961-1990. Figure 1 shows that the two summers, though with similarly large temperature anomalies, have very different “backstories”: 2006 saw a spring with reasonably high precipitation anomalies jointly with only slightly elevated North Atlantic sea surface temperatures (sst). In contrast, the spring of 2007 was exceptionally dry with considerable positive sst anomalies already present in the coastal waters. As we will see below, the two summers also had very different atmospheric circulations, specifically with respect to blocking highs. Due to the similar temperature anomalies, attribution studies based on temperature statistics alone will inevitably reach similar conclusions for both years. Nevertheless, it is possible that European summer temperatures with antecedent droughts or strong atmospheric blocking respond very differently to climate change.

To take such factors into account, we need a conditional attribution approach that asks, for example, “which temperatures would we have observed if the spring and blocking conditions of 2007 had occurred in a pre-industrial climate?” Such questions are addressed by the so-called storyline [32, 25, 26] and analogue [34] 2007] methods of attribution. For a recent overview of different conditional attribution methods, we refer to [31]; their common thread is that the properties of an event are studied conditional on some aspects of the climate state that embedded the event in space and time.

The introduction of storylines, in particular, sparked a contentious debate in the scientific literature [19], which is also of interest from the perspective of the philosophy of science [20, 7]. Following this debate, the current best practice is to clearly state which form of conditionality was used, if any, and to combine multiple attribution approaches to gain a more comprehensive understanding of the event. This is advocated by [31] and put into practice by [27], [1], [18], [22], [5], [3] and others. These studies present attribution results with different degrees of conditioning and discuss their similarities and differences on a case-by-case basis.

However, conditioning can have a significant impact on the outcome of the attribution statement, depending on the specific role a conditioning variable plays within the causal chain in the climate system. Figure 2 shows three archetypal roles that a variable can play. First, a variable can function as a covariate (panel a) within the climate system without itself being influenced by climate change or altering the response of the target variable or event to climate change. The second role shown in panel (b) is that of a mediator [12], through which climate change influences the target variable or event. The conditioning variable represents a process that establishes the causal link between changing greenhouse gas concentrations and the target variable or event. In a third role (c), a variable can act as a moderator that modulates the effect of anthropogenic emissions without representing a necessary condition of the event.

In this article, these three archetypes serve as simplified conceptual models to understand the effects of different conditions on a fundamental level. We will address the following research questions:

1. What is the impact of conditioning on a variable that accounts for a significant portion of the variability in the event of interest?

2. How do the results change depending on the causal relationship between event, conditions and forcing?

As a concrete example, we also ask

3. Did climate change affect the summers of 2006 and 2007 differently, if we take their different stories into account?

To address these research questions systematically, we use an intermediate approach that incorporates conditioning variables into the traditional statistical event attribution framework (Sect. 2). For the attribution of the summers 2006 and 2007, we use random effects model (Sect. 3) that allows us to include and combine information from multiple climate models and cover effects of uncertain observations or event definitions. Section 4 summarises the data used in our experiment before we attribute the two summers with and without conditions in Sect. 5. The results are discussed in terms of possible causal relationships in Sect. 6, where we consider a simplified “toy”-model of mediator and moderator conditions. Our main concluding remarks are given in Sect. 7.

Risk-based conditional attribution

This section and Sect. 3 introduce the theoretical underpinnings and practical implementation of risk-based conditional attribution. A list of the main mathematical symbols used is given in Table 3 in Appendix A.

PR	probability ratio
LR	likelihood ratio (ratio of probability densities)
$S_1, S_0$	scenario with and without climate change
$\mathbf{X}, \mathbf{Y}$	climate state vector in observations and models
$\mu_{t,i}$	mean climate state time $t$ under scenario $S_i$
$\mathbf{V}$	natural variability
$\epsilon$	random observation error
$b_{t,i}$	climate model bias at time $t$ under scenario $S_i$
$\boldsymbol{\Sigma}_{V,i}$	covariance matrix of natural variability under $S_i$
$\boldsymbol{\Sigma}_j$	covariance matrix of natural variability for climate model $j$
$\boldsymbol{\Sigma}_{b,t,i}$	cov. matrix of model biases at time $t$ under $S_i$
$\hat{\boldsymbol{\Sigma}}_{\mu,t,i}$	error cov. matrix of the estimated mean state at time $t$ under $S_i$
$\boldsymbol{\Sigma}_{\epsilon}$	observation error cov. matrix

Risk-based attribution studies estimate the probability ratio (or risk ratio) for an event $E$

where $S_1$ and $S_0$ are scenarios representing a world with and without climate change, respectively. For example, $E$ could be defined as an observed temperature $X$ exceeding some high threshold $u$. A probability ratio $\operatorname{PR}(X > u) = 2$ would then indicate that the risk of temperatures at least as extreme as $u$ has doubled due to climate change.

We can limit the scope of our attribution study by studying the probability of $E$ conditional on some set of conditions $C$. The outcome is a conditional probability ratio

For instance, if $E$ is again a temperature extreme and $C$ describes the presence of atmospheric blocking, then $\operatorname{PR}(E \mid C) = 3$ would imply that the probability of reaching such extreme temperatures under blocking conditions in both scenarios has tripled due to climate change.

By definition, conditional and unconditional probabilities are linked via $\mathbb{P}(A \mid B) = \mathbb{P}(A, B)/\mathbb{P}(B)$. Inserting this in (2), we find a simple expression for the conditional probability ratio:

This is a Bayes theorem of probability ratios, which states that the conditional PR is given by the ratio of the joint PR of $E$ and $C$ and the PR of the conditions $C$ alone. From this result, we can draw some immediate conclusions: firstly, if $E$ and $C$ are independent, we can factorise $\operatorname{PR}(E, C) = \operatorname{PR}(E)\operatorname{PR}(C)$, resulting in $\operatorname{PR}(E \mid C) = \operatorname{PR}(E)$. In the general case, where the conditions affect the event, the conditional $\operatorname{PR}(E \mid C)$ will differ from $\operatorname{PR}(E)$. When the conditions themselves experience no climate change, we have $\operatorname{PR}(C) = 1$ and the conditional $\operatorname{PR}(E \mid C)$ is equal to the joint $\operatorname{PR}(E, C)$ of event and conditions. Otherwise, the conditional $\operatorname{PR}(E \mid C)$ is smaller than its joint counterpart when $\operatorname{PR}(C) > 1$, i.e., when the conditions themselves become more likely under climate change.

Attribution method

For a real-world case, both the event $E$ and the conditions $C$ depend on the climate state random vector $\mathbf{X}$, which could, in general, contain any number of physical variables at various points in space and time. We will describe the observed climate in year $t$ under scenario $S_i$ as

where $\mu_{t,i}$ is the mean state, $V_i$ denotes the natural variability and $\epsilon$ represents the observational error. Here, we have included the simplifying assumption that $V_i$ depends only on the scenario and not on time, while the observation error $\epsilon$ depends on neither of those factors.

Since we have only one observation of $X_{t,i}$ in the factual world $S_1$, and no observations of the counterfactual world $S_0$, we will estimate the necessary probabilities from ensembles of climate models. The random vector of their modelled climate state for year $t$ under scenario $i$ will be denoted as $Y_{t,i,j}$, where the additional index $j$ enumerates the various climate models. In the model world, there are no uncertain measurements, but the model mean can differ from $\mu_{t,i}$ by a bias $b_{t,i}$, i.e.,

When multiple model ensembles are considered, each model $j$ is assumed to have its own realisation of the bias, which could depend on time and scenario as well. Each ensemble member in turn contains an independent realisation of the natural variability $V_{i,j}$ for that model and scenario. Equation (5) then corresponds to a random effects model (see for example Von Storch and Zwiers [2002] chapter 9.2.8). [33]

In general, $X$ and $Y$ could contain the full state of the climate in space and time. To make the problem tractable and easier to interpret, we consider four-dimensional climate states

where the components represent northern-hemispheric summer temperatures, summer blocking, spring precipitation and spring sst for a particular year. Exact definitions are given in Sect. 4.

Event definition

We can formally relate an event $E$ to the observed climate state $X$ via

where $\Omega$ is the set of all possible states. The first factor in the integral contains the event definition, which is the same under both scenarios. The second term is the probability density function (pdf) of the observations we would make under scenario $S_i$.

Univariate attribution studies often define the event as a threshold exceedance (i.e., as a Heaviside step function)

Inserted into a one-dimensional version of (7), this simply results in

where $F$ is the cumulative distribution function (cdf) of $X$.

To compute conditional probability ratios in (3), we need the joint probability of event and conditions, as well as the probability for the conditions alone. Defining these in terms of component-wise threshold exceedances is not generally appropriate: in our example case in (6), the second component of $X$ represents a blocking frequency. Defining $\mathcal{C}$ as a threshold exceedance would mean conditioning on all situations at least as strongly blocked as during the event, which is not what we want. To avoid this complication, we formally replace the probability in the event definition by a Dirac delta function at the observed value $x_o$:

(7) then simplifies to $f(x_o\mid S_i)$ and the probability ratio becomes a ratio of (multivariate) probability densities. To avoid confusion and abuse of notation, we will refer to this as quantity as the likelihood ratio

Throughout this manuscript, we will summarise our attribution results in terms of LR, which is the natural extension of PR to events that are not threshold exceedance.

Gaussian random effects model

To evaluate $f(\mathbf{x})$, we describe the natural variability $V_i$ in (4) as a zero mean Gaussian random variable with covariance matrix $\Sigma_{V,i}$. For the mean for each year and scenario, we insert an estimator $\hat{\mu}_{t,i}$ from our climate models. This estimator is also a Gaussian random variable with error covariance matrix $\Sigma_{\hat{\mu},t,i}$ depending on the year and scenario. Lastly, the observation errors $\epsilon$ are assumed to be zero-mean Gaussian with covariance $\Sigma_{\epsilon}$. Assuming that variability, bias and observation error are independent of each other, the desired density is then multivariate normal

The four unknown parameters in (12) are estimated as follows: for each climate model $j = 1,\ldots,J$, year $t$ and scenario $i$, we obtain the ensemble means $\bar{\mu}_{t,i,j}$. For each model, the natural variability covariance matrix $\Sigma_j$ is estimated from a long pre-industrial control run. From these quantities, we estimate $\hat{\mu}_{t,i}$ and $\Sigma_{\hat{\mu},t,i}$ via maximum likelihood using the mixmeta software package [24]. Following [8], the log-likelihood corresponding to the model in (5) reads

where $\Sigma_{b,t,i}$ is the unknown covariance matrix of the model bias $b_{t,i}$. Due to the small number of climate models used (see Table 2), we limit $\Sigma_{b,t,i}$ to have a diagonal structure, neglecting the correlations between biases in different variables. (13) is then iteratively optimised with respect to $\hat{\mu}_{t,i}$ and $\Sigma_{b,t,i}$. The error covariance of $\hat{\mu}_{t,i}$ is given by

model	hist-nat	historical	piControl (# years)
canesm5	42	35	451
cesm2	1	10	1200
cnrm-cm6-1	3	23	213
hadgem3-gc31-ll	4	4	500
ipsl-cm6a-lr	6	23	115
miroc6	3	10	500
mri-esm2-0	5	5	200

To evaluate (12), we also need an estimate of the observed natural variability covariance matrix $\Sigma_{V,i}$, which is not a parameter of the random effects model. We can define a natural estimator as the constant covariance matrix $\hat{\Sigma}_{V,i}$ by which all $\Sigma_j$ could be replaced such that (14) is fulfilled. Thus inserting $\hat{\Sigma}_{V,i}$ into (14), we have

which can be re-arranged to obtain

This is inserted as our best estimate of $\Sigma_{V,i}$ in (12). To check that the covariance matrix defined in this way is positive semi-definite (and thus valid), we have to verify that

From (14), we recognise that $J\Sigma_{\hat{\mu},t,i}$ is the harmonic mean over matrices $(\Sigma_j + \Sigma_{b,t,i})_{j=1,\ldots,J}$. Since each $\Sigma_j$ is positive semi-definite by definition, (17) just expresses that fact that the harmonic mean (like any other mean) is monotone increasing in the sense of the Loewner order: adding some positive semi-definite $\Sigma_j$ to each element will not decrease the mean [2].

We note that

1. If the bias term is negligible, (16) simplifies to the arithmetic mean of the individual covariance matrices $\Sigma_j$.

2. Due to $\Sigma_{b,t,i}$, the resulting covariance matrix of $X$ can differ between scenarios and years, even though we assume a single, fixed $\Sigma_j$ per model.

3. The number of ensemble members is intentionally not used in estimating the overall mean and covariance. We thereby avoid giving extra credence to models from centres with more computational resources (“ensemble of opportunity” problem, see [30]). In particular, an unequal number of members per scenario (see Table 2) would artificially modify the probability ratios.

Lastly, we need an estimate of the observation error matrix $\Sigma_{\epsilon}$. We obtain a rough estimate of this quantity from the discrepancies between four global reanalysis datasets (see Sect. 4). Because of this small sample size, we cannot rely on the naive estimates for any individual year. Instead, we compute the six unique time series of differences (as a pre-whitening operation) between pairs of datasets and estimate $\Sigma_{\epsilon}$ from the pooled sample of these differences (analogous to the ensemble dressing approach in [23]).

With all terms in (12) estimated, we can evaluate the individual and joint likelihood ratio LR for our variables of interest. According to (3), the conditional LRs can then simply be computed as the ratio between joint and unconditional LRs.

Bootstrap confidence intervals

To contextualise our estimates of conditional and unconditional LR, we need to assess their uncertainty resulting from the randomly realised model biases and member states. To this end, we perform a parametric bootstrap based on the model given in (5): for each year $t$, scenario $i$ and model $j$, we draw a realisation of $b_{t,i}$ using the estimated $\Sigma_{b,t,i}$. Then, the same number of members as in the original ensemble (see Table 2 below) is sampled from $\mathcal{N}(\cdot|b_{t,i},\Sigma_j)$. With these resampled ensembles, we repeat the same estimation procedure as before to obtain 100 bootstrap samples of the various LRs.

Data

As mentioned in the introduction, our real-world example revisits the classic case of European summer temperatures, studied in the very first event attribution study [29]. More specifically, we consider spatial means of daily mean temperature in the IPCC AR6 region 17 (central Europe, see Figure 1), averaged over the Northern Hemisphere (NH) summer months June to August as our event variable. To represent the dynamical conditions during each summer, we include the mean frequency of blocking or subtropical ridge activity, averaged over the same region and timespan. Some details about the blocking definition following Sousa et al. [28] are given in Appendix B. Two further variables are added to represent the pre-existing conditions at the beginning of each summer season: first, the central European precipitation total during the preceding NH spring (March to May), and second, the NH spring mean sst averaged over the European region (Figure 1). For 2006 and 2007, we have presented these spring conditions in Figure 1. All variables, their regions, units and abbreviations are summarised in Table 1.

	full name	spatial aggregation	time aggregation	unit
tas	daily mean temperature	AR6 reg. 17	JJA mean	°C
bl	presence of block or ridge	AR6 reg. 17	JJA sum	days / summer
pr	precipitation	AR6 reg. 17	MAM mean	$mm$ / month
sst	sea surface temperature	AR6 reg. 16,17,19	MAM mean	°C

Since we are interested in the interplay between these variables, we need a representation of the observed state in which all four are physically consistent. We therefore employ ERA5 reanalysis data [13] as our best estimate of the observed NH spring and summer conditions. While ERA5 does provide a 10-member ensemble used during data assimilation, the variance across members is miniscule in recent years for the large-scale averages we study (not shown). Since further uncertainty may result from the assimilation system and choice of observations, we instead obtain additional data from three global reanalyses: NCEP2 [16], MERRA2 [9] and JRA3Q [17]. The error covariance matrix is estimated from the common time period 1980-2024, following the approach described in Sect. 3.2.

For the climate state under scenarios $S_{0}$ and $S_{1}$, we rely on CMIP6 simulations of the historic and hist-nat scenario, respectively [6, 10]. Our choice of models is mainly limited by the need for daily geopotential height fields for the blocking and ridge algorithm. Table 2 lists the included models and the used number of members from each scenario. In addition to these scenario runs, we use the respective piControl simulations to compute the natural variability covariance matrix of each model $\Sigma_{j}$.

All variables are converted to anomalies by subtracting the average values from the climate reference period 1961-1990 for each model and scenario as well as the reanalyses, separately. To first order, this step effectively removes the model biases in the event and especially in the conditions, and focuses our analysis on the climate change signals alone.

Figure 3 displays the ERA5 anomalies of sst, bl and tas. As expected, anomalously hot NH summers in Europe are associated both with strong blocking and positive antecedent sst anomalies. Many of the most recent years show positive anomalies in both tas and sst, indicating an observed climate change signal in both variables. In contrast, a shift towards stronger or weaker blocking is not readily apparent. The two summers of 2006 and 2007 to be discussed are encircled in Figure 3.

Attribution of two European summers

Figure 3 shows the strong contrast between our two selected NH summers from Figure 1: the summer of 2006 (tas anomaly of $+1.4\,^{\circ}\mathrm{C}$) has one of the strongest blocking anomalies on record (ten days more than average), with relatively low antecedent sst and above average spring precipitation. In contrast, the following summer of 2007 was slightly warmer than 2006, but without anomalous blocking activity. It was preceded by a relatively dry spring with exceptionally high regional sst anomalies.

For both years the ensembles means and piControl covariance matrices of our models are combined by the random-effects model described in Sect. 3. Before discussing the resulting distributions, we can look at the correlation matrices corresponding to our estimates of natural variability $\hat{\Sigma}_{V,i}$. For the historical scenario in 2006, we find

As expected, NH summer blocking is a strong predictor of temperatures in central Europe, explaining over half of the variance. Precipitation shows a slight anti-correlation with tas, reflecting a weak link between dry NH springs and hot NH summers. The correlation between tas and sst has a slightly larger magnitude and opposite sign, indicating that pre-existing positive sst anomalies in the European oceans favour a warm NH summer in central Europe. We find almost identical matrices for the year 2007, as well as the hist-nat scenario in both years. While this is not guaranteed to be the case, it is hardly surprising since our computation of anomalies has removed much of the inter-model variance, resulting in near-zero between model covariances.

We now turn our attention to the resulting conditional and unconditional temperature distributions for the two years shown in Figure 4. To better understand the effects of the three condition variables, Figure 5 shows their distributions together with the ERA5 values in the two years.

For 2006, Figure 4 shows that the observed temperature value is close to the unconditional expectation under the historical scenario (red), but far in the tail of the hist-nat distribution (blue). In terms of the likelihood ratio, $LR(\mathrm{tas}) \gg1$. This result is flipped completely when we condition the distribution on blocking: firstly, $f(\mathrm{tas}|\mathrm{bl})$ is considerably narrower than the unconditional pdf. In addition, the highly unusual blocking activity (Figure 5 a) shifts the distribution under both scenarios to higher temperatures. As a result, the 2006 value is actually more typical for hist-nat than for the historical scenario, which implies $LR(\mathrm{tas}|\mathrm{bl}) < 1$!

The positive precipitation anomaly (Figure 5 b) has hardly any impact on the temperature distribution, but sst shows a significant effect: due to the strong sst warming (Figure 5 c), the small observed anomaly is rather unusual for the historical scenario (panel d). Consequently the distribution of $tas$ conditional on $sst$ is shifted to lower values for the hist scenario. If we consider this condition in isolation, it also reduces the likelihood ratio, but not as strongly as blocking. However, when we combine all three conditions (top row of Figure 4), the effects of strong blocking and lower $sst$ compensate each other. As a result, the 2006 temperature is close to the intersection of the conditional pdfs, i.e.,

For the NH summer 2007 the attribution result is very different (right part of Figure 4). Here, conditioning on blocking also reduces the variance of the temperature distribution, but does not lead to strong shift since the observed blocking value is close to the climatological mean (Figure 5 a). As a result, the difference between the two pdfs is massively enhanced, resulting in a much stronger conditional attribution statement (Figure 7 below).

The effect of precipitation conditions is again hardly visible. While the anomaly has a similar magnitude but opposite sign from 2007, Figure 5 b appears to indicate a slightly stronger climate change signal in the precipitation distribution itself. However after inspection of the underlying model timeseries, we interpret this as a random fluctuation in the hist-nat simulation. Due to the relatively small sample size for that scenario (Table 2) and the inherent uncertainty of precipitation in coarsely resolved climate simulations, we will not interpret this phenomenon any further.

The rather high $sst$ anomaly in 2007 is close to normal for the historical scenario but far too large for hist-nat (Figure 5 c). The corresponding conditional temperature distribution $f(tas \mid sst)$ is therefore shifted towards its historical counterpart (reduced LR). When we condition $tas$ on the “whole story” (top right curves in Figure 4), the effects again compensate slightly. However, the variance reduction is clearly dominant, meaning

The outcome of our conditional and unconditional attribution experiments for 2006 and 2007 is quantitatively summarised in terms of likelihood ratios in Figure 6. We observe that all conditional and unconditional LRs are exponential functions of the temperature anomaly (linear on the logarithmic axis). The unconditional curves are nearly identical for the two years, resulting in LR $\approx10$ for the observed values. The curves for the conditional LR are steeper and have opposite offsets for the two years. As a result, the confidence interval for 2006 includes LR $= 1$, while the conditional LR lies around 100 in 2007. Hence, when we control for the dynamics and the precoditions in 2006, there is no evidence for an increased event likelihood. In 2007, conditioning on the same variables leads to an increase in LR by orders of magnitude.

For completeness, Figure 7 summarises the LR for the observed values given all or no conditions (as in Figure 6) and conditional on each variable individually. As expected from the previous discussion, we see a massive effect from the reduced variance when conditioning on bl. This is slightly compensated by the other two conditions, with sst in particular resulting in nearly identical LR for both years.

Lastly, Figure 7 also shows the effect of the observation errors $\Sigma_{\epsilon}$. The difference in attribution outcomes with and without uncertain observations (triangles and circles in Figure 7) is very limited in our example. The overall effect, seen most clearly in the extremely uncertain LR(tas|bl), tends to bring the ratios closer to unity by ever so slightly reducing the signal to noise ratio. We expect this effect to become much more relevant for smaller spatio-temporal scales and data-sparse regions outside of Europe.

Role of the conditioning variable

We now want to understand and classify the effects of conditioning seen in Sect.5 in terms of the different archetypal causal relations, denoted as pure covariate, mediator, and moderator (Figure 2). To this end, we consider a simple model for a univariate target variable $X$ and a univariate conditioning variable $M$ with the following assumptions:

• All process variables are normally distributed with variances independent of the scenario.

• $M$ has an additive effect on $X$, which may vary depending on the scenario.

• Other processes affecting $X$ are represented by a random variable $\varepsilon$ which is independent of $M$.

The corresponding model is represented by

where $\varepsilon$ is a zero-mean Gaussian random variable, $i \in\{0,1\}$ again denotes the scenario $S_i$, and $\mu_X$, $\Delta x$, $\beta_1$, and $\beta_2$ are constant parameters.

represents our condition variable, where $\eta$ follows zero-mean Gaussian distribution independent of $\epsilon$, and $\Delta m$ is a constant climate change signal in $M$. For simplicity, we have also assumed that $M$ has mean zero in the counterfactual scenario. It follows, that both the unconditional distribution of $X_i$ and the conditional distribution $X_i|M_i = m$ are also Gaussian, and the results of our investigation are fully determined by the conditional and unconditional means and variances. They are given by

With this simplified model, we now study the three archetypal causal relationships between climate change, event and conditions, shown in Figure 2. For these examples, it is convenient to keep the unconditional distributions at the same location in all three cases, thus mirroring the situation in Sect. 5 (i.e., one unconditional distribution, and multiple possible conditions). This is achieved by re-parametrising the model in terms of the unconditional expectations $E[X_0]$, $E[X_1]$ instead of $\mu_x$, $\Delta x$. The conditional expectation values then read

where we have tacitly inserted $i^2=i$ since $i \in\{0,1\}$. For our example distributions shown in Figure 8, we set $E[X_0]=0$, $E[X_1]=2$, $\operatorname{Var}[\varepsilon]=\operatorname{Var}[\eta]=1$ and consider a condition value of $m=1$. The parameters $\beta_1$, $\beta_2$, $\Delta m$ are chosen to reflect the different causal situations:

• For the pure covariate, $\beta_1=1$, $\beta_2=0$, $\Delta m=0$.

• For the pure mediator, $\beta_1=1$, $\beta_2=0$, $\Delta m=1$.

• For the pure moderator, $\beta_1=0$, $\beta_2=1$, $\Delta m=0$.

Conditions act as a pure covariate

In this setting, the climate change scenario $S_i$ and the conditions $M$ act on $X$ completely independently (Figure 2 a). In our simple model this is realised if there is no climate change signal in the conditions ($\Delta m=0$) and no interaction term between $S_i$ and $M$ (thus $\beta_2=0$). Consequently, the conditional distributions are shifted by $\beta_1m$ in both scenarios, and the conditioning reduces the variance by $\beta_1^2\operatorname{Var}[\eta]$. The resulting pdfs are shown in Figure 8.

This is approximately the situation we observed in Sect. 5 when conditioning tas on blocking activity, which is largely unaffected by climate change: the distribution under both scenarios becomes narrower and receives a shift. Since in our toy model the variance remains unaffected by climate change, we have a simple equation for the log-likelihood ratio of the conditional probability densities (see Appendix C), consisting of three components:

We note that the unconditional $\operatorname{LR}(x)$ is recovered by $\beta_1=0$. The first term $(i)$ in (28) represents the anomaly of $x$ with respect to the conditional mean averaged over the two scenarios. As we have seen in Figure 6, $\log\operatorname{LR}$ scales linearly with $x$ and is shifted by $m$. The second factor $(ii)$ is the climate change signal of $X$. Factor $(iii)$ describes the reduction of $X$’s variance due to conditioning on $M$, which is always positive. If we assume a positive shift $\Delta x>0$ due to climate change, then $\operatorname{LR}>1$ holds if and only if the event under investigation $x$ is unusually large even conditioned on $m$. This explains how conditioning on $M$ can modify and even flip the attribution statement. For example if we control for very favourable conditions $\beta_1m\gg0$, we may find that $x$ was actually unusually small for the climate change scenario $S_1$ and more typical of the counterfactual climate state $S_0$. (28) also explains the general increase in LR we observe when conditioning on a strong predictor like blocking: the larger $\beta_1$, the larger is the part of the variance explained by the conditions, leading to a larger LR due to term (iii).

Conditions act as a pure mediator

If a variable is acting as a mediator, it is causally located between $S$ and $X$ (see Figure 2b) and forms part of the causal pathway through which $X$ is changing [14]. In our simplified case, this implies a non-zero climate change signal in $M$ ($\Delta m \ne0$). If $M$ is purely a mediator, the interaction term is absent ($\beta_2 = 0$).

The reduction of the variances is then the same as for a pure covariate as defined in Sect. 6.1. The changes in conditional and unconditional expectation value ((27)) are now related via

Assuming, for simplicity, that $\beta_1 \Delta m > 0$, i.e., conditions shift in a way that favours large $X$, the conditional climate change signal is smaller than the unconditional one, because we have removed the part of the effect that is mediated by the conditioning variable $M$.

The likelihood ratio from (28) is modified in two opposing ways: in term (i), $m$ is replaced by $m - \Delta m/2$, accounting for the fact that $E[M]$ now depends on the scenario. For $\beta_1 \Delta m > 0$ this effect would increase LR. Conversely, term (ii) is reduced by $\beta_1 \Delta m$, thus reducing LR. For very unusual conditions ($m \gg\Delta m$), the effect in (i) can be neglected and the net effect of conditioning on a mediator will be a reduction of LR. In Sect. 5, we have observed a similar effect, where part of the tas signal was removed by conditioning on sst. However, it is important to note that this analysis cannot determine whether sst actually mediates the relationship between climate change and central European summer temperatures. This distinction will be discussed in Sect. 6.4 below.

Conditions act as a pure moderator

While a mediator is part of the physical mechanism through which climate change modifies $X$, a moderator modulates the strength of the effect (Figure 2c). A simple example is when the climate change signal exhibits a seasonal cycle, in which case the season acts as a moderator of climate change. In our model ((21)), this is represented by an interaction term between $S_i$ and $M$, i.e., $\beta_2 \ne0$. In the simplest case, shown in Figure 2c), all other causal connections are absent, i.e., $\Delta m = 0$, $\beta_1 = 0$. The conditional and unconditional expectation changes are then related as

In addition, the moderator introduces a scenario dependence to the unconditional variance

which is removed, when we condition on $M$. These two effects are shown in the top row of Figure 8. Here, conditioning on a moderator increases the separation between distributions and affects the variances differently under the two scenarios, because $M$ contributes more variance under $S_1$.

While changing variances between scenarios were theoretically possible in our experiments (through the between model covariance, see Sect. 3.2), they did not occur in Sect. 5. This suggests that possible moderation effects were neglected in our analysis. They might be represented in a more sophisticated statistical model, where the natural variability covariance matrix is permitted to differ between scenarios.

Correlation vs Causation

A statistical model like (21) does not encode any information on the actual causal relations between the processes it represents. The parameter regimes we discussed above are consistent with the three causal pictures from Figure 2, but they would be equally consistent with different, possibly more complex situations, where for example $M$ and $X$ have other common drivers, or some of the directions of causality are reversed. Based on our results in Sect. 5 alone, we cannot determine whether there is a physical pathway linking NH spring sst to NH summer temperatures, or whether the link is purely statistical.

In a strict sense, the true causal structure can only be determined by an experiment in which the value of the conditions is prescribed. In the notation of causal theory proposed by Pearl [21], such experimentally constrained probabilities are denoted as $\mathbb{P}(X \mid do(m))$, which is generally not equal to $\mathbb{P}(X \mid m)$. While a controlled experiment is impossible with the real climate system, it can be realised in silico using weather or climate model simulations. For the effect of the scenario forcing $S$ itself, this was indeed realised by the historical and hist-nat simulations we used. The conditions under which these experiments allow for a causal interpretation are clarified in [11]. In principle, the analogous experiment can be run for the effect of the weather or seasonal conditions $M$, for example by fixing the large-scale flow through data assimilation. This is the principle of the storyline school of climate change attribution [26]. While an in-depth discussion of the different styles of conditional attribution is beyond the scope of this study, we can already see that this experimental style of conditioning can easily have effects that differ drastically from purely statistical conditioning: fixing a condition with no causal link from $M$ to $X$ will have no effect in a controlled storyline experiment. Conversely, it may heavily influence the outcome of our conditional statistics. This is particularly relevant in light of the third popular style of conditional attribution –flow analogues– which also relies on statistical rather than experimental conditions.

Concluding remarks

Many recent attribution studies have begun to incorporate several complementary methodologies with varying kinds of conditioning. In this study, we have asked for the systematic effects of including such conditions. The NH summers of 2006 and 2007 serve as a striking example case: when we control for the NH summer blocking and, to a lesser degree, the conditions of the preceding NH spring sst, the probability ratio for tas in 2007 increases by orders of magnitude. In contrast, the summer of 2006 was not unusually warm once we account for the extremely favourable blocking circulation, leading to a conditional likelihood ratio close to 1. This illustrates that conditional and unconditional attribution approaches can differ systematically, and even lead to seemingly contradictory conclusions.

Specifically, we have seen two main effects of conditioning. Firstly, removing part of the variability through conditioning leads to a better distinction between scenarios, resulting in more confident attribution statements. This reflects one of the original motivations for the storyline approach [32], namely the ability to analyse very rare and thus highly uncertain phenomena. In our example, the blocking index of Sousa et al. [28] has proven to be a very strong predictor of NH summer temperatures in central Europe. By conditioning on this index, we explain around half of the natural variability and thereby greatly increase the separation between the two temperature distributions. Secondly, the conditions shift our expectations of what is “normal” under each scenario. In our example, the exceptional blocking in 2006 caused the temperature distributions in both climates to shift towards higher temperatures. Consequently, the observed anomaly was actually more typical of the counterfactual climate, rather than being attributed to climate change.

As discussed in Sect. 6, both of these effects can occur when our conditions act like a pure covariate, i.e., independently of climate change. The situation is significantly more complex when we allow for conditions that are involved in the causal chain, for example acting as moderators or mediators of climate change: conditioning on a mediator will “regress” out part of the mean signal. Conversely, a moderator condition may contribute additional variance under the factual scenario, which is removed when the condition is fixed.

Although our conditional attribution experiments were based on low-dimensional normal distributions, many of the results can be applied to the conditional attribution frameworks that are commonly used in practice. In the case of flow analogues, for example, the condition variable describes how close the observed dynamical state of the atmosphere is to the target state. Depending on the presence of dynamical trends, this could be a moderator or mediator. As discussed in Sect. 6.4, the situation is fundamentally different for so-called dynamical storylines [26], where the conditions are experimentally enforced. We have argued that such experiments could be used to test whether the link between event and conditions is in fact causal or purely statistical. In the presence of dynamical trends, we should therefore not be surprised if analogues and storylines produce different results.

In this study, we considered climate change in the expectation values of normal distributions only. While this assumption is reasonable for seasonal averages and large geographic regions, it will break down when we move to more localised and short-lived events, such as individual storms or extreme rainfall events. Although the general “Bayes theorem of probability ratios” ((3)) still holds, the conditioning may have further consequences beyond those considered here. One example is the moderator effect, which would manifest as changes in the covariance structure under climate change. In the real world, the effective causal structure may be more complex with variables acting as both moderators and mediators, for example as moderators on shorter timescales and as mediators on longer timescales. Evidence for such effects could potentially be collected by applying models like (21) to data in a moderation and mediation analysis [14]. Carefully controlled experiments similar to the storyline approach may then help in distinguishing between pure correlations an genuine causal pathways.

List of symbols

Blocking definition

We want to summarise the presence of atmospheric blocking by a single index which is strongly correlated with surface temperatures, can be localised over the region of interest and does not automatically contain a thermodynamic trend (for example due to the expansion of the atmosphere). These criteria are met by the relatively recent approach of Sousa et al. [28] based on daily 500 hPa geopotential height fields. Here, the effect of thermal expansion is removed by considering only spatial gradients and anomalies with respect to a short time window. A key novelty of this blocking definition is the identification of subtropical ridges in addition to classic Rex- or Omega-type blocks. Grid points are considered instantaneously blocked if either (a) the geostrophic flow in 500 hPa is reversed (Rex or Omega, similar to Davini et al. [4]) or (b) a ridge with geopotential height above the hemispheric mean in the preceding 15 days extends north of the edge of the subtropical belt. Contiguous blocking objects are identified in the daily fields, filtered for sizes $\ge500.000\ \mathrm{km}^2$ and then tracked over time to ensure a lifetime of at least four days. For the details of the algorithm, we refer to Sousa et al. [28]. The resulting daily binary blocking fields are then averaged over the study region and each summer to obtain a yearly timeseries of blocking activity.

Gaussian LR in case of constant variances

In general, the logarithm of the normal pdf under scenario $S_i$ can be written as

Assuming equal variances under the two scenarios, we find for log-likelihood ratio

The outcome of our attribution study is thus determined by the variance of the distribution, the climate change signal in the mean and the anomaly of the event with respect to the average mean state.

Noting that the conditional distribution of $X_i|M_i=m$ is normal, we obtain Eq. (28) for the pure covariate case by inserting the expressions for the conditional means and variances (Eqs. (27) and (26)) into Eq. (37).

Acknowledgements

This research was funded within the BMBF project ClimXtreme phase 1 under grant number 01LP1902A and ClimXtreme phase 2 under grant number 01LP2323A. This work used resources of the Deutsches Klimarechenzentrum (DKRZ) granted by its Scientific Steering Committee (WLA) under project ID 1152. We are especially grateful to Etor Lucio Eceiza at DKRZ for his invaluable technical support. We furthermore thank Pedro Sousa for advice on his blocking definition and Ieda Pscheidt-Willems for her contributions in the early phases of this work.