2.1 Satellite Remote Sensing Data

We used six long-term and three short-term satellite remote sensing datasets in this study, covering four types of vegetation indicators - NDVI, LAI, SIF and VOD (Table 1). These remote sensing products include: (a) VIP15 NDVI product (1981–2014), which has 0.05° and 15-day resolutions and is developed by harmonizing the observations of Advanced Very High Resolution Radiometer (AVHRR) from 1981 to 1999 and Moderate Resolution Imaging Spectroradiometer (MODIS) C5 from 2000 to 2014 (Didan et al., 2015); (b) GIMMS NDVI3g (1981–2015), which has 15-day and 1/12° resolutions, and was produced by aggregating daily AVHRR surface reflectance (Pinzon & Tucker, 2014); (c) GIMMS LAI3g product (1981–2015), which was further produced by GIMMS NDVI3g product using a neural network algorithm (Zhu et al., 2013); (d) PKU GIMMS NDVI (1982–2020), which is a new version of GIMMS NDVI product produced by a machine learning model incorporating Landsat images (Li et al., 2023); (e) GLASS LAI product (1981–2018), which has 8-day temporal resolution and 0.05° spatial resolution, and was reconstructed by combing AVHRR LAI from 1981 to 1999 and MODIS LAI from 2000 to 2018 using a bidirectional long short-term memory (Bi-LSTM) model (Ma & Liang, 2022); (f) GLOBMAP LAI (1982–2019) dataset at a spatial resolution of ∼0.07°, covering the period from 1982 to 2019, has half-month (1982–2000) and 8-day (2001–2019) temporal resolutions, and was produced by establishing a pixel-level AVHRR Simple Ratio (SR)-MODIS LAI relationship (Liu et al., 2012); (g) MOD13C1 NDVI (2000–2020) product (C61), which has a 0.05° spatial resolution and a 16-day temporal resolution (Didan & Munoz, 2019); (h) OCO2 SIF product (2000–2020) at resolutions of 0.05° and 4 days, which was generated from MODIS surface reflectance and OCO2 SIF data using a neural network approach (Zhang et al., 2018); and (i) VOD (1987–2017) dataset, which was produced by merging observations from multiple microwave sensors at daily temporal resolution and 0.25° spatial resolution (Moesinger et al., 2020). We standardized these satellite remote sensing products to a 0.5° spatial resolution by using pixel aggregation (PA) method and to a monthly temporal interval by using maximum value composite (MVC) method (Ma et al., 2022; Tian et al., 2021). Considering the high similarity among GIMMS-version datasets (Figure S1 in Supporting Information S1), we only used GIMMS NDVI3g, along with other three independent long-term satellite remote sensing products (i.e., VIP15 NDVI, GLASS LAI and GLOBMAP LAI) for the main analysis.

2.2 LAI Estimates

We used monthly LAI estimates of 12 DGVMs from TRENDY v9 (Friedlingstein et al., 2022; Sitch et al., 2015) to check the consistency of the changes in the vegetation IAV derived from satellite remote sensing products and model measurements. The 12 DGVMs include CABLE, CLASSIC, CLM5, ISAM, ISBA, JULES, LPJ, LPX, ORCHIDEE, SDGVM, VISIT and YIBs (Table S1 in Supporting Information S1), which were driven by monthly CRU or 6-hourly CRU-JRA gridded climate datasets, as well as dynamic atmospheric CO₂ concentrations. These models provide LAI estimates under four different scenarios, namely no change (S0), varying CO₂ only (S1), varying CO₂ and climate (S2), and varying CO₂, climate and land cover change (S3). We analyzed CV trends from LAI estimates of S3 scenario in this study. To match the period of data availability with satellite remote sensing observations, we used LAI estimates from 1982 to 2014.

2.3 Quantifying the Vegetation IAV Changes

We used the coefficient of variation (i.e., CV), the ratio of the standard deviation to the mean, to indicate vegetation IAV for each dataset. The use of CV is meant to remove the differences in magnitude between datasets (i.e., NDVI, LAI, SIF and VOD) to ensure the comparability among the results. To further quantify the CV changes, we first defined the growing season as the period when the mean daily air temperature is above zero (Jiang et al., 2017; Smith et al., 2019) for each vegetated pixel (Figure S2 in Supporting Information S1). We adopted a “methodological growing season” (Körner et al., 2023) as it allows for easy extraction and reproducibility without incurring uncertainties from more complex definitions of phenology. Meanwhile, this definition effectively eliminates the frozen period when there is no vegetation growth, and aligns well with the phenological dates extracted by MODIS NDVI (Leeper et al., 2021). We obtained yearly composites by summing all values within the growing season (Piao et al., 2020; Zhu et al., 2016). Subsequently, we calculated the CV for every 10-year moving window and assigned the CV value to the middle year of the window. The Theil-Sen method was used to estimate the CV trend (i.e., slope) (Wang et al., 2019), and its significance (i.e., p value) was determined using a two-tailed Student’s t-test (Jiang et al., 2019; Xu et al., 2022). Lastly, the non-parametric Mann-Kendall test was used to detect whether a significant monotonic increasing or decreasing trend exists (Jiang et al., 2019).

2.4 Classifying the Levels of Consistency of Remote Sensing-Based CV Trends

Considering potential inconsistency of CV trends across four long-term satellite remote sensing products (i.e., VIP15 NDVI, GIMMS NDVI3g, GLASS LAI and GLOBMAP LAI), we classified the levels of consistency of CV trends for each vegetated pixel using the following criterion (Kause et al., 2022; Xu et al., 2022): (a) “virtually certain” (CER) means the signs of CV trends derived from all four long-term satellite products are the same and significant (Mann-Kendall test, p ≤ 0.1); (b) “likely” (LIK) if the signs of the CV trends are the same and significant in any three satellite products; (c) “about likely as not” (ALN) if the signs of the CV trends are the same and significant in any two satellite products; (d) “possibly” (POS) if only one satellite product yields significant CV trend, and others are insignificant; (e) “no change” (NOC) if no significant CV changes were detected in all four satellite products; and (f) “conflicting” (CON) if the observed CV trends are conflicting with each other (i.e., significant positive and negative CV trends were detected simultaneously across the four products). We further assigned “+” and “−” signs to the consistency levels that refer to the direction of CV trends (i.e., positive and negative trends, respectively). Overall, the levels of consistency were classified into four types, positive, that is, CER (+), LIK (+), ALN (+) and POS (+), negative, that is, CER (−), LIK (−), ALN (−) and POS (−), no change (NOC), and conflicting (CON), respectively. Based on previous studies, NDVI and LAI are both indicators of canopy structure and greenness, and they are strongly related (Wang et al., 2022; Zeng, Hao, et al., 2022). Therefore, we deem that there is a clear physiological and statistical basis to treat them equally. In addition, we used both indicators to augment the number of available long-term datasets for our analysis, as some datasets (e.g., VIP15) only contain NDVI, while others (e.g., GLOBMAP) only include LAI.

2.5 Identifying Factors Driving the IAV Changes

Based on the findings from previous studies (Baldocchi et al., 2016; Berdugo et al., 2022; Forkel et al., 2016; Gray et al., 2014; Luo & Keenan, 2022; Piao et al., 2015; Zhu et al., 2016), we investigated several factors that may serve as potential drivers for the IAV changes. These factors are aridity index (AI), mean annual air temperature (MAT), mean annual precipitation (MAP), land use and land cover change (LUCC), mean monthly LAI, the trend of LAI (LAI_trend) and nitrogen deposition (N deposition). AI was obtained from the Global Aridity Index dataset at a 30 arc-second resolution from 1970 to 2000 (Zomer et al., 2022), in which lower AI values mean drier conditions. MAT and MAP data were extracted from the gridded CRU-JRA V2.1 dataset (6-hr and 0.5° resolutions) (Harris et al., 2020). LUCC intensity was obtained using a 300-m ESA CCI land cover type product with 22 land cover classes from 1992 to 2019. We first determined whether land cover change occurred in each 300-m pixel from 1992 to 2019, and then quantified the LUCC intensity as the proportions of land cover change in each 0.5° spatial grid. The mean annual LAI was the average of GLASS LAI and GLOBMAP LAI from 1982 to 2014, and the LAI_trend (i.e., slope) was estimated using Mann-Kendall and Theil-Sen method for each vegetated pixel. Most global vegetated regions show significant greening trends (i.e., LAI_trend > 0) over the past decades (Chen, Chi, et al., 2019; Zhu et al., 2016). This implies that terrestrial ecosystems have accumulated more green leaves for carbon fluxes and in response to climate change, which may induce IAV changes. The N deposition was the average of the nitrogen deposition dataset (2° × 2.5° grid resolution) from 1984 to 2016 (Ackerman et al., 2019). All the above datasets were resampled into a 0.5° spatial grid.

We first assessed the significance of each factor in driving the vegetation IAV by using a two-sample t-test method, in which we examined whether the high and low values of each factor caused statistically different CV trends. For each factor, we used different criteria for the binary classification of high and low values, that is, the 50% for LUCC, the zero for LAI_trend, and the global average for other drivers. After that, we selected the four factors (i.e., AI, MAT, LAI_trend and N deposition; see Results) that have significant impacts on the CV trend and quantified their respective impacts using a multiple regression method: