Climate Prediction of Satellite-Based Spring Eurasian Vegetation Index (NDVI) using Coupled Singular Value Decomposition (SVD) Patterns

: Satellite-based normalized di ﬀ erence vegetation index (NDVI) data are widely used for estimating vegetation greenness. Seasonal climate predictions of spring (April–May–June) NDVI over Eurasia are explored by applying the year-to-year increment approach. The prediction models were developed based on the coupled modes of singular value decomposition (SVD) analyses between Eurasian NDVI and climate factors. One synchronous predictor, the spring surface air temperature from the NCEP’s Climate Forecast System (SAT-CFS), and three previous-season predictors (winter (December–January–February) sea-ice cover over the Barents Sea (SICBS), winter sea surface temperature over the equatorial Paciﬁc (SSTP), and winter North Atlantic Oscillation (NAO) were chosen to develop four single-predictor schemes: the SAT-CFS scheme, SICBS scheme, SSTP scheme, and NAO scheme. Meanwhile, a statistical scheme that involves the three previous-season predictors (i.e., SICBS, SSTP, and NAO) and a hybrid scheme that includes all four predictors are also proposed. To evaluate the prediction skills of the schemes, one-year-out cross-validation and independent hindcast results are analyzed, revealing the hybrid scheme as having the best prediction skill. The results indicate that the temporal correlation coe ﬃ cients at 92% of grid points over Eurasia are signiﬁcant at the 5% signiﬁcance level in the hybrid scheme, which is the best among all the schemes. Furthermore, spatial correlation coe ﬃ cients (SCCs) of the six schemes are signiﬁcant at the 1% signiﬁcance level in most years during 1983–2015, with the averaged SCC of the hybrid scheme being the highest (0.60). The grid-averaged root-mean-square-error of the hybrid scheme is 0.04. By comparing the satellite-based NDVI value with the independent hindcast results during 2010–2015, it can be concluded that the hybrid scheme shows high prediction skill in terms of both the spatial pattern and the temporal variability of spring Eurasian NDVI.


Introduction
Unequivocal global warming since the mid-19th century has been revealed. In the Northern Hemisphere, the last three decades were the warmest of the last 1400 years [1]. The increased surface air temperature could prompt a longer growing season, increased photosynthetic activity and respiration, and an increase in productivity of terrestrial vegetation [2][3][4][5][6][7]. Meanwhile, terrestrial vegetation interacts with the climate through the surface albedo, hydrological processes, roughness, carbon exchange, and so on [8][9][10][11][12][13][14]. In Eurasia, the rising temperature has led to the northward extension of the forest treeline into tundra areas [15,16]. The changes in Eurasian vegetation not only have important effects on regional climate, such as the East Asian atmospheric circulation, Eurasian temperature and Remote Sens. 2019, 11, 2123 3 of 21 This study aims to develop prediction models for spring Eurasian NDVI based on current-season and previous-season climate factors. The models are expected to predict spring (April-May-June-averaged) NDVI over Eurasia (50 • -65 • N, 50 • -120 • E) in February. A new prediction method based on SVD paired modes of the predictand and the predictors was applied, which is expected to improve the prediction skill of spring Eurasian NDVI. The rest of this paper is organized as follows. In Section 2, the materials and methods used in this study are introduced. The predictors and their physical linkages to Eurasian NDVI are illustrated in Section 3. Section 4 reports the verification results of the prediction models. Finally, a conclusion is given in Section 5.

NDVI and Climate Data
The NDVI data used in this study were from the third-generation NDVI dataset (NDVI3g) derived from the advanced very high-resolution radiometer instruments onboard the National Oceanic and Atmospheric Administration (NOAA) series of meteorological satellites [37,66]. As chlorophyll absorbs red light and the mesophyll leaf structure scatters near-infrared light, the vegetation greenness, or vigor, can be represented by the combination of the earth surface reflectance of radiation at red and near-infrared wavelengths [67][68][69]. The NDVI was calculated based on the red spectral band (0.58-0.68 µm) and the near-infrared spectral band (0.725-1.10 µm) reflectance [NDVI = (NIR − RED)/(NIR + RED)]. The NDVI values of vegetation regions are greater than 0.1 [70,71]. The original global AVHRR data of the Global Area Coverage format at 4 km spatial and daily temporal resolution were processed (including sampling, calibration, noise removal, mapping, gap treatment, etc.) by the Global Inventory Modeling and Mapping Studies (GIMMS) archiving from 1981 through to the present day [72]. GIMMS produces the NDVI3g dataset at a 1/12 • × 1/12 • spatial resolution and 1/24 a year temporal resolution. Then, the NDVI3g dataset was reprojected into a monthly dataset at a 1.0 • × 1.0 • spatial resolution to suppress data noise.
The first spatial mode of the empirical orthogonal function (EOF) of Eurasian NDVI shows large values of approximately (50 • -65 • N, 50 • -120 • E) (Figure 1a), which means significant variations of NDVI. The time series of the first mode of EOF shows the characteristic of interannual variability ( Figure 1b). Therefore, our prediction of Eurasian NDVI (NDVIEA) was focused on the region of (50 • -65 • N, 50 • -120 • E). Figure 1c indicates that Eurasian NDVI experiences a rapid increase from April to June. Considering the importance of spring vegetation to the climate system, the focus is on the prediction of spring NDVIEA.
Monthly climate data were derived from the European Centre for Medium-Range Weather Forecast's interim reanalysis dataset (ERA-Interim) at the 0.75 • × 0.75 • spatial resolution during 1982-2015 [73], including the variables of monthly sea level pressure, zonal and meridional wind, geopotential height, 2 m air temperature, and soil temperature (level 1; 0.07 m). The monthly sea-ice cover and sea surface temperature data from the Met Office Hadley Centre observational datasets at the 1.0 • × 1.0 • spatial resolution during 1982-2015 were also used [74]. The satellite-derived monthly snow water equivalent data derived from the National Snow and Ice Data Center are available at 25 km Equal-Area Scalable Earth Grids from 1982 to 2007 [75]. The monthly snow cover data in the 89 × 89 cells Northern Hemisphere grid during 1982-2015 were derived from the Rutgers University Global Snow Lab and transformed to a 2.5 • × 2.5 • spatial resolution in this study [76]. The monthly NAO index and Niño3 index during 1982-2015 were derived from the publicly available NOAA-CPC database [77,78].
The spring monthly surface air temperature at a 1.0 • × 1.0 • spatial resolution during 1982-2015 were derived from CFSv2 [79]. CFSv2 makes retrospective forecasts using initial states from the Climate Forecast System Reanalysis from 1982 to 2010 and onward to calibrate subsequent operational real-time predictions. The nine-month forecast surface air temperature data initialized in March during 1982-2015 of April to June was used in this study.

Methods
The year-to-year increment approach was applied in the prediction and analysis processes: where ( ) is a variable in the current year; ( − 1) is the variable one year before ( ); and ( ) is the DY of the variable. On the one hand, the DY can amplify the interannual signal for both the predictand and the predictors compared to the anomalies of variables. On the other hand, the prediction based on the year-to-year increment approach can obtain the interdecadal variability of variables by employing ( − 1) [61]. Therefore, the year-to-year increment approach was applied to achieve higher prediction skill. The prediction approach based on the coupled spatiotemporal variability patterns was used. Firstly, SVD analyses are used to extract the pairs of spatial patterns of the predictand and predictor that are optimally correlated [58]: In the formulae, ( , ) and ( , ) are the DYs of the predictor and predictand, respectively; represents the number of all the SVD modes; ( ) and ( ) represent the singular vector of the th mode of the predictor and predictand, respectively; ( ) and ( ) indicate the corresponding time expansion coefficients of the th mode. Then, the predicted time expansion coefficients for each SVD mode of the predictand ( ) are obtained using linear regression of ( ) of the predictor:

Methods
The year-to-year increment approach was applied in the prediction and analysis processes: where V(t) is a variable in the current year; V(t − 1) is the variable one year before V(t); and V DY (t) is the DY of the variable. On the one hand, the DY can amplify the interannual signal for both the predictand and the predictors compared to the anomalies of variables. On the other hand, the prediction based on the year-to-year increment approach can obtain the interdecadal variability of variables by employing V(t − 1) [61]. Therefore, the year-to-year increment approach was applied to achieve higher prediction skill. The prediction approach based on the coupled spatiotemporal variability patterns was used. Firstly, SVD analyses are used to extract the pairs of spatial patterns of the predictand and predictor that are optimally correlated [58]: In the formulae, X(t, x) and Y(t, x) are the DYs of the predictor and predictand, respectively; m represents the number of all the SVD modes; U i (x) and R i (x) represent the singular vector of the ith mode of the predictor and predictand, respectively; A i (t) and B i (t) indicate the corresponding time expansion coefficients of the ith mode. Then, the predicted time expansion coefficients for each SVD mode of the predictand B i (t) are obtained using linear regression of A i (t) of the predictor: Here, β i is the linear regression coefficient of the ith mode, and ε i is the residual. Then, the predictand can be obtained by the predicted time expansion coefficients B i (t) and the singular vector: in which Y k (t, x) is the DY of the predicted variable based on the kth single predictor. In this study, more than one predictor was used in the statistical and hybrid models. The predictors are weighted using the multiple linear regression method: where Y DY (t, x) is the DY of the predicted variable based on multiple-predictor schemes; K is the total number of predictors; a k (t, x) is the multiple linear regression coefficients. The prediction process is based on a grid-to-grid scale, which means every grid carries a different weighting. Finally, the DY of the predicted variable is added to the observed variable in the previous year to obtain the predicted variable Y in the current year: in which Y Obs (t − 1, x) is the observed variable in the previous year and Y(t, x) is the final predicted variable in the current year.
To verify the prediction skill of the models, the one-year-out cross-validation method was used for the period 1982-2015 [80]. One year in the time span of the predictand was deleted, and the remaining years were used as the training period to develop the prediction model. The cross-validation was performed 34 times during 1982-2015, and the results were compared to the NDVI3g data to estimate the prediction skill. Meanwhile, an independent hindcast for the period 2010-2015 was also performed. The independent hindcast scheme was developed based on the previous years for the NDVI in the following year (i.e., the prediction in 2010 was based on the training period of 1982-2009; the prediction in 2015 was based on the training period of 1982-2014). This procedure was repeated six times.
The temporal and spatial correlation coefficients (TCCs and SCCs, respectively) between the cross-validated results and reanalysis data were calculated to evaluate the prediction skills. Furthermore, the root-mean-square-errors (RMSEs) are calculated: in which N is the number of years during the study period; y i and y oi represent the NDVI3g data and the prediction schemes' outputs, respectively. In the verification of the physical linkages between spring NDVIEA and climate factors, composite analyses were explored based on the criterion that the normalized DY of the indexes (including sea-ice index, Niño3 index and NAO index) were greater than 0.5 or less than −0.5. Note that all the composite analysis results were based on the variables in year-to-year-increment form. The linear correlations used in this study were in the form of Pearson correlation. Furthermore, the wave activity flux for diagnosing the propagation of Rossby waves was calculated using the following formula [81]: Remote Sens. 2019, 11, 2123 6 of 21 in which p = (pressure/1000hPa) and z = −H ln p (H is a constant scale height); φ and λ are the latitude and longitude; U = (U, V, 0) T represents the basic flow; a is Earth's radius; ψ represents the geostrophic streamfunction; N 2 is the buoyancy frequency squared; C U is a vector that indicates the phase propagation in the zonal direction; and M is the wave-activity pseudomomentum.

Results
In this section, the chosen predictors and the reasons for choosing them are introduced. The physical mechanisms of the relationships between spring Eurasian NDVI and the predictors are explained. Then, the prediction models applying the coupled SVD patterns based on the year-to-year increment method were developed for the spring Eurasian NDVI. The one-year-out cross-validation and independent hindcasts were performed to determine and compare the accuracy of the models.

Predictors
To develop an effective prediction model, the predictors should be physically significant and statistically significantly correlated with the predictand. Based on investigation and analyses of climate factors related to NDVIEA, one current-season predictor (spring surface air temperature) and three previous-season predictors (winter sea-ice cover, sea surface temperature, and NAO) were taken into consideration.

Spring Surface Air Temperature over Eurasia from CFSv2
Classical statistical seasonal climate prediction is based on the premise that the physical linkages between the predictand and the previous-season predictors remain stable in the coming prediction period. However, current-season climate factors may also closely relate to the predictand. If current-season climate factors available in a dynamic climate prediction system (e.g., CFSv2) can be used as synchronous predictors, the prediction skill could be improved. A prediction model that combines previous-season and current-season predictors is called a hybrid dynamical-statistical model. This a hybrid prediction approach has been used previously in the prediction of East Asian summer monsoon [55], winter temperature [60], spring dust weather frequency [57], and so on. Considering the good performance of the hybrid approach in seasonal climate prediction, this study attempted to include current-season predictors derived from CFSv2 in the prediction schemes.
Previous studies have found that spring NDVIEA is mainly associated with the condition of spring temperature and precipitation [49,82]. The growing season of Eurasian vegetation starts in March or April. The higher spring temperature is associated with an earlier and longer growing season [2,50,83]. Previous studies have also indicated that vegetation is more related to temperature than precipitation in northern middle and high latitudes [48,84]. Thus, this study attempted to use the spring surface air temperature as a predictor. As observational or reanalysis climate data for spring were unavailable in advance, CFSv2's real-time prediction results of spring surface air temperature data were used. Figure 2 shows the first paired mode of the SVD analysis of the DY of spring surface air temperature and NDVIEA. Note that the spring surface air temperature was derived from CFSv2 data initialized in March (over (40 The SAT-CFS field shows positive anomalies and the NDVIEA field shows positive anomalies over most of Eurasia, except for a negative center near (55 • -65 • N, 70 • -110 • E) (Figure 2a,b). The first paired mode explains 58.3% of the total variance, and the correlation coefficient between the time series of expansion coefficients is 0.65 (significant at the 1% significance level; Figure 2c). Meanwhile, the prediction skill of CFSv2 was tested ( Figure 3). The NDVIEA shows significant positive correlations with the Eurasian surface air temperature from ERA-Interim ( Figure 3a). Furthermore, the correlation coefficients between the surface air temperature of CFSv2 and ERA-Interim are significant over most parts of Eurasia ( Figure 3b). SAT-CFS not only describes important physical linkages to NDVIEA, but is also statistically significantly correlated with NDVIEA. Meanwhile, CFSv2 shows good real-time prediction skill for Eurasian spring surface air temperature. Therefore, SAT-CFS was considered as a predictor. time prediction skill for Eurasian spring surface air temperature. Therefore, SAT-CFS was considered as a predictor.   time prediction skill for Eurasian spring surface air temperature. Therefore, SAT-CFS was considered as a predictor.

Winter Sea-Ice Cover over the Barents Sea
The decrease in the Arctic sea-ice cover and its effects on weather and climate, especially over middle and high latitudes in the Northern Hemisphere, has become a hot topic [85][86][87][88]. Ji and Fan [53] found that the Arctic sea-ice cover is positively correlated with NDVIEA. Figure 4 depicts the leading paired mode of the SVD analysis of the DY of winter sea-ice cover over the Barents Sea, (78 • -85 • N, 25 • -130 • E), abbreviated as SICBS), and spring NDVIEA. The leading SVD mode shows that the positive SICBS anomalies and the overall positive NDVIEA anomalies are optimally correlated. The leading paired mode explains 63.3% of the total variance, and the correlation coefficient between the time series of expansion coefficients is 0.67 (significant at the 1% significance level; Figure 4c).
The possible mechanism underlying the linkage between SICBS and NDVIEA has been explored by Ji and Fan [53]. An increase in SICBS might induce a Rossby wave via anomalous adiabatic heating. This anomalous Rossby wave can propagate southeastward and southward, exciting anticyclonic circulation anomalies over Eurasia. Meanwhile, anomalous westerly winds bring warm and moist air from the Atlantic to Eurasia, helping to increase the winter surface air temperature south of 70 • N [53]. Therefore, the corresponding snow cover is less over Eurasia. These negative snow cover anomalies can persist from winter to spring to influence spring climate conditions that are relevant to NDVIEA. As less snow cover corresponds to a lower snow albedo, the earth surface could be heated owing to the receipt of more solar radiation. Meanwhile, less snow cover and less snowmelt might also influence surface temperatures via a reduction in latent heat [89][90][91]. Corresponding to an increase in winter SICBS, the spring soil temperature over Eurasia shows positive anomalies, which is favorable for higher spring NDVIEA. Therefore, Eurasian snow cover acts as a bridge between winter SICBS and spring NDVIEA, with increased winter SICBS potentially inducing increased spring NDVIEA via its positive effect on spring Eurasian soil temperatures. Thus, winter SICBS was used as a predictor of spring NDVIEA. The decrease in the Arctic sea-ice cover and its effects on weather and climate, especially over middle and high latitudes in the Northern Hemisphere, has become a hot topic [85][86][87][88]. Ji and Fan [53] found that the Arctic sea-ice cover is positively correlated with NDVIEA. Figure 4 depicts the leading paired mode of the SVD analysis of the DY of winter sea-ice cover over the Barents Sea, (78°-85°N, 25°-130°E), abbreviated as SICBS), and spring NDVIEA. The leading SVD mode shows that the positive SICBS anomalies and the overall positive NDVIEA anomalies are optimally correlated. The leading paired mode explains 63.3% of the total variance, and the correlation coefficient between the time series of expansion coefficients is 0.67 (significant at the 1% significance level; Figure 4c).
The possible mechanism underlying the linkage between SICBS and NDVIEA has been explored by Ji and Fan [53]. An increase in SICBS might induce a Rossby wave via anomalous adiabatic heating. This anomalous Rossby wave can propagate southeastward and southward, exciting anticyclonic circulation anomalies over Eurasia. Meanwhile, anomalous westerly winds bring warm and moist air from the Atlantic to Eurasia, helping to increase the winter surface air temperature south of 70°N [53]. Therefore, the corresponding snow cover is less over Eurasia. These negative snow cover anomalies can persist from winter to spring to influence spring climate conditions that are relevant to NDVIEA. As less snow cover corresponds to a lower snow albedo, the earth surface could be heated owing to the receipt of more solar radiation. Meanwhile, less snow cover and less snowmelt might also influence surface temperatures via a reduction in latent heat [89][90][91]. Corresponding to an increase in winter SICBS, the spring soil temperature over Eurasia shows positive anomalies, which is favorable for higher spring NDVIEA. Therefore, Eurasian snow cover acts as a bridge between winter SICBS and spring NDVIEA, with increased winter SICBS potentially inducing increased spring NDVIEA via its positive effect on spring Eurasian soil temperatures. Thus, winter SICBS was used as a predictor of spring NDVIEA.

Winter El Niño
El Niño is an important aspect of the climate system that can induce changes not only in tropical weather and climate, but also subtropical and high-latitude climate via atmospheric teleconnection [92][93][94][95]. According to Li et al. [96], spring vegetation greenness over Eurasia is closely linked to El Niño events. Figure 5 depicts the first paired mode of the SVD analysis of the sea surface temperature over the equatorial Pacific and NDVIEA. The leading SVD mode of sea surface temperature shows an El Niño-like pattern, with positive sea surface temperature anomalies over the eastern equatorial Pacific. Meanwhile, as shown in Figure 5b, negative spring NDVI anomalies east of 60 • E are associated with El Niño events. The variance explained by the first paired mode is 76.1%, and the correlation coefficient between their time series of expansion coefficients is 0.68 (significant at the 1% significance level). Therefore, El Niño events may have negative influences on spring NDVIEA.

Winter NAO
The NAO is a prominent atmospheric circulation pattern characterized by a large-scale seesaw pattern of sea level pressure between the Icelandic low and Azores high that dictates the climate variability of North America and Eurasia [100]. Li et al. [23] explored the linkage between late winter (January-February-March) NAO and spring Eurasia NDVI and found that positive-phase late winter NAO can lead to higher than average NDVIEA. In the first paired mode of the SVD analysis of sea level pressure (20°-80°N, 90°W-40°E) and NDVIEA, the sea level pressure shows a similar pattern to positive-phase NAO, and NDVIEA shows overall positive anomalies ( Figure 6). The variance explained by the first paired mode is 62.5%, and the correlation coefficient between their time series of expansion coefficients is 0.58 (significant at the 1% significance level). The NAO might therefore be an important factor of influence on NDVIEA.
In positive-phase winter NAO years, there are cyclonic circulation anomalies over the northern part of the North Atlantic and anticyclonic circulation anomalies over the southern part. In turn, due to the propagation of Rossby waves originating from the Novaya Zemlya, anticyclonic circulation anomalies are excited in a region centered near Lake Baikal [23]. Therefore, the corresponding The mechanism responsible for this linkage was investigated by Li et al. [96]. The mature phase of El Niño is in winter, when the eastern Pacific shows positive sea surface temperature anomalies and the western Pacific negative anomalies. In the following spring, much of the tropical Indian Ocean warms up according to the Indian Ocean capacitor effect proposed by Xie et al. [97]. Therefore, the anomalous easterly wind that blows from the cold waters of the western Pacific to the relatively warmer waters of the North Indian Ocean strengthens in spring [98]. In turn, the sea surface temperature over the North Indian Ocean cools in spring, and this cooling persists because of positive wind-evaporation-sea surface temperature feedback [99]. The anomalously suppressed convection that results from the cooling sea surface temperature over the North Indian Ocean excites Rossby waves and can affect the Asian climate [96]. Rossby waves originating from the Bay of Bengal propagate northeastward to the Kamchatka peninsula via southern China and Japan, which excites anomalous cyclonic circulation over a large region east to 80 • E [96]. This, in turn, affects spring NDVIEA via anomalous low surface air temperatures associated with the northeasterly and northwesterly winds in the west of the anomalous cyclonic circulation centered near the Kamchatka peninsula. As El Niño events indeed have negative influences on spring NDVIEA, the winter equatorial Pacific sea surface temperature (over (20 • S-20 • N, 140E • -70 • W), abbreviated as SSTP) was considered as a predictor.

Winter NAO
The NAO is a prominent atmospheric circulation pattern characterized by a large-scale seesaw pattern of sea level pressure between the Icelandic low and Azores high that dictates the climate variability of North America and Eurasia [100]. Li et al. [23] explored the linkage between late winter (January-February-March) NAO and spring Eurasia NDVI and found that positive-phase late winter NAO can lead to higher than average NDVIEA. In the first paired mode of the SVD analysis of sea level pressure (20 • -80 • N, 90 • W-40 • E) and NDVIEA, the sea level pressure shows a similar pattern to positive-phase NAO, and NDVIEA shows overall positive anomalies ( Figure 6). The variance explained by the first paired mode is 62.5%, and the correlation coefficient between their time series of expansion coefficients is 0.58 (significant at the 1% significance level). The NAO might therefore be an important factor of influence on NDVIEA.
Remote Sens. 2019, 11, x FOR PEER REVIEW 11 of 22 relative to spring NDVIEA. Meanwhile, the correlation coefficients between the time series of expansion coefficients corresponding to the first SVD mode between the predictand and predictors are significant at the 1% significance level. Therefore, prediction schemes were developed based on these predictors.

Cross-Validation
Four single-predictor hindcast schemes based on each predictor were explored (i.e., the SAT-CFS scheme, SICBS scheme, SSTP scheme, and NAO scheme). Meanwhile, the NDVIEA results from the three previous-season single-predictor schemes (i.e., SICBS, SSTP and NAO) were considered in a statistical scheme via multiple linear regression, in which different stations have different weights. In positive-phase winter NAO years, there are cyclonic circulation anomalies over the northern part of the North Atlantic and anticyclonic circulation anomalies over the southern part. In turn, due to the propagation of Rossby waves originating from the Novaya Zemlya, anticyclonic circulation anomalies are excited in a region centered near Lake Baikal [23]. Therefore, the corresponding anomalous westerly winds from the wet and warm Atlantic can lead to increases in surface air temperature and decreases in snow water over Eurasia. In the following spring, the corresponding thinner snow cover melts rapidly and Eurasian soil temperatures show positive anomalies because of the lower surface albedo and reduced latent heat from the snowmelt. Therefore, NDVIEA shows higher than normal values in spring following positive-phase winter NAO. Thus, the sea level pressure over (20 • -80  . Each of the four predictors represent a physical connection relative to spring NDVIEA. Meanwhile, the correlation coefficients between the time series of expansion coefficients corresponding to the first SVD mode between the predictand and predictors are significant at the 1% significance level. Therefore, prediction schemes were developed based on these predictors.

Cross-Validation
Four single-predictor hindcast schemes based on each predictor were explored (i.e., the SAT-CFS scheme, SICBS scheme, SSTP scheme, and NAO scheme). Meanwhile, the NDVIEA results from the three previous-season single-predictor schemes (i.e., SICBS, SSTP and NAO) were considered in a statistical scheme via multiple linear regression, in which different stations have different weights. Furthermore, a hybrid scheme was developed based on the one synchronous predictor scheme from CFSv2 and the three previous-season predictor schemes (i.e., SAT-CFS, SICBS, SSTP, and NAO).
At each grid point, the temporal correlation coefficients (TCCs) for the SAT-CFS scheme, SICBS scheme, SSTP scheme, and NAO scheme are shown in Figure 7. Of all the 1130 grid points, the TCCs between the cross-validation and satellite data of NDVI are statistically significant at the 5% significance level at 686, 608, 819 and 624 grid points, which account for 61%, 54%, 73% and 55% of the total grid points, respectively (Table 1). The averaged TCCs for the single-predictor schemes are 0.38, 0.36, 0.43 and 0.36, respectively, which are all significant at the 5% significance level. The TCCs between the NDVI3g data and the cross-validation results of the statistical scheme and the hybrid scheme are significant at the 5% significance level at 977 and 1033 grid points, which account for 87% and 92% of the total grid points, respectively (Table 1, Figure 7). The averaged TCCs for them are 0.51 (significant at the 1% significance level) and 0.55 (significant at the 0.1% significance level), respectively. In the hybrid scheme, the TCCs at almost all the grid points (96% of the grid points) are significant at the 10% significance level. Meanwhile, the averaged temporal RMSEs at each grid point were calculated (Table 1). The single-predictor schemes show relatively higher temporal RMSE (0.07, 0.08, 0.07 and 0.08, respectively), while the hybrid scheme and the statistical scheme show the least error (0.04, Table 1).
The time series of spatial correlation coefficients (SCCs) between the satellite NDVI3g data and the results of the schemes are also presented (Figure 8). For the single-predictor scheme, the SCCs in most of the years are significant at the 1% significance level (0.08) during 1983-2015, except for 1991 and 2015 in the NAO scheme and 1996 in the SICBS scheme. The SCCs of the statistical scheme and the hybrid scheme are significant at the 1% significance level in all the years throughout 1983 to 2015. The multi-year averaged SCC for each scheme is presented and the hybrid scheme shows the highest averaged SCC (0.60, Table 1). Therefore, the hybrid scheme, which combines the synchronous spring predictor with the previous-season predictors, shows the highest prediction skill in terms of the spatiotemporal distribution and amplitude of the temporal variability.

Independent Hindcast
Due to the good performance of the hybrid scheme in the cross-validation results, independent hindcasts for the period 2010-2015 were performed based on the hybrid scheme to further verify its prediction skill. Figure

Discussion
Vegetation plays an important role in the interaction between ecosystems and the atmosphere, and has gained much attention in recent years. To represent the vegetation dynamics at the global scale, the dynamic global vegetation models (DGVMs) have been developed [101]. Although these DGVMs are able to roughly simulate the distribution of global vegetation types and the equilibrium of ecosystems, it is hard for them to predict the evolution of terrestrial vegetation under the current global climate change situation [102,103]. Therefore, statistical climate predictions based on satellitebased vegetation data, applying local weather and climate predictors, have been developed [40][41][42]. For instance, Fu, et al. [41] developed a multiple linear regression model to predict NDVI over the Heihe River Basin in China based on local precipitation, temperature, and soil moisture in each grid cell. The model showed high prediction skill over the regions where the NDVI was highly correlated with local precipitation, temperature, and soil moisture. However, the prediction skill showed a need for improvement over the regions where the NDVI was not significantly related to the local predictors. Given that large-scale climate patterns provide a conceptual framework and a broader understanding of observed changes in the local physical environment, large-scale climate indices seem to be better predictors of ecological processes than local climate [45][46][47]. Therefore, This study focused on developing a climate prediction model based on the coupled spatiotemporal patterns between the Eurasian NDVI and large-scale as well as local climate factors. The coupled spatiotemporal patterns between NDVI and the different predictors enable the prediction models to capture the spatiotemporal variability of NDVI effectively, which improves the prediction skill

Discussion
Vegetation plays an important role in the interaction between ecosystems and the atmosphere, and has gained much attention in recent years. To represent the vegetation dynamics at the global scale, the dynamic global vegetation models (DGVMs) have been developed [101]. Although these DGVMs are able to roughly simulate the distribution of global vegetation types and the equilibrium of ecosystems, it is hard for them to predict the evolution of terrestrial vegetation under the current global climate change situation [102,103]. Therefore, statistical climate predictions based on satellite-based vegetation data, applying local weather and climate predictors, have been developed [40][41][42]. For instance, Fu, et al. [41] developed a multiple linear regression model to predict NDVI over the Heihe River Basin in China based on local precipitation, temperature, and soil moisture in each grid cell. The model showed high prediction skill over the regions where the NDVI was highly correlated with local precipitation, temperature, and soil moisture. However, the prediction skill showed a need for improvement over the regions where the NDVI was not significantly related to the local predictors. Given that large-scale climate patterns provide a conceptual framework and a broader understanding of observed changes in the local physical environment, large-scale climate indices seem to be better predictors of ecological processes than local climate [45][46][47]. Therefore, This study focused on developing a climate prediction model based on the coupled spatiotemporal patterns between the Eurasian NDVI and large-scale as well as local climate factors. The coupled spatiotemporal patterns between NDVI and the different predictors enable the prediction models to capture the spatiotemporal variability of NDVI effectively, which improves the prediction skill compared with multiple linear regression models, especially over the regions where the NDVI is not significantly related to the local predictors [59,[104][105][106][107].
In this work, cross-validation and independent hindcast results indicate that the year-to-year increment schemes are able to predict NDVIEA effectively, especially the hybrid model. The hybrid scheme performs better than the single-predictor schemes and the statistical scheme because of its comprehensive consideration of the previous-season signals of the sea and atmosphere and the synchronous atmospheric signal derived from the CFSv2 dynamic climate prediction model. The skillful prediction of NDVIEA uncovered in this study can be seen as an important reference for the evolution of vegetation dynamics, and is helpful for our understanding of the Eurasian climate and environment.
Nevertheless, there is still room for improvement in terms of the prediction skill over some regions where the TCCs are not significant. This improvement could be achieved in future work in several respects. Firstly, with the development of dynamic climate prediction models, it may be possible to incorporate more synchronous climatic predictors, such as precipitation and soil moisture, to further improve the prediction skill. Meanwhile, other processes that are not directly linked with climate can also influence the Eurasian NDVI, such as land use, land cover change, pests, and so on [108][109][110]. How to incorporate these factors into the NDVI prediction models needs to be investigated in future work.
Eurasia is the largest continent, covering more than 30% of Earth's total land area. This study focused on the region where the first spatial mode of the EOF pattern of NDVI shows the largest values in the middle to high latitudes, which means significant variations of NDVI. However, other regions in the middle to high latitudes over Eurasia may be affected by different climate factors. For example, Li et al. [96] found that the NDVI in Western Europe is strongly related to Central Pacific Niño events, whereas the NDVI in East Russia is significantly influenced by Eastern Pacific Niño events. Gong and Shi [111] found that large-scale climate signals could exert different influences on the interannual variations of Eurasian spring vegetation greenness from region to region. For example, for the NDVI over Eastern Europe, the Southern Oscillation dominates over other factors, whereas the Eurasian pattern is the most important atmospheric system for West Siberia. Therefore, future work is needed to expand the area of the prediction model.

Conclusions
Seasonal climate prediction schemes of spring vegetation over Eurasia were developed using the year-to-year increment approach. The prediction schemes are based on the coupled spatiotemporal patterns between NDVIEA and relevant climate factors. First, the local spring surface air temperature was found to be a dominant factor for NDVIEA. However, observational or reanalysis data for this factor are unavailable in advance. Thus, the surface air temperature derived from the CFSv2 dynamic climate prediction model was applied as a synchronous predictor (SAT-CFS). Meanwhile, the winter sea-ice cover over the Barents Sea (referred to as SICBS), the winter sea surface temperature over the equatorial Pacific (referred to as SSTP), and winter NAO, all have effects on NDVIEA via different teleconnections. Therefore, they too were considered, as three previous-season predictors.
Four single-predictor schemes (i.e., the SAT-CFS scheme, SICBS scheme, SSTP scheme, and NAO scheme) were proposed, as well as a statistical scheme comprising SICBS, SSTP and NAO, and a hybrid scheme that included all four predictors. The one-year-out cross-validation and independent hindcasts were performed to evaluate their prediction skills. For the cross-validation results, in general, all the schemes show high skill in predicting NDVIEA. For all six schemes, the TCCs of more than 50% of the grid points are significant at the 5% significance level. Meanwhile, the SCCs of the six schemes are significant at the 1% significance level in almost all years during 1983-2015. Among them, the hybrid scheme performs better than the other schemes, with higher TCCs (0.55) and SCCs (0.60). Meanwhile, the RMSE of the hybrid scheme is the least amongst the schemes (0.04). The independent hindcast of the hybrid scheme was performed during 2010-2015, and the results compared with the NDVI3g data. It can be concluded that the hybrid scheme can capture both the spatial and the temporal variability of NDVIEA.