A Timescale Decomposed Threshold Regression Downscaling Approach to Forecasting South China Early Summer Rainfall

An important task for flood and drought management is the provision of accurate rainfall prediction. The rainfall distribution in China is generally characterized by a "southern flood and northern drought" pattern, due to the weakening of the East Asian summer monsoon after the late 1970s (Nitta and Hu, 1996; Wang, 2001; Gong and Ho, 2002). Accompanied by the northward seasonal march of the East Asia summer monsoon, abundant rainfall first appears over South China (SC) in mid-May (Tao, 1987; Lau and Weng, 2001; Ding et al., 2008; Wu et al., 2012), and leads to the peak of annual rainfall in early summer (June) over the region (Fig. 1a). Because SC is one of China's largest economic zones with a large population, extreme flood events in early summer often cause a large number of human casualties and

considerable economic loss (Chan and Zhou, 2005; Zhou et al., 2006). Thus, improving seasonal forecasting skill for early summer rainfall over SC is of great importance for disaster prevention and mitigation.

Figure 1.  (a) The annual climatological cycle of the rainfall at the stations of Guangzhou, Hong Kong and Macao, as well as the 3-station average, during 1910-2011. (b) Detrended correlation of the 3-station averaged rainfall with the 160-station rainfall from CMA in June during 1951-2011. The stations at Guangzhou, Hong Kong and Macao are indicated by the three red dots. The grey shading indicates CCs significant at the 99% confidence level based on the t-test, presenting the SC region. (c) Power spectrum for June rainfall. The peak over the red dashed line indicates the confidence level is greater than 80% against red noise.

There are several approaches to forecasting rainfall. One is using numerical models (e.g. Collischonn et al., 2005; Aligo et al., 2009). The raw model prediction is dynamically meaningful, but still has low skill for many reasons; for instance, because of the model resolution (Martin, 1999), subgrid processes (Grotch and MacCracken, 1991) and parameterization schemes (Eitzen and Randall, 1999). Hence, precipitation modeling is regarded as one of the most difficult issues in climate modeling, resulting in low skill on the basis of the pure numerical approach. Another approach to forecasting rainfall is based on statistical downscaling models (Liu and Fan, 2012a, 2014; Sun and Chen, 2012). A statistical downscaling method combining the GCM-simulated and observed information is developed, which shows much better predictability for global precipitation forecasting (Sun and Chen, 2012). In addition, pure statistical downscaling models are also used to forecast rainfall. The method is generally based on an empirical observed relationship between the large-scale climate anomalies and local rainfall fluctuations. There are various methods that can be used to develop statistical downscaling models, including multiple linear regression (Wilby, 1998), principle components (Li and Smith, 2009), and singular value decomposition (Zhu et al., 2008; Liu and Fan, 2012b). Other more sophisticated methods include partial least squares (PLS) regression (Zhang et al., 2011; Wu et al., 2013).

Pure statistical downscaling models show reasonable skill in predicting regional rainfall. For example, (Sahai et al., 2003) made optimum use of global SST for Indian summer monsoon rainfall prediction nine months in advance. Rainfall may contain variabilities on various timescales, with low-frequency interdecadal variability associated with pronounced wetting or drying trends (Ding et al., 2008) and high-frequency interannual variability related to severe floods or droughts (Huang et al., 2006). (Guo et al., 2012) described a timescale decomposition (TSD) approach to statistically downscale the late summer rainfall over North China, which made use of two distinct downscaling models corresponding to the interannual and interdecadal rainfall variability, respectively. However, the TSD approach by (Guo et al., 2012) neglected the possibility that the interannual rainfall variability may be modulated by an interdecadal climatic background.

Recently, an increasing number of studies have suggested that the interannual relationships of large-scale climate anomalies with local or remote climate fluctuations are not stationary (e.g. Torrence and Webster, 1999; Wang, 2002; Wu and Wang, 2002; Gao et al., 2006; Sun and Wang, 2012; Chen et al., 2013a; Chen et al., 2015a, 2015b; Cao et al., 2015, 2016). In particular, the modulation of a decadal-scale coupled ocean atmospheric mode named the Pacific Decadal Oscillation (PDO) (Mantua et al., 1997; Mantua and Hare, 2002) has attracted more attention in the last decade (Gershunov and Barnett, 1998; Power et al., 1999; Chan and Zhou, 2005; Wang et al., 2008; Mao et al., 2011; Chen et al., 2013b; Duan et al., 2013). For example, it has been found that the typical influences of ENSO on the North American climate are strong and consistent only during preferred phases of the PDO (Gershunov and Barnett, 1998). The ENSO-East Asian summer monsoon relationship during 1962-77 has been found to be significantly different from that during 1978-93 (Wu and Wang, 2002), which is actually consistent with a phase transition from negative to positive PDO in the late-1970s. In the SC region, the interannual rainfall variation has also been suggested to be modulated by the PDO (Chan and Zhou, 2005; Mao et al., 2011; Duan et al., 2013). It has been proposed that incorporating information on the PDO could improve the long-term predictability of early summer SC rainfall changes (Chan and Zhou, 2005; Zhou et al., 2006). Therefore, a statistical model that takes into account the modulation of the PDO may lead to better skill in seasonal rainfall prediction over SC. This naturally raises two key questions: How can the modulation effect of the decadal-scale coupled oceanic-atmospheric mode of the PDO be incorporated when developing a statistical downscaling model to forecast the early summer rainfall over SC? And can the forecasting skill be improved?

The aim of this paper is to describe a TSD threshold regression (TSDTR) downscaling approach for the statistical forecasting of South China early summer rainfall (SCESR). The TSDTR downscaling approach includes two distinct regression downscaling models, called the interannual model (IAM) and interdecadal model (IDM), respectively linking SST patterns prior to early summer to the interannual and interdecadal variability of SCESR. On the interannual timescale, the IAM, based on a threshold PLS regression model, is fitted to interannual components of SCESR and March SST patterns by taking account of the modulation of negative and positive phases of the PDO. On the interdecadal timescale, the IDM, based on a PLS regression model, is employed to fit the relationship between SCESR and preceding November SST patterns. The total rainfall prediction is obtained by the sum of the outputs from both the IAM and IDM.

The rest of this paper is arranged as follows: Section 2 introduces the datasets used in this work. The TSDTR downscaling approach is described in section 3. Results of the TSDTR approach in forecasting SCESR are presented in section 4. Section 5 is a summary and discussion.

Long-term reliable observed rainfall data are derived from the monthly rain gauge datasets from three stations: Hong Kong (22.11^°N, 114.14^°E), Macao (22.16^°N, 113.35^°E) and Guangzhou (23.08^°N, 113.16^°E). These rainfall data have been used by previous studies (Chan and Zhou, 2005; Mao et al., 2011; Duan et al., 2013) and are chosen here because of their long-term reliability since 1910. The maximum rainfall of these three stations occurs in June (Fig. 1a), which purely belongs to the South China Sea summer monsoon rainfall (Tao, 1987; Wang et al., 2004). Another set of monthly rainfall data derived from 160 Chinese meteorological stations provided by the China Meteorological Administration (CMA) is also used in this study, available from 1951. The simultaneous correlation pattern of three-station averaged rainfall in June with the 160-station rainfall from CMA during the period 1951-2011 shows significant positive correlation over the region of SC (Fig. 1b). Thus, it is reasonable to use the long-term, high-quality three-station (Hong Kong, Macao and Guangzhou) averaged rainfall in June to represent the SCESR in this study.

The PDO index, defined as the leading EOF of SST anomalies in the North Pacific Ocean poleward of 20^°N after removing the global warming signal (Mantua et al., 1997; Mantua and Hare, 2002), is taken from the website of the Joint Institute for the Study of the Atmosphere and Ocean (http://jisao.washington.edu/pdo/PDO.latest), which contains data from 1900 onward. An 11-year running mean of wintertime (October-March averaged) PDO index is used to obtain the interdecadal variability of the PDO (Mantua and Hare, 2002; Mao et al., 2011) over the period 1910-2011 (Fig. 2a). Note that 10 years are missing, caused by the 11-year running mean, and thus the period for the forecasting experiment in this work is chosen as 1915-2006. Monthly SST data are extracted from NOAA ERSST.v3b (Smith and Reynolds, 2004; Smith et al., 2008) on 2.0^°× 2.0^° grid, which is available from the year 1854 at http://www.esrl.noaa.gov/psd/data/gridded/.Monthly atmospheric data are obtained for the period from 1948 from the NCEP-NCAR reanalysis products on 2.5^°× 2.5^° grid (Kalnay et al., 1996), including horizontal winds and specific humidity at different levels (also available at http://www.esrl.noaa.gov/psd/data/gridded/).

Figure 2.  (a) The standardized time series of wintertime (October-March) PDO index and their 11-year running mean during 1915-2006. As an example, the year 1915 denotes October-December 1914 and January-March 1915. (b) Total rainfall (Rain$_\rm T$, units: mm) in June during 1915-2006, which is decomposed into an interannual component (Rain$_\rm A$) with variation less than 11 years and an interdecadal component (Rain$_\rm D$) with variation longer than 11 years.

Two primary peaks with periods of about 4 years and 17 years exist in the SCESR, as determined by spectrum analysis (Fig. 1c), indicating apparent interannual and interdecadal variability. Therefore, the observed time series of total rainfall (Rain _T) is decomposed into interannual (variation less than 11 years, denoted as Rain _A) and interdecadal (variation longer than 11 years, denoted as Rain _D) components by an 11-year high-pass and low-pass filter (Fig. 2b). That is, Rain _T= Rain _A+ Rain _D.

SST anomalies have been found to be important in influencing the East Asian climate (Huang and Wu, 1989; Wang et al., 2000; Lau and Weng, 2001; Wu et al., 2012), have a relatively long "memory" (Sahai et al., 2003), and can therefore be regarded as a preceding predictor for SCESR. This motivates us to explore the relationships between SCESR and associated SST patterns on both interannual and interdecadal timescales. As such, we also decompose the SST field into interannual (SST _A) and interdecadal (SST _D) components. To reveal the dominant SST patterns associated with SCESR variations, we employ the PLS regression method. Specifically, PLS embodies the well-known concept of partial correlation, as it seeks the predictors Z, which are linear combinations of the factors X, being referred to as latent vectors or PLS components, and maximizes the variance explained in Y and the correlation between X and Y (Haenlein and Kaplan, 2004; Smoliak et al., 2010). Unlike EOF analysis, which identifies major patterns explaining SST variations, using PLS regression we can find the PLS components of SST variations that best explain the covariance between SST variations and the SCESR variations. In other words, it reveals the dominant SST patterns that not only account for most of the SST variations but are also closely related to the SCESR variability.

Figure 3.  Schematic plot of the TSDTR approach to forecasting the SCESR.

Figure 3 shows a schematic plot of the TSDTR approach to forecasting the SCESR based on the PLS regression method. On the interannual timescale, in order to incorporate the modulation effect of the PDO on the interannual variability of SCESR (Mao et al., 2011; Duan et al., 2013), a threshold PLS regression model is calibrated for the relationship between Rain _A and SST _A under the positive and negative phase of the PDO. On the interdecadal timescale, a standard PLS regression model is calibrated for the relationship between Rain _D and associated SST _D patterns. The total rainfall prediction is obtained by the sum of the outputs \(\hat{R}ain_\rm A\)(t) and \(\hat{R}ain_\rm D\)(t) from both the IAM and IDM (Fig. 3). To test the performance of the TSDTR approach in forecasting the SCESR, the study period 1915-2006 (N=92) is separated into a calibration period [1915-84 (n=70)] and validation period (1985-2006). The two PLS-based regression downscaling models in the TSDTR approach are calibrated by using the calibration data in 1915-84, and the forecasting skill of the TSDTR is tested by the independent validation data in 1985-2006. To obtain the forecasted values in the independent validation period, we use the running forecasting method based on the calibrated models. That is, for t=n+1,n+2,…,N (=92), when the observed preceding SST(t) data and wintertime PDO are available over the validation period 1985-2006, we add these new SST data to those in the training period, and then decompose the combined SST data into interannual SST _A(t) and interdecadal SST _D(t) components. The forecasted value of interannual [interdecadal] rainfall can be estimated by using SST _A(t) [SST _D(t)] and the calibrated IAM [IDM]. The performance of the TSDTR downscaling approach is assessed through the correlation coefficient (CC) between predicted and observed values and the RMSE (Zhang et al., 2011; Guo et al., 2012). The uncertainty of the forecast is indicated by the spread of bootstrapping prediction intervals [see appendix in (Li and Smith, 2009)].

This section provides details of the TSDTR downscaling approach to building the IAM and IDM based on the PLS regression for the relationship between the preceding month SST and the SCESR on the interannual and interdecadal timescales, respectively.

The interannual variation of SCESR is remarkably different under different PDO phases. (Chan and Zhou, 2005) demonstrated that the interannual relationship between ENSO and SCESR is modulated by the phase of the PDO. (Mao et al., 2011) indicated that the interannual SCESR variations are remarkably different under different PDO phases, based on a comparison in two typical epochs of 1958-76 and 1980-98. Furthermore, it was suggested by (Duan et al., 2013) that the predictability of interannual SCESR variability is modulated by different PDO-phase backgrounds, and they also indicated that the relationship between interannual rainfall and preceding Pacific SST anomalies experienced a robust interdecadal change due to the PDO's modulation through the so-called "seasonal footprinting mechanism".

Figure 4.  Lead-lag correlation of Rain$_\rm A$ with SST$_\rm A$ in the Pacific and Indian oceans in the (a-c) negative phase and (d-f) positive phase of the PDO. The contour lines present the CCs of $\pm0.3$, statistically significant at the 95% confidence level.

Figure 5.  First two leading modes of March SST$_\rm A$ in the Pacific and Indian oceans from PLS regression analyses in the (a, b) negative phase and (c, d) positive phase of the PDO. The first percentage value is the total variance of Rain$_\rm A$ explained by the SST$_\rm A$ mode, and the second value is the total SST$_\rm A$ variance explained by the same mode.

To further explore the predictability of interannual SCESR variability modulated by different PDO-phase backgrounds, Fig. 4 shows the lag-correlation patterns of interannual SCESR with preceding interannual monthly (May, April and March) SST _A over the Pacific and Indian oceans in negative and positive phases of the PDO, denoted by PDO(-) and PDO(+), respectively. It is evident that the interannual relationship between Rain _A and SST _A exhibits very different structures during the positive and negative PDO phases. A striking difference in Fig. 4 is that there is a traditional eastern Pacific warming ENSO-like correlation pattern between Rain _A and SST _A during the PDO(+) phase, but a central Pacific warming ENSO-like pattern is more pronounced during the PDO(-) phase. Such correlation patterns may lead the SCESR by up to 3 months from March to May. This result is consistent with the result of (Duan et al., 2013), who compared 1955-76 PDO(-) and 1977-98 PDO(+). The result supports the fact that the interannual SCESR variability and its relationship with ENSO are modulated by the PDO phases (Chan and Zhou, 2005; Zhou et al., 2006; Mao et al., 2011). Therefore, the interannual relationship between large-scale SST anomalies and the SCESR is not stationary due to the modulation of the PDO. As such, in order to incorporate the modulation effect of the PDO on the interannual variability of SCESR, we need to consider the nonstationary relationship between Rain _A and SST _A in the Pacific and Indian oceans. To this end, a threshold PLS regression model (Fig. 3) is employed to establish the relationship between Rain _A and the leading SST _A modes, which best explain the co-variations of Rain _A and SST _A under the PDO(-) and PDO(+), separately. After a series of tests using SST _A predictors in each preceding month, including the preceding May, April, March, February, January, December and November, it is found that SST _A conditions in the preceding March tend to yield the best results. Thus, we concentrate on reporting the dominant SST _A modes in March over the Pacific and Indian oceans that influence Rain _A.

Figure 5 shows the first two leading interannual modes of March SST (SST _A-I and SST _A-II) represented by PLS loadings in PDO(-) and PDO(+), respectively. The significance of each mode is reflected by two numbers: the first is the percentage of Rain _A variance explained by the SST-PLS mode; the other is the percentage of SST _A variance explained by the same SST-PLS mode. In PDO(-), the first two leading modes together explain 73.2% of the total variance of Rain _A. The first dominant mode (Fig. 5a) during the negative phase of the PDO is, by and large, an El Niño Modoki-like pattern in the tropical Pacific region, with positive loadings in the central Pacific and negative loadings in the western and eastern Pacific. Another significant positive loading occurs over the South-central Pacific region around New Zealand. The tropical Indian Ocean is characterized by weak positive loadings. As shown in Fig. 6a with respect to the corresponding patterns of precipitation and vertically integrated moisture flux to the SST _A-I mode, notable southwesterly moisture transportation to the SC region is observed in the troposphere, which is accompanied by a significant low-level anomalous anticyclone located in the South China Sea and Philippine Sea (not shown). At the same time, southward penetration of northerly moisture transportation from a large part of North and Northeast China also exists and leads to significant moisture convergence over the SC region. Therefore, abundant interannual rainfall is received in this region (Fig. 6a).

The second dominant SST mode (Fig. 5b) in PDO(-) has strong negative loadings in the central-eastern tropical Pacific and weaker negative loadings in the tropical Indian Ocean, indicating a La Niña-like pattern. Signals can also be found in the midlatitudes, with anomalous warming in the North Pacific region and South-central Pacific region around New Zealand. The pattern in the tropical and midlatitude North Pacific resembles the anomalous SST features of negative PDO phase. Correspondingly, one can see notable southerly moisture transportation over the majority of eastern China (Fig. 6b), which is accompanied by an anomalous low-level anticyclone located in the north to northeast of the Philippine Sea (not shown). The SC region is under the control of southerly moisture transport, but the moist flows are strong enough to advance more northward. Thus, distribution of anomalous precipitation in China shows weak above-normal interannual rainfall over SC (Fig. 6b).

Figure 6.  Correlations of the detrended scores of March SST$_\rm A$ modes in Fig. 5 with high-pass filtered June precipitation (color shading) and vertically integrated moisture flux between 1000 hPa and 300 hPa (vectors) in the (a, b) negative phase and (c, d) positive phase of the PDO.

Figure 7.  (a) PLS regression model results for forecasting Rain$_\rm A$. PDO negative (positive) years are marked with "o" ("*"). The black line is the time series of observed Rain$_\rm A$; the red line is the reconstructed Rain$_\rm A$ during 1915-84; and the blue line is the hindcasted Rain$_\rm A$ during 1985-2006. (b) PLS regression model results for forecasting Rain$_\rm D$. The black line is the time series of observed Rain$_\rm D$; the red line is the reconstructed Rain$_\rm D$ during 1915-84; and the blue line is the hindcasted Rain$_\rm D$ during 1985-2006. (c) Model results for forecasting the total rainfall. The black line is the time series of observed Rain$_\rm T$; the red line is the reconstructed Rain$_\rm T$ during 1915-84; and the blue line is the hindcasted Rain$_\rm T$ during 1985-2006. The blue shading in each of the panels represents the upper and lower bands of the bootstrapping 95% confidence intervals for the hindcasted values.

In PDO(+), the first two leading modes together explain 58.3% of the total variance of interannual rainfall. The leading interannual SST modes (SST _A-I and SST _A-II) dramatically change in contrast with those in PDO(-). The first dominant mode (Fig. 5c) tends to suggest a strengthened traditional El Niño-like signal in influencing the interannual early summer rainfall over SC during the positive phase of the PDO, with very strong positive loadings in the tropical eastern Pacific. Another positive loading occurs in the tropical Indian Ocean. In addition, the positive loading in the South-central Pacific region around New Zealand in PDO(-) is replaced by negative loading. Meanwhile, the pattern in the midlatitude North Pacific shares similar features with SST anomalies during the positive phase of the PDO. The circulation responses to this mode (Fig. 6c) show notable southwesterly moisture transportation from the South China Sea to the southeast coast of China, which leads to above-normal precipitation over the SC region. The second dominant interannual SST mode (SST _A-II) in PDO(+) has negative loadings in the central-eastern tropical Pacific that display a La Niña-like pattern (Fig. 5d), but has reduced in this mode compared with that in PDO(-) (Fig. 5b). The Indian Ocean is occupied by weak negative loadings. In addition, a warming anomaly occurs in Peru's inshore waters. Figure 6d shows the corresponding interannual circulation responses. The vertically integrated moisture flux pattern indicates that the main moisture source region for interannual rainfall over SC in this mode is from the Bay of Bengal. Through the westerly flows carrying warm and wet air, the SC region tends to receive above-normal interannual rainfall.

Based on the above analysis, the different leading interannual SST modes related to the interannual early summer rainfall over SC in PDO(-) and PDO(+) may indicate that different PDO "backgrounds" modulate the connection between the interannual early summer SC rainfall and SST anomalies over the Pacific and Indian oceans. Therefore, when building a statistical forecasting model for interannual rainfall with preceding SST conditions, one needs to consider such nonlinearity. As previously mentioned, a threshold PLS regression model (Fig. 3) is employed to establish the relationship between Rain _A and the leading SST _A modes by taking into account the modulation of PDO(-) and PDO(+).

Figure 7a compares the interannual variation of predicted and observed Rain _A values from the threshold PLS regression model (i.e. the IAM in Fig. 3) based on the first two leading interannual SST modes, which is derived using observed data over the training period of 1915-84 (n=70). The results indicate that Rain _A can be reconstructed with significant skill for the training period (Fig. 7a), with a CC₀ of 0.87 and RMSE₀ of 54.45 mm between the predicted and observed Rain _A values (Model IV in Table 1). Here, CC₀ and RMSE₀ are the CC and RMSE between the observed rainfall and reconstructed rainfall in the calibration period. To obtain the forecasted values [\(\hat{R}ain_\rm A\)(t)] in the independent validation period of 1985-2006, we use the running forecasting method based on the calibrated IAM. The forecasted value of interannual rainfall [\(\hat{R}ain_\rm A\)(t)] can be estimated by using SST _A(t) and the calibrated IAM (i.e. the calibrated threshold PLS regression model based on the preceding wintertime PDO phases). Figure 7a shows the forecasted values [\(\hat{R}ain_\rm A\)(t)] in the validation period of 1985-2006 and their uncertainty in terms of bootstrapping 95% confidence intervals.

Figure 8.  First two leading modes of November SST$_\rm D$ in the western-central North Pacific from PLS regression analyses. The first percentage value is the total variance of Rain$_\rm D$ explained by the SST$_\rm D$ mode and the second number is the total SST$_\rm D$ variance explained by the same mode.

Next we develop an IDM for the interdecadal variability of the SCESR. Previous studies have noted that temporal variations of the SCESR exhibit an interdecadal oscillation related to the PDO, with more dry (wet) years during periods of positive (negative) PDO index (Chan and Zhou, 2005; Duan et al., 2013). (Chan and Zhou, 2005) also mentioned that the effect of the PDO is more important than that of ENSO in the control of SCESR, although such a conclusion might be premature given our limited understanding of the PDO and ENSO, as well as their possible interaction. Because the rainfall and PDO index have a close negative relationship on the interdecadal timescale, and the most visible climatic fingerprints of the PDO exist in the North Pacific (Mantua et al., 1997; Mantua and Hare, 2002), it is supposed that the most important interdecadal SST signal influencing the interdecadal rainfall variability may exist in the midlatitude North Pacific region. With the aim to develop a PLS-based regression downscaling model with preceding interdecadal SST conditions to forecast the interdecadal rainfall variability, a series of tests using SST _A predictors in each preceding month, including the preceding May, April, March, February, January, December and November, are conducted, and the results show that interdecadal SST conditions over the midlatitude western-central North Pacific [(20^°-50^°N, 120^°E-140^°W)] in the preceding November give the best results. Thus, we focus on reporting the leading interdecadal SST modes in November over the western-central Pacific region and the corresponding atmospheric circulation responses to the SST modes influencing the interdecadal SC rainfall variability. To this end, the PLS regression model (Fig. 3) is employed to establish the relationship between Rain _D and the leading SST _D modes using the training data in the period 1915-84.

Figure 8 shows the first two leading interdecadal SST modes (SST _D-I and SST _D-II) represented by PLS loadings. These two leading modes together explain 74.1% of the variance of Rain _D, with the majority contributed from the first leading mode (58.1%). The SST _D-I mode (Fig. 8a) is characterized by significant positive loadings that almost cover the entire western-central North Pacific, with two maximum anomalous SST centers located in the Yellow Sea and midlatitude central North Pacific located at about (35^°N, 160^°W). However, these positive loadings in the western-central North Pacific shrink in the second mode (SST _D-II); instead, negative loadings widely expand (Fig. 8b). The SST _D-I mode resembles the characteristics of the PDO negative phase to some extent, with a typical warm SST anomaly pattern in the central North Pacific (Mantua et al., 1997; Mantua and Hare, 2002). This is therefore consistent with previous results in which persistent wet conditions over the SC are related to periods of negative PDO index (Chan and Zhou, 2005; Zhou et al., 2006). In the second dominant SST mode (SST _D-II), the shrinking of positive loadings may imply a phase transformation from negative to positive PDO.

To examine the corresponding atmospheric circulation anomalies associated with the above two dominant SST modes on the interdecadal timescale, we show the correlation patterns of the scores of SST modes with the vertically integrated moisture flux between 1000 hPa and 300 hPa (Fig. 9). The anomalous precipitation distributions in China are also displayed in Fig. 9. For the SST _D-I mode, remarkable southwesterly moisture flux enters into the SC region, which would lead to persistently wet conditions in the region (Fig. 9a). As a result, a significant positive correlation with interdecadal rainfall variability over the SC region ultimately occurs (Fig. 9a). The atmospheric circulation responses to the second mode (SST _D-II) are basically similar to those for the first mode, only with small differences in magnitude, thus also exerting a beneficial influence on above-normal interdecadal rainfall variability over the SC region (Fig. 9b). The SST anomalies in the North Pacific associated with the PDO can persist for 20-30 years, from winter to summer (Mantua et al., 1997), and thus seem to be important in the control of the persistently wet or dry conditions in early summer over SC. Note that since the anomalies associated with persistently wet and dry conditions generally tend to have opposite polarities, the reverse is true for interdecadal below-normal rainfall variability.

Figure 7b compares the interdecadal variation of predicted and observed Rain _D values from the PLS regression model (i.e. the IDM in Fig. 3) based on the first two leading modes (SST _D-I and SST _D-II), which is derived using observed data over the training period of 1915-84 (n=70). It is evident that the IDM can reproduce the time series of interdecadal rainfall variability in the training period well, with a CC₀ of 0.88 and RMSE₀ of 62.67 mm between the predicted and observed values of Rain _D. To obtain the forecasted interdecadal values [\(\hat{R}ain_\rm D\)(t)] in the independent validation period of 1985-2006, we use the running forecasting method based on the calibrated IDM. The forecasted value of interdecadal rainfall [\(\hat{R}ain_\rm D\)(t)] can be estimated by using SST _D(t) and the calibrated IDM. Figure 7b shows the forecasted values [\(\hat{R}ain_\rm D\)(t)] in the validation period of 1985-2006 and their uncertainty in terms of bootstrapping 95% confidence intervals.

It is straightforward to forecast the values of the SCESR (i.e. Rain _T) by summing up the forecasted values \(\hat{R} ain_\rm A\)(t) and \(\hat{R}\rm ain_\rm D \)(t) from the IAM and IDM models (i.e. the TSDTR downscaling approach). Figure 7c shows the performance of forecasting the SCESR. In general, the performance in the training period is maintained in the subsequent validation period. For example, compared to the observed climatological rainfall of 343.80 mm during the training period of 1915-84 and 367.70 mm during the validation period of 1985-2006, the statistical model provides a preproduction of 343.80 mm and 351.15 mm, respectively. The TSDTR downscaling approach provides reasonable forecasting skill by using preceding SST as the only predictor. Table 1 shows the forecasting skill of total rainfall by showing the associated CC ( CC₀=0.88) and RMSE ( RMSE₀=62.67 mm) between the downscaled and observed rainfall for the training period of 1915-84. This skill is maintained reasonably well in the validation period with a CC₁ of 0.56 and RMSE₁ of 139.47 mm (Table 1, Model IV). Here, CC₁ and RMSE₁ are the CC and RMSE between the observed rainfall and predicted rainfall in the validation period. The bootstrapping 95% confidence intervals associated with the forecasted values are shown in Fig. 7c, indicating the uncertainty of the TSDTR downscaling approach. All these results indicate that the TSDTR downscaling approach can achieve some reasonable rainfall forecasting skill over SC.

Figure 9.  Correlations of the detrended scores of November SST$_\rm D$ modes in Fig. 8 with low-pass filtered June precipitation (color shaded) and vertically integrated moisture flux between 1000 hPa and 300 hPa (vectors).

To see how well the TSDTR downscaling approach performs in forecasting the SCESR, compared to other simple methods, we examine three different methods for forecasting the SCESR. We analyze the use of the previous 5 years' SCESR average (P5YSA), a single PLS model, and a TSD-PLS model. The P5YSA is a computationally convenient and frequently used technique to forecast a regular time series, which is carried out without recourse to a formal statistical model. The single PLS model is calibrated by directly regressing the SCESR onto the SST field, while the TSD-PLS model is calibrated by two steps; namely, decomposing the rainfall and SST field into interannual and interdecadal components, and then applying the PLS method to them. All these three methods do not integrate the modulation of the PDO on the interannual relationship between SCESR and SST. Note that the first two dominant SST patterns are used in both the PLS and TSD-PLS models for a fair comparison with the TSDTR downscaling approach. Table 1 summarizes how the forecasting skill of the TSDTR downscaling approach compares with the three simpler models. The CC between the forecasted and observed SCESR values in the validation period for P5YSA, PLS and TSD-PLS are CC₁=0.13, 0.06 and 0.46, respectively, which is smaller than that (0.56) from the TSDTR downscaling approach. In addition, the RMSE between the forecasted and observed SCESR values in the validation period for P5YSA, PLS and TSD-PLS are RMSE₁=167.22, 183.81 and 148.06 mm, respectively——larger than the RMSE ( RMSE₁=139.47 mm) for the TSDTR downscaling approach. Therefore, the more complex TSDTR downscaling approach performs better than the simpler P5YSA, PLS and TSD-PLS models in forecasting the SCESR. We thus emphasize that the application of the newly proposed TSDTR method can improve regional rainfall prediction skill by incorporating the PDO modulation effect, even when the skill of the pure PLS model is low.

In summary, the rationale behind the TSDTR downscaling approach is that it allows a model parameter change for the unstable relationship between the interannual rainfall variability and preceding SST conditions in different PDO phases. As a result, the TSDTR downscaling approach can further improve rainfall forecasting skill by considering the modulation effect of decadal-scale coupled oceanic-atmospheric modes on interannual climate variability, compared to other simpler methods (e.g. P5YSA, PLS and TSD-PLS), without considering the preceding SST conditions in different PDO phases.

This paper puts forward the TSDTR downscaling approach to forecast SCESR by using preceding SST over the Indian and Pacific oceans through modeling the interannual and interdecadal rainfall variability via an IAM and IDM. The IAM is built upon a threshold PLS regression by taking account of the modulation effect of the PDO on the interannual variability of the SCESR. Through the IAM, the first two leading March SST modes (an El Niño Modoki-like pattern and a La Niña-like pattern) are linked to the interannual rainfall in the negative phase of the PDO, while the first two leading modes (a strengthened traditional El Niño-like pattern and a La Niña-like pattern) are linked to the interannual variability of the SCESR in the positive phase of the PDO. The first two leading interdecadal SST modes are linked to the interdecadal component of the SCESR via the IDM based on the PLS regression method. The forecasted total rainfall is obtained by summing up the values of the two forecasted components from the IAM and IDM.

In applying the TSDTR approach to forecasting SCESR, it is found that the interannual relationship between the preceding Pacific-Indian Ocean SST and SCESR experiences interdecadal changes, displaying remarkable differences in the negative and positive phases of the PDO. In particular, the traditional ENSO-like anomalous SST pattern is robust only in the PDO positive phase. For interdecadal rainfall variability, the dominant interdecadal SST patterns over the North Pacific resemble some characteristics of the PDO. When the North Pacific SST is persistently warming (cooling), the interdecadal variability of SCESR tends to be below-normal (above-normal). This is consistent with previous studies in which persistently wet (dry) conditions over SC in early summer are related to periods of negative (positive) PDO index (Chan and Zhou, 2005; Duan et al., 2013).

A central implication of this study is that regional rainfall forecasting skill can be improved via the TSDTR approach, wherein the unstable relationship between interannual rainfall and large-scale variability, such as SST, can be taken into the process of model development. For a statistical downscaling model, it is important to select the predictors that have the most stable relationships with observed rainfall (Sun and Chen, 2012): If the predictor has a stable relationship with the predictand, the predictability of the statistical downscaling scheme developed using this predictor will be stable; otherwise, a predictor with an unstable relationship with the predictand could result in unstable predictability of the statistical downscaling scheme. As the TSDTR approach considers the modulation effect of decadal-scale coupled oceanic-atmospheric modes on interannual climate variability, it may provide a new perspective to improve climate prediction. In this paper, we demonstrate that the TSDTR downscaling forecasting skill is superior to three other simpler methods (P5YSA, PLS and TSD-PLS), without considering the preceding SST conditions in different PDO phases. The results suggest that the TSDTR approach has a higher predictive capability SCESR.

In spite of the improvement in forecasting skill, the TSDTR downscaling forecasting model for SCESR encounters a problem in that it fails to predict some extreme flood events, such as in 2001. This failure is attributable to the interannual rainfall model, because the preceding interannual spring SST condition used as the predictor in 2001 displays no significant anomalies (not shown), indicating that the SST anomaly itself is not the only factor that can exert an important influence on the interannual SC rainfall. Thus, research gaps exist insofar as other potential factors need to be considered in the development of the TSDTR model. We expect in the future to include some of these other factors [e.g. the atmospheric circulation variability over the extratropics and tropics, and soil moisture conditions (e.g. Thompson and Wallace, 1998; Zhao et al., 2007; Chen et al., 2014)], with the aim to improve the statistical forecasting skill further.

PLS regression

PLS regression is usually described as a two-staged approach. Following \cite{Butler2000},the first stage is to produce a sequence of $k\le m$ PLS components $Z_i$ of length $n$ to be included in the regression. However, unlike principal components, which are only formed by accounting for the maximum amount of joint variability of ${X}$, the PLS component $Z_i$($i=1,2,\ldots,k$) is chosen to explain as much as possible the covariance between ${X}$ and $Y$. That is,$$Z_i={Xc}_i \ \ (A1)$$ is chosen to maximize $Z'_iY$ subject to the loadings ${ c}_i$ that $\|{ c}_i\|=1$ and that $Z_i$ is orthogonal to the space spanned by the basis $\{Z_1,Z_2,\ldots,Z_{i-1}\}$. The second stage of PLS regression is to regress $Y$ on the PLS components $Z_i$($i=1,2,\ldots,k$), which gives a simple linear model $$ Y=\sum_{i=1}^k {b_iZ_i}+\varepsilon_k~,(A2)$$ where $\varepsilon_k$ is the appropriate error term. The parameter $b_i$ is derived by minimizing the sum of squares,$$\left(Y-\sum_{i=1}^k{b_iZ_i}\right)'\left(Y-\sum_{i=1}^k{b_iZ_i}\right)=(Y-X\beta_{{\rm PLS}})'(Y-X\beta_{{\rm PLS}})~,(A3)$$ where $\beta_{{\rm PLS}}=\sum_{i=1}^k{b_ic_i}$ is the PLS parameter vector. From \cite{Helland1988},we have $$\hat{\beta}_{{\rm PLS}}=\sum_{i=1}^k\hat{\gamma}_i({X}'{X})^{i-1}{X}'Y~,$$ where the parameter $\hat{\gamma}_i$ can be estimated by using the PLSR1 algorithm. Predictions $\hat{y}_{x'}$ of future responses can then be made by $$ \hat{y}_{x'}=\overline{y}+\sum_{j-1}^k\hat{\beta}_{{\rm PLS},j}(x'_j-\overline{x}_j)~.$$

Note that the selected $Z_i$($i=1,2,\ldots,k$) is obtained by accounting for the maximum amount of the covariance between ${X}$ and $Y$. Thus, PLS components are obtained by not only accounting for the variance of explanatory variables ${X}$, but also the variances in the predicant $Y$.