Was the extreme Northern Hemisphere greening in 2015 predictable?

The year 2015 was, at the time, the warmest since 1880, and many regions in the Northern Hemisphere (NH) registered record breaking annual temperatures. Simultaneously, a remarkable and widespread growing season greening was observed over most of the NH in the record from the Moderate Resolution Imaging Spectroradiometer (MODIS) normalized difference vegetation index (NDVI). While the response of vegetation to climate change (i.e. the long term trend) is assumed to be predictable, it is still unclear whether it is also possible to predict the interannual variability in vegetation activity. Here, we evaluate whether the unprecedented magnitude and extent of the greening observed in 2015 corresponds to an expected response to the 2015 climate anomaly, or to a change in the sensitivity of NH vegetation to climate. We decompose NDVI into the long-term and interannual variability components, and find that the Pacific Decadal Oscillation (PDO) and the Atlantic Multidecadal Oscillation (AMO) explain about half of NDVI interannual variability. This response is in addition to the long-term temperature and human-induced greening trend. We use a simple statistical approach to predict the NDVI anomaly in 2015, using the PDO and AMO states as predictors for interannual variability, and temperature and precipitation trends for the long-term component. We show that the 2015 anomaly can be predicted as an expected vegetation response to temperature and water-availability associated with the very strong state of the PDO in 2015. The link found between climate variability patterns and vegetation activity should contribute to increase the predictability of carbon-cycle processes at interannual time-scales, which may be relevant, for instance, for optimizing land-management strategies.


Introduction
The sustained increasing vegetation activity trend (greening) in the Northern Hemisphere (NH) has been a prominent feature in satellite observations since the 1980s and is consistently simulated by models [1][2][3]. The trend in vegetation greenness has been linked to increasing growing season length at high latitudes [1] and enhancemed terrestrial CO 2 uptake in northern ecosystems [2,4]. The greening pace has been associated with asymmetric effects of climate trends in vegetation activity [5] or variations in the climate forcing [6]. It has also been shown that regional greening trends are further attributed to land use change, land management, CO 2 fertilization, and nitrogen deposition [3,7].
Northern ecosystems also present strong interannual variability (IAV), liked with climate variability [8][9][10] and extremes [11,12], although the underlying mechanisms are less well understood [2,8]. Luo et al [13] have pointed that while the response of ecosystems to climate change is intrinsically predictable, the predictability of IAV is largely unknown. Previous studies have linked IAV in vegetation activity and the carbon balance of ecosystems with atmospheric circulation variability [10,14,15]. In fact, Hallet et al [16] proposed that climate indices relating to atmospheric circulation patterns (or teleconnections) may be better predictors of ecological variability, because they integrate the co-variation between the different drivers of ecological activity. Bastos et al [10] have described how changes in sea-level pressure in the North Atlantic associated with two such patterns influence heat and moisture advection towards Europe, ultimately driving variations in vegetation activity.
In 2015, the Moderate resolution Imaging Spectroradiometer (MODIS) shows widespread record greening during the growing season in NH mid-and high-latitudes (figure 1(a)). The normalized difference vegetation index (NDVI) has been proven to be a good surrogate of vegetation activity and is widely used in studies of vegetation phenology, productivity, and in disturbance monitoring [17,18]. Over the 16 yr period the mean NH NDVI has steadily increased, but the broadening of the distribution ( figure 1(b)) also indicates increasing spatial heterogeneity. NDVI in 2015 clearly stands out as an anomaly to the 2000-2014 distribution, extending to values 3-4 standard deviations above the reference mean (s ref ). At hemispheric scale, 2015 produced the largest absolute NDVI value on record, and ranked highest for 38% of the NH pixels (figures 1(a) and (c)).
While abnormally high growing season temperature was registered at hemispheric scale [19], climate anomalies in 2015 were spatially heterogeneous, with some regions at high-latitudes actually experiencing temperatures below the previous 15 yr average (figure 2). By contrast, southern and central Eurasia registered precipitation deficits combined with warming, resulting in strong dry conditions [20]. In spite of low summer precipitation, soil moisture was exceptionally high in central USA, mostly due to positive anomalies in rainfall and snow during previous seasons.
Some of the factors influencing the greening trend in the NH such as nitrogen deposition, CO 2 fertilization and land-use change [3], are unlikely to produce the rapid greening needed to explain the 2015 anomaly. The unprecedented magnitude and extent of the greening in 2015 may be due to (i) the acceleration of the climatechange related trend [6], (ii) a response to an extreme climate anomaly in 2015, or (iii) changes in the sensitivity of NH vegetation to the climate forcing [8].
Vegetation presents distinct responses to climate at different time-scales, and the drivers of the long-term greening and the interannual anomalies are not necessarily the same [5,7,8,21]. Therefore, we use principal component analysis (PCA) to decompose NDVI into its dominant spatio-temporal patterns: (a) a non-monotonic trend and (b) anomalies due to IAV and extremes. We then identify the corresponding climate-related drivers by evaluating the sensitivity of NDVI components to temperature, light and soil moisture, as well as to several teleconnection indices. We finally develop simple statistical models to evaluate the predictability of each NDVI component and of the 2015 anomaly.

Vegetation greenness
The NDVI is calculated from contrasting surface reflectance in the near-infrared (r nir ) and red (r red ) bands and is related to the amount of photosynthetically active radiation (400 to 700 nm) absorbed by vegetation [22].
Here we use the NDVI from the Terra MODIS Collection 6 standard product [23], spanning from 2000 to 2015 with 16 d temporal composite at 0.05°lat/lon. The Collection 6 is the latest version of the MODIS product, which has been improved in various ways from its predecessor (i.e. Collection 5): Level 1 B calibration, snow/cloud detection, aerosol retrieval/correction, polarization correction, [24,25]. Furthermore, a refined temporal compositing approach has enabled it to provide less biased radiometric observation of the Earth's surface for the last 16 yr [26]. All selected data are processed from MODIS land surface reflectance data and thoroughly corrected for atmospheric effects. To obtain valid observations for surface vegetation, we applied a strict quality control strategy to minimize non-vegetative signals (i.e. snow and clouds). Land cover was assessed by the MCD12C1 for the year 2007 (MODIS Land Cover Type: IGBP classes, available from 2001-2012 only) [27], and urban, permanent snow and ice, and sparsely vegetated or barren areas were excluded. Three quality control strategies with increasing strictness in order to further filter out low quality pixels were compared. Results were consistent between the quality control strategies; here we show the least conservative one as it keeps the largest number of pixels.
Data were first aggregated to monthly time step and growing season NDVI was calculated by averaging Jun-Sep observations. The growing season standardized NDVI anomaly (henceforth simply NDVI) is defined as the departure from the long-term Environ. Res. Lett. 12 (2017) 044016 climatology for 2000-2014; this procedure avoids the influence of the abnormal greening in 2015. The data were resampled to the coarser resolution of climate data (0.75°lat/lon and selected for the latitudes above 30°N. NDVI values were further compared with the ones derived from MODIS Aqua for the common period (2003-2015, figure S1).

Climate variability patterns
We select three main coupled atmosphere-ocean variability patterns influencing global climate anomalies: the El-Niño/Southern Oscillation (ENSO) [29], the Pacic Decadal Oscillation (PDO) [30] and the Atlantic Multidecadal Oscillation (AMO) [31]. We further analysed two patterns of atmospheric variability in the North Atlantic, influencing both Eurasia and North America: the North Atlantic Oscillation (NAO), the East-Atlantic Pattern (EA), which have been shown to modulate the interannual variability of ecosystem activity in Europe [10].
The PDO and AMO indices correspond to leading modes of Sea Surafce Temperature (SST) variability in the North Pacific [30] and the North Atlantic [32] oceans, respectively, and are provided by NOAA's Physical Sciences Division (www.ncdc.noaa.gov/tele connections/). The Multivariate ENSO Index (MEI) [33], NAO and EA indices are provided by NOAA Climate Prediction Centre (www.cpc.ncep.noaa.gov/).
Monthly values of the five indices between December 1999 and December 2015 were selected.

NDVI decomposition
In order to separate the long-term and interannual components of NDVI, we performed a rotated principal component analysis (rPCA) on NDVI. The rPCA is a multivariate statistical technique that allows decomposing a data set into a given number of uncorrelated components explaining most of the variance in decreasing order. The components explaining most of the variance are referred to as leading modes. This transformation allows the reconstruction of the original time-series from all, or only a given set of principal-components. More details are given in SI.
The leading mode presents high spatial (r ¼ 0.96) and temporal coherence (r ¼ 0.97) with the linear trend of NDVI, but it further captures year-to-year variations of the greening pace [3,6]. We reconstruct the long-term NDVI component (NDVI PC1 , i.e. 'trend') by retaining only the first component of NDVI and reversing the PCA transformation. The other components (2À16) were similarly used to reconstruct the IAV component (NDVI IAV ). The contribution of each component (trend and IAV) to NDVI in 2015 can then be calculated for each pixel (figure S3).

NDVI driving factors and prediction
In order to evaluate the drivers of NDVI PC1 (i.e. trend) we calculate its correlation with T, PAR and P, as well as with their principal-components (calculated as described for NDVI). NDVI PC1 correlates significantly with T and P and their leading modes, as well as one component (fifth) of PAR (figure S4). A set of multiple linear regression models (MLRM) using these variables as predictors of NDVI PC1 were trained at hemispheric and pixel scale during 2000-2014. The best and most parsimonious fit (according to Akaike's information criterion, AIC) was found for the MLRM using the leading modes of T (first) and P (first two) as predictors (figure S5), consistent with the stronger correlations found. The MLRM was trained at hemispheric and pixel-scale for 2000-2014 and used to predict NDVI PC1 in 2015 (figure S6).
We calculated the correlation between the different teleconnections in winter (DJF, as atmospheric circulation is stronger in winter) and during the growing season with the five leading NDVI components (54% of the variance) and NDVI IAV (table S1). All leading components except the first (trend) and NDVI IAV present significant correlation with at least one teleconnection, consistent with their prominent influence in global and hemispheric climate anomalies and consequent vegetation response [10,14]. A set of MLRMs were trained using combinations of T, P, PAR or, alternatively, the five teleconnection indices. The winter PDO and AMO were found to provide the best fit (lowest AIC). An MLRM using these indices as predictors was trained during 2000-2014 at hemispheric and pixel level (figures S7(a) and (b)) and then used to predict NDVI IAV ( figure 3).
The 2015 NDVI field was then reconstructed based on the predictions of the two MLRMs (figure 4 (a)) and the corresponding model residuals were evaluated for each pixel ( figure 4(b)). The same procedure was repeated for shorter training periods (figure S8).

Disturbance and extremes
To understand the poor skill of climate-based prediction in the pixels with residuals greater than ±1s ref , we investigate the occurrence of climate extremes and/or disturbances.
The occurrence of climate extremes is assessed separately for positive and negative residuals by evaluating the distribution of T, P, PAR and SWC Environ. Res. Lett. 12 (2017) 044016 anomalies ranking highest or lowest in 2015 (figure 4 (c) and figure S9). As it is convenient to distinguish regions with water-limited and energy-limited regimes [34], we perform the analysis separately for the arid, transitional and humid regions defined in [35].
The occurrence of fire was assessed by the GFED4s burned-area fraction product [36], downloaded at 0.25°lat/lon spatial resolution from http://globalfire data.org. Monthly data were re-gridded to the common 0.75°lat/lon resolution and aggregated for January-September.

NDVI decomposition
The leading mode of NDVI (figure S2) explains 22% of the variance and corresponds to a generalized and spatially coherent greening trend over most of the NH (figures 1(c) and (d)), excepting the region north of the Black and Caspian seas (browning trend). The other four leading components present variability at interannual to decadal time-scales (figure S2), and all present significant relationships with atmospheric circulation patterns (table S1).
The decomposition of NDVI into its long-term ( figure 1(d)) and interannual variability (figure 1(e)) components indicates different contributions of the trend and IAV to 2015 NDVI (figure S3). While NDVI PC1 explains most (83%) of the NH average NDVI anomaly in 2015 ( figure 1(a)), it explains only 57% of the pixel-scale anomalies and does not capture anomalies above 2s ref ( figure 1(d)). The strongest pixel-scale anomalies (extending to ± 3s ref , figure 1(e)) are mostly captured by NDVI IAV . NDVI PC1 explains most of the 2015 values in the North American and east-Asian tundra regions, as well as in eastern Europe, parts of central USA and east Asia (figure S3), regions dominated by managed ecosystems [3]. On the other hand, the biggest anomalies, such as the high NDVI observed in central Eurasia and southern USA and the negative anomalies (browning) in Europe and western North-America, are mostly explained by NDVI IAV . The cancelling regional greening and browning anomalies explains why NDVI IAV does not contribute to the hemispheric average NDVI as much as the trend.

MLR model fit
The leading modes of growing season temperature and precipitation were found to be the best climaterelated predictors of the spatio-temporal patterns of 2000-2014 NDVI PC1 , as they filter out anomalies due to natural climate variability (figures S4 and S5). The MLRM based on the leading components of T and P (explaining 22% and 27% of T and P variance, respectively) trained for 2000-2014 (R 2 ¼ 0.63 figure S6(a)), predicts the strong hemispheric greening in 2015 and its spatial pattern, with most residuals smaller than ± 0.5s ref , but predominantly negative (figures S6(b) and (c)), of À0.15s ref for the NH average.
Hemispheric NDVI IAV is predominantly high before 2008, followed by a period of lower but highly variable values ( figure 3(a)). The winter states of the PDO and the AMO ( figure 3(b)), explain about half of the variance in NDVI IAV between 2000-2014 and allow NH NDVI IAV in 2015 to be predicted with a very small difference (0.02s ref ). Most of the predicted pixel-scale anomalies ( figure 3(c)) remain within ±1s ref of observations (75%), although some pixels still present errors above ± 2s ref .
The coefficients of the model fit at pixel-scale (NDVI IAV sensitivity to PDO and AMO) are shown in figures. S7(a) and S7(b). The spatial pattern of NDVI sensitivity to AMO is consistent with the spatial configuration of the second leading mode of NDVI, adding robustness to the temporal correlations found (table S1). A consistent pattern is also found between pixel-scale NDVI IAV sensitivity to PDO and the second and third leading modes of NDVI. While the strong sensitivity of NDVI to PDO in southeastern USA is reflected in the second component, the browning in central Europe accompanied by greening in central Eurasia (with positive PDO) is more evident in the third.

Predicting 2015 NDVI anomalies
The reconstruction of NDVI in 2015 as the sum of predicted NDVI PC1 and NDVI IAV by the two MLRMs described above allows reproducing most of the spatio-temporal variability in NDVI over the 16 yr MODIS record and accurately predict the strong greening in 2015, with 74% of the residuals below 1s ref (figures 4(a) and (b)). The variables used as predictors of NDVI PC1 (leading modes of T and P) and NDVI IAV (PDO and AMO) appear to be robust, because they still present predictive skill when trained for shorter periods (figure S8). Negative residuals (underestimates of greening) are generally found in regions experiencing climate extremes favourable for vegetation growth: extreme wetness in arid and transitional regions which are water-limited [2,34,35] as, for example, in Turkey; and extreme warming and light in humid regions (energy-limited) as in far-east Russia.

Evaluating the model residuals
Positive residuals (overestimates) are generally associated with extreme heat and dryness in waterlimited regions, e.g. western USA and southern Environ. Res. Lett. 12 (2017) 044016 Eurasia. NDVI was also overestimated in areas experiencing conditions generally favourable for vegetation growth but associated with fire occurrence, e.g. energy-limited regions with extremely high T and PAR (e.g. eastern Europe) or water-limited regions registering extreme wetness (south of Caspian Sea). Overall, 69% (81%) of the pixels with residuals above 1s ref (2s ref ) are associated with burned areas.

Discussion
The leading modes of temperature and precipitation (associated with trends) generally provide a good fit for NDVI PC1 and allow capturing the main patterns of the trend-related 2015 NDVI. This is consistent with the predominance of climate-related trends in NH vegetation [1,2,3,5]. The magnitude of predicted NDVI PC1 is systematically underestimated, as the climate-based MLRM does not include the effect of CO 2 fertilization, nitrogen deposition or structural changes related with land-use change that further influence the greening trend.
The significant correlation found between the winter PDO and AMO and the components associated with IAV is likely due to their modulation of the multiannual variability pattern that dominates NDVI IAV during the 16 year period ( figure 3(a) and (b)), and consistent with the long time-scales of variability of these patterns. In winter, atmospheric circulation variability is generally stronger, and may influence vegetation in spring and summer through phenology and delayed physical effects (e.g. snow cover or soil moisture) [10,21].
The PDO is a decadal basin-wide Pacific SST variability pattern with spatial similarities to the ENSO pattern but with stronger expression in the North Pacific [30]. The winter of 2015 registered the highest PDO value (warmer than average SST along the North American western coast) in the 16 yr record ( figure 3(b)), likely amplified by the development of El Niño conditions in the autumn of 2014 [37]. The effect of the PDO in NDVI IAV may be generally linked to its influence on the temperaturemoisture feedback [38]: greening is associated with cooler and cloudier conditions and enhanced soil moisture in response to positive PDO (e.g. southeastern USA), while browning occurs in regions where positive PDO leads to warmer and drier growing season with increased radiation; these conditions increase soil water depletion (e.g. central Europe, part of the interior lowlands and western North America).
The AMO is a multi-decadal variability pattern in North Atlantic SST, likely linked to natural changes in the thermohaline circulation [31]. Winter AMO remained generally positive (warmer conditions) during the 16 year period, peaking in 2004-2007 and registering a low value in 2015 ( figure 3(b)). The AMO influences NDVI especially in higher latitudes, with greening response to high AMO in the eastern edge of Russia and Alaska linked to higher temperature and radiation (the latter mostly in Alaska). High AMO induces browning at the very high latitudes in eastern Eurasia and eastern Canada, due to cooling and cloudier conditions (figure S7). The opposite response is expected for low AMO.
The climate anomalies observed in 2015 (figure 2) are consistent with a very strong positive phase of the PDO with low AMO: the sharp drought/wet patterns and abnormal warming over most mid-latitude regions and central Eurasia; this explains the greening/browning patterns found in NDVI IAV (figure S7). The MLRM using PDO and AMO indices generally captures the patterns of NDVI IAV in 2015: the greening over central Eurasia and southern USA as well as the browning in central Europe and western North America ( figure 3(c)). The statistical approach used here is based solely on climate variables (T and P for NDVI PC1 and teleconnection indices for NDVI IAV ). The poor skill in predicting NDVI in 2015 may be due to: (i) non-climate related trends, associated for instance with land-use/ land-cover change, changes in nutrient availability, CO 2 fertilization or management practices [3], (ii) the occurrence of climate extremes, explained neither by the trends in Tand P (predictors of NDVI PC1 ), nor by AMO and PDO states (predictors of NDVI IAV ); (iii) the occurrence of disturbances during the growing season or during the preceding months. Since the highest residuals (figure 4(b)) are found in regions dominated by large IAV (figure S3), the processes described in (i) unlikely explain the high residuals in 2015.
In those regions where NDVI predicted by our linear model was too high, generally the 'potential' (climaterelated) greening was not observed either because of fire occurrence or leaf wilting in response to extreme heat and soil dryness. In regions where the model underestimated the greening, extremely favourable climate for vegetation growth was observed. These extremes may be partly linked with the influence of other teleconnections, in particular the 2014/15 El Niño which may have reinforced PDO anomalies [37,39].
Still, not all pixels with high residuals can be explained by climate extremes or fire, pointing to other processes not directly linked with climate (e.g. land use/land cover change, forest demography, pests). Also, in water-limited regions, vegetation may react more positively to wetness than predicted by a linear model, consistent with the strong and non-linear sensitivity of the predominant grassland and shrub ecosystems in these regions to climate variability, in particular under high soil moisture conditions [9,40,41].

Conclusion
Here we show that the record greening observed throughout the Northern Hemisphere in 2015 can be explained by the very strong state of the PDO embedded into a long-term greening trend, rather than changes in the greening pace [6] or in vegetation sensitivity to climate at the hemispheric-scale [8].
We further show that understanding the links between low-frequency climate variability patterns and NDVI allows predicting, to some extent, variability in vegetation activity at interannual time scales. The poor predictive skill found in some pixels is generally linked to extreme climate anomalies not controlled by PDO or AMO (used here to predict IAV) and the occurrence of fires and may further be due to other disturbances and possible non-linear responses that could not be captured with a simple linear model of IAV.
Since the AMO and PDO are potentially predictable a few years in advance [42], our results show that IAV in ecosystem activity might after all be potentially predictable for certain biomes [13]. The link found between AMO and PDO indices and growing season NDVI may be used to develop simple statistical models  The positive extremes of radiation and temperature are generally linked with reduced cloudiness and negative extremes in precipitation and soil moisture, which were not represented to improve readability. Arid (dark shade), transitional (light shade), and humid regions (white) as in [35] shown in (c).
Environ. Res. Lett. 12 (2017) 044016 to predict growing season vegetation productivity in the forthcoming seasons, although it might not perform so well for years with weaker PDO and likely needs to be re-evaluated for periods with negative AMO. Although simple, this approach may prove useful for governments, land managers and farmers.
However, it should be noted that the relationships found may not be stationary, as teleconnections may interact [10,39,43] or respond to anthropogenic influences [44]. Due to the long variability time-scales of these patterns (in the order of several decades), the full depiction of their influence on global ecosystems is still not possible. Yet, this may be feasible in the near future, as the length of the period of continuous global ecosystem monitoring increases.