# Partitioning uncertainty components in climate projections using smoothing splines

Guillaume Evin

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

## To cite this version:

Guillaume Evin. Partitioning uncertainty components in climate projections using smoothing splines. 2022. hal-03720621

HAL Id: hal-03720621

https://hal.inrae.fr/hal-03720621

Preprint submitted on 12 Jul 2022

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Partitioning uncertainty components in climate projections using smoothing splines

Guillaume Evin¹

¹Univ. Grenoble Alpes, INRAE, UR ETGR, Grenoble, France

## ABSTRACT

A critical issue in climate change impact studies is the assessment of uncertainties associated with future projections. Various methods have been proposed for partitioning uncertainty sources, usually based on an Analysis of Variance (ANOVA). In this paper, we show how Smoothing-Spline ANOVA approaches (SS-ANOVA) can be used to estimate the total uncertainty and its partition in climate projection ensembles. A Bayesian framework is proposed to handle heteroscedastic and autocorrelated residual errors between the climate change responses and the main additive effects modelled with cubic smoothing splines.

Keywords: Smoothing-Spline ANOVA, climate projection ensemble

## 1 INTRODUCTION

A critical issue in climate change studies is the estimation of uncertainties in projections along with the contribution of the different uncertainty sources, including scenario uncertainty, the different components of model uncertainty, and internal variability (see, e.g., Hawkins and Sutton, 2009). Scenario uncertainty is related to the possible evolution of greenhouse gas emissions, which is usually accounted for using Representative Concentration Pathways (RCP) scenarios (van Vuuren et al., 2011). Model uncertainty corresponds to the dispersion between the different climate responses obtained with different models for the same forcing configuration. Model uncertainty can concern Global Climate Models (GCMs), regional downscaling models such as regional climate models (RCMs) and/or statistical downscaling methods, and impact models (e.g. agricultural or hydrological models). Internal variability is due to the chaotic variability of the climate (Deser et al., 2012). Estimating and partitioning uncertainties in future climate projections is first intended to help evaluating the significance of estimated changes for adaptation purposes. Besides, it highlights the most important uncertainty sources for a better allocation of future research efforts.

Over the recent years, uncertainty in climate projections has been mostly explored and partitioned based on Multiscenarios Multimodel Multimember Ensembles (MMEs) of transient climate projections. Various methods have been proposed for this, most of them based on an Analysis of Variance (ANOVA) of projections available for the future time window considered (Hingray et al., 2007; Yip et al., 2011; Paeth et al., 2017). Available methods are usually applied to different future time windows in turn (Jacob et al., 2014; Northrop and Chandler, 2014; Reintges et al., 2017) which may lead to temporal fluctuations of uncertainty estimates. As reported by Hingray et al. (2019), such variations are likely to be due to a lack of robustness of the analysis. In most cases, the climate response of each simulation chain is indeed expected to evolve gradually with time. As a consequence, the different uncertainty components of a given MME should also be smooth signals. A robust estimation of all uncertainty components can thus be proposed by assuming a gradual evolution of climate responses. This is considered by Geinitz et al. (2015), who introduce an explicit model of the signals (climate responses and ANOVA effects) using simple trend models. The analysis framework described in the present study follows the same direction and extends the Bayesian approach of Evin et al. (2019). More specifically, it promotes the use of a Bayesian Smoothing-spline ANOVA (SS-ANOVA, see, e.g. Cheng and Speckman, 2012) to partition the different sources of uncertainty. The Bayesian framework handles heteroscedastic and autocorrelated residual errors between the climate change responses and the main additive effects modelled with cubic
smoothing splines.

Section 2 details the methodology proposed for partitioning climate change uncertainties. Section 4 illustrates the SS-ANOVA approach with an application to a MME of regional temperature projections obtained from the CMIP5-EUROCORDEX experiment (Jacob et al., 2014; Vautard et al., 2021). Section 5 concludes.

# 2 METHODOLOGICAL FRAMEWORK

We consider that a multi-model ensemble of climate projections $y_i$ where $i = 1, ..., n$ indicates a specific simulation chain. Each chain spans a period of $T$ years, this vector of years being denoted as $\mathbf{x}^{\text{year}} = \{x_1^{\text{year}}, ..., x_T^{\text{year}}\}$.

## 2.1 Climate response and climate change response

The proposed methodology assumes that a climate response $\phi_i$ represents the forced response of each simulation chain, as a result of climate change. Polynomial trend models or cubic splines have been used to obtain $\phi_i$ in past studies (Hawkins and Sutton, 2009; Hingray and Saïd, 2014; Evin et al., 2019). In this study, the smooth climate response $\phi_i$ is modeled as a cubic smoothing spline (Green and Silverman, 1993) capable of representing a great variety of possible evolution, including non-monotonic trends (see, e.g. Evin et al., 2019).

Most climate impact studies quantify uncertainty sources from change variables, obtained as differences between a future and a reference period. The rationale behind this transformation is that changes provided by climate projections with respect to a reference period are deemed to be more informative that raw projected values which are subject to various biases. Here, we consider a change variable defined in terms of absolute differences with respect to a reference year, i.e. $\phi_i^* = \phi_i^* - \phi_i^*$, where c indicates the reference year, which can be for example representative of the pre-industrial climate (e.g. 1850) or, as done in this study, of the current climate (c = 1999).

## 2.2 SS-ANOVA of climate change responses

Let us denote $\phi^*$ the concatenated vector of $n$ climate change responses $\phi_i^*$, $i = 1, ..., n$ of length $n \times T$. Following Cheng and Speckman (2012), we consider the following smoothing-spline analysis of variance (SS-ANOVA) decomposition:

$$ \phi^* = \sum_{e=1}^{E} \theta_e + \xi_{RV}, \qquad (1) $$

where $\theta_e$, $e = 1, ..., E$ are the main effects of the ANOVA, i.e. the evolution of the grand mean change response, shared by all simulation chains, and the effects considered as major sources of uncertainty, e.g. the evolution of the climate model (GCM/RCM) effects, main effects related to the emission scenarios. In this study, the SS-ANOVA approach is applied within a Bayesian framework and considers that the main effects are represented with cubic smoothing splines which minimize the residual error $\xi_{RV}$, their “smoothness” being introduced via Gaussian priors including smoothing parameters $\lambda_e > 0$, $e = 1, ..., E$, which control the trade-off between small residual errors and smoothness of the main effects.

Because the climate change responses are themselves smooth functions of the year, the residual variability $\xi_{RV}$ is considerably autocorrelated. Moreover, it also usually increases as a function of the horizon, due to the greater variability of the climate projections for larger horizons. Smoothing spline models with correlated errors have been investigated in Wang (1998) who proposes to consider the following model for the ANOVA errors:

$$ \xi_{RV} \sim N(0, \delta_{RV} W_{\rho}), \qquad (2) $$

where $\delta_{RV}$ is a parameter representing the variance of $\xi_{RV}$ and $W_{\rho} = VC_{\rho}V$ is the variance-covariance matrix of $\xi_{RV}$, decomposed into a diagonal matrix $V$ specifying the weights that represent the heteroscedasticity, and $C_{\rho}$ is the matrix of correlations between error terms. Wang (2011) discusses possible parametric
choices for $V$ and $C_\rho$. In this study, in a simple manner, we prescribe the choice of $V$ in order to introduce a linear evolution of the standard deviation of $\xi_{RV}$ with the horizon:

$$V = \mathbb{1}_n \otimes \text{diag}(\sqrt{1/T}, \sqrt{2/T}, \dots, 1), \qquad (3)$$

where $\otimes$ denotes the kronecker product and $\text{diag}(\sqrt{1/T}, \sqrt{2/T}, \dots, 1)$ is a diagonal matrix specifying the weights given to the horizons $x^{\text{year}}$, such that the variance of residual errors for the horizon $x_1^{\text{year}}$ is $\delta_{RV}/T$, and $\delta_{RV}$ for the last horizon $x_T^{\text{year}}$.

A simple autoregressive model of order 1 is considered for the errors corresponding to one climate change response, and are considered independent otherwise, such that:

$$C_{\rho} = \mathbb{1}_n \otimes C_{\rho}^{\text{KMS}}, \qquad (4)$$

where $C_\rho^{\text{KMS}}$ is the Kac-Murdock-Szegö matrix (Kac et al., 1953) of size $T$. This symmetric Toeplitz matrix defines the cross-correlations between the error terms corresponding to a climate projection, and can be written as follows:

$$C_{\rho}^{\text{KMS}} = \begin{pmatrix} 1 & \rho & \rho^2 & \dots & \rho^{T-1} \\ \rho & 1 & \rho & \dots & \rho^{T-2} \\ \rho^2 & \rho & 1 & \dots & \vdots \\ \vdots & \vdots & \vdots & \ddots & \rho \\ \rho^{T-1} & \rho^{T-2} & \dots & \rho & 1 \end{pmatrix}.$$

# 3 BAYESIAN INFERENCE

The different unknown quantities {$\phi_1, \dots, \phi_n, \delta_{IV,1}, \dots, \delta_{IV,n}, \theta_1, \dots, \theta_E, \lambda_1, \dots, \lambda_E, \delta_{RV}, \rho$} are estimated within a Bayesian framework for the following reasons. First, several contributions show that natural Bayesian interpretations can be obtained for the application of smoothing splines (Kimeldorf and Wahba, 1971; Wahba, 1990; Speckman and Sun, 2003) or SS-ANOVA models (Reich et al., 2009; Cheng and Speckman, 2012). Moreover, the methods often applied to select the smoothing parameters (generalized maximum likelihood, generalized cross-validation, and unbiased risk) generally assume independent observations (Wang, 1998). In the context of autocorrelated and heteroscedastic errors, the application of these methods (Wang, 2011) can become unstable and our attempts with the R package assist failed (i.e. the estimates did not converge). Finally, one important advantage is that data augmentation (Tanner and Wong, 1987) is easily implemented within a Bayesian framework, which can be used to treat incomplete designs of climate experiments (Déqué et al., 2007; Evin et al., 2019; Christensen and Kjellström, 2021).

## 3.1 Priors

Without loss of generality, and in agreement with the application shown in Section 4, let us consider here a MME obtained for one emission scenario, and several GCMs and RCMs. For each element $j = 1, \dots, n \times T$ of $\phi^*$, the concatenated vector of $n$ climate change responses $\phi_i^*$, the corresponding “horizon” $x_j^{\text{year}}$, type of GCM $x_j^{\text{GCM}}$ or RCM $x_j^{\text{RCM}}$ can be used to partition the different sources of uncertainties. We aim at decomposing the different climate change response using the SS-ANOVA framework given in Eq. (1) where:

* $\theta_1$ is the main effect corresponding to the grand mean change response, shared by all simulation chains,

* $\theta_2$ and $\theta_3$ are the main effects corresponding to the evolution of the GCM and RCM as a function of the horizon,

with the corresponding priors:

$$\theta_e | \delta_{RV}, \lambda_e \sim N(0, \frac{\xi_{RV}}{\lambda_e} \Sigma_e), \qquad (5)$$
where the covariance matrices $\Sigma_e$ are referred to as “reproducing kernels”, i.e. specific covariance matrices which respect the properties of reproducing kernel Hilbert space (Wahba, 1990; Gu, 2013) and are obtained from the possible values of the predictors. Each element of $\Sigma_e$, $e = 1, 2, 3$ corresponds to the pairs $\{j, k\}$ of indices in $\phi^*$ and are obtained as:

* $k_C(x_j^{\text{year}}, x_k^{\text{year}})$ for $\theta_1$,

* $k_C(x_j^{\text{year}}, x_k^{\text{year}}) \times k_D(x_j^{\text{GCM}}, x_k^{\text{GCM}})$ for $\theta_2$,

* $k_C(x_j^{\text{year}}, x_k^{\text{year}}) \times k_D(x_j^{\text{RCM}}, x_k^{\text{RCM}})$ for $\theta_3$.

where $k_C$ and $k_D$ are functions specifying the reproducing kernels for continuous and discrete predictors, respectively. Many choices are possible for the reproducing kernels, depending on the nature of the predictors (continuous or discrete) but do not seem to have an important impact on the results in this study. Following Cheng and Speckman (2012), for a discrete predictor $\mathbf{x}$ which can take $J$ different values, we consider the following reproducing kernel function $k_D$ between two elements $x_i$ and $x_j$ of $\mathbf{x}$:

$$k_D(x_i, x_j) = \frac{J-1}{J} \mathbb{1}_{x_i=x_j} - \frac{1}{J} \mathbb{1}_{x_i \neq x_j},$$

where $\mathbb{1}$ is the indicator function. For a continuous predictor which has been scaled in $[0, 1]$, we consider the following reproducing kernel function $k_C$:

$$k_C(x_i, x_j) = \min(x_i, x_j)^2 \{3\max(x_i, x_j) - \min(x_i, x_j)\} / 6.$$

Since the reproducing kernel $\Sigma_{CR,i}$ is not necessarily positive definite, Cheng and Speckman (2012) advise to take the following spectral decomposition:

$$\Sigma_e = Q_e D_e Q_e',$$

where $Q_e$ is the $nT \times r_e$ matrix of eigenvectors corresponding to the $r_e$ nonzero eigenvalues of $\Sigma_e$ and $D_e$ is the diagonal matrix of $r_e$ nonzero eigenvalues. Let us consider the following vector:

$$v_e | \delta_{RV}, \lambda_e \sim \mathcal{N}(0, \frac{\delta_{RV}}{\lambda_e} D_e), e = 1, \dots, E,$$

such that $\theta_e = Q_e v_e$ has the prior (5).

Gamma priors are considered for the smoothing parameters $\lambda_e$, $e = 1, \dots, E$:

$$\lambda_e \sim \text{Gamma}(1/2, 2b_\lambda),$$

the corresponding mean being a fixed hyperparameter $b_\lambda$ and the variance is $b_\lambda^2 / 2$ (i.e. shape parameter equals to $1/2$ and scale parameter equals to $2b_\lambda$).

Finally, noninformative priors for $\delta_{RV}$ and $\rho$ are chosen:

$$\delta_{RV} \sim 1/\delta_{RV}, \qquad (6)$$

$$\rho \sim \text{Unif}[-0.999, 0.999], \qquad (7)$$

where $|\rho| < 1$ constraints errors $\xi_{RV}$ to follow a stationary process.

## 3.2 Conditional distributions and Gibbs sampling

The Bayesian framework proposed in this study is implemented with the Gibbs sampling strategy. Indeed, the conditional distributions are easily expressed for all quantities. However, one remaining issue is that the estimation of the smoothing parameters $\lambda_e$ and the autocorrelation parameter $\rho$ are strongly dependent, as discussed in Wang (2011). Moreover, the full conditional distributions are numerically intractable and would require multiple numerical matrix inversions for each step of the Markov chain Monte Carlo (MCMC) sampling. We rely here on a pragmatic multi-step approach, where we first infer all unknown quantities ignoring the heteroscedasticity and autocorrelation of the residual errors by taking $W_\rho = \mathbb{1}_{n \times T}$ in Eq. 2. We then infer the heteroscedasticity / autocorrelation parameters $\delta_{RV}$ and $\rho$ conditionnally on the parameters estimated in the first step. Finally, in a third step, we reestimate the SS-ANOVA effects with these heteroscedasticity / autocorrelation parameters.
### 3.2.1 Step 1

In this first step, we first estimate {$\theta_1^{(1)}, \dots, \theta_E^{(1)}, \lambda_1, \dots, \lambda_E, \delta_{RV}$} assuming that $W_\rho = \mathbb{1}_{n \times T}$. As indicated above, we actually infer $v_e$ and use the equality $\theta_e^{(1)} = Q_e v_e$. The full conditional distributions are obtained as follows:

$$ \delta_{RV}^{(1)} | \phi^*, \theta_1, \dots, \theta_E, \lambda_1, \dots, \lambda_E \sim IGa \left( \frac{n \times T + \sum_e r_e}{2}, \frac{1}{2} \bar{\theta}_e^{(1)'} \bar{\theta}_e^{(1)} + \frac{1}{2} \sum_e \lambda_e v_e' \mathbf{D}_e^{-1} v_e \right), \quad (8) $$

$$ \lambda_e | \delta_{RV}^{(1)}, v_1, \dots, v_E \sim Ga \left( \frac{r_e + 1}{2}, \frac{1}{2} \delta_{RV}^{(1)} v_e' \mathbf{D}_e^{-1} v_e + \frac{1}{2b_e} \right), \quad (9) $$

$$ v_e | \phi^*, \lambda_e, \theta_1, \dots, \theta_E, \delta_{RV}^{(1)} \sim N \left( \lambda_e \mathbf{C}_e^{(1)} \bar{\theta}_e^{(1)}, \delta_{RV}^{(1)} \mathbf{C}_e^{(1)} \right), \quad (10) $$

where

$$ \theta_e^{(1)} = Q_e v_e, \quad (11) $$

$$ \bar{\theta}_e^{(1)} = \phi^* - \sum_k \theta_k^{(1)}, \quad (12) $$

$$ \tilde{\theta}_e^{(1)} = \phi^* - \sum_{k \neq e} \theta_k^{(1)}, \quad (13) $$

$$ \mathbf{C}_e^{(1)} = (\mathbb{1}_{r_e} + \lambda_e \mathbf{D}_e^{-1})^{-1}. \quad (14) $$

We remind that we first obtain the climate change response $\phi^*$ using cubic smoothing splines. Following the Gibbs sampling strategy, we sample iteratively 6000 draws from these conditional distributions, the first iteration being sampled from the prior distribution if the quantity has not been sample yet. The averages of the last $M = 5000$ draws are retained as point estimates for $\theta_e^{(1)}, e = 1, \dots, E$ and denoted as $\hat{\theta}_e^{(1)}, e = 1, \dots, E$. Similarly, we obtain estimates $\hat{\lambda}_e, e = 1, \dots, E$ for the smoothing parameters.

### 3.2.2 Step 2

In the second step, we estimate $\delta_{RV}$ and $\rho$ by taking into account the autocorrelation and heteroscedasticity of the residual errors:

$$ \delta_{RV}^{(2)} | \dots \sim IGa \left( \frac{n \times T}{2}, \frac{1}{2} \left( \phi^* - \sum_k \hat{\theta}_k^{(1)} \right)' \mathbf{W}_\rho^{-1} \left( \phi^* - \sum_k \hat{\theta}_k^{(1)} \right) \right), \quad (15) $$

$$ \log p(\rho) | \dots \sim -\frac{n \times T}{2} \log \delta_{RV}^{(2)} - \frac{1}{2} \log |W_\rho| - \frac{1}{2} \delta_{RV}^{(2)} (\phi^* - \sum_k \hat{\theta}_k^{(1)})' W_\rho^{-1} (\phi^* - \sum_k \hat{\theta}_k^{(1)}). \quad (16) $$

where $|W_\rho| = (1 - \rho^2)^{(n-1) \times T}$ (Kac et al., 1953). We rely on the Metropolis-Hastings algorithm to sample the full conditional distribution of $\rho$. Similarly than in the first step, the full conditional distributions are sampled iteratively and the last $M = 5000$ draws are retained in order to obtain point estimates $\hat{\rho}$ and $\hat{\delta}_{RV}^{(2)}$ as the average of these $M$ draws.

### 3.2.3 Step 3

Finally, we re-estimate the SS-ANOVA effects taking into account heteroscedastic and autocorrelated errors. The full conditional distribution is:

$$ \theta_e^{(3)} | \dots \sim N \left( \hat{\lambda}_e C_e^{(3)} \mathbf{W}_\rho^{-1} \bar{\theta}_e^{(3)}, \hat{\delta}_{RV}^{(2)} C_e^{(3)} \right), \quad (17) $$

with

$$ \tilde{\theta}_e^{(3)} = \phi^* - \sum_{k \neq e} \theta_k^{(3)}, \quad (18) $$

$$ C_e^{(3)} = (\mathbf{W}_\rho^{-1} + \hat{\lambda}_e \Sigma_e^{-1})^{-1}. \quad (19) $$

$C_e^{(3)}$ is obtained using a spectral decomposition of $\mathbf{W}_\rho^{-1} + \hat{\lambda}_e \Sigma_e^{-1}$, where $\mathbf{W}_\rho^{-1}$ is obtained analytically, and where $\Sigma_e^{-1}$ has already been obtained in step 1 using a spectral decomposition.
# 4 APPLICATION TO AN ENSEMBLE OF REGIONAL CLIMATE PROJECTIONS

## 4.1 Climate projections

For illustration, the SS-ANOVA approach is applied to transient climate simulations of regional winter temperature. The MME used in this study is composed of $n = 20$ simulations obtained from the CMIP5-EUROCORDEX experiment (Jacob et al., 2014; Vautard et al., 2021) for all combinations of 4 GCMs and 5 RCMs (see Table 1). Simulation chains are composed of historical runs for the periods 1971-2005, and of future runs for the period 2006-2099 obtained with the emission scenario RCP8.5 (van Vuuren et al., 2011). Our application focuses on mean temperature in winter averaged over the large Central Europe (CEU) region (land and sea points) considered in the IPCC SREX report (Seneviratne et al., 2012) which cover most of European countries above the 45th parallel north at the exception of Norway, Sweden, Finland, Denmark and United Kingdom.

<table><caption>**Table 1.** Characteristics of the scenarios of the EURO-CORDEX climate projection ensemble.</caption><thead><tr><th>GCM</th><th>Member</th><th>RCM</th><th>Version</th></tr></thead><tbody><tr><td>CNRM-CM5</td><td>r1ilp1</td><td>RACMO22E</td><td>v2</td></tr><tr><td>CNRM-CM5</td><td>r1ilp1</td><td>CCLM4-8-17</td><td>v1</td></tr><tr><td>CNRM-CM5</td><td>r1ilp1</td><td>HIRHAM5</td><td>v2</td></tr><tr><td>CNRM-CM5</td><td>r1ilp1</td><td>RCA4</td><td>v1</td></tr><tr><td>CNRM-CM5</td><td>r1ilp1</td><td>REMO</td><td>v1</td></tr><tr><td>EC-EARTH</td><td>r12ilp1</td><td>RACMO22E</td><td>v1</td></tr><tr><td>EC-EARTH</td><td>r12ilp1</td><td>CCLM4-8-17</td><td>v1</td></tr><tr><td>EC-EARTH</td><td>r12ilp1</td><td>HIRHAM5</td><td>v1</td></tr><tr><td>EC-EARTH</td><td>r12ilp1</td><td>RCA4</td><td>v1</td></tr><tr><td>EC-EARTH</td><td>r12ilp1</td><td>REMO</td><td>v1</td></tr><tr><td>HadGEM2-ES</td><td>r1ilp1</td><td>RACMO22E</td><td>v2</td></tr><tr><td>HadGEM2-ES</td><td>r1ilp1</td><td>CCLM4-8-17</td><td>v1</td></tr><tr><td>HadGEM2-ES</td><td>r1ilp1</td><td>HIRHAM5</td><td>v2</td></tr><tr><td>HadGEM2-ES</td><td>r1ilp1</td><td>RCA4</td><td>v1</td></tr><tr><td>HadGEM2-ES</td><td>r1ilp1</td><td>REMO</td><td>v1</td></tr><tr><td>MPI-ESM-LR</td><td>r1ilp1</td><td>RACMO22E</td><td>v1</td></tr><tr><td>MPI-ESM-LR</td><td>r1ilp1</td><td>CCLM4-8-17</td><td>v1</td></tr><tr><td>MPI-ESM-LR</td><td>r1ilp1</td><td>HIRHAM5</td><td>v1</td></tr><tr><td>MPI-ESM-LR</td><td>r1ilp1</td><td>RCA4</td><td>v1a</td></tr><tr><td>MPI-ESM-LR</td><td>r1ilp1</td><td>REMO</td><td>v1</td></tr></tbody></table>

## 4.2 Climate change responses

Climate change responses $\phi_i^*$, $i = 1, \dots, n = 20$ are obtained using cubic smoothing splines (implemented by the function `smooth.spline` in R software, R, 2017) to each simulation chain $i$. A high smoothing parameter (`spar` argument of `smooth.spline` equals to 1) is chosen in order to avoid including decennial variability into these fitted forced responses. In this study, the year $c = 1999$ is retained as the reference year and absolute changes are considered.

Figure 1 illustrates the climate change responses $\phi_i^*$ for this MME. The absolute temperature changes are equal to 0 in 1999 by construction, and gradually increase up to between 4°C and 6°C in 2099. The different simulations for each GCM clearly show that some GCMs lead to higher temperature changes (CNRM-CM5 and HadGEM2-ES) than the other (EC-EARTH, MPI-ESM-LR).

## 4.3 ANOVA and SS-ANOVA decompositions

The Bayesian SS-ANOVA approach is applied by sampling 6,000 MCMC draws of all unknown quantities, i.e. the smoothing spline effects $\theta_e^*$, the residual variability $\delta_{RV}$, the autocorrelation of the residual errors $\rho$ and the smoothing parameters $\lambda_e$. Informative priors are used for the smoothing parameters, the hyperparameter $b_\lambda$ being fixed to $10^{-6}$. Convergence is reached quickly, and a burn-in period composed
<figure><img src="image_3.png" /><figcaption>Figure 1. Absolute temperature change (climate change responses $\phi_i^*$) for the n = 20 simulation chains of the MME between each year of the period 2000-2099 and the reference year c = 1999, with one color by GCM.</figcaption></figure>

of the first 1,000 draws appears to be reasonable. The last 5,000 draws are retained as representative draws from the different posterior distributions.

The SS-ANOVA approach is compared to a simple ANOVA approach where the following linear model can be expressed as follows for a simulation chain i:

$$ \phi_i^*(t) = \mu(t) + \alpha_r(t) + \beta_g(t) + \xi_{RV}(t), \qquad (20) $$

where $\mu(t)$ is mean climate change response, $\alpha_r(t)$ is the main RCM effect for the RCM model r used in the simulation chain i, and $\beta_g(t)$ is the main GCM effect for the GCM model g used in the simulation chain i. The residual terms $\xi_{RV}(t)$ are assumed to be independent and identically distributed (i.i.d.) over all GCMs and RCMs, and to follow normal distributions, with mean 0 and variance $\sigma_{RV}^2(t)$. This linear model is applied to each time step in turn and does not exploit the time dependency of the climate change response between consecutive time steps. This simple ANOVA approach is close to the methodology proposed by Hawkins and Sutton (2009), the only difference being a cubic spline used here instead of 4-order polynomial function in (Hawkins and Sutton, 2009). It is implemented here with the function 1m in R software and is referred to as ANOVA-TI hereafter to indicate the time independence of this approach.

## 4.4 Main effects

Figure 2 compares the GCM and RCM main effects for the ANOVA-TI and SS-ANOVA approaches. Concerning the SS-ANOVA approach, median estimated effects are indicated by thick lines. The uncertainty related to the estimation of each individual effect is indicated by 95% credible intervals obtained from the corresponding posterior distributions. Mean estimated effects obtained with the two approaches are very similar. These results highlight the discrepancies between two groups of GCMs CNRM-CM5 and HadGEM2-ES versus EC-EARTH and MPI-ESM-LR, the former leading to higher temperature changes than the latter.

## 4.5 Smoothing parameters and parameters related to the residual variability

Figure 3 shows the posterior distributions of the smoothing parameters $\lambda_e$ and of the parameters $\delta_{RV}$ and $\rho$ of the residual variability. The posterior distributions of the smoothing parameters indicated well-identified parameters, the smoothing parameter $\lambda_3$ related to the grand ensemble mean being clearly lower than $\lambda_1$ and $\lambda_2$, indicating that the smooth effect $\theta_{\lambda_3}$ for the grand mean is clearly smoother than the two other main effects. The autocorrelation of the residual parameters is rather high, with a mode
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Comparison of the GCM (top plots) and RCM (bottom plots) effects for the ANOVA-TI (left plots) and SS-ANOVA (right plots) approaches. These main effects can be interpreted as absolute temperature changes (°C) attributable to each individual climate model, compared to the year 1999. For the SS-ANOVA approach, the uncertainty related to the estimation of each individual effect is represented by colored intervals (95% credible intervals) around median values (thick lines), obtained from the posterior distributions.</figcaption></figure>

around 0.9, indicating as expected very autocorrelated residual errors due to the autocorrelation in climate
change responses **phi**\*

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Posterior distributions of (a) the smoothing parameters $\lambda_e$, $e = 1, \dots, 3$, (b) the variance of the residual variability at time $T$, $\delta_{RV}$ and (c) the autocorrelation parameter $\rho$ of the residual errors.</figcaption></figure>
## 4.6 Residual variance

Figure 4 shows the standard deviation of the residual errors obtained with benchmark approach ANOVA-TI and the SS-ANOVA approach proposed in this study. The ANOVA-TI approach applies a linear model for each year and does not assume a particular form of evolution for the standard deviation of the residual errors, whereas this evolution is linear for the SS-ANOVA approach. The residual variability obtained with ANOVA-TI seems to increase linearly until 2060 and becomes constant, which creates a discrepancy with the SS-ANOVA approach.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Standard deviation of the residual errors obtained from the approaches ANOVA-TI (linear ANOVA model, time independence) and SS-ANOVA.</figcaption></figure>

# 5 DISCUSSION AND CONCLUSION

In this study, we follow the “time-series” approach (Hingray and Saïd, 2014) which consists in extracting the climate response using a prescribed trend model. The climate responses are estimated using flexible cubic smoothing splines, the degree of smoothing being fixed *a priori* in order to obtain smooth and gradual climate responses. The climate responses are used to compute anomalies between future and reference climate periods, so-called climate change responses. This study aims at developing a statistical framework based on a smoothing spline analysis of variance (SS-ANOVA) to estimate mean expected changes and main climate change effects from an ensemble of climate change responses. We do not discuss here the assumptions concerning the extraction of the climate response, the degree of flexibility for this climate trend being highly subjective. Another aspect which is not discussed here is the internal variability resulting from the natural variability of the climate, often defined as the difference between the climate response and the raw climate projections. This source of uncertainty is irreducible and can be dominant for some indicators (e.g. seasonal precipitation, see Evin et al., 2021, in Europe).

Our Bayesian analysis considers noninformative priors for most inferred quantities (see Section 2). Informative priors are used for the smoothing parameters $\lambda_e$ using a Gamma distribution parametrized with the hyperparameter $b_\lambda = 10^{-6}$. Larger hyperparameter values lead to constrained estimates, and the advantage of using flexible smoothing splines is lost (results not shown). More advanced elicitation methods could be considered, for example using an “effective degrees of freedom” (Hastie et al., 2009; Cheng and Speckman, 2012).

The “time-series” approach applied here can be applied to any climate ensemble and does not require multiple members for each element of the simulation chain (e.g. scenario/GCM/RCM). As such,
it has often been applied to large MMEs which embrace all available climate models with only one
realization is available for some scenario/model combinations, typically CORDEX ensembles (Evin et al.,
2021) or some CMIP ensembles (Lehner et al., 2020). Because the climate change responses evolve
smoothly, the climate change responses are highly autocorrelated. The SS-ANOVA framework proposed
in this study formalize the representation of the autocorrelation of the residual errors, together with their
heteroscedasticity. However, the proposed framework has several limitations:

* We rely on the three-step inference procedure in order to infer the different parameters of the statistical framework, the dependence between some of the parameters is thus ignored,

* Simple forms are taken for the representation of the autocorrelation (AR1 model) and heteroscedasticity (linear evolution of the standard deviation of the residual errors). The linear form for the heteroscedasticity might not be adapted for other applications.

The added complexity of the SS-ANOVA model in comparison of a simple linear model (ANOVA-TI) puts into question the relevance of this framework.

CODE AND DATA AVAILABILITY

The code used for this paper is available through the R-package qualypsoss, which can be down-
loaded from https://CRAN.R-project.org/package=qualypsoss. The results shown in
this study are obtained with the script available at https://github.com/guillaumeevin/
QUALYPSOSS/blob/master/demo/appli_CEU_DJF_tas.r. The code used to apply the ap-
proach "ANOVA / TI" is available through the R-package QUALYPSO, which can be downloaded from
https://CRAN.R-project.org/package=QUALYPSO. The ensemble of climate projections
exploited in this study is attached to both packages.

REFERENCES

Cheng, C.-I. and Speckman, P. L. (2012). Bayesian smoothing spline analysis of variance. *Computational Statistics & Data Analysis*, 56(12):3945-3958.

Christensen, O. B. and Kjellström, E. (2021). Filling the matrix: an ANOVA-based method to emulate regional climate model simulations for equally-weighted properties of ensembles of opportunity. *Climate Dynamics*.

Deser, C., Phillips, A., Bourdette, V., and Teng, H. (2012). Uncertainty in climate change projections: the role of internal variability. *Climate Dynamics*, 38(3-4):527-546.

Déqué, M., Rowell, D. P., Lüthi, D., Giorgi, F., Christensen, J. H., Rockel, B., Jacob, D., Kjellström, E., Castro, M. d., and Hurk, B. v. d. (2007). An intercomparison of regional climate simulations for Europe: assessing uncertainties in model projections. *Climatic Change*, 81(1):53–70.

Evin, G., Hingray, B., Blanchet, J., Eckert, N., Morin, S., and Verfaillie, D. (2019). Partitioning Uncertainty Components of an Incomplete Ensemble of Climate Projections Using Data Augmentation. *Journal of Climate*, 32(8):2423–2440.

Evin, G., Somot, S., and Hingray, B. (2021). Balanced estimate and uncertainty assessment of European climate change using the large EURO-CORDEX regional climate model ensemble. *Earth System Dynamics*, 12(4):1543–1569. Publisher: Copernicus GmbH.

Geinitz, S., Furrer, R., and Sain, S. R. (2015). Bayesian multilevel analysis of variance for relative comparison across sources of global climate model variability. *International Journal of Climatology*, 35(3):433–443.

Green, P. J. and Silverman, B. W. (1993). *Nonparametric Regression and Generalized Linear Models: A roughness penalty approach*. Chapman and Hall/CRC, New York.

Gu, C. (2013). *Smoothing Spline ANOVA Models*. Springer Series in Statistics. Springer-Verlag, New York, 2 edition.

Hastie, T., Tibshirani, R., and Friedman, J. (2009). *The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Second Edition*. Springer Series in Statistics. Springer-Verlag, New York, 2 edition.

Hawkins, E. and Sutton, R. (2009). The Potential to Narrow Uncertainty in Regional Climate Predictions. *Bulletin of the American Meteorological Society*, 90(8):1095–1107.
Hingray, B., Blanchet, J., Evin, G., and Vidal, J.-P. (2019). Uncertainty component estimates in transient climate projections. *Climate Dynamics*, 53(5):2501–2516.

Hingray, B., Mezghani, A., and Buishand, T. A. (2007). Development of probability distributions for regional climate change from uncertain global mean warming and an uncertain scaling relationship. *Hydrol. Earth Syst. Sci.*, 11(3):1097–1114.

Hingray, B. and Saïd, M. (2014). Partitioning Internal Variability and Model Uncertainty Components in a Multimember Multimodel Ensemble of Climate Projections. *Journal of Climate*, 27(17):6779–6798.

Jacob, D., Petersen, J., Eggert, B., Alias, A., Christensen, O. B., Bouwer, L. M., Braun, A., Colette, A., Déqué, M., Georgievski, G., Georgopoulou, E., Gobiet, A., Menut, L., Nikulin, G., Haensler, A., Hempelmann, N., Jones, C., Keuler, K., Kovats, S., Kröner, N., Kotlarski, S., Kriegsmann, A., Martin, E., Meijgaard, E., Moseley, C., Pfeifer, S., Preuschmann, S., Radermacher, C., Radtke, K., Rechid, D., Rounsevell, M., Samuelsson, P., Somot, S., Soussana, J.-F., Teichmann, C., Valentini, R., Vautard, R., Weber, B., and Yiou, P. (2014). EURO-CORDEX: new high-resolution climate change projections for European impact research. *Regional Environmental Change*, 14(2):563–578.

Kac, M., Murdock, W. L., and Szegö, G. (1953). On the Eigen-values of Certain Hermitian Forms. *Journal of Rational Mechanics and Analysis*, 2:767–800. Publisher: Indiana University Mathematics Department.

Kimeldorf, G. and Wahba, G. (1971). Some results on Tchebycheffian spline functions. *Journal of Mathematical Analysis and Applications*, 33(1):82–95.

Lehner, F., Deser, C., Maher, N., Marotzke, J., Fischer, E. M., Brunner, L., Knutti, R., and Hawkins, E. (2020). Partitioning climate projection uncertainty with multiple large ensembles and CMIP5/6. *Earth System Dynamics*, 11(2):491–508. Publisher: Copernicus GmbH.

Northrop, P. J. and Chandler, R. E. (2014). Quantifying Sources of Uncertainty in Projections of Future Climate. *Journal of Climate*, 27(23):8793–8808.

Paeth, H., Vogt, G., Paxian, A., Hertig, E., Seubert, S., and Jacobeit, J. (2017). Quantifying the evidence of climate change in the light of uncertainty exemplified by the Mediterranean hot spot region. *Global and Planetary Change*, 151:144–151.

R (2017). R: A language and environment for statistical computing. https://www.r-project.org/.

Reich, B. J., Storlie, C. B., and Bondell, H. D. (2009). Variable Selection in Bayesian Smoothing Spline ANOVA Models: Application to Deterministic Computer Codes. *Technometrics*, 51(2):110–120.

Reintges, A., Martin, T., Latif, M., and Keenlyside, N. S. (2017). Uncertainty in twenty-first century projections of the Atlantic Meridional Overturning Circulation in CMIP3 and CMIP5 models. *Climate Dynamics*, 49(5-6):1495–1511.

Seneviratne, S. I., Nicholls, N., Easterling, D., Goodess, C. M., Kanae, S., Kossin, J., Luo, Y., Marengo, J., Mc Innes, K., Rahimi, M., Reichstein, M., Sorteberg, A., Vera, C., Zhang, X., Rusticucci, M., Semenov, V., Alexander, L. V., Allen, S., Benito, G., Cavazos, T., Clague, J., Conway, D., Della-Marta, P. M., Gerber, M., Gong, S., Goswami, B. N., Hemer, M., Huggel, C., Van den Hurk, B., Kharin, V. V., Kitoh, A., Klein Tank, A. M. G., Li, G., Mason, S., Mc Guire, W., Van Oldenborgh, G. J., Orlowsky, B., Smith, S., Thiaw, W., Velegrakis, A., Yiou, P., Zhang, T., Zhou, T., and Zwiers, F. W. (2012). Changes in climate extremes and their impacts on the natural physical environment. In *Managing the Risks of Extreme Events and Disasters to Advance Climate Change Adaptation*, pages 109–230. Cambridge University Press. Num Pages: 122.

Speckman, P. L. and Sun, D. (2003). Fully Bayesian Spline Smoothing and Intrinsic Autoregressive Priors. *Biometrika*, 90(2):289–302.

Tanner, M. A. and Wong, W. H. (1987). The Calculation of Posterior Distributions by Data Augmentation. *Journal of the American Statistical Association*, 82(398):528–540.

van Vuuren, D. P., Edmonds, J., Kainuma, M., Riahi, K., Thomson, A., Hibbard, K., Hurtt, G. C., Kram, T., Krey, V., Lamarque, J.-F., Masui, T., Meinshausen, M., Nakicenovic, N., Smith, S. J., and Rose, S. K. (2011). The representative concentration pathways: an overview. *Climatic Change*, 109(1-2):5.

Vautard, R., Kadygrov, N., Iles, C., Boberg, F., Buonomo, E., Bülow, K., Coppola, E., Corre, L., van Meijgaard, E., Nogherotto, R., Sandstad, M., Schwingshackl, C., Somot, S., Aalbers, E., Christensen, O. B., Ciarlo, J. M., Demory, M.-E., Giorgi, F., Jacob, D., Jones, R. G., Keuler, K., Kjellström, E., Lenderink, G., Levavasseur, G., Nikulin, G., Sillmann, J., Solidoro, C., Sørland, S. L., Steger, C., Teichmann, C., Warrach-Sagi, K., and Wulfmeyer, V. (2021). Evaluation of the Large EURO-CORDEX Regional Climate Model Ensemble. *Journal of Geophysical Research: Atmospheres*,
126(17):e2019JD032344. eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2019JD032344.

Wahba, G. (1990). *Spline Models for Observational Data*. CBMS-NSF Regional Conference Series in Applied Mathematics. Society for Industrial and Applied Mathematics.

Wang, Y. (1998). Smoothing Spline Models with Correlated Random Errors. *Journal of the American Statistical Association*, 93(441):341–348. Publisher: [American Statistical Association, Taylor & Francis, Ltd.].

Wang, Y. (2011). Spline Smoothing with Heteroscedastic and/or Correlated Errors. In *Smoothing Splines*. Chapman and Hall/CRC. Num Pages: 24.

Yip, S., Ferro, C. A. T., Stephenson, D. B., and Hawkins, E. (2011). A Simple, Coherent Framework for Partitioning Uncertainty in Climate Predictions. *Journal of Climate*, 24(17):4634–4643.