# Was the extreme Northern Hemisphere greening in 2015 predictable?

Ana Bastos, Philippe Ciais, Taejin Park, Jakob Zscheischler, Chao Yue, Jonathan Barichivich, Ranga B Myneni, Shushi Peng, Shilong Piao, Zaichun Zhu

► To cite this version:

Ana Bastos, Philippe Ciais, Taejin Park, Jakob Zscheischler, Chao Yue, et al. Was the extreme Northern Hemisphere greening in 2015 predictable?. Environmental Research Letters, 2017, 12 (4), 10.1088/1748-9326/aa67b5. [https://doi.org/10.1088/1748-9326/aa67b5](https://doi.org/10.1088/1748-9326/aa67b5). [hal-01583964](https:// HAL Id: )

https://hal.science/hal-01583964

Submitted on 16 Sep 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
LETTER • OPEN ACCESS

# Was the extreme Northern Hemisphere greening in 2015 predictable?

To cite this article: Ana Bastos et al 2017 *Environ. Res. Lett.* **12** 044016

View the [article online](#) for updates and enhancements.

## Related content

- Atmospheric teleconnection influence on North American land surface phenology
  Matthew P Dannenberg, Erika K Wise, Mark Janko et al.

- Vegetation anomalies caused by antecedent precipitation in most of the world
  C Papagiannopoulou, D G Miralles, W A Dorigo et al.

- Satellite view of seasonal greenness trends and controls in South Asia
  Sangeeta Sarmah, Gensuo Jia and Anzhi Zhang

## Recent citations

- Youngwook Kim et al.

- Evaluation of the MODIS collections 5 and 6 for change analysis of vegetation and land surface temperature dynamics in North and South America
  Ethan Heck et al.

- Understanding the Land Carbon Cycle with Space Data: Current Status and Prospects
  Jean-François Exbrayat et al.
<figure><img src="image_2.png" /></figure>

Original content from
this work may be used
under the terms of the
Creative Commons
Attribution 3.0 licence.

Any further distribution
of this work must
maintain attribution to
the author(s) and the
title of the work, journal
citation and DOI.

<figure><img src="image_3.png" /></figure>

LETTER

# Was the extreme Northern Hemisphere greening in 2015 predictable?

Ana Bastos<fnref n="1" />,<fnref n="6" />, Philippe Ciais<fnref n="1" />, Taejin Park<fnref n="2" />, Jakob Zscheischler<fnref n="3" />, Chao Yue<fnref n="1" />, Jonathan Barichivich<fnref n="4" />, Ranga B Myneni<fnref n="2" />, Shushi Peng<fnref n="5" />, Shilong Piao<fnref n="5" /> and Zaichun Zhu<fnref n="5" />

E-mail: ana.bastos@lsce.ipsl.fr

Keywords: MODIS, NDVI, PDO, AMO, climate variability

Supplementary material for this article is available online

Abstract

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.

1. 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–3]. The trend in vegetation greenness has been linked to increasing growing season length at high

latitudes [1] and enhanced terrestrial CO₂ 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₂ fertilization, and nitrogen deposition [3, 7].

<footnote n="1">¹ Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="2">² Department of Earth and Environment, Boston University, Boston, Massachusetts 02215, United States of America</footnote>

<footnote n="3">³ Institute for Atmospheric and Climate Science, ETH Zurich, Zurich, Switzerland</footnote>

<footnote n="4">⁴ Instituto de Geografia, Pontificia Universidad Católica de Valparaíso, Valparaíso, Chile</footnote>

<footnote n="5">⁵ Sino-French Institute for Earth System Science, College of Urban and Environmental Sciences, Peking University, Beijing 100871, People's Republic of China</footnote>

<footnote n="6">⁶ Author to whom any correspondence should be addressed.</footnote>
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Northern Hemisphere 'greening' in 2015. (a) spatial distribution of growing season NDVI from MODIS Terra, calculated as the Jun–Sep standardized anomalies relative to 2000–2014 (NDVI in 2000–2014 standard-deviation units $\sigma_{ref}$). The regions where NDVI in 2015 ranks 1 st in the 2000–2015 period are highlighted with markers. (b) Pixel-scale distribution of NDVI anomalies in 2015 (green) compared with 2000–2014 (warm colours, one probability density function (PDF) per year, with light red corresponding to 2000 and dark red to 2014); (c) Time-series of observed NH NDVI (bold green line), and the time-series corresponding to the long-term trend (NDVI<sub>PC1</sub>, dashed green) and interannual variability (NDVI<sub>IAV</sub>, dotted green) components calculated from rPCA decomposition (SI). The resulting fit between 2000–2014 and corresponding 95% confidence interval (black line and grey interval respectively) calculated by the two MLRMs and the predicted anomaly in 2015 (dark grey line). (d), (e) the spatial distribution of NDVI<sub>PC1</sub> and NDVI<sub>IAV</sub> for 2015, respectively. All NDVI values are in $\sigma_{ref}$ units. For a comparison of NDVI anomalies from MODIS Terra and Aqua see figure S1 available at stacks.iop.org/ERL/12/044016/mmedia.</figcaption></figure>

Northern ecosystems also present strong interan-
nual variability (IAV), liked with climate variability
[8–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 predict-
able, the predictability of IAV is largely unknown.
Previous studies have linked IAV in vegetation activity
and the carbon balance of ecosystems with atmo-
spheric circulation variability [10, 14, 15]. In fact,
Hallet *et al* [16] proposed that climate indices relating
to atmospheric circulation patterns (or teleconnec-
tions) 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 Spec-
troradiometer (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 devia-
tions above the reference mean ($\sigma_{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
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 2.** Climate anomalies during 2015. (a) Temperature (T, in °C); (b) precipitation (P, in mm.d⁻¹) and (c) soil water content down to 3 m (SWC, in %). Climate variables are from the ERA-Interim Reanalysis, and anomalies are calculated as departure of the 2000–2014 average. The regions where 2015 anomalies were extreme (either ranking as the highest or lowest in the 2000–2015 period) are highlighted with markers.</figcaption></figure>

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₂ 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 climate-change 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.

## 2. Methods

### 2.1. Vegetation greenness

The NDVI is calculated from contrasting surface reflectance in the near-infrared ($\rho_{nir}$) and red ($\rho_{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
<figure><img src="image_8.png" /></figure>

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

## 2.2. Climate data

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.

## 2.3. 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 Pacific 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 Surface 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/teleconnections/). 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.

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

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

## 2.6. Disturbance and extremes

To understand the poor skill of climate-based prediction in the pixels with residuals greater than $\pm 1\sigma_{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
<figure><img src="image_9.png" /></figure>

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.

# 3. Results

## 3.1. 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 $2\sigma_{ref}$ (figure 1(d)). The strongest pixel-scale anomalies (extending to $\pm 3\sigma_{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.

## 3.2. MLR model fit

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 MLR model 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 $\pm 0.5\sigma_{ref}$, but predominantly negative (figures S6(b) and (c)), of $-0.15\sigma_{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\sigma_{ref}$). Most of the predicted pixel-scale anomalies (figure 3(c)) remain within $\pm 1\sigma_{ref}$ of observations (75%), although some pixels still present errors above $\pm 2\sigma_{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.

## 3.3. 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 $1\sigma_{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).

## 3.4. Evaluating the model residuals

High residuals (above $\pm 1\sigma_{ref}$) are dominated by negative values (i.e. 2015 NDVI is underestimated) in pixels where 2015 NDVI was exceptionally high ($2\sigma_{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
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 3.** Predicting interannual variability in NDVI. (a) NH average $NDVI_{IAV}$ observed (green), corresponding MLRM fit for 2000–2014 (black) and predicted for 2015 (dashed) and 95% confidence interval (grey). The equation of the NH MLR fit ($NDVI_{IAV} = aPDO + bAMO$), $R^2$ and *p*-value for the 2000–2014 fit and the deviation between prediction and observation in 2015 are presented in the plot. (b) Time-series of AMO and PDO during DJF between 2000–2015, the best predictors of $NDVI_{IAV}$ in 2000–2014, and their corresponding correlation with $NDVI_{IAV}$; (c) Spatial distribution of the corresponding prediction, calculated from the coefficients *a* and *b*, calculated at pixel level and shown in figure S6.</figcaption></figure>

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\sigma_{ref}$ ($2\sigma_{ref}$) are associated with burned areas.

## 4. 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₂ 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 multi-annual 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 temperature-moisture 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)).
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 4.** Analysis of the prediction residuals. (a) Pixel-scale predicted NDVI for 2015, calculated as the sum of the two predictions for the interannual variability and long-term components (figure 3 and figure S5), to be compared with figures 1(a) and (b) the corresponding residuals (predicted minus observed). Negative (positive) residuals in green (purple) indicate that the observed NDVI in 2015 was higher (lower) than the climate-based prediction. (c) Analysis of disturbances or climate extremes in the regions where residuals exceed ±1σ<sub>ref</sub>: occurrence of fires during Jan–Sep 2015 (black squares), highest (upward triangles) or lowest (downward pointing triangles) values of temperature (red), precipitation (dark blue), soil moisture (cyan) and photosynthetically active radiation (yellow) in 2015. 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).</figcaption></figure>

The statistical approach used here is based solely on
climate variables (T and P for NDVI<sub>PC1</sub> and tele-
connection indices for NDVI<sub>IAV</sub>). 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<sub>2</sub>
fertilization or management practices [3], (ii) the
occurrence of climate extremes, explained neither by the
trends in T and P (predictors of NDVI<sub>PC1</sub>), nor by AMO
and PDO states (predictors of NDVI<sub>IAV</sub>); (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’ (climate-
related) 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 under-
estimated 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].

**5. 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 predict-
able 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
<figure><img src="image_14.png" /></figure>

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.

## Acknowledgments

This work is supported by the Commissariat à l'Énergie Atomique et aux Énergies Alternatives (CEA). This work was partly supported by the European Research Council Synergy grant ERC-2013-SyG-610028 IMBALANCE-P. This work was partially funded by the NASA Earth Science Division (Grant No. NNX14AP80A and NNX14AI71G). We thank Sonia Seneviratne and Ashley Ballantyne for helpful comments while preparing the manuscript.

## References

[1] Myneni R B, Keeling C, Tucker C, Asrar G and Nemani R 1997 Increased plant growth in the northern high latitudes from 1981 to 1991 *Nature* **386** 698–702

[2] Nemani R R *et al* 2003 Climate-driven increases in global terrestrial net primary production from 1982 to 1999 *Science* **300** 1560–63

[3] Zhu Z *et al* 2016 Greening of the earth and its drivers *Nat. Clim. Change* **6** 791–5

[4] Barichivich J *et al* 2013 Large-scale variations in the vegetation growing season and annual cycle of atmospheric CO₂ at high northern latitudes from 1950 to 2011 *Glob. Change Biol.* **19** 3167–83

[5] Peng S *et al* 2013 Asymmetric effects of daytime and nighttime warming on northern hemisphere vegetation *Nature* **501** 88–92

[6] de Jong R, Verbesselt J, Schaepman M E and de Bruin S 2012 Trend changes in global greening and browning: contribution of short-term trends to longer-term change *Glob. Change Biol.* **18** 642–55

[7] Los S 2013 Analysis of trends in fused AVHRR and MODIS NDVI data for 1982–2006: indication for a CO₂ fertilization effect in global vegetation *Glob. Biogeochem. Cycles* **27** 318–30

[8] Piao S *et al* 2014 Evidence for a weakening relationship between interannual temperature variability and northern vegetation activity *Nat. Commun.* **5** 5018

[9] Seddon A W, Macias-Fauria M, Long P R, Benz D and Willis K J 2016 Sensitivity of global terrestrial ecosystems to climate variability *Nature* **531** 229–32

[10] Bastos A *et al* 2016 European land CO₂ sink influenced by NAO and East-Atlantic pattern coupling *Nat. Commun.* **7** 10315

[11] Liu G, Liu H and Yin Y 2013 Global patterns of NDVI- indicated vegetation extremes and their sensitivity to climate extremes *Environ. Res. Lett.* **8** 025009

[12] Zscheischler J, Mahecha M D, Harmeling S and Reichstein M 2013 Detection and attribution of large spatiotemporal extreme events in earth observation data *Ecol. Inform.* **15** 66–73

[13] Luo Y, Keenan T F and Smith M 2015 Predictability of the terrestrial carbon cycle *Glob. Change Biol.* **21** 1737–51

[14] Gong D Y and Ho C H 2003 Detection of large-scale climate signals in spring vegetation index (normalized difference vegetation index) over the northern hemisphere *J. Geophys. Res.* **108** 4498

[15] Bastos A, Running S W, Gouveia C and Trigo R M 2013 The global NPP dependence on ENSO: La Niña and the extraordinary year of 2011 *J. Geophys. Res. Biogeosci.* **118** 1247–55

[16] Hallett T B *et al* 2004 Why large-scale climate indices seem to predict ecological processes better than local weather *Nature* **430** 71–5

[17] Pettorelli N *et al* 2005 Using the satellite-derived NDVI to assess ecological responses to environmental change *Trends Ecol. Evol.* **20** 503–10

[18] Gouveia C *et al* 2008 The North Atlantic Oscillation and European vegetation dynamics *Int. J. Climatol.* **28** 1835–47

[19] NOAA 2016 State of the Climate: Global Analysis for Annual 2015, (NOAA National Centers for Environmental Information), Technical report. Retrieved on June 9, 2016 from (www.ncdc.noaa.gov/sotc/global/201513)

[20] Orth R, Zscheischler J and Seneviratne S I 2016 Record dry summer in 2015 challenges precipitation projections in central europe *Nat. Sci. Rep.* **6** 28334

[21] Barichivich J *et al* 2014 Temperature and snow-mediated moisture controls of summer photosynthetic activity in northern terrestrial ecosystems between 1982 and 2011 *Remote Sens.* **6** 1390–431

[22] Tucker C J 1979 Red and photographic infrared linear combinations for monitoring vegetation *Remote Sens. Environ.* **8** 127–50

[23] Didan K 2015 MOD13C1 MODIS/Terra Vegetation Indices 16-Day L3 Global 0.05Deg CMG V006 NASA EOSDIS Land Processes DAAC (https://doi.org/10.5067/MODIS/MOD13C1.006)

[24] Lyapustin A *et al* 2014 Scientific impact of MODIS C5 calibration degradation and c6+ improvements *Atmos. Meas. Tech.* **7** 4353–65

[25] Vermote D 2013 Surface reflectance over Land in Collection 6. *MODIS Science Team Meeting*, 15–17 April 2013

[26] Didan K, Munoz B, Solano R and Huete A 2015 MODIS Vegetation Index Users Guide (MOD13 series). *Vegetation Index and Phenology Lab, The University of Arizona* pp 1–38

[27] WWW-MCD12C1 MCD12C1 Combined MODIS Land Cover Type Yearly L3 Global 0.05Deg CMG V051 NASA EOSDIS Land Processes DAAC (https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mcd12c1) (Accessed: 1 February 2016)

[28] Dee D P *et al* 2011 The Era-interim Reanalysis: configuration and performance of the data assimilation system *Q. J. R. Meteorol. Soc.* **137** 553–97

[29] Rasmusson E M and Wallace J M 1983 Meteorological aspects of the El Niño/Southern Oscillation *Science* **222** 1195–202

[30] Mantua N J, Hare S R, Zhang Y, Wallace J M and Francis R C 1997 A pacific interdecadal climate oscillation with impacts on salmon production *Bull. Amer. Meteor. Soc.* **78** 1069–79

[31] Knight J R, Folland C K and Scaife A A 2006 Climate impacts of the Atlantic multidecadal oscillation *Geophys. Res. Lett.* **33** L17706

[32] Leathers D J, Yarnal B and Palecki M A 1991 The pacific/ North American teleconnection pattern and United States climate. Part I: regional temperature and precipitation associations *J. Clim.* **4** 517–28

[33] Wolter K and Timlim M S 1998 Measuring the strength of ENSO events: how does 1997/98 rank? *Weather* **9** 315–24

[34] Seneviratne A *et al* 2010 Investigating soil moisture climate interactions in a changing climate: a review *Earth-Sci. Rev* **99** 125–61
<figure><img src="image_15.png" /></figure>

[35] Greve P et al 2014 Global assessment of trends in wetting and drying over land *Nat. Geosci.* **7** 716–21

[36] Randerson J, Chen Y, Werf G, Rogers B and Morton D 2012 Global burned area and biomass burning emissions from small fires *J. Geophys. Res. Biogeosci.* **117** G04012

[37] Di Lorenzo E and Mantua N 2016 Multi-year persistence of the 2014/15 North Pacific marine heatwave *Nat. Clim. Change* **6** 1042–47

[38] Seneviratne S I, Luthi D, Litschi M and Schar C 2006 Land-atmosphere coupling and climate change in Europe *Nature* **443** 205–9

[39] Wang S, Huang J, He Y and Guan Y 2014 Combined effects of the Pacific decadal oscillation and El Niño-southern oscillation on global land dry-wet changes *Nat. Sci. Rep.* **4** 6651

[40] Sharp E D, Sullivan P F, Steltzer H, Csank A Z and Welker J M 2013 Complex carbon cycle responses to multi-level warming and supplemental summer rain in the high arctic *Glob. Change Biol.* **19** 1780–92

[41] Myers-Smith I H et al 2015 Climate sensitivity of shrub growth across the tundra biome *Nat. Clim. Change* **5** 887–91

[42] Branstator G and Teng H 2012 Potential impact of initialization on decadal predictions as assessed for CMIP5 models *Geophys. Res. Lett.* **39** L12703

[43] Zhang R and Delworth T L 2007 Impact of the Atlantic multidecadal oscillation on North Pacific climate variability *Geophys. Res. Lett.* **34** L23708

[44] Collins M et al 2010 The impact of global warming on the tropical Pacific Ocean and El Niño *Nat. Geosci.* **3** 391–7