# Partitioning Uncertainty Components of an Incomplete Ensemble of Climate Projections Using Data Augmentation

Guillaume Evin, Benoit Hingray, Juliette Blanchet, Nicolas Eckert, Samuel Morin, Deborah Verfaillie

► To cite this version:

Guillaume Evin, Benoit Hingray, Juliette Blanchet, Nicolas Eckert, Samuel Morin, et al. Partitioning Uncertainty Components of an Incomplete Ensemble of Climate Projections Using Data Augmentation. Journal of Climate, 2019, 32 (8), pp.2423-2440. 10.1175/jcli-d-18-0606.1. hal-02327459

HAL Id: hal-02327459

https://hal.science/hal-02327459

Submitted on 31 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Partitioning Uncertainty Components of an Incomplete Ensemble of Climate Projections Using Data Augmentation

<figure><img src="image_2.png" /></figure>

## GUILLAUME EVIN

Univ. Grenoble Alpes, Irstea, UR ETGR, Grenoble, France

## BENOIT HINGRAY AND JULIETTE BLANCHET

Univ. Grenoble Alpes, CNRS, IRD, Grenoble INP, IGE, Grenoble, France

## NICOLAS ECKERT

Univ. Grenoble Alpes, Irstea, UR ETGR, Grenoble, France

## SAMUEL MORIN AND DEBORAH VERFAILLIE

Univ. Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM, CEN, Grenoble, France

(Manuscript received 13 September 2018, in final form 4 February 2019)

## ABSTRACT

The quantification of uncertainty sources in ensembles of climate projections obtained from combinations of different scenarios and climate and impact models is a key issue in climate impact studies. The small size of the ensembles of simulation chains and their incomplete sampling of scenario and climate model combinations makes the analysis difficult. In the popular single-time ANOVA approach for instance, a precise estimate of internal variability requires multiple members for each simulation chain (e.g., each emission scenario-climate model combination), but multiple members are typically available for a few chains only. In most ensembles also, a precise partition of model uncertainty components is not possible because the matrix of available scenario/models combinations is incomplete (i.e., projections are missing for many scenario-model combinations). The method we present here, based on data augmentation and Bayesian techniques, overcomes such limitations and makes the statistical analysis possible for single-member and incomplete ensembles. It provides unbiased estimates of climate change responses of all simulation chains and of all uncertainty variables. It additionally propagates uncertainty due to missing information in the estimates. This approach is illustrated for projections of regional precipitation and temperature for four mountain massifs in France. It is applicable for any kind of ensemble of climate projections, including those produced from ad hoc impact models.

## 1. Introduction

A critical issue in climate change impact studies is the assessment of uncertainties associated with future projections. A comprehensive quantification of the different sources of uncertainty is necessary in order to obtain the range of possible future changes, their significance,

and their value for adaptation. The identification of the most important uncertainty sources is also important for the allocation of research and development resources (Northrop and Chandler 2014).

Uncertainty sources include scenario uncertainty, model uncertainty, and climate internal variability. Scenario uncertainty corresponds to the uncertain future of greenhouse gas emissions. It is now usually evaluated with climate projections obtained for different representative concentration pathways (RCPs), obtained from various socioeconomic projections (van Vuuren et al. 2011). Model uncertainty is related to the imperfections of the models used to obtain the climate projections (climate models and subsequent impact

Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JCLI-D-18-0606.s1.

Corresponding author: Guillaume Evin, guillaume.evin@irstea.fr
models when sectorial declinations such as ecology, water resources, hydropower, and so on are considered). Indeed, many different models can be used to represent the Earth system, leading to different climate responses for the same forcing configuration. For a given emission scenario for instance, different climate responses are generally obtained at the scale of the planet from different general circulation models (GCMs) and, for a given scenario/GCM simulation, different climate responses are expected at a regional scale from different regional climate models (RCMs). In addition, for a given GCM/RCM simulation, an ensemble of impact models can be used to represent the uncertainty regarding sectorial processes such as the hydrological response at the catchment scale, snowpack evolution, etc. Model uncertainty is typically estimated using different simulation chains, combining for instance different climate, downscaling, and impact models (Stocker et al. 2013). Climate internal variability originates from the chaotic and nonlinear nature of the climate system (Deser et al. 2012). Unlike scenario and model uncertainty, which could potentially be reduced if our estimates of future emissions and our knowledge and representation of geophysical processes were improved, climate internal variability is purely stochastic and is irreducible (Hawkins and Sutton 2009; Lafaysse et al. 2014; Fatichi et al. 2016).

These different uncertainty sources are typically estimated from multiscenario multimodel ensembles (MMEs) of transient climate projections produced for some control and future periods. Various methods have been proposed for this, usually based on an analysis of variance (ANOVA). An important limitation of most existing ANOVA models applied to MMEs is due to the nature of climate projection datasets available for the analysis. Typically, given the high computational costs of GCMs and RCMs, a large number of possible scenario/GCM/RCM combinations may be missing from the MMEs. For instance, in the EURO-CORDEX dataset (Jacob et al. 2014; Kotlarski et al. 2014), 7 RCMs have been used to produce high-resolution projections over the European domain, from 5 CMIP5 GCMs and for 3 RCP scenarios (Taylor et al. 2012). From the $7 \times 5 \times 3 = 105$ possible scenario/GCM/RCM combinations, only 17 were actually available for the analysis of Jacob et al. (2014). The number is increasing progressively but the matrix of GCM/RCM combinations will probably not fill up entirely by the end of this generation of the EURO-CORDEX dataset.

As indicated above, different approaches based on ANOVA analyses have been proposed to deal with these MMEs. The single-time approach consists in applying the ANOVA on projections available for one given projection lead time (Hingray et al. 2007; Yip et al. 2011; Paeth et al. 2017), independently from the

projections available for other lead times. Provided that multiple members (i.e., replicates) are available for each simulation chain, the climate response of the chain is estimated from the multimember mean and its internal variability is estimated from the intermember dispersion. Model uncertainty components are then estimated from the dispersion between the climate responses of the different simulation chains. In the single-time approach, multiple runs are required for each simulation chain, for both the estimation of the chain's climate response and internal variability. However, in the great majority of available MMEs, multiple GCM runs are only available for a small number of simulation chains. This multiple-run constraint often leads scientists to discard single-run GCM projections from the analysis. Only 6 GCMs out of 32 GCMs, for instance, are considered for this reason in Bracegirdle et al. (2014). Another approach to circumvent this constraint is to avoid separating the variations of the climate variable in the climate response from the fluctuations due to internal variability (e.g., Giuntoli et al. 2015). In this case, however, the contributions of model uncertainty and internal variability to the total uncertainty cannot be evaluated anymore. Another option that has been considered in the literature is to increase the size of MMEs using weather generator models (e.g., Lafaysse et al. 2014; Fatichi et al. 2016).

Alternatively, in the time series approach, the temporal variation of the climate response of any given chain is assumed to be gradual and smooth, the high-frequency variations of the time series being due to internal variability alone (e.g., Hawkins and Sutton 2009; Hingray and Saïd 2014; Reintges et al. 2017). For each simulation chain, the climate response is estimated with a trend model (e.g., polynomial functions) and the internal variability of the chain, considered constant throughout the time series, is estimated by the variance of the deviations from the climate response. As for the single-time approach, for any projection lead time, model uncertainty components are estimated from the dispersion between the climate responses. Thanks to the discrimination of the raw projections between internal variability and climate response, the time series approach can be applied when only a single member is available for some (or even all) simulation chains. It is thus not subject to the multiple-run constraint. In many MMEs configurations, especially when MMEs are noisy or when a small number of members is available, estimates of uncertainty components obtained with a time series approach are expected to be more precise than with the single-time approach (Hingray et al. 2019).

For both single-time and time series approaches, classical ANOVA methods cannot be applied with missing data (i.e., when several scenario/GCM/RCM combinations are
not available). A drastic solution would be to retain only scenario/GCM/RCM runs for which all combinations are available. However, this would lead to a dramatic waste of information, reducing the dataset to very few combinations (e.g., 2 RCMs × 4 GCMs for the EURO-CORDEX dataset; see Verfaillie et al. 2018). Missing-data methods have been developed to reduce biases and increase the efficiency of the estimation in such cases. They use, for instance, an ad hoc data reconstruction algorithm (Déqué et al. 2007) or likelihood-based techniques (Little and Rubin 2014). Under the assumption that missing combinations follow the same decomposition of variance as the data actually at hand, Bayesian inference and data augmentation techniques (Tanner and Wong 1987) can be applied and present the advantage to propagate missing-data uncertainty. As an example, Tingley (2012) applies a Bayesian ANOVA method to infer climate anomalies and to handle missing spatial locations.

In the present study, we introduce a method allowing for the characterization of multiple uncertainty sources, namely scenario, GCM, and RCM uncertainties, as well as internal variability, from an incomplete dataset of climate projections. The impact of missing scenario/ GCM/RCM combinations is explicitly treated using a Bayesian ANOVA framework and data augmentation. We show how missing combinations influence the uncertainty of the estimated ANOVA effects, and we compare these estimates to those obtained with a simple empirical estimation approach. This approach, named Quasi-Ergodic Analysis of Climate Projections Using Data Augmentation (QUALYPSO), is illustrated using MMEs of precipitation and temperature projections in four different French mountain massifs located in the Pyrenees and the Alps mountain ranges, at an elevation of 1500 m (Verfaillie et al. 2017, 2018).

Section 2 presents the pretreatments applied to projections, and details the Bayesian framework used for the estimation of the climate response and for the partition of the different uncertainty components. An application of the QUALYPSO method to MMEs of precipitation and temperature projections is presented in section 3. Section 4 discusses the estimation uncertainty and presents a comparison of Bayesian and direct estimates. Section 5 summarizes the results and provides some perspectives.

## 2. Methodology

We consider a given MME composed of different simulation chains where each chain corresponds to a given GCM/RCM combination for a given emission scenario. Let us note $Y_{i,j,k}(t)$, the raw climate projections of the studied climate variable (e.g., temperature, precipitation) obtained with a given chain for a given year $t$, where $i$ is the

index related to the emission scenario, $j$ is the index related to the GCM, and $k$ is the index related to the RCM.

### a. Estimation of the climate response

For any chain, the time evolution of the climate response to greenhouse gas emissions is assumed to be gradual and smooth, the higher-frequency variations of the time series being due to internal variability of the studied variable alone (the so-called quasi-ergodicity assumption; see Hingray and Saïd 2014). The climate response of a particular simulation chain for a given emission scenario is then estimated from the long-term trend of the climate projections. The internal variability corresponds to the variance of the deviations of projections from the estimated climate response.

For any chain and each year $t$, the raw outputs of the studied climate variable can be expressed as

$$Y_{i,j,k}(t) = \phi_{i,j,k}(t) + \eta_{i,j,k}(t), \quad (1)$$

where $\phi_{i,j,k}(t)$ is the climate response and $\eta_{i,j,k}(t)$ is the deviation from the climate response as a result of internal variability.

The climate response $\phi_{i,j,k}(t)$ is often obtained by fitting an analytical trend model to the raw data $Y_{i,j,k}(t)$. A linear trend or second-order to fourth-order polynomial trend models are for instance used by Hawkins and Sutton (2009), Hingray and Saïd (2014), Bracegirdle et al. (2014), and Vidal et al. (2016). In what follows, the climate response of each chain is obtained using cubic smoothing splines (see section 3b).

### b. Estimation of the climate change response

As in most climate impact studies, uncertainty sources are quantified from change variables, obtained as the difference of climate projections between a future and a reference period. The change variable can be defined in terms of absolute changes (e.g., for temperature)

$$Y_{i,j,k}^*(t) = Y_{i,j,k}(t) - \phi_{i,j,k}(c), \quad (2)$$

or in terms of relative changes (e.g., for precipitation)

$$Y_{i,j,k}^*(t) = Y_{i,j,k}(t)/\phi_{i,j,k}(c) - 1, \quad (3)$$

where $\phi_{i,j,k}(c)$ is the value of the climate response estimated from the trend model for the control year $c$ (e.g., $c = 1990$ in what follows).

The change variable $Y_{i,j,k}^*(t)$ can be split up into

$$Y_{i,j,k}^*(t) = \phi_{i,j,k}^*(t) + \eta_{i,j,k}^*(t), \quad (4)$$

where $\phi_{i,j,k}^*(t)$ is the climate change response of the scenario/GCM/RCM combination and $\eta_{i,j,k}^*(t)$ is the
deviation from the climate change response for this scenario/GCM/RCM combination, as a result of internal variability. When absolute changes are considered, these components read

$$ \phi_{i,j,k}^{*}(t) = \phi_{i,j,k}(t) - \phi_{i,j,k}(c), \quad (5) $$

$$ \eta_{i,j,k}^{*}(t) = Y_{i,j,k}(t) - \phi_{i,j,k}(t), $$

and when relative changes are considered they read as

$$ \phi_{i,j,k}^{*}(t) = \phi_{i,j,k}(t)/\phi_{i,j,k}(c) - 1, \quad (6) $$

$$ \eta_{i,j,k}^{*}(t) = [Y_{i,j,k}(t) - \phi_{i,j,k}(t)]/\phi_{i,j,k}(c). $$

### c. Decomposition of the climate change response

For each time t, we assume that the climate change response function $\phi_{i,j,k}^*(t)$ of any scenario/GCM/RCM combination can be expressed as

$$ \phi_{i,j,k}^{*}(t) = \mu(t) + \alpha_{i}(t) + \beta_{j}(t) + \gamma_{k}(t) + \xi_{i,j,k}(t), \quad (7) $$

where

* $\mu(t)$ is the ensemble mean climate change response,

* $\alpha_i(t)$ is the main effect of emission scenario i [i.e., the mean deviation of RCP scenario i from $\mu(t)$],

* $\beta_j(t)$ is the main effect of GCM j,

* $\gamma_k(t)$ is the main effect of RCM k, and

* $\xi_{i,j,k}(t) = \phi_{i,j,k}^*(t) - \mu(t) - \alpha_i(t) - \beta_j(t) - \gamma_k(t)$ corresponds to residual terms. For each year t, $\xi_{i,j,k}(t)$ are assumed to be independent and identically distributed (i.i.d.) over all scenarios, GCMs, and RCMs, and to follow normal distributions, with mean 0 and variance $\sigma^2(t)$.

Section 2d describes how the different parameters $\mu(t)$, $\alpha_i(t)$, $\beta_j(t)$, $\gamma_k(t)$, and $\sigma^2(t)$ are estimated. To ensure identifiability of the parameters, the vectors of RCP, GCM, and RCM effects are subject to sum-to-zero constraints for each time t:

$$ \sum_{i=1}^{I} \alpha_{i}(t) = 0, \sum_{j=1}^{J} \beta_{j}(t) = 0, \text{ and } \sum_{k=1}^{K} \gamma_{k}(t) = 0. \quad (8) $$

### d. Bayesian inference

Let $\{i, j, k\}^o$ denote the ensemble of $n_o$ available runs for year t, and $\phi^{*o}(t)$ denote the corresponding vector of climate change responses. Similarly, let $\{i, j, k\}^m$ denote the ensemble of $n_m$ missing climate scenarios, which are required in order to complete the set of all scenario/ GCM/RCM possible combinations, and $\phi^{*m}(t)$, the

vector of these $n_m$ missing climate change responses. For each year, the unknown quantities to infer include the different effects of Eq. (7), $\mu, \alpha = \{\alpha_i\}, i = 1, \dots, I, \beta = \{\beta_j\}, j = 1, \dots, J, \gamma = \{\gamma_k\}, k = 1, \dots, K$, the variance of residual terms $\sigma^2$, and the climate change response of missing scenario/GCM/RCM combinations $\phi^{*m}$, where the year index t is removed in order to lighten the equations.

In this study, statistical inference of all missing quantities $\theta = \{\phi^{*m}, \mu, \alpha, \beta, \gamma, \sigma^2\}$ is performed using Bayesian methods [see, e.g., monographs by Robert (1994) and Gelman et al. (2013)]. A Bayesian inference combines available information (here the available runs $\phi^{*o}$), represented by the likelihood function, with prior knowledge about the parameters of the likelihood. An advantage of Bayesian methods is that missing quantities that are not part of the likelihood can also be introduced if they can ease the inference. These so-called latent variables are treated as additional parameters and are also part of the inference procedure. When these latent variables can be interpreted as missing data (here the missing runs $\phi^{*m}$), this procedure is called “data augmentation” (Tanner and Wong 1987). With Bayesian methods, the inference consists in obtaining the joint posterior distribution of all missing quantities which is computed using Bayes’ theorem:

$$ \begin{align} \Pi(\theta | \phi^{*o}) &\propto \Pi(\phi^{*o} | \theta) \times \Pi(\theta) \tag{9} \\ &= \underbrace{\Pi(\phi^{*o} | \theta)}_{\text{likelihood}} \times \underbrace{\Pi(\phi^{*m} | \mu, \alpha, \beta, \gamma, \sigma^2, \phi^{*o})}_{\text{Missing data given param. and obs. data}} \nonumber \\ &\quad \times \underbrace{\Pi(\mu, \alpha, \beta, \gamma, \sigma^2)}_{\text{joint prior}} \tag{10} \end{align} $$

where the symbol $\propto$ in expression (9) indicates that the normalizing constant of the posterior is not known. Concerning expression (10):

* The first term is the likelihood of available climate change responses $\phi^{*o}$ given the unknowns. As we assume that residual errors $\xi_{i,j,k}(t)$ are i.i.d. normal in Eq. (7), $\phi^{*o}$ does not depend on the missing values $\phi^{*m}$. Its elements are independent and normally distributed with mean $\mu + \alpha_i + \beta_j + \gamma_k$ and variance $\sigma^2$, such that

$$ \begin{align} \Pi(\phi^{*o} | \theta) &= \Pi(\phi^{*o} | \phi^{*m}, \mu, \alpha, \beta, \gamma, \sigma^2) \nonumber \\ &= \prod_{\{i,j,k\}^o} \Pi(\phi_{i,j,k}^* | \mu, \alpha_i, \beta_j, \gamma_k, \sigma^2). \tag{11} \end{align} $$

* The second term is the distribution associated with the missing climate change responses, given the ANOVA
model parameters. As indicated in the introduction, we assume that the additive model (7) is valid for missing and nonmissing chains. Similarly to available climate change responses $\phi^{*o}$, the missing elements $\phi_{i,j,k}^{*m}$ are independent and normally distributed with mean $\mu + \alpha_i + \beta_j + \gamma_k$ and variance $\sigma^2$. This leads to

$$ \Pi(\boldsymbol{\phi}^{*m} | \mu, \alpha, \beta, \gamma, \sigma^2, \boldsymbol{\phi}^{*o}) = \Pi(\boldsymbol{\phi}^{*m} | \mu, \alpha, \beta, \gamma, \sigma^2) \\ = \prod_{\{i,j,k\}^m} \Pi(\phi_{i,j,k}^* | \mu, \alpha_i, \beta_j, \gamma_k, \sigma^2). \qquad (12) $$

* The third term is the joint prior distribution of the ANOVA model parameters, for which we consider independent priors, a classical nonrestrictive assumption:

$$ \Pi(\mu, \alpha, \beta, \gamma, \sigma^2) = \Pi(\mu) \times \Pi(\alpha) \times \Pi(\beta) \times \Pi(\gamma) \times \Pi(\sigma^2). \qquad (13) $$

The posterior distribution $\Pi(\boldsymbol{\theta}|\boldsymbol{\phi}^{*o})$ cannot be expressed or sampled directly. In this case, a classical solution is to generate a large number of draws from this distribution using a Markov chain Monte Carlo (MCMC) algorithm (Gilks et al. 1995; Robert and Casella 2004). Here, we obtain draws from the posterior distribution using a Gibbs sampling strategy (Casella and George 1992), which consists in sampling iteratively the conditional distributions of all unknown quantities given all other quantities. These so-called full conditional distributions are easily obtained in our case, which motivates the use of the Gibbs algorithm, and are specified in the appendix, along with the definition of the prior distributions and the choice of the hyperparameters. The prescribed hyperparameters lead to vague conjugate priors for the parameters. With regards to more generic approaches such as the Metropolis-Hastings algorithm (Metropolis et al. 1953), Gibbs sampling ensures an optimal exploration of the posterior distribution, speeding up convergence and, hence, limiting the number of required iterations.

In this study, the posterior distributions of all unknown quantities (parameters and missing climate projections) are sampled sequentially using the Gibbs algorithm. After a burn-in period of 2000 samples, 50 000 draws are retained. Several tests have been proposed to test the convergence of multiple chains (Gelman and Rubin 1992; Brooks and Gelman 1998). In our case, a visual assessment of the chains (not presented here) shows that convergence is reached very quickly, after 100 iterations.

### e. Total uncertainty and uncertainty components

Following expression (4), if $\phi_{i,j,k}^*(t)$ and $\eta_{i,j,k}^*(t)$ are assumed to be independent, the total variance of the change variable $Y_{i,j,k}^*(t)$ is

$$ \mathrm{Var}[Y_{i,j,k}^*(t)] = \mathrm{Var}[\phi_{i,j,k}^*(t)] + \mathrm{Var}[\eta_{i,j,k}^*(t)], \quad (14) $$

where $\mathrm{Var}[\phi_{i,j,k}^*(t)]$ is the total uncertainty in the climate change response and $\mathrm{Var}[\eta_{i,j,k}^*(t)]$ is the uncertainty associated to internal variability of the change variable. For each time $t$, the variance of the climate change response $\phi_{i,j,k}^*(t)$ is the sum of the variance of the different uncertainty components in Eq. (7):

$$ \mathrm{Var}[\phi_{i,j,k}^*(t)] = \mathrm{Var}[\alpha_i(t)] + \mathrm{Var}[\beta_j(t)] \\ + \mathrm{Var}[\gamma_k(t)] + \mathrm{Var}[\xi_{i,j,k}(t)], \quad (15) $$

where $\mathrm{Var}[\alpha_i(t)]$, $\mathrm{Var}[\beta_j(t)]$, and $\mathrm{Var}[\gamma_k(t)]$ are related to the dispersion between the RCP, GCM, and RCM effects, respectively, and $\mathrm{Var}[\xi_{i,j,k}(t)]$ is the residual variability. From the retained MCMC sample generated with our Gibbs algorithm, we can obtain draws of the first three components as

$$ \mathrm{Var}^{(d)}[\alpha_i(t)] = \frac{1}{I} \sum_i [\alpha_i^{(d)}(t)]^2, \qquad (16) $$

$$ \mathrm{Var}^{(d)}[\beta_j(t)] = \frac{1}{J} \sum_j [\beta_j^{(d)}(t)]^2, \qquad (17) $$

$$ \mathrm{Var}^{(d)}[\gamma_k(t)] = \frac{1}{K} \sum_k [\gamma_k^{(d)}(t)]^2, \qquad (18) $$

where $\alpha_i^{(d)}(t)$, $\beta_j^{(d)}(t)$, and $\gamma_k^{(d)}(t)$ are the $d$th draws from the posterior distribution of $\alpha_i(t)$, $\beta_j(t)$, and $\gamma_k(t)$, respectively. Variance estimates $\mathrm{Var}[\alpha_i(t)]$, $\mathrm{Var}[\beta_j(t)]$, and $\mathrm{Var}[\gamma_k(t)]$ are then obtained as the mean of these draws. The last component, $\mathrm{Var}[\xi_{i,j,k}(t)]$, is the residual variability and corresponds directly to the quantity $\sigma^2(t)$, which is part of the Bayesian inference. The mean of the draws from the posterior distribution of $\sigma^2(t)$ is thus a natural point estimate of this variance.

Internal variability of each simulation chain corresponds to the variance of $\eta_{i,j,k}^*(t)$, that is, the deviations of $Y_{i,j,k}^*(t)$ from the chain's climate change response [see Eq. (4)]. Here, we assume that internal variability is constant over the available period (see Fig. 3 for a visual assessment of this assumption). The internal variability component of the MME is then estimated as the multichain mean of the variance over time, for each chain, of $\eta_{i,j,k}^*(t)$:

$$ \widehat{\mathrm{Var}}(\eta_{i,j,k}^*) = \frac{1}{n_o} \sum_{\{i,j,k\}^o} \frac{1}{n_{i,j,k}} \sum_t \eta_{i,j,k}^*(t)^2, \qquad (19) $$
<figure><img src="image_3.png" /><figcaption>FIG. 1. Map of France with the locations of the four mountain massifs.</figcaption></figure>

where $n_{i,j,k}$ corresponds to the length of the time series available for the chain {$i, j, k$}.

### 3. Illustration of the Bayesian ANOVA

#### a. Data

The QUALYPSO method is applied to transient climate simulations of annual precipitation and annual mean temperature available for four French mountain massifs located in the Pyrenees (Pays-Basque and Cerdagne-Canigou) and in the Alps (Mont-Blanc and Haut-Var Haut-Verdon) (see Fig. 1). Total precipitation (rainfall and snow) refers to the mass of water per unit

area (kg m⁻²), which can also be expressed as a water flux [i.e. volume per area and per unit of time (mm)]. These four massifs are a subset of the geographical delineation of French mountain regions for climatological purposes (Durand et al. 2009) and were chosen to maximize the contrast in meteorological and climatological conditions. Mont-Blanc is typical of the internal mountain ranges of the northern French Alps, near the Swiss and Italian borders. Haut-Var Haut-Verdon is typical of the southern French Alps. Pays-Basque is the westernmost massif of the French Pyrenees, at a short distance from the Atlantic Ocean, while Cerdagne-Canigou is the easternmost massif, close to the Mediterranean Sea. These climate projections include 26 simulations obtained with *I* = 2 RCP emission scenarios (RCP4.5 and RCP8.5; see van Vuuren et al. 2011) and 13 different combinations of *J* = 5 GCMs and *K* = 6 RCMs (see Table 1). Simulation chains are composed of historical runs for the periods 1950–2005 (for 12 chains), 1970–2005 (for 8 chains), or 1981–2005 (for 6 chains), and of future runs for the period 2006–2100 (2006–99 for the 6 chains starting in 1981). They have been produced from EUR-11-CORDEX projections (Jacob et al. 2014; Kotlarski et al. 2014) for different elevation bands of the French massifs, using the Adaptation of RCM Outputs to Mountain (ADAMONT) statistical adjustment method (Verfaillie et al. 2017, 2018) and the Système d’analyse fournissant des renseignements atmosphériques à la neige (SAFRAN) meteorological reanalysis as an observation dataset (Durand et al. 2009). In this paper, we consider projections obtained at 1500 m, a typical elevation representative of the mountain environment, allowing comparison between all the four massifs used here, given that this elevation level is common to the four massifs chosen.

#### b. Climate responses

The climate response function of each simulation chain is estimated with cubic smoothing splines (de Boor 1978), using the entire time series (i.e., starting in 1950,

<table><caption>TABLE 1. Combination of available GCM and RCM climate projections with scenarios RCP4.5 and RCP8.5. More information about the climate models can be found in Verfaillie et al. (2018). (Expansions of acronyms are available online at http://www.ametsoc.org/PubsAcronymList.)</caption><thead><tr><th rowspan="2">RCM</th><th colspan="5">GCM</th></tr><tr><th>CNRM-CM5</th><th>EC-EARTH</th><th>IPSL-CM5A-MR</th><th>HadGEM2-ES</th><th>MPI-ESM-LR</th></tr></thead><tbody><tr><td>CCLM 4.8.17</td><td>X</td><td>X</td><td></td><td>X</td><td>X</td></tr><tr><td>ALADIN 53</td><td>X</td><td></td><td></td><td></td><td></td></tr><tr><td>WRF 3.3.1.F</td><td></td><td></td><td>X</td><td></td><td></td></tr><tr><td>RACMO 2.2E</td><td></td><td></td><td></td><td>X</td><td></td></tr><tr><td>REMO 2009</td><td></td><td></td><td></td><td></td><td>X</td></tr><tr><td>RCA 4</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>FIG. 2. Climate projections $Y_{i,j,k}(t)$ (dotted lines) and climate responses $\phi_{i,j,k}(t)$ (plain lines) for (left) mean annual temperature (°C) and (right) annual precipitation (mm or kg m⁻²) for the Haut-Var Haut-Verdon massif, for scenarios (top) RCP4.5 and (bottom) RCP8.5, and for each available combination of climate models GCM/RCM. The different colors correspond to the different GCMs (the different lines of a given color correspond to the different RCMs driven by the same GCM). Vertical dotted lines indicate the start of future runs in 2006.</figcaption></figure>

1970, and 1981 and ending in 2100). The fit is obtained with the function smooth.spline in R software (https://www.r-project.org/). Further details can be found in Hastie and Tibshirani (1990, chapter 3.4). To obtain smooth trends, a large smoothing parameter is chosen ($\lambda = 1$). The corresponding number of degrees of freedom (knots) is 3.7 (79), 3.5 (72), and 3.35 (68) for runs starting in 1950, 1970, and 1981 (and ending in 2100), respectively. Figure 2 shows climate responses $\phi_{i,j,k}(t)$ obtained for mean annual temperature (°C) and annual precipitation (mm or kg m⁻²) for the Haut-Var Haut-Verdon massif. Mean annual temperature is obviously expected to increase for all scenario/GCM/RCM configurations. For annual precipitation, contrasted climate responses are simulated. Similar conclusions hold for the other massifs (results not shown).

Figure 3 shows the deviations $\eta_{i,j,k}^*(t)$ from the climate change response, as a result of internal variability, for the Haut-Var Haut-Verdon massif. Here, internal variability corresponds to the combined effect of low-frequency internal variability and higher-frequency interannual variability, because annual scale data are used. These signals are roughly homoscedastic (e.g., no trend in the mean deviation, roughly constant variability of the deviations in time), which supports the

assumption of a constant internal variability over the simulation period. Similar results are obtained for the other massifs (not shown).

### c. Mean climate change response and posterior of the main effects

For each scenario/GCM/RCM combination and each year $t$, we obtain the climate change response in terms of absolute changes for temperature [Eq. (5)] and relative changes for precipitation [Eq. (6)], compared to the control year $c = 1990$. For each year $t = 1990, \dots, 2011$, we then obtain the posterior distributions of the mean climate change response function for each RCP scenario [$\mu(t) + \alpha_i(t)$], as well as the main effects of the different GCMs and RCMs. These posterior distributions are presented in Fig. 4 for the Haut-Var Haut-Verdon massif. Roughly similar results are obtained for the other massifs (see Figs. S1–S3 in the online supplemental material). This figure also presents the 95% credible intervals representing the uncertainty in the estimates, obtained from the distribution of the 50 000 draws. As a result of the smooth climate change response $\phi_{i,j,k}^*(t)$ extracted independently for each simulation chain, the mean climate change response function of each RCP and the main effect of each climate model evolve slowly and smoothly
<figure><img src="image_5.png" /><figcaption>FIG. 3. Deviations from the climate change response $\eta_{i,j,k}^*(t)$ as a result of internal variability for (left) mean annual temperature (absolute changes; °C) and (right) annual precipitation (relative changes; unitless) for the Haut-Var Haut-Verdon massif, for each available combination of climate models and GCMs and RCMs, and for scenarios (top) RCP4.5 and (bottom) RCP8.5. The different colors correspond to the different GCMs (the different lines of a given color correspond to the different RCMs driven by the same GCM).</figcaption></figure>

with time. For the considered MME, temperature is expected to increase up to 3°C (6°C) at the end of the century for scenario RCP4.5 (RCP8.5). Precipitation is expected to decrease up to 10% between 1990 and 2100 for scenario RCP8.5. For scenario RCP4.5 the mean expected precipitation change is less than a few percent regardless the projection lead time.

For both variables, the climate model effects (GCMs or RCMs) logically start from zero for the control period and tend to increase with the projection lead time (GCMs effects actually tend to stabilize during the second part of the century). For the Haut-Var Haut-Verdon massif, for both variables, the effects due to RCMs are of the same order (of magnitude) as those due to GCMs.

For temperature, the main effects of the different GCMs and/or RCMs in 2100 are less than ±1°C. The most important GCM effects are obtained for HadGEM2-ES and for CNRM-CM5, which produce slightly larger and smaller warming than the other GCMs, respectively. Similarly, RCA4 warms up to 0.8°C more than the other RCMs and WRF 3.3.1.F up to 0.8°C less. Whatever the simulation chain, these effects are much less than the mean expected warming.

For precipitation, the main effects of GCMs and/or RCMs are up to ±10% in 2100. The most important

effects are obtained for GCMs CNRM-CM5 (+10%) and IPSL-CM5A-MR (−10%) and for RCM WRF 3.3.1.F (+10%). These effects are roughly the same as the mean expected change. Similar results are obtained for the other massifs (see Figs. S1–S3). An exception is the Mont-Blanc massif, the alpine region with the highest elevations and altitude gradients, for which the dispersion between the RCM effects is larger than that the dispersion between the GCM effects.

Figure 5 shows the standard deviation $\sigma(t)$ of the residual terms $\xi_{i,j,k}(t)$, as a function of time, for absolute temperature changes (°C), and relative precipitation changes (unitless), for the Haut-Var Haut-Verdon massif. As indicated above, the standard deviation is related to the magnitude of the errors $\xi_{i,j,k}(t)$, due to the imperfect representation of the climate change response as a sum of main effects. For both temperature and precipitation, $\sigma(t)$ increases as a function of time. For precipitation, as a result of the stabilization of the main effects around midcentury, $\sigma(t)$ is almost constant during the period 2030–70.

### d. Total uncertainty and uncertainty components

As mentioned previously, different types of uncertainty contribute to the total uncertainty of the climate
<figure><img src="image_6.png" /></figure>

FIG. 4. Decomposition of the effects contributing to the variance of the climate projections for absolute temperature changes (°C), and relative precipitation changes (unitless), compared to year 1990, for the Haut-Var Haut-Verdon massif. These results are available for the other massifs in Figs. S1–S3. Posterior distribution of (a) the mean climate change response function for each RCP scenario $\mu(t) + \alpha_i(t)$, (b) the GCM effects $\beta_i(t)$, and (c) RCM effects $\gamma_k(t)$, as a function of time. The uncertainty on the estimation of each individual effect is represented by the colored bounds around the expected value (95% credible intervals).

projections, namely uncertainties related to RCP scenarios, GCMs, RCMs, residual variability, and internal variability. For each RCP scenario, the evolution of the corresponding total variance is presented in Figs. 6 and 7. We also present the fraction of total variance explained by each source of uncertainty.

For temperature, the main source of uncertainty is internal variability (in orange in Figs. 6 and 7) until midcentury, which is then exceeded by the uncertainty due to the RCP scenario. For all massifs, differences between the results obtained with RCP4.5 and RCP 8.5 are clear at the end of the century and scenario uncertainty grows up to 60% of the total uncertainty in 2100. Except for the Mont-Blanc massif, the part of total variance due to RCM and GCM uncertainty is no more

than 30% and decreases down to less than 20% at the end of the century.

For precipitation, total uncertainty is large regardless of the projection lead time and mainly due to internal variability. Climate model uncertainty (GCMs and RCMs) increases with lead time but always remains less than 30% of total uncertainty (except for the Mont-Blanc massif). Uncertainty related to RCP scenarios remains small and even negligible compared to climate models uncertainty, especially for the Mont-Blanc massif. For this massif, the large fraction related to the RCM uncertainty is mainly explained by the strong departure of the RCM WRF 3.3.1.F, which is much “colder” and “wetter” than the other RCMs (see Figs. S1–S3). The large dispersion between RCM effects
<figure><img src="image_7.png" /><figcaption>FIG. 5. Standard deviation $\sigma(t)$ of the residual terms $\xi_{i,j,k}(t)$ as a function of time for absolute temperature changes (°C) and relative precipitation changes (unitless) for the Haut-Var Haut-Verdon massif. The posterior mean (plain line) and 95% credible intervals (gray interval) are shown.</figcaption></figure>

for both temperature and precipitation highlights the contribution of the regional model, which is determinant in mountainous regions.

Note that for both variables, the fraction of total variance due to the residual variability (in yellow in Figs. 6 and 7) is small. This confirms that first-order terms in the ANOVA model (7) explain the main part of the variability of the climate change responses, so that further interaction terms are not necessarily needed in the ANOVA model.

### e. Significance of changes

The colored intervals in Figs. 6a, 6b, 7a, and 7b represent the total uncertainty, and correspond to the confidence interval of possible future changes at the 90% confidence level. A significant climate change realization is expected when zero (i.e., no change) lies outside this confidence interval.

The time of emergence of a significant warming (Giorgi and Bi 2009), defined here as the first future lead time for which zero lies outside the confidence interval, is found within the 2030–40 period for both emission scenarios, with the exception of the Mont-Blanc massif with the scenario RCP4.5. For precipitation, no significant change is observed, even at the end of the century. Precipitation that may be experienced for a given future period could therefore be higher or lower than what has been observed for the control period. This is mainly due to the large internal variability for this variable.

## 4. Discussion

### a. Internal variability

In our application, the contribution of internal variability to total variability of the ensemble is large, for

both variables and all regions. Indeed, in our case, internal variability includes a high-frequency (interannual) variability at the annual scale and a variability corresponding to lower frequencies (e.g., decadal scales), which could possibly be separated (Solomon et al. 2011; Seitola 2016). However, as shown by Vidal et al. (2016), for instance, a smaller contribution would have been obtained by applying QUALYPSO to time-slice averages instead of annual values (e.g., 20-yr mean precipitation instead of annual precipitation), as done in many studies (see, e.g., Yip et al. 2011; Lafaysse et al. 2014). Furthermore, it is important to note that spatial aggregation also removes an important part of the variability, which explains that regional projections lead to a larger internal variability than global means (Kendon et al. 2008; Hawkins and Sutton 2009).

### b. Uncertainty of the estimation

The main output of Bayesian methods is the posterior distributions obtained for the different unknown quantities to be estimated. The distribution obtained for each quantity provides a direct assessment of the uncertainty associated to the estimation. We focus on uncertainty related to the estimation of ANOVA model parameters (i.e., mean climate change response, scenarios, and climate model effects). Figure 8 represents the standard deviation of the posterior distribution obtained for the mean climate change response of each RCP scenario, for the GCM and RCM effects. These standard deviations are related to the width of the credible intervals provided in Fig. 4 (the width of the 95% is about 4 times the standard deviation when posterior distributions are Gaussian). For all effects and for both temperature and precipitation, the estimation uncertainty roughly increases with time and directly follows, by construction,
<figure><img src="image_8.png" /><figcaption>FIG. 6. Total uncertainty and fraction of total variance explained by each source of uncertainty for mean annual temperature changes (°C) compared to year 1990 as a function of time. Mean climate change response (white curve) for scenarios (a) RCP4.5 and (b) RCP8.5 [$\mu(t) + \alpha_i(t)$] and corresponding 90% confidence interval. The width of the overall colored interval representing the different sources of uncertainty corresponds to $\mu(t) + \alpha_i(t) \pm 1.645\sqrt{\text{Var}[Y_{i,t,k}^*(t)]}$. This interval is subdivided into subintervals whose width is proportional to the fraction of total uncertainty variance explained by the different sources of uncertainties. (c) Fraction of total variance explained by each source of uncertainty.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 7. As in Fig. 6, but for relative changes of annual precipitation (unitless).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 8. Uncertainty of the estimation (standard deviation of the posterior distribution) for (a) the mean climate change response for each RCP scenario $\mu(t) + \alpha_i(t)$, (b) the GCM effects $\beta_j(t)$, and (c) the RCM effects $\gamma_k(t)$ as a function of time for absolute temperature changes (°C) and relative precipitation changes (unitless) for the Haut-Var Haut-Verdon massif. Note that different line widths are used to show overlapping curves and do not represent different uncertainties.</figcaption></figure>

the evolution of the standard deviation $\sigma(t)$ obtained for the residual errors of the ANOVA model [see Fig. 5 and Eqs. (A2), (A9), (A11), and (A13) in the appendix].

The magnitude of the estimation uncertainty depends mostly on the size and setup of the dataset available for the estimation. For climate model effects for instance, it depends on the number of available RCM/GCM combinations. As an illustration, the estimation uncertainty of the GCM effect is similar for all GCMs, except for IPSL-CM5A-MR for which it is significantly larger. Conversely to the other GCMs, IPSL-CM5A-MR was not used to drive the RCM CCLM 4.8.17, for which the effect is well estimated since it is available for four GCMs. Similarly, the same estimation uncertainty is obtained for the RCMs ALADIN 53, RACMO 2.2E,

and REMO 2009. These three RCMs share the same configuration. They are driven by only one GCM, which has also been used to run the RCMs CCLM 4.8.17 and RCA 4. The estimation uncertainty is lower for RCMs CCLM 4.8.17 and RCA 4 as they are driven by 4 and 5 GCMs, respectively. The largest estimation uncertainty obtained for the RCM WRF 3.3.1.F is finally due to the fact that it was driven by the GCM IPSL-CM5A-MR only, which has the largest estimation uncertainty among all GCM effects.

As illustrated above, missing scenario/RCM/GCM combinations logically determine the estimation uncertainty of the different parameters of the analysis (climate model effects, mean climate change response, etc.). The data augmentation approach accounts for this
<figure><img src="image_11.png" /><figcaption>FIG. 9. Mean climate change response for each RCP scenario $\mu_i(t) = \mu(t) + \alpha_i(t)$ as a function of time for absolute temperature changes (°C) and relative precipitation changes (unitless) for the Haut-Var Haut-Verdon massif. Bayesian estimates $\mu_i^{\text{bay}}(t)$ and direct empirical estimates $\mu_i^{\text{emp}}(t)$ of $\mu_i(t)$ are shown. For the Bayesian estimates, 95% credible intervals are shown.</figcaption></figure>

and propagates the uncertainty due to missing data in the analysis.

### c. Comparison of Bayesian and direct estimates

In numerous climate impact studies, the mean climate change response is estimated by the empirical mean of all available climate experiments. When many scenario/ RCM/GCM combinations are missing, this empirical estimate might be significantly different from the mean response that would be obtained with the complete ensemble.

Let $\mu_i(t) = \mu(t) + \alpha_i(t)$ denote the mean climate change response for each RCP scenario. Figure 9 compares the Bayesian estimates $\mu_i^{\text{bay}}(t)$ of $\mu_i(t)$ (already shown in Fig. 4a) to the empirical estimates $\mu_i^{\text{emp}}(t) = (1/n_{i,o}) \sum_{\{j,k\}^{i,o}} \phi_{i,j,k}^*(t)$, where $\{j,k\}^{i,o}$ denotes the ensemble of $n_{i,o}$ available runs for year $t$ and scenario $i$. At the end of the century, regardless of the scenario, direct empirical estimates lead to a difference on the order of 0.1°C for temperature and 1% for precipitation.

When the mean response is estimated by the empirical mean of available scenario/GCM/RCM combinations, each available climate experiment has the same weight in the estimation. However, some climate models are more represented than other ones. For example, RCMs CCLM 4.8.17 and RCA 4 concern 8 and 10 climate experiments, respectively, among the $n_o = 26$ climate experiments. For temperature changes, this logically leads to a slight overestimation of the mean climate change response (the effect of CCLM 4.8.17 is roughly zero for the whole century; RCA 4 presents the largest positive effect, reaching up to +0.8°C at the end of the century; see Fig. 4).

To better illustrate the impact of overrepresented climate models on the mean estimated trends, if direct estimates are used, let us further consider a synthetic

MME composed of 5 GCMs and 5 RCMs. For each GCM $j$ and each RCM $k$, climate projections $\tilde{Y}_{j,k}(t)$ are generated according to the following formula:

$$ \tilde{\phi}_{j,k}(t) = \tilde{\mu}(t) + \tilde{\beta}_j(t) + \tilde{\gamma}_k(t) + \tilde{\eta}(t), \quad (20) $$

where the mean climate change response is constant and equal to 0 [$\tilde{\mu}(t) = 0$] (i.e., there is no trend) and GCM effects are proportional to time $t$ and are equal to $\tilde{\beta}_j(t) = a_j \times t$, with $t = 0.01, 0.02, \dots, 1$ and $a_1 = 1, a_2 = 0.5, a_3 = 0, a_4 = -0.5$, and $a_5 = -1$. RCM effects are also proportional to time $t$ and are equal to $\tilde{\gamma}_k(t) = b_k \times t$, with $t = 0.01, 0.02, \dots, 1$ and $b_1 = 0.5, b_2 = 0.25, b_3 = 0, b_4 = -0.25$, and $b_5 = -0.5$, and $\tilde{\eta}(t)$ is a random noise mimicking internal variability and residual variability, which follows a Gaussian distribution with mean 0 and a standard deviation of 0.3.

This configuration corresponds more or less to the MMEs presented previously for the temperature projections (see Fig. 4). In the present case, however, we consider a different set of available experiments, as indicated in Table 2. Clearly, this set is unbalanced in the sense that GCM1 and RCM1 are overrepresented. This situation is quite similar to real-world applications for which one or two GCMs/RCMs are more used than other existing climate models. Figure 10a shows the generated climate projections $\tilde{\phi}_{j,k}(t)$ and the corresponding smooth signals $\phi_{j,k}(t)$ extracted using cubic splines.

Bayesian estimates of the mean climate change response $\mu_{\text{bay}}(t)$, obtained with QUALYPSO and direct empirical estimates $\mu_{\text{emp}}(t) = (1/n_o) \sum_{\{j,k\}} \phi_{j,k}^*(t)$, obtained for each time $t$ as the mean of estimated climate change responses, are presented in Fig. 10b. Bayesian estimates are logically unbiased in contrast to empirical ones, which significantly overestimate, in the present
<table><caption>TABLE 2. Combinations of available GCM and RCM climate projections for the synthetic experiment.</caption><thead><tr><th rowspan="2">RCM</th><th colspan="5">GCM</th></tr><tr><th>GCM1</th><th>GCM2</th><th>GCM3</th><th>GCM4</th><th>GCM5</th></tr></thead><tbody><tr><td>RCM1</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td></tr><tr><td>RCM2</td><td>X</td><td>X</td><td></td><td></td><td></td></tr><tr><td>RCM3</td><td>X</td><td></td><td>X</td><td></td><td></td></tr><tr><td>RCM4</td><td>X</td><td></td><td></td><td>X</td><td></td></tr><tr><td>RCM5</td><td>X</td><td></td><td></td><td></td><td>X</td></tr></tbody></table>

configuration, the true mean climate change response $\tilde{\mu}(t) = 0$ (solid line). In contrast to our real application, for which QUALYPSO and empirical approaches lead to close results due to a more balanced projection set, this synthetic example clearly shows the added value of our approach.

## 5. Conclusions and outlooks

This work presents the development and application of a Bayesian approach, named QUALYPSO, which uses data augmentation to assess the different sources of uncertainty in incomplete multiscenario multimodel ensembles (MMEs) of climate experiments. In a first step, the climate response of each available simulation chain is estimated with a trend model (i.e., cubic splines here) fitted to raw climate projections. Climate change responses can then be obtained and residuals from the climate change response are used to estimate the internal variability of the chain. The other uncertainty components of the projections (scenario uncertainty and climate model uncertainty) are estimated with a Bayesian ANOVA model applied to the climate change responses of available simulation chains. The ANOVA model provides an estimate of the mean climate change response of the MME, as well as an estimate of the main effects of the different emission scenarios and climate models (GCMs and RCMs). The different parameters of the ANOVA model and the missing quantities associated to the missing simulation chains are jointly estimated using data augmentation techniques. For illustration, we apply QUALYPSO to MMEs of climate projections (mean annual temperature and total annual precipitation) produced for four different French masses at 1500-m elevation. Projections are available for 13 GCM/RCM combinations and two emission scenarios.

QUALYPSO presents many advantages over more classical estimation approaches. Along with the estimation of missing data, it provides an assessment of the estimation uncertainty and adequately propagates the uncertainty due to missing scenario/GCM/RCM combinations. With the explicit treatment of missing climate

experiments, it is then expected to produce unbiased estimates of all parameters, in contrast to direct empirical estimates, typically obtained as the average over all available projections. The Bayesian approach also exploits all available climate experiments, avoiding a dramatic loss of information when standard single-time or time series approaches are applied (in this case, a classical solution is to select a complete subset of climate experiments). The QUALYPSO methodology can be applied to any kind of climate variable and any kind of MMEs of climate projections. Additional effects due to some impact model (e.g., a hydrological model or a snow model; see Vidal et al. 2016; Verfaillie et al. 2018) or to some other explanatory factor (e.g., spatial effects; see Geinitz et al. 2015; Tingley 2012) could also be easily included in the analysis.

A key step of QUALYPSO is the estimation of the climate response for each simulation chain. In the present analysis, the climate responses are modeled by

<figure><img src="image_12.png" /><figcaption>FIG. 10. (a) Synthetic climate change projections $\phi_{i,j}(t)$ (dotted lines) and estimated trends $\phi_{i,j}(t)$ (plain lines) for a synthetic ensemble of climate experiments. Each time series corresponds to one GCM/RCM combination. All times series with a same color correspond to GCM/RCM combinations with the same GCM. (b) Mean response: true mean response function $\tilde{\mu}(t)$, equal to 0 (plain horizontal line), Bayesian estimates $\mu_{\text{bay}}(t)$ of $\tilde{\mu}(t)$ (gray interval), and direct empirical estimates $\mu_{\text{emp}}(t)$ (dotted line). For the Bayesian estimates, 95% credible intervals are shown.</figcaption></figure>
cubic splines. Many alternatives are possible and have been considered in previous works. For example, Hawkins and Sutton (2009) and Hingray and Saïd (2014) apply fourth-order polynomial functions to temperature changes. Cubic splines are expected to provide a flexible and robust fit of the climate response. This estimation can also be improved using multiple runs of each chain, if available (e.g., Kendon et al. 2008; Deser et al. 2012; Hingray et al. 2019). In all cases, an improvement of the method would be to consider the estimated climate response of each chain as uncertain, and to account for this uncertainty in the Bayesian ANOVA approach.

Acknowledgments. This work was supported by the French National program LEFE (Les Enveloppes Fluides et l'Environnement) and by the French National Research Agency, via the CDP-Trajectories project, in the framework of the “Investissements d'avenir” program (ANR-15-IDEX-02). This study also benefited from the outcomes of the ADAMONT project (GICC program and ONERC) and of the Interreg project POCTEFA/Clim'Py. Irstea, IGE and CNRM/CEN are part of Labex OSUG@2020. Deborah Verfaillie's work has been funded by the European project EUCP (H2020-SC5-2016-776613). We also thank Matthieu Lafaysse (Univ. Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM, CEN, Grenoble, France), Samuel Somot, Aurélien Ribes and Michel Déqué (CNRM, Météo-France, CNRS, Université de Toulouse, Toulouse, France) and Wilfried Thuiller (LECA, Grenoble, France) for interesting discussions on this subject. We also thank the editor and two anonymous reviewers for their constructive comments, which helped us to improve the manuscript.

# APPENDIX

## Full Conditional Distributions for the Decomposition of Climate Change Response

As indicated in section 2, QUALYPSO applies a Bayesian approach in order to infer the different unknown quantities (parameters of the ANOVA model and missing scenario/RCM/GCM combinations). The joint posterior distribution of all unknown quantities is sampled using the Gibbs algorithm, which necessitates to sample iteratively the full conditional posterior distributions of all unknown quantities. In this section, we specify conjugate priors for the parameters, we provide the full conditional posteriors, and we propose values for the hyperparameters, following Tingley (2012). In the remainder of this section, $A|.$ denotes the unknown

quantity A conditionally on all other variables. The Bayesian ANOVA approach is applied independently for each time step $t$. However, in the following equations, the year index $t$ is removed in order to lighten the equations.

### a. Mean climate change response function $\mu$

As indicated above, we assume that the climate change response $\phi_{i,j,k}^*$ is the sum of fixed effects $\mu$, $\alpha_i$, $\beta_j$, $\gamma_k$, and residual terms $\xi_{i,j,k} \sim N(0, \sigma^2)$. As a consequence, the normal distribution is a conjugate prior for the mean response function $\mu$:

$$ \mu \sim N(\mu_0, \sigma_{\mu}^2). \quad (A1) $$

The conditional posterior is likewise normal:

$$ \mu| \sim N(V_{\mu}, \psi_{\mu}, V_{\mu}), \quad (A2) $$

where $\psi_{\mu} = (1/\sigma^2) \sum_{\{i,j,k\}} \phi_{i,j,k}^* + (\mu_0/\sigma_{\mu}^2)$, and $V_{\mu} = \{[(I \times J \times K)/\sigma^2] + (1/\sigma_{\mu}^2)\}^{-1}$.

Note that the expression for $\psi_{\mu}$ considers all climate change responses $\phi_{i,j,k}^*$, including climate change responses that are not directly available and are estimated.

### b. Variance of the residual terms $\sigma^2$

For the variance parameter $\sigma^2$, a natural conjugate prior is the inverse-gamma distribution:

$$ \sigma^2 \sim \text{IGa}(\kappa, \nu), \quad \text{and} \quad (A3) $$

$$ p(\sigma^2) \propto (\sigma^2)^{-(\kappa+1)} \exp(-\nu/\sigma^2). \quad (A4) $$

It follows that the conditional posterior distribution remains inverse-gamma distributed:

$$ \sigma^2 \sim \text{IGa} \left[ \frac{I \times J \times K}{2} + \kappa, \frac{1}{2} \sum_{\{i,j,k\}} (\phi_{i,j,k}^* - \mu - \alpha_i - \beta_j - \gamma_k)^2 + \nu \right]. \quad (A5) $$

### c. Effect of the GCM $\beta$

We first present the treatment of GMC effects in order to illustrate the use of Helmert contrasts with more than $I = 2$ effects. The sum-to-zero constraint implies that there are only $J - 1$ free parameters in the vector $\beta$. Following Tingley (2012), a $J$ by $J - 1$ matrix $Q_J$ can be used to make the correspondence between the constrained $\beta$ vector and an unconstrained $\beta^*$ vector of $J - 1$ i.i.d. normal variables, such that the linear transform $\beta = Q_J\beta^*$ has the correct covariance matrix. A possible choice for $Q_J$ is
$$Q_j = Q_j^*(Q_j^T Q_j)^{-1/2}, \quad (A6)$$

where $Q_j^*$ is composed of columns of Helmert contrasts:

$$Q_j^* = \begin{pmatrix} 0 & 0 & 0 & -4 \\ 0 & 0 & -3 & 1 \\ 0 & -2 & 1 & 1 \\ -1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 \end{pmatrix}, \text{ for } J = 5. \quad (A7)$$

We first specify the prior and the conditional posterior for the vector $\boldsymbol{\beta}^*$ of $J-1$ parameters. The prior for $\boldsymbol{\beta}^*$ given $\sigma_{\boldsymbol{\beta}}^2$ is normal:

$$\boldsymbol{\beta}^* | \sigma_{\boldsymbol{\beta}}^2 \sim N(\mathbf{0}_{J-1}, \sigma_{\boldsymbol{\beta}}^2 \mathbf{1}_{J-1}), \quad (A8)$$

where $\mathbf{1}_{J-1}$ represent a $J-1$ by $J-1$ identity matrix and $\mathbf{0}_{J-1}$ is a $J-1$ null vector. The full conditional posterior for $\boldsymbol{\beta}^*$ is likewise normal:

$$\boldsymbol{\beta}^* | \sim N(V_{\boldsymbol{\beta}*} \boldsymbol{\psi}_{\boldsymbol{\beta}*}, V_{\boldsymbol{\beta}*} \mathbf{1}_{J-1}), \quad (A9)$$

where $\boldsymbol{\psi}_{\boldsymbol{\beta}*} = (1/\sigma^2) \mathbf{Q}_j^T \sum_{\{i,k\}} (\boldsymbol{\phi}_{i..,k}^* - \mu \mathbf{1}_J)$ and $V_{\boldsymbol{\beta}*} = \{[(I \times K)/\sigma^2] + (1/\sigma_{\boldsymbol{\beta}}^2)\}^{-1}$. The expression for $V_{\boldsymbol{\beta}*}$ indicates that the estimation uncertainty of the GCM affects decreases as a function of the number of RCP scenarios and RCMs.

### d. Effect of the RCM $\gamma$

Similar expressions are obtained for $\gamma$. The prior for $\gamma^*$ given $\sigma_{\gamma}^2$ is normal:

$$\gamma^* | \sigma_{\gamma}^2 \sim N(\mathbf{0}_{K-1}, \sigma_{\gamma}^2 \mathbf{1}_{K-1}), \quad (A10)$$

which leads to a normal conditional distribution for $\gamma^*$:

$$\gamma^* | \sim N(V_{\gamma*} \boldsymbol{\psi}_{\gamma*}, V_{\gamma*} \mathbf{1}_{K-1}), \quad (A11)$$

where $\boldsymbol{\psi}_{\gamma*} = (1/\sigma^2) \mathbf{Q}_j^T \sum_{\{i,j\}} (\boldsymbol{\phi}_{i..,j}^* - \mu \mathbf{1}_K)$ and $V_{\gamma*} = \{[(I \times J)/\sigma^2] + (1/\sigma_{\gamma}^2)\}^{-1}$.

### e. Effect of the RCP $\alpha$

We specify a normal prior for $\alpha^*$ given $\sigma_{\alpha}^2$:

$$\alpha^* | \sigma_{\alpha}^2 \sim N(\mathbf{0}_{I-1}, \sigma_{\alpha}^2 \mathbf{1}_{I-1}), \quad (A12)$$

which leads to a normal conditional distribution for $\alpha^*$:

$$\alpha^* | \sim N(V_{\alpha*} \boldsymbol{\psi}_{\alpha*}, V_{\alpha*} \mathbf{1}_{I-1}), \quad (A13)$$

where $\boldsymbol{\psi}_{\alpha*} = 1/\sigma^2 \mathbf{Q}_j^T \sum_{\{j,k\}} (\boldsymbol{\phi}_{..,j,k}^* - \mu \mathbf{1}_I)$ and $V_{\alpha*} = \{[(J \times K)/\sigma^2] + (1/\sigma_{\alpha}^2)\}^{-1}$.

### f. The missing values $\boldsymbol{\phi}^{*m}$

Following Eq. (7), the full conditional posterior for each missing climate change response $\boldsymbol{\phi}_{i,j,k}^*$ is normal:

$$\boldsymbol{\phi}_{i,j,k}^* | \sim N(\mu + \alpha_i + \beta_j + \gamma_k, \sigma^2). \quad (A14)$$

The mean of the posterior distribution is the sum of the fixed effect and its variance corresponds to the variance of the residual terms $\xi_{i,j,k}$, namely $\sigma^2$.

### g. Hyperparameters

The standard choices for the hyperparameters are the following:

* $\mu_0$: mean of all available climate change responses,

* $\sigma_{\mu}^2$: 16 times the variance of all available climate change responses,

* $\kappa = 0.5$ and $\nu$ is equal to half the variance of a direct estimate of residual terms (called “estimated residual variance” in Tingley (2012)), and

* $\sigma_{\alpha}$, $\sigma_{\beta}$, and $\sigma_{\gamma}$ are set to 16 times the variance of all available climate change responses.

## REFERENCES

* Bracegirdle, T. J., J. Turner, J. S. Hosking, and T. Phillips, 2014: Sources of uncertainty in projections of twenty-first century westerly wind changes over the Amundsen Sea, West Antarctica, in CMIP5 climate models. *Climate Dyn.*, **43**, 2093–2104, https://doi.org/10.1007/s00382-013-2032-1.

* Brooks, S. P., and A. Gelman, 1998: General methods for monitoring convergence of iterative simulations. *J. Comput. Graph. Stat.*, **7**, 434–455, https://doi.org/10.1080/10618600.1998.10474787.

* Casella, G., and E. I. George, 1992: Explaining the Gibbs sampler. *Amer. Stat.*, **46**, 167–174, https://doi.org/10.1080/00031305.1992.10475878.

* de Boor, C., 1978: *A Practical Guide to Splines*. Springer-Verlag, 392 pp.

* Déqué, M., and Coauthors, 2007: An intercomparison of regional climate simulations for Europe: Assessing uncertainties in model projections. *Climatic Change*, **81** (Suppl. 1), 53–70, https://doi.org/10.1007/s10584-006-9228-x.

* Deser, C., A. Phillips, V. Bourdette, and H. Teng, 2012: Uncertainty in climate change projections: The role of internal variability. *Climate Dyn.*, **38**, 527–546, https://doi.org/10.1007/s00382-010-0977-x.

* Durand, Y., M. Laternser, G. Giraud, P. Etchevers, B. Lesaffre, and L. Mérindol, 2009: Reanalysis of 44 yr of climate in the French Alps (1958–2002): Methodology, model validation, climatology, and trends for air temperature and precipitation. *J. Appl. Meteor. Climatol.*, **48**, 429–449, https://doi.org/10.1175/2008JAMC1808.1.

* Fatichi, S., and Coauthors, 2016: Uncertainty partition challenges the predictability of vital details of climate change. *Earth’s Future*, **4**, 240–251, https://doi.org/10.1002/2015EF000336.

* Geinitz, S., R. Furrer, and S. R. Sain, 2015: Bayesian multilevel analysis of variance for relative comparison across sources of
global climate model variability. *Int. J. Climatol.*, **35**, 433–443,
https://doi.org/10.1002/joc.3991.

Gelman, A., and D. B. Rubin, 1992: Inference from iterative simulation using multiple sequences. *Stat. Sci.*, **7**, 457–472, https://doi.org/10.1214/ss/1177011136.

—, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin, 2013: Bayesian Data Analysis. 3rd ed. CRC Press, 675 pp.

Gilks, W. R., S. Richardson, and D. Spiegelhalter, 1995: *Markov Chain Monte Carlo in Practice*. CRC Press, 512 pp.

Giorgi, F., and X. Bi, 2009: Time of emergence (TOE) of GHG-forced precipitation change hot-spots. *Geophys. Res. Lett.*, **36**, L06709, https://doi.org/10.1029/2009GL037593.

Giuntoli, I., J.-P. Vidal, C. Prudhomme, and D. M. Hannah, 2015: Future hydrological extremes: The uncertainty from multiple global climate and global hydrological models. *Earth Syst. Dyn.*, **6**, 267–285, https://doi.org/10.5194/esd-6-267-2015.

Hastie, T. J., and R. J. Tibshirani, 1990: *Generalized Additive Models*. CRC Press, 352 pp.

Hawkins, E., and R. Sutton, 2009: The potential to narrow uncertainty in regional climate predictions. *Bull. Amer. Meteor. Soc.*, **90**, 1095–1107, https://doi.org/10.1175/2009BAMS2607.1.

Hingray, B., and M. Saïd, 2014: Partitioning internal variability and model uncertainty components in a multimember multimodel ensemble of climate projections. *J. Climate*, **27**, 6779–6798, https://doi.org/10.1175/JCLI-D-13-00629.1.

—, A. Mezghani, and T. A. Buishand, 2007: Development of probability distributions for regional climate change from uncertain global mean warming and an uncertain scaling relationship. *Hydrol. Earth Syst. Sci.*, **11**, 1097–1114, https://doi.org/10.5194/hess-11-1097-2007.

—, J. Blanchet, G. Evin, and J.-P. Vidal, 2019: Uncertainty component estimates in transient climate projections. Precision of estimators in the single time and time series approaches. *Climate Dyn.*, https://doi.org/10.1007/s00382-019-04635-1.

Jacob, D., and Coauthors, 2014: EURO-CORDEX: New high-resolution climate change projections for European impact research. *Reg. Environ. Change*, **14**, 563–578, https://doi.org/10.1007/s10113-013-0499-2.

Kendon, E. J., D. P. Rowell, R. G. Jones, and E. Buonomo, 2008: Robustness of future changes in local precipitation extremes. *J. Climate*, **21**, 4280–4297, https://doi.org/10.1175/2008JCLI2082.1.

Kotlarski, S., and et al, 2014: Regional climate modeling on European scales: A joint standard evaluation of the EURO-CORDEX RCM ensemble. *Geosci. Model Dev.*, **7**, 1297–1333, https://doi.org/10.5194/gmd-7-1297-2014.

Lafaysse, M., B. Hingray, A. Mezghani, J. Gailhard, and L. Terray, 2014: Internal variability and model uncertainty components in future hydrometeorological projections: The Alpine Durance basin. *Water Resour. Res.*, **50**, 3317–3341, https://doi.org/10.1002/2013WR014897.

Little, R. J. A., and D. B. Rubin, 2014: *Statistical Analysis with Missing Data*. 2nd ed. Wiley, 408 pp.

Metropolis, N., A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, 1953: Equation of state calculations by fast computing machines. *J. Chem. Phys.*, **21**, 1087–1092, https://doi.org/10.1063/1.1699114.

Northrop, P. J., and R. E. Chandler, 2014: Quantifying sources of uncertainty in projections of future climate. *J. Climate*, **27**, 8793–8808, https://doi.org/10.1175/JCLI-D-14-00265.1.

Paeth, H., G. Vogt, A. Paxian, E. Hertig, S. Seubert, and J. Jacobeit, 2017: Quantifying the evidence of climate change in the light of uncertainty exemplified by the Mediterranean hot spot region. *Global Planet. Change*, **151**, 144–151, https://doi.org/10.1016/j.gloplacha.2016.03.003.

Reintges, A., T. Martin, M. Latif, and N. S. Keenlyside, 2017: Uncertainty in twenty-first century projections of the Atlantic meridional overturning circulation in CMIP3 and CMIP5 models. *Climate Dyn.*, **49**, 1495–1511, https://doi.org/10.1007/s00382-016-3180-x.

Robert, C., 1994: *The Bayesian Choice: A Decision-Theoretic Motivation*. Springer-Verlag, 436 pp.

—, and G. Casella, 2004: *Monte Carlo Statistical Methods*. 2nd ed. Springer-Verlag, 649 pp.

Seitola, T., 2016: Decomposition of the 20th century climate variability. Ph.D. thesis, Finnish Meteorological Institute, 45 pp.

Solomon, A., and Coauthors, 2011: Distinguishing the roles of natural and anthropogenically forced decadal climate variability: Implications for prediction. *Bull. Amer. Meteor. Soc.*, **92**, 141–156, https://doi.org/10.1175/2010BAMS2962.1.

Stocker, T., and Coauthors, 2013: Technical summary. *Climate Change 2013: The Physical Science Basis*. T. F. Stocker et al., Eds., Cambridge University Press, 33–115.

Tanner, M. A., and W. H. Wong, 1987: The calculation of posterior distributions by data augmentation. *J. Amer. Stat. Assoc.*, **82**, 528–540, https://doi.org/10.1080/01621459.1987.10478458.

Taylor, K. E., R. J. Stouffer, and G. A. Meehl, 2012: An overview of CMIP5 and the experiment design. *Bull. Amer. Meteor. Soc.*, **93**, 485–498, https://doi.org/10.1175/BAMS-D-11-00094.1.

Tingley, M. P., 2012: A Bayesian ANOVA scheme for calculating climate anomalies, with applications to the instrumental temperature record. *J. Climate*, **25**, 777–791, https://doi.org/10.1175/JCLI-D-11-00008.1.

van Vuuren, D. P., and Coauthors, 2011: The representative concentration pathways: An overview. *Climatic Change*, **109**, 5, https://doi.org/10.1007/s10584-011-0148-z.

Verfaillie, D., M. Déqué, S. Morin, and M. Lafaysse, 2017: The method ADAMONT v1.0 for statistical adjustment of climate projections applicable to energy balance land surface models. *Geosci. Model Dev.*, **10**, 4257–4283, https://doi.org/10.5194/gmd-10-4257-2017.

—, M. Lafaysse, M. Déqué, N. Eckert, Y. Lejeune, and S. Morin, 2018: Multi-component ensembles of future meteorological and natural snow conditions for 1500 m altitude in the Charreuse mountain range, northern French Alps. *Cryosphere*, **12**, 1249–1271, https://doi.org/10.5194/tc-12-1249-2018.

Vidal, J.-P., B. Hingray, C. Magand, E. Sauquet, and A. Ducharne, 2016: Hierarchy of climate and hydrological uncertainties in transient low-flow projections. *Hydrol. Earth Syst. Sci.*, **20**, 3651–3672, https://doi.org/10.5194/hess-20-3651-2016.

Yip, S., C. A. T. Ferro, D. B. Stephenson, and E. Hawkins, 2011: A simple, coherent framework for partitioning uncertainty in climate predictions. *J. Climate*, **24**, 4634–4643, https://doi.org/10.1175/2011JCLI4085.1.