An evaluation of experimental decadal predictions using CCSM4

Abstract

This study assesses retrospective decadal prediction skill of sea surface temperature (SST) variability in initialized climate prediction experiments with the Community Climate System Model version 4 (CCSM4). Ensemble forecasts initialized with two different historical ocean and sea-ice states are evaluated and compared to an ensemble of uninitialized coupled simulations. Both experiments are subject to identical twentieth century historical radiative forcings. Each forecast consists of a 10-member ensemble integrated over a 10-year period. One set of historical ocean and sea-ice conditions used for initialization comes from a forced ocean-ice simulation driven by the Coordinated Ocean-ice Reference Experiments interannually varying atmospheric dataset. Following the Coordinated Model Intercomparison Project version 5 (CMIP5) protocol, these forecasts are initialized every 5 years from 1961 to 1996, and every year from 2000 to 2006. A second set of initial conditions comes from historical ocean state estimates obtained through the assimilation of in-situ temperature and salinity data into the CCSM4 ocean model. These forecasts are only available over a limited subset of the CMIP5 recommended start dates. Both methods result in retrospective SST prediction skill over broad regions of the Indian Ocean, western Pacific Ocean and North Atlantic Ocean that are significantly better than reference skill levels from a spatio-temporal auto-regressive statistical model of SST. However the subpolar gyre region of the North Atlantic stands out as the only region where the CCSM4 initialized predictions outperform uninitialized simulations. Some features of the ocean state estimates used for initialization and their impact on the forecasts are discussed.

Appendices

Appendix 1: Details of the two initialization methods

The ocean–sea-ice states used to initialize the HDInit coupled model forecasts came from an integration of the nominal 1$^\circ$ resolution POP2 ocean model forced with the Coordinated Ocean-Ice Reference Experiment sinterannually varying atmospheric data set (CORE-II; Large and Yeager 2009; Danabasoglu 2014) from 1948 to 2007. Because this is a single-realization ocean and sea-ice state estimate, ensembles for the HDInit experiment were generated by initializing the atmosphere from time-staggered atmospheric states from the NoInit twentieth century integration. This set of prediction experiments is described in detail in Yeager et al. (2012).

For the DAInit experiment, the ocean was initialized using an ensemble of historical ocean states obtained from the assimilation of in-situ temperature and salinity observations into the POP2 ocean model. Details of the data assimilation method, data sources, and atmospheric forcing used to generate the ocean state estimate can be found in Karspeck (2013). The ensemble assimilation system was initialized on Jan 1 1998, and hydrographic observations were assimilated daily using an ensemble adjustment Kalman filter (Anderson 2001). The Data Assimilation Research Testbed (DART; Anderson et al. 2009) software was used to implement the ensemble filter. The assimilation was run continuously until Jan 1 2006, supplying initial states for Jan 1 of years 2000, 2001, 2002, 2003, 2004, 2005, and 2006. For each of the start dates Jan 1 1975, 1980, 1985, 1990, and 1995 the filter was reinitialized one year preceding each of the forecast start dates, such that one year of in situ observations were assimilated prior to the ocean state estimate being used for forecast initialization. For the DAInit experiment a single atmospheric state taken from a CAM4 solution forced with observed SST and sea-ice was used. The sea-ice state for the DAInit predictions was from a CORE-II experiment where strong restoring to climatological salinity was employed. Sea-ice states for the two experiments were chosen for best compatibility with their respective ocean states. Similarly, while the initial atmospheric state is thought to be unimportant for prediction on this time-scale, all atmospheric initial states were taken at approximately coordinated start-dates for best agreement with the ocean state.

Note that all of the uncertainty in the initial-conditions in the DAInit experiment comes from the ensemble of ocean states, while the HDInit experiment relies on the initial uncertainty in the atmosphere to generate forecast spread.

Appendix 2: Drift correction method

We adopt the cross-validated method of drift correction recommended for CMIP5 experiments outlined in CLIVAR (2011). At each point on the ocean grid, the model drift $d$, corresponding to the $n$ start years $j$ at forecast lead year $t$ is computed as

$$\begin{aligned} d_{j,t} = \frac{1}{n-1}\sum ^n_{i=1} (\langle X_{i,t}\rangle - O_{i,t}),\;\;i \ne j \end{aligned}$$

where $\langle X_{i,t} \rangle$ is the annual and ensemble mean SST forecasts initialized on Jan 1 of year $i$ and $O_{i,t}$ is the HADISST from the corresponding year. Drift-corrected forecasts are then computed as the difference between the raw forecast members and the model drift, i.e. $X_{j,t}-d_{j,t}$.

As formulated above, the forecast at any particular start date does not contribute to the drift estimate that is subsequently applied. While cross-validation is intended to reduce the likelihood of overestimating skill, Gangsto et al. (2013) demonstrate that when the sample size is low cross-validation can sometimes result in a negatively biased skill score.

Appendix 3: A bootstrapped single-ensemble member correlation coefficient and an approximate formula for its expected value

In this appendix we describe a bootstrapping procedure for building a distribution of single-member correlation coefficients that accounts for the spread in the ensemble forecasts. This distribution is used as the skill-score in the body of the paper. We also give an approximate formula for the expected value of the distribution. Note that the specifics of this bootstrap were developed for this assessment. The reader is referred to Efron and Tibshirani (1993) for a general description of bootstrapping and its strengths and limitations.

Bootstrapped samples of correlation coefficients associated with an initialized forecast experiment are generated in the following way: For each of the $n$ start dates, one of the 10 forecast ensemble members is chosen at random. The SST averaged over a box-region (or simply taken at a grid-point) is then computed for the 2–4 or 6–9 years lead times. Using this set of $n$ single-member forecasts, a correlation coefficient is computed with (1). These ensemble members are then returned to the pool of ensembles associated with their start date, so that they may be randomly chosen again. This process is repeated 2,000 times, providing samples that can be used to build a bootstrapped distribution of $r$ associated with the initialized prediction experiment.

Here we derive an approximate formula for the expected value of the distribution of the single-member anomaly correlation coefficient, $r$. We can treat $X_{j,\tau }$ from (1) as a Gaussian random variable, with mean assumed to be the forecast ensemble mean, and variance assumed to be the forecast ensemble variance (denoted $\sigma ^2_{j,\tau }$). This implies that $r$, which is a function of $X_{j,\tau }$, is also a random variable. Of course, the sampling distribution of $r$ is not normally distributed, but we can use standard definitions of the statistical expectance of functions of normally distributed variables (Ross 2007) to derive the expected value of the distribution of $r$.

We make the simplifying assumption that $s^2_X$ in (1), which is a random variable by virtue of its functional dependence on $X_{j,\tau }$, can be replaced by a constant $\tilde{s^2_X}$ that is defined as the expected value, i.e.,

$$\begin{aligned} \tilde{s^2_X} \equiv E(s^2_X) = s^2_{\langle X \rangle } + \frac{1}{n}\sum _{j=1}^n \sigma ^2_{i,\tau } = s^2_{\langle X \rangle }+ \overline{\sigma ^2_{\tau }}. \end{aligned}$$

Here $\langle \cdot \rangle$ denotes an ensemble average and the overbar is an average over start dates. Using $\tilde{s^2_X}$, we can approximate the expected value of $r$ as:

$$\begin{aligned} E(r) \approx \frac{1}{n-1}\sum _{j=1}^n \frac{(\; \langle X_{j,\tau }\rangle -\overline{\langle X_{\tau }\rangle }\;) \;\; (O_{j,\tau }-\overline{O_{\tau }})}{s_O \sqrt{s^2_{\langle X \rangle }+ \overline{\sigma ^2_{\tau }}}}. \end{aligned}$$

(2)

We have expressed the expected value of $r$ as a function of the forecast ensemble statistics to highlight the relationship between this single-member based correlation metric and the correlation score as it is commonly computed using the ensemble mean. Due to the forecast ensemble spread, the expected value of $r$ will be smaller than what would be computed if $X_{j,\tau }$ were an ensemble mean forecast at each start date. In the limit that the average ensemble spread ($\overline{\sigma ^2_{\tau }}$) is much smaller than the time-variance of the ensemble mean ($s^2_{\langle X \rangle }$), the expected value of $r$ will be the correlation coefficient as it is commonly computed using the ensemble mean. The approximate relationship in (2), can be in error due to the assumption that $s^2_X$ can be replaced by $\tilde{s^2_X}$, and also because bootstrapping with replacement yields values of $r$ that are not independent samples. In a separate test of the utility of (2), we found that for our experimental setup (2) underestimates the bootstrapped mean value of $r$ by about 3 %. It is thus a reasonable approximate formula.

Appendix 4: Reference statistical forecasts

To create a benchmark for what are meaningful levels of skill, we generate two reference forecasts. Both reference forecasts used in this study are based on representing global SST variability as a first order auto-regressive (AR1) process that includes both spatial and temporal auto-correlations in the statistical model. Modeling the auto-correlation in both time and space is important because we are evaluating SST in regionally and temporally averaged boxes and because there is long-term autocorrelation in the observations and the forecasts.

The multivariate spatial patterns are computed as the leading 15 Empirical Orthogonal Functions (EOF) of the annually averaged HADISST anomaly relative to the 1870–1960 time mean. The leading 15 modes can be used to account for about 88 % of the annual, coarse-resolution global SST variance from 1870 to 1960. We only use the data up until 1960 to avoid overlap with our forecast evaluation period. The multivariate AR1 process for the vector $v$ of time-varying Principal Components (PC) associated with the spatial EOF patterns can be expressed:

$$\begin{aligned} v_{t+1} = Gv_t + \epsilon , \;\;\;\; \epsilon \sim \mathcal {N}(0,\Sigma ^2). \end{aligned}$$

(3)

The vector random variable $\epsilon$ is normally distributed with mean zero and variance $\Sigma ^2$. It is time-uncorrelated and uncorrelated across modes (i.e. $\Sigma ^2$ is a diagonal matrix). The least squares estimate of $G$ can be computed as $G = C_1C_0^{-1}$, where the one year lag covariance is $C_1 = \langle v_{t+1} v_t^T\rangle$ and the lag-zero covariance is similarly $C_0 = \langle v_{t} v_t^T \rangle$. Here $\langle \cdot \rangle$ refers to an average over all years $t$ from 1870 to 1960. The error variance associated with each EOF/PC mode (the elements of $\Sigma ^2$), can be computed directly from the residual (i.e., $v_{t+1}-Gv_t$). Statistical realizations of SST variability that share the spatio-temporal statistical properties of the observed SST anomaly can be generated by integrating (3) forward in time to form PC time series that are then multiplied by the associated EOF spatial patterns.

Samples of correlation coefficients associated with the null hypothesis of no correlation are generated in the following way: two independent 50-year realizations of (3) are generated with randomly generated initial conditions at $t_0$. The SST averaged over a box-region (or at a grid-point) is then computed. Four year averages are computed at intervals set by the prediction experiment. For consistency with these experiments, this are nine 4-year averages at five year intervals (associated with the 1961–2000 start dates) followed by six 4-year averages at one year interval (associated with the 2001–2006 start-dates).⁷ A sample correlation coefficient is then computed with these averaged and subsampled time-series. This process is repeated 2,000 times, providing samples that can be used to build an empirical distribution of $r$ associated with ‘no-skill.’

The AR1 process described above can also be used to make initialized statistical predictions. When used in this way, an AR1 model is similar to the linear inverse model (or LIM) benchmark forecasts outlined by Newman (2013). To make LIM forecasts, the HADISST from the start years of the HDInit and DAInit experiments are projected onto the leading 15 EOFs, and a single ensemble member is forecasted from the resulting PC using (3). These ensemble members have identical initial conditions and are different only because of the random variations in $\epsilon$. Just as in the HDInit and DAInit experiments, correlation coefficients are computed between the statistical forecast and the corresponding observations of HADISST. This process is repeated 2,000 times. Sample correlation coefficients from these LIM predictions are used just as for the HDInit and DAInit coupled model predictions to build an empirical distribution of $r$ associated with LIM.

Appendix 5: Computing the probability that one set of forecasts outperforms another

The probability that one forecast outperforms another, which are tabulated in Tables 1 and 2, are computed in the following way. A normal kernel density function (we use the Matlab ‘ksdensity’ function; Bowman and Azzalini 1997) is used to generate smooth (but discrete) probability density functions from the bootstrapped correlation samples. We then perform a set of discrete integrals to compute the probability $p$ that correlation coefficients from distribution $f$ (e.g., from the HDInit and DAInit experiments) will exceed correlation coefficients from distribution $g$ (e.g, from the ‘no-skill’ or LIM statistical forecasts). This can be written:

$$\begin{aligned} p=\sum _{r=-1}^{1}g(r)\sum _{r'=r}^{1} f(r') \Delta r' \Delta r\;\;\;\;\; (r,r') \in [-1:\Delta r:1] \end{aligned}$$

where we have made use of the fact that correlation coefficients can only range from $-1$ to $1, \Delta r$ and $\Delta r'$ are (identical) bin-widths that define the discrete values of $r$ and $r'$, and $f$ and $g$ are expressed as discrete probability density functions.

In the case of the NoInit experiments, where there are only five samples, a distribution is not computed. Instead, we assume that each ensemble member is equally likely to occur with a probability of 20 %. The probably of outperforming the NoInit experiment can then be written:

$$\begin{aligned} p=\sum _{r} 0.20\sum _{r'=r}^{1} f(r') \Delta r',\;\; r \in [r_1, r_2 ... r_5], \;\;r' \in [-1:\Delta r:1] \end{aligned}$$

where $r_1$ through $r_5$ are the correlation coefficients computed for each of the NoInit ensemble members.