Seasonal prediction of winter haze days in the north central North China Plain

Recently, the winter (December–February) haze pollution over the north central North China Plain (NCP) has become severe. By treating the year-to-year increment as the predictand, two new statistical schemes were established using the multiple linear regression (MLR) and the generalized additive model (GAM). By analyzing the associated increment of atmospheric circulation, seven leading predictors were selected to predict the upcoming winter haze days over the NCP (WHDNCP). After cross validation, the root mean square error and explained variance of the MLR (GAM) prediction model was 3.39 (3.38) and 53 % (54 %), respectively. For the final predicted WHDNCP, both of these models could capture the interannual and interdecadal trends and the extremums successfully. Independent prediction tests for 2014 and 2015 also confirmed the good predictive skill of the new schemes. The predicted bias of the MLR (GAM) prediction model in 2014 and 2015 was 0.09 (−0.07) and −3.33 (−1.01), respectively. Compared to the MLR model, the GAM model had a higher predictive skill in reproducing the rapid and continuous increase of WHDNCP after 2010.


Introduction
In recent years, the north central North China Plain (NCP; 34-43 • N, 114-120 • E) has suffered from increasingly severe winter (December-February) haze pollution (Ding and Liu, 2014), particularly after persistent heavy fog and haze in January 2013 (Zhang et al., 2014;Zhao et al., 2014).Af-ter 2000, the combined effects of a rapid increase in total energy consumption and the influence of climate change intensified the haze pollution in central north China (Wang and Chen, 2016).In conditions of heavy and slowly varying pollutant emissions, the fine particles in the atmosphere reach their saturation levels easily, and the climate conditions become another critical contributor of haze.Some new climatic studies should be helpful for diagnosing seasonal predictors of winter haze days over the NCP (WHD NCP ).The East Asian winter monsoon (EAWM) has a significantly negative relationship with WHD NCP (Yin et al., 2015a, b;Q. Li et al., 2015).By weakening EAWM circulations, negative sea surface temperature (SST) anomalies over the subtropical western Pacific could significantly intensify WHD NCP (Yin and Wang, 2016).Furthermore, the decline of preceding autumn (September-November) Arctic sea ice (ASI) has led to favorable environments for haze with stable atmosphere and greatly intensified haze pollution in eastern China (Wang et al., 2015).Although recent studies on the changes in WHD NCP and their associated mechanisms are new and still insufficient, they support the possibility of seasonal prediction.
The climate variables in east Asia showed obvious characteristics of tropospheric biennial oscillation, based on which a new interannual increment approach was applied for shortterm climate prediction (Wang et al., 2000(Wang et al., , 2012)).This new approach treated the year-to-year increment of a variable, i.e., the difference between the current and previous year (DY), as the predictand.Because the DY approach utilized the observed information from the previous year and the features of biennial oscillation, the interannual variation and interdecadal trend could be captured well.In addition, the signals (i.e., variance) of the predictors and predictand were both amplified (Huang et al., 2014) and thus of benefit to improve the prediction skill.If the predictive objects (Y ), e.g., haze days, were cross influenced by socioeconomic factors and climatic conditions, the predictand could be represented by Y = YS+YC, where YS and YC were the slowly varying socioeconomic and climatic components, respectively.
where the subscripts t and t − 1 indicate the current and previous years, respectively.
Commonly, the difference in pollutant emissions between current and previous year was very small, resulting in (YS t − YS t−1 ) ≈ 0, so DY ≈ (YC t − YC t−1 ).To some extent, the WHD NCP DY reflected the fluctuations caused by climate variability.After adding the predicted WHD NCP DY to the observed WHD NCP of the previous year, the interdecadal and socioeconomic components were contained in the final prediction.In prior studies, the DY approach has been used to explore the prediction of summer rainfall in China (Fan et al., 2008), heavy winter snow activity in northeast China (Fan and Tian, 2013), summer Asian-Pacific Oscillation (Huang et al., 2014) and winter North Atlantic Oscillation (Tian and Fan, 2015).Furthermore, some variables cross influenced by socioeconomic and climatic factors were predicted successfully using the DY approach, e.g., rice production in northeast China (Zhou and Wang, 2014) and the discoloration day for Cotinus coggygria leaves in Beijing (Yin et al., 2014).Considering the seriously negative impact of winter haze and the substantial need to predict WHD NCP , we made it the goal of this study to apply the DY approach to the seasonal prediction of WHD NCP .
The data and methods employed were introduced in Sect. 2. Section 3 described the predictors and associated circulations.We applied the DY approach to build the prediction models for WHD NCP in Sect. 4. In this section, the statistical models were built based on multiple linear regression (MLR) and a generalized additive model (GAM).Then, leave-one-out cross-validation and independent tests were performed to assess the statistical schemes of WHD NCP prediction.

Data sets and methods
Monthly atmospheric data, such as geopotential height and surface air temperature (SAT), were derived from the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) global reanalysis data set with a horizontal resolution of 2.5 • ×2.5 • from 1979 to 2016 (Kalnay et al., 1996).The monthly mean extended reconstructed SST data sets with a horizontal resolution of 2 • × 2 • from 1979 to 2016 were obtained from the National Oceanic and Atmospheric Administration (NOAA) (Smith et al., 2008).ASI extent was calculated from the ASI concentration data, downloaded from the Hadley Centre with a horizontal resolution of 1 • ×1 • from 1979 to 2016 (Rayner et al., 2003).The monthly gridded soil moisture data from 1979 to 2016 were downloaded from NOAA's Climate Prediction Center, with a horizontal resolution of 0.5 • × 0.5 • (Huug et al., 2003).The monthly Antarctic Oscillation (AAO) indices from 1979 to 2016 were also obtained from the Climate Prediction Center (Mo, 2000).
China ground observations from 39 NCP stations, collected by the National Meteorological Information Center of China four times per day from 1979 to 2016, were used to reconstruct the climatic WHD data (Yin and Wang, 2016).
Here, haze was defined as visibility less than a certain threshold and relative humidity less than 90 %.After excluding other weather phenomena affecting visibility, a day with haze at any time was defined as a haze day.Site WHD data were converted into grids after Cressman interpolation (Cressman, 1959), and then the WHD NCP was computed as the mean value of the gridded data.
In this study, the statistical models were built based on MLR and GAM methods.The MLR approach, a modeldriven method, was ultimately expressed as a linear combination of K predictors (x i ) that could generate the least error for prediction of ŷ (Wilks, 2011).With coefficients β i , intercept β 0 and residual ε, the MLR formula could be described as follows: (1) The GAM approach was more advanced and was developed from MLR and the generalized linear model (Hastie and Tibshirani, 1990).This method was particularly effective at handling the complex nonlinear and non-monotonous relationships between the predictand and the predictors, whose expressions were replaced by unspecified smooth functions (s).Similar to the generalized linear model, the dependent variable in GAM could have different probability distributions, such as Gaussian, Poisson and binomial, any of which could be transferred by the link function (g).The GAM was data driven rather than model driven.The resulting fitted values did not come from an a priori model that was adopted by the MLR and generalized linear model.The rationale behind fitting a nonparametric model was that the structure of data should be examined first to choose an appropriate smooth function for each predictor; i.e., the GAM allowed the data to determine the shape of the smooth function (Yee and Mitchell, 1991).The GAM could be written in the form (2) The normalized data sets from 1979 to 2013 were trained as the basic samples to fit the models, and those from 2014 to 2015 were treated as test data for independent prediction.Thereafter, the root mean standard error (RMSE), mean absolute error (MAE) and explained variance were calculated for evaluation by simple fitting and leave-one-out cross validation.

The predictors and associated circulations
To choose the DY predictors, the correlated DY atmospheric circulations were identified, as shown in Fig. 1.The positive phase of the east Atlantic/west Russia (EA/WR) and west Pacific (WP; Barnston and Livezey, 1987) patterns and the negative phase of the Eurasia (EU; Wallace and Gutzler, 1981) pattern were obvious, and we took the anticyclone circulation over north China as an intermediary that led to a more stable atmosphere to analyze the associated physical process.The positive anomaly over the NCP could confine the particles within a thinner boundary layer by suppressing vertical movement and, together with the cyclone, they could induce an easterly to weaken the East Asian jet stream (EAJS), producing weaker cold air.Meanwhile, the water vapor transportation was also enhanced by an anomalous southeaster in the lower troposphere (figure omitted), creating favorable conditions for more WHD NCP than in the previous year.
The pivotal local anticyclone over the NCP was the most important contributor; we therefore speculated that preautumn SAT DY around the NCP should be effective to impact WHD NCP DY.There were significantly negative correlations between WHD NCP DY and pre-autumn SAT DY from the Sea of Japan to the Stanovoy Range (35-65 • N, 130-140 • E), the area mean of which was selected as predictor x 1 (Fig. 2).The correlation coefficient (CC) between WHD NCP DY and predictor x 1 was −0.47, exceeding the 99 % confidence level.The features of negative EU and positive WP pattern could be identified clearly and the anomalous cyclone over south China and the South China Sea was significant in the circulations associated with predictor x 1 (× −1) (Fig. 3).Although the associated land-air interaction, especially in the DY field, was complicated and still unclear, according to the analysis of Fig. 1, the horizontal and vertical diffusion of pollutant particles would be restricted efficiently.
The pre-autumn SST anomalies and their associated winter SST of the Pacific could influence WHD NCP significantly via the air-sea interaction (Yin and Wang, 2016).Figure 4 shows the CC between WHD NCP DY and pre-autumn Pacific SST DY.The most significant CC distributed around the Gulf of Alaska (36-56 • N, 130-170 • W), and the area-averaged SST DY here was defined as predictor x 2 , whose CC with WHD NCP DY was 0.47 (above the 99 % confidence level).Chen and Wang (2015) found that the severe winter haze events in north China were closely related with the weaker and northward EAJS.The positive SST DY around the Gulf of Alaska could induce an obviously anomalous cyclone over eastern China and the adjacent ocean, and the stimulated easterly weakened the core of EAJS.Furthermore, there was a significantly anomalous southerly at the high latitude that restricted the cold activities from their source region and intensified the haze pollution over the NCP (Fig. 5).
Prior studies have documented that the triple SST pattern was a dominant mode of the northern Atlantic SST in autumn (Czaja and Frankignoul, 1999).When the pre-autumn SST anomalies were distributed in a "+ − +" pattern from south to north, the subsequent EAWM was stronger, and the surface temperature of north China was lower (Shi, 2014).Xiao et al. (2015) proved the SST anomalies over the North Atlantic from summer to the following winter exhibit a significant relationship with winter haze days on both decadal and interannual timescale.Similarly, the CC between WHD NCP DY and pre-autumn SST DY of the Atlantic was distributed  in a "− + −" pattern (Fig. 6).The area-averaged SST DY of the northern center was defined as predictor x 3 , whose CC with WHD NCP DY was −0.50, passing the 99 % confidence test.The most obvious DY atmospheric circulations related to predictor x 3 (× −1) were the positive WP pattern, whose south center linked with a subtropical high (Fig. 7).The continental high and marine low was both weakened by the anomalous geopotential height from the lower to the middle layer that led to weaker EAWM and weaker cold air.The pressure gradient over the east coast of China also resulted in significant southerly anomalies, indicating smaller surface wind and more moisture and resulting in more WHD NCP .
ASI decreased dramatically with significant variance and was a significant contributor influencing WHD in eastern China (Wang et al., 2015;Wang and Chen, 2016).The CC between pre-autumn ASI DY and WHD NCP DY was calculated (Fig. 8) and was significantly positive around the Beaufort Sea (73-78 • N, 130-165 • W).The area-averaged extent of ASI DY of the Beaufort Sea was selected as the fourth predictor (x 4 ), and its CC with WHD NCP DY was  0.37 (above a 95 % confidence level).A positive center of geopotential height at 500 hPa was located over the Central Siberian and Mongolian plateaus, and negative centers were distributed zonally from southern China to the subtropical Pacific (Fig. 9).Thus, the EAJS was weakened by the induced easterly and shifted northward, which illustrated less cold activities over the NCP (Yang et al., 2002) and generated more haze days.
Following SST, the soil moisture was another important factor for seasonal prediction (Guo et al., 2007).The WHD NCP was closely correlated with the moisture conditions due to the hygroscopicity of the atmospheric particles   (Yin et al., 2015a).Thus, the questions with respect to soil moisture were whether pre-summer or autumn soil moisture would be effective for seasonal prediction of WHD NCP DY.The area-averaged pre-autumn soil moisture DY of the Bohai Economic Rim (35-42 • N, 117-127 • E), defined as predictor x 5 , showed a significantly negative correlation with WHD NCP DY, i.e., the CC was −0.59, exceeding a 99 % confidence test (Fig. 10).The CC between predictor x 5 and geopotential height at 500 hPa was distributed in a similar way as in Fig. 1.The positive EA/WR and WP phases and the negative EU phase was obvious and led to more WHD NCP than in the previous year (Fig. 11).Being specific to local circulations, the cyclone over south China and the anticy-  The area-averaged soil moisture DY in the east of Mongolia was defined as predictor x 6 , whose CC with WHD NCP DY was 0.41 (above a 95 % confidence level).The negative EU pattern could be recognized from the associated atmospheric circulation with predictor x 6 (Fig. 13).The anomalous geopotential height was distributed zonally at high latitude indicating that the meridional circulations that transported cold air were weak.The positive high over the NCP could confine the vertical motion and the vertical diffusion  of atmospheric particles and intensify the haze pollution over the NCP.
Recently, some studies documented that AAO could affect the east Asian climate through cross-equatorial flow, e.g., the Somali jet stream (Fan and Wang, 2004, 2006, 2007a, b).After the late 1990s, global sea level pressure and geopotential height at 300 hPa in boreal January were characterized by the concurrence of the Aleutian Low and the negative phase of the AAO (F.Li et al., 2015).We investigated the relationship between WHD NCP DY and geopotential height at 850 hPa in the Southern Hemisphere and found that the distribution was remarkably similar to that of the negative phase of the AAO (Fig. 14).Furthermore, the CC between the September-October AAO DY and WHD NCP DY was −0.54, exceeding a 99 % confidence test.As shown in Fig. 15, the positive phases of the EA/WR and WP patterns were closely correlated with the negative phase of AAO and were responsible for more WHD NCP than in the previous year.The anomalous  anticyclone over the NCP and adjacent ocean not only led to stable atmosphere but also resulted in small wind and high humidity.Hence, the September-October mean AAO index was selected as the last predictor (x 7 ) to forecast the interannual increment of WHD NCP .

The prediction models and validations
In total, seven DY predictors (x 1 -x 7 ) were chosen to build the seasonal prediction model (SPM) for WHD NCP DY.Among the predictors were 21 types of pair combinations, of which only 5 pairs presented significant linear correlation.Thus, the multicollinearity would not be a problem when modeling with the MLR approach.Although the linear correlation between the predictand and each predictor was significant, the nonlinear interaction would also affect the WHD NCP and should be taken into account.In this section, seasonal prediction models were established using MLR (SPM MLR ) and GAM (SPM GAM ) and validated in detail.
The WHD NCP DY showed obvious features of biennial oscillation (Fig. 16), illustrating the DY approach was suitable for its prediction.The SPM MLR of WHD NCP DY was as follows: DY×10 = −2.774x 1 +2.582x 2 −1.631x 3 +2.528x 4 − 2.229x 5 + 2.555x 6 − 1.812x 7 .After leave-one-out cross validation, the RMSE CV of SPM MLR was 3.39 days, and the CC between fitted and observed WHD NCP DY was 0.73, accounting for 53 % of the total variance (Table 1).The percentage of the same sign (i.e., meaning the mathematical sign of the fitted and observed WHD NCP DY was the same) was 79.4 %.The SPM MLR showed good ability to predict the negative and minimum WHD NCP DY but did not adequately capture the continuous positive value after 2011 (Fig. 16a).The fitted WHD NCP DY from 2011 to 2013 varied similarly with the one before 2010 and did not reflect the rapid rising trend after 2010.As an independent prediction test, the predicted bias, i.e., the predicted value minus the measurement, in 2014 was 0.09, illustrating good performance, but the bias in 2015 was larger, i.e., −3.33.
We also applied the GAM approach to build a prediction model that would contain the nonlinear relationship with smooth functions.The SPM GAM of WHD NCP DY was as follows: DY × 10 = −2.164s(x 1 ) + 2.036s(x 2 ) − 1.721x 3 +2.588s(x 4 )−2.157s(x 5 )+2.187x 6 −2.506x 7 .During the simple fitting, the SPM GAM performed very well.The RMSE was 1.56 days, and the CC between the fitted and observed WHD NCP DY was 0.95.The SPM GAM could fit the minimum (in 2003) and maximum (in 2013), and show the trend well, indicating an advantage to process the nonlinear relationship.After cross validation, the performance of SPM GAM decreased dramatically, meaning that its stability was worse than that of SPM MLR .The RMSE CV of SPM GAM was 3.38 days and the CC between fitted and observed WHD NCP DY was 0.74, accounting for 54 % of the total variance (Table 1).The percentage of the same sign of SPM GAM results was 73.5 %, which was close to the re-  sult from SPM MLR .The SPM GAM also showed good ability to predict the negative and minimum WHD NCP DY and better performance to fit the maximum in 2013 (Fig. 16b).The predicted bias in 2014 and 2015 was −0.07 and −1.01, respectively, and the results were slightly better than those from SPM MLR .The CC between the bias of SPM MLR and SPM GAM from 1980 to 2013 was 0.83 (above a 99.99 % confidence level).If the SPM MLR performed well in some years, the SPM GAM also showed good ability in these years, and vice versa.We speculated that the reason was that some useful factors were not diagnosed and included here.
After adding the predicted WHD NCP DY to the observed information in the previous year, the predicted WHD NCP in the current year was obtained.For example, the predicted WHD NCP DY in 2012 was added to the measured WHD NCP in 2011, and the result was the final predicted WHD NCP in 2012.In Fig. 17, the simulated WHD NCP anomaly was fitted by cross validation from 1980 to 2013 and predicted in 2014 and 2015.For SPM MLR and SPM GAM , the CC between the original (detrended) observed and simulative WHD NCP was 0.89 (0.87) and 0.90 (0.88), respectively.Both of these September-October AAO index DY −0.54 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q −10 prediction models could capture the interannual and interdecadal trend and the extremums.The percentage of the same sign of the anomalies from the two models was 100 %, meaning these two models could predict the sign of the WHD NCP anomaly successfully.The SPM GAM could simulate the abrupt rising trend in 2010 better than SPM MLR , which was important for the prediction of recent years.

Conclusions and discussions
In this paper, we treated the WHD NCP DY as the predictand and built two prediction models using the MLR and GAM approach.In the DY atmospheric circulation, the positive phases of the EA/WR and WP patterns and the negative phase of the EU pattern intensified the haze pollution by inducing positive anomalies over the NCP and the Sea of Japan.Finally, seven leading predictors were selected and were listed in Table 2.
After cross validation, the RMSE CV and explained variance of SPM MLR (SPM GAM ) was 3.39 (3.38) and 53 % (54 %).The percentage of the same sign of these two prediction models was also similar, i.e., more than 73 %.The WHD NCP DY increased rapidly and persistently after 2010, and the SPM GAM could capture this trend better.For the final predicted WHD NCP , both of these two prediction models could capture the interannual and interdecadal trends and the extremums.The percentage of the same sign of the anomalies from two models was 100 %, and the SPM GAM simulated the abrupt increase in 2010 better than SPM MLR .The predicted bias of SPM MLR (SPM GAM ) in 2014 and 2015 was 0.09 (−0.07) and −3.33 (−1.01), respectively.Both of these models performed well in the independent tests, but the biases of SPM GAM were slightly smaller.The consistency of these two models might indicate that, after including plentiful predictors, the linear relationship dominated the WHD NCP DY prediction.Actually, the studies about the associated physical mechanism, i.e., how the external forcings influenced haze pollution, were new and still insufficient.In this paper, the underlying physical process was presented mostly from the way that the associated circulations impacted the WHD NCP DY.Thus, the physical mechanisms via which the external forcings stimulated such anomalous circulations were worthy of further study.
Although these two statistical models performed well during most of the past 3 decades and could predict the WHD NCP in 2014 and 2015 with small biases, they showed disadvantages when simulating the rapid rising trend after 2010.The large abrupt change was a common challenge to the statistical models, including the DY approach, so the numerical model should be introduced into the prediction of haze pollution.At the same time, if the SPM MLR performed well in some years, the SPM GAM also showed good ability in these years, and vice versa.One possible reason could be that some useful factors, most notably the human activities, were not included here.There was no doubt that the human activities, especially the energy consumption, were the fundamental drivers for the increase of haze pollution.In this paper, we simply assumed that the difference in pollutant emissions between current and previous years was very small and that the socioeconomic component of WHD NCP varied slowly.This assumption could support the seasonal prediction of haze days in most of the years but was still a compromise.In certain years, especially the recent years, this pollutant emission proportion varied rapidly, which needed to be taken into account.The preceding autumn energy consumption should be a good choice but is difficult to measure, and its DY could be introduced into the developed models directly to improve the predictive skill.

Figure 1 .
Figure 1.The CC between WHD NCP DY and geopotential height at 500 hPa (Z500) in winter from 1980 to 2013.The white curves indicate that the CC exceeded the 95 % confidence level.A and C represent anticyclone and cyclone, respectively.

Figure 2 .
Figure 2. The CC between WHD NCP DY and SAT DY in autumn from 1980 to 2013.The shades indicate that the CC exceeded the 95 % confidence level, and the rectangle represents the selected region (35-65 • N, 130-140 • E) of predictor x 1 .

Figure 3 .
Figure 3.The CC between predictor x 1 (× −1) and Z500 DY in winter from 1980 to 2013.The white curves indicate that the CC exceeded the 95 % confidence level.A and C represent anticyclone and cyclone, respectively.

Figure 4 .
Figure 4.The CC between WHD NCP DY and Pacific SST DY in autumn from 1980 to 2013.The shades indicate that the CC exceeded the 95 % confidence level, and the rectangle represents the selected region (36-56 • N, 130-170 • W) of predictor x 2 .

Figure 5 .
Figure 5.The CC between predictor x 2 and wind vector DY at 200 hPa in winter from 1980 to 2013.The shades indicate that the CC between the zonal wind DY and x 2 exceeded the 95 % confidence level.

Figure 6 .
Figure 6.The CC between WHD NCP DY and Atlantic SST DY in autumn from 1980 to 2013.The shades indicate that the CC exceeded the 95 % confidence level, and the rectangle represents the selected region (50-70 • N, 30-65 • W) of predictor x 3 .

Figure 7 .
Figure 7.The CC between predictor x 3 (× −1) and Z500 DY (shade)/850 hPa wind DY (arrows) in winter from 1980 to 2013.The dots indicate that the CC with meridional wind exceeded the 95 % confidence level.A and C represent anticyclone and cyclone, respectively.

Figure 8 .
Figure 8.The CC between WHD NCP DY and ASI DY in autumn from 1980 to 2013.The shades indicate that the CC exceeded the 95 % confidence level, and the box represents the selected region (73-78 • N, 130-165 • W) of predictor x 4 .

Figure 9 .
Figure 9.The CC between predictor x 4 and Z500 DY in winter from 1980 to 2013.The white curves indicate that the CC exceeded the 95 % confidence level.

Figure 10 .
Figure 10.The CC between WHD NCP DY and soil moisture DY in autumn from 1980 to 2013.The shades indicate that the CC exceeded the 95 % confidence level, and the rectangle represents the selected region (35-42 • N, 117-127 • E) of predictor x 5 .

Figure 11 .
Figure 11.The CC between predictor x 5 (× −1) and Z500 DY in winter from 1980 to 2013.The white curves indicate that the CC exceeded the 95 % confidence level.A and C represent anticyclone and cyclone, respectively.

Figure 12 .
Figure 12.The CC between WHD NCP DY and soil moisture DY in summer from 1980 to 2013.The shades indicate that the CC exceeded the 95 % confidence level, and the rectangle represents the selected region (48-52 • N, 115-125 • E) of predictor x 6 .

Figure 13 .
Figure 13.The CC between predictor x 6 and Z500 DY in winter from 1980 to 2013.The white curves indicate that the CC exceeded the 95 % confidence level.A and C represent anticyclone and cyclone, respectively.

Figure 14 .
Figure 14.The CC between WHD NCP DY and September-October Z850 DY from 1980 to 2013.The white curves indicate that the CC exceeded the 95 % confidence level.

Figure 15 .
Figure 15.The CC between predictor x 7 (× −1) and Z500 DY in winter from 1980 to 2013.The white curves indicate that the CC exceeded the 95 % confidence level.

Figure 16 .
Figure 16.The temporal variation of measured (black) WHD NCP DY, MLR (red, a) and GAM (red, b) cross-validation fitted WHD NCP DY from 1980 to 2013.The results for 2014 and 2015 represent the measured (black square) and predicted (red hollow circle) WHD NCP DY.

Figure 17 .
Figure 17.The temporal variation of measured (black) WHD NCP anomaly from 1980 to 2015, MLR (blue) and GAM (red) simulative WHD NCP anomaly, which was composed of cross fitted series from 1980 to 2013 and predicted values in 2014 and 2015.

Table 1 .
The RMSE, MAE, CC and explained variance (EV) of MLR and GAM models, and predicted bias for 2014 and 2015.The subscripts S and CV indicated simple and cross-validation fitting.

Table 2 .
The predictors and their meaning.CC indicates the correlation coefficient between predictor and WHD NCP DY from 1980 to 2013.