Increasing interannual variability of global vegetation greenness

Despite the long-term greening trend in global vegetation identified in previous investigations, changes in the interannual variability (IAV) of vegetation greenness over time is still poorly understood. Using Global Inventory Modeling and Mapping Studies normalized difference vegetation index (NDVI) third generation data and corresponding meteorological data from 1982 to 2015, we studied the changes and drivers of the IAV of vegetation greenness as indicated by the coefficient of variation of vegetation greenness at a global scale. Dry and high-latitude areas exhibited high NDVI variability whereas humid areas exhibited relatively low NDVI variability. We detected an increase in the global IAV of vegetation greenness over time using a 15 year moving window. Spatially, we observed significant increases in the IAV of vegetation greenness in greater than 45% of vegetated areas globally and decreases in 21%. Our comparison of ecological models suggests good performance in terms of simulating spatial differences in vegetation variability, but relatively poor performance in terms of capturing changes in the IAV of vegetation greenness. Furthermore, the dominant climate variables controlling changes in the IAV of vegetation greenness were determined spatially using principal component regression and partial least squares regression. The two methods yielded similar patterns, revealing that temperature exerted the biggest influence on changes in the IAV of vegetation greenness, followed by solar radiation and precipitation. This study provides insights into global vegetation variability which should contribute to an understanding of vegetation dynamics in the context of climate change.


Introduction
Terrestrial ecosystems are important global carbon sinks that can offset a substantial proportion of anthropogenic CO 2 emissions (Beer et al 2010, Pan et al 2011, Ahlström et al 2015a. Terrestrial vegetation strongly regulates atmospheric CO 2 concentrations in terms of accumulation, seasonality, and interannual variability (IAV) (Keenan et al 2016, Yuan et al 2018). Hence, monitoring the drivers of ecosystem dynamics is crucial for forecasting future climate change.
Satellite observations reveal a greening trend on Earth's surface since the 1980s .
However, this does not necessarily indicate that greening improves ecosystem functioning because longterm averages ignore IAV which is a critical indicator of ecosystem vulnerability or resilience to external disturbances (Thornton et al 2014, Smith et al 2014, Ray et al 2015, Sloat et al 2018. Low IAV in vegetation greenness indicates that vegetation growth is stable and not easily disturbed by external factors (Thornton et al 2014), whereas high IAV indicates unstable vegetation growth that has undergone large disturbances and remains vulnerable to future disturbance. For example, low IAV of crop yields indicates that food supplies and farm incomes are stable (Reidsma et al 2010, Ray et al 2013. In ecosystems, indicators of vegetation greenness and ecosystem production have been used to study IAV (Ahlström et al 2015a, Zhang et al 2016, Jiang et al 2017, Hu et al 2018, but it is unclear how IAV changes over the long term. A global survey of how the IAV of ecosystem productivity changes over time is critical for deepening our understanding of ecosystem resilience to climate change. Climate change is expected to trigger more frequent climate extremes, including droughts and heatwaves (IPCC 2013). Such changes in climate variability can have profound impacts on ecosystem productivity (Smith et al 2014). Ray et al (2015) analyzed the effects of climate variability on crop yields and found that one-third of global crop yield variability could be explained by climate variability. Le Hourou et al (1988) found that precipitation variability influences variability in ecosystem productivity. Sloat et al (2018) found an increasing influence of precipitation variability on productivity in grazing lands globally. However, it is presently unclear how the IAV of ecosystem productivity responds to climate change.
Here, we used the satellite-based normalized difference vegetation index (NDVI) (Rouse et al 1973) as a proxy for ecosystem productivity to examine temporal-spatial patterns in IAV from 1982 to 2015 and explored the potential climatic drivers of the IAV of vegetation greenness.

Data sets
The NDVI is typically used to detect vegetation greenness. Here we used the Global Inventory Modeling and Mapping Studies (GIMMS) NDVI third generation data for the time series 1982-2015, which has a spatial resolution of 1/12 degree and a temporal resolution of 15 d. The maximum value composite (MVC) method (Holben 1986) was adopted to aggregate the biweekly data into monthly intervals. This method can largely remove the contaminations of cloud and atmospheric noises (Holben 1986). Barren regions with an annual mean NDVI less than 0.1 were excluded from the analysis (Piao et al 2005. Climate data from the Climatic Research Unit (CRU) Time Series 4.01 with a 0.5°resolution of monthly temperature (TMP) and precipitation (PRE) from 1901 to 2016 were used in this study (Harris et al 2014). Monthly shortwave radiation (SWD) data with a resolution of 0.5°was obtained from the CRU-National Centers for Environmental Prediction (CRU-NCEP) v7 data set.
To evaluate the performance of different Dynamic Global Vegetation Models (DGVMs) in capturing the IAV of vegetation growth, we compared gross primary productivity (GPP) simulations from eight DGVMs archived in the Trends in Net Land-Atmosphere Exchange (TRENDY-v6) project from 1982 to 2015 with observations of NDVI variability (Sitch et al 2015). We selected the following eight ecological models: CABLE, DLEM, ISAM, LPJ-wsl, OCN, ORCHI-DEE, VEGAS, and VISIT.
In addition, global land cover data with a spatial resolution of 0.05°was obtained from MODIS Land Cover Climate Modeling Grid (MCD12C1) Version 6 from 2001 to 2015. We used the International Geosphere and Biosphere Programme land cover classifications to select nine biomes (figure S1(a) is available online at stacks.iop.org/ERL/14/124005/mmedia and table S1), including evergreen needleleaf forests (ENF), deciduous needleleaf forests (DNF), evergreen broadleaf forests (EBF), deciduous broadleaf forests (DBF), mixed forests (MF), savannas, shrublands, grasslands, and croplands. We used only areas (pixels) with a constant vegetation type from 2001 to 2015 to prevent land cover changes from impacting our results.
Global climate zones were determined based on Köppen-Geiger classifications (Kottek et al 2006) obtained from http://koeppen-geiger.vu-wien.ac.at/. We focused on five major climate zones: equatorial, arid, warm temperate, snow, and polar, as shown in figure S1(b) and

Coefficient of variation
Here, we used the coefficient of variation (CV) of NDVI (NDVI_ CV ) and the CV of GPP (GPP_ CV ) to represent in the IAV of vegetation growth (Fang et al 2001, Schucknecht et al 2013, Sloat et al 2018. CV is an absolute value that measures the degree of dispersion of a variable. Its size is affected not only by the degree of variable dispersion, but also by the average of the variable's value. The formula is generally expressed as: where, SD represents the standard deviation and MN represents the mean of the studied variables. Both the annual growing season and annual maximum NDVI_cv were determined here.
To test the sensitivity of CV of studied variable to SD and MN, a simple linear regression (Yang et al 2016) of CV against SD and MN was performed. The sensitivity was determined according the slope of CV to SD and MN, and a larger slope indicates a higher sensitivity. Here, the growing season for land plants was determined from the GIMMS NDVI3g data set using a Savitzky-Golay filter that effectively smoothed the original curve and somewhat reduced the noise in the time series data (Chen et al 2004, Jönsson and Eklundh 2004, Cong et al 2012. Figure S2 shows the global patterns in surface phenology from 1982 to 2015, average start day and length of the growing season.

Mann-Kendall (MK) statistical test
To study the dynamics of NDVI_ CV and GPP_ CV , we applied a 15 yr moving time window to the NDVI/ GPP sequence from 1982 to 2015. A CV value representing the IAV of vegetation growth over a specific 15 yr period was calculated and assigned to the middle year of each window. This approach generated a series of 20 CV values with the first value representing the first 15 yr period from 1982 to 1996 and the 20th CV value representing the last 15 yr period from 2001 to 2015. The trend in window position values was detected using the Mann-Kendall (MK) non-parametric test method (Mann 1945). The advantage of the method is that the series does not need to follow a specific distribution and is seldom disturbed by abnormal values, so the calculation is straightforward and convenient (Hamed 2008). We applied the trendfree pre-whitening technique to the studied time series (Yue et al 2002, Erasmi et al 2014 before performing the trend test because the autocorrelation in the time series has been found to influence the significance of the Mann-Kendall trend test (Yue et al 2002). The strength of the time series trend was estimated using the Theil-Sen method (Fensholt et al 2012, Schucknecht et al 2013). This trend estimator calculates the slope of every single data pair in a time series and uses the median slope to characterize a trend in the data (Sen 1968). The significance of the CV trend was assessed at a p-value less than 0.05. In addition, we also examined NDVI_ CV trends determined from 11, 13, 17, and 19 yr windows to test the potential influence of window length on our results.

Principal component regression (PCR) and partial least squares regression (PLSR)
Several factors were selected to explore the drivers of NDVI_ CV , including the mean values of climate variables (TMP_ mean , PRE_ mean , and SWD_ mean ) and the CVs of climate variables (TMP_ CV , PRE_ CV , and SWD_ CV ) corresponding to each 15 yr window. Regression analysis between NDVI_ CV and the six potential drivers were performed to determine the main drivers. To reduce the uncertainty caused by the analysis method, two regression methods, the principal component regression (PCR) and the partial least squares regression (PLSR), were used in this study. Compared with traditional multivariate linear regression, these two methods can effectively overcome the multicollinearity problem existing among two or more explanatory variables. PCR is based on principal component analysis (Jolliffe 1982), which decomposes explanatory variables into several principal components. Instead of regressing the original dependent variable on the explanatory variables directly, the principal components of the explanatory variables are used for regression analysis (Jolliffe 1982). The PLSR is similar to PCR in that it integrates the basic functions of multiple linear regression analysis, canonical correlation analysis, and principal component analysis (Chun andKeleş 2010, Mehmood et al 2012). Variable importance in the projection (VIP) is the PLSR standard for judging the importance of independent variables in modeling. When VIP is greater than or equal to 1, its corresponding independent variables have strong explanatory significance to the dependent variables (Chong and Jun 2005). All variables were standardized before calculation.
We noted that data processing was performed in MATLAB 2014a environment, and global maps were drawn using ArcGIS 10.3. Figure 1(a) shows the spatial distribution of the IAV of NDVI (NDVI_ CV ) from 1982 to 2015. NDVI_ CV demonstrates high spatial heterogeneity. We observed relatively low IAV of vegetation greenness in the equatorial and snow climate zones (figures 1(b), S1), including tropical rainforests near the equator and boreal forests between 50°and 60°N latitude. In contrast, we observed high CV values in the vegetated lands of arid and polar climates zones, including western and southwestern Asia, Australia, northern and southern Africa, southwestern North America, and areas above 60°N. At the biomes scale (figure 1(c)), evergreen broadleaf forests exhibited the lowest IAV of vegetation greenness, whereas grasslands demonstrated the highest IAV.

Dynamics of variability of vegetation greenness
The trend of NDVI_ CV sequences generated from the 15 yr moving window from 1982 to 2015 was detected for each grid cell, as shown in figure 2. Significant vegetation IAV, indicated by annual maximum NDVI (figure 2(a)) decreases, were observed in more than 21% of vegetated lands scattered across Europe, the Indian peninsula, the sahel region, the southwest of America and the north of Asia. This indicates increasing stability of local ecosystems. Simultaneously, significant increases in vegetation variability occurred in northeastern China, western Asia, the southwestern United States, much of Australia, and central Africa, accounting for more than 45% of vegetated areas globally. The trend in growing season NDVI_ CV (figure 2(b)) roughly agreed with that of annual maximum NDVI_ CV . We observed decreases in growing season NDVI_ CV in 30% of vegetated areas and increases in 40% of vegetated areas globally. Figure 2(c) shows the probability density function of the trends in annual maximum NDVI_ CV and growing season NDVI_ CV . We observed a higher proportion of increases than decreases. In addition, growing season NDVI_ CV demonstrated more increasing trends than did annual maximum NDVI_ CV . Note that the NDVI_ CV trends determined using different moving time windows present a similar pattern with that of the 15 yr window (figure S3), suggesting a basic reliability of NDVI_ CV change built here.
NDVI_ CV changes can be caused by variations in standard deviation or the mean value of NDVI, so we examined trends in these two values and tested the sensitivity of NDVI_ CV to them (figure S4). The mean values of annual maximum NDVI displayed increasing trends in most regions (figure S4(a)), with decreasing trends in northeast China, the northern of North America, and the Congo rainforest. Trends in NDVI standard deviation were similar to those of NDVI_ CV ( figure S4(b)). Sensitivity tests suggest that NDVI_ CV was more sensitive to the standard deviation of NDVI than to the mean value of NDVI in nearly all continents (figure S4(c)), suggesting that standard deviation exerted a dominant influence over NDVI_ CV trends.
We further explored trends in regional mean annual maximum NDVI_ CV by dividing global vegetated lands into different regions (figure 3). Global mean NDVI_ CV demonstrated a significant increasing trend and the slope of NDVI_ CV in the Southern Hemisphere (SH) was larger than that in the Northern Hemisphere (NH). For different climate zones, the largest increase in NDVI_ CV was observed in the vegetated lands of polar and arid climate zones. In contrast, vegetation in humid and warm temperate climate zones presented the lowest increase in variability. At the biome scale, significant increases in NDVI_ CV were found for the majority of vegetation types, including shrublands, grasslands, evergreen broadleaf forests, deciduous needleleaf forests, and savannas. Simultaneously, NDVI_ CV exhibited decreasing trends in mixed forests and croplands, but the trends were not statistically significant.

Performance of ecosystem models in simulating vegetation variability
To examine whether the observed IAV of vegetation growth can be captured by ecological models, we studied the GPP_ CV simulated from TRENDY-models  1(a)), producing high CV values in drier regions and low values in humid and high-latitude regions. We further examined the patterns of GPP_ CV determined using a 15 yr moving widow, as shown in figure 4(b). The models only partially reproduced the spatial patterns of NDVI_ CV trend (figures 2(a) and (b)), including the increasing CV in northeastern China, western Asia, southern North America, and the Congo. However, opposite trends were observed in northern Asia, eastern Europe, and northern Australia, suggesting that the performance of current ecosystem models in modelling the IAV of vegetation growth still needs to be improved.

Drivers of variability of vegetation greenness
To explore the drivers of NDVI_ CV , we first examined the partial correlations between NDVI_ CV and climate variability and mean state over different periods ( figure S5). There was a strong positive correlation between TMP_ CV and NDVI_ CV in Europe and Australia, and a large continuous negative correlation in southern North America, central South America, and southwestern and northeastern Asia. TMP_ mean presented more extensive positive impacts on NDVI_ CV (22% of global vegetated lands) than did TMP_ CV (16% of global vegetated lands), mainly concentrating in relatively humid regions such as the eastern South America, southern China, and the Congo rainforest. Negative impacts between TMP_ mean and NDVI_ CV were also observed over 17.3% of global vegetated areas, such as Europe and the Amazon rainforest. We found that PRE_ CV exerted a strong positive influence on NDVI_ CV for 20.5% of global vegetated lands, whereas PRE_ mean exerted a strong negative influence for 19.1% of global vegetated lands. The influences from precipitation were scattered over both arid and humid regions. SWD_ CV exerted a comparable negative influence on NDVI_ CV (18.8% of global vegetated lands) with a positive influence (17.8% of global vegetated lands); but the SWD_ mean had more positive impacts on NDVI_ CV than negative impacts. Generally speaking, climate mean states had more extensive impacts on NDVI_ CV than did climate variability. In addition, all factors except PRE_ mean and SWD_ CV exerted more positive than negative impacts on NDVI_ CV .
Next, PLSR and PCR were used to determine which climate factors had the strongest effects on the IAV of vegetation greenness. The two regression methods generated comparable results and similar spatial patterns, as shown in figures 5 and S6. Temperature states exerted the most extensive influence on NDVI_ CV , accounting for 36.5% (indicated by PLSR) and 35.5% (indicated by PCR) of global vegetated areas with a stronger influence from TMP_ mean than from TMP_ CV . Dominant temperature impacts were observed in southern Asia, eastern and northern North America, northern and central South America, and central Africa. Solar radiation states presented the second largest influence on NDVI_ CV variations, accounting for 34% (by PLSR) and 34.3% (by PCR) of global vegetated lands. These regions were scattered across western Asia, western North America, southern South America, and western Australia. Slightly more regions were controlled by SWD_ mean than by SWD_ CV . PRE states dominating NDVI_ CV were identified in 29.5% (by PLSR) and 30.2% (by PCR) of global vegetated lands, including northwest Asia, northern and western Europe, northern North America, and southern Africa. We also determined the climatic factors that most influenced NDVI_ CV variations in each climate zone, as shown in figure S7. We observed large divergences in the factors controlling NDVI_ CV changes in different climate zones. TMP_ mean displayed a dominant positive correlation with NDVI_ CV changes in the equatorial climate zone, but a negative relationship in the snow climate zone. NDVI_ CV changes in the arid climate zone were driven by PRE_ mean , by SWD_ CV in the warm temperate zone, and by SWD_ mean in the polar climate zone.  lands of the arid and polar climate zones. This pattern of IAV of vegetation greenness is roughly consistent with previous studies (Zhao et al 2018). At the biome scale, the IAV of vegetation greenness in grasslands was highest, whereas variability in evergreen broadleaf forests was lowest. Furthermore, nearly all forests displayed relatively low variability, which agrees with previous research (Fang et al 2001, Knapp and Smith 2001). The primary reason for this low variability is that forests are more resistant to climate anomalies than other vegetation types (Hirota et al 2011). Without considering disturbances from human activities, vegetation variability should be determined primarily by climate conditions (Wu et al 2015). In the current study, we also calculated the spatial distributions of climate factor variability (PRE, TMP, and SWD), as shown in figure S8. Temperature variability presents a clear latitudinal distribution, gradually decreasing with decreasing latitude. The distribution of solar radiation variability has great spatial heterogeneity, with high variability in Europe, southern South America, and western Australia, and lowest in the equatorial regions. The distribution of precipitation variability is highly consistent with the pattern of vegetation variability, with high variability in drier regions and low variability in humid regions, suggesting that precipitation variability is a strong determinant of vegetation variability distribution. Previous studies have suggested a heavy influence of precipitation variability on ecosystem productivity in global pastures (Sloat et al 2018). In addition, precipitation extremes were found to determine the distribution of global vegetation extremes (Liu et al 2013, Ahlström et al 2015a).

Discussion
In addition to determining the spatial patterns of the IAV of vegetation greenness trends from 1982 to 2015, we also analyzed the dynamics of vegetation variability over the long term. We found increasing variability in global vegetation greenness, suggesting that the IAV of vegetation growth is increasing. This increase in variability of vegetation greenness was found to be mainly due to the increase in the standard deviation other than the mean value of NDVI, implying that the increase variability should not be attributed to the long-term greening over continents . . Therefore, the observed increase of vegetation variability would be expected to have profound impacts on global carbon cycle. As previous studies have found, vegetation in arid zones and polar regions were more susceptible to disturbances from the external environment (IPCC 2007, Lioubimtseva andHenebry 2009). This also explains the increase in extreme events such as droughts and heat waves in recent years (Perkins et al 2012). Propastin et al (2010) found that the vegetation classes 'grassland' and 'closed shrubland' exhibited the largest proportion of land area containing both moderate and high vulnerability to the external environment. Similarly, our study demonstrates that grassland and shrub variability increased most severely. Propastin et al (2010) also suggested that forest areas were significantly less affected by climate than the other land cover types. The current study confirms these findings by showing that changes in the IAV of growth of evergreen broadleaf forests and deciduous broadleaf forests were the most stable of the vegetation types. Furthermore, despite large uncertainties, satellite NDVI measurements and multi-model estimates of GPP demonstrated a consistent spatial distribution of vegetation variability over the last 34 years. However, we also noted that these TRENDY models did not fully capture the change trend of vegetation variability as NDVI from 1982 to 2015, suggesting a need for improved ecological models. A recent study also suggested that current ecological models underestimate the IAV of global ecosystem production associated with water availability conditions (Humphrey et al 2018).
Numerous factors can influence the IAV of vegetation greenness, including fire, deforestation, and land use changes but large scale changes are driven by climatic factors (Osborne and Wheeler 2013). This study reveals the complex influence of climatic factors on the IAV of vegetation greenness dynamics. The most influential impacts on the IAV of vegetation greenness were from temperature, followed by solar radiation and precipitation. Vegetation growth in the northern high latitudes is usually limited by temperature (Nemani et al 2003, Seddon et al 2016). Interestingly, we observed strong impacts from precipitation and solar radiation on changes in the IAV of vegetation greenness. Piao et al (2014) reported a weakening influence of vegetation change on the IAV of vegetation growth in high latitudes due to disturbances from drought and heatwaves. In addition, warming has caused spring to come early in these regions, which may result in increased water demand by plants during the summer months (Buermann et al 2013). Hence, precipitation would be expected to exert a stronger influence on the IAV of annual maximum NDVI. We also observed a dominant influence of temperature on changes in the IAV of vegetation growth in areas where radiation is traditionally thought to control vegetation growth such as tropical regions. Previous studies have found a strong coupling between tropical temperature changes and changes in the IAV of atmospheric CO 2 growth rate, and that this coupling was primarily regulated by the influences of temperature changes on tropical ecosystem productivity (Wang et al 2013, Wang et al 2014. High temperature tends to promote the vapor pressure deficits, which will induce the closure of plant stomata and thereby the reduction of photosynthesis rate and productivity ( . High temperatures can trigger drought conditions by increasing vapor pressure deficits and increasing evapotranspiration even in normal rainfall conditions (Corlett 2011, Silva et al 2013. In some arid regions, such as southern South America, southern Africa, and Australia, precipitation and solar radiation exert a blended influence on changes in the IAV of vegetation greenness. It is easy to understand the large impact of precipitation on changes in the IAV of vegetation growth because water is the main limiting factor for vegetation growth in these regions (Nemani et al 2003). A possible explanation for the observed impact of solar radiation on changes in the IAV of vegetation growth in these regions may be that radiation is a primary energy component that influences soil water conditions by changing atmospheric water demand (Dai 2011).
Despite the long-term greening of Earth's surface, the IAV of vegetation greenness is also increasing, which suggests that the stability of ecosystems is decreasing. This change in the IAV of vegetation greenness should have profound impacts on the dynamics of ecosystem carbon sinks. Considering the strong coupling between ecosystem productivity and atmospheric CO 2 concentrations, the atmospheric CO 2 growth rate would also likely be affected (Richardson et al 2010, Graven et al 2013, Piao et al 2017. It should be noted that there are some limitations to this study. First, the uncertainties associated with remote sensing data cannot be ignored. For example, the NDVI has been suggested to suffer from saturation problem in areas with dense vegetation (Morton et al 2014). In addition, it may also contain errors produced by the shift or degradation of sensors and by the noise signals from atmosphere and ground (Atzberger et al 2013, Jiang et al 2013, Zeng et al 2013. Second, this study provides only a preliminary investigation of the climatic drivers of changes in the IAV of vegetation greenness, but the mechanisms of these drivers is still unclear. Although this study focused on the role of climate, other factors that influence the IAV of vegetation growth such as fires, land use and CO 2 fertilization, are also important. Therefore, further analysis involving these factors is needed to fully understand the mechanisms behind the IAV of vegetation growth.

Conclusions
This study provides insights into global vegetation variability using satellite-based estimates of vegetation greenness for the period 1982-2015. On a global scale, we observed a high correlation between the spatial patterns for precipitation variability and vegetation variability, with high precipitation variability corresponding to high vegetation variability. In addition, this study reveals increasing IAV in annual maximum vegetation in more than 45% of vegetated areas globally, whereas only 21% exhibited decreasing IAV. Regression analysis suggests that temperature changes explain the largest proportion of changes in the IAV of vegetation greenness, followed by solar radiation and precipitation.
All codes for data processing can be obtained by contacting the corresponding author.