Was the extreme Northern Hemisphere greening in 2015 predictable?

The NDVI is calculated from contrasting surface reflectance in the near-infrared (ρ_nir) and red (ρ_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 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).

The climate variables used in this study were obtained from ERA-Interim reanalysis [28] (http://apps.ecmwf.int/datasets/) as monthly means of daily values for the period 2000–2015, at 0.75°lat/lon resolution, for latitudes above 30°N. Seasonal anomalies (Jun-Sep) were calculated with reference to the period 2000–2014 for: 2 m air temperature (T, °C), average daily precipitation (P, mm.d⁻¹), photosynthetically active radiation (PAR, %), and volumetric soil water content (depth 0−2.89 m, SWC, %). Climate anomalies in 2015 are shown in figure 2.

Zoom In Zoom Out Reset image size

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 (http://www.ncdc.noaa.gov/teleconnections/). The Multivariate ENSO Index (MEI) [33], NAO and EA indices are provided by NOAA Climate Prediction Centre http://www.cpc.ncep.noaa.gov/).

Monthly values of the five indices between December 1999 and December 2015 were selected.

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).

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).

Zoom In Zoom Out Reset image size

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).

Zoom In Zoom Out Reset image size

To understand the poor skill of climate-based prediction in the pixels with residuals greater than ±1σ_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 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://globalfiredata.org. Monthly data were re-gridded to the common 0.75°lat/lon resolution and aggregated for January–September.

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 2σ_ref (figure 1(d)). The strongest pixel-scale anomalies (extending to ±3σ_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.

The leading modes of growing season temperature and precipitation were found to be the best climate-related 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² = 0.63 figure S6(a)), predicts the strong hemispheric greening in 2015 and its spatial pattern, with most residuals smaller than ±0.5σ_ref, but predominantly negative (figures S6(b) and (c)), of −0.15σ_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.02σ_ref). Most of the predicted pixel-scale anomalies (figure 3(c)) remain within ±1σ_ref of observations (75%), although some pixels still present errors above ±2σ_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.

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 1σ_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).

High residuals (above ±1σ_ref) are dominated by negative values (i.e. 2015 NDVI is underestimated) in pixels where 2015 NDVI was exceptionally high (2σ_ref or more) and mostly explained by IAV (figure S3). We analyse the occurrence of climate extremes or fire disturbance in 2015 as a possible explanation for poor-model skill in (figure 4(c), S9).

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.

Positive residuals (overestimates) are generally associated with extreme heat and dryness in water-limited regions, e.g. western USA and southern 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 1σ_ref (2σ_ref) are associated with burned areas.