An Approach for Improving Short-Term Prediction of Summer Rainfall over North China by Decomposing Interannual and Decadal Variability

North China (NC) is a highly populated region, and anomalous climatic conditions have a substantial impact on agriculture, industry and society. Located within the northern extension of the East Asian summer monsoon, NC exhibits typical mid-latitudinal monsoon characteristics with summer (June-August) mean rainfall accounting for more than 60% of the annual total. Rainfall anomalies often result in severe floods or drought. Hence, predicting summer rainfall over NC is of considerable importance for climate scientists, and continues to be a challenging topic up to the present day (Chen, 1999).

During recent decades, many studies have been devoted to the prediction of rainfall ahead of seasons (Barnston and He, 1996; Chen, 2000; Chu et al., 2008). Scientists have developed a variety of approaches to achieve this, including traditional experience-based statistics, physical-mathematical statistics, numerical dynamic model simulations, and the dynamical-statistical joined method (Zeng et al., 2003). The latter method not only obtains advantages through use of state-of-the-art dynamical models in simulating large-scale atmospheric circulation variability, but also overcomes the disadvantages involved in simulating regional-scale rainfall; thus, it is widely considered the most effective approach (Fuentes and Heimann, 1996; Li and Chen, 1999; Chen et al., 2003; Zeng et al., 2003).

However, owing to its relatively coarse resolution, it is difficult for the dynamical approach to predict station-scale rainfall. Thus, researchers have developed various methods to resolve the model's grid scale to the smaller station-scale, and the process is called "downscaling". Such downscaling methods can be categorized into two types: dynamical and statistical. The former involves nudging a higher-resolved regional model into a global model, and although it has been shown to have the potential to simulate extreme events (Díez et al., 2005; Wang et al., 2011), it consumes a huge amount of computing resource. For this reason, in comparison, the latter (the statistical downscaling approach) is more widely used by the community (von Storch et al., 1993; Murphy, 1999).

Various statistical downscaling techniques have been developed. As far as the source of predictors is concerned, researchers extract valuable information from the outputs of a single dynamical model (Hewitson and Crane, 1996) or a set of dynamical models (Liu and Fan, 2012), even from both the preceding-term observations and the model outputs (Lang, 2011). As far as the prediction equation construction is concerned, an efficient approach called the "interannual increment approach" was proposed by Fan et al. (Fan et al., 2008; Fan and Wang, 2009, 2010). They predicted the year-to-year increment of an object variable (predictand) like surface air temperature or rainfall, rather than the variable itself. This approach highlights the anomalies of the predictand and incorporates the significant quasi-biennial oscillation feature of the East Asian climate, thus exhibiting a great advantage. Recently, (Wang and Fan, 2009) developed another new approach called the "tropics analog" method, in which they considered the primary origination of the predictability of the East Asian summer monsoon from the tropics and made predictions by using the analogue of the current factor field to historical records. The "tropics analog" method has also been used to substantially improve the short-term prediction of East Asian summer rainfall (Wang and Fan, 2009; Fan et al., 2011).

However, the current status of downscaling is not optimistic, and the prediction level is unsatisfactory. In particular, the skill level is unstable, with significant year-to-year or decadal variation. Previous studies have illustrated that one important factor responsible for this instability is the inability to separate interannual and decadal signals when constructing the prediction model (Chen and Zhao, 1998; Shukla, 2007; Chen et al., 2008; Qian and Lu, 2010). It is well known that rainfall over NC exhibits significant interannual and decadal variability. On the interannual scale, it features a quasi-biennial period together with an almost five-year period (Chen, 1999; Dai et al., 2003), while on the decadal scale it experienced a substantial shift around the late 1970s from a wetter period that began in the 1950s to a drier period thereafter (Lu, 2003; Ma, 2007).

The key for improving short-term climate prediction is to identify the link between predictors and predictands at different background timescales. Considering this, (Chen et al., 2008) and (Qian and Lu, 2010) suggested that both kinds of signal from decadal and interannual timescales should be utilized. Following this suggestion, we developed an approach for improving short-term prediction of summer rainfall over NC by decomposing interannual and decadal variabilities. We describe and test this new approach in the present paper, which is organized as follows. Section 2 describes the datasets and methodology used. An analysis of the results is presented in section 3, followed by a summary and discussion of the major findings in section 4.

Three sets of outputs from ensemble hindcast experiments from three global coupled ocean-atmosphere models involved in the DEMETER (Development of a European Multi-model Ensemble System for Seasonal-to-Interannual Prediction) project were used for constructing the factors for downscaling.The three models were developed by the European Centre for Medium-Range Weather Forecasts (ECMWF), the United Kingdom Meteorological Office (UKMO), and the Centre National de Recherches Me'te'orologiques (CNRM) respectively. Each set consisted of one ensemble with nine members, and each was integrated beginning from the initial days in May until the end of August for each individual year during the period 1960-2001. A brief summary of the models is provided in Table 1, and further details can be found in (Palmer et al., 2004).

Different kinds of observational variables including sea surface temperature (SST) and atmospheric circulation variables were used to construct the predictors. The SST data used were from the Hadley Centre Global Sea Ice and Sea Surface Temperature (HadISST) dataset (Rayner et al., 2003). The atmospheric circulation data used were from the ERA-40 Reanalysis (Upplal et al., 2004) dataset. Observed monthly rainfall from over 160 stations within mainland China were used for validation and obtained from the National Meteorological Information Center, China Meteorological Administration. The rainfall dataset spanned the period from January 1951 to the present day.

We decomposed the rainfall over NC into the two distinct timescales through a Fourier function transformation filter. The decadal components of NC rainfall were defined as the variability with period greater or equal to nine years. The potential predictors were filtered too. Then, we chose the predictors for decadal/interannual rainfall components and calibrated the decadal/interannual predictive equations, respectively. Finally, the downscaled decadal and interannual rainfall components were combined to obtain the total predicted rainfall. This method is called the "decomposing timescale downscaling" (DTD) scheme.

When building the DTD model for the decadal component, two basin-scale SST oscillation indices describing the Pacific Decadal Oscillation (PDO) and the Atlantic Multidecadal Oscillation (AMO) were used as predictors, since previous studies have demonstrated their substantial influences (Ma and Fu, 2003; Ma and Shao, 2006; Lu et al., 2006; Li and Bates, 2007; Wang et al., 2009). The statistical downscaling prediction model was built through a linear multivariate fit, together with training and validation processes, for each station. The training period was 1960-94, for which a cross-validation-based regression was used. Further details about the method can be found in (Kim et al., 2004).

When building the DTD model for the interannual component, predictors were selected based on significant connections between observed rainfall with large-scale atmospheric circulation variables in the dynamical model's outputs or observational dataset. To reduce noise, an EOF (Empirical Orthogonal Function) analysis was applied to the predictors and predictands, and the first ten leading EOFs were retained to obtain the principal components. Then, an SVDA (singular value decomposition analysis) was used to obtain the coupled connection between regional rainfall over NC and the large-scale atmospheric circulation pattern. The coupled patterns derived from the SVDA could be expressed mathematically as follows:

Here, x and t represent spatial and temporal grid, respectively. m is the total number of SVD modes, and U_i(x) and R_i(x) are the singular vectors of predictor and predictand, respectively, in the i mode.Then, a downscaling transfer function was used: where P_j(t, x) indicates the downscaled prediction, S_i(t) and K_i(t) represent the time expansion coefficient for predictand and predictors of the i mode, and n is the total number of SVD modes retained. In this study, the leading six coupled modes were retained. Thus, individual statistical models at each station for the interannual component could be built by using the training period and the pathway similar to the above decadal component model.

For comparison, a parallel downscaling scheme was constructed by using unfiltered raw data, which we refer to as the "no-decomposition downscaling" (NDD) scheme. The same predictors as in the above interannual model were used, except for that they were derived directly from the model's raw outputs or observations (unfiltered). In other words, the NDD scheme does not distinguish the interannual and decadal components for both the predictands and the predictors.

An independent-sample validation was carried out for the 7-yr period from 1995 to 2001. This period was not included when training the above models. Several variables were used to quantify the skills of the predictions/hindcasts. The first was the anomaly correlation coefficient (ACC), which was calculated as the correlation between predicted and observed rainfall anomaly time series and reflected the model's ability to yield rainfall anomalies at each station. The second was the spatial correlation coefficient (CC) between the predicted and observed rainfall anomaly pattern. This described the model's ability to predict the spatial pattern of rainfall. The third was the root-mean-square error (RMSE), which was calculated as follows: where y_t is hindcasted from the models or the downscaling, y _ot represents the observation, and N and T represent the amount of stations within the NC region and forecast time length, respectively. The smaller the RMSE, the better the predictions. In addition, a comparison of the historical evolution of NC regionally-averaged rainfall with observations was also employed to assess the skill.

The selection of appropriate predictors is one of the most important steps. Two requirements should be satisfied before a variable is selected. It has to be well simulated by the models, and have a stable correlation with the predictand. In other words, only those variables with a high correlation coefficient with the predictand can be possible candidates (Kang et al., 2007). Whether a particular variable is selected should be based on previous studies and operational experience.

Water vapor transported from the South China Sea (SCS), the Bay of Bengal (BOB) and the subtropical western Pacific plays an important role in the East Asian summer monsoon (Huang and Wu, 1989; Huang et al., 1998). An increase in vapor flux transportation from the SCS and BOB can result in increased summer rainfall over NC (Zhao et al., 2002; Li et al., 2002; Zhou and Yu, 2005). Because vapor transportation achieves a maximum in the lower troposphere, the 850-hPa meridional wind was selected as one potential factor.

Figure 1.  Correlation coefficients between the interannual components of observed NC summer rainfall and 850-hPa meridional wind in the three models and their MME. The interval is 0.1, and only values >0.3(<-0.3) are displayed, as solid (dashed) lines. Shading indicates significance at the 95% confidence level.

Figure 1 displays the simultaneous correlation of observed NC summer rainfall with simulated 850-hPa meridional wind in each individual model, as well as the multi-model ensemble (MME) of all three models. Significant positive correlations can be seen over the BOB, Indochina and the Maritime Continent in nearly all three models and their MME. These confirm that the increased NC rainfall is associated with the anomalous vapor transportation from the BOB and the southern SCS. It also shows an opposite correlation of NC rainfall with South China. Besides, negative correlations can be seen in the equatorial central-western Indian Ocean. Thus, we selected the averaged meridional wind as one predictor, calculated over various domains for the different models. For CNRM, the domain was (25^°S-20^°N, 65^°-110^°E), which covers the central-western Indian Ocean and BOB. For ECMWF and the MME, the domain was (7.5^°-30^°N, 110^°-135^°E) and (10^°S-30^°N, 90^°-130^°E), respectively, which covers the SCS and the East China Sea. And finally, for UKMO, the domain was (15^°S-30^°N, 170^°-120^°W), a sub-region of the tropical central-eastern Pacific Ocean.

Figure 2.  Correlation coefficients of the interannual components of observed NC summer rainfall with the simultaneous (a) 500-hPa geopotential heights and (b) the preceding winter (Dec-Feb) SST. The interval is 0.1. Shading over the oceans indicates significance at the 95% confidence level.

Figure 3.  Schematic diagram of the statistical downscaling approaches used in this study: (a) for the decomposing timescale downscaling (DTD) scheme, and (b) for the no-decomposition downscaling(NDD) scheme. Here, Y, Y₁ and Y₂ represent the raw predictand and interdecadal, interannual component, respectively. a₀, a₁, a₂, a₃ is coefficient of different predictors.

Figure 4.  Spatial distribution of temporal correlation coefficients of station summer rainfall time series of the models' raw predictions with observations. Shading indicates significance at the 95% confidence level.

Figure 5.  Statistical significance of ACCs between observations and the models' hindcast rainfall for the 42-yr period (1960-2001) in the DTD scheme (left column) and the NDD scheme (right column). Smaller (bigger) black circles represent stations with significance exceeding the 95% (99%) confidence level. Blank boxes indicate no significance.

Figure 6.  ACCs of hindcast rainfall with observations during the training period 1960-94. Panels (a-d) represent the interannual component in individual models and the MME, and (e) represents the decadal component in the MME. Light and dark shading represent significance at the 95% and 99% confidence levels, respectively.

The western Pacific subtropical high is another key circulation system influencing the East Asian summer monsoon (e.g., Zhu and Yang, 2003). A previous study on summer rainfall over NC and 500-hPa geopotential height anomalies suggested a high degree of correlation between them (Lau and Weng, 2002; Wang et al., 2004). This is confirmed in Fig. 2a, and a similar correlation could be seen when the simulated 500-hPa heights were used (not shown). Thus, for all the models and the MME, the averaged 500-hPa height over the equatorial region (10^°S-20^°N, 60^°W-120^°E) was used as another predictor.

In addition, SST is an important factor influencing summer rainfall over NC. Generally, in the year following an El Niño event, less rainfall occurs over NC, along with more rainfall over the Yangtze River valley and Northeast China (Zhang and Huang, 1998; Yang and Lau, 2004). Tropical Indian Ocean SST, including the dipole mode and the basin-scale mode, also have a substantial influence (e.g., Li et al., 2008; Xie et al., 2009). Additionally, SST in the SCS in early winter has been found to significantly correlate with rainfall over NC in the following summer (He et al., 2003). Thus, tropical SST may provide potential predictability. Figure 2b shows the correlation of NC summer rainfall with preceding winter SST. A significant positive correlation can be seen in the equatorial western Pacific, whereas opposite correlations are seen in the central-eastern equatorial Pacific and the whole basin of the tropical Indian Ocean. These correlations reflect the lagged impact of El Niño events on East Asian summer climate (Zhang et al., 1997, 1999; Wang et al., 2000; Huang et al., 2000, 2004).

On the decadal timescale, most previous studies have shown NC experienced a relatively wet period from the 1950s and then a dry period since the 1970s (Chen, 1999; Lu, 2003). This decadal variability is highly correlated with the PDO and AMO (Ma, 2007; Wang et al., 2009). The warmer (colder) phases of the PDO are coupled with colder (warmer) SST of the westerly drift region. The sea level pressure anomaly is negative (positive) over Siberia and north of Japan, but positive (negative) over most of southern China, which generates an anomalous northwesterly over NC and reduces water vapor transportation from the oceans. The positive AMO phase leads to a stronger East Asian summer monsoon with enhanced rainfall over NC, along with a delayed withdrawal of the Indian summer monsoon. Hence, the preceding spring AMO and annual PDO indices were chosen as the decadal predictors. The PDO index is determined by projecting SST on the leading EOF of monthly SST on the Pacific north of 20^°N (Mantua and Hare, 2002), and can be downloaded from http://jisao.washington.edu./pdo/PDO.latest. The AMO index is based on Enfield et al. (2001), and can be downloaded from http:// www.esrl.noaa.gov/psd/data/correlation/amon.us.long.data.

All the above chosen predictors are summarized in Table 2. Through a step-by-step training process, a statistical model was built for each station, and the schematic diagram for these models is illustrated in Fig.3. Thus, the interannual rainfall component could be hindcasted by the individual models and the decadal rainfall component could be hindcasted by the observed AMO and PDO indices. The downscaled interannual and decadal rainfall components were combined to yield the total predicted rainfall.

To provide a reference for comparison, the raw prediction in the models was analyzed. Figure 4 shows a comparison of the distribution of ACCs between the observed and the modeled raw prediction. As can be seen, nearly all the models show little skill, with the correlation coefficients ranging from -0.25 to 0.1. At some stations, the modeled raw prediction is even opposite to the observed. Also, the MME exhibits an improvement relative to any individual model, albeit an insignificant one. This illustrates that the models' abilities to predict summer rainfall over NC is weak, and the scheme of the MME is favorable for improving the seasonal prediction.

When the DTD scheme is used, the level of prediction significantly improves. From the left column of Fig. 5, one can see that the ACCs are positive at all stations, and are over 0.4 at most of them, which is significant at the 99% level. The regional-averaged ACC is over 0.5 for all three models and their MME. Comparing with the model's raw prediction, the ACCs increase from -0.15 to 0.57 near the northern rim for CNRM, while for ECMWF and UKMO the ACC achieves a value of 0.46 at the stations near the domain center from -0.25.

To understand the individual contributions from different timescale components, a comparison of ACCs between observations and the hindcasts for the two distinct timescales of the downscaling models is given in Fig.6.For both timescales, all ACCs are positive with the maximum over 0.6, much higher than the models' raw predictions. At some stations near the northern rim, the correlation even reaches 0.75 for the interannual timescale (significant at the 99% confidence level). Meanwhile, near the center, the maximum correlation is close to 0.79 for the decadal timescale. This means that the source of improvement near the northern rim is primarily from the interannual component, while from the decadal component for the improvement in the central region. This can also be seen from a quantitative comparison of the relative contribution rates from the two timescales, which were calculated through the follow equation: Here, ŷ represents the predicted total, while y₁ and y₂ represent the interannual and decadal rainfall, respectively. The relatively more important component at each station is displayed in the left column of Fig. 7. Clearly, the contribution rate of decadal components is higher in the central station, while near the northern part the interannual components prevail.

Figure 7.  A comparison of the two different timescale components (left column) or the various predictors (right column) in the relative contributions to the improvements of the hindcasted rainfall. In the left column, solid black circles (blank boxes) indicate that the decadal (interannual) component has the greater contribution. In the right column, solid black circles, blank boxes and triangles represent that SST, 850-hPa meridional wind and 500-hPa height has the most important contribution, respectively.

Figure 8.  Evolution (solid line marked with crosses) of spatial correlation coefficients of the DTD hindcasts for (a-c) individual models and the (d) MME with the observed rainfall interannual component. For comparison, the evolution (dashed lines marked with boxes) of a similar correlation but from the models' raw predictions is included. The value in the upper-left corner of each panel is the correlation coefficient between the two evolution curves.

As a comparison, the ACCs in the NDD scheme are displayed in the right column of Fig. 5. Relative to the models' raw predictions, the ACCs are higher than all of them, as well as their MME, for some stations (cp. Fig. 4). This suggests that NDD downscaling does indeed improve the level of prediction to some extent. However, the ACCs from the NDD scheme are significantly smaller than in the DTD scheme. In particular, the number of stations with a significant improvement in skill is much less for all the models. This indicates that the DTD scheme can indeed significantly improve the skill relative to both the models' raw predictions and the NDD scheme.

The above results can also be seen from the temporal evolution of the spatial anomaly CCs between the observations and hindcasts. For the DTD scheme (Fig. 8), a considerable consistency can be seen in all three models and their MME. The mean spatial anomaly CC during the training period (1960-94) is 0.58, 0.57, 0.58 and 0.57 for CNRM, ECMWF, UKMO and MME, respectively, while these values are 0.17, 0.21, 0.21 and 0.18 in their raw predictions. We also calculated the significance of the differences between the two sets of CCs, and obtained a t value of 2.81, 3.02, 2.68 and 2.59 for ECMWF, UKMO, CNRM and MME, respectively-all of them being significant at the 95% confidence level. This again suggests the DTD scheme's improvement relative to the models' raw predictions. The modest correlation between the two evolution curves of DTD scheme's hindcast and the raw predictions illustrates that the improvement is not uniform from year to year. For the validation period (1995-2001), the correlation is also higher than the models' raw predictions.

Figure 9.  Differences of RMSE in the models' raw predictions minus that in the hindcast after DTD downscaling.

Figure 10.  Evolution of regional-averaged summer rainfall anomalies over NC during 1960-2001. The short dashed curve indicates the observed, while the solid curve represents the values from the DTD scheme. The dashed (solid) line indicates the trend in the observations (DTD scheme predictions). The gray vertical line indicates the division for the period used for training the downscaling models and for validation as independent samples. The correlation coefficient between the observed and DTD-downscaled rainfall are plotted in the top right of each panel.

When NDD was conducted, the improvement after the statistical downscaling was modest (not shown). The mean spatial CC between the predicted and observed rainfall anomalies during the period 1960-2001 is 0.27, 0.31, 0.3 and 0.32 for CNRM, ECMWF, UKMO and MME, respectively. This suggests that the spatial distribution of rainfall can be better predicted using the NDD scheme than the models' raw predictions, which is in agreement with its good performance in capturing the temporal evolution of rainfall at some stations. However, the efficiency of the NDD scheme is still not as high as the DTD scheme (cp. Fig. 8).

As mentioned above, the bias of prediction or hindcast relative to the observations can be measured by the RMSE. Figure 9 presents the difference in RMSE between the models' raw predictions and the DTD scheme's downscaling prediction. As can be seen, all values are positive for the three models and their MME, illustrating that the DTD scheme improves the level of prediction consistently. The maximum of the differences is located near the northern rim, indicating the largest significance there, in agreement with that shown in Fig. 6.

In addition, a further evaluation was conducted by comparing the time series of the observed and predicted regional mean rainfall anomalies (Fig. 10), which was defined as the mean of 21 stations within the domain (Fig. 5a). In the observations, a negative (drying) trend can be seen during the whole period. A similar negative trend exists in the DTD hindcasts, particularly for UKMO and CNRM. This consistency in trend is in agreement with their substantial correlation with a coefficient of around 0.7, which is significant at the 99% confidence level. For the independent sample validation period (1995-2001), the ratio with the same sign of anomaly in the predicted rainfall with the observed is 100%. This result is higher than that reported by (Fan et al., 2008), which was 40% (three instances when the predictions were consistent with observations within eight years) achieved by adapting the year-to-year increment technique to construct a statistical downscaling scheme without timescale decomposition (cp. their Fig.6 with Fig.10 here). This indicates a considerable advantage when timescale decomposition is used.

It is extremely valuable to examine which component of predictors used plays the primary role in improving the ACCs in the DTD scheme. Similar to the above calculation regarding the contribution of the two distinct timescale components, we calculated the relative contribution rates of the various predictors to the hindcast interannual rainfall component. The right column of Fig. 7 displays the predictor with the greatest contribution over the individual stations. From this, it can be seen that SST is the most important factor for all the models and their MME, since at nearly half of the stations it explains the greatest contribution rate. This suggests that the decomposition of SST may be the most important factor for DTD downscaling.

In the present reported work, a statistical downscaling technique was developed to predict the summer rainfall over NC. First, both the predictors and predictands were decomposed into interannual and decadal timescale components. Then, two individual models were built for the two timescales using multivariate linear regression. Finally, the predicted rainfall for the two timescales was combined to represent the final prediction. Validations of the predictions based on independent samples suggested a substantial improvement in skill in station-based temporal correlation, ACCs and RMSE, relative to the conventional without-timescale decomposition and the models' raw predictions. Moreover, the present approach reproduces the observed linear trend of rainfall over NC.

The key to the success of this approach is the selection of different potential predictors for the distinct timescale predictands. In other words, it considers the different timescale predictands are influenced by different physical processes. As in conventional statistical downscaling, potential predictors are selected based on the correlations between the observed rainfall component and the large-scale circulation from the model's simulations or observations, and also the correlations with the lower boundary conditions, such as SST. Thus, it conveys information on large-scale atmospheric circulation to local regional rainfall, considering the dynamical model's advantage in simulating large-scale circulation and avoiding the model's disadvantage in simulating regional or local rainfall. For the interannual timescale, three large-scale variables (850-hPa meriondinal wind, 500-hPa geopotential heights and tropical SST) were selected as the potential predictors. For the decadal timescale, the PDO and the AMO were calibrated to the prediction equation.

Previous studies have speculated that one potential approach to improving short-term climate predictions is to incorporate decadal signals (e.g., Chen et al., 2008; Qian and Lu, 2010). This study set out to investigate that hypothesis, and the results showed significant improvements in comparison with traditional statistical downscaling approaches without timescale decomposition. Nonetheless, more validations using state-of-the-art techniques are expected in the future.