Seasonal prediction and predictability of the Asian winter temperature variability

Efforts have been made to appreciate the extent to which we can predict the dominant modes of December–January–February (DJF) 2 m air temperature (TS) variability over the Asian winter monsoon region with dynamical models and a physically based statistical model. Dynamical prediction was made on the basis of multi-model ensemble (MME) of 13 coupled models with the November 1 initial condition for 21 boreal winters of 1981/1982–2001/2002. Statistical prediction was performed for 21 winters of 1981/1982–2001/2002 in a cross-validated way and for 11 winters of 1999/2000–2009/2010 in an independent verification. The first four observed modes of empirical orthogonal function analysis of DJF TS variability explain 69 % of the total variability and are statistically separated from other higher modes. We identify these as predictable modes, because they have clear physical meaning and the MME reproduces them with acceptable criteria. The MME skill basically originates from the models’ ability to capture the predictable modes. The MME shows better skill for the first mode, represented by a basin-wide warming trend, and for second mode related to the Arctic Oscillation. However, the statistical model better captures the third and fourth modes, which are strongly related to El Niño and Southern Oscillation (ENSO) variability on interannual and interdecadal timescales, respectively. Independent statistical forecasting for the recent 11-year period further reveals that the first and fourth modes are highly predictable. The second and third modes are less predictable due to lower persistence of boundary forcing and reduced potential predictability during the recent years. In particular, the notable decadal change in the monsoon–ENSO relationship makes the statistical forecast difficult.


Your article is published under the Creative
Commons Attribution license which allows users to read, copy, distribute and make derivative works, as long as the author of the original work is cited. You may selfarchive this article on your own website, an institutional repository or funder's repository and make it publicly available immediately.

Introduction
During the 2009/2010 winter season, many parts of Asia experienced conspicuous climate anomalies concurrent with the Central Pacific El Niño and strong negative phase of the Arctic Oscillation (AO) in which frequent severe cold surges (outbreaks of cold air) and record-breaking heavy snowfall were observed. During the first few days of January 2010, North China and Korea dropped to 50-and 70-year record low temperature and snowfall, respectively. These events resulted in heavy losses from the viewpoint of agriculture, transportation, and other socioeconomic activities. The cold winter over the most parts of Asia seems against the global warming trend projected by the contemporary climate models (e.g. Lee and Wang 2012b). Although economic and social influences are equally affected by salient variability in the Asian winter monsoon (AWM) and that in the summer, prediction of the AWM Wu et al. 2011;Sohn et al. 2011) has received less attention than that of the summer counterpart which has turned out to be the most challengeable in current climate models and observation shown by numerous studies (e.g. Kang et al. 2002;Ha et al. 2005;Kang and Shukla 2006;Yang et al. 2008;Wang et al. 2004Wang et al. , 2008Wang et al. , 2009Lee et al. 2010Lee et al. , 2011aLee and Wang 2012a and many others). The extent to which we can predict dominant modes of the AWM variability with dynamical and with physically based empirical models is still unclear.
The identification of the dominant modes of AWM variability is a crucial step in understanding its predictability. Wang et al. (2010) identified two dominant modes of air temperature variability over the East Asian winter monsoon (EAWM) region, which were called the northern and southern mode respectively. The former represents a cold winter in northern East Asia due to cold air intrusion from northeastern Siberia, and the latter is characterized by deepening of the East Asian trough and strengthening of the Mongolian High. They showed that nine of ten existing EAWM circulation indices (Ji et al. 1997;Cui and Sun 1999;Lu and Chan 1999;Chen et al. 2000;Li and Zeng 2002;Jhun and Lee 2004;Wu et al. 2006) basically describe the southern modes and one index (Guo 1983) describes both modes equally. It was further demonstrated that the two dominant modes can explain a significant amount of temperature variability over the entire Asian region.
Several factors have been determined to influence interannual and interdecadal variability of the AWM. Many observational and modeling studies have shown that the decrease of snow cover over the Eurasian continent during the previous fall season results in weakening of the EAWM (Walland and Simmonds 1997;Watanabe and Nitta 1999;Clark and Serreze 2000;Jhun and Lee 2004). The importance of autumn Arctic sea ice to abnormal AWM climate was suggested (Wu et al. 2011;Li and Wu 2012). It is found that the joint action of western Pacific Subtropical High and Siberian High affects cold wave frequency in China (Ma et al. 2012). Moreover, El Niño-Southern Oscillation (ENSO) plays an important role in regulating the interannual variability of the AWM (Zhang et al. 1996;Wang et al. 2000;Chan and Li 2004;Chang et al. 2004;Wang et al. 2010). On the other hand, several studies have shown that AO and North Atlantic Oscillation (NAO) are related to the AWM on a decadal time scale (Gong et al. 2001;Jhun and Lee 2004;Wu et al. 2006;Li and Bates 2007). Wu et al. (2009) indicated that the Southern Hemisphere annular mode is also linked to the AWM variability.
Despite of the aforementioned studies, how to determine predictability of AWM variability on a seasonal time scale is still elusive. Better understanding of the origins and predictability sources of the AWM may contribute to the improvement of its seasonal prediction capability. Wang et al. (2007) and Lee et al. (2011a) suggested the concept of predictable mode analysis (PMA) because the multi-model ensemble (MME) seasonal prediction skill is based essentially on the major modes. With this method, attainable potential climate predictability can be quantified using the fractional variance of the ''predictable'' leading modes determined by examination of observation and hindcast results of the models. The purpose of this study is to explore the attained prediction skill and achievable potential predictability of the AWM variability by using the PMA approach and to understand the range for potential improvement of dynamical and statistical model predictions. The specific questions to be addressed include (1) how to determine the current level of dynamic prediction skill for the AWM seasonal mean 2 m air temperature (hereafter, TS) anomaly and identify the source of the dynamical forecast skill; (2) how to empirically estimate the potential predictability of the AWM; and (3) how to optimally construct a statistical model with physically based predictors for the AWM TS prediction. These points will be respectively addressed in Sects. 3, 4 and 5, respectively, following a brief discussion of the datasets, models, and methodology used in this study. Section 6 summarizes our major findings.
2 Models, data, and methodology 2.1 Coupled models' hindcast data Thirteen coupled models used for the present study were adopted from Asia-Pacific economic cooperation climate center (APCC)/climate prediction and its application to society (CliPAS) and development of European multimodel ensemble system for seasonal-to-interannual prediction (DEMETER) projects; hindcast datasets are detailed by Wang et al. (2009) and Lee et al. (2010). A brief summary of these models is presented in Table 1. The models show a large range of model horizontal resolution and number of ensemble members, from 6 to 15. Flux correction was not applied to any of the selected coupled models and these models have retrospective hindcasts with the November 1 initial condition that yield 1-month-lead winter DJF forecasts for the common period of 1981/1982-2001/2002. MME prediction was made using simple average of the ensemble means of the 13 coupled models.

Observed data
The observed TS data were obtained from the National Centers for Environmental Prediction (NCEP) reanalysis version 2 (NCEP R2; Kanamitsu et al. 2002). The monthly Niño 3.4 SST index was calculated using improved Extended Reconstructed SST Version 3 (ERSST V3) data (Smith et al. 2008). All data used in this study were normalized by their own standard deviations. Wang et al. (2007) and Lee et al. (2011a) suggested a method for determining predictable modes from empirical orthogonal function (EOF) modes of observation and stateof-the-art climate model prediction. The correlation matrix is used for EOF analysis because dynamical models more effectively capture observed dominant modes with this approach than with that of the covariance matrix for the AWM TS anomalies. Two basic criteria were considered for the determination: (1) In observation, predictable modes should explain a large part of the total variability with physical interpretations and should be statistically separated from other higher modes; and (2) the climate prediction models should be capable of predicting these major modes.

Predictable mode analysis
To represent climate model capability in predicting EOF modes, a skill score is defined in terms of the pattern correlation coefficient (PCC) skill for eigenvector (EV) and temporal correlation coefficient (TCC) skill for the principal component (PC) time series of each mode (i). The skill score for each mode is calculated by The skill score ranges from 0 to 1, respectively indicating no skill at all and perfect forecast. It should be noted that we reordered the EOF modes of the MME prediction according to the skill score because the order of the predicted EOF mode is not necessarily the same as its observed counterpart. To reorder the predicted EOF modes, the skill score for the first observed mode was first calculated against all of the predicted modes, and the predicted mode with the best skill score in the first observed mode was taken as the first predicted mode. Other predicted modes were determined by repeating this process.
From the determination of predictable modes, the total field (TS) is decomposed into the predictable and unpredictable parts. The predictable part (TS pred ), as a function of longitude (lon), latitude (lat), and time (t), is reconstructed by the linear combination of the predictable EOF modes defined by where k i is the eigenvalue of ith mode, and N is the total number of predictable modes. The unpredictable part (TS unpred ) is then calculated by subtracting the predictable part from the total field. We consider the correlation coefficient between the original TS and TS pred as a measure of attainable potential predictability assuming that predictable PCs are accurately predicted (Lee et al. 2011a; Lee and Wang 2012a).

Statistical model
We attempted to determine the extent to which we can predict the predictable dominant modes of AWM variability with a physically based statistical model. In this study, the statistical model, named predictable mode forecast model (PMFM), has two steps: (1) prediction of the predictable EOF PCs of DJF TS over the AWM region using physically based and optimally selected predictors; and (2) reconstruction of the TS field over the entire AWM region from predicted PCs with determined EV and eigenvalues during the training period. A simple multiple linear regression model was used for prediction of each PC defined by where STM_PC i is the predicted PC for the ith EOF mode; t, the time for forecast target; M, the total number of predictors for each PC; Pred ij , the jth predictor for the ith PC normalized by its own standard deviation with time lag s (1 month in this study); and a ij , the coefficient of the jth predictor for the ith PC. The most important factor for accurate prediction of PCs in the PMFM is the determination of optimal predictors for each predictand. The details of predictors for each PC will be addressed in Sect. 4.  1979/1980-1998/1999 (1989/1990-2008/2009) DJF. The predicted DJF TS was determined by the following equation: We considered the correlation coefficient between the original TS and STM_TS training as a measure of attained fitting skill and that between the original TS and STM_TS forecast as a measure of attained forecast skill for the PMFM.

Evaluation of AWM TS prediction
In this section, we first examine the AWM prediction skills of the 13 coupled models and their MME predictions. The strength of the AWM is measured by TS over the Asia-Western North Pacific (WNP) region (Eq-60°N, 60°E-140°E) in the period of DJF. Figure 1 shows the TCC skills of the 1-month-lead seasonal prediction for DJF TS obtained from the 13 coupled models  and their MME ( Fig. 1a) with November 1 initial condition over the AWM region for the 21-year period of 1981/1982-2001/2002. Although most coupled models accurately predict anomalies over the tropical ocean due to persistent SSTs 1 month in advance, prediction of the DJF TS anomaly over the continental AWM region is difficult. For the hindcast period, the TCC skill averaged over the entire domain (Eq-60°N, 60°E-140°E) for individual models ranges from 0.16 to 0.46 with the averaged skill of 0.33. In comparison, the all-model MME is capable of predicting the AWM 1 month in advance with the domainaveraged TCC skill of 0.51 for the 21 years, which is better than that observed for any individual coupled model and significantly higher than the averaged skill of all individual models.

Prediction of major EOF modes
The first four major EOF modes of AWM TS variability are investigated in observation and in a 1-month-lead MME prediction. Figure 2 shows spatial distribution of the EV and PC time variation in addition to the percentage variance accounted for by each EOF mode, shown at the top of each panel of the spatial structure. Moreover, the MME PCC skill for EV and the TCC skill for PC are presented for each mode. A comparison of spatial structure and temporal coefficient indicates that the second, third, and fourth predicted modes respectively correspond to the fourth, second, and third observed modes, suggesting that the MME tends to highly overestimate the variance of the observed fourth mode and underestimates that of the observed second and third modes. Thus, the predicted modes were reordered following the procedure introduced in Sect. 2.3.
The observed first mode represents a domain-wide warming trend across the period examined and accounts for 32.1 % of the total observed variability. This mode is the same as the second EOF mode of the EAWM (southern mode) shown in Wang et al. (2010). The MME captures its spatial pattern and temporal variation 1 month in advance with high fidelity. The PCC skill for EV 1 is 0.95, and the TCC skill for PC 1 is 0.83.
The second observed EOF mode features a north-south seesaw pattern with an interannual timescale, accounting for 16.6 % of the total observed variability. This mode resembles the first mode of the EAWM (northern mode) shown in Wang et al. (2010). The MME accurately captures the spatial and temporal characteristics of the second observed mode but underestimates its variance. The PCC skill for EV 2 is 0.72, and the TCC skill for PC 2 is 0.50.
The observed third mode displays a sandwich pattern characterized by a cold TS anomaly over the mid-latitude Asian continent with a warm TS anomaly over the highlatitude Asian continent and the Tropics and vice versa. This mode has mainly interannual variability and accounts The TCC skill for 1-month lead DJF prediction of 2 m TS obtained from 13 coupled models and their MME with November 1 initial condition over the AWM region for the period of The numbers in the left upper corners indicate the averaged correlation skill over the AWM region for 11.2 % of the total observed variability. The MME reasonably predicts the spatial and temporal characteristics, although its accuracy is less than that of the first two observed modes. The percentage variance is significantly underestimated, and spatial errors are exhibited over most of the oceanic region. In addition, the strong positive phase of the observed PC in DJF 1982DJF /1983DJF and 1997DJF /1998 is not well captured. The fourth observed mode exhibits prominent decadal variability and three variability centers over the WNP (negative anomaly), northwestern China (positive), and Kazakhstan (negative). Interestingly, the second predicted mode most resembles the fourth observed mode. However, the MME significantly overestimates its percentage variance and had remarkable spatial errors over most of midlatitude Asia. Nonetheless, the MME offers some skill in predicting the fourth observed mode with a PCC skill of 0.36 and a TCC skill of 0.39. indicates PCCs between the observed and predicted EV. The value of TCC indicates temporal correlation coefficients between the observed and the predicted PC time series predictable modes were identified for the AWM TS variability. Figure 3 shows a scatter diagram between the observed percentage variance (ordinate) and the MME skill score (abscissa), defined by Eq. (1), for each EOF mode. The first two observed EOF modes are statistically well separated from higher modes according to the rule of North et al. (1982) and are predicted with high fidelity by the MME of the current coupled models. The skill scores were 0.89 and 0.6 for the first and second EOF modes, respectively. The third and fourth modes, to a certain extent, are separated from other higher modes and captured by the MME with lower fidelity than the first and second. Based on the result, we defined the first four modes as predictable modes for DJF TS over the AWM region. These modes account for approximately 69 % (83 %) of the total variance of DJF TS in observation (MME prediction).

Prediction skill and attainable potential predictability
This subsection describes the current level of prediction skill attained by the 1-month-lead MME prediction and attainable potential predictability from the first four predictable modes. The total predicted and observed TS were decomposed into predictable (TS pred ) and unpredictable (TS unpred ) parts as described in Sect. 2.3. Reconstruction from PCs was achieved using Eq. (2). Figure 4a, b compare TCC skills for the DJF TS obtained from TS pred and TS unpred components, respectively, of the MME prediction.
For the calculation, the observed total field was used for both predictable and unpredictable cases. The skills of the dynamical model predictions made by the four predictable modes are essentially comparable to the original prediction using all empirical modes (Fig. 1a). On the contrary, the contribution of all residual modes to seasonal prediction skill is minimal, except for those over the Indochina peninsula. This result clearly indicates that the current coupled model prediction skill for the DJF AWM TS essentially originates from the skill for prediction of the first four modes.
To improve the 1-month-lead MME prediction, we applied statistical postprocessing to the prediction. Since the MME prediction exhibited errors in capturing the spatial pattern of the observed EVs, the predicted EVs were replaced by the observed values during the reconstruction procedure for the predictable modes of the MME prediction. The statistical postprocessing potentially improved the dynamical prediction over most of the AWM region shown in Fig. 4c. However, if the statistical postprocessing was applied to the MME prediction with a cross-validated approach, no significant improvement was observed (not shown).
From the conventional perspective, potential predictability can be defined by the fractional variance of the predictable part. In this case, 69 % of the total observed variability was potentially predictable over the AWM region. In present study, however, the realizable potential predictability was estimated by the TCC between the observed total field and the observed TS pred to facilitate comparison with the MME prediction skill, as mentioned in Sect. 2.3. Figure 4d indicates that the variability of the DJF TS is highly predictable, particularly over the EAWM region, the tropical oceans, and a large region of northern Asia. The predictability is relatively low over northern India. The area-averaged TCC skill is 0.82, which corresponds to 69 % of the total observed variability in the linear regression.

Sources of predictability
The results presented in Figs. 3 and 4 suggest that the first four leading observed EOF modes are predictable. To support this postulation, we further examined the sources of predictability of the modes. Figure 5 shows TCC between the DJF TS anomaly and each mode for the 21-year period of 1981/1982-2001/2002. The first observed mode, which exhibits a warming trend over most of the AWM, significantly correlates with SST warming over the WNP, North Indian Ocean, and subtropical North Atlantic Ocean and with continental TS warming over North America, middle Africa and the Maritime Continent (Fig. 5a). In addition, this mode significantly correlates with Niño 3.4 SST anomalies observed in the previous September-October-November (SON) with a TCC of 0.4 (Fig. 6). Fig. 3 The percentage variances that are accounted for by the observed first eight EOF modes (ordinate) and the combined forecast skill score for the EV and PC for each mode (abscissa) for DJF TS over the AWM region. The first four major modes capture about 69 % of the total observed interannual variability The second observed mode shows strong positive correlation with the DJF TS over most of northern Asia and Europe (Fig. 5b). This mode exhibits no relationship with the simultaneous Niño 3.4 SST anomaly but may be slightly related to the developing ENSO during following June-July-August (JJA), shown in Fig. 6. It is noted that this mode has significant correlation with the simultaneous AO with a TCC of 0.66 for the 21-year period of DJF 1981DJF /1982DJF -2001DJF /2002. A significant positive correlation (0.68) remains for the recent 20-years period of DJF 1990of DJF /1991of DJF -2009of DJF /2010. ENSO regulates the third observed EOF mode (Fig. 5c) on the interannual time scale. The lead-lag relationship between the seasonal Niño 3.4 SST anomaly and the third observed PC clearly indicates that the third mode is observed during the mature phase of ENSO (Fig. 6). During the mature phase of El Niño, the cold TS anomaly tends to occur over most of mid-latitude Eurasia, and the warm anomaly occurs over most of northern Eurasia and the Equatorial and South Indian Ocean.
The fourth observed mode is related to the decadal variability of the North Pacific Ocean variability (Fig. 5d) (a) (b) (c) (d) Fig. 4 The MME's TCC skill for DJF TS using a reconstructed field from the first four modes, b reconstructed field from higher modes, and c statistical correction using the first four modes. d The realizable potential predictability from the first four predictable EOF modes.
Solid ( and has a prolonged relationship with the Niño 3.4 SST anomaly from the previous JJA to the following DJF (Fig. 6), indicating that the decadal component of ENSO and the North Pacific Ocean variability regulate the mode. The fourth observed PC is closely related to the long-term TS variability over the eastern Tibetan Plateau and northwest China, where the TS variability shows a weaker correlation with the observed PC 1 than that over the other Asian region indicated in Fig. 5a. It is further noted that PC 4 has no relationship (nearly zero) with the AO for the 21-year period of DJF 1981DJF /1982DJF -2001DJF /2002; however, a significant negative correlation of -0.4 is evident for the recent 20-year period of DJF 1990DJF /1991DJF -2009DJF /2010. Interestingly, the second predicted PC, which resembles the fourth observed PC, has a significant relationship with the predicted ENSO on both interannual and interdecadal time scales (not shown).

Statistical prediction
This section demonstrates the achievable prediction skill of the AWM using the physically based statistical model (predictable modes forecast model; PMFM) introduced in Sect. 2.4. In this approach, the first four predictable PCs are first predicted with optimally selected predictors, and the prediction of the DJF TS over the entire AWM is achieved from TS reconstruction using the observed EVs and statistically predicted PCs from Eq. (5). The selection of predictors plays a crucial role on improving the forecast skill of the PMFM.

Sources of predictability
To facilitate comparison with the 1-month-lead MME prediction, the preceding September-October mean (SO) TS was used for selecting predictors for the first four predictable PCs. We first investigated the correlation coefficients between the SO mean TS over the globe and each PC time series of the DJF AWM TS shown in Fig. 7. A comparison between Figs. 5 and 7 reveals remarkable persistent characteristics of lower-boundary forcings related to the predictable PCs, except for the second PC, which indicates an existing perspective on the statistical forecast. During the SO prior to a weakened AWM (i.e. basinwide warming of the TS), pronounced warming is observed in the north Indian Ocean and in the subtropical-mid-latitude North Atlantic Ocean (Fig. 7a) that persists through the following DJF (Fig. 5a). In accordance with the noteworthy relationship, the predictors for the PC 1 (Pred 1j ) are defined by a 1j Pred 1j ðtÞ ¼ 0:5 r 1j CTS 1j ; j ¼ 1; 2 CTS 1j ¼ ½CORðlon; latÞ Á SO TSðlon; lat; tÞ Rj if CORðlon; latÞ j j ! 0:38; where COR is the correlation coefficient between SO TS and the first PC as a function of longitude and latitude, r is the standard deviation of CTS time series, square brackets indicate area averaging over region R1 (Eq-25°N, 60°E-130°E) and region R2 (Eq-60°N, 80°W-10°W) shown in Fig. 7a. The SO TS anomaly, which has a correlation coefficient larger than ?0.38 or smaller than -0.38 (95 % confidence level), is averaged over the selected region. For the forecast purpose, the correlation coefficient was obtained during training period. Contrary to that in PC 1 , outstanding precursors were not observed for PC 2 (Fig. 7b) over the Northern Hemisphere. In the simultaneous relationship, the positive AO is strongly concurrent with the warm (cold) TS anomaly over the northern (subtropical) AWM. Because of the chaotic nature of the AO, however, no significant relationship of PC 2 with the SO AO was observed. However, useful precursors were detected over the Barents and Kara seas, some where region R1 is 35°N-80°N, 20°E-160°E, and region R2 is 30°N-70°N, 180°-80°W. We have determined that the precursors for PC 2 are not optimal; therefore, further study is needed to find more effective predictors.
Since the positive phase of PC 3 tends to occur during the mature phase of El Niño, its correlation field with SO TS exhibits a clear progressing El Niño pattern over the tropical Indo-Pacific Ocean (Fig. 7c). In addition, a positive correlation is observed over Russia and Kazakhstan that may be related to SO snow cover over the region. The predictors for PC 3 (Pred 3j ) are defined by where region R1 is 20°S-20°N, 40°E-100°E; R2 is 20°S-20°N, 100°E-160°E; R3 is 20°S-20°N, 170°W-60°W; and R4 is 30°N-70°N, 40°E-90°E. The positive phase of PC 4 is strongly related to the cold SST anomaly over the WNP and the warm SST anomaly over the central and eastern Pacific (Fig. 7d) that persist through the following DJF. The predictors for PC 4 (Pred 4j ) are defined by where region R1 is Eq-50°N, 80°E-160°E, and R2 is 10°S-60°N, 170°E-120°W.

Cross-validated statistical forecast
To compare with the 1-month-lead MME prediction, we first attempted statistical forecasting for the first four PCs using the aforementioned predictors with a cross-validated approach for the 21-year period of 1981/1982-2001/2002. All coefficients and standard deviations of predictors were determined during the training period excluding the forecast year. Figure 8 shows the observed and empirically predicted PCs for the cross-validated method (Fcst). As a control forecast, a fitting result (Fit) that used all 21-year information for the statistical forecast is also displayed. In Fig. 8, r indicates the correlation coefficient between the observed and predicted PCs. Because of the overfitting problem of the statistical model, the cross-validated skill dropped in comparison to the fitted skill. However, the first and fourth PCs that exhibit trend and interdecadal variability, respectively, are highly predictable with the current empirical approach. The empirical model is also able to predict the third PC with less fidelity than the first and fourth PCs. It is noted that current dynamical MME has better skill for the first and second PCs than the PMFM whereas the PMFM relatively better captures the third and fourth PCs, which are strongly related to ENSO variability on interannual and interdecadal timescales, respectively, than the MME. The statistical prediction for the entire AWM TS was next performed using the reconstruction method from the fitted (Eq. 4) and predicted PCs (Eq. 5) introduced in Sect. 2.4. Figure 9a, b show the TCC skill for statistically fitted (STM_TS training ) and cross-validated prediction (STM_TS forecast ) for the 21-year period. The fitted skill is as good as the dynamical prediction skill with spatial bias correction (Fig. 4c), and the cross-validated prediction skill is comparable to the 1-month-lead MME prediction reconstructed from the first four PCs (Fig. 4a). While the dynamical model has better skill in predicting the TS over the subtropical AWM, the statistical model more effectively predicts the TS over the mid-latitude EAWM region. In addition, the spatial distribution of the statistically fitted skill is similar as that of the attainable potential predictability shown in Fig. 4d, indicating that the selected predictors for each PC are capable of capturing most of interannual-to-interdecadal variability of the AWM TS. The persistent forecast skill was also compared with the statistical and dynamical forecast (Fig. 9c). The persistent forecast was obtained with the assumption that the October mean TS persists through the following DJF season. A comparison of Fig. 9b, c indicates that the PMFM exhibits significantly better skill than the persistent forecast except over some parts of Central Asia and the Arabian Sea, where large persistence was observed.
The spatial distributions of the dynamical (Fig. 4a) and statistical forecast skills (Fig. 9b) are demonstrated to be, to some degree, independent and complementary to each other. The dynamical forecast has significant skill over most of the oceanic region, while the statistical forecast has skill over the mid-latitude EAWM region. From these results, we simply averaged the dynamical and statistical forecasts and investigated their resultant forecast skill. Figure 9d shows that a simple average of two forecasts can improve the DJF TS forecast over most of the AWM region. The area-averaged forecast skill for the combined forecast was 0.59, which is larger than that of the dynamical (0.53) and statistical (0.51) skills.

Independent forecast
Although the statistical forecast utilized the cross-validated approach, overfitting problems remain inevitable. Thus, we performed independent forecasting for the 11-year period of DJF 1999DJF /2000DJF -2009DJF /2010 to confirm whether the suggested methodology is actually useful. The independent forecast procedure is detailed in Sect. 2.4. Figure 10 shows the observed, fitted, and independently predicted PCs. Since the training period is used for the It is obvious that the independent forecast skill drops further than the fitted skill from the cross-validated forecast in the previous section. Nonetheless, the trend mode (PC 1 ) and interdecadal mode (PC 4 ) are highly predictable with TCC skills of 0.76 and 0.65, respectively, for the recent 11 years. However, two interannual modes of PC 2 and PC 3 are not well predicted. For the AO-related second PC (northern mode in Wang et al. 2010), the remarkable positive phases in 2001/2002, 2003/2004, and 2006/2007 are totally missed in the PMFM, while negative phases are relatively well predicted. The low performance of the ENSO-related third PC may be attributed to the interdecadal change in the ENSO-monsoon relationship. The correlation between the third PC and DJF Niño 3.4 SST index is 0. 56 for 1979/1980-1999/2000 and 0.37 for 1989/1990-2009/2010, which indicates a weakening of the relationship during the recent decade. This notable decadal change in the monsoon-ENSO relationship likely causes difficulties in statistical forecasting. Figure 11 shows the TCC skill for the DJF TS over the entire AWM by the reconstructed field from the fitted and independently predicted PCs along with the aforementioned potential predictability during 1999/2000-2009/2010. During the recent 11 years the TS over the subtropical AWM is more statistically predictable than that over the northern AWM. The area-averaged TCC skill is 0.43 for the 11-year period, which is less than the skill of 0.51 for the 21-year period of 1981/1982-2001/2002. It is important to note that the TS during the recent decade is less predictable than that of 1981/1982-2001. The attainable potential predictability is 0.71 over the AWM, which is  1981/1982-2001/2002 (0.82). In addition, the persistence of TS was also reduced from 0.24 (Fig. 9c) to 0.12 (Fig. 11d). Considering the lower persistence and potential predictability during the recent decades, the independent statistical forecast appears to be valid and useful.

Summary
The 1-month-lead seasonal forecast skills of current dynamical and statistical models were investigated on the DJF TS anomaly over the AWM region and were compared with attainable potential predictability obtained from the ''predictable'' leading modes for the 21-year period of 1981/1982-2001/2002. The 1-month-lead dynamical prediction for the DJF AWM TS was achieved by the MME of thirteen coupled models with the November 1 initial condition that were included in DEMETER and CliPAS projects. The MME is capable of predicting the AWM TS 1 month in advance with the domain-averaged TCC skill of 0.51 for the 21-year period, which is significantly higher than the 0.33 averaged skill of all individual models and better than individual model skills, which ranged from 0.16 to 0.46.
The first four observed modes are identified as predictable modes because they explains 69 % of the total observed variability with physical interpretations and are statistically separated from other higher modes. Moreover, they are also well predicted by the current climate prediction models, to some extent. The MME effectively captures the first and second observed EOF modes with high fidelity, which are characterized by a domain-wide warming trend and by AO-related interannual variability, respectively. The observed third and the fourth modes are respectively associated with interannual and interdecadal components of ENSO and are reasonably captured by the MME, although less faithfully than the first two modes. The MME highly overestimates percentage variance of the fourth observed mode but far underestimates that of the second observed mode.
The observed total field of the DJF AWM TS was decomposed into predictable and unpredictable components. Considering the assumption that the first four predictable modes are perfectly predicted, attainable potential predictability can be quantified by the TCC between total TS and reconstructed TS from the predictable modes. The area-averaged attainable TCC skill is 0.82, and this corresponds to 69 % of the total observed variability in linear regression. The MME tends to more accurately predict areas of high potential predictability except over that above 50°N of the AWM region. It is important to note that the MME forecast skill essentially originates from the ability of dynamical models to capture the first four predictable modes.
A statistical prediction with the PMFM was achieved using predictors obtained from the SO mean TS that significantly correlated with the predictable modes of the DJF AWM TS. During 1981During /1982During -2001During /2002, the cross-validated statistical forecast shows significantly high skills for the first and fourth observed modes, which exhibits trend and interdecadal variability, respectively. It is noted that the current dynamical MME shows better skills for the first and second PCs whereas the PMFM more effectively captures the third and fourth PCs, which are strongly related with ENSO variability on interannual and interdecadal timescales, respectively, compared to the MME. Since the dynamical and statistical predictions are complementary, a simple composite of two predictions is capable of improving the forecast skill of the DJF TS over most of the AWM region. The area-averaged forecast skill for the combined forecast is 0.59 which is larger than that of 0.53 dynamical and 0.51 statistical skills.
Independent statistical forecasting for the DJF AWM TS was achieved for the 11-year period of 1999/2000-2009/ 2010 to validate the statistical model's usefulness for realtime forecasting during recent years. The area-averaged TCC skill of the statistical forecast is 0.43 for the recent (a) (b) (c) (d) Fig. 11 The TCC skill for DJF TS for the period of 1999/2000-2009/ 2010 DJF using a realizable potential predictability measure, b fitting procedure, c independent forecast of 1-month lead statistical prediction, and d persistent forecast using October mean TS. Solid (dashed) line indicates a correlation coefficient of 0.6 (0.4). The numbers in the left upper corners indicate averaged correlation skill over the AWM region 11 years, which is less than that of 0.51 for the 21-year period of 1981/1982-2001/2002. It is important to note that the attainable potential predictability is 0.71 over the AWM, which is less than that of 0.81 for the 21-year period, and this result indicates that the TS during the recent decade was potentially less predictable. In addition, the persistent skill of the TS from October to DJF of the 21-year period is also reduced from 0.24 to 0.12 in the recent 11 years. Considering the lower persistence and potential predictability during the recent decade, the independent statistical forecast skill appears to be valid and useful. The interannual variability is less predicted than the interdecadal variability and warming trend by the current statistical forecast. The statistical forecast shows difficulty in capturing the positive phase (warming over the northern AWM region) but well predicts the negative phase of the second observed PC. The low performance of the ENSOrelated third PC may be attributed to the interdecadal change in the monsoon-ENSO relationship. The correlation coefficient between the third PC and the Niño 3.4 SST index dropped to 0.37 for the recent 21 years from 0. 56 for 1979/1980-1999/2000. This notable decadal change in the monsoon-ENSO relationship likely makes difficulty in statistical forecast.
It should be emphasized that the forecast skill demonstrated here for the AWM TS is a baseline and attained skill for both dynamical and statistical models. During the recent two decades, climate models have been improved remarkably; extensive effort will be devoted to further improvement. The statistical forecast can be further improved with the incorporation of better predictors than only the SO mean TS.