A Possible Approach for Decadal Prediction of the PDO

The Pacific Decadal Oscillation (PDO) is a leading mode of decadal sea surface temperature variability in the North Pacific. Skillful PDO prediction can be beneficial in many aspects because of its global and regional impacts. However, current climate models cannot provide satisfied decadal prediction of the PDO and related decadal variability of sea surface temperature. In this study, we propose a new approach, i.e., the increment method, to predicting the PDO. A series of validations demonstrate that the increment method is effective in improving decadal prediction of PDO and it can well capture the phase change of PDO with high accuracy. The prediction processes include three steps. First, a five-year smoothing is performed; second, effective preceding predictors for PDO are selected, with all predictors and predictands in the form of a three-year decadal increment (DI); third, the prediction model is set up for PDO three-year decadal increment (DI_PDO), and PDO prediction can be obtained by adding the predicted DI_PDO to the observed PDO three years ago. This new method can also be applied for decadal climate prediction of other modes (e.g., Atlantic multidecadal oscillation) and predictands (e.g., sea surface temperature).


Introduction
With great socioeconomic and environmental significance, decadal/near-term climate prediction covering the period up to 30 yr has attracted lots of attention recently, especially for the policy-and-decision makers (Meehl et al., 2009). It has been addressed as a major issue in the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. The Decadal Climate Prediction Project (Boer et al., 2016) has been found to make contribution to the sixth Coupled Model Intercomparison Project (Eyring et al., 2016) and the World Climate Research Program (WCRP) Grand Challenge on Near Term Climate Prediction.
Pacific Decadal Oscillation (PDO), which is one of the dominant internal variability in the climate system, shows great influence on the other components of the cli-mate system (Newman et al., 2016). For instance, its regime shifts have been suggested to contribute to the decadal climate variation over East and South Asia (Yu et al., 2015;Zhu et al., 2015;Fan and Fan, 2017), Australia (Arblaster et al., 2002) and North America (Mc-Cabe et al., 2004, and to the recent global-warming hiatus (Kosaka et al., 2013). Besides, the agriculture, water resources, and the fisheries can also be affected by the PDO (Mantua et al., 1997;Miller et al., 2004).
To predict the internal decadal climate variability is necessary but full of challenges. Currently, using initialized climate models to produce decadal climate prediction is the main technical approach. Successful strategy for initialization of the model is one of the biggest technical challenges affecting the quality of decadal climate prediction Zhou and Wu, 2017;Mehta et al., 2019). The initial shock is a long-standing problem in the model initialization He et al., 2017). In a case study, Mochizuki et al. (2010Mochizuki et al. ( , 2012 show that the initialization in their climate model can improve the predictive skills for the PDO to some extent. Generally, the initialization shows higher improved predictive skill ability in North Atlantic (e.g., Smith et al., 2010;Yang et al., 2013) than in North Pacific (Guemas et al., 2012;Kim et al., 2012;Newman, 2013;Meehl and Teng, 2014;Wang and Miao, 2018), because of the inherent sensitivity to initial state uncertainty  and the uncertainty in the mechanism for internally generated decadal climate variability in the Pacific . Moreover, the regime shift prediction for the PDO is sometimes considered to be challengeable because it may be largely and randomly forced (Newman et al., 2016). Some statistical/statistical-dynamical method are used to predict the surface air temperature in the following decades (e.g., Fu et al., 2011, Luo andLi, 2014), but none of them for the PDO. Therefore, the current decadal prediction of PDO is far from being successful.
It is worth noticing that the recent interannual increment approach seems to provide a new clue in the climate prediction . The predictand is the interannual increment of a variable (difference between the current year and the preceding year) instead of the traditional variable anomalies from the climatology. The final predicted variable is produced by adding the predicted interannual increment to the observed variable in the preceding year. Since the preceding observed variable contains decadal-interdecadal signals (Wang et al., 2010) and the interannual increment of the variable in the numerical climate model can reduce the model system bias (Huang et al., 2014), the interannual increment approach has made some achievements in the interannual climate prediction (e.g. Fan et al., 2008;Fan, 2009Fan, , 2010Fan and Wang, 2009;Tian and Fan, 2015).
In the traditional decadal prediction steps, the predictand in each following year is usually predicted firstly, and then the decadal variability is obtained through the filter or running mean. In this study, we will directly use the decadal variability of the PDO as the predictand. The increment approach is used to build the statistical model to improve the effective samples and obtain more information from the previous observation. Section 2 gives the data and method in details. Sections 3 and 4 show the prediction of the decadal variability for the PDO index and the SSTA over North Pacific, respectively. A summary is presented in Section 5.

Data and methods
In this paper, we focus on the decadal prediction of the PDO during winter season (December-January-February, DJF), because many of the physical processes in the North Pacific have undergone a substantial long-term intensification during winter (Mantua et al., 1997;Deser and Phillips, 2006;Yeh et al., 2011;. Five-year running mean is used to obtain the decadal variability, and the yearly mark represents the middle year of the 5-yr mean period. The decadal PDO index is defined as the leading principal component of the empirical orthogonal function of 5-yr running mean DJF SSTAs over the North Pacific (20°-70°N), in which the SSTAs are the anomalies from the climatological annual cycle after removing the global mean SSTs (Mantua et al., 1997).
Four types of monthly mean datasets are involved to predict the decadal PDO, including: sea surface temperature from NOAA Extended Reconstructed SST v3b on a horizontal resolution of 2° x 2° (Smith et al., 2008), sea surface height from SODA (Simple Ocean Data Assimilation) v2.1.6 on a horizontal resolution of 0.5° x 0.5° ( Carton and Giese, 2008), sea level pressure from NOAA-CIRES (Cooperative Institute for Research in Environmental Sciences) 20th Century Reanalysis version 2c on a horizontal resolution of 2° x 2° (Compo et al., 2011), and sea ice concentration on a horizontal resolution of 1° x 1° from the Met Office Hadley Centre (Rayner et al., 2003).
After obtaining the decadal PDO index, the difference of the PDO between the current year and the three years before is calculated as the 3-yr decadal increment (DI) of the PDO (namely DI_PDO) [see Eq. (1)]. For example, the DI_PDO in 1906 is the difference between PDO in 1906 and PDO in 1903.
where PDO i and PDO i−3 represent the decadal PDO index in current year and three years before the current year, respectively. DI_PDO i indicates the decadal increment of the PDO in current year. The DI_PDO is predicted firstly through a statistical model with three predictors, which are also in the DI form (DI_Predictors). DI_Predictors are calculated as the 3-yr decadal increment of the 5-yr running mean for the detrended seasonal mean predictors. The final predicted PDO is derived by adding the predicted DI_PDO to the observed PDO three years ago by using Eq. (2) below, represents the the statistical model predicted DI_PDO in the target year. is the decadal PDO index from the observation three years before the target year. PDO i indicates the final predicted PDO in the target year.
In order to avoid the artificial built-in skill, we use two modulated validation methods to verify the predictive skill of the statistical model. (1) Cross-validation with five years left out (Michaelsen, 1987): the predictand in a target year is predicted with the forecast model by using the training samples in all the years except the five years that are centered around the target year. This process is repeated for all the other target years to obtain the crossvalidated re-forecast for the whole prediction period. The first/last three years are predicted by leaving the first/last five years out from the prediction period.
(2) Independent hindcast: the training period is the forward-rolling 67 yr, which is from 72 (69) yr before the target year to 6 (3) yr before the target year in predictors (predictand) to avoid the prediction model using any information from the prediction period. For example, the DI_PDO in 1975(1976 is predicted by the DI_Predictors in 1972(1973 based upon the forecast model built up using the period of 1903-1969/1906-1972 (1904-1970/1907-1973) in DI_Predictors/DI_PDO.
The moving t test (MTT) with a 9-yr moving window is used to detect the abrupt change points of the decadal variation of the PDO. The student's t-test is used to examine the statistical significance, in which the effective sample size N* is computed following Bretherton et al. (1999) and Ding et al. (2012) as follows, where N is the number of available time steps, and r 1 and r 2 are the lag-1 autocorrelation of the two correlated variables, respectively. Figure 1a shows the spatial patterns of the decadal DJF global SST associated with the decadal PDO index (the first EOF of decadal DJF SSTAs over the North Pacific). This mode accounts for 27% of North Pacific decadal SST variance. A positive PDO features negative SSTAs over central and western North Pacific and positive SSTAs over eastern North Pacific and eastern tropical Pacific. The decadal DJF PDO index during 1901-2009, as well as the DI_PDO, is shown in Fig. 1b.

Decadal prediction of the PDO
The DI_PDO is treated as the predictand. The empirical statistical model is built to predict the DI_PDO using DI_Predictors. This empirical statistical model is built through three steps. First, the potential preceding DI_Predictors are derived from the following three main processes driving the PDO as summarized by Newman et al. (2016): a) fluctuations in the Aleutian low; b) teleconnections from the tropics; and c) midlatitude ocean dynamics and coupled variability. Second, the DI_Predictors have to lead the DI_PDO at least three years to avoid  the DI_Predictors using any information from the prediction period, since we use 5-yr running mean to derive the decadal variability. Third, we use stepwise regression to identify the final DI_Predictors, which should be important and less dependent, to build the empirical statistical model. The fundamental rules in the stepwise regression are to select the DI_Predictor that is most significantly correlated to the predictand, and remove those predictors that are significantly related to this DI_Predictor, thus the DI_Predictors selected this way are independent to each other. We use the 99% confidence level for Fisher's F test to select the final DI_Predictors.
Based on the above steps, three leading DI_Predictors for three years are selected, which are autumn DI of Aleutian low (DI_AL), winter DI of sea ice over of Greenland Sea (DI_SIC), and spring DI of sea surface height over central Pacific (DI_SSH), respectively [Eq.
(4)]. Physically, the SST variability over most of the North Pacific is driven primarily by the atmospheric forcing (Smirnov et al., 2014). Generally, the atmospheric variations lead the SST variations that are resulted from the physical processes of the flux-driven SSTA pattern or Ekman transports (Davis, 1976;Deser and Timlin, 1997). The preceding Aleutian low (AL) has been suggested to be one of the important atmospheric forcing to drive the PDO variability (Schneider and Cornuelle, 2005;Newman et al., 2016) (Fig. 2a). On the other hand, the SIC anomalies over Greenland Sea is connected to the PDO variability through influencing the Arctic Oscillation anomalies ( Fig. 2b) (Lindsay and Zhang, 2005;Sun and Wang, 2006). Additionally, the sub-Arctic frontal zone (SAFZ) over the western Pacific has the large SST variations associated with the PDO (Nakamura and Kazmin, 2003). The preceding thermal capacity over SAFZ can drive the PDO variability through thermodynamic response ( Fig. 2c) (Qiu, 2003;Newman et al., 2016). though three of the DI_Predictors are all significantly correlated with DI_PDO, it is inevitable that these relationships have decadal shifts since the predicted period is more than 100 years (Fig. 3). It is noted that both of the relationships in DI_AL and DI_SSH become weakened around the period of 1960-1980, which may influence the predictive effect in the statistical model during this period (Figs. 3a, b). The moderate explained variances from the three DI_Predictors for the DI_PDO may also hint this decadal shift in the correlations. Figure 4a shows the predicted DI_PDO during 1906-2009 in the cross-validation with 5 years left out. As we can see, the limited skill for the predicted DI_PDO appears over the period around 1960-1980, which may be attributed to the above mentioned weakened relationships in DI_AL and DI_SSH. Excluding the period of 1960-1980, the predicted DI_PDO shows consistency in  The predicted DI_PDO is added with the PDO three years ago in the observation to obtain the final predicted PDO (Fig. 4c). The inconsistency between the predicted and observed PDO during 1960PDO during -1980 can also be seen in Fig. 4c. In the rest of the prediction period, the final predicted PDO fits the observation very well in amplitude and variability and generally successfully captures the phase change at each time. The CC between the ob-served and final predicted PDO reaches 0.82 during 1906-2009 (Fig. 4c). Actually, compared to regime shifts in the observation, the errors for each regime shift in prediction is less than or equal to two years (Fig. 4e). Overall, the increment method employed in this study shows effective skill for the decadal prediction of the PDO.
The predictive skill of the statistical model is further investigated by examining the independent hindcast for predicting the DI_PDO during 1975-2009. The magnitude of the predicted DI_PDO during 1975-1980 is lower than 0.5, which is far less than the observation, indicating a limited predictive skill. After 1980, the statistical model displays a high predictive skill for the DI_PDO. The CC between the prediction and observation of the DI_PDO during 1975-2009 can reach 0.78 (significant at the 99% confidence level) (Fig. 4b). The final predicted PDO during 1975PDO during -2009 has the CC of 0.81 with the observation (Fig. 4d). There are two regime shifts in observed PDO during 1975PDO during -2009PDO during , which are 1988PDO during /1989PDO during and 1999PDO during /2000 The final predicted PDO has the first regime shift at 1990/1991, which is two years later than the observation. The second regime shift in the final predicted PDO is at 1999/2000, which is precisely consistent with the observation (Fig.  4f). Overall, the statistical model combined with the increment method presents a high skill for the decadal prediction of the PDO, including the prediction for the regime shift.
In order to identify the effect of the increment method, the statistical prediction model in its original form (i.e., without the increment method) is used to directly predict PDO [Eq. (5)]. According to the correlation map of the PDO and the 3-yr leading predictors (Fig. 5), the key areas for the AL and SSH move slightly to the highest correlation regions compared to the key areas for DI_AL and DI_SSH. The key areas for the SIC move to the north of Kara Sea, since the region in Greenland Sea is not significant any more.
where AL indicates the September-October-November   show certain skill in predicting the PDO with significant correlation coefficient between the prediction and the observation (Figs. 6a, c). However, the cross-validated PDO shows that the extreme negative period of PDO index around 1950 has been missed (Fig. 6a). In the independent hindcast for the PDO during 1975-2009, the magnitude of the prediction is much lower than the observation and the predicted PDO keeps in the negative phase after the late 1980s, which is incorrect ac-cording to the observation. Generally, the predictive skill of the statistical predictive model without the increment method is lower than that with the increment method.

Prediction for the decadal variability of the SSTA over North Pacific
The statistical predictive model combined with the increment method is also employed to predict the decadal 1900 1920 1940 1960 19080 2000 1975 1980 1985 1990 1995 1900 1920 1940 1960 19080 2000 1975 1980 1985 1990 1995 2000 2005 1900 1920 1940 1960 19080 2000 1975 1980 1985 1990 1995 For the DI_SSTA during 1975-2009, a significantly high predictive skill appears over the Kuroshio-Oyashio Extensions to the central Pacific and the region north of 50°N (Fig. 7a). Although the predictive skill for the final predicted SSTA is not impeccable (Fig. 7b), the final predicted SSTA successfully demonstrates a regime shift at 1999 (Fig. 7c), which is a challenge for current dynamic numerical models (Newman et al., 2016). Moreover, we further investigate the predictive skill of this statistical model combined with the increment method to produce the spatial pattern of the PDO. Figure 7d gives the CC between the final predicted PDO index (Fig. 4d) and the final predicted SSTA over the North Pacific (Fig. 7b) in the independent hindcast. Corresponding to the positive PDO, there are negative SSTAs appearing over the western and central North Pacific and positive SSTAs appearing over the eastern North Pacific, resembling the observation (Fig. 7e). Actually, the pattern correlation coefficients between predictions and observations reaches 0.95. The statistical model combined with the increment method shows a high predictive skill for the regime shift of the SSTA over North Pacific and the spatial pattern of the PDO.

Conclusions and discussion
In this study, the decadal variability of the PDO is predicted by a statistical model combined with the increment method. Particularly, all the involved data have been detrended and applied with a 5-yr running meaning to maximize the decadal variability. The 3-yr increment of the PDO (DI_PDO) is predicted firstly by the statistical model with DI_AL, DI_SSH, and DI_SIC at preceding three years. The predicted DI_PDO is then added with the PDO three years ago in the observation to derive the final predicted PDO. The statistical model combined with the increment method shows a high predictive skill for the decadal variability of the PDO, especially for its regime shifts.
In addition to the fact that the increment method can take advantage of the previous observation, another benefit is that the increment method can increase the effective sample size after we apply the 5-yr running mean to the meta data. For example, for the correlation between DJF PDO during 1906-2009and DJF SIC over Greenland Sea during 1903, the effective sample size is 8, but it is 19 in correlation between DI_PDO and DI_SIC. The increment method provides a possibility for directly focusing on the decadal variability prediction of the PDO.  Huang, Y. Y., and H. J. Wang 69 1900Wang 69 1920Wang 69 1940Wang 69 1960Wang 69 19080 2000Wang 69 1975Wang 69 1980Wang 69 1985Wang 69 1990Wang 69 1995Wang 69 2000Wang 69 2005Wang 69 1900Wang 69 1920Wang 69 1940Wang 69 1960Wang 69 19080 2000Wang 69 1975Wang 69 1980Wang 69 1985Wang 69 1990Wang 69 1995     Our results provide a promising clue in current decadal climate prediction, especially for the regime shift prediction. This new method can be applied for decadal climate prediction on other modes (e.g., Atlantic multidecadal oscillation) and predictands (e.g., sea surface temperature). Additionally, the increment method may be applicable to the current decadal prediction products of climate models. The DI of the variable may be more predictable in the numerical model than the original form of the variable, since the DI form can partly reduce the system bias in the climate model. We can add the DI form prediction in climate models to the observed value at previous year to improve the model's prediction. Also, we can build a dynamical-statistical model using the DI form prediction in climate models and previous observations, to provide more realistic decadal climate prediction.