DatasetsVegetation cover data<\/p>\n
The Vegetation Continuous Fields version 1 product (VCF5KYR) was adopted to provide global fine-scale fractional vegetation cover information14<\/a>. This product has a 0.05\u00b0\u2009\u00d7\u20090.05\u00b0 spatial resolution and a yearly interval, combining the advantages of multiple satellite sensors with very high spatial resolutions. The product is available from 1982 to 2016 and features the percentage of tree cover (tall vegetation \u2265\u20095\u2009m in height, abbreviated as \\({{\\rm{TC}}}\\)), short vegetation (\\({{\\rm{SV}}}\\)), and bare ground (\\({{\\rm{BG}}}\\)), respectively for land pixels, with a precision of 1%, representing the vegetation composition of each pixel during the peak growing season14<\/a>. The sum of the proportions of the three land-cover types for each pixel is 100%. The spatial pattern of global vegetation cover and the dominant changes in vegetation cover type are shown in Supplementary Fig.\u00a03<\/a>. Note that the years 1994 and 2000 were excluded from analyses due to a lack of VCF5KYR data (https:\/\/lpdaac.usgs.gov\/products\/vcf5kyrv001\/<\/a>). The broad vegetation cover classification was adopted in order to maximize the utility of the available data while ensuring the statistical significance and robustness of the results. Despite the inherent resolution limitation, our approach provides valuable information on VCC at a scale that is consistent with the LAI data.<\/p>\n LAI data<\/p>\n Our analysis involved the examination of four LAI products covering the satellite era, including the Advanced Very High-Resolution Radiometer (AVHRR) Global Inventory Monitoring and Modeling Studies (GIMMS) LAI3g dataset3<\/a>,49<\/a>, the National Oceanic and Atmospheric Administration (NOAA) Climate Data Record (CDR)50<\/a>, the Global Land Surface Satellite (GLASS) AVHRR LAI V50 product51<\/a>, and the GLOBMAP LAI V3 (ref52<\/a>.), all of which aim to provide consistent and reliable attributions of global LAI change. Note that while LAI serves as a valuable proxy of vegetation greenness11<\/a>,53<\/a>, it may be subject to saturation effects in areas with dense vegetation, similar to other vegetation indices54<\/a>,55<\/a>. In addition, estimating LAI from satellite imagery comes with several well-known limitations2<\/a>,56<\/a>,57<\/a>,58<\/a>,59<\/a>. These include spatial resolution constraints, saturation in densely vegetated areas, orbital drift, sensor replacements, retrieval algorithm uncertainties, and disturbances from land-cover heterogeneity, cloud contamination, and background reflectance. Despite these challenges, our data-driven analysis offers critical insights into the impact of VCC on LAI, allowing us to uncover broad trends and patterns across large spatial scales.<\/p>\n Our findings are based on the AVHRR GIMMS LAI3g dataset3<\/a>,49<\/a> (see Robustness Analysis), which provides a comprehensive and continuous time series of LAI during the period from 1982 to 2015. Generated based on a Feed-Forward Neural Network (FFNN) model, it has a spatial resolution of 1\/12\u00b0, is a bi-monthly composite, and has global coverage. It is one of the longest series among current LAI datasets and has demonstrated good performance in representing vegetation dynamics3<\/a>,10<\/a>,11<\/a>. We first interpolated the GIMMS LAI3g into 0.05\u00b0 grid spacing by using linear interpolation to fit the VCF5KYR land-cover data and averaged LAI into a monthly scale by taking the temporal mean of all valid records. We then retrieved the maximum monthly LAI that represents the peak vegetation growth35<\/a> for each pixel in each year. We calculated the annual mean LAI simultaneously to discuss the uncertainties associated with the proxy of vegetation greenness (Supplementary Fig.\u00a010<\/a>; Robustness Analysis). The quantitative attributions of long-term LAI changes were determined for the period from 1982 to 2015 to align with the availability of VCF5KYR and LAI products.<\/p>\n To provide a robust assessment associated with satellite-based products, we utilized LAI data from three alternative long-term satellite products to examine attributions of LAI change: NOAA CDR50<\/a>, which spans from July 1981 to the present with daily frequency and a 0.05\u00b0 resolution (Supplementary Fig.\u00a011<\/a>); GLASS51<\/a>, which offers 8-day LAI data with a spatial resolution of 0.05\u00b0 and covers the period from 1981 to 2018 (Supplementary Fig.\u00a012<\/a>); GLOBMAP52<\/a>, which offers biweekly LAI data during 1981\u20132000 and 8-day LAI data during 2001\u20132020 at 8\u2009km resolution (Supplementary Fig.\u00a013<\/a>). All LAI datasets were interpolated to a 0.05\u00b0 grid spacing and aggregated into monthly scales for analyses.<\/p>\n Disentanglement of the vegetation greenness<\/p>\n As shown in Supplementary Fig.\u00a01<\/a>, the main concept of our framework contains two parts. Developing such a framework is motivated by the lack of LAI observations for specific vegetation types, as existing LAI data are grid-based at a resolution of ~\u20095 kilometers, with each grid cell representing a single value that averages across multiple vegetation types. This means the lack of information on individual vegetation types for LAI, making it difficult to quantify how sub-grid scale vegetation cover changes contribute to overall LAI variations.<\/p>\n The first step is elucidating the proportion of different vegetation cover types from satellite-observed LAI based on an improved Bounded-Variable Least Squares (BVLS) model60<\/a> (Eqs.\u00a01<\/a>\u20133<\/a>). In our study, we extended the traditional BVLS approach by not only defining value ranges for the three greenness variables (\\({{{\\rm{\\beta }}}}_{{{\\rm{TC}}}}\\), \\({{{\\rm{\\beta }}}}_{{{\\rm{SV}}}}\\), \\({{{\\rm{\\beta }}}}_{{{\\rm{BG}}}}\\)) but also imposing additional constraints on their relative magnitudes (Eq.\u00a03<\/a> and Supplementary Fig.\u00a019<\/a>). These improvements ensured that the estimated parameters adhered to physically meaningful relationships, enhancing the stability and reliability of the results, especially under conditions of noise or multicollinearity. High-resolution satellite-observed LAI data (e.g., Supplementary Fig.\u00a02<\/a>) and the VCF5KYR product14<\/a> (Supplementary Fig.\u00a03<\/a>) during 1982\u20132015 are used as inputs for the statistical model. To overcome the spatial homogeneity issue raised by the coarse precision of VCF data, the percentage of TC, SV, and BG, as well as the LAI data, all at 0.05\u00b0\u2009\u00d7\u20090.05\u00b0 resolution, were first aggregated to 0.1\u00b0\u2009\u00d7\u20090.1\u00b0 cells for each year by taking the average of all valid records. These 0.1\u00b0 cells were then used as input to the BVLS regression. We used 3\u00b0\u2009\u00d7\u20093\u00b0 grid spacing as the analysis unit and retained only the grids and their inside cells with more than 20% valid LAI and vegetation cover observations at their original resolution. For detailed discussion regarding the size of cells and grids, please refer to Robustness Analysis.<\/p>\n We decomposed the observed LAI to TC, SV, and BG at the 3\u00b0\u2009\u00d7\u20093\u00b0 grid level for each year (Supplementary Fig.\u00a01<\/a>). The following regression was established for each valid 0.1\u00b0\u2009\u00d7\u20090.1\u00b0 cell to account for the quantitative relationship between LAI and vegetation cover in the current 3\u00b0\u2009\u00d7\u20093\u00b0 grid:<\/p>\n $$\\left\\{\\begin{array}{c}{{{\\rm{LAI}}}}_{1}={{{\\rm{\\beta }}}}_{{{\\rm{TC}}}}\\cdot {{{\\rm{TC}}}}_{1}+{{{\\rm{\\beta }}}}_{{{\\rm{SV}}}}\\cdot {{{\\rm{SV}}}}_{1}+{{{\\rm{\\beta }}}}_{{{\\rm{BG}}}}\\cdot {{{\\rm{BG}}}}_{1}+{{{\\rm{\\varepsilon }}}}_{1}\\\\ {{{\\rm{LAI}}}}_{2}={{{\\rm{\\beta }}}}_{{{\\rm{TC}}}}\\cdot {{{\\rm{TC}}}}_{2}+{{{\\rm{\\beta }}}}_{{{\\rm{SV}}}}\\cdot {{{\\rm{SV}}}}_{2}+{{{\\rm{\\beta }}}}_{{{\\rm{BG}}}}\\cdot {{{\\rm{BG}}}}_{2}+{{{\\rm{\\varepsilon }}}}_{2}\\\\ \\ldots \\\\ {{{\\rm{LAI}}}}_{{{\\rm{i}}}}={{{\\rm{\\beta }}}}_{{{\\rm{TC}}}}\\cdot {{{\\rm{TC}}}}_{{{\\rm{i}}}}+{{{\\rm{\\beta }}}}_{{{\\rm{SV}}}}\\cdot {{{\\rm{SV}}}}_{{{\\rm{i}}}}+{{{\\rm{\\beta }}}}_{{{\\rm{BG}}}}\\cdot {{{\\rm{BG}}}}_{{{\\rm{i}}}}+{{{\\rm{\\varepsilon }}}}_{{{\\rm{i}}}}\\\\ \\ldots \\\\ {{{\\rm{LAI}}}}_{{{\\rm{n}}}}={{{\\rm{\\beta }}}}_{{{\\rm{TC}}}}\\cdot {{{\\rm{TC}}}}_{{{\\rm{n}}}}+{{{\\rm{\\beta }}}}_{{{\\rm{SV}}}}\\cdot {{{\\rm{SV}}}}_{{{\\rm{n}}}}+{{{\\rm{\\beta }}}}_{{{\\rm{BG}}}}\\cdot {{{\\rm{BG}}}}_{{{\\rm{n}}}}+{{{\\rm{\\varepsilon }}}}_{{{\\rm{n}}}}\\end{array}\\right.$$<\/p>\n \n (1)\n <\/p>\n where \\({{\\rm{n}}}\\) is the number of valid cells with a maximum of 900; \\({{{\\rm{LAI}}}}_{{{\\rm{i}}}}\\) is the satellite-observed LAI value at cell \\({{\\rm{i}}}\\); \\({{{\\rm{TC}}}}_{{{\\rm{i}}}}\\), \\({{{\\rm{SV}}}}_{{{\\rm{i}}}}\\), and \\({{{\\rm{BG}}}}_{{{\\rm{i}}}}\\) are the percentage of TC, SV, and BG at cell \\({{\\rm{i}}}\\), respectively, and their sum equals 100%; \\({{{\\rm{\\beta }}}}_{{{\\rm{TC}}}}\\), \\({{{\\rm{\\beta }}}}_{{{\\rm{SV}}}}\\), and \\({{{\\rm{\\beta }}}}_{{{\\rm{BG}}}}\\) represent the estimated LAI parameter (in m2 m\u22122) of the current 3\u00b0\u2009\u00d7\u20093\u00b0 grid for assumed 100% coverage of TC, SV, and BG, respectively. Our improved BVLS model solves the linear least-squares regression with upper and lower bounds on the parameters (i.e., \\({{{\\rm{\\beta }}}}_{{{\\rm{TC}}}}\\), \\({{{\\rm{\\beta }}}}_{{{\\rm{SV}}}}\\), and \\({{{\\rm{\\beta }}}}_{{{\\rm{BG}}}}\\)), given as:<\/p>\n $${\\min }_{{{\\rm{l}}}\\le {{\\rm{x}}}\\le {{\\rm{u}}}}{{\\rm{||Ax}}}-{{\\rm{b}}}{||}_2$$<\/p>\n \n (2)\n <\/p>\n where \\({{\\rm{A}}}\\) is an \\({{\\rm{n}}}\\)-by-3 matrix of observed vegetation cover, \\({{\\rm{b}}}\\) is an \\({{\\rm{n}}}\\)-by-1 vector of observed LAI. Specifically, the parameters are constrained by the following inequality with upper and lower (non-negative) limits:<\/p>\n \\({{{\\rm{\\beta }}}}_{\\max }\\) was defined individually for each grid in each year based on the observation. This inequality was verified using a stratified random-sample approach (details refer to Supplementary Fig.\u00a019<\/a>). The significance of F-statistics (p-value) and the coefficient of determination (R2) for the BVLS model were calculated to show the goodness of fit in estimating vegetation greenness parameters. The F-test, with a significance level of p-value\u20092 reflects how well the model explains the pixel-level (0.1\u00b0) LAI variation using the grid-level averaged greenness parameters (\\({{\\rm{\\beta }}}\\)) within a 3\u00b0\u2009\u00d7\u20093\u00b0 grid. Note that the R2 value could be negative, indicating that the improved BVLS model with constraints is not appropriate for the samples.<\/p>\n The improved BVLS model was established and solved in each 3\u00b0\u2009\u00d7\u20093\u00b0 grid for each year. The models with low confidence were removed based on the model parameters (i.e., a p-value of F-statistics >\u20090.05 and\/or negative R2). Grids were retained only if they had at least 16 years of credible BVLS models during the study period. Overall, there are 1544 valid grids left after the model quality assessment, accounting for 96.0% of all grids. Among them, 1299 grids (80.8%) have valid results in all 32 years (the years 1994 and 2000 were excluded from analyses due to a lack of VCF data). The model performs better (i.e., exhibits higher R2 values) in regions with LAI increase hotspots and\/or significant vegetation cover change (e.g., East Asia, South Asia, Europe, Sahel, eastern South America, eastern Australia), while it shows lower R2 in areas with dense, uniform vegetation (e.g., the Amazon and Congo forests, and high-latitude tundra), likely due to limited variability in vegetation composition (Supplementary Fig.\u00a04<\/a>). The invalid grids are mainly located in coastal areas and tropical rainforests (Supplementary Fig.\u00a04<\/a>). Based on the improved BVLS statistical model and the high-resolution satellite-observed LAI and VCF data, we successfully disentangled the \\({{{\\rm{\\beta }}}}_{{{\\rm{TC}}}}\\), \\({{{\\rm{\\beta }}}}_{{{\\rm{SV}}}}\\), and \\({{{\\rm{\\beta }}}}_{{{\\rm{BG}}}}\\) parameters for each 3\u00b0\u2009\u00d7\u20093\u00b0 grid and each year (Supplementary Fig.\u00a05<\/a>).<\/p>\n Attribution of Earth\u2019s LAI response to vegetation cover change<\/p>\n In the second step of our framework, we reconstructed annual LAI values by using the parameters retrieved from the BVLS model. For each 3\u00b0\u2009\u00d7\u20093\u00b0 grid in each year, the LAI value (\\({{\\rm{\\sigma }}}{{\\rm{LAI}}}\\); in m2 m\u22122) was calculated as:<\/p>\n $${{\\rm{\\sigma }}}{{\\rm{LAI}}}={{{\\rm{\\beta }}}}_{{{\\rm{TC}}}}\\cdot {{\\rm{TC}}}+{{{\\rm{\\beta }}}}_{{{\\rm{SV}}}}\\cdot {{\\rm{SV}}}+{{{\\rm{\\beta }}}}_{{{\\rm{BG}}}}\\cdot {{\\rm{BG}}}$$<\/p>\n \n (4)\n <\/p>\n where \\({{\\rm{TC}}}\\), \\({{\\rm{SV}}}\\), and \\({{\\rm{BG}}}\\) are the mean proportion (%) of tree cover, short vegetation, and bare ground averaged from the valid cells in the current grid. For regional and global analyses, the grid-level estimations of LAI were aggregated by taking the area-weighted average approach, with the upper and lower ranges of the 95% confidence interval being estimated. The long-term trends of reconstructed LAI (refer to as \\({{\\rm{\\delta }}}{{\\rm{LAI}}}\\)) were derived by using linear regression during the study period at the grid, regional and global scales. The upper and lower ranges of \\({{\\rm{\\delta }}}{{\\rm{LAI}}}\\) were estimated based on the standard error of the mean. Generally, the reconstructed \\({{\\rm{\\delta }}}{{\\rm{LAI}}}\\) performs very well in reproducing the trend and variation of the observed LAI dynamics across the global, continental, and regional scales (details refer to Validation).<\/p>\n To attribute the reconstructed LAI trends to drivers related to VCC, we differentiated Eq. (4<\/a>) with respect to the year (denoted as \\({{\\rm{y}}}\\)), as follows:<\/p>\n $${{\\rm{\\delta }}}{{\\rm{LAI}}}=\t\\frac{\\partial {{\\rm{\\sigma }}}{{\\rm{LAI}}}}{\\partial {{\\rm{y}}}}=\\frac{\\partial {{{\\rm{\\beta }}}}_{{{\\rm{TC}}}}}{\\partial {{\\rm{y}}}}\\cdot \\overline{{{\\rm{TC}}}}+\\frac{\\partial {{\\rm{TC}}}}{\\partial {{\\rm{y}}}}\\cdot \\overline{{{{\\rm{\\beta }}}}_{{{\\rm{TC}}}}}+\\frac{\\partial {{{\\rm{\\beta }}}}_{{{\\rm{SV}}}}}{\\partial {{\\rm{y}}}}\\cdot \\overline{{{\\rm{SV}}}} \\\\ \t+ \\frac{\\partial {{\\rm{SV}}}}{\\partial {{\\rm{y}}}}\\cdot \\overline{{{{\\rm{\\beta }}}}_{{{\\rm{SV}}}}}+\\frac{\\partial {{{\\rm{\\beta }}}}_{{{\\rm{BG}}}}}{\\partial {{\\rm{y}}}}\\cdot \\overline{{{\\rm{BG}}}}+\\frac{\\partial {{\\rm{BG}}}}{\\partial {{\\rm{y}}}}\\cdot \\overline{{{{\\rm{\\beta }}}}_{{{\\rm{BG}}}}}$$<\/p>\n \n (5)\n <\/p>\n where the overbar denotes the multi-year average value of the variable. On the right-hand side of Eq. (5<\/a>), the six terms represent the changes in reconstructed LAI associated with the changes in the greenness of tree cover (\\({{{\\rm{\\beta }}}}_{{{\\rm{TC}}}}\\)), the proportion of tree cover (\\({{\\rm{TC}}}\\)), the greenness of short vegetation (\\({{{\\rm{\\beta }}}}_{{{\\rm{SV}}}}\\)), the proportion of short vegetation (\\({{\\rm{SV}}}\\)), the greenness of bare ground (\\({{{\\rm{\\beta }}}}_{{{\\rm{BG}}}}\\)) and the proportion of bare ground (\\({{\\rm{BG}}}\\)), respectively.<\/p>\n The reconstructed \\({{\\rm{\\delta }}}{{\\rm{LAI}}}\\) was further attributed to vegetation cover change (\\({{\\rm{\\delta }}}{{\\rm{VC}}}\\)), given by:<\/p>\n $${{\\rm{\\delta }}}{{\\rm{VC}}}=\\frac{\\partial {{\\rm{TC}}}}{\\partial {{\\rm{y}}}}\\cdot \\overline{{{{\\rm{\\beta }}}}_{{{\\rm{TC}}}}}+\\frac{\\partial {{\\rm{SV}}}}{\\partial {{\\rm{y}}}}\\cdot \\overline{{{{\\rm{\\beta }}}}_{{{\\rm{SV}}}}}+\\frac{\\partial {{\\rm{BG}}}}{\\partial {{\\rm{y}}}}\\cdot \\overline{{{{\\rm{\\beta }}}}_{{{\\rm{BG}}}}}$$<\/p>\n \n (6)\n <\/p>\n and the contributions of VCC were calculated as the ratio of \\({{\\rm{\\delta }}}{{\\rm{VC}}}\\) to \\({{\\rm{\\delta }}}{{\\rm{LAI}}}\\). Similarly, the attributions of the changes in the proportion of tree cover (\\(\\frac{\\partial {{\\rm{TC}}}}{\\partial {{\\rm{y}}}}\\cdot \\overline{{{{\\rm{\\beta }}}}_{{{\\rm{TC}}}}}\\)), short vegetation (\\(\\frac{\\partial {{\\rm{SV}}}}{\\partial {{\\rm{y}}}}\\cdot \\overline{{{{\\rm{\\beta }}}}_{{{\\rm{SV}}}}}\\)) and bare ground (\\(\\frac{\\partial {{\\rm{BG}}}}{\\partial {{\\rm{y}}}}\\cdot \\overline{{{{\\rm{\\beta }}}}_{{{\\rm{BG}}}}}\\)) were calculated at regional and global scales. The remaining part (\\({{\\rm{\\delta }}}{{\\rm{ID}}}\\)) represents the LAI changes in constant vegetation, which is largely influenced by indirect factors, as given by:<\/p>\n $${{\\rm{\\delta }}}{{\\rm{ID}}}=\\frac{\\partial {{{\\rm{\\beta }}}}_{{{\\rm{TC}}}}}{\\partial {{\\rm{y}}}}\\cdot \\overline{{{\\rm{TC}}}}+\\frac{\\partial {{{\\rm{\\beta }}}}_{{{\\rm{SV}}}}}{\\partial {{\\rm{y}}}}\\cdot \\overline{{{\\rm{SV}}}}+\\frac{\\partial {{{\\rm{\\beta }}}}_{{{\\rm{BG}}}}}{\\partial {{\\rm{y}}}}\\cdot \\overline{{{\\rm{BG}}}}$$<\/p>\n \n (7)\n <\/p>\n The interannual variations of the VCC-induced LAI (\\({{{\\rm{\\sigma }}}{{\\rm{VC}}}}_{{{\\rm{y}}}}\\); in m2 m\u22122) could be reconstructed using a variable-controlling approach61<\/a> by fixing the greenness of tree cover, short vegetation and bare ground as multi-year average values:<\/p>\n $${{{\\rm{\\sigma }}}{{\\rm{VC}}}}_{{{\\rm{y}}}}=\\overline{{{{\\rm{\\beta }}}}_{{{\\rm{TC}}}}}\\cdot {{{\\rm{TC}}}}_{{{\\rm{y}}}}+\\overline{{{{\\rm{\\beta }}}}_{{{\\rm{SV}}}}}\\cdot {{{\\rm{SV}}}}_{{{\\rm{y}}}}+\\overline{{{{\\rm{\\beta }}}}_{{{\\rm{BG}}}}}\\cdot {{{\\rm{BG}}}}_{{{\\rm{y}}}}$$<\/p>\n \n (8)\n <\/p>\n The anomalies of \\({{\\rm{\\sigma }}}{{\\rm{VC}}}\\) during the study period (e.g., Fig.\u00a03a<\/a>) were then obtained by subtracting the multi-year mean from the resulting time series. The variable-controlling approach was also applied spatially to calculate the contribution of LAI change within specific continents, countries, or regions to global \\({{\\rm{\\delta }}}{{\\rm{LAI}}}\\) (Fig.\u00a02c<\/a> and Supplementary Table\u00a02<\/a>).<\/p>\n Estimating potential effect of stable cropland and managed forest<\/p>\n The global map of cropland extent62<\/a> spanning from 2003 to 2019, with a spatial resolution of 0.00025\u00b0, was utilized to delineate the extent of stable cropland without VCC signals. The original data were aggregated into 3\u00b0 grid spacing using the area-weighted average approach, generating a global map showing the proportion of stable cropland (Supplementary Fig.\u00a018a<\/a>). Assuming that the LAI change of these areas is solely driven by human land-use practices, we estimated the potential contribution of cropland management (\\({{\\rm{\\delta }}}{{\\rm{SC}}}=\\frac{\\partial {{{\\rm{\\beta }}}}_{{{\\rm{SV}}}}}{\\partial {{\\rm{y}}}}\\cdot {{\\rm{SC}}}\\)) globally and regionally, where \\({{\\rm{SC}}}\\) represents the proportion of stable cropland at the grid level. Note that this assumption may lead to an overestimation of the contribution by disregarding the background LAI change from indirect factors like CO2 fertilization, climate change and nitrogen deposition.<\/p>\n Validation<\/p>\n Our statistical models were validated for their ability to reproduce the characteristics of \\({{\\rm{\\sigma }}}{{\\rm{LAI}}}\\) and \\({{\\rm{\\delta }}}{{\\rm{LAI}}}\\) through comparisons with satellite-derived observations. At the grid level, the reconstructed annual \\({{\\rm{\\sigma }}}{{\\rm{LAI}}}\\) shows high correlation coefficients (r\u2009>\u20090.99 for over 95.5% of valid 3\u00b0\u2009\u00d7\u20093\u00b0 grids; Supplementary Fig.\u00a06a<\/a>) and minimal bias (2\u2009m\u22122 for over 95.3% of grids; Supplementary Fig.\u00a06b<\/a>). Only slight overestimations occur in semi-arid regions, such as Central Asia and Sahel, while minor underestimations are mainly observed in tropical areas. For the reconstructed \\({{\\rm{\\delta }}}{{\\rm{LAI}}}\\), more than 90% of the grids exhibit a bias of less than \u00b1\u20090.01\u2009m2\u2009m\u22122 decade\u22121, with particularly low bias in grids with significant LAI changes (Supplementary Fig.\u00a07<\/a>). These results indicate that while the pixel-level fit of BVLS model may be moderate (R2 value of ~\u20090.5 in Supplementary Fig.\u00a04<\/a>), our framework still captures grid-level signals effectively in terms of magnitude and interannual variability, as evidenced by high correlation coefficients and low bias metrics. Note that our attribution analysis focus on LAI dynamics at the 3\u00b0\u2009\u00d7\u20093\u00b0 grid-level, and applying the estimated grid-average \\({{\\rm{\\beta }}}\\) values to disentangle and attribute LAI at finer scales (e.g., 0.1\u00b0\u2009\u00d7\u20090.1\u00b0 pixels) may introduce uncertainty.<\/p>\n Overall, the results show high-quality performance globally with a high correlation coefficient (r\u2009=\u20090.997, p\u20091a and Supplementary Table\u00a01<\/a>) between the reconstructed annual LAI and the observations. On a continental scale, our reconstruction is robust in capturing the temporal variation and trend of the observed LAI, with correlation coefficient values varying from 0.977 to 0.999; and relative bias ranged between \u2212\u20090.9% and 6.7% (Fig.\u00a01c<\/a> and Supplementary Table\u00a01<\/a>). In evaluating the model performance of the largest countries and regions, we show that our model had strong capability across most countries\/regions, with correlation coefficient values larger than 0.97 and relative biases of less than 6% (Supplementary Table\u00a01<\/a>). The model is relatively weak in reproducing the trend of observed LAI change in Argentina and the Democratic Republic of the Congo. This difficulty likely stems from the low coefficient of determination for the BVLS model that was induced by the spatial homogeneity of vegetation cover in those grids.<\/p>\n The robustness analyses further confirmed the quality of our statistical models by demonstrating similar results when using the annual mean LAI (r\u2009=\u20090.990, p\u200910), the NOAA CDR LAI product (r=\u20090.993, p\u200911), the GLASS LAI product (r\u2009=\u20090.968, p\u200912), and the GLOBMAP LAI product (r\u2009=\u20090.990, p\u200913). For the three alternative LAI products, the reconstructed \\({{\\rm{\\delta }}}{{\\rm{LAI}}}\\) also exhibits a low bias of less than \u00b1\u20090.01\u2009m2\u2009m\u22122 decade\u22121 in 94.9%, 94.4%, and 81.6% of the grids, with particularly low bias in grids exhibiting significant LAI changes (Supplementary Fig.\u00a020<\/a>).<\/p>\n Robustness analysisGrid spacing<\/p>\n Initially, we conducted our analysis using a grid spacing of 0.5\u00b0\u2009\u00d7\u20090.5\u00b0 as the analysis unit, with each grid containing up to 100 samples at the original 0.05\u00b0\u2009\u00d7\u20090.05\u00b0 resolution. However, during the application of the improved BVLS models to those samples, we encountered approximately 40% of grids that produced models with low confidence, as indicated by a p-value of F-statistics >\u20090.05 and\/or negative R2 (Supplementary Table\u00a03<\/a>). This can be attributed to the homogeneity of vegetation cover fraction caused by the relatively low precision (1%) of the VCF product14<\/a>. At a local scale (i.e., 0.5\u00b0\u2009\u00d7\u20090.5\u00b0 grid), many 0.05\u00b0\u2009\u00d7\u20090.05\u00b0 cells shared the same fraction of tree cover and\/or short vegetation due to the landscape similarities, yet exhibited different satellite-observed LAI values.<\/p>\n To tackle the issue of spatial homogeneity, we implemented a strategy involving spatial averaging of vegetation cover fractions and an increase in the number of samples at the grid level. We conducted experiments with various settings for grids and cells (Supplementary Table\u00a03<\/a>). Initially, the percentages of TC, SV, and BG, along with the LAI values at 0.05\u00b0\u2009\u00d7\u20090.05\u00b0 resolution, were aggregated to cells ranging from 0.1\u00b0\u2009\u00d7\u20090.1\u00b0 to 0.5\u00b0\u2009\u00d7\u20090.5\u00b0 for each year. This involved calculating the average of all valid records within each cell. Subsequently, the BVLS models were established and solved at the grid level, with grid size ranging from 1\u00b0\u2009\u00d7\u20091\u00b0 to 5\u00b0\u2009\u00d7\u20095\u00b0. Across different combinations of grid and cell sizes, the reconstructed contributions of \\(\\delta {VC}\\) consistently ranged from 18.0% to 23.5%, underscoring the non-negligible role of VCC in global LAI increasing (Supplementary Table\u00a03<\/a>). Importantly, the robustness of our findings was not affected by the choice of analysis unit or spatial average window. Finally, we present the results based on a 3\u00b0\u2009\u00d7\u20093\u00b0 grid spacing with a 0.1\u00b0\u2009\u00d7\u20090.1\u00b0 spatial average window, which yielded the highest ratio of valid grids (96.0%) and a low bias (1.3%) in the reconstruction of \\({{\\rm{\\delta }}}{{\\rm{LAI}}}\\).<\/p>\n Maximum monthly LAI and annual mean LAI<\/p>\n To assess the robustness of our proxy indicator of vegetation greenness, we investigated the uncertainties associated with different indices, comparing the results derived from annual mean LAI to those from maximum monthly LAI (Supplementary Fig.\u00a08<\/a>). The reconstructed global trend in annual mean LAI (0.030\u2009\u00b1\u20090.007\u2009m2\u2009m\u22122 decade\u22121; t(30)\u2009=\u20094.18, p\u2009p\u2009=\u20090.02) can be attributed to \\({{\\rm{\\delta }}}{{\\rm{VC}}}\\), confirming the substantial impact of VCC. Furthermore, the spatial distribution of \\({{\\rm{\\delta }}}{{\\rm{VC}}}\\) demonstrates a high level of consistency (Fig.\u00a01<\/a>\u2009d versus Supplementary Fig.\u00a08b<\/a>). We report the results based on the maximum monthly LAI due to the following considerations: (1) it represents the peak vegetation greenness during the local peak growing season, a period identified by the vegetation fields data14<\/a>; (2) it can better capture the interannual changes in the greenness of croplands and plantations with different cropping cycles, making it suitable for a global study35<\/a>; (3) it is less affected by missing data and can be easily retrieved without subjectively defining a growing season or dealing with multiple growing seasons; and (4) the results have higher statistical significances compared with those derived from annual mean LAI.<\/p>\n LAI products<\/p>\n We performed the analysis using NOAA CDR50<\/a>, GLASS51<\/a> and GLOBMAP52<\/a> as alternative products to the GIMMS LAI3g49<\/a> to test the uncertainties associated with different satellite-based datasets. Similarly, the reconstructed \\({{\\rm{\\delta }}}{{\\rm{LAI}}}\\) based on these three datasets captures the observed magnitude and temporal variation of global vegetation greenness satisfactorily (Supplementary Figs.\u00a011<\/a>\u201313<\/a>). Our findings consistently indicate that \\({{\\rm{\\delta }}}{{\\rm{VC}}}\\) accounts for 14.6\u2009\u00b1\u20095.5% (t(30)\u2009=\u20092.88, p\u2009=\u20090.007), 34.5\u2009\u00b1\u200910.4% (t(30)\u2009=\u20093.30, p\u2009=\u20090.003) and 28.6\u2009\u00b1\u200911.8% (t(30)\u2009=\u20092.44, p\u2009=\u20090.021) of \\(\\delta {LAI}\\) based on the NOAA CDR, the GLASS and the GLOBMAP products, respectively, reaffirming the validity of the importance of \\({{\\rm{\\delta }}}{{\\rm{VC}}}\\) in our main analysis. However, we observed an abrupt shift of LAI value around the year 2000 for these products (Supplementary Figs.\u00a011c<\/a>, 12c<\/a>, 13c<\/a>). The shift might be caused by limitations associated with sensor degradation and retrieval processes2<\/a>,56<\/a>,57<\/a>,58<\/a>,59<\/a>,63<\/a>. Given the potential impact of this issue on the estimation of long-term LAI change during 1982\u20132015 (e.g., the widespread LAI increasing signals across the globe in Supplementary Fig.\u00a011a<\/a>), we reported the assessment derived from the GIMMS LAI3g product while using the results from the NOAA CDR, GLASS, and GLOBMAP products as robustness tests.<\/p>\n Reporting summary<\/p>\n