A Bayesian hierarchical approach for spatial analysis of climate model bias in multi-model ensembles

Abstract

Coupled atmosphere–ocean general circulation models are key tools to investigate climate dynamics and the climatic response to external forcings, to predict climate evolution and to generate future climate projections. Current general circulation models are, however, undisputedly affected by substantial systematic errors in their outputs compared to observations. The assessment of these so-called biases, both individually and collectively, is crucial for the models’ evaluation prior to their predictive use. We present a Bayesian hierarchical model for a unified assessment of spatially referenced climate model biases in a multi-model framework. A key feature of our approach is that the model quantifies an overall common bias that is obtained by synthesizing bias across the different climate models in the ensemble, further determining the contribution of each model to the overall bias. Moreover, we determine model-specific individual bias components by characterizing them as non-stationary spatial fields. The approach is illustrated based on the case of near-surface air temperature bias in the tropical Atlantic and bordering regions from a multi-model ensemble of historical simulations from the fifth phase of the Coupled Model Intercomparison Project. The results demonstrate the improved quantification of the bias and interpretative advantages allowed by the posterior distributions derived from the proposed Bayesian hierarchical framework, whose generality favors its broader application within climate model assessment.

Appendix: Choice of priors

In this section, we provide details of prior and hyperparamters choices. All priors are approximately non-informative. For the error variances $\{\sigma ^2_{\varepsilon ,j}: j=1\dots , 6\}$, we assign the uniform prior on the standard deviation scale $\sigma _{\varepsilon ,j}\sim \hbox {Unif}(a, b)$ by choosing $a=0$ and $b=10^2$ for each j independently. Accordingly, the error variances, which are proportional to $(b-a)^2$, are very large so that the priors are approximately non-informative. For $\{\tau ^2_j: j=1\dots , 6\}$, we use a Half-Cauchy (HC) prior, which is a conditionally conjugate family of a half t distribution (Gelman 2006). The Half t distribution corresponds to the absolute value of a Student-t distribution centered at zero, whose probability distribution is proportional to

$$\begin{aligned} \left( 1+\frac{1}{\hbox {df}}\left( \frac{\tau _j}{\theta }\right) ^2\right) ^{-\left( \hbox {df}+1\right) /2} \end{aligned}$$

(8)

with two parameters: degrees of freedom $\hbox {df}$ and scale parameter $\theta$. We obtain the proper HC probability distribution for $\tau _j$ as a special case of (8) by setting $\hbox {df}=1$,

$$\begin{aligned} p(\tau _j)\propto \left( \theta ^2+\tau ^2_j\right) ^{-1},\qquad j=1,\ldots ,6 \end{aligned}$$

we specify priors for $\tau _j$ as $\tau _j\sim \hbox {HC}(\theta )$, independently for each j. Large but finite value of the scale parameter $\theta$ represents an approximately non-informative prior distribution. In the limit $\theta$ $\rightarrow$ $\infty$ this becomes a uniform prior density on $p(\tau _j)$. For our analysis, we set $\theta =30$. To choose a prior for the $p\times p$ covariance matrix $\mathbf{G}$, the variances $G_1,\ldots , G_p$ and the pair-wise covariances $G_{kl}:k,l=1,\ldots ,p$ must be explicitly specified. One way to achieve this is to use the separation technique (Gelman and Hill 2006; O’Malley and Zaslavsky 2008)

$$\begin{aligned} \mathbf{G}=\varGamma \mathbf Q \varGamma \end{aligned}$$

where $\varGamma$ is the diagonal matrix with diagonal elements $\omega ^2_1,\ldots , \omega ^2_p$ and $\mathbf{Q}$ is new $p \times p$ covariance matrix. The role of the new parameters $\omega ^2_k$ and $\mathbf{Q}$ is to derive appropriately scaled priors for the variances and pair-wise covariances related to $\mathbf{G}$. We assign proper uniform prior on $\omega ^2_k\sim \hbox {Unif}(0,10^{2})$ independently for each k. The covariance component $\mathbf{Q}$ is given the inverse Wishart distribution $\hbox {IW}(p+1,\mathbf{I}_p)$. The two parameters degrees of freedom $p+1$ and the identity matrix $\mathbf{I}_p$ fully determine the distribution. The variances and pair-wise covariances associated to $\mathbf{G}$ are then obtained as $G_k=\omega ^2_k Q_p$ and $G_{kl}=\omega _k\omega _lQ_{kl}$. To make inference, we require the standard deviations $|G_k|^{1/2}$ and correlations $\rho _{kl}$

$$\begin{aligned} |G_k|^{1/2}=|\omega _k|\sqrt{Q_k}\quad \hbox {and}\quad \rho _{kl}=\frac{G_{kl}}{|G_k|^{1/2}|G_l|^{1/2}}, \quad k,l=1,\ldots , p \end{aligned}$$