Global vegetation greenness interannual variability and its evolvement in recent decades

The interannual variability (IAV) of global vegetation greenness needs careful assessment as it relates to the stability of the climate, conservation of biodiversity, sustainable ecosystem services, and global food security. Here, we investigated the spatial feature and temporal evolvement of global vegetation greenness interannual variability from 1982 to 2015 using the Global Inventory Modeling and Mapping Studies (GIMMS) Normalized Difference Vegetation Index third generation (NDVI3g) data. Generally, regions with herbaceous and short woody plants had larger IAV of vegetation greenness than those with tall woody biome types (7.9% versus 2.9%). On average, all the biomes displayed increasing IAV of vegetation greenness from 1982 to 2015, with notable increases over northern high latitudes (0.135%/year), Eastern Europe (0.037%//year), and Central Australia (0.231%/year). Croplands in China and India experienced decreasing IAV of vegetation greenness (–0.037%/year for China and –0.004%/year for India). The changing IAV of vegetation greenness had implications for climatic, environmental, and anthropogenic changes that influence vegetation dynamics. Some note-worthy factors include climate warming, the CO2 fertilization effect, agricultural practice improvement, cropland abandonment, and China’s Grain-for-Green Program.


Introduction
The spatial and temporal characteristics of global terrestrial vegetation dynamics are of great importance to the Earth's climate system (Sellers et al 1996), with implications and consequences for the global carbon cycle (Quéré et al 2016, Piao et al 2022, water cycle (Wang and Dickinson 2012), energy balance (Schaaf et al 2002), land use land cover change (Hansen et al 2013, Song et al 2018, and ecosystem services (Costanza et al 1997). There are two categories of vegetation dynamics: intra-annual and interannual. Satellite remote sensing (Friedl et al 2002, Huete et al 2002, Pinzon and Tucker 2014 is indispensable for monitoring both the intra-annual (Zhang et al 2006, Bi et al 2015 and interannual vegetation dynamics (Myneni et al 1997, Hilker et al 2014, Poulter et al 2014, Mao et al 2016, Chen et al 2019b, Berner and Goetz 2022 at regional, continental, and global scales. Regions of substantial changes in vegetation dynamics include northern high latitudes (Goetz et al 2005, Fraser et al 2014, Piao et al 2017, ecological restoration regions (Feng et al 2016, Song et al 2018, Chen et al 2019b, agricultural practice improving areas (Chen et al 2019b), abandoned croplands (Kolecka 2021), and tropical forests (Hansen et al 2013, Hilker et al 2014, Zhang et al 2016. The multi-decadal trends of global vegetation dynamics and their drivers have been investigated extensively (Myneni et al 1997, Nemani et al 2003, Chen et al 2019b, Piao et al 2020, Berner and Goetz 2022, Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. Shen et al 2022). With more frequent extreme climatic events (i.e., droughts, heat waves, fires, and storms) (Fischer and Knutti 2014, Jolly et al 2015, Chang et al 2020, Chiang et al 2021, Power et al 2021, several researchers have also looked into the vegetation interannual variability (IAV) and how it responds to the changing climate. Jiang et al (2017) compared four long-term Leaf Area Index (LAI) data sets between 1982 and 2011 and found that Global Land Surface Satellite (GLASS) LAI betrayed the largest interannual variability of global LAI, almost 4.5 times that of the GLOBMAP LAI. Chen et al (2019a) investigated the global vegetation greenness interannual variability and the trends of the variability from 1982 to 2015 using the Global Inventory Modeling and Mapping Studies (GIMMS) NDVI3g data and found different spatial patterns of global vegetation interannual variability than those in (Jiang et al 2017). Furthermore, Chen et al (2019a) found increasing interannual variability of vegetation greenness in more than 45% of vegetated areas globally, and only 21% of vegetated land exhibited decreasing interannual variability of vegetation greenness.
There are many ways to define annual vegetation greenness: annual maximum greenness (Chen et al 2019a), growing season mean greenness (Myneni et al 1997), growing season integrated greenness , and annual mean greenness (Chen et al 2019b). Each of them has both strengths and limitations. Annual maximum greenness can eliminate the impact of atmospheric corruption most thoroughly but cannot incorporate vegetation greenness all year round. Growing season mean and growing season integrated greenness describes vegetation photosynthetic activity most directly, but the definition of growing seasons at the global scale can be tricky. Annual mean greenness can incorporate vegetation greenness all year round, which is suitable for globalscale studies; nevertheless, it might include some non-vegetation signals.
The discrepant magnitudes and spatial extent of global vegetation interannual variability found in (Jiang et al 2017) and (Chen et al 2019a) could be due to the following reasons: different time spans of the study periods, different satellite sensors, different atmospheric correction algorithms, different vegetation parameter retrieval algorithms, different annual vegetation greenness calculating procedures, and different interannual variability assessment methods. Jiang et al (2017) used growing season mean LAI, while Chen et al (2019a) used annual maximum NDVI as well as growing season mean NDVI and found the two variables generated similar results, albeit with differences in some regions. Jiang et al (2017) used the standard deviation of the detrended time series of annual LAI to assess vegetation greenness interannual variability, which accounted for the long-term vegetation greenness trends. Still, it hindered the inter-region or inter-biome comparability of the calculated variability due to the different levels of vegetation greenness among land cover types. Chen et al (2019a) used a metric called Coefficient of Variation (CV) to assess the vegetation greenness interannual variability, which suits the inter-region and inter-biome comparison because the metric is calculated with normalization by the annual mean vegetation greenness (equation (1)).
where CV is the Coefficient of Variation, SD is the standard deviation, and Mean the mean of the time series. Still, this metric is meddled by long-term trends of vegetation greenness. A better method would be to combine the merits of both approaches used in (Jiang et al 2017) and (Chen et al 2019a).
Considering the limitations of the previous research, we assessed the interannual variability of global vegetation greenness with improved methods in order to further assess the global vegetation greenness interannual variability from a different perspective. We also explored factors that drove the spatial and temporal characteristics of the interannual variability of vegetation greenness. Our study aims to investigate: (1) the spatial distribution of the interannual variability of global vegetation greenness from 1982 to 2015; (2) the temporal evolvement of the interannual variability of global vegetation greenness in recent decades, especially that of biome types directly influenced by human activities, such as croplands; and (3) the climatic, environmental, and anthropogenic causes of the spatial and temporal features of vegetation greenness interannual variability.

Data and methods
2.1. Global NDVI data and land cover data We used the GIMMS NDVI3g data (Tucker et al 2005, Pinzon and in our investigation. The GIMMS NDVI3g data were derived from the remote observations of a series of Advanced Very High Resolution Radiometer (AVHRR) sensors onboard the National Oceanic and Atmospheric Administration (NOAA) platforms. We got the GIMMS NDVI3g data from the Big Earth Data Platform for Three Poles. The temporal resolution of the data is half monthly, the spatial resolution is 1/12°, the time span is from 1982 to 2015, and there are cloud and snow contamination flags in the vegetation index data. The data have been widely used in global and regional vegetation dynamics studies (Piao et al 2020).
We used the International Geosphere-Biosphere Program (IGBP) land cover classification data from the Collection 6 MODIS land cover product MCD12C1 (Sulla-Menashe et al 2019) to identify the global land cover types at the spatial resolution of 0.05°(figure S1). We excluded the non-vegetated land cover types in our investigation.

Defining the annual vegetation greenness
We identified the snow-contaminated pixels in the GIMMS NDVI3g data as missing data, filled the missing values using the Savitzky-Golay filter, and resampled the data to 0.05°using the nearest neighbor method to match the spatial resolution of the MODIS land cover type data. Chen et al (2019a) used the annual maximum vegetation greenness index to indicate the yearly vegetation status (i.e., the peak growing season vegetation status). Annual maximum NDVI can neither fully incorporate the vegetation greenness dynamics all year round nor reflect long-term vegetation phenological changes, which are not negligible in the era of climate change. In contrast, the annual mean vegetation greenness incorporates the vegetation dynamics all year round, which means that the changes in phenology are implied in annual mean vegetation greenness (Chen et al 2019b). We defined annual vegetation greenness as an indicator that is representative of vegetation greenness all year round, including all the phenological stages, so that the changes in vegetation greenness in any part of the year can be reflected in our definition of annual vegetation greenness. Hence, we used annual mean NDVI to indicate the annual vegetation greenness level. We calculated the annual vegetation greenness by averaging the 12 monthly NDVI (the average of the two half-monthly NDVI in a month). The spatial patterns of multi-year annual mean NDVI from 1982 to 2015 for the global vegetated lands are as in figure S2.

Calculation of the interannual variability
The Interannual Variability (IAV) was calculated according to equation (2).
where IAV is the interannual variability, SD detrended is the standard deviation of the detrended time series data, in which the trend was estimated using the Ordinary Least Squares method, and Mean is the long-term mean of the time series data. Our method here differs from the methods used in (Jiang et al 2017) and (Chen et al 2019a); in essence, we combined the merits of both studies. Jiang et al (2017) accounted for the long-term trends of time series data by detrending, and Chen et al (2019a) made the metric of variability comparable among regions/ biomes by normalization with the long-term mean of the time series data. Equation (2) not only excludes the impacts of the long-term trends of time series but also normalizes the variability so that IAV is comparable among biome and land cover types.

Detecting the changes in vegetation greenness interannual variability
For each 0.05°vegetated grid identified by the MODIS land cover product, we constructed a moving window of 15 years for a total of 34 years of NDVI data from 1982-2015, calculated IAV for each 15-year moving window, and assigned the calculated IAV to the middle year of each 15-year window. Hence, nominally, we got 20 IAVs for each vegetated pixel from 1989 to 2008. We calculated the trend of the 20 IAV time series using the Ordinary Least Squares (OLS) method and tested the significance with the t-test at the significance level of 0.05 (α = 0.05).
In this way, we assessed the changes in vegetation greenness interannual variability as well as the implications and impacts. We also discussed the climatic, environmental, and anthropogenic implications of the changing IAVs.

Relationship between amplitudes of vegetation greenness interannual variability and the corresponding trends
The IAV of vegetation greenness was related to the corresponding vegetation greenness IAV trends ( figure 3). The spatial distribution of the combination of the amplitudes of vegetation greenness interannual variability and the corresponding trends is as in figure 3(a). Northern high latitudes and semi-arid regions in Central Australia were dominated by vegetation greenness IAVs above the global average and their increasing trends, suggesting the vegetation greenness IAVs in these regions tend to get larger in the era of climate change (figure S3). Southeastern Asian tropical rainforests and the monsoon regions in Asia were dominated by vegetation greenness IAVs below the global average and decreasing trends. Unlike northern high latitudes and semi-arid regions, Southeast Asian tropical forests were getting more stable, opposite to Congolese rainforests (figure S4). Croplands in Eastern Europe and Northeastern China had vegetation greenness IAVs below the global average and increasing trends, suggesting the croplands in these regions were becoming less and less stable or were abandoned (Kolecka 2021), implying increasing climatic risk for these agricultural lands  future warming and limited adaptation potential for wheat and barley in Europe, which agrees with our results for croplands in Europe ( figure S5).
The correlation coefficient of the latitudinal profiles for vegetation greenness IAV and its trends was 0.87 and statistically significant (α = 0.05) ( figure 3(b)). The relationship between vegetation greenness IAV and its trends for the 12 global biome types is as in figure 3(c). Biome types with larger IAVs tend to become more responsive to climatic variations in the era of climate change. All the biome types except Evergreen Broadleaf Forests, on average experienced increasing vegetation greenness IAV in the past decades.

Comparison with previous vegetation greenness IAV studies
The spatial pattern of global vegetation greenness IAV in this study ( figure 1(a)) differs from that in (Chen et al 2019a). In all the regions with high vegetation greenness IAV except northern Eurasia (e.g., Australia, central Asia, western North America), the vegetation greenness IAV in (Chen et al 2019a) is higher than in this study because (Chen et al 2019a) did not detrend the time series of annual vegetation greenness, hence the IAV is higher due to the increasing trends of vegetation greenness. As for northern Eurasia, the method in (Chen et al 2019a) did not consider the substantial variation of growing season length in northern high latitudes (Potter and Alexander 2020); hence the estimated IAV of vegetation greenness is lower in (Chen et al 2019a) than in this study. Because we used annual mean NDVI to represent the annual vegetation greenness, the variation of growing season length was incorporated into the annual vegetation greenness metric. Biased estimation of vegetation greenness IAV could conflict with conclusions based on other lines of evidence, such as carbon budgeting, atmospheric data, and forest inventories (Walker et al 2021).
The spatial pattern of the trends of global vegetation greenness IAV in this study agrees with that in (Chen et al 2019a) for most parts of the world. For example, open shrublands have the largest range of IAV, suggesting that they are most heterogeneous globally, probably due to the varying climatic factors or species composition. Open shrublands generally experienced the strongest increases in vegetation greenness IAV, which might be due to the expansion of shrubland in northern high latitudes due to climate warming (Fraser et al 2014), conversion of deciduous needleleaf forests to evergreen needleleaf forests in Siberia (He et al 2017), as well as the increase in leaf area in central Australia caused by increasing atmospheric CO 2 concentrations (Zhang et al 2022). Open shrublands distribute in many regions of the world with varying climate conditions, such as northern high latitudes, central Australia, southwest United States, etc (figure S1), with these distinctly different climates and species compositions, open shrublands also exhibited large variations in the trends of vegetation greenness IAV.
The open shrublands in Australia and northern high latitudes experienced notable increases in IAV of vegetation greenness, which might be related to the expansion of shrublands in these regions (McManus et al 2012, Frost et al 2014, Saintilan and Rogers 2015. The Southeast Asian tropical forests experienced widespread decreases in IAV of vegetation greenness, the tropical forests in Congo experienced increases in IAV of vegetation greenness, and the Amazonian forests displayed mixed trends in IAV of vegetation greenness ( figure  S4). The discrepant trends of tropical forests' greenness IAV need further investigation in particular, because of their important roles in the global carbon cycle. The increases in the vegetation greenness IAV of croplands in Europe, Asia, and North America suggest the potential regional vulnerability of food security (figure S5) (Kolecka 2021, Ortiz-Bobea et al 2021. Large portions of croplands in China and India experienced significant decreases in cropland greenness IAV, mainly due to improved agricultural management practices such as irrigation and fertilizer use in the two developing countries (Chen et al 2019b). China and India have increased the crop harvest frequency (Ray and Foley 2013), and India has increased the area of irrigated croplands (Ambika et al 2016).
Croplands in Eastern Europe and Northeastern China had vegetation greenness IAVs below the global average and increasing trends, suggesting the croplands in these regions were becoming less and less stable or were abandoned (Kolecka 2021), implying increasing climatic risk for these agricultural lands (Lobell et al 2011) (figure S5). For example, (Moore and Lobell 2014) found high adaptation potential for maize to future warming and limited adaptation potential for wheat and barley in Europe, which agrees with our results for croplands in Europe ( figure S5). However, for the croplands in China, India, and western Europe, (Chen et al 2019a) detected increasing trends of vegetation greenness IAV, rather than decreasing trends. This contradicts the fact that the agricultural practices in these regions improved in recent decades, especially in China and India, which should help to stabilize agricultural production against climate disturbances (Chen et al 2019b). Again, the increasing IAV of croplands in (Chen et al 2019a) could be due to the lack of consideration of vegetation greenness dynamics out of the peak growing season and the lack of consideration of the evolving cropland phenology. Because we used annual mean NDVI to represent the annual vegetation greenness, the evolving cropland phenology was incorporated into the annual vegetation greenness metric. In (Chen et al 2019a), the trends of vegetation greenness IAV showed more extreme positive or negative trends when using annual maximum NDVI than when using growing season mean NDVI, which implies that vegetation greenness of peak growing seasons is more interannually variational than that of whole growing seasons.
We used the annual mean NDVI to represent annual vegetation greenness, which is not indicative of peak growing season vegetation greenness. To learn more about the IAV of global vegetation peak growing season greenness, readers can refer to this research (Chen et al 2019a).

The climatic driving factor
The air temperature was one of the major driving forces of vegetation dynamics, particularly over northern high latitudes. The vegetation dynamics of the northern high latitudes have been investigated extensively (Myneni et al 1997, Macias-Fauria et al 2012, Bi et al 2013, Ju and Masek 2016, Welp et al 2016, Berner et al 2020, Liu et al 2021a, 2021b in the era of climate change because of the amplified warming there . Macias-Fauria et al (2012) suggested Eurasian Arctic greening revealed the potential for structurally novel ecosystems. The increasing interannual variability in Asian Arctic vegetation greenness could imply that the shrub genera were shifting toward more sensitive ones to climate variability (Liu et al 2021a(Liu et al , 2021b. Besides air temperature, precipitation was another important climatic factor affecting vegetation greenness IAV. Al-Yaari et al (2021) investigated the interannual variability of biomass over the contiguous United States and found the greatest interannual variability in shrubs and grasslands, with forests the least variable; and the IAV of biomass was strongly correlated with annual precipitation. This agrees with Peng et al.'s results in China (Peng et al 2021) and agrees with our results (figure 1). The large range of grassland greenness IAV found in our study ( figure 1(b)) suggests the wide range of precipitation characteristics over grasslands.

The environmental driving factor
Biodiversity plays an important role in the vegetation greenness IAV of terrestrial ecosystems. Oliveira et al (2022) found biodiversity mediates ecosystem sensitivity to climate variability-regions with greater plant diversity exhibit lower sensitivity to temperature variability at the interannual and seasonal scales and lower sensitivity to precipitation variability at the interannual scales. The increasing interannual variability of vegetation greenness might be related to decreasing biodiversity, particularly in China's Grain-for-Green Program (GFGP) regions, where GFGP forests are mostly monocultures or compositionally simple mixed forests (Hua et al 2016).
The increasing atmospheric CO 2 concentration is another important environmental factor driving increasing global vegetation greenness IAV. Zhang et al (2022) found increasing sensitivity of dryland vegetation greenness to precipitation because of the increase of leaf area in global drylands caused by increasing atmospheric CO 2 (i.e., the CO 2 fertilization effect). This mechanism could explain the increasing vegetation greenness IAV we found in arid and semi-arid regions (figure 2).

The anthropogenic driving factor
Anthropogenic factors, such as armed conflicts, also contribute to the changes in vegetation greenness IAV. Tao et al (2022) found increasing and widespread vulnerability of intact tropical rainforests in the Americas, Africa, and Asia, suggesting the capacity of intact rainforests to withstand future droughts is limited. This agrees with our results, especially in the Congolese rainforests. However, our results in Indonesia of Southeast Asia suggest an increase in the resilience of rainforests there. The divergent trends in the resilience of tropical forests in Africa and Southeast Asia mandate further investigation. A possible reason might be anthropogenic. The anarchic conditions in the Congo basin might be the reason for the deteriorating condition of Congolese rainforests (Tamm 2019).
Cropland abandonment also can contribute to the changes in vegetation greenness IAV and, at the same time, jeopardize food security. Asseng et al (2015) suggested that wheat yield could become more variable with rising temperatures, and this agrees with the increasing vegetation greenness IAV of wheat planting croplands in Eastern Europe. The increasing IAV of vegetation greenness in Eastern Europe could be caused by agricultural land abandonment (Kolecka 2021). With no agricultural management, the croplands converted to grasslands and shrublands, which were more sensitive to climate variability, and had more vegetation greenness IAV than croplands ( figure 1(b)). Challinor et al (2014) not only projected the decrease in crop yields in the second half of the 21st century due to climate change but also projected the increase in IAV of crop yields, which might be more noteworthy since it relates to the steadiness of food supply. The potential of agricultural management for increasing food supply is not limitless, and when agrarian management is not adequate to deal with more extreme climate variability, the IAV of cropland greenness will probably increase accordingly and jeopardize food security. Yue et al (2020) found that managed land accounted for 30%-45% of the IAV of global land carbon sink from 1959 to 2015. We found that the vegetation greenness IAV of croplands in Eastern Europe and Northeastern China had increased from 1989 to 2008. The increase of vegetation greenness IAV in Eastern Europe was probably due to cropland abandonment (Kolecka 2021); the increase of vegetation greenness IAV in Northeastern China might be caused by the conversion of croplands to shrublands (Yu et al 2021).
The Chinese ecological restoration project, such as the Grain for Green program (Geng et al 2020) was another cause for the increase in vegetation greenness IAV. When croplands converted to shrublands, because the vegetation greenness IAV of shrublands is greater than that of croplands ( figure 1(b)), the ecological restoration regions experienced increasing vegetation greenness IAV. And the mechanisms proposed in (Zhang et al 2022) could explain the continuously increasing vegetation greenness IAV over the Chinese ecological restoration regions-the CO 2 fertilization effect leads to increased water use efficiency of plants in semi-arid regions (Fatichi et al 2016), which decreases the sensitivity of vegetation dynamics to variation of precipitation.

Conclusion
We investigated the spatial feature and temporal evolvement of global vegetation greenness interannual variability in recent decades from 1982 to 2015 using the GIMMS NDVI3g data. Generally, regions with herbaceous biome types had larger IAV of vegetation greenness than those with tall woody biome types. It's note-worthy that croplands had smaller IAV of vegetation greenness than shrublands and grasslands, probably due to agricultural management practices. On average, all the biomes displayed increasing IAV of vegetation greenness from 1982 to 2015, with prominent increases over northern high latitudes, Eastern Europe, and Central Australia. These changes had implications for climatic, environmental, and anthropogenic changes, which influence vegetation dynamics. Some note-worthy changes include climate warming, the CO 2 fertilization effect, agricultural practice improvement, cropland abandonment, and China's Grain for Green Program. In particular, large portions of croplands in China and India experienced significant decreases in cropland greenness IAV, mainly due to improved agricultural management practices such as irrigation and fertilizer use in the two developing countries. The increasing IAV of vegetation greenness in China's Grain-for-Green Program (GFGP) regions was probably related to decreasing biodiversity, because GFGP forests were mostly monocultures or compositionally simple mixed forests.

Data availability statement
No new data were created or analysed in this study.

Declaration of Interests
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Funding
This work was financially supported by the Chinese Academy of Engineering (Grant No. 2022-XY-139) and Fundamental Research Funds for the Central Universities (Grant No. lzujbky-2022-sp13).