Improving Multi-model Ensemble Probabilistic Prediction of Yangtze River Valley Summer Rainfall

The Yangtze River valley (YRV), the most populous region and the primary agricultural and industrial region in China, lies in the East Asian monsoon region. The strong interannual variability of East Asian summer monsoon leads to frequent occurrence of floods and droughts over the YRV, which causes huge economic loss and affects the freshwater resource of the country (Huang et al., 2003). Although seasonal prediction of YRV summer rainfall is important, the complexity of related tropical and mid-latitude processes and their interactions makes the task quite challenging and with large uncertainty (Tao and Chen, 1987; Wang, 2001; Xiao, 2010; Huang et al., 2012).

General circulation models (GCMs) have played an increasingly important role in seasonal prediction in China since 1990 due to their solid physical foundation, adaptability to the shifts in climate, and ability to handle a wide range of linear and nonlinear interactions (Zeng et al., 1990; Ding et al., 2004; Li et al., 2005; Luo et al., 2013). However, uncertainty in dynamical seasonal prediction, which comes from errors and simplifications in the model as well as uncertainty in the initial conditions, is inevitable (Huang, 1993; Palmer et al., 2004). Correspondingly, predictions are essentially probabilistic, and thus expected to be expressed in a probability density function (PDF) manner (Gneiting, 2008; Doblas-Reyes et al., 2009; Lavers et al., 2009). The multi-model ensemble (MME), whose product is probabilistic, has the potential to combine the strengths of individual models and reduce/estimate the uncertainty in dynamical seasonal prediction (WCRP, 2005; WCRP, 2008; Doblas-Reyes et al., 2009; Alessandri et al., 2011).

Many prediction schemes have been developed for summer rainfall prediction in China to extract prediction information from MME datasets (Feng and Fu, 2007; Ke, 2007; Ke et al., 2009; Wang and Fan, 2009; Sun and Chen, 2012). All these prediction schemes focused on single-value (deterministic) prediction, and thus could not provide uncertainty information for YRV summer rainfall prediction. On the other hand, there have been four commonly-used PDF ensemble schemes: the equal-weighted ensemble scheme with raw model outputs; the equal-weighted ensemble scheme with calibrated model outputs; the multiple linear regression ensemble scheme; and the Bayesian scheme. The first three of these schemes evolved from their deterministic prediction versions (Gneiting et al., 2005; Weigel et al., 2009; Yatagai et al., 2014), while the Bayesian scheme was born for PDF prediction and has in general been found to perform well (Raftery et al., 1997; Coelho et al., 2004; Luo et al., 2007; Wang et al., 2012). Recently, (Li, 2011) assessed the performance of the four schemes in the seasonal PDF prediction of regional summer rainfall over East China based on rainfall hindcasts from the ENSEMBLES. Results showed that all of their PDF predictions of YRV summer rainfall were close to the climatology forecast (the benchmark of PDF prediction), suggesting that the four schemes are not adapted to the PDF prediction of YRV summer rainfall.

ENSEMBLES is an EU-funded seasonal-to-annual MME project (Weisheimer et al., 2009; Alessandri et al., 2011), which provides multi-initial-condition ensemble hindcast datasets from five coupled atmosphere-ocean general circulation models (CGCMs) for research of prediction uncertainty. It succeeds the Development of a European Multimodel Ensemble System for Seasonal to Interannual Prediction (DEMETER) project (Palmer et al., 2004) with improvements in the CGCMs, including better representation of sub-grid scale physical processes and boundary forcing, increased resolution, and more widespread use of assimilation for ocean initialization. (Li et al., 2012) evaluated the predictability of the western North Pacific (WNP) summer monsoon index proposed by (Wang and Fan, 1999) in ENSEMBLES. According to the recommendation of (Wang et al., 2008), the reversed Wang and Fan index can represent well the interannual variability of East Asian summer monsoon circulation. The evaluation results showed that with a one-month lead, the WNP summer monsoon index was highly predictable, with a correlation coefficient of 0.68 between observations and the MME for 1960-2005 (Li et al., 2012). Unfortunately, for summer rainfall over continental East Asia, there were only scattered regions where rainfall was skillfully predicted; the correlation coefficient between observations and the MME-simulated summer rainfall at almost all stations over the YRV was less than 0.3 for the period 1979-2005 (Lu et al., 2012; Li et al., 2012).

In the present study, a PDF prediction scheme is developed for YRV summer rainfall based on ENSEMBLES multi-model hindcasts from 1960 to 2002. The scheme is similar to the Bayesian ensemble prediction scheme in (Coelho et al., 2004), but with major differences in ensemble members (section 3.1) and modeling the variance of the likelihood function (section 3.2). In the present study, the optimized ensemble members are used rather than the model output of target variable used by (Coelho et al., 2004). Also, the revised variance estimation of the likelihood function is a linear function of ensemble spread with a constant item, rather than a constant multiple of ensemble spread as in (Coelho et al., 2004). The remainder of the paper is structured as follows. Section 2 introduces the data and data pre-processing. Section 3 describes the proposed PDF prediction scheme for YRV summer rainfall. Section 4 presents its cross-validation skill, including calibration, resolution, and accuracy. In addition, the performance of the new scheme is compared with those of the equal-weighted ensemble scheme with raw model outputs and the Bayesian scheme from (Coelho et al., 2004). Further discussion and a summary of the key conclusions are presented in section 5.

Summer (June-July-August; JJA) rainfall observations at 160 Chinese rain gauge stations are provided by the National Climate Center, China Meteorological Administration. Summer 500-hPa geopotential height (Z500) observations are from the European Centre for Medium-Range Weather Forecasts (ECMWF) 40-year Reanalysis (ERA-40, Uppala et al., 2005) at a 2.5^° spatial resolution for the period 1960-2002.

The simulations of summer rainfall and Z500 are archived by the ENSEMBLES stream 2 dataset (Weisheimer et al., 2009; Alessandri et al., 2011). The dataset comprises outputs of five CGCMs from the UK Met Office (UKMO), M\'et\'eo-France (MF), ECMWF, Leibniz Institute of Marine Sciences at Kiel University (IFM-GEOMAR), and Euro-Mediterranean Centre for Climate Change (CMCC-INGV) in Bologna. For each CGCM, there are nine ensemble members from different initial ocean conditions with wind stress and SST perturbation. The hindcasts are from 1960 to 2005 and have a spatial resolution of 2.5^° of latitude and longitude. For each year, 7-month-long seasonal forecasts start on the first day of February, May, August, and November. The JJA rainfall data initiated from May 1 are used in the present study.

The observed time series of YRV summer rainfall is derived as follows. First, five adjacent stations along the Yangtze River with standard deviation (D _s) of summer rainfall larger than 2.5 mm d^-1 are selected (Fig. 1a), and their average rainfall is set as the reference time series. Second, the correlation coefficients between the reference time series and rainfall at each of the 160 Chinese stations are calculated. As shown in Fig. 1b, there are 27 stations with significant positive correlation at a 0.05 significance level according to the nonparametric bootstrap test with 1000 random samplings (Efron, 1979). Finally, the average rainfall of the 27 stations is regarded as the observed YRV summer rainfall. Note that, even though it has a D _s larger than 2.5 mm d^-1, Fuyang station (32.93^°N, 115.83^°E) is excluded from the calculation of the reference time series given that it is located in the region of the Huaihe River valley. Some previous studies pointed out that interannual variations of summer rainfall and associated circulation anomalies between the YRV and Huaihe River valley were significantly different (Ma et al., 2011; Hong and Liu, 2012).

Figure 1.  Spatial distribution of (a) standard derivation (D _s) of summer rainfall and (b) correlation between the summer rainfall at 160 Chinese stations and the reference time series. Areas with D _s>2.5 mm d^-1 in (a) and correlation exceeding the 0.05 significance level in (b) are shaded. The black dots are the five stations whose rainfall are used to derive the reference time series in (a) and the representative stations for the YRV in (b). The contour interval is 0.5 mm d^-1 in (a) and 0.2 in (b). Contour levels less than 2 mm d^-1 are omitted in (a).

The model rainfall outputs are bilinearly interpolated to the 27 stations in Fig. 1b, and then the rainfall averaged over the stations is the YRV summer rainfall from ENSEMBLES.

The PDF ensemble prediction scheme includes two components: the optimization of ensemble members, and prediction with the Bayesian scheme. They are described in sections 3.1 and 3.2, respectively. In the present study, a "leave-one-out" cross-validation framework (Wilks, 1995) is adopted for building and evaluating the PDF ensemble prediction model to prevent over-fitting and wasting data. For a certain target year, only the data in other years are used as the training set to build the prediction model. The forecast skill obtained within this framework is called the cross-validation skill.

(Li, 2011) showed that, for YRV summer rainfall prediction, the likelihood function in the Bayesian ensemble prediction scheme was barely capable of extracting prediction information from ENSEMBLES rainfall outputs, and did not show its impressive power in ENSO prediction in (Coelho et al., 2004). One important reason is that the CGCMs in ENSEMBLES simulate YRV summer rainfall poorly, like other models (Zhu et al., 2008; Wang et al., 2009; Sun and Chen, 2012; Luo et al., 2013). To increase the prediction information in ensemble members, we optimize ensemble members before going through the Bayesian ensemble scheme.

It is well-known that CGCMs can reasonably simulate large-scale circulation, such as 500-hPa geopotential height (Z500), but fail to reproduce mid-latitude rainfall over land (Huth, 1999; Wang et al., 2009). Moreover, synchronous Z500 in many regions is significantly correlated with YRV summer rainfall (Fig. 2). Therefore, ensemble members of YRV summer rainfall regressed from CGCM-simulated Z500 may hold more prediction information. In the present study, the ith (i=1,2,3,4,5) ensemble member M_i (i.e. the ith regressed YRV summer rainfall) is the predictand of a simple linear regression with the ith selected Z500 simulation factor C_i as the predictor:

\begin{equation} \label{eq1} \bm{M}_i=k_i\bm{C}_i+\bm{l}_i, \end{equation}

where the regression coefficients k_i and constant vector l_i are estimated from the training set. The five Z500 factors are selected via three steps. First, we calculate the multi-initial-condition ensemble mean of summer Z500 for each CGCM. Second, correlations between simulated and observed Z500 and between simulated Z500 and observed YRV summer rainfall are calculated. For simulated and observed Z500, 10^{°× 10°} domains are used to avoid uncertainty caused by a single grid and to improve the stability of selected circulation factors. Third, to increase the prediction information contained in ensemble members, we select five domains in which the correlations between simulated Z500 and observed YRV summer rainfall are the highest from the domains where Z500 is simulated skillfully (correlation between observed and simulated Z500 is significant at the 95% confidence level, so selected Z500 factors are not spurious). Also required is that the overlap between any two of them in the same model is no more than 1/4, to decrease the interdependence of the selected Z500 factors.

Figure 2.  Locations of the domains of the selected Z500 factors for different target years (boxes) and correlation between the observed YRV summer rainfall and observed synchronous 500-hPa geopotential height (Z500, shading). Shaded areas from light to dark color denote significant correlation at the 90% and 95% confidence levels.

Figure 2 shows that the domains of the selected Z500 factors are located in the tropical Pacific. These optimal circulation factors are simulated skillfully (correlation coefficients are ∼0.6, significant at the 99.9% confidence level), consistent with the results of (Wu and Kirtman, 2005), who reported that the large-scale circulation fields in the tropics could generally be simulated reliably by CGCMs. (Wang et al., 2009) pointed out that the high skill of CGCM simulations for the large-scale circulation in the tropics is because CGCMs could capture the influence of ENSO teleconnection. The ENSO-induced heating excites rapidly propagated equatorial Kelvin and Rossby waves in the eastern Pacific. As a result, the air in the entire tropics warms up and 500-hPa geopotential height rises up during El Ni\ no. The opposite is also true during La Ni\ na. In addition, simulated Z500 in the selected domains is significantly correlated with observed YRV summer rainfall (correlation coefficients are ∼0.6, significant at 99.9% confidence level). This relationship is also presented between the observed YRV summer rainfall and the observed Z500 in the tropical Pacific (see Fig. 2). The possible physical mechanisms are discussed later in section 5.

Observed, modeled, and regressed time series of YRV summer rainfall can pass the Jarque-Bera test of normal distribution (Jarque and Bera, 1980) at the 0.05 significance level, so all ensemble members and target variable satisfy the Bayesian modeling assumption of normality. The Bayesian scheme includes three steps: selection of the prior distribution; modeling of the likelihood function; and determination of the posterior distribution (Coelho et al., 2004).

First, the climatology forecast is adopted as the prior distribution in the present study:

\begin{equation} \theta_t\sim{\rm N}(\mu_{\rm oc},\sigma_{\rm oc}^2) , \end{equation}

where N means normal distribution, and the mean μ _oc and variance σ _oc² are estimated by the sample mean and variance of observations in the training set.

Then, the likelihood \(p(\overlineM_t|\theta_t)\) is modeled by performing a simple linear regression with the average of ensemble members \(\overline\bmM\) as the predictand and observations θ as the predictor:

\begin{equation} \overline{M}_t|\theta_t\sim{\rm N}(a\theta_t+b,cs_t^2+d) , \end{equation}

where a and b are regression coefficients. The variance of the likelihood function is a linear function of ensemble spread s_t² with a constant d, rather than the constant times of ensemble spread as in (Coelho et al., 2004). We add a layer of leave-one-out framework to estimate the regression coefficients c (≥0 required) and d with the residual variance of

\begin{equation} \overline{\bm{M}}=a\bm{\theta}+\bm{b} , (4)\end{equation}

as the predictand and ensemble spread s_t² as the predictor. In Eq. (3), the constant d represents the independent part of residual variance from ensemble spread. When c is equal to 0, d is equal to the variance of residuals estimated from the training set, similar to that used by (Luo et al., 2007) and (Yuan et al., 2013). Otherwise, the dependent part of residual variance on the spread rate can capture the information from the ensemble spread in the target year. For YRV summer rainfall, the dependence of residual variance on ensemble spread is small for both modeled and regressed YRV summer rainfall, so the uncertainty cannot be estimated well if constant multiples of ensemble spread are used for YRV summer rainfall prediction. Note that YRV summer rainfall prediction is not a special case for low spread-error relationship. (Whitaker and Loughe, 1998) demonstrated that even for a perfect ensemble (one in which all sources of prediction uncertainty are sampled correctly), correlation between ensemble spread and error need not be high.

Finally, from Bayesian theorem, the posterior distribution is also normally distributed when the prior distribution and the likelihood function are both normal distributions (Lee, 1997). The normal posterior distribution for the target year t is

\begin{equation} \theta_t|\overline{M}_t\sim{\rm N}(\mu_t,\sigma_t^2) , \end{equation}

where the mean μ_t and variance σ_t² are

\begin{equation} \dfrac{1}{\sigma_t^2}=\dfrac{1}{\sigma_{\rm oc}^2}+\dfrac{a^2}{cs_t^2+d} ,\\ \label{eq4} \dfrac{\mu_t}{\sigma_t^2}=\dfrac{\mu_{\rm oc}}{\sigma_{\rm oc}^2}+\dfrac{a^2}{cs_t^2+d}\left(\dfrac{\overline{M}_t-b}{a}\right). \end{equation}

Equations (6) and (7) state that the precision of the posterior distribution 1/σ_t² is the sum of the precisions of the prior distribution 1/σ _oc² and the ensemble system a²/cs_t²+d, while the posterior ensemble mean μ_t is equal to the weighted average of the prior mean and ensemble mean.

This section presents the cross-validation skill of the new PDF ensemble prediction scheme, including calibration and sharpness in section 4.1, and accuracy (i.e. combined assessment of calibration and sharpness) in section 4.2. Furthermore, its skill is compared with the equal-weighted ensemble scheme (EE) and the Bayesian scheme from (Coelho et al., 2004) with the climatology forecast as the prior distribution (Bayes). The EE and Bayes use the multi-initial-condition ensemble mean of the YRV summer rainfall outputs from five CGCMs as ensemble members.

Calibration and sharpness are two desirable properties of probabilistic prediction; and PDF prediction aims to maximize its sharpness subject to calibration (Raftery et al., 2005; Gneiting et al., 2005, 2007). Calibration is a measure of the statistical consistency between the predictive distribution and observations, while sharpness refers to the concentration of the predictive distributions. Following (Raftery et al., 2005) and Gneiting et al. (2005, 2007), in the present study, we use probability integral transform (PIT) histogram and the average width of the central 95% prediction interval (95% PI) to assess the calibration and sharpness, respectively.

Figure 3.  Probability integral transform (PIT) histogram for (a) the scheme developed in the present study (New), (b) the equal-weighted ensemble (EE) scheme, and (c) the Bayesian scheme in (Coelho et al., 2004) with the climatology forecast as the prior distribution (Bayes).

The PIT value is the predictive cumulative distribution at an observation value (Dawid, 1984; Gneiting et al., 2007). If the predictive PDF is calibrated, the PIT values should be uniformly distributed when the sample size is infinite. As shown in Fig. 3, histograms of the new prediction scheme and Bayes are nearly flat, while that of EE (i.e. the raw ensemble) is obviously hump-shaped. This indicates that the new prediction scheme and Bayes are well-calibrated, while the raw ensemble is over-dispersive. To quantify the uniformity of the PIT histogram, Squared Bias (SB) proposed by (Glahn et al., 2009) is used here. The smaller the SB, the better calibrated the predictive PDF is. As shown in Table 1, the SB of the new scheme is 0.23 mm d^-1, smaller than that of EE (SB=1.26 mm d^-1) and Bayes (SB=0.37 mm d^-1), indicating that the prediction of the new scheme is the best calibrated among the three schemes.

For normally distributed variables, the 95% PI is equal to 2z_0.975σ, where z_0.975 denotes the 0.975 quantile of the standard normal distribution and σ is the standard deviation of the predictive PDF. The smaller the 95% PI, the sharper the predictive PDF is. As shown in Table 1, the 95% PI of the new scheme is 1.78 mm d^-1, smaller than that of EE (3.40 mm d^-1) and Bayes (2.26 mm d^-1), indicating that the PDF prediction of the new scheme is shaper than EE and Bayes.

The temporally averaged continuous ranked probability score (CRPS) is an appropriate skill score for the combined assessment of calibration and sharpness (Gneiting et al., 2007). The CRPS is the integral of the Brier scores at all possible threshold values of the predictand, and refers to the discrepancy between the cumulative distribution functions of predictions and observations. The smaller the CRPS, the more skillful the prediction system is. The climatology forecast is the benchmark of the PDF prediction. The PDF prediction is skillful when it outperforms the climatology forecast. As shown in Fig. 4, the CRPS of the new scheme is 0.57, smaller than that (0.69) for the climatology forecast. Therefore, the new scheme is skillful. The new scheme performs better than EE and Bayes whose CRPSs are 0.67. The Bayes scheme is not better than the EE, though it is sharper and more calibrated. The relative CRPS of the new scheme's prediction to the climatology forecast is 0.83, smaller than the lowest CRPSs of summer rainfall PDF prediction in all regions of East China using the four commonly-used PDF ensemble schemes [0.92 for the northern part and 0.96-0.98 for the other regions (Li, 2011)].

Figure 4.  Temporally averaged continuous ranked probability score (CRPS) for the PDF prediction of different schemes and climatology prediction (the benchmark of PDF prediction).

Both the optimization of ensemble members and the revision of variance modeling of the likelihood function contribute to this higher skill of the new prediction scheme as compared to the Bayes scheme of (Coelho et al., 2004). When only the optimization of ensemble members is used, the calibration and sharpness of the PDF prediction are improved (SB=0.29 mm d^-1; 95% PI=1.90 mm d^-1), and the CRPS is decreased from 0.67 in Bayes to 0.62. The reason is that the correlation between the optimized ensemble members and the target variable is higher (Fig. 5a) and more stable (Fig. 5b). The correlation coefficients increase to ∼0.6 from ∼0.2 between the raw ensemble members and the target variable.

Figure 5.  (a) Average and (b) coefficient of variability (CV) of leave-one-out correlation coefficients (Cor) between observed YRV summer rainfall and the five ensemble members (Member) and ensemble average (Avg) before (Raw) and after (Optimized) the optimization of ensemble members.

Meanwhile, the coefficients of variability (CVs), which measure the instability of correlation, decrease to ∼0.02 from ∼0.15. Because the average of ensemble members is stably and highly correlated to the predictand, the statistical significance of Eq. (4) is improved after the optimization; and correspondingly, the modeling of the likelihood function can extract more prediction information from ensemble members. When the revision of variance modeling of the likelihood function is also included, the calibration, sharpness, and accuracy of the PDF prediction skill are improved further. SB decreases from 0.29 mm d^-1 to 0.23 mm d^-1, the 95% PI decreases from 1.90 mm d^-1 to 1.78 mm d^-1, and the CRPS decreases from 0.62 to 0.57 (Table 1).

In this study, we develop an MME PDF prediction scheme for YRV summer rainfall based on hindcasts from ENSEMBLES. The newly-developed scheme is similar to the Bayesian ensemble scheme in (Coelho et al., 2004), but with major differences in ensemble members and the variance modeling of the likelihood function. The optimized ensemble members are regressed YRV summer rainfall, which are obtained by a simple linear regression with selected factors of synchronous 500-hPa geopotential height as predictors, rather than raw model rainfall outputs. The revised variance estimation of the likelihood function is a linear function of ensemble spread with a constant, rather than constant multiple of ensemble spread as in (Coelho et al., 2004). The constant in the linear function represents the independent part of the variance of the likelihood function from the ensemble spread, while the regressed coefficient represents the time-varying dependent part of the ensemble spread. The cross-validation prediction skill of 1960-2002 YRV summer rainfall is used to assess its performance. Results show that PDF prediction of the new scheme is well-calibrated, sharp, and skillful, and obviously outperforms the EE and Bayes schemes.

In the present study, we select five Z500 simulation factors for each target year that are located in the tropical Pacific and highly correlated with observed YRV summer rainfall. Similarly, correlation is also significant for observed Z500 in the tropical Pacific and YRV summer rainfall (Fig. 2). This relationship between tropical circulation and YRV summer rainfall is probably attributable to the effect of ENSO. On the one hand, tropospheric temperature increases over the tropics following El Ni\ no (Yulaeva and Wallace, 1994; Sobel et al., 2002), which can persists into the following summer mostly due to SST warming over the tropical Indian and Atlantic oceans (Kumar and Hoerling, 2003; Lau et al., 2005; Xie et al., 2009). The tropospheric warming increases Z500 in the tropics. On the other hand, rainfall is significantly enhanced in the YRV during the summer following El Ni\ no (e.g. Huang and Wu, 1989; Huang et al., 2004), associated with a lower-tropospheric anticyclonic anomaly over the WNP (Wang et al., 2000, 2003; Li et al., 2007) due to the local air-sea interaction (Wang et al., 2000; Lau and Nath, 2003) and the capacitor effect of tropical Indian Ocean warming (Xie et al., 2009). In short, during the summer following an El Ni\ no event, rainfall increases in the YRV and Z500 rises up in the tropics, consistent with the significantly positive correlation between them identified in this study.

We attempt to select optimal factors from simulated 850-hPa and 200-hPa zonal and meridional velocities, sea level pressure, and SST in ENSEMBLES rather than the Z500 only. The PDF prediction skill of YRV summer rainfall is lower than that when the Z500 only is used, because the correlation between the selected optimal factors and observed YRV summer rainfall is less stable. We can obtain higher prediction skill only if the locations and variables of the optimal factors are fixed and identified based on the whole period. However, in such a case, artificial skill is introduced (i.e. not the true cross-validation skill) because information in the target year is used when we identify the optimal factors for the target year. Other schemes to extract information from these fields need to be tested or developed in the future.

There are some potential applications of the new prediction scheme. Although it is developed for PDF prediction of YRV summer rainfall in the present study, the optimization of ensemble members by seeking a statistical relationship between well simulated large-scale circulation fields and rainfall provides an alternative approach to the Bayesian PDF prediction of rainfall in other regions or on grid levels that are simulated poorly by GCMs. In addition, there have been many simple and effective statistical models developed for the seasonal prediction of summer rainfall over China (Liao and Zhao, 1992; Chen, 2000; Li and Zeng, 2008; Fan et al., 2008; Wu et al., 2009). Besides calibration and combing predictions from GCMs, the Bayesian scheme with revised variance modeling of the likelihood function may be suitable for extracting prediction information from multiple statistical models. The new scheme with the two revisions may be suitable for making final hybrid predictions (i.e. combining outputs from statistical and dynamic models) to initiate operational PDF prediction of summer rainfall in China.