Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Trends in Global Vegetation Activity and Climatic Drivers Indicate a Decoupled Response to Climate Change

  • Antonius G. T. Schut ,

    tom.schut@wur.nl

    Affiliations Plant Production Systems Group, Wageningen University, Droevendaalsesteeg 1, 6708 PB Wageningen, The Netherlands, PBL Assessment Agency for the Environment, Antonie van Leeuwenhoeklaan 9, 3721 MA Bilthoven, The Netherlands

  • Eva Ivits,

    Current address: European Environment Agency, Kongens Nytorv 6, 1050 Copenhagen, Denmark

    Affiliation Joint Research Centre, Via Enrico Fermi 2749, I - 21027 Ispra, Italy

  • Jacob G. Conijn,

    Affiliation Plant Research International, Wageningen UR, Droevendaalsesteeg 1, 6708 PB Wageningen, The Netherlands

  • Ben ten Brink,

    Affiliation PBL Assessment Agency for the Environment, Antonie van Leeuwenhoeklaan 9, 3721 MA Bilthoven, The Netherlands

  • Rasmus Fensholt

    Affiliation Section of Geography, Department of Geosciences and Natural Resource Management, Faculty of Science, University of Copenhagen, Oster Voldgade 10, 1350 Copenhagen, Denmark

Abstract

Detailed understanding of a possible decoupling between climatic drivers of plant productivity and the response of ecosystems vegetation is required. We compared trends in six NDVI metrics (1982–2010) derived from the GIMMS3g dataset with modelled biomass productivity and assessed uncertainty in trend estimates. Annual total biomass weight (TBW) was calculated with the LINPAC model. Trends were determined using a simple linear regression, a Thiel-Sen medium slope and a piecewise regression (PWR) with two segments. Values of NDVI metrics were related to Net Primary Production (MODIS-NPP) and TBW per biome and land-use type. The simple linear and Thiel-Sen trends did not differ much whereas PWR increased the fraction of explained variation, depending on the NDVI metric considered. A positive trend in TBW indicating more favorable climatic conditions was found for 24% of pixels on land, and for 5% a negative trend. A decoupled trend, indicating positive TBW trends and monotonic negative or segmented and negative NDVI trends, was observed for 17–36% of all productive areas depending on the NDVI metric used. For only 1–2% of all pixels in productive areas, a diverging and greening trend was found despite a strong negative trend in TBW. The choice of NDVI metric used strongly affected outcomes on regional scales and differences in the fraction of explained variation in MODIS-NPP between biomes were large, and a combination of NDVI metrics is recommended for global studies. We have found an increasing difference between trends in climatic drivers and observed NDVI for large parts of the globe. Our findings suggest that future scenarios must consider impacts of constraints on plant growth such as extremes in weather and nutrient availability to predict changes in NPP and CO2 sequestration capacity.

Introduction

Quantification of trends in ecosystem productivity is essential to understand observed and expected responses to environmental changes. Recent easing of climatic constraints have increased global vegetation productivity [1], although clear regional differences are present with both stronger positive but also negative changes in net primary production (see reproduced figure in Dent et al. [2]. For the northern latitudes, warming during the photosynthetic active period resulted in an increase in photosynthetic activity [3] and productivity [4]. Increasing atmospheric CO2 concentration is another contributing factor to increased productivity, estimated to be about equally important as climate [5]. Increases in productivity are limited in areas with nutrient [6] or water availability constraints and increasing drought frequency [7,8].

Studies on climate records covering long time-series with multiple decades indicate that trends are rarely monotonic, but include phases with different trends [911]. For vegetated areas that strongly respond to e.g. climatic drivers, short term fluctuation [12] and thus transitions in trends are to be expected in long time-series of vegetation indices, such as the Normalised Difference Vegetation Index (NDVI), as they are directly related to productivity [13,14]. Various authors have identified trend breaks in observed (wheat) yields at global and regional scales, with yields reaching a plateau in recent years after a long period of increasing yields [1519]. These transitions have been observed in GIMMS3g NDVI for the 1982–2008 period [20]. In all areas with a significant trend in vegetation activity (15% of the global land area), trend transitions were more common (74%) than monotonic trends (26%), trend breaks and a change from greening to browning occurred more frequently in the 2000s than in the 1980s [20].

Methodological choices affect conclusions about trends in specific regions. For example in the Sahel zone, the fraction of pixels with significant positive or negative trends depends on the method used to determine a seasonal or annual value [21] and the definition of the time-period that best reflects a season [22]. A wide range of choices have been used including: annual NDVI sums [23]; sums of months around the peak of season [22]; maximum LAI values [24]; residual trends after regression with rainfall amounts based on seasonal sums for dryland areas [21,25]; or integrals based on variable start and end points (either annual or growing season) determined from phenology [26,27] or soil thawing [3]. Also the selected trend detection method can make an important difference. Compared to a simple linear regression, fitting polynomials can accommodate gradual changes [10,28] under the assumption that residuals are randomly distributed and only include “white noise” [10]. The BFAST approach is specifically designed to determine discontinuities in remote sensing datasets [29].

Explanation of observed trends in NDVI datasets is notoriously difficult as trends in coarse spatial resolution datasets such as GIMMS3g-NDVI are difficult to verify [30]. Attempts to attribute observed trends derived from the GIMMS3g-NDVI dataset to e.g. land degradation resulted in controversy about methods and interpretation of trends. See for example comments of Wessels [31] on work of Bai et al. [23] and the response of Dent et al. [2]. Extreme events such as droughts clearly influence ecosystem responses [32], further complicating the interpretation. Only few studies validated NDVI-trends with field observations of changes in biomass, for obvious reasons. Dardel et al. [22] found that two-monthly mean NDVI were moderately strong related to herbaceous mass, although changes in vegetation composition can have a strong influence [33].

A better understanding of the influence of climatic, biophysical and human induced drivers of change is needed to better predict the vegetation response in scenarios for the future. The occurrence of a negative trend in NDVI time-series with a positive trend in climatic drivers indicate a potential decoupled vegetation response. A full understanding of this potential decoupling is essential. Firstly, it will affect our estimation of the amount of terrestrial carbon sequestration under future climate scenarios [34]. Secondly, it enables a better understanding of observed ecosystem productivity with methodologies based on vegetation indices [3537]. Thirdly, the human influence on ecosystems needs to be better understood to initiate relevant and targeted interventions. There is a considerable component of spatially correlated variation in NDVI datasets that is not related to temperature or rainfall, but that may be linked to land-use [28] or other human activity [38]. Although water availability is one of the major constraints in productivity of dryland areas, responses of vegetation indices to changes in rainfall are non-linear e.g. due to differences in rain use-efficiency [21,39]. This indicates that all climatic factors need to be weighted [40,41] or considered in combination, e.g. by using plant growth models.

This study aims to understand better a possible decoupling of observed changes in vegetation activity (as an approximation for biomass productivity) from an expected change derived from models responding to climatic drivers only. To this end, global trends in the GIMMS3g NDVI dataset covering the years 1982 to 2010 were compared with trends in modeled water-limited biomass production. The objective was to map regions where the vegetation response was decoupled from the climatic drivers and showed a negative trend in greening combined with a positive trend in modeled water-limited biomass production, indicating a decreasing productivity that may be associated with land use change or land degradation processes. To assess uncertainty due to methodological choices, trends in NDVI metrics derived from the GIMMS3g NDVI dataset (including seasonal sums, seasonal integrals and maximum NDVI), linear and segmented trend detection methods and correlations of spatial variation with MODIS-NPP and modeled water-limited biomass production were compared.

Materials and Methods

Remote sensing datasets

A 29 year GIMMS3g dataset derived from the series of AVHRR sensors was provided with bimonthly NDVI from 1982 up to and including 2010 [42]. The AVHRR channel 1 and 2 data, used to compute the GIMMS3g NDVI are calibrated as suggested by Vermote and Kaufman [43], and the derived NDVI is further adjusted using the technique of Los [44]. No atmospheric correction is applied to the GIMMS data except for periods (1982–1984 and 1991–1994) with stratospheric aerosols from volcanic origins [42]. A satellite orbital drift correction is performed [45], minimizing effects of orbital drift by removing common trends between time series of solar zenith angle and NDVI. The original spatial resolution of 8 km has been resampled to 5 arcminutes. Various global datasets with annual values derived from the GIMMS3g were compared: the sum of monthly means of NDVI on an annual (calendar) basis (AS-NDVI), between April and October (AO-NDVI) including the growing season in the Northern Hemisphere; between January and June (JJ-NDVI) reflecting the growing season of monsoonal climates in the southern hemisphere; between June and December (JD-NDVI) reflecting the growing season of monsoonal climates in the northern hemisphere; the large integral of the growing season (LI-NDVI); and the annual maximum NDVI (Max-NDVI) value derived with Phenolo software, for a full description see Ivits et al. [27].

The Large Integral was computed using a Savitsky–Golay filter available in the TIMESAT software [46]. These Savitsky-Golay functions were fitted to the GIMMS3g NDVI dataset using the following parameters: a seasonal parameter of 0.5; 2 envelope iterations; an adaptation strength of 1; a Savitzky–Golay window size of 2; and amplitude of the season start and season end of 20%.

The MOD17A3 dataset [47], based on the MODIS sensor, provides estimates of the annual net primary production (NPP) with an original spatial resolution of 1 km, here referred to as MODIS-NPP. MODIS-NPP is determined using NDVI and climatic data, relates reasonably well to flux-tower measurements and shows no overall bias [7,36,47]. MODIS-NPP grids for the period 2000–2010 were resampled to 5 arcminutes, exactly matching the spatial extent of grid cells in the GIMMS3g NDVI dataset.

Ecosystems dataset

Kier et al. [48] provided detailed maps of ecoregions and 14 biomes on a global scale. These biomes were overlayed with continents, resulting in 98 unique continent-biome combinations. Land-cover was derived from the GLC2000 dataset [49], including 18 land-cover classes.

Modelled total biomass weight

The influence of changes in daily radiation, temperature, humidity and rainfall on productivity was evaluated with LINPAC [50], a crop growth model based on the light-use-efficiency approach [51,52]. The model has been used to assess rainfed potential (i.e. water limited) biomass production at local to global scales. LINPAC computes the annual total biomass weight (TBW) produced for annual and perennial vegetation and is calibrated and validated using datasets of four perennial plant species from seven continents [53], and is also used in a comparative study for calculating maize production at four global sites [54]. The start and end of a growing season was determined by temperature and soil moisture conditions. In the calculations for this study the effects of changing CO2 concentrations during three decades up to 2010 have not been taken into account.

Climate data covering the years 1981–2010 were derived from the CRUTS3.10 monthly 30 arcminute resolution gridded dataset, available for each year [55,56], with weather variables recalculated to daily values as input for the crop model. Daily rainfall was computed by a random generator in combination with the number of rainy days per month to allow days without rain. For a more detailed description and references see [57]. Crop water availability is determined using a simple soil water module including one soil compartment, considering run-off, evapotranspiration and water percolation beyond rooting depth. Thickness of the soil compartment equals rooting depth up to a maximum as determined by soil limitations or maximum rooting depth of the plant species. Soil information was derived from the ISRIC-WISE database for soil properties [58], and the digital soil map of the world with a 5 arcminute resolution [59].

In this work, we quantify the influence of changes in climate on total biomass weight (TBW) of the vegetation, allowing some error in the absolute value for biomass production for a particular year and location. For each pixel and year, TBW for arable land use was determined using model parameters for a spring wheat (temperate) or a grain maize (tropical) crop. For other land uses the TBW was determined using Miscanthus (Miscanthus spp, a C4 species), with parameters adapted to allow growth in cold or short-season climates. A sensitivity analysis using Reed canary grass (Phalaris arundinacea) as C3 perennial species instead of Miscanthus showed only a minor effect on trend values [57]. The combined TBW including all land uses was calculated per pixel based on weights for the relative area of crop land and other land-use as determined by Erb et al. [60]. Ultimately, the global TBW map was produced at a spatial resolution of 5 arcminutes for the years 1981 up to and including 2010.

Comparison of trends

Trends were estimated with a simple linear regression (SLR), a piecewise regression (PWR) and a Thiel-Sen estimator for slopes, the latter based on the median of slopes from all pair-wise combinations. The Thiel-Sen medium trend estimator is particularly useful when there are extreme values in the data-series [61]. For this, the MATLAB® software package was used, and the function Thiel_Sen_Regress was downloaded from Matlab central (http://www.mathworks.se/matlabcentral/). In a broader context, sophisticated trend analyses have been proposed [911]. However, the assumption that residuals are randomly distributed, so-called “white noise”, needs to be assessed. Considering the large amount of pixels, a visual pixel-by-pixel evaluation of residuals as recommended is impractical. Also, detection of a large number of different trends in a time-series, as possible in e.g. the BFAST approach [29], each covering a limited number of years is not meaningful to understand responses to long-term processes as these are strongly influenced by cyclic weather patterns. Therefore, a trend estimate should preferably cover at least one full El Nino/La Nina cycle of 4–7 years to identify long term changes rather than responses to short-term climatic variations.

In a piece-wise regression, the number of segments in the PWR can vary, but for this study only two segments were used. This was done in order to restrict the influence of events with a relative short period of influence on landscape greenness. Therefore, each segment covered at least 6 full years (equivalent to one full El Nino/La Nina cycle), i.e. the year of the breakpoint was restricted to the interval between 1 January 1988 and 31 December 2004. Events that occur within a time-series can gradually change the NDVI or cause an abrupt change. For example widespread damaging fires, logging of trees and/or severe droughts all cause a sudden but strong change in annual NDVI values and a discontinuous time-series. A PWR can accommodate this by fitting segments that are discontinuous, i.e. end points of a fitted regression line within a segment are not required to be starting points of the next segment. In this work we are more interested in gradual changes associated with long-term processes. Therefore, trends in the PWR segments were allowed to change in direction but two sequential segments were forced to pass through the same breakpoint, i.e. creating a broken stick where the two parts are still connected. This means that a single event in a stable (e.g. evergreen forest) environment, such as a damaging widespread fire with a recovery within a few years, does not trigger the start of a new segment. Such an event will affect overall average greenness, but will not strongly affect the overall trend. However, an increasing frequency or fires affecting the ability of the landscape to recover to its previous greenness or a decreasing frequency allowing the landscape to better recover will have an effect on the overall trend.

The PWR breakpoint was estimated by splitting the time-series iteratively in two segments. This process was repeated using first relatively large and then refined steps to allow optimal positioning of the breakpoint within a year. Slope values of a segment was set to zero if not different from 0 (t-test, with p < 0.1). The optimum position of the breakpoint within the time series was determined by maximizing the coefficient of determination (CD) of the PWR including both segments. A segmented regression is only meaningful when truly better than simple linear regression as it needs to compensate for the loss of one degree of freedom. This means that for our shortest time-series (1983–2010) with 28 observations, the PWR-R2 needs to be at least 0.034 larger than the SLR-R2 to compensate, i.e. for a model with a SLR-R2 value of 0.1 the adjusted R2 for PWR is also larger than the adjusted R2 for SLR. To remain conservative, a pixel was therefore only counted as segmented when the CD was at least 0.1 better than the SLR-R2 value.

A composite map was created combining PWR trend estimates for segmented pixels and SLR trends for all other land pixels. When the CD of the PWR or the R2 of the SLR was below 0.1, the trend values of both segments were considered insignificant and set to 0. All trend values were converted to relative values by expressing them as a percentage of the mean to enable direct comparisons between areas with low and high NDVI values. Relative values were computed as either SLR-trends (and Thiel-Sen median trends) divided by the overall mean including all years or as PWR-trend for the first or second segment divided by the mean for that specific segment respectively.

Trend values are not reliable for areas with a large bare ground fraction or many years without vegetative growth. Others have e.g. used 0.02 standard deviations in NDVI to mask out these areas [21], based on the assumption that bare areas do not vary in annual NDVI. We found that productive areas with evergreen vegetation were also occasionally below this threshold. To this end, we have defined an upper NDVI threshold of 0.3 for bare areas based on the mean AS-NDVI over the full time-series and combined this with a threshold of 0.02 standard deviations to mask bare areas, this mask was used for the figure comparing trends and breakpoints for the three metrics only.

Comparison of annual measures of productivity

The NDVI metrics derived from the GIMMS3g dataset were related to the MODIS-NPP estimates. For each NDVI metric, the mean was determined per pixel for the years that overlapped with the MODIS-NPP dataset (2001–2010). Means for these same years were also determined for the MODIS-NPP dataset. Pixels without TBW (due to unproductive years) or NDVI (e.g. no season determined in TIMESAT) values in more than half of the years in the time-series were excluded. Simple second-order polynomial relationships between means of NDVI per pixel as independent and means of MODIS-NPP per pixel as dependent variables were established. Firstly, one single (global) relationship including all land pixels was fitted. Secondly, continent-biome stratifications were created and separate relationships were fitted for all land pixels within each continent-biome combination. All residuals and all fitted values from these continent-biome specific relationships were combined. Finally, the overall fraction of explained variation (R2) and root mean squared error (RMSE) values were computed using these combinations. The same was done for stratifications based on land-cover-continent combinations.

Using temporal averages of pixel values in a similar manner, NDVI metrics were also related to TBW values, however overlapping years in these datasets were 1982–2010. This provided statistics to quantify differences in NDVI metrics and influence of stratification (using biomes or land cover) on the accuracy of relationships between NDVI metrics and MODIS-NPP or TBW.

Determining decoupled vegetation response to climatic drivers of growth

In the comparisons between TBW and NDVI trends only productive areas are considered, defined here as the area including all pixels where a TBW value was calculated in at least half of the years in the time-series. Trends in NDVI metrics were considered decoupled from its climatic drivers when the direction of the trend in NDVI was negative, despite a positive TBW trend. Pixels with a positive NDVI trend and a neutral (including all pixels with a weak trend and SLR-R2 or PWR-CD < 0.1) or negative TBW trend are referred to as diverging, to differentiate these greening areas from browning areas. Pixels showing no or a very weak (SLR-R2 or PWR-CD<0.1) trend in NDVI were excluded and not counted as decoupled or diverging.

Trends in TBW were most often not segmented and (for reasons of simplicity) maps were created showing classes combining either positive or negative SLR-trends in TBW with monotonic (SLR, greening or browning) and segmented NDVI trends (PWR, greening or browning in both segments; greening followed by browning; or browning followed by greening) using the composite NDVI maps as explained above. Fractions of pixels in categories with various combinations of positive and negative TBW and these NDVI trends were determined.

For latitudinal comparisons of trends in NDVI and TBW, only productive pixels with a mean NDVI (1982–2010) above 0.1 were considered. For each pixel row, latitudinal means were determined within one degree latitude.

Results

Comparison of trends

Pixels with weak SLR-trends (R2 value < 0.3) occurred much more frequently than pixels with a moderately strong SLR-trend (R2 value > 0.6) (Fig 1). When compared to the SLR, the PWR increased the fraction of explained variation in the NDVI time-series, most notably in the lower R2-range, and PWR increased the frequency of pixels with a moderate to strong relationship. For more than 50% of the pixels in the AS-NDVI metric with weak or moderate SLR-trends (SLR-R2 value < 0.6), the fraction of explained variation improved with at least 0.1 when using the PWR (not shown). The frequency distributions of AS-NDVI and AO-NDVI were quite similar in each trend estimation method. The frequency of pixels with a very weak trend (SLR-R2 and PWR-CD < 0.1) was larger for the LI-NDVI metric when compared to the other two NDVI metrics shown in Fig 1.

thumbnail
Fig 1. Cumulative frequency distribution of the coefficient of determination for Simple Linear Regression (SLR) and piece-wise regression (PWR) trends in land pixels for the Annual Sum (AS-NDVI), Large Integral (LI-NDVI) and April-October (AO-NDVI) GIMMS3g NDVI metrics.

https://doi.org/10.1371/journal.pone.0138013.g001

The improvement of PWR over the SLR trend method varied between the NDVI metrics (evidenced by difference in the fraction of land pixels where the PWR-CD was at least 0.1 larger than the SLR-R2, Table A in S1 File), in particular south of 60°N where differences between LI-NDVI and the other NDVI metrics were largest (Fig 2). The fraction of segmented pixels varied with latitude but were lowest near the equator (between 20°S and 10°N) for all datasets (Fig 2).

thumbnail
Fig 2.

a, Latitudinal totals of land area and area with segmented trends for the Annual Sum (AS-NDVI), Large Integral (LI-NDVI) and April-October (AO-NDVI) GIMMS3g NDVI metrics. b, Fraction of land areas with segmented trends in these three datasets.

https://doi.org/10.1371/journal.pone.0138013.g002

The TBW dataset showed mostly monotonic trends with only for 6% of the pixels a segmented trend (Table A S1 File). The SLR and Thiel-Sen trends did not differ much across the NDVI metrics, with more or less comparable proportions of pixels with a positive or negative trend (Table B in S1 File). Most of the pixels did not show strong TBW trends, the R2 remained below 0.1 for 72% of the productive area (i.e. a neutral trend). From all pixels in productive areas, 23% showed a positive SLR-trend in TBW indicating more favorable climatic conditions and 5% a negative SLR-trend was found for only 5% (Table B in S1 File). For TBW, the PWR did not improve the fraction of explained variation strongly, only 5% of the TBW pixels showed a segmented TBW trend (Table B in S1 File).

Globally, the difference between NDVI metrics in the number of pixels with a positive or negative SLR and Thiel-Sen trend was small (except Max-NDVI). The differences were more pronounced when looking at pixels with a segmented PWR trend, although fractions of pixels for summed NDVI metrics were similar (e.g. 14–18% with monotonic negative trending pixels and 38–40% for monotonic positive trending pixels, Table B in S1 File).

The overall differences between NDVI metrics in the fraction of pixels within a specific SLR-R2 or PWR-CD range (see also Fig 1) were even more pronounced for biomes with evergreen vegetation, where the frequency of low CD values was much higher for LI-NDVI than for the other two metrics (compare legend numbers 1–4 in S1a–S1c Fig. a-). The differences between metrics resulted in distinct differences in estimated regional trend patterns (Fig 3). The number of unmasked pixels with a CD or R2 value above 0.1 were much more frequent for the AS-NDVI than AO-NDVI or LI-NDVI, in particular in the south-western areas of West Africa where evergreen vegetation occurs frequently. The percentage of annual change (e.g. Senegal) and the direction of the trends (e.g. western Mali) also differed between the NDVI metrics for particular areas. Spatially, there was more similarity and less scatter in the estimated year of the breakpoint for the AO-NDVI and AS-NDVI metrics than for the LI-NDVI metric (Fig 3).

thumbnail
Fig 3. Differences in the estimated year of the breakpoint (panels a, c, e) and most recent trends (panels b, d, f) between the AO-NDVI (panels a, b), AS-NDVI (panels c, d) and LI-NDVI (panels e, f) metrics.

Trends shown are the percentage change per year. Maps are combinations of SLR trends and the trend in the second segment, the PWR segment 2 trend was used when PWR explained at least 0.1 more variation than SLR. Pixels in bare areas or with a low fraction (< 0.1) of explained variation are shown as white.

https://doi.org/10.1371/journal.pone.0138013.g003

Comparison of annual measures of productivity

The spatial variation in mean NPP within the MOD17A3 dataset that can be explained varied between NDVI metrics (Table 1). Separate relationships for each biome or land-use type improved the fraction of explained variation and reduced RMSE values for all NDVI metrics (R2 of 0.72–0.80) and for TBW (R2 of 0.74–0.75) when compared with a single relationship (Table 1). Differences between land-use and biome specific relationships were relatively small. However, there were clear differences in the relationships across continent-biome combinations (Table 2) and land cover types. The differences in explained NPP variation between biomes or land cover types were large, e.g. for AS NDVI ranging from 22–74% with respect to biomes (Table 2) and from 12–63% for land-cover types (data not shown).

thumbnail
Table 1. Descriptive statistics for second order polynomial relationships between various metrics derived from GIMMS3G NDVI-metrics and NPP (g C m-2 yr-1 MOD17A3).

Comparisons of one generic and biome/land cover (LC) specific relationships per continent.

https://doi.org/10.1371/journal.pone.0138013.t001

thumbnail
Table 2. Means and ranges of R2 values across continents for (quadratic) relationships between NDVI and MODIS-NPP per biome.

Relationships were developed for each continent separately, shown here for the annual sum (AS), large seasonal integral (LI) and April-October (AO) NDVI metrics. For each biome, the overall best performing GIMMS3g NDVI metric (including January-June, JJ) is listed.

https://doi.org/10.1371/journal.pone.0138013.t002

On a global scale, spatial variation in NDVI correlated strongly to spatial variation in TBW (Table C in S1 File). Within biomes, between 46–67% of variation in NPP and 29–47% of variation in TBW was explained, when spatially combining the best NDVI metric per continent-biome (Table D in S1 File).

There were clear differences between biomes (Table 2) and continents (Table D in S1 File) in the explained fraction of variation when comparing continent-biome combinations. Relationships were in general weak for tropical and subtropical moist broadleaf forest and strongest for biomes in more Northern areas. Moderately strong relationships were found between NDVI metrics and MODIS-NPP for most biome-continent combinations, but the best performing NDVI metric varied per biome (Table 2).

Determining decoupled vegetation response to climatic determinants of growth

In total, 76% of all land pixels were productive, defined in this study as pixels where a TBW value was determined in more than half of the years. Between 17–36% of all pixels in these productive areas showed a decoupled vegetation response with a negative NDVI trend combined with a positive or neutral TBW trend (Table 3). However, this occurred on all continents in small areas, for example in north-east China (S2 Fig) or between the center and north of Argentina (S3 Fig). This area in Argentina borders to a wider region in South America with negative trends in NDVI metrics, although to the north of Argentina this negative trend is combined with a negative trend in TBW indicating worsening climatic conditions (e.g. Paraguay and the far north of Argentina). In Africa, improving climatic conditions with positive TBW trends combined with negative NDVI trends occurred in a relatively small area in the north of Morocco (S4 Fig) and in a wider area between 10 and 20°S, in particular in Angola, Tanzania, Zambia and Mozambique (S5 Fig), and on Madagascar.

thumbnail
Table 3. Proportions of productive pixels (with a TBW value calculated for more than half of the years) with decoupled, diverging or following or trends.

Decoupled trends combine a negative NDVI trend with a positive or neutral TBW trend, diverging trends combining a positive NDVI trend with a negative TBW trend, and pixels with following trends have NDVI trends in line with TBW trends. A proportion of 0.23, 0.72 and 0.05 of all productive pixels show positive, neutral or negative TBW trends respectively.

https://doi.org/10.1371/journal.pone.0138013.t003

From all productive pixels with a positive TBW trend (23%, thus excluding the pixels without a TBW trend) 17–33% were decoupled, equivalent to 4–8% of all pixels in productive areas (Table 3). This decoupled trend, mostly with a negative trend in recent times, occurred in wider areas in the Sahel zone and around 50 degrees latitude in Europe and Asia (S2 and S4 Figs, please note that the observed patterns differ between metrics). For the remaining pixels with a positive TBW trend, a fraction of 38–47% showed positive trends, equivalent to 9–11% of all productive pixels (Table 3). The areas showing a greening trend in line with climatic drivers after an declining period (see e.g. S3 and S5 Figs: central-south Angola, Ethiopia, Mexico and western edge of Brazil), may indicate some recovery.

From all productive pixels, 5% showed a negative TBW trend. From these pixels, 30–36% were diverging, indicating recovery, equivalent to 1–2% of all pixels in productive areas (Table 3). This was observed near the Loess plateau in northeast China (S1 Fig, both LI-NDVI and AS-NDVI) and the agricultural areas in northern India (S2 Fig, only AS-NDVI) and southwest Australia (S7 Fig).

The fraction of pixels with decoupled or diverging TBW and NDVI trends varied only slightly between metrics on a global scale, with exception of Max-NDVI (Table 3). However, differences between NDVI metrics were important when studying specific geographical regions, as illustrated in Fig 3 (Sahel) and Fig 4 (e.g. Belarus and Ukraine).

thumbnail
Fig 4.

a, Global map showing combinations of SLR-trends in TBW and PWR-trends in the two segments (S1 and S2 for the two segments) for the AS-NDVI metric. The “+” sign in the legend indicates a neutral or positive trend, a “–”sign a negative trend. For white areas, the temporal trend explained less than 10% of the observed variation in AS-NDVI or TBW. A continuing trend indicates that PWR did not improve upon a SLR trend and the later trend value was used; (b): as in (a) but now for the LI-NDVI metric. See supplementary materials for more detailed figures for separate regions.

https://doi.org/10.1371/journal.pone.0138013.g004

On a global scale, a decoupling between changes in vegetation activity and modelled productivity was visible for all datasets (Table 3 and compare Fig 5c and 5d). Generally, NDVI trends in the first PWR segment were positive and more neutral for the second PWR segment, whereas TBW trends remained positive in both PWR segments for most latitudes (Fig 5c and 5d). A negative TBW trend was found for latitudinal means between 25–40°S and 40–50°N. The year when a breakpoint between PWR segments was determined showed a latitudinal pattern, somewhat similar to mean NDVI values, with breakpoints identified in later years for latitudes with lower mean NDVI values (S8a and S8b Fig).

thumbnail
Fig 5.

a, Latitudinal means of NDVI metrics based on AS-NDVI and AO-NDVI for the years 1982–2010 (for legend see Fig 5c); b, mean year of breakpoint for pixels where PWR regression improved upon SLR; c, mean trends in three NDVI metrics for the first (S1) and second segment (S2) of PWR. Pixel trends were based on PWR where PWR improved the coefficient of determination with at least 0.1 and SLR trends otherwise; (d) as in (c) but now for TBW.

https://doi.org/10.1371/journal.pone.0138013.g005

Discussion

This work shows a clear difference in the area with monotonic and segmented trends and direction of these trends when comparing climate-constrained, modelled biomass production with observed NDVI for large parts of the globe. A decoupled NDVI-trend, with decreasing vegetation activity opposing a positive or neutral trend in climatic drivers, was observed for 17–36% of the pixels in productive areas depending on the NDVI metric used. Only a small proportion (6–9%) of these productive areas were trending downwards in NDVI over the full period from 1982–2010, and the majority of these areas showed a transition of trends, with first a positive and then a negative trend in NDVI. These findings are effected by uncertainties in the GIMMS3g dataset, as discussed elsewhere [62]. Most of the pixels with a decoupled trend had a neutral TBW trend, aligning with the overall dominance of neutral trends in TBW (72%). This could be due to the relatively short period (29 years) where impacts of changes in climatic conditions on TBW may be limited in comparison to year to year weather variation.

Observed changes in CO2 concentrations and temperature are major drivers of increasing landscape greenness over all latitudes, whereas the influence of changes in the amount of precipitation varies strongly on regional scales [4]. For the northern hemisphere above 50°N, strong greening kept pace with the increased temperatures when averaging over latitudes [3]. Greening rates in tropical mountains increase with altitude but greening reversed to browning in the mid-1990s at lower elevations [63]. We found that on a global scale, only 38–47% of pixels with a positive TBW were indeed greening up (as observed with NDVI), equivalent to 9–11% of all productive pixels. The large fraction of pixels with a positive (and often segmented) NDVI trend, in particular in agricultural areas, aligns with a positive trend in agricultural productivity and yield, most probably resulting from a combination of improved climatic conditions, increasing nutrient applications and soil fertility [64], rising CO2 concentrations and genetic improvements [65] and changes in seeding times [19]. Plateauing crop yields have been reported for various parts of the world [1518], for time intervals comparable to our study that can be linked to biophysical constraints to crop productivity [66]. Hence, the positive and often segmented NDVI trend in the 1981–2010 period is reasonably well understood for cropland areas considering the strong relationships between yield and NDVI [67,68].

Segmented trends caused by transitions from greening to browning in non-agricultural areas have been associated with changes in rainfall patterns and temperature, resulting in larger net evapotranspiration deficits [63,69] and decreased NPP [7]. Modelled biomass integrates these effects of temperature, rainfall amounts and distribution and shifts in season lengths [50,53]. For the majority of productive areas, the trend in modelled biomass was absent or very weak (R2 < 0.1), about 24% showed an increase in modelled production and only 5% a decrease. Mao et al. [4] found that modelled change in landscape greenness strongly varied with latitude, increasing in wet climates between 10°S and 15°N and above 50 °N but decreasing in drier climates between 10–40 °S and 10–30 °N. We also found large differences in TBW trends between latitudes, although TBW trends were positive over a wide range of latitudes, with a negative TBW trend in latitudes from 25–40°S and 40–50°N [4]. First signs of carbon sink saturation are reported for Europe [70,71]. Higher temperatures and longer growing seasons resulted in more biomass production, whereas climate change has not resulted in an increase of the land area under drought in the northern hemisphere [72]. This in contrast to the southern hemisphere [72], where for latitudes from 25–40 °S evapotranspiration deficits are increasing for large areas in all three continents [32]. In general, rising CO2 concentrations increase landscape greenness [4], rendering our estimates of changes in biomass conservative as CO2 concentration increases that accelerate growth were not included in our study. The modelled biomass production could capture the global spatial variability in MODIS-NPP reasonably well (R2 values of 0.63–0.75), with similar fractions of explained variation when compared to AS-NDVI (R2 values of 0.72–0.79). Part of the unexplained variation may arise from spatial variation in MODIS-NDVI as TBW was calculated with only climatic variables, in contrast to MODIS-NPP that also uses MODIS-NDVI.

Decoupling between positive climatic drivers and a negative NDVI response, predominantly in the last PWR segment, may be understood when considering constraints on landscape greenness other than climate. The vegetation present within the landscape may not be able to respond to changes in climate, e.g. a drought induced productivity declines of white spruce in Alaska [73] leading to a biome shift [69]. Changes in severity and area of burning and/or fire frequency [74] may be another important component affecting NDVI trends. Nutrient constraints are prevalent in wide areas on the globe [6]. Changes in nitrogen deposition may, therefore, increase or decrease NDVI [4]. Changes in rainfall amounts and distribution influence nutrient constraints in dryland areas, e.g. Delgado-Baquerizo et al. [75] found that aridity controls N and P availability. The intensity of rainfall affects the depth of the wetted root zone, determining the amount of nutrients in the soil profile that are available for plant uptake [76]. Nutrient constraints may also result in CO2 acclimation [77], restricting tree stem growth response to increased CO2 availability in temperate [78] and tropical forests [79], evidenced by a deteriorating mineral nutrition of trees in Europe [80].

Human influences are more difficult to assess. Land use change (e.g. due to increasing cropping areas) may result in a positive [81] or a negative trend in NDVI, the latter might occur when mean NDVI for natural vegetation is higher than for crops, such as in e.g. Argentina [82]. In this area, changes in the observed start of season indicate conversion of natural vegetation into crops [83]. However, on a global scale land use change is not expected to affect NDVI as strongly as the changes in climate or CO2 [4]. Results of Seaquist et al. [84] suggest that “demographic and agricultural pressures in the Sahel are unable to account for differences between simulated and observed vegetation dynamics, even for the most densely populated areas”. Firstly, the human influence may be masked by changes in natural vegetation. It is suggested that native vegetation in the Sahel has been in a long-term recovery period after the 1970–1980 drought [85], coinciding with increasing rain-use efficiency [21,86]. These increasing rainfall use efficiencies may be understood when considering limited access to nutrients without water sufficiently wetting the soil during drought periods [76], resulting in an increased nutrient availability in wetter periods. Secondly, soil degradation may lead to land abandonment and eventually lead to an increase in NDVI. All NDVI metrics showed greening in the peanut-basin of Senegal. This area has undergone land-abandonment, with agricultural land-use dropping from 80% to 67% between 1980 and 2000, resulting in increasing tree cover [87], and a higher NDVI as the natural vegetation has a stronger NDVI amplitude than crops [88]. Secondly, soil fertility decline on agricultural farms may be most visible on a scale much smaller than the course spatial resolution of the GIMMS3g dataset as according to Giller et al. [89] “differences in soil fertility between fields within a single farm may be as wide as those found between agroecological zones”.

Differences due to methodological choices strongly affected trend estimates and the fraction of pixels with a positive or negative trend in specific regions strongly differed between NDVI metrics. For tropical and sub-tropical broadleaf forest biomes, the explained fractions of variation in NPP or TBW were small for a few continents (e.g. North and South America). For these biomes the response of GIMMS3g-NDVI differed strongly from MODIS [90]. The NDVI is likely not the best index for these biomes with large amounts of green biomass due to saturation, and different results are obtained when other indices (e.g. the enhanced vegetation index, EVI) are used [13,14]. However, MODIS-EVI underestimates GPP in African landscapes [37], suggesting that selection of a vegetation index is not trivial. The differences between NDVI metrics indicate that regional assessments based on trends in a single NDVI metric, derived from e.g. the GIMMS3g dataset, needs to be interpreted with the highest caution as these trends are affected by a large degree of uncertainty. For global studies, it is critical to be able to select the appropriate NDVI metric for the purpose, but this is not a trivial choice and influences the trend estimate and direction of trends on regional scales. Further study is needed to better understand which NDVI-metrics are best suited to particular biomes or eco-regions. In the absence of metrics with long-term biomass observations suited for comparisons with GIMMS3g data, a comparison with independent data (e.g. MODIS or other sensors) is recommended [62].

Discontinuous trends in GIMMS3g time series result in differences between observed trends for particular epochs. NDVI trends were segmented for 32–48% of the pixels, confirming findings of others [28,91] and in line with the type of trends observed in climatic datasets [911]. The frequent presence of discontinuous trends can explain differences in trend estimates when a slightly different time-period is used [22]. This highlights the importance of these long-term datasets, enabling to differentiate between external factors influencing trends and long-term cyclic patterns in landscape greenness. The presence of discontinuous trends also indicates that there is uncertainty about the linearity and on the duration of a particular trend, requiring extra caution when extrapolating trends into the future. The approach presented here provides an improvement of our understanding of the vegetation responses to changes in climate and provides better means to monitor and assess the human influence on ecosystem productivity.

The observed segmented and often decoupled trends in ecosystem productivity opposing climatic drivers may have been caused by a combination of natural and human factors including land use change, declining soil fertility, nutrient limitations, biome shifts and effects of fire, but disentangling these remains challenging. This work stresses the importance of a better understanding of current limitations to plant growth and suggests that future scenarios must consider the influence of droughts and extreme weather events, nutrient availability and other constraints on predicted changes in NPP.

Supporting Information

S1 Fig. Cumulative frequency distribution of Piecewise linear regression (PWR) coefficient of determination (CD) values within biomes for the AS-NDVI (a), LI-NDVI (b) and AO-NDVI (c) metrics.

Legend numbers indicate cumulative frequency distributions of the following biomes: 1 Tropical and subtropical moist broadleaf forests; 2 Tropical and subtropical dry broadleaf forests; 3 Tropical and subtropical coniferous forests; 4 Temperate broadleaf and mixed forests; 5 Temperate conifer forests; 6 Boreal forests/taiga; 7 Tropical and subtropical grasslands, savannas and shrub lands; 8 Temp. grasslands, savannas and shrub lands; 9 Flooded grasslands and savannas; 10 Montane grasslands and shrub lands; 11 Tundra; 12 Mediterranean forests, woodlands and shrub; 13 Deserts and xeric shrub lands; 14 Mangroves.

https://doi.org/10.1371/journal.pone.0138013.s001

(TIF)

S2 Fig. South-East Asia region.

The map shows combinations of SLR-trends in TBW and PWR-trends in the two segments (S1 and S2 for the two segments) for (a) the AS-NDVI dataset and (b) the LI-NDVI dataset. The “+” sign in the legend indicates a neutral or positive trend, a “–”sign a negative trend. For white areas, the temporal trend explained less than 10% of the observed variation in AS-NDVI or TBW. A continuing trend indicates that PWR did not improve upon a SLR trend and the later trend value was used.

https://doi.org/10.1371/journal.pone.0138013.s002

(TIF)

S3 Fig. South and central America region.

The map shows combinations of SLR-trends in TBW and PWR-trends in the two segments (S1 and S2 for the two segments) for (a) the AS-NDVI dataset and (b) the LI-NDVI dataset. The “+” sign in the legend indicates a neutral or positive trend, a “–”sign a negative trend. For white areas, the temporal trend explained less than 10% of the observed variation in AS-NDVI or TBW. A continuing trend indicates that PWR did not improve upon a SLR trend and the later trend value was used.

https://doi.org/10.1371/journal.pone.0138013.s003

(TIF)

S4 Fig. Wider Europe region.

The map shows combinations of SLR-trends in TBW and PWR-trends in the two segments (S1 and S2 for the two segments) for (a) the AS-NDVI dataset and (b) the LI-NDVI dataset. The “+” sign in the legend indicates a neutral or positive trend, a “–”sign a negative trend. For white areas, the temporal trend explained less than 10% of the observed variation in AS-NDVI or TBW. A continuing trend indicates that PWR did not improve upon a SLR trend and the later trend value was used.

https://doi.org/10.1371/journal.pone.0138013.s004

(TIF)

S5 Fig. Africa.

The map shows combinations of SLR-trends in TBW and PWR-trends in the two segments (S1 and S2 for the two segments) for (a) the AS-NDVI dataset and (b) the LI-NDVI dataset. The “+” sign in the legend indicates a neutral or positive trend, a “–”sign a negative trend. For white areas, the temporal trend explained less than 10% of the observed variation in AS-NDVI or TBW. A continuing trend indicates that PWR did not improve upon a SLR trend and the later trend value was used.

https://doi.org/10.1371/journal.pone.0138013.s005

(TIF)

S6 Fig. North America.

The map shows combinations of SLR-trends in TBW and PWR-trends in the two segments (S1 and S2 for the two segments) for (a) the AS-NDVI dataset and (b) the LI-NDVI dataset. The “+” sign in the legend indicates a neutral or positive trend, a “–”sign a negative trend. For white areas, the temporal trend explained less than 10% of the observed variation in AS-NDVI or TBW. A continuing trend indicates that PWR did not improve upon a SLR trend and the later trend value was used.

https://doi.org/10.1371/journal.pone.0138013.s006

(TIF)

S7 Fig. Australia and surrounding region.

The map shows combinations of SLR-trends in TBW and PWR-trends in the two segments (S1 and S2 for the two segments) for (a) the AS-NDVI dataset and (b) the LI-NDVI dataset. The “+” sign in the legend indicates a neutral or positive trend, a “–”sign a negative trend. For white areas, the temporal trend explained less than 10% of the observed variation in AS-NDVI or TBW. A continuing trend indicates that PWR did not improve upon a SLR trend and the later trend value was used.

https://doi.org/10.1371/journal.pone.0138013.s007

(TIF)

S8 Fig.

(a) Latitudinal means of NDVI metrics based on AS-NDVI and AO-NDVI for the years 1982–2010; (b) Mean year of breakpoint for pixels where PWR regression improved upon SLR; mean trends in (c) three NDVI metrics and (d) TBW for the first (S1) and second segment (S2) of PWR. Pixel trends were based on PWR where PWR improved the coefficient of determination with at least 0.1 and SLR trends otherwise.

https://doi.org/10.1371/journal.pone.0138013.s008

(TIF)

S1 File. Supplementary material including tables.

https://doi.org/10.1371/journal.pone.0138013.s009

(PDF)

Acknowledgments

The authors thank the NASA Global Inventory Modelling and Mapping Studies (GIMMS) group for producing and sharing the AVHRR GIMMS3g NDVI dataset. The MODIS GPP/NPP project is kindly acknowledged for producing and sharing the MOD17A3 dataset. Authors gratefully acknowledge the contributions of M. Cherlet, M. Bakkenes, Z. Bai and anonymous reviewers to this work.

Author Contributions

Conceived and designed the experiments: AS EI JC BB. Analyzed the data: AS EI JC. Contributed reagents/materials/analysis tools: RF. Wrote the paper: AS EI JC BB RF.

References

  1. 1. Nemani RR, Keeling CD, Hashimoto H, Jolly WM, Piper SC, Tucker CJ, et al. (2003) Climate-driven increases in global terrestrial net primary production from 1982 to 1999. Science 300: 1560–1563. pmid:12791990
  2. 2. Dent DL, Bai Z, Schaepman M, Olsson L (2009) Response to Wessels: Comments on 'Proxy global assessment of land degradation'. Soil Use Manag. 25: 93–97.
  3. 3. Xu L, Myneni RB, Chapin FS, Callaghan TV, Pinzon JE, Tucker CJ, et al. (2013) Temperature and vegetation seasonality diminishment over northern lands. Nat. Clim. Change 3: 581–586.
  4. 4. Mao J, Shi X, Thornton PE, Hoffman FM, Zhu Z, Myneni RB. (2013) Global latitudinal-asymmetric vegetation growth trends and their driving mechanisms: 1982–2009. Remote Sens. 5: 1484–1497.
  5. 5. Los SO (2013) Analysis of trends in fused AVHRR and MODIS NDVI data for 1982–2006: Indication for a CO2 fertilization effect in global vegetation. Global Biogeochem. Cycles 27: 318–330.
  6. 6. Fisher JB, Badgley G, Blyth E (2012) Global nutrient limitation in terrestrial vegetation. Global Biogeochem. Cycles 26: GB3007.
  7. 7. Zhao M, Running SW (2010) Drought-induced reduction in global terrestrial net primary production from 2000 through 2009. Science 329: 940–943. pmid:20724633
  8. 8. Vicente-Serrano SM, Gouveia C, Camarero JJ, Beguería S, Trigo R, López-Moreno JI et al. (2013) Response of vegetation to drought time-scales across global land biomes. Proc. Natl. Acad. Sci. USA 110: 52–57. pmid:23248309
  9. 9. Tomé AR, Miranda PMA (2004) Piecewise linear fitting and trend changing points of climate parameters. Geophys. Res. Lett. 31: L02207 02201–02204.
  10. 10. Visser H (2004) Estimation and detection of flexible trends. Atmos. Environ. 38: 4135–4145.
  11. 11. Bates BC, Chandler RE, Bowman AW (2012) Trend estimation and change point detection in individual climatic series using flexible regression methods. J Geophys Res 117: D16106.
  12. 12. Ivits E, Cherlet M, Sommer S, Mehl W (2013) Addressing the complexity in non-linear evolution of vegetation phenological change with time-series of remote sensing images. Ecol. Indicators 26: 49–60.
  13. 13. Huete A, Saleska S. Remote sensing of tropical forest phenology: issues and controversies; 2010; Alice Springs, 13–17 September 2010. Available: http://www.scribd.com/doc/37456541/15arspc-Submission-136.
  14. 14. Huete AR, Didan K, Shimabukuro YE, Ratana P, Saleska SR, Hutyra LR, et al. (2006) Amazon rainforests green-up with sunlight in dry season. Geophys. Res. Lett. 33: L06405.
  15. 15. Brisson N, Gate P, Gouache D, Charmet G, Oury F-X, Huard F (2010) Why are wheat yields stagnating in Europe? A comprehensive data analysis for France. Field Crops Res. 119: 201–212.
  16. 16. Lobell DB, Ortiz-Monasterio JI, Asner GP, Matson PA, Naylor RL, Falcon WP (2005) Analysis of wheat yield and climatic trends in Mexico. Field Crops Res. 94: 250–256.
  17. 17. Lobell DB, Schlenker W, Costa-Roberts J (2011) Climate Trends and Global Crop Production Since 1980. Science 333: 616–620. pmid:21551030
  18. 18. Maltais-Landry G, Lobell DB (2012) Evaluating the Contribution of Weather to Maize and Wheat Yield Trends in 12 US Counties. Agron. J. 104: 301–311.
  19. 19. Lobell DB, Ortiz-Monasterio JI, Sibley AM, Sohu VS (2013) Satellite detection of earlier wheat sowing in India and implications for yield trends. Agric. Sys. 115: 137–143.
  20. 20. De Jong R, Verbesselt J, Zeileis A, Schaepman ME (2013) Shifts in global vegetation activity trends. Remote Sens. 5: 1117–1133.
  21. 21. Fensholt R, Rasmussen K, Kaspersen P, Huber S, Horion S, Swinnen E (2013) Assessing land degradation/recovery in the african sahel from long-term earth observation based primary productivity and precipitation relationships. Remote Sens. 5: 664–686.
  22. 22. Dardel C, Kergoat L, Hiernaux P, Mougin E, Grippa M, Tucker CJ (2014) Re-greening Sahel: 30 years of remote sensing data and field observations (Mali, Niger). Remote Sens. Environ. 140: 350–364.
  23. 23. Bai ZG, Dent DL, Olsson L, Schaepman ME (2008) Proxy global assessment of land degradation. Soil Use Manag. 24: 223–234.
  24. 24. Cook BI, Pau S (2013) A global assessment of long-term greening and browning trends in pasture lands using the GIMMS LAI3g dataset. Remote Sens. 5: 2492–2512.
  25. 25. Wessels KJ, Prince SD, Malherbe J, Small J, Frost PE, VanZyl D (2007) Can human-induced land degradation be distinguished from the effects of rainfall variability? A case study in South Africa. J. Arid Environ. 68: 271–297.
  26. 26. Ivits E, Cherlet M, Mehl W, Sommer S (2013) Ecosystem functional units characterized by satellite observed phenology and productivity gradients: A case study for Europe. Ecol. Indicators 27: 17–28.
  27. 27. Ivits E, Cherlet M, Tóth G, Sommer S, Mehl W, Vogt J, et al. (2012) Combining satellite derived phenology with climate data for climate change impact assessment. Global Planet. Change 88–89: 85–97.
  28. 28. Jamali S, Seaquist J, Eklundh L, Ardö J (2014) Automated mapping of vegetation trends with polynomials using NDVI imagery over the Sahel. Remote Sens. Environ. 141: 79–89.
  29. 29. Verbesselt J, Hyndman R, Newnham G, Culvenor D (2010) Detecting trend and seasonal changes in satellite image time series. Remote Sens. Environ. 114: 106–115.
  30. 30. de Jong R, de Bruin S, Schaepman M, Dent D (2011) Quantitative mapping of global land degradation using earth observations. Int. J. Remote Sens. 32: 6823–6853.
  31. 31. Wessels KJ (2009) Comments on 'Proxy global assessment of land degradation' by Bai et al. (2008). Soil Use Manag. 25: 91–92.
  32. 32. Ivits E, Horion S, Fensholt R, Cherlet M (2014) Global Ecosystem Response Types Derived from the Standardized Precipitation Evapotranspiration Index and FPAR3g series. Remote Sens. 6: 4266–4288.
  33. 33. Mbow C, Fensholt R, Rasmussen K, Diop D (2013) Can vegetation productivity be derived from greenness in a semi-arid environment? Evidence from ground-based measurements. J. Arid Environ. 97: 56–65.
  34. 34. IPCC-WGIAR5 (2013) Working group I contribution to the IPCC 5th Assessment Report Climate Change 2013: the Physical Basis. Online: IPCC. Report. www.ipcc.ch/report/ar5/wg1/
  35. 35. Running SW, Nemani RR, Heinsch FA, Zhao M, Reeves M, Hashimoto H (2004) A continuous satellite-derived measure of global terrestrial primary production. Bioscience 54: 547–560.
  36. 36. Turner DP, Ritts WD, Cohen WB, Gower ST, Running SW, Zhao M, et al. (2006) Evaluation of MODIS NPP and GPP products across multiple biomes. Remote Sens. Environ. 102: 282–292.
  37. 37. Sjöström M, Ardö J, Arneth A, Boulain N, Cappelaere B, Eklundh L, et al. (2011) Exploring the potential of MODIS EVI for modeling gross primary production across African ecosystems. Remote Sens. Environ. 115: 1081–1089.
  38. 38. de Jong R, Schaepman ME, Furrer R, de Bruin S, Verburg PH (2013) Spatial relationship between climatologies and changes in global vegetation activity. Global Change Biol. 19: 1953–1964.
  39. 39. Hein L, De Ridder N (2006) Desertification in the Sahel: a reinterpretation. Global Change Biol. 12: 751–758.
  40. 40. Stephens DJ, Lyons TJ (1998) Rainfall-yield relationships across the Australian wheatbelt. Aust. J. Agric. Res. 49: 211–223.
  41. 41. Stephens DJ, Lyons TJ, Lamond MH (1989) A simple model to forecast wheat yield in Western Australia. J. R. Soc. West. Aust. 71: 77–81.
  42. 42. Tucker CJ, Pinzon JE, Brown ME, Slayback DA, Pak EW, Mahoney R, et al. (2005) An extended AVHRR 8-km NDVI dataset compatible with MODIS and SPOT vegetation NDVI data. Int. J. Remote Sens. 26: 4485–4498.
  43. 43. Vermote E, Kaufman YJ (1995) Absolute calibration of AVHRR visible and near-infrared channels using ocean and cloud views. Int. J. Remote Sens. 16: 2317–2340.
  44. 44. Los SO (1998) Estimation of the ratio of sensor degradation between NOAA AVHRR channels 1 and 2 from monthly NDVI composites. IEEE Trans. Geosci. Remote Sens. 36: 206–213.
  45. 45. Pinzon JE, Brown ME, Tucker CJ (2005) Hilbert-Huang transform and its applications. In: Huang NE, editor. Interdisciplinary mathematical sciences v. 5. Singapore; Hackensack NJ; London: World Scientific Publishing Co. pp. 311.
  46. 46. Jonsson P, Eklundh L (2004) TIMESAT—a program for analyzing time-series of satellite sensor data. Comp. Geosci. 30: 833–845.
  47. 47. Zhao M, Heinsch FA, Nemani RR, Running SW (2005) Improvements of the MODIS terrestrial gross and net primary production global data set. Remote Sens. Environ. 95: 164–176.
  48. 48. Kier G, Mutke J, Dinerstein E, Ricketts TH, Küper W, Kreft H, et al. (2005) Global patterns of plant diversity and floristic knowledge. J. Biogeogr. 32: 1107–1116.
  49. 49. Bartholomé E, Belward AS (2005) GLC2000: A new approach to global land cover mapping from earth observation data. Int. J. Remote Sens. 26: 1959–1977.
  50. 50. Conijn JG, Querner E, Rau ML, Hengsdijk H, Kuhlman T, Meijerink G, et al. (2011) Agricultural Resource Scarcity and Distribution: a Case Study of Crop Production in Africa. Wageningen: Plant Research International. Report 380. Available: http://edepot.wur.nl/165556
  51. 51. Bouman BAM, Van Keulen H, Van Laar HH, Rabbinge R (1996) The 'School of de Wit' crop growth simulation models: A pedigree and historical overview. Agric. Sys. 52: 171–198.
  52. 52. Van Ittersum MK, Leffelaar PA, Van Keulen H, Kropff MJ, Bastiaans L,Goudriaan J (2003) On approaches and applications of the Wageningen crop models. Eur. J. Agron. 18: 201–234.
  53. 53. Jing Q, Conijn SJG, Jongschaap REE, Bindraban PS (2012) Modeling the productivity of energy crops in different agro-ecological environments. Biomass Bioenergy 46: 618–633.
  54. 54. Bassu S, Brisson N, Durand J-L, Boote K, Lizaso J, Jones JW, et al. (2014) How do various maize crop models vary in their responses to climate change factors? Global Change Biol. 20: 2301–2320.
  55. 55. CRU (2008) Climate Research Unit (CRU) University of East Anglia. CRU Time Series (CRUTS) high resolution gridded datasets. Available: http://catalogue.ceda.ac.uk/search?search_term=CRU&return_obj=ob&search_obj=ob: NCAS British Atmospheric Data Centre.
  56. 56. Harris I, Jones PD, Osborn TJ, Lister DH (2014) Updated high-resolution grids of monthly climatic observations—the CRU TS3.10 Dataset. International Journal of Climatology 34: 623–642.
  57. 57. Conijn JG, Bai ZG, Bindraban PS, Rutgers B (2013) Global changes in net primary productivity, affected by climate and abrubt land use changes since 1981. Towards mapping global soil degradation. Wageningen, the Netherlands: ISRIC. Report 2013/01. Available: http://www.isric.org/sites/default/files/Report%202013_01_total_SEC.pdf
  58. 58. Batjes NH (1996) Development of a world data set of soil water retention properties using pedotransfer rules. Geoderma 71: 31–52.
  59. 59. FAO (1996) Digital soil map of the world and derived soil properties, version 3.5, November 1995, derived from the FAO/UNESCO soil map of the world, original scale 1:5 000 000. Rome: Food and Agricultural Organisation.
  60. 60. Erb K-H, Gaube V, Krausmann F, Plutzar C, Bondeau A, Haberl H, et al. (2007) A comprehensive global 5 min resolution land-use data set for the year 2000 consistent with national census data. J. Land Use Sci. 2: 191–224.
  61. 61. Fernandes R, G Leblanc S (2005) Parametric (modified least squares) and non-parametric (Theil–Sen) linear regressions for predicting biophysical parameters in the presence of measurement errors. Remote Sens. Environ. 95: 303–316.
  62. 62. Guay KC, Beck PSA, Berner LT, Goetz SJ, Baccini A, Buermann W (2014) Vegetation productivity patterns at high northern latitudes: A multi-sensor satellite data assessment. Global Change Biol. 20: 3147–3158.
  63. 63. Krishnaswamy J, John R, Joseph S (2014) Consistent response of vegetation dynamics to recent climate change in tropical mountain regions. Global Change Biol. 20: 203–215.
  64. 64. Huang Y, Sun W (2006) Changes in topsoil organic carbon of croplands in mainland China over the last two decades. Chin. Sci. Bull. 51: 1785–1803.
  65. 65. Rijk B, van Ittersum M, Withagen J (2013) Genetic progress in Dutch crop yields. Field Crops Res. 149: 262–268.
  66. 66. Grassini P, Eskridge KM, Cassman KG (2013) Distinguishing between yield advances and yield plateaus in historical crop production trends. Nature Comm. 4: 1–11.
  67. 67. Schut AGT, Stephens DJ, Stovold RGH, Adams M, Craig RL (2009) Improved wheat yield and production forecasting with a moisture stress index, AVHRR and MODIS data. Crop Pasture Sci. 60: 60–70.
  68. 68. Mkhabela MS, Bullock P, Raj S, Wang S, Yang Y (2011) Crop yield forecasting on the Canadian Prairies using MODIS NDVI data. Agric. For. Meteorol. 151: 385–393.
  69. 69. Beck PSA, Juday GP, Alix C, Barber VA, Winslow SE, Sousa EE, et al. (2011) Changes in forest productivity across Alaska consistent with biome shift. Ecol. Lett. 14: 373–379. pmid:21332901
  70. 70. Ruiz-Benito P, Madrigal-González J, Ratcliffe S, Coomes DA, Kändler G, Lehtonen A, et al. (2014) Stand Structure and Recent Climate Change Constrain Stand Basal Area Change in European Forests: A Comparison Across Boreal, Temperate, and Mediterranean Biomes. Ecosystems 17: 1439–1454.
  71. 71. Nabuurs GJ, Lindner M, Verkerk PJ, Gunia K, Deda P, Michalak R, et al. (2013) First signs of carbon sink saturation in European forest biomass. Nat. Clim. Change 3: 792–796.
  72. 72. Damberg L, AghaKouchak A (2014) Global trends and patterns of drought from space. Theor. Appl. Clim. 117: 441–448.
  73. 73. Wilmking M, Juday GP, Barber VA, Zald HSJ (2004) Recent climate warming forces contrasting growth responses of white spruce at treeline in Alaska through temperature thresholds. Global Change Biol. 10: 1724–1736.
  74. 74. Giglio L, Randerson JT, Van Der Werf GR (2013) Analysis of daily, monthly, and annual burned area using the fourth-generation global fire emissions database (GFED4). Journal of Geophysical Research: Biogeosciences 118: 317–328.
  75. 75. Delgado-Baquerizo M, Maestre FT, Gallardo A, Bowker MA, Wallenstein MD, Quero JL, et al. (2013) Decoupling of soil nutrient cycles as a function of aridity in global drylands. Nature 502: 672–676. pmid:24172979
  76. 76. Breman H, Groot JJR, Van Keulen H (2001) Resource limitations in Sahelian agriculture. Global Environ. Change 11: 59–68.
  77. 77. Gunderson CA, Wullschleger SD (1994) Photosynthetic acclimation in trees to rising atmospheric CO2: A broader perspective. Photosynthesis Res. 39: 369–388.
  78. 78. Bader MKF, Leuzinger S, Keel SG, Siegwolf RTW, Hagedorn F, Schleppi P, et al. (2013) Central european hardwood trees in a high-CO2 future: Synthesis of an 8-year forest canopy CO2 enrichment project. J. Ecol. 101: 1509–1519.
  79. 79. van der Sleen JP (2014) Environmental and physiological drivers of tree growth: a pan-tropical study of stable isotopes in tree rings. Wageningen: PhD Thesis Wageningen University.
  80. 80. Jonard M, Fürst A, Verstraeten A, Thimonier A, Timmermann V, Potocic N, et al. (2015) Tree mineral nutrition is deteriorating in Europe. Global Change Biol. 21: 418–430.
  81. 81. Pouliot D, Latifovic R, Olthof I (2009) Trends in vegetation NDVI from 1 km AVHRR data over Canada for the period 1985–2006. Int. J. Remote Sens. 30: 149–168.
  82. 82. Volante JN, Alcaraz-Segura D, Mosciaro MJ, Viglizzo EF, Paruelo JM (2012) Ecosystem functional changes associated with land clearing in NW Argentina. Agric., Ecosyst. Environ. 154: 12–22.
  83. 83. van Leeuwen W, Hartfield K, Miranda M, Meza F (2013) Trends and ENSO/AAO Driven Variability in NDVI Derived Productivity and Phenology alongside the Andes Mountains. Remote Sens. 5: 1177–1203.
  84. 84. Seaquist JW, Hickler T, Eklundh L, Ardö J, Heumann BW (2009) Disentangling the effects of climate and people on Sahel vegetation dynamics. Biogeosciences 6: 469–477.
  85. 85. Bégué A, Vintrou E, Ruelland D, Claden M, Dessay N (2011) Can a 25-year trend in Soudano-Sahelian vegetation dynamics be interpreted in terms of land use change? A remote sensing approach. Global Environ. Change 21: 413–420.
  86. 86. Fensholt R, Rasmussen K (2011) Analysis of trends in the Sahelian 'rain-use efficiency' using GIMMS NDVI, RFE and GPCP rainfall data. Remote Sens. Environ. 115: 438–451.
  87. 87. Tappan GG, Sall M, Wood EC, Cushing M (2004) Ecoregions and land cover trends in Senegal. J. Arid Environ. 59: 427–462.
  88. 88. Bobée C, Ottlé C, Maignan F, De Noblet-Ducoudré N, Maugis P, Lézine AM, et al. (2012) Analysis of vegetation seasonality in Sahelian environments using MODIS LAI, in association with land cover and rainfall. J. Arid Environ. 84: 38–50.
  89. 89. Giller KE, Tittonell P, Rufino MC, van Wijk MT, Zingore S, Mapfumo P, et al. (2011) Communicating complexity: Integrated assessment of trade-offs concerning soil fertility management within African farming systems to support innovation and development. Agric. Sys. 104: 191–203.
  90. 90. Fensholt R, Proud SR (2012) Evaluation of Earth Observation based global long term vegetation trends—Comparing GIMMS and MODIS global NDVI time series. Remote Sens. Environ. 119: 131–147.
  91. 91. de Jong R, Verbesselt J, Schaepman ME, de Bruin S (2012) Trend changes in global greening and browning: Contribution of short-term trends to longer-term change. Global Change Biol. 18: 642–655.