Vegetation cover change as a growing driver of global leaf area index dynamics

Domain: earth_science
Published: 2025-10-20
Editorial outcome: major_revision

Abstract
The ongoing four-decade increase in global leaf area index (LAI), inferred from satellite observations, suggests enhanced carbon sequestration and evaporative cooling with potential benefits that may help mitigate climate warming. However, the role of vegetation cover change in driving these trends remains debated, raising concerns about the potential effects of reforestation and deforestation on ecosystem-climate interactions. Here, we develop a data-driven framework combining satellite-based observations to quantify the contributions of vegetation cover change to LAI dynamics across space and time. We find that vegetation cover change explains 18.1 ± 5.9% of the observed LAI increase since the 1980s, with pronounced contribution from increased tree cover primarily occurring in the northern hemisphere (33.8 ± 3.6%), partially offset by deforestation mainly taking place in the southern hemisphere (−15.7 ± 3.1%). Increases are most prominent in forestation regions such as China and Europe, while the sensitivity, defined as LAI change per unit percentage of tree cover gain, is lower than declines induced by deforestation in tropical areas. Our findings reveal that vegetation cover change is already, and increasingly, shaping global LAI dynamics to a greater extent than previously recognized, with important implications for future LAI trajectory projections, model development, and climate mitigation policies.

Introduction
A worldwide increase in leaf area index (LAI) has been reported based on satellite observations over the past four decades1–6. This emerging signal can partially alleviate global warming by promoting vegetation productivity and enhancing the land carbon sink7–9, as well as sustaining plant-driven evaporative cooling and accelerating the hydrological cycle10–12. Consequently, increased LAI assumes a crucial role in climate mitigation and ecosystem resilience, drawing considerable interest in global change science13.

Vegetation cover change (VCC), encompassing facets including reforestation, deforestation, and cropland expansion, has the potential to directly enhance or suppress LAI increases by altering vegetation extent and function2,14. While various indirect climatic and ecological factors4, including CO2 fertilization15,16, climate change (e.g., temperature, water availability, solar radiation)7,17, nitrogen deposition18, and agricultural management practices4,19, have been extensively studied as drivers of regional and global LAI increase, the specific contribution and impact of VCC on observed LAI changes remain contested2. As the world intensifies commitment to forest restoration and conservation as nature-based solutions to combating climate change20, unraveling the complexities of VCC-induced LAI increase becomes ever more critical. The urgency in doing so is underscored by reports revealing the persistent challenge of reducing deforestation despite previous efforts21. Enhanced understanding of the role of VCC in driving LAI changes can aid refining parameterization of vegetation processes in Earth system models, thereby improving the accuracy of projected vegetation dynamics and associated ecosystem benefits22,23.

Previous modeling studies have suggested that CO2 fertilization is the dominant contributor to the global rise in LAI3,15, while the effects of VCC were estimated to account for only 4% of the observed trend3. Conversely, studies utilizing field and satellite observations have emphasized a distinct overlap between croplands and observed LAI increases, particularly prominent in China and India4,19, leading to inferences that human agricultural management could potentially explain more than 33% of global LAI changes4. However, the influence of VCC on LAI change was blended among land management changes in their analyses. Much of the conflicting results and uncertainties may arise from the lack of comprehensive statistical exploration that can effectively separate the contributions of individual factors19,24, as well as limitations in models that inadequately represent fine-scale land-use practices and parameterize physical processes22,25. Specifically, large-scale forest restoration and tree planting26–30, along with commodity-driven deforestation21,31,32, have undoubtedly reshaped global vegetation cover14. Yet, their contribution to Earth’s LAI changes remains unclear, as only a few studies have investigated the influence of VCC over limited regions33,34. To date, there is a substantial lack of quantitative, data-driven attribution of LAI changes in response to VCC, raising concerns about the impact of VCC-related policies and actions for ecosystem monitoring and climate mitigation.

To tackle these challenges, we develop a framework based on statistical models that utilize satellite-observed LAI data and the vegetation continuous fields (VCF) product14 as inputs, covering the satellite era (“Methods”, Supplementary Figs. 1–5). The maximum monthly LAI35 is used because the VCF data identify the vegetation composition at the local peak growing season14. First, we decompose the greenness parameters (β) within 3°×3° grids in each year for three different land-cover categories—tree cover (βTC), short vegetation (βSV), and bare ground (βBG)—using LAI and VCF observations at the original resolution. The greenness parameters (β) represent a grid-averaged LAI value under the assumption of 100% coverage of each vegetation type within the grid cell. This is achieved using an improved Bounded-Variable Least Squares approach, which not only defines value ranges but also imposes constraints on the relative magnitudes of these three variables (Eqs. 1–3). Next, we reconstruct LAI (denoted as σLAI) across 3°×3° grid, regional and global scales based on the decomposed greenness parameters (β) and the proportion of each land-cover type, then calculate the trend of reconstructed LAI (δLAI) (Eqs. 4–6).

To attribute the drivers of LAI changes more precisely, we categorize them into two distinct groups: (1) LAI change signals related exclusively to VCC (δVC), and (2) LAI trends induced by indirect factors in areas where vegetation cover remained unchanged, calculated as the residual (δID = δLAI − δVC) (Supplementary Fig. 1). The statistical models are validated by comparing the reconstructed σLAI and δLAI with observations across various scales (Supplementary Figs. 6, 7 and Supplementary Table 1). In parallel, we demonstrate the robustness of our results by comparing against the annual mean LAI that is used to represent vegetation greenness, using other alternative long-term satellite products, and at variable grid spacings ranging from 0.5°×0.5° to 5°×5° (“Methods”). Importantly, our approach focuses on distinguishing LAI changes exclusively associated with VCC, providing solid, quantitative observational evidence of the LAI changes in response to VCC, rather than isolating specific natural or anthropogenic factors as in previous studies.

Results
Globally, there has been a significant increase in satellite-observed LAI maximum monthly LAI of 0.093 m2 m−2 decade−1 (standard error of the mean (SEM) = 0.011 m2 m−2 decade−1; statistical significance was assessed using a two-sided t test: t(32) = 8.12, p < 0.001) since the early 1980s. In agreement, our statistical model reconstructed a consistent trend of 0.094 ± 0.012 m2 m−2 decade−1 (mean ± SEM; t(30) = 8.04, p < 0.001), demonstrating its robustness in capturing the temporal variation and trend of the observed LAI changes (σLAI and δLAI in Fig. 1a; r = 0.997; Relative Bias = 1.3%). At the grid-level, the reconstructed annual σLAI shows strong correlations (r > 0.99) and minimal bias (< ±0.03 m²m⁻²) in over 95% of the grids (Supplementary Fig. 6). Overall, the reconstructed δLAI demonstrates high accuracy, particularly in areas with significant LAI changes (Supplementary Fig. 7).

The reconstructed global trend is further attributed to individual and combined drivers associated with VCC (Fig. 1b). VCC contributed 18.1 ± 5.9% of the global LAI increases (0.017 ± 0.005 m2 m−2 decade−1; t(30) = 3.09, p < 0.01), with the main contribution coming from increased tree cover (25.8%), which was partially offset by reduced presence of short vegetation (−7.4%) and bare ground (−0.3%). Despite the decline in global bare ground cover14, its contribution to global LAI changes remains minimal due to its notably lower LAI parameters (βBG) compared with tree cover and short vegetation (Supplementary Fig. 5). Focusing on regions where significant signals of LAI change were observed4,6, such as China, Europe, India, Sahel, and Brazil, we find that δVC can reach a mean of up to 0.030 ± 0.006 m2 m−2 decade−1 (t(30) = 5.16, p < 0.001). This accounts for 20.6 ± 4.0% of the LAI change in these areas (Supplementary Fig. 8).

At the continental level, Europe stands out with the highest magnitude of areal mean δLAI, where δVC dominates, accounting for 54.2% (Fig. 1c, d). Asia plays the most dominant role in global LAI trends, contributing 33.0% to overall changes, with δVC explaining 38.3% of regional variation and 12.5% of global shifts (Fig. 1d and Supplementary Table 2). VCCs in North America and Africa affect their regional LAI variations by 23.9% and 14.0%, respectively. In contrast, in South America and Oceania, the contributions of δVC are negative (−35.9% and −20.8%, respectively), indicating a decline in LAI due to forest loss. However, this signal is outweighed by the enhanced growth of existing vegetation at scales of tens to hundreds of kilometers (Fig. 1d versus Supplementary Fig. 9). Generally, changes of tree cover fraction are the primary factor driving δVC, except for Africa, where increases in short vegetation fraction play a more prominent role.

To ensure robustness of our findings, we tested our model under different settings (Methods). In addressing uncertainties related to the representation of vegetation greenness, we conducted additional tests using the annual mean LAI instead of the maximum monthly LAI (Supplementary Fig. 10). To assess uncertainties associated with different LAI datasets, we performed analyses using three alternative long-term satellite products (Supplementary Figs. 11–13). These tests reaffirmed the validity of the non-negligible role of δVC in our main analysis. In addition, the effect of varying spatial analysis scales showed that the proportions of δVC consistently fell within a similar range, reinforcing the reliability of our findings (Supplementary Table 3).

The magnitude of δVC exhibits a noticeable latitudinal gradient, with a contrast between the two hemispheres (Fig. 1e). The overall positive δVC globally results from LAI increases that mainly occur in the northern hemisphere (33.8 ± 3.6%; 0.031 ± 0.003 m2 m−2 decade−1, t(30) = 9.43, p < 0.001), compared with LAI decreases due to forest loss primarily in the southern hemisphere (−15.7 ± 3.1%; −0.014 ± 0.003 m2 m−2 decade−1, t(30) = −4.96, p < 0.001) (Fig. 2a). The temporal variations are generally consistent with major forestation or conservation projects, including the forest recovery following land abandonment in Eastern Europe and European Russia since the late 1980s (ref. 27), the Grain for Green Program and the Natural Forest Protection Program implemented around 2000 in China26,29, and the Bonn Challenge which started in 2011 (ref. 28). From 2003 to 2015, the VCC-induced LAI (σVC) reveals a substantial enhancement, consistent with the previously reported reversal of terrestrial biomass resulting from forest expansion mainly in Asia and Europe36.

Spatially, clusters of positive δVC are primarily found in recently afforested areas, spanning regions in China, Europe, the Middle Africa Climatic Transition Zone (MACTZ), Eastern Siberia, and the eastern United States, as well as in areas that experienced agricultural intensification, such as India and northern China (Fig. 1d and Supplementary Fig. 3b). The contribution of VCC could reach up to 29–73% in those areas (Fig. 2b). Notably, in China, δVC contributes to 72.5% of the nationwide δLAI, affirming the primary role of VCC in shaping long-term shifts in vegetation greenness of China33. This is likely in response to national programs such as the Great Green Wall of China, initiated in 1978 and is expected to continue until 2050, making it the largest afforestation initiative in human history29. India exhibits the highest magnitude of δVC, due to both the positive effects of agriculture expansion and tree cover gain, along with a substantial decrease in bare ground (Fig. 2b, c). The contributions of δVC from Russia, China, MATCZ, and the European Union to the global δLAI exceed those of India, and these four regions collectively account for 17.1% of Earth’s LAI increase (Fig. 2c and Supplementary Table 2).

In contrast, negative δVC predominates in Brazil, Argentina, Madagascar, and the Miombo woodlands in south-central Africa (Fig. 1d). Specifically, in Brazil and Argentina, δVC contributed to δLAI of about −50% (Fig. 2b). Some LAI increases in these areas largely resulted from indirect factors4 such as CO2 fertilization, climate change and nitrogen deposition3,7 (δID; Supplementary Fig. 9a). These indirect increases may have offset the deforestation-induced decrease in LAI at a regional scale over the past four decades (Fig. 1d versus Supplementary Fig. 9b). However, the persistence of these benefits amidst escalating disturbances cannot be guaranteed. If deforestation continues and warming intensify, the current buffering effect may weaken, leading to accelerated LAI decline and the potential disruption of tropical forest ecosystems, which are renowned for their biodiversity and high productivity37,38. Further assessment of the sensitivity of δVC to vegetation cover changes across continents reveals that the sensitivities to tree cover losses in South America (−0.306 ± 0.010 m2 m−2 decade−1 per 10% tree cover loss, t(151) = 31.24, p < 0.001) and Oceania (−0.360 ± 0.015, t(42) = 24.22, p < 0.001) are notably higher than those linked with tree cover gains in Asia (0.264 ± 0.008 m2 m−2 decade−1 per 10% tree cover gain, t(423) = 31.65, p < 0.001) and Europe (0.218 ± 0.009, t(166) = 24.73, p < 0.001) (Supplementary Fig. 14). This emphasizes the detrimental consequences of losing carbon-dense forests in tropical regions, meanwhile highlighting the urgent need for effective forest protection and management39,40.

The contribution of VCC to LAI trends has increased progressively over time (Fig. 3a). Analyzing the time series of δVC with expanding windows since 1982, we observe a significant upward trend at a rate of 0.017 ± 0.001 m2 m−2 decade−2 (t(16) = 12.36, p < 0.001) from the period 1982–1997 to 1982–2015. This increase is driven by a reduction in LAI decreases (LAI↓) from forest loss (0.012 ± 0.001 m2 m−2 decade−2; t(16) = 17.23, p < 0.001) and an amplification in LAI gains (LAI↑) from vegetation expansion (0.005 ± 0.001 m2 m−2 decade−2; t(16) = 5.33, p < 0.001). Concurrently, δVC holds an increasingly prominent role in Earth’s LAI changes, following an onlinear growth pattern over time. Globally, the fraction of total LAI change explained by δVC increased from around −10% before 2003 to a peak of 20.6% in 2009, followed by a slowdown in the 2010s (Fig. 3a). This fluctuation may reflect a stalling of LAI increases by indirect factors during the 2000s (ref. 24). At the continental level, the increased magnitude and contribution of δVC from both enhanced LAI↑ and reduced LAI↓ are also evident in Asia, Europe, North America, and South America (Supplementary Figs. 15, 16). Meanwhile, a declining trend in vegetation gain in Africa has led to a steady weakening of δVC. Overall, our findings highlight the increasing role of VCC in driving LAI changes, particularly prominent in the 2000s.

The spatially heterogeneous trends in δVC across countries and regions reflect diverse VCC processes and geographical constraints (Supplementary Fig. 17). In areas where increases in tree cover expansion dominated VCC, such as the Loess Plateau in China, Carpathian Mountains in Europe, and eastern United States, δVC increased substantially in the 2000s, reaching 0.10~0.20 m2 m−2 decade−1. In contrasts, water-limited savannas in Central Africa exhibited a gradual long-term decline in δVC during the past four decades (Fig. 3b, c), suggesting a potential future leveling-off of vegetation greenness under increasing water constraints in newly afforested areas like the Loess Plateau41,42. In northwest India, where short vegetation expansion dominates VCC, a slight decreasing trend of δVC exists (Fig. 3b, c). In regions that experienced deforestation, such as Cerrado and Gran Chaco in South America and Madagascar, the significant negative δVC exhibits a decreasing trend, reflecting a reduction in VCC-induced LAI↓ during 1982–2015 (Fig. 3d).

The influence of VCC, which reshaped the vegetation structure worldwide over the last four decades, has significant policy implications32,39,43. The non-negligible role of VCC-related actions in LAI increase is primarily driven by afforestation in the northern hemisphere26–30. However, on a global scale, this effect is partially offset by the LAI declines associated with widespread tropical deforestation and degradation14,21,31. Some regions, such as southeastern Amazonia, have transitioned from being carbon sinks to carbon sources due to forest loss and disturbance37,38. Deforestation in areas like the Cerrado and Gran Chaco resulted in negative δVC values, directly countering the benefits of indirect factors and leading to a shift from LAI increase to decrease (Supplementary Figs. 2 and 9).

Apart from forest change dynamics, cropland management, such as irrigation, fertilization, and species transformation, has also affected the LAI variability without altering the vegetation cover type. A recent study4 highlighted that human land-use practices might contribute to more than one-third of the observed Earth’s LAI increases, particularly in areas with extensive croplands, notably China (25%) and India (6.8%). Building on their conclusion that over one-third of global LAI increase is linked to direct land-use practices, our approach provides a more precise partitioning of LAI change signals. We estimate that cropland management without vegetation cover change (δSC) accounts for no more than 14.2 ± 2.0% (“Methods”; Supplementary Fig. 18a, b). This lower contribution compared to VCC is evident in China (15.6% of δSC versus 72.5% of δVC) and India (37.3% versus 43.5%), underscoring the necessity to distinguish the long-term LAI change signals induced by land-cover transition from those caused by land management alone (Supplementary Fig. 18c, d and Supplementary Table 3). Generally, agricultural intensification (management) predominated in developed regions, whereas agricultural-related LAI increase was primarily attributed to cropland expansion (VCC) in less developed regions44.

Concerns have persisted for decades regarding how climatic drivers affect global LAI dynamics2, motivated by negative feedback loops within the biosphere responding to climate change13. Recent findings point to a weakening or saturation trend in the observed LAI increase, along with elevated CO2 concentration and warming in this century24,45, as well as the adverse effects of extreme events and decreasing water availability on vegetation photosynthesis42,46. However, relying solely on satellite-observed LAI may obscure the true effects of those drivers, as VCC can partially mask the climatic signals of LAI change at the local-to-regional scale. Our results suggest that VCC may play an increasingly important role than previously thought, implying that the impacts of indirect factors on LAI increases may be overestimated. Hence, there is an urgent need to carefully exclude VCC impacts from the broader LAI changes and revisit the underlying mechanisms governing natural vegetation dynamics. Unfortunately, the ongoing intensification of ecosystem disturbances worldwide creates uncertainty in future LAI dynamics, with potential reductions driven not only by increasingly frequent extreme events, but unchecked deforestation, agriculture expansion, urbanization, and expanding energy infrastructure1,21,43,47.

Our findings contradict the commonly held belief in ecosystem models that the contribution of vegetation cover change to LAI dynamics (δVC) is minor1,2. This is evident in the notable difference between the model-projected δVC of 4% (ref. 3) and our higher observational attribution of 18.1 ± 5.9%. Ecosystem models likely underestimate the influence of VCC due to their omission or under-representation of the complex patterns of land-use practices and critical physiological processes2,4,22,25, which lead to an overemphasis on the benefits of indirect factors. While VCC’s influence on LAI dynamics is evident, with a majority of VCC in this study linked to direct human activities14, our aim is not to distinguish between anthropogenic and natural changes due to the coarse resolution of observations and the complexity of their interactions. Instead, our data-driven assessment highlights the role of VCC-related nature-based solutions in shaping LAI dynamics, both historically and in future projections. These solutions—such as forest preservation, ecological restoration, and afforestation40—support climate mitigation goals while offering co-benefits like terrestrial carbon uptake, biodiversity conservation, heat stress reduction for humans and livestock, improved ecosystem services, and greater resilience to future climate impacts26,39,48.

Methods
Datasets
Vegetation cover data. The Vegetation Continuous Fields version 1 product (VCF5KYR) was adopted to provide global fine-scale fractional vegetation cover information14. This product has a 0.05° × 0.05° 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 ≥ 5 m in height, abbreviated as TC), short vegetation (SV), and bare ground (BG), respectively for land pixels, with a precision of 1%, representing the vegetation composition of each pixel during the peak growing season14. 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. 3. Note that the years 1994 and 2000 were excluded from analyses due to a lack of VCF5KYR data (https://lpdaac.usgs.gov/products/vcf5kyrv001/). 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.

LAI data. 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,49, the National Oceanic and Atmospheric Administration (NOAA) Climate Data Record (CDR)50, the Global Land Surface Satellite (GLASS) AVHRR LAI V50 product51, and the GLOBMAP LAI V3 (ref 52.), 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,53, it may be subject to saturation effects in areas with dense vegetation, similar to other vegetation indices54,55. In addition, estimating LAI from satellite imagery comes with several well-known limitations2,56–59. 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.

Our findings are based on the AVHRR GIMMS LAI3g dataset3,49. We first interpolated the GIMMS LAI3g into 0.05° 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 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. 10; 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.

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, which spans from July 1981 to the present with daily frequency and a 0.05° resolution (Supplementary Fig. 11); GLASS51, which offers 8-day LAI data with a spatial resolution of 0.05° and covers the period from 1981 to 2018 (Supplementary Fig. 12); GLOBMAP52, which offers biweekly LAI data during 1981–2000 and 8-day LAI data during 2001–2020 at 8 km resolution (Supplementary Fig. 13). All LAI datasets were interpolated to a 0.05° grid spacing and aggregated into monthly scales for analyses.

Disentanglement of the vegetation greenness
As shown in Supplementary Fig. 1, 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 ~5 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.

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 (Eqs. 1–3). In our study, we extended the traditional BVLS approach by not only defining value ranges for the three greenness variables (βTC, βSV, βBG) but also imposing additional constraints on their relative magnitudes (Eq. 3 and Supplementary Fig. 19). 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. 2) and the VCF5KYR product14 (Supplementary Fig. 3) during 1982–2015 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° × 0.05° resolution, were first aggregated to 0.1° × 0.1° cells for each year by taking the average of all valid records. These 0.1° cells were then used as input to the BVLS regression. We used 3° × 3° 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.

We decomposed the observed LAI to TC, SV, and BG at the 3° × 3° grid level for each year (Supplementary Fig. 1). The following regression was established for each valid 0.1° × 0.1° cell to account for the quantitative relationship between LAI and vegetation cover in the current 3° × 3° grid:

LAI_i = β_TC · TC_i + β_SV · SV_i + β_BG · BG_i + ε_i (1)

where n is the number of valid cells with a maximum of 900; LAI_i is the satellite-observed LAI value at cell i; TC_i, SV_i, and BG_i are the percentage of TC, SV, and BG at cell i, respectively, and their sum equals 100%; β_TC, β_SV, and β_BG represent the estimated LAI parameter (in m2 m−2) of the current 3° × 3° 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., β_TC, β_SV, and β_BG), given as:

min ||Ax − b||_2, l ≤ x ≤ u (2)

where A is an n-by-3 matrix of observed vegetation cover, b is an n-by-1 vector of observed LAI. Specifically, the parameters are constrained by the following inequality with upper and lower (non-negative) limits:

0 < β_BG ≤ β_SV ≤ β_TC ≤ β_max (3)

β_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. 19). 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 < 0.05, is implemented to determine whether the model disentangles the LAI parameters better than random chance. While the R2 reflects how well the model explains the pixel-level (0.1°) LAI variation using the grid-level averaged greenness parameters (β) within a 3° × 3° grid. Note that the R2 value could be negative, indicating that the improved BVLS model with constraints is not appropriate for the samples.

The improved BVLS model was established and solved in each 3° × 3° grid for each year. The models with low confidence were removed based on the model parameters (i.e., a p-value of F-statistics > 0.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., exhibit 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. 4). The invalid grids are mainly located in coastal areas and tropical rainforests (Supplementary Fig. 4). Based on the improved BVLS statistical model and the high-resolution satellite-observed LAI and VCF data, we successfully disentangled the β_TC, β_SV, and β_BG parameters for each 3° × 3° grid and each year (Supplementary Fig. 5).

Attribution of Earth’s LAI response to vegetation cover change
In the second step of our framework, we reconstructed annual LAI values by using the parameters retrieved from the BVLS model. For each 3° × 3° grid in each year, the LAI value (σLAI; in m2 m−2) was calculated as:

σLAI = β_TC · TC + β_SV · SV + β_BG · BG (4)

where TC, SV, and 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 δLAI) were derived by using linear regression during the study period at the grid, regional and global scales. The upper and lower ranges of δLAI were estimated based on the standard error of the mean. Generally, the reconstructed δ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).

To attribute the reconstructed LAI trends to drivers related to VCC, we differentiated Eq. (4) with respect to the year (denoted as y), as follows:

δLAI = ∂σLAI / ∂y = (∂β_TC / ∂y) · TC + (∂TC / ∂y) · β_TC + (∂β_SV / ∂y) · SV + (∂SV / ∂y) · β_SV + (∂β_BG / ∂y) · BG + (∂BG / ∂y) · β_BG (5)

where the overbar denotes the multi-year average value of the variable. On the right-hand side of Eq. (5), the six terms represent the changes in reconstructed LAI associated with the changes in the greenness of tree cover (β_TC), the proportion of tree cover (TC), the greenness of short vegetation (β_SV), the proportion of short vegetation (SV), the greenness of bare ground (β_BG) and the proportion of bare ground (BG), respectively.

The reconstructed δLAI was further attributed to vegetation cover change (δVC), given by:

δVC = (∂TC / ∂y) · β_TC + (∂SV / ∂y) · β_SV + (∂BG / ∂y) · β_BG (6)

and the contributions of VCC were calculated as the ratio of δVC to δLAI. Similarly, the attributions of the changes in the proportion of tree cover (∂TC · β_TC), short vegetation (∂SV · β_SV) and bare ground (∂BG · β_BG) were calculated at regional and global scales. The remaining part (δID) represent the LAI changes in constant vegetation, which is largely influenced by indirect factors, as given by:

δID = (∂β_TC / ∂y) · TC + (∂β_SV / ∂y) · SV + (∂β_BG / ∂y) · BG (7)

The interannual variations of the VCC-induced LAI (σVC_y; in m2 m−2) could be reconstructed using a variable-controlling approach61 by fixing the greenness of tree cover, short vegetation and bare ground as multi-year average values:

σVC_y = β_TC · TC_y + β_SV · SV_y + β_BG · BG_y (8)

The anomalies of σVC during the study period (e.g., Fig. 3a) 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 δLAI (Fig. 2c and Supplementary Table 2).

Estimating potential effect of stable cropland and managed forest
The global map of cropland extent62 spanning from 2003 to 2019, with a spatial resolution of 0.00025°, was utilized to delineate the extent of stable cropland without VCC signals. The original data were aggregated into 3° grid spacing using the area-weighted average approach, generating a global map showing the proportion of stable cropland (Supplementary Fig. 18a). Assuming that the LAI change of these areas is solely driven by human land-use practices, we estimated the potential contribution of cropland management (δSC = ∂SC / ∂y · β_SV) globally and regionally, where 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.

Validation
Our statistical models were validated for their ability to reproduce the characteristics of σLAI and δLAI through comparisons with satellite-derived observations. At the grid level, the reconstructed annual σLAI shows high correlation coefficients (r > 0.99 for over 95.5% of valid 3° × 3° grids; Supplementary Fig. 6a) and minimal bias (< ±0.03 m2 m−2 for over 95.3% of grids; Supplementary Fig. 6b). 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 δLAI, more than 90% of the grids exhibit a bias of less than ±0.01 m2 m−2 decade−1, with particularly low bias in grids with significant LAI changes (Supplementary Fig. 7). These results indicate that while the pixel-level fit of BVLS model may be moderate (R2 value of ~0.5 in Supplementary Fig. 4), 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° × 3° grid-level, and applying the estimated grid-average β values to disentangle and attribute LAI at finer scales (e.g., 0.1° × 0.1° pixels) may introduce uncertainty.

Overall, the results show high-quality performance globally with a high correlation coefficient (r = 0.997, p < 0.001) and low relative bias (1.3%) (Fig. 1a and Supplementary Table 1) 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 −0.9% and 6.7% (Fig. 1c and Supplementary Table 1). 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 1). 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.

The robustness analyses further confirmed the quality of our statistical models by demonstrating similar results when using the annual mean LAI (r = 0.990, p < 0.001; Supplementary Fig. 10), the NOAA CDR LAI product (r = 0.993, p < 0.001; Supplementary Fig. 11), the GLASS LAI product (r = 0.968, p < 0.001; Supplementary Fig. 12), and the GLOBMAP LAI product (r = 0.990, p < 0.001; Supplementary Fig. 13). For the three alternative LAI products, the reconstructed δLAI also exhibits a low bias of less than ±0.01 m2 m−2 decade−1 in 94.9%, 94.4%, and 81.6% of the grids, with particularly low bias in grids exhibiting significant LAI changes (Supplementary Fig. 20).

Robustness analysis
Grid spacing. Initially, we conducted our analysis using a grid spacing of 0.5° × 0.5° as the analysis unit, with each grid containing up to 100 samples at the original 0.05° × 0.05° 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 > 0.05 and/or negative R2 (Supplementary Table 3). This can be attributed to the homogeneity of vegetation cover fraction caused by the relatively low precision (1%) of the VCF product14. At a local scale (i.e., 0.5° × 0.5° grid), many 0.05° × 0.05° cells shared the same fraction of tree cover and/or short vegetation due to the landscape similarities, yet exhibited different satellite-observed LAI values.

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 3). Initially, the percentages of TC, SV, and BG, along with the LAI values at 0.05° × 0.05° resolution, were aggregated to cells ranging from 0.1° × 0.1° to 0.5° × 0.5° 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° × 1° to 5° × 5°. Across different combinations of grid and cell sizes, the reconstructed contributions of δVC consistently ranged from 18.0% to 23.5%, underscoring the non-negligible role of VCC in global LAI increasing (Supplementary Table 3). 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° × 3° grid spacing with a 0.1° × 0.1° spatial average window, which yielded the highest ratio of valid grids (96.0%) and a low bias (1.3%) in the reconstruction of δLAI.

Maximum monthly LAI and annual mean LAI. 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. 8). The reconstructed global trend in annual mean LAI (0.030 ± 0.007 m2 m−2 decade−1; t(30) = 4.18, p < 0.001) agrees well with the observation; and 97.3% of all grids are retained after the quality control. Of the observed global LAI increase, 29.7% ± 12.7% (t(30) = 2.37, p = 0.02) can be attributed to δVC, confirming the substantial impact of VCC. Furthermore, the spatial distribution of δVC demonstrates a high level of consistency (Fig. 1d versus Supplementary Fig. 8b). 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; (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; (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.

LAI products. We performed the analysis using NOAA CDR50, GLASS51 and GLOBMAP52 as alternative products to the GIMMS LAI3g49 to test the uncertainties associated with different satellite-based datasets. Similarly, the reconstructed δLAI based on these three datasets captures the observed magnitude and temporal variation of global vegetation greenness satisfactorily (Supplementary Figs. 11–13). Our findings consistently indicate that δVC accounts for 14.6 ± 5.5% (t(30) = 2.88, p = 0.007), 34.5 ± 10.4% (t(30) = 3.30, p = 0.003) and 28.6 ± 11.8% (t(30) = 2.44, p = 0.021) of δLAI based on the NOAA CDR, the GLASS and the GLOBMAP products, respectively, reaffirming the validity of the importance of δVC in our main analysis. However, we observed an abrupt shift of LAI value around the year 2000 for these products (Supplementary Figs. 11c, 12c, 13c). The shift might be caused by limitations associated with sensor degradation and retrieval processes2,56–59,63. Given the potential impact of this issue on the estimation of long-term LAI change during 1982–2015 (e.g., the widespread LAI increasing signals across the globe in Supplementary Fig. 11a), 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.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This paper identifies an important knowledge gap in our understanding of how trends of changes in vegetation density/cover (what the authors refer to as ‘greening’) are associated with changes in land cover (as opposed to being driven directly by increased [CO2] and/or shifts in climate). The authors then devise a well thought through and very comprehensive approach to testing this questions using long-term satellite observations at a global scale. Overall I found the paper to be well written, easy to follow (despite the complex methodology) and statistically rigorous.\nI do however have a couple of concerns about the analysis, the conclusions that can be drawn from it and the way the authors define greening and browning.\nFirstly, the definitions of land cover used by the authors are really broad (forest, short vegetation, bare ground) and this to me really limits (i) the ability to identify the drivers behind the observed changes (e.g., human induced land use change vs climate-induced changes) and (ii) disentangle these effects from the very processes the authors aim to tease apart (e.g., how do you know that some of the changes in land cover aren’t driven by those in [CO2] and if they are can we really attribute the ‘greening’ effects to land cover change as opposed to the primary effects of changes in [CO2]?). Part of the reason for this the resolution of the analysis is quite coarse (~5km) which makes it really hard to be any more detailed in the attribution and qualification of land cover change. I know this is a fundamental limit of the data and not something the authors can do much about, but it does limit the ability of the authors to test their overarching question.\nSecond, the paper focuses heavily on the contrast between vegetation greening and browning. But what is actually being estimated are changes in LAI. This might sound like semantics, but I do think these are actually very different things (e.g., a loss of LAI in a high LAI system doesn’t mean vegetation browning in the way that most people think of this, which is vegetation die-off after something like a wildfire or drought). I would suggest reframing the paper around LAI change rather than greening vs browning. Linked to this, however, I do think the authors need to better acknowledge the well-know limitations of estimating LAI from satellite imagery. To convince readers that the patterns you are detecting really do translate to changes in LAI could you compare some of your current-day estimates of LAI with those derived from GEDI, which should be much closer to the truth? How well do the two compare, especially in high LAI systems?\nOther comments\nFig. 1: I’m always partial to a schematic diagram, but I really struggle to see what this adds here and I found the mix of drawing and equations quite challenging to navigate. Suggest removing from the main text. The other figures are great!\nL71-4: are these really a result of statistical limitations that don’t account for non-linear interactions? Is it not just that no one has really tried to address this question yet? If you stick with this statement I would suggest better justifying it.\nL104: what is this increase in tree cover replacing?\nL105: isn’t this minimal change in bare ground that you find in contrast to the results of recent global studies showing that bare ground cover has decreased globally, likely due to CO2 fertilization driving greater water use efficiency in arid environments?", "reviewer_2": "Aiming to evaluate the role of land cover change (LCC) in contributing to earth’s greening phenomenon observed by satellite remotely sensed leaf area index (LAI) increasing, the authors decomposed the trend of LAI into LCC-based signals including changes in tree cover, short vegetation, and bare ground. Then, the authors claimed that 18% of greening since 1980s can be attributed to LCC. The topic is interesting, but I think the current version lacks many details on the analysis and the discussion. My major concerns are:\nFirst, I may not agree that the consistency of reconstructed trend of σLAI and (cid:1)LAI can be used as evidence that the model is well representative (L99-101). Fig 2a shows mean global average, not grid-specific model performance. Assume a naïve model that only predicts the mean LAI for each grid, the same result of consistency can also be achieved. Particularly, R^2=0.997 is suspectable. Plus, SI Fig. 3 shows that average R^2 is about 0.45, so the model on average explains 45% of LAI variance, which suggests that other major factors are lacking in the model.\nSecond, related to the above problem, without controlling for other major factors you mentioned in L46-48, it may not be sound to claim that LCC alone contributes 18% to LAI greening, especially since this (cid:1)LAI was estimated by the model not observation and you are not comparing the magnitude of LCC with other factors.\nThird, it seems confusing to me that how this study solves the problems described in L60-76. I would appreciate more details and discussions on previous studies and how and why this analysis is better. Otherwise, it would just be another method derived result, not necessarily a better one, and the consistency/inconsistency with previous results would be meaningless. This is where I think is very important to clarify before talking about any results.\nSome specific comments:\n1. L79. “high-resolution” turns out to be 0.05° (~5 km) and then aggregated to 3° (~300 km)? Although AVHRR is named “high-resolution”, that was due to some historical issue. I’d recommend be more specific here.\n2. L83. Can you describe some details here about what the improved Bounded-Variable Least Squares approach is and how it solves the challenges you mentioned in the previous paragraph?\n3. L90-91. This explanation is a bit abrupt because you did not say max LAI was used before this sentence.\n4. L162-166. SI Fig. 1 does not talk about these factors, e.g., CO2, climate change, right? How did you get this conclusion?\n5. L227-228. Can you be more specific? Why do you think previous efforts underestimated LCC? And why do you think your method is better?\n6. Fig. 1. “soving” should be “solving”.\n7. L330-332. I don’t think I understand the “substantial lack of LAI observations”, can you give more descriptions?\n8. L370-371. SI Fig. 4 does not support “we successfully disentangled…” as they are basically model estimates, which do not say whether the model is good or bad.", "reviewer_3": "Reports of global greening are widespread and credible, but the contribution of land cover change has been uncertain, and create substantial confusion on the projected efficacy of natural climate solutions. By using vegetation continuous fields data to statistically disentangle land cover change from change within a land cover class (e.g. increases in forest productivity that may be attributable to CO2 fertilization), Zeng et al. provide credible data on how, when, and where land cover changes have influenced global greening.\nThe clear methods, robust analysis of multiple datasets at varying resolutions, disaggregated reporting of quantitative results (including in helpful detailed supplemental information), and the effective graphics and writing are major strengths of the manuscript. The analysis is convincing and comprehensive in its presentation and quality. These strengths are important because, as the authors rightly detail in the concluding section, the results of this analysis are vital for guiding natural climate solutions.\nIn my limited suggestions to improve the paper, I only found the few minor grammatical or stylistic errors below.\nLine 132: should be \"Contrasting effects...\"\nLine 226-227: Tone down to read: \"This underestimation likely originates...\"\nLine 231: I assume the 60% statistic is from the cited study, but it is unclear why that statistic from another study directly applies to this study. Perhaps avoid the confusion of that precise statistic and tone down to read: \"A majority of the LCC in this study is linked to direct human activities...\"\nLine 232: This sentence is awkward to read. Starting with \"Besides,\" is perhaps the problem, but the whole sentence is oddly phrased. Reword.\nLine 238: Change to the plural \"ecosystem services\"\nLine 271: \"recent failure\" is singular and time-specific, but deforestation is pervasive and chronic. Tone down to \"failures to limit deforestation\"."}, "author_response": "Response to the reviewers (NCOMMS-24-44410-T)\nReviewer #1 (Remarks to the Author):\nReviewer #1 General Comments: This paper identifies an important knowledge gap in our\nunderstanding of how trends of changes in vegetation density/cover (what the authors refer to\nas ‘greening’) are associated with changes in land cover (as opposed to being driven directly\nby increased [CO2] and/or shifts in climate). The authors then devise a well thought through\nand very comprehensive approach to testing this question using long-term satellite observations\nat a global scale. Overall, I found the paper to be well written, easy to follow (despite the\ncomplex methodology) and statistically rigorous.\nI do however have a couple of concerns about the analysis, the conclusions that can be drawn\nfrom it and the way the authors define greening and browning.\n[Response] We deeply thank you for your thoughtful and comprehensive review of our\nmanuscript. We greatly appreciate your recognition of the importance of our study in addressing\nthe existing knowledge gap, as well as your positive feedback on our methodology, statistical\nrigor, and the clarity of our writing. In response to your suggestions, we have carefully revised\nthe manuscript, providing additional descriptions and discussions, and have reframed the terms\nof greening and browning throughout the text. Detailed responses to each of your comments,\nalong with the corresponding revisions, are provided below. We believe that your insightful\nsuggestions have greatly enhanced the robustness and clarity of our study.\nSpecific comments:\n[Reviewer #1 Comment 1]: Firstly, the definitions of land cover used by the authors are really\nbroad (forest, short vegetation, bare ground) and this to me really limits (i) the ability to identify\nthe drivers behind the observed changes (e.g., human induced land use change vs climate-\ninduced changes) and (ii) disentangle these effects from the very processes the authors aim to\ntease apart (e.g., how do you know that some of the changes in land cover aren’t driven by those\nin [CO ] and if they are can we really attribute the ‘greening’ effects to land cover change as\n2\nopposed to the primary effects of changes in [CO ]?). Part of the reason for this the resolution\n2\nof the analysis is quite coarse (~5km) which makes it really hard to be any more detailed in the\n1\nattribution and qualification of land cover change. I know this is a fundamental limit of the data\nand not something the authors can do much about, but it does limit the ability of the authors to\ntest their overarching question.\n[Response] Thank you for your thoughtful feedback. We appreciate your valuable insights, and\nbased on your comments, we have made efforts to clarify the core concepts of our statistical\nframework and the main objectives of our analysis. Below, we provide a detailed explanation\nand corresponding revisions.\nOur research distinguishes itself from previous studies, such as those by Zhu et al. (2016) and\nChen et al. (2019), by offering quantitative evidence on how, when, and where vegetation cover\nchanges have influenced global LAI trends—an issue that has not been thoroughly explored in\nprior work. In our statistical framework, we have aimed to classify the drivers of LAI changes\nmore precisely by creating two distinct categories: (1) LAI change signals related exclusively\nto vegetation cover change (VCC)—a direct and observable consequence of global change\n(𝜹𝑽𝑪), and (2) LAI trends induced by indirect factors in areas where land-cover remains\nunchanged, calculated as the residual (𝜹𝑰𝑫 = 𝜹𝑳𝑨𝑰−𝜹𝑽𝑪). We fully agree with you that\ndisentangling the effects of natural versus anthropogenic factors on LAI and/or VCC dynamics\nis challenging, especially given the coarse resolution of observations and the complex\ninteractions involved. Therefore, rather than dwelling on the difficulties of determining whether\nVCC is driven by natural disturbances or human activities, we concentrate on distinguishing\nLAI changes directly associated with VCC using reliable observations within current limitations.\nMeanwhile, we emphasize the potential of human VCC-related actions in shaping LAI\ndynamics and highlight the promising effects of forestation programs on the maintenance of\necosystem services.\nWe have elaborated on this distinction in the Introduction and Implications sections of the\nrevised manuscript, emphasizing that our findings provide the first quantitative, observational\nevidence of the LAI changes response to LCC and offer more detailed insights on this topic\n(Page 4, Lines 87-107; Page 10, Lines 275-284 in the clean version of revised manuscript).\nas the coarse spatial resolution of the data in our analysis, constrain the ability to precisely\nattribute LAI changes to various individual factors at a finer scale. The broad land cover\nclassification used in our study was adopted in order to maximize the utility of the available\ndata while ensuring the statistical significance and robustness of the results. Despite the inherent\nresolution limitation, this approach still provides valuable information on VCC at a scale that\nis consistent with the LAI data. We believe that, given these constraints, the current framework\nprovides a meaningful contribution to understanding the influence of VCC on LAI changes,\nparticularly at regional-to-global scales. We have discussed this distinction in detail in the\nMethods section (Page 11, Lines 298-302).\nAdditionally, we have included a revised version of our conceptual framework (Fig. R1, also\nSupplementary Fig. 1 in the revised manuscript) to clarify our approach and better convey the\nmain objectives of our analysis. Once again, we sincerely thank you for your constructive\nfeedback, which has significantly helped to improve the clarity and focus of our manuscript.\nFig. R1. The conceptual framework for decoding the contribution of vegetation cover\n3\nchange (𝜹𝑽𝑪 ) on the satellite-observed LAI changes, along with the three individual\ndrivers. Two main steps are involved in the framework: 1) decompose the greenness parameters\nof three different land-cover types, i.e., tree cover (TC), short vegetation (SV), and bare ground\n(BG), from satellite-observed LAI based on an improved Bounded-Variable Least Squares\n(BVLS) model, and 2) reconstruct annual LAI values using the parameters and attribute the\nlong-term changes in reconstructed LAI to vegetation cover change.\n[Reviewer #1 Comment 2]: Second, the paper focuses heavily on the contrast between\nvegetation greening and browning. But what is actually being estimated are changes in LAI.\nThis might sound like semantics, but I do think these are actually very different things (e.g., a\nloss of LAI in a high LAI system doesn’t mean vegetation browning in the way that most people\nthink of this, which is vegetation die-off after something like a wildfire or drought). I would\nsuggest reframing the paper around LAI change rather than greening vs browning. Linked to\nthis, however, I do think the authors need to better acknowledge the well-known limitations of\nestimating LAI from satellite imagery. To convince readers that the patterns you are detecting\nreally do translate to changes in LAI could you compare some of your current-day estimates of\nLAI with those derived from GEDI, which should be much closer to the truth? How well do the\ntwo compares, especially in high LAI systems?\n[Response] We sincerely appreciate your feedback, which has greatly enhanced the clarity and\naccuracy of our work. While our use of the terms greening and browning aligns with numerous\nprevious studies (e.g., Zhu et al., 2016; Forzieri et al., 2017; Chen et al., 2019; Piao et al., 2020;\nCortés et al., 2021), we agree with you that LAI changes do not always equate to the traditional\nunderstanding of greening and browning. In response, we have revised the manuscript to focus\nmore explicitly on LAI changes, updating the terminology throughout the text, figures, and\nsupplementary materials for greater precision.\nWe have also included an acknowledgement of the limitations of in our Datasets section: “Note\nthat while LAI serves as a valuable proxy of vegetation greenness11,53, it may be subject to\nsaturation effects in areas with dense vegetation, similar to other vegetation indices54,55. In\naddition, estimating LAI from satellite imagery comes with several well-known limitations2,56-\n4\n59. These include spatial resolution constraints, saturation in densely vegetated areas, orbital\ndrift, sensor replacements, retrieval algorithm uncertainties, and disturbances from land-cover\nheterogeneity, cloud contamination, and background reflectance. Despite these challenges, our\ndata-driven analysis offers critical insights into the impact of VCC on LAI, allowing us to\nuncover broad trends and patterns across large spatial scales.” (Pages 11-12, Lines 309-317).\nTo ensure the robustness of our key findings, we employed three alternative long-term LAI\ndatasets to examine attributions of LAI change (Supplementary Figs. 11-13; Page 12, Lines\n333-340).\nFurthermore, we compared the LAI products used in this study with data from the Global\nEcosystem Dynamics Investigation (GEDI), which offers high-resolution lidar-based\nobservations of Earth’s three-dimensional structure starting from April 18, 2019. Due to the\ntemporal limitations of the GIMMS LAI3g dataset, we selected two alternative LAI datasets—\nGLASS and GLOBMAP at a spatial resolution of 0.05°—for the comparison. We compared the\nLAI data from these two products (for the period 2019–2020) with GEDI-derived metrics,\nspecifically: 1) the mean canopy height (L3 level, in meters; Dubayah et al., 2021), 2) the mean\naboveground biomass density (L4B level, in Mg ha-1; Dubayah et al., 2022), and 3) the total\nplant area index (L2B level, in m2 m-2; Burns et al., 2024). Both variables are provided at a 1\nkm × 1 km resolution, covering a latitude range from 51.6°S to 51.6°N.\nDespite discrepancies between the physical meanings of the LAI variable and the GEDI\nvegetation structure metrics (plant area index, canopy height or aboveground biomass density),\nwe found that the spatial patterns derived from GEDI vegetation structure metrics and satellite-\nbased LAI products are generally consistent (Fig. R2). This is particularly true in regions with\nhigh LAI, such as the Amazon, tropical Africa, Southern China, and Indonesia (Figs. R2a-c).\nAt the pixel level, the comparison between LAI data and GEDI-derived metrics shows strong\ncorrelation coefficients (r=0.71~0.82) (Fig. R3), although some differences were observed,\nespecially in areas of high LAI for the GLASS data.\n5\nLAI data. a, GEDI Level 3 gridded canopy height estimates at a 1 km × 1 km resolution (source:\nGEDI03_rh100_mean_2019108_2021104_002_02.tif). b, GEDI Level 4B gridded mean\naboveground biomass density estimates at a 1 km × 1 km resolution (source:\nGEDI04_B_MW019MW223_02_002_02_R01000M_MU.tif). c, GEDI Level 2B gridded total\nplant area index during 2019-2020 at a 1 km × 1 km resolution (Burns et al., 2024). d, Annual\nmean maximum monthly LAI at 0.05° resolution derived from the GLASS LAI data for 2019-\n2020. e, Annual mean maximum monthly LAI at 0.05° resolution derived from the GLOBMAP\nLAI data for 2019-2020.\nsatellite-based LAI data. The density is the number of valid pixel pairs between datasets within\ndifferent intervals. The grey lines represent the best fit using ordinary least squares regression.\nDue to the time-period limitations, we are currently unable to compare the long-term LAI\nchanges calculated in this study. However, we believe the comparison with GEDI, which is\nconsidered closer to the “truth” due to its lidar-based observations, strengthens the confidence\nthat satellite-derived LAI estimates remain a valuable variable for detecting vegetation structure\n(Fang et al., 2019). Finally, we sincerely appreciate your helpful suggestions, and we have\nreframed the manuscript to focus more explicitly on LAI changes to avoid any potential\nsemantic misunderstanding of greening or browning.\n[Reviewer #1 Comment 3]: Fig. 1: I’m always partial to a schematic diagram, but I really\nstruggle to see what this adds here and I found the mix of drawing and equations quite\nchallenging to navigate. Suggest removing from the main text. The other figures are great!\n[Response] Thank you for highlighting this issue. We have relocated the original Fig. 1 to the\nsupplementary information file as Supplementary Fig. 1. Additionally, we have revised this\nsupplementary figure to offer a more detailed breakdown of our motivation and methodology,\nthereby enabling readers to follow our approach more closely (Fig. R1 in the response letter).\n[Reviewer #1 Comment 4]: L71-4: are these really a result of statistical limitations that don’t\naccount for non-linear interactions? Is it not just that no one has really tried to address this\nquestion yet? If you stick with this statement, I would suggest better justifying it.\n[Response] We sincerely appreciate your insightful comment regarding the potential sources of\nconflicting results here. It is indeed that the lack of extensive exploration into this question has\ncontributed to the current gaps in understanding. We have revised the text to more accurately\nreflect this nuance as follows: “Much of the conflicting results and uncertainties may arise from\nthe lack of comprehensive statistical exploration that can effectively separate the contributions\nof individual factors19,24, as well as limitations in models that inadequately represent fine-scale\nland-use practices and parameterize physical processes22,25.” (Page 3, Lines 72-76). “So far,\n8\nthere is a substantial lack of quantitative, data-driven attributions of LAI changes in response\nto VCC, raising concerns about the impact of VCC-related policy and actions on LAI dynamics.”\n(Page 3, Lines 79-81).\n[Reviewer #1 Comment 5]: L104: what is this increase in tree cover replacing?\n[Response] The global VCF product used in our analysis classifies land pixels into three\ncategories: tree cover (TC), short vegetation (SV), and bare ground (BG), based on the\npercentage of vegetation composition, with their sums totaling 100% (Song et al., 2018).\nTherefore, the increase in tree cover is attributed to the decrease in short vegetation and bare\nground. We have involved the above information into the revision (Page 11, Lines 291-295).\n[Reviewer #1 Comment 6]: L105: isn’t this minimal change in bare ground that you find in\ncontrast to the results of recent global studies showing that bare ground cover has decreased\nglobally, likely due to CO fertilization driving greater water use efficiency in arid environments?\n2\n[Response] The vegetation continuous fields (VCF) dataset from Song et al. (2018), which we\nused in our study, also indicates a global decrease in bare ground cover by 1.16 million km2,\nrepresenting a 3.1% reduction relative to 1982 levels. However, despite this decline in bare\nground cover, its contribution to global LAI changes over the past four decades remains minimal\n(-0.3%) in our analysis. This is primarily because the LAI values associated with bare ground\n(𝛽 <0.6 m2 m-2) are significantly lower than those of tree cover (𝛽 ) and short vegetation\n𝐵𝐺 𝑇𝐶\n(𝛽 ) (Fig. R4, also Supplementary Fig. 5 in the revised manuscript). As a result, even though\n𝑆𝑉\nthe extent of bare ground has decreased, its overall impact on LAI trends is limited. We\nappreciate your insightful comment and have clarified this point in the revised manuscript (Page\n5, Lines 124-126).\ntree cover (𝜷 ), short vegetation (𝜷 ), and bare ground (𝜷 ) during the period of 1982-\n𝑻𝑪 𝑺𝑽 𝑩𝑮\n2015 at 3° × 3° grid spacing.\nFinally, thank you for reviewing our paper and providing useful comments and suggestions. We\nhave acknowledged this in the paper: “We acknowledge the anonymous reviewers for detailed\nand helpful comments on the original manuscript.” (Page 21, Lines 567-568).\nReferences\nBurns, P., Hakkenberg, C. R. & Goetz, S. J. Multi-resolution gridded maps of vegetation\nstructure from GEDI. Sci. Data 11, 881 (2024).\nChen, C., et al. China and India lead in greening of the world through land-use management.\nNat. Sustain. 2, 122–129 (2019).\nCortés, J. et al. Where are global vegetation greening and browning trends significant? Geophys.\nRes. Lett. 48, e2020GL091496 (2021).\nDubayah, R.O., S.B. Luthcke, T.J. Sabaka, J.B. Nicholas, S. Preaux, and M.A. Hofton. 2021.\nGEDI L3 Gridded Land Surface Metrics, Version 1. ORNL DAAC, Oak Ridge, Tennessee,\nUSA. https://doi.org/10.3334/ORNLDAAC/1865\nDuncanson, and J.R. Kellner. 2022. GEDI L4B Gridded Aboveground Biomass Density,\nVersion 2. ORNL DAAC, Oak Ridge, Tennessee, USA.\nhttps://doi.org/10.3334/ORNLDAAC/2017\nFang, H., Baret, F., Plummer, S., & Schaepman-Strub, G. An overview of global leaf area index\n(LAI): Methods, products, validation, and applications. Rev. Geophys. 57, 739–799 (2019).\nForzieri, G., Alkama, R., Miralles, D. G. & Cescatti, A. Satellites reveal contrasting responses\nof regional climate to the widespread greening of Earth. Science 356, 1180–1184 (2017).\nPiao, S. et al. Characteristics, drivers and feedbacks of global greening. Nat. Rev. Earth Environ.\n1, 14–27 (2020).\nSong, X.-P. et al. Global land change from 1982 to 2016. Nature 560, 639–643 (2018).\nZhu, Z. et al. Greening of the Earth and its drivers. Nat. Clim. Change 6, 791–795 (2016).\nReviewer #2 General Comments: Aiming to evaluate the role of land cover change (LCC) in\ncontributing to earth’s greening phenomenon observed by satellite remotely sensed leaf area\nindex (LAI) increasing, the authors decomposed the trend of LAI into LCC-based signals\nincluding changes in tree cover, short vegetation, and bare ground. Then, the authors claimed\nthat 18% of greening since 1980s can be attributed to LCC. The topic is interesting, but I think\nthe current version lacks many details on the analysis and the discussion.\n[Response] Thank you for your thoughtful review of our manuscript. We appreciate your\nrecognition of the importance of our study in quantitatively decomposing the contribution of\nvegetation cover changes to global LAI trends. This perspective was also acknowledged as\ntimely and interesting by Reviewer #1 and Reviewer #3. In response to your insightful\ncomments, we have thoroughly revised the Introduction and Implications sections to more\nclearly highlight the novelty of our research. We have also refined the motivation of the study\nto provide a clearer framework for understanding the critical role of vegetation cover change\n(VCC) in explaining global LAI changes. Additionally, we have included more detailed\nvalidation analyses to enhance the reliability and clarity of our model reconstructions. We\nbelieve these revisions have significantly improved the quality and completeness of our study,\nand we hope you will find that all the issues raised have been adequately addressed in the\nrevised manuscript.\nSpecific comments:\n[Reviewer #2 Comment 1]: First, I may not agree that the consistency of reconstructed trend\nof σLAI and 𝛿LAI can be used as evidence that the model is well representative (L99-101). Fig\n2a shows mean global average, not grid-specific model performance. Assume a naïve model\nthat only predicts the mean LAI for each grid, the same result of consistency can also be\nachieved. Particularly, R2=0.997 is suspectable. Plus, SI Fig. 3 shows that average R2 is about\n0.45, so the model on average explains 45% of LAI variance, which suggests that other major\nfactors are lacking in the model.\n[Response] We sincerely appreciate your thoughtful feedback and apologize for any confusion\n12\nin our original explanation. It is important to clarify that in our Bounded-Variable Least Squares\n(BVLS) model, the R2 value represents the model’s goodness of fit in estimating vegetation\ngreenness parameters. Specifically, it quantifies the proportion of variance explained by BVLS\nwhen decomposing the LAI signal at the 0.05° × 0.05° pixels. However, this R2 does not directly\nreflect the accuracy of the final reconstructed LAI (𝜎𝐿𝐴𝐼 and 𝛿𝐿𝐴𝐼). Therefore, using it to\nquestion the reliability of our reconstructed LAI trends may not be appropriate. To ensure the\nrobustness of our results, we conducted a series of validations of the reconstructed LAI,\ndemonstrating its high accuracy (Figs. R5-R7; Supplementary Table 1). Additionally, in our\nanalysis, we applied a filtering process based on the BVLS R2 values, excluding grid cells where\nthe model fit was not statistically significant (Pages 14-15, Lines 393-397 in the clean version\nof revised manuscript). This step helps to further refine the reliability of our estimates.\nTo validate our statistical models, we conducted extensive comparisons of reconstructed 𝜎𝐿𝐴𝐼\nand 𝛿𝐿𝐴𝐼 with satellite-derived observations at grid, regional, and global scales. While our\nmain text presents the global-scale validation, we also assessed model performance at finer\nspatial scales, including countries, regions, and continents. These evaluations demonstrated\nstrong agreement, with R2 values exceeding 0.97 and relative biases below 6% (Supplementary\nTable 1). Additionally, robustness analyses confirmed the reliability of our models, as similar\nresults were obtained using alternative indices and LAI products (Supplementary Figs. 10-13).\nIn response to your concern, we have now provided more detailed analyses of grid-specific\nmodel performance. At the grid level, the reconstructed annual 𝜎𝐿𝐴𝐼 exhibits high correlation\ncoefficients (𝑟>0.99 for over 95.5% of grids; Fig. R5a) and minimal bias (<±0.03 m2 m-2 for\nover 95.3% of grids; Fig. R5b). Slight overestimations occur mainly in semi-arid regions such\nas Central Asia and the Sahel, while minor underestimations are observed in tropical areas.\nSimilarly, the reconstructed 𝛿𝐿𝐴𝐼 shows overall high accuracy at the grid level, with more\nthan 90% of grids exhibiting a bias of less than ±0.01 m2 m-2 decade-1, and particularly low bias\nin grids with significant LAI changes (Fig. R6). For the three alternative LAI products, the\nreconstructed 𝛿𝐿𝐴𝐼 also exhibits a low bias of less than ±0.01 m2 m-2 decade-1 in 94.9%, 94.4%,\nand 81.6% of the grids (Fig. R7).\nof our analysis. In the revision, we have incorporated your suggestions for additional details,\nincluding in the Main text (Page 4, Lines 99-101; Pages 4-5, Lines 112-118), Methods (Page 14,\nLines 385-391), and Validation (Pages 17-18, Lines 458-488), which have enriched the depth\nand breadth of our analysis.\nFig. R5 (also Supplementary Fig. 6 in the revised manuscript). Validation of the model’s\nreproduction of 𝝈𝑳𝑨𝑰 at grid level. a, Spatial pattern of Pearson correlation coefficient (r)\nbetween the reconstructed annual LAI (𝜎𝐿𝐴𝐼) and the observations from 1982 to 2015 at 3° ×\n14\n3° grid spacing, with the percentage of grids falling in each bin of r values labeled above the\ncolor bar. The correlation coefficients are statistically significant at P <0.001 for all valid grids.\nb, Spatial pattern of mean bias of the reconstructed 𝜎𝐿𝐴𝐼 compared with satellite-derived\nobservations.\nFig. R6 (also Supplementary Fig. 7). Validation of the model’s reproduction of 𝜹𝑳𝑨𝑰 at\ngrid level. Spatial pattern of the bias in the reconstructed 𝜹𝑳𝑨𝑰 at 3° × 3° grid spacing. Dots\nrepresent statistically significant LAI trend at P <0.01 derived from satellite observations. The\ninset displays the frequency distribution of the bias in the reconstructed 𝜹𝑳𝑨𝑰.\ngrid level. Same as Fig. R6, but for the model reproduced 𝜹𝑳𝑨𝑰 derived from a, the CDR data,\nb, the GLASS data, and c, the GLOBMAP data.\n16\n[Reviewer #2 Comment 2]: Second, related to the above problem, without controlling for other\nmajor factors you mentioned in L46-48, it may not be sound to claim that LCC alone contributes\n18% to LAI greening, especially since this 𝛿LAI was estimated by the model not observation\nand you are not comparing the magnitude of LCC with other factors.\n[Response] Thank you for your valuable comment. Based on your concerns and feedbacks, we\nhave refined the explanation of our statistical approach to provide a clearer framework that\nbetter articulates our research objectives. Specifically, we have clarified:\n1) How we control for other factors to isolate the effect of vegetation cover change (VCC).\nIn our statistical framework, we aimed to precisely isolate the contributions of VCC to LAI\nchanges by dividing the analysis into two distinct categories: (1) LAI change signals\ndirectly related to VCC (denoted as 𝜹𝑽𝑪), and (2) LAI trends in areas where vegetation\ncover remained unchanged, representing the influence of other indirect factors (calculated\nas the residual 𝜹𝑰𝑫 = 𝛿𝐿𝐴𝐼 −𝛿𝑉𝐶 ). These other factors, including CO fertilization,\n2\nclimate change (e.g., temperature, water availability, solar radiation), and nitrogen\ndeposition, are closely linked to the LAI change of existing vegetation, rather than being\nprimarily driven by VCC (Zhu et al., 2016; Chen et al., 2019; Piao et al., 2020). The effects\nof these factors are thus generally captured in the 𝜹𝑰𝑫 component.\nIn our study, we specifically focused on quantifying the contribution of VCC to long-term\nchanges in LAI (𝜹𝑽𝑪 ) using a statistical attribution framework that incorporates all\navailable and reliable satellite observations. This approach ensures that our analysis\naccurately reflects the contribution of VCC—a direct and observable consequence of global\nchange. Meanwhile, we emphasize the potential of human VCC-related actions in shaping\nLAI dynamics and highlight the promising effects of forestation programs on the\nmaintenance of ecosystem services.\nTo better clarify this point, we have revised our conceptual framework (Fig. R1) to clearly\nillustrate this approach and better communicate the objectives of our analysis. We have also\nadded a new figure (Fig. R8) that directly displays the spatial pattern of 𝜹𝑰𝑫. Additionally,\n17\nthe Introduction has been updated to more thoroughly explain our motivation and the\ndesign of our statistical framework (Pages 3-4, Lines 66-107), while the Implications\nsection now includes a discussion associated with isolating the influence of LCC (Page 10,\nLines 275-284).\nFig. R8 (also Supplementary Fig. 9). Spatial pattern of reconstructed LAI change. a,\nSpatial pattern of reconstructed 𝜹𝑰𝑫, representing LAI trends induced by indirect factors in\nareas where vegetation cover remains unchanged, calculated as the residual (𝜹𝑰𝑫 = 𝛿𝐿𝐴𝐼 −\n𝛿𝑉𝐶). b, Spatial pattern of reconstructed 𝜹𝑳𝑨𝑰, mirroring the pattern show in Supplementary\nFig. 2 but at a 3° × 3° grid spacing. Dots represent statistically significant trends of LAI changes\n(P <0.01).\n18\n2) How 𝛿LAI was reconstructed by the statistical approach to closely mirror the observations.\nAs outlined in Fig. R1, our statistical framework consists of two main steps. The first step\ninvolves decomposing the greenness parameters (𝛽 , 𝛽 , 𝛽 ) for three different land-\n𝑇𝐶 𝑆𝑉 𝐵𝐺\ncover types using satellite-derived LAI products and the Vegetation Continuous Fields\n(VCF) data (Equations 1-3 in the Methods section). The second step involves reconstructing\nthe annual LAI values (𝜎𝐿𝐴𝐼) using these parameters and attributing long-term changes in\nthe reconstructed LAI (𝛿𝐿𝐴𝐼) to vegetation cover change (𝛿𝑉𝐶) (Equations 4-7).\n𝜎𝐿𝐴𝐼 = 𝛽 ∙𝑇𝐶 +𝛽 ∙𝑆𝑉 +𝛽 ∙𝐵𝐺 (4)\n𝑇𝐶 𝑆𝑉 𝐵𝐺\n𝜕𝜎𝐿𝐴𝐼 𝜕𝛽 𝜕𝑇𝐶 𝜕𝛽 𝜕𝑆𝑉 𝜕𝛽 𝜕𝐵𝐺\n𝛿𝐿𝐴𝐼 = = 𝑇𝐶 ∙𝑇̅̅̅𝐶̅+ ∙𝛽̅̅̅̅̅+ 𝑆𝑉 ∙̅𝑆̅̅𝑉̅+ ∙𝛽̅̅̅̅̅+ 𝐵𝐺 ∙̅𝐵̅̅𝐺̅+ ∙𝛽̅̅̅̅̅\n𝜕𝑦 𝜕𝑦 𝜕𝑦 𝑇𝐶 𝜕𝑦 𝜕𝑦 𝑆𝑉 𝜕𝑦 𝜕𝑦 𝐵𝐺\n(5)\n𝛿𝑉𝐶 = 𝜕𝑇𝐶 ∙𝛽̅̅̅̅̅+ 𝜕𝑆𝑉 ∙𝛽̅̅̅̅̅+ 𝜕𝐵𝐺 ∙̅𝛽̅̅̅̅ (6)\n𝑇𝐶 𝑆𝑉 𝐵𝐺\n𝜕𝑦 𝜕𝑦 𝜕𝑦\n𝛿𝐼𝐷 = 𝜕𝛽𝑇𝐶 ∙𝑇̅̅̅𝐶̅ + 𝜕𝛽𝑆𝑉∙𝑆̅̅̅𝑉̅+ 𝜕𝛽𝐵𝐺 ∙𝐵̅̅̅𝐺̅ (7)\n𝜕𝑦 𝜕𝑦 𝜕𝑦\nWhile it is true that the greenness parameters are derived using a statistical model (the\nBVLS model), and that this disentanglement process (the first part of the framework) may\nintroduce some residuals and estimation biases, both the reconstructed 𝜎𝐿𝐴𝐼 and the\n𝛿𝐿𝐴𝐼 (derived from a differentiated of 𝜎𝐿𝐴𝐼 with respect to the year) are underpinned by\nsolid statistical and physical principles (the second part of the framework). The derivations\nfrom Equations 4 to 7 ensure the reliability of our approach, confirming that we accurately\nisolate the contribution of VCC while controlling other factors influencing LAI (isolating\n𝛿𝑉𝐶 from 𝛿𝐿𝐴𝐼). The Methods section has been updated to more thoroughly explain the\ndesign of our statistical framework (Pages 12-16, Lines 342-445).\nTo validate the accuracy of our reconstructed 𝜎𝐿𝐴𝐼 and 𝛿𝐿𝐴𝐼, we compared them with\nsatellite-derived observations, as shown in Figs. R5-7 (response to Reviewer #2 Comment\n1). The comparison demonstrates that the discrepancies between the reconstructed LAI and\nthe satellite data are minimal and well within acceptable margins, further showing that our\n19\napproach provides a robust reconstruction of LAI that is highly consistent with actual\nsatellite observations.\n3) How we compare the magnitude of 𝜹𝑽𝑪 with other factors.\nWe fully acknowledge the importance of other contributing factors, such as CO\n2\nfertilization, climate change, and nitrogen deposition, and their impact on LAI trends.\nHowever, disentangling the individual effects of these factors is challenging to implement\ngiven the scope and limitations of the current data. Therefore, our focus is on accurately\nestimating the contribution of VCC vs other indirect factors to LAI changes, while\nrecognizing the difficulties in separating the effects of other drivers (Pages 9-10, Lines 254-\n268).\nHere we have also provided the spatial patterns (Fig. R8a) and trends (Fig. R9) illustrating\nthe contributions of other factors to global LAI changes (𝛿𝐼𝐷 ). LAI changes while\nvegetation cover remains constant are primarily influenced by indirect factors. LAI changes\nfor tree cover ( 𝜕𝛽𝑇𝐶 ∙𝑇̅̅̅𝐶̅ ), short vegetation ( 𝜕𝛽𝑆𝑉∙𝑆̅̅̅𝑉̅ ), and bare ground ( 𝜕𝛽𝐵𝐺 ∙𝐵̅̅̅𝐺̅ )\n𝜕𝑦 𝜕𝑦 𝜕𝑦\n(Equation 7) account for 27.9%, 54.8%, and -0.8% of 𝛿𝐿𝐴𝐼, respectively. The combined\neffect of these factors leads to a significant increase in 𝛿𝐼𝐷 at a rate of 0.076 ±0.014 m2\nm-2 decade-1 (P <0.001), explaining 81.9 ±15.1% of the LAI change. The results highlight\nthe dominant influence of indirect factors on global LAI changes over the past four decades,\nwhile emphasizing the significant and growing contribution of VCC-related actions to the\nincrease in LAI. We hope these clarifications address your concerns. Thank you for your\ninsightful comments, which have helped to improve the clarity and rigor of our work.\n**P <0.01; *P <0.05; n.s., not significant. Error bars show the standard deviation.\n[Reviewer #2 Comment 3]: Third, it seems confusing to me that how this study solves the\nproblems described in L60-76. I would appreciate more details and discussions on previous\nstudies and how and why this analysis is better. Otherwise, it would just be another method\nderived result, not necessarily a better one, and the consistency/inconsistency with previous\nresults would be meaningless. This is where I think is very important to clarify before talking\nabout any results.\n[Response] We thank the reviewer for the feedback, which has prompted us to re-evaluate and\nstrengthen our arguments. Importantly, our research distinguishes itself from previous works,\nincluding Zhu et al. (2016) and Chen et al. (2019), by providing the first quantitative,\nobservational-based evidences on how, when, and where vegetation cover changes (VCC)\nhave influenced global LAI changes—a crucial issue that has not been explored before. Below\nwe provide detailed accounts on the novelty of our research in comparison to previous works.\n1) We provide the most reliable observational data-driven evidence.\nZhu et al. (2016) analyzed the key drivers of global LAI increase using three long-term\nsatellite LAI records and ten global ecosystem models. They concluded that CO\n2\nfertilization is a major factor driving global vegetation greening, with the effects of VCC\n21\ncontributing only 4% to the observed trend. However, their findings were based on the\nimplicit assumption that the models accurately capture the physiological sensitivity of\nvegetation to environmental changes, whereas the model sensitivity is not entirely correct\n(Bonan & Doney, 2018; Jiang et al., 2019). Moreover, several studies have pointed out that\nthe models used in previous attribution analyses often provided simplified representations\nof the complex, evolving patterns of land-use practices (Prestele et al., 2017; Chen et al.,\n2019; Piao et al., 2020). As a result, these models likely under-represented or overlooked\ncritical physiological processes, leading to an underestimation of VCC’s influence and an\noveremphasis on the benefits of indirect factors.\nTo address this gap, our research presents a comprehensive, data-driven statistical analysis\nthat provides reliable evidence of how, when, and where VCC contributes to global-to-\nregional LAI changes. Unlike previous work, our approach focuses on disentangling VCC\nfrom other factors using actual observational data rather than relying on model projections.\nThis is a significant advancement, as highlighted in our response to Reviewer #2 Comment\n2.\n2) We quantitatively attribute LAI changes in response to VCC.\nIn a subsequent high-impact study, Chen et al. (2019) utilized field and satellite observations\nto highlight a distinct overlap between croplands and observed LAI increases, particularly\nin regions like China and India. They inferred that human agricultural management could\nexplain more than 33% of global LAI changes. However, their study was qualitative in\nnature, and we have reservations regarding their assumption that all changes in LAI within\nVCC areas can be solely attributed to land management, without accounting for the potential\ninfluence of indirect factors. This gap highlights the need for a quantitative investigation\nthat can more precisely attribute LAI changes to VCC.\nIn contrast to previous approaches, such as those used by Chen et al. (2022), which subtract\nlocal signals from nearby non-VCC areas using “Paired Land Use Experiment”, or the\nassumption in Chen et al. (2019) that all LAI changes in VCC areas are solely due to land\n22\nmanagement, our study employs a novel, data-driven statistical attribution framework. This\napproach disentangles greenness parameters across different vegetation cover types and\nthen reconstructs LAI (𝛿𝐿𝐴𝐼 ) while attributing the changes to VCC (𝛿𝑉𝐶 ), explicitly\nexcluding the greening effects from areas that have not undergone VCC (Equations 4-7).\nBy utilizing all available long-term LAI products and dynamic VCF data, our framework\nprovides the first robust, quantitative analysis that disentangles the contributions of VCC\nfrom mixed signals across various scales. The grid-specific validations further demonstrate\nthat our approach provides a robust reconstruction of LAI that is highly consistent with\nactual satellite observations (Figs. R5-7; response to Reviewer #2 Comment 1).\nIn response to your insightful comments, we have revised the manuscript to more clearly\narticulate the motivation behind our analysis and highlight the novelty of our approach in\nrelation to previous studies (Pages 3-4, Lines 66-107). Specifically, we have emphasized the\nkey distinctions between our work and that of Zhu et al. (2016) and Chen et al. (2019), with\ndetailed explanations provided in the revised manuscript (Page 9, Lines 239-252; Page 10, Lines\n270-275). We sincerely appreciate your constructive feedback, which has greatly contributed to\nenhancing the clarity and depth of our research by addressing important gaps and highlighting\nthe unique aspects of our approach.\n[Reviewer #2 Comment 4]: L79. “high-resolution” turns out to be 0.05° (~5 km) and then\naggregated to 3° (~300 km)? Although AVHRR is named “high-resolution”, that was due to\nsome historical issue. I’d recommend be more specific here.\n[Response] Thank you for your suggestion. We acknowledge that the term high-resolution may\nbe ambiguous in this context. In our study, we first decomposed greenness parameters using\ndata at the original 0.05° × 0.05°-pixel scale, which aligns with the highest available resolution\nof the AVHRR-based LAI products. To enhance clarity, we have revised the text to explicitly\nstate the resolutions and removed the term “high-resolution” (Page 4, Lines 87-89).\n[Reviewer #2 Comment 5]: L83. Can you describe some details here about what the improved\nBounded-Variable Least Squares approach is and how it solves the challenges you mentioned\n23\nin the previous paragraph?\n[Response] The Bounded-Variable Least Squares (BVLS) approach is an optimization method\ndesigned to solve least-squares problems while incorporating constraints on variable values\n(Stark & Parker, 1995). In our study, we extended the traditional BVLS approach by not only\ndefining value ranges for the three greenness variables (𝛽 , 𝛽 , 𝛽 ) but also imposing\n𝑇𝐶 𝑆𝑉 𝐵𝐺\nadditional constraints on their relative magnitudes (0 < 𝛽 ≤ 𝛽 ≤ 𝛽 ≤ 𝛽 ; Equation 3;\n𝐵𝐺 𝑆𝑉 𝑇𝐶 𝑚𝑎𝑥\nSupplementary Fig. 19). These improvements ensured that the estimated parameters adhered to\nphysically meaningful relationships, enhancing the stability and reliability of the results,\nespecially under conditions of noise or multicollinearity.\nUsing this improved BVLS approach, we derived the three greenness variables (Equations 1-3)\nand reconstructed 𝛿𝐿𝐴𝐼 across 3° × 3° grid, regional and global scales. This allowed us to\nquantify the attributions of LAI changes in response to vegetation cover change (𝛿𝑉𝐶) in a\nconsistent and systematic manner using a differentiated equation framework (Equations 4-6).\nWe have clarified this in the revised manuscript to better explain the advantages and\nimplementation of our improved BVLS approach (Page 4, Lines 87-94; Page 13, Lines 350-\n357; Page 14, Lines 375-385).\nAdditionally, we have included a revised version of our conceptual framework (Fig. R1, also\nSupplementary Fig. 1) to clarify our improved BVLS approach and better convey the main\nobjectives of our analysis.\n[Reviewer #2 Comment 6]: L90-91. This explanation is a bit abrupt because you did not say\nmax LAI was used before this sentence.\n[Response] To improve clarity, we have revised the manuscript by introducing the use of\nmaximum LAI earlier in the paragraph, making it clear from the beginning why it was chosen.\nSpecifically, we now mention the use of maximum monthly LAI in the second sentence when\ndescribing the methodology, ensuring a smoother logical flow (Page 4, Lines 85-87).\n[Reviewer #2 Comment 7]: L162-166. SI Fig. 1 does not talk about these factors, e.g., CO2,\n24\nclimate change, right? How did you get this conclusion?\n[Response] Thank you for pointing this out, and we apologize for the confusion. Our intention\nwas to illustrate that in regions with severe deforestation, such as Brazil, Argentina, Madagascar,\nand the Miombo woodlands in the southern part of continental Africa, the negative impact of\nvegetation cover change on LAI trends (𝛿𝑉𝐶 ) can be masked by the positive influence of\nindirect factors such as CO₂ fertilization and climate change in areas with constant vegetation\n(i.e., where land cover remains unchanged) (𝛿𝐼𝐷 ). This can be observed by comparing the\nprevious version of Fig. S1 with Fig. 2d (now revised as Supplementary Fig. 2 and Fig. 1d).\nTo better clarify this point, we have added a new figure (Fig. R8, also Supplementary Fig. 9 in\nthe revised manuscript) that directly displays the spatial pattern of 𝛿𝐼𝐷 (Fig. R8a, representing\nthe influence of indirect factors on the LAI trend) and the reconstructed 𝛿𝐿𝐴𝐼 (Fig. R8b). It\ncan be observed that in these regions, the magnitude of 𝛿𝐼𝐷 tend to be high. As a result, the\noverall LAI trend appears to be positive, despite severe deforestation having occurred.\nIn response, we have updated the corresponding description and figure references in the text:\n“Some LAI increases in these areas, largely resulted from indirect factors4 such as CO\n2\nfertilization, climate change and nitrogen deposition3,7 (𝛿𝐼𝐷; Supplementary Fig. 9a), may have\noffset the deforestation-induced decrease in LAI at a regional scale over the past four decades\n(Fig. 1d versus Supplementary Fig. 9b). However, the persistence of these benefits amidst\nescalating disturbances cannot be guaranteed.” (Page 7, Lines 184-188)\n[Reviewer #2 Comment 8]: L227-228. Can you be more specific? Why do you think previous\nefforts underestimated LCC? And why do you think your method is better?\n[Response] Thank you for point this out. We have revised this sentence to clarify the potential\nunderestimation associated with the results from modeling studies and have included relevant\nreferences: “Ecosystem models likely underestimate the influence of VCC due to their omission\nor underrepresentation of the complex patterns of land-use practices and critical physiological\nprocesses in ecosystem models2,4,22,25, which lead to an overemphasis on the benefits of indirect\nfactors.” (Page 10, Lines 272-275). In contrast, our approach, which utilizes observational data\n25\nwithin a statistical framework that accounts for the complexity dynamics of VCC, provides a\nmore accurate estimate of its contribution to global LAI changes. A more detailed explanation\nis provided in the response to Reviewer #2, Comment 3.\n[Reviewer #2 Comment 9]: Fig. 1. “soving” should be “solving”.\n[Response] We have corrected the spelling error in this figure.\n[Reviewer #2 Comment 10]: L330-332. I don’t think I understand the “substantial lack of LAI\nobservations”, can you give more descriptions?\n[Response] Thank you for raising this concern. We have revised this sentence to clarify current\nchallenges and our motivation: “Developing such a framework is motivated by the lack of LAI\nobservations for specific vegetation types, as existing LAI data are grid-based at a resolution of\n~5 kilometers, with each grid cell representing a single value that averages across multiple\nvegetation types. This means the lack of information on individual vegetation types for LAI,\nmaking it difficult to quantify how sub-grid scale vegetation cover changes contribute to overall\nLAI variations.” (Pages 12-13, Lines 344-348). We have also incorporated this challenge into\nthe revised Supplementary Fig. 1 (Fig. R1) to clarify our approach and better convey the main\nobjectives of our analysis.\n[Reviewer #2 Comment 11]: L370-371. SI Fig. 4 does not support “we successfully\ndisentangled…” as they are basically model estimates, which do not say whether the model is\ngood or bad.\n[Response] Thank you so much for your comment. As mentioned in our Response to your\nComment 1, we have now provided more detailed analyses of grid-specific model performance\n(Figs. R5-7). The comparison demonstrates that the discrepancies between the reconstructed\nLAI and the satellite data are minimal and well within acceptable margins, further showing that\nour approach provides a robust reconstruction of LAI that is highly consistent with actual\nsatellite observations.\nFinally, thank you again for reviewing our paper and providing valuable suggestions. We have\n26\nacknowledged this in the paper: “We acknowledge the anonymous reviewers for detailed and\nhelpful comments on the original manuscript.” (Page 21, Lines 567-568).\nReferences\nBonan, G. B. & Doney, S. C. Climate, ecosystems, and planetary futures: The challenge to\npredict life in Earth system models. Science 359, eaam8328 (2018).\nChen, C., et al. China and India lead in greening of the world through land-use management.\nNat. Sustain. 2, 122–129 (2019).\nChen, T. et al. Land management explains the contrasting greening pattern across China-Russia\nborder based on Paired Land Use Experiment approach. J. Geophys. Res. Biogeosciences\n127, e2021JG006659 (2022).\nJiang, M., Caldararu, S., Zaehle, S., Ellsworth, D. S. & Medlyn, B. E. Towards a more\nphysiological representation of vegetation phosphorus processes in land surface models.\nNew Phytol. 222, 1223–1229 (2019).\nPiao, S. et al. Characteristics, drivers and feedbacks of global greening. Nat. Rev. Earth Environ.\n1, 14–27 (2020).\nPrestele, R. et al. Current challenges of implementing anthropogenic land-use and land-cover\nchange in models contributing to climate change assessments. Earth Syst. Dynam. 8, 369–\n386 (2017).\nStark, P. B. & Parker, R. L. Bounded-variable least-squares: an algorithm and applications.\nComput. Stat. 10, 129–129, (1995).\nZhu, Z. et al. Greening of the Earth and its drivers. Nat. Clim. Change 6, 791–795 (2016).\nReviewer #3 General Comments: Reports of global greening are widespread and credible,\nbut the contribution of land cover change has been uncertain, and create substantial confusion\non the projected efficacy of natural climate solutions. By using vegetation continuous fields\ndata to statistically disentangle land cover change from change within a land cover class (e.g. increases in forest productivity that may be attributable to CO fertilization), Zeng et al. provide\n2\ncredible data on how, when, and where land cover changes have influenced global greening.\nThe clear methods, robust analysis of multiple datasets at varying resolutions, disaggregated\nreporting of quantitative results (including in helpful detailed supplemental information), and\nthe effective graphics and writing are major strengths of the manuscript. The analysis is\nconvincing and comprehensive in its presentation and quality. These strengths are important\nbecause, as the authors rightly detail in the concluding section, the results of this analysis are\nvital for guiding natural climate solutions.\nIn my limited suggestions to improve the paper, I only found the few minor grammatical or\nstylistic errors below.\n[Response] We are truly grateful for your encouraging feedback on our manuscript. We greatly\nappreciate your high compliments on the broad interest of our paper and detailed feedback on\nour manuscript. We have carefully addressed the grammatical and stylistic suggestions you\nprovided. These revisions have been incorporated into the revised manuscript to further enhance\nits clarity and readability. Details of the revision are provided below.\nSpecific comments:\n[Reviewer #3 Comment 1]: Line 132: should be “Contrasting effects...”\n[Response] Thanks, we have modified this sub-title in the revision (Page 6, Line 153 in the\nclean version of revised manuscript).\n[Reviewer #3 Comment 2]: Line 226-227: Tone down to read: “This underestimation likely\noriginates...”\n[Response] We have revised this sentence accordingly (Page 10, Line 272-275).\n28\n[Reviewer #3 Comment 3]: Line 231: I assume the 60% statistic is from the cited study, but it\nis unclear why that statistic from another study directly applies to this study. Perhaps avoid the\nconfusion of that precise statistic and tone down to read: “A majority of the LCC in this study\nis linked to direct human activities...”\n[Response] Thank you for pointing this out. We agree that the specific statistic may cause some\nconfusion regarding its direct applicability to our study. Follow your suggestions, the revised\nsentence now reads: “While VCC’s influence on LAI dynamics is evident, with a majority of\nVCC in this study linked to direct human activities14” (Page 10, Lines 275-276).\n[Reviewer #3 Comment 4]: Line 232: This sentence is awkward to read. Starting with “Besides,”\nis perhaps the problem, but the whole sentence is oddly phrased. Reword.\n[Response] Thank you for your valuable feedback. We have thoroughly rewritten the\nImplications section to better highlight the novelty of our research and the critical policy\nimplications of the growing influence of vegetation cover change (VCC) on LAI dynamics. As\npart of these revisions, this original sentence has been removed.\n[Reviewer #3 Comment 5]: Line 238: Change to the plural “ecosystem services”\n[Response] We have revised the text as recommended (Page 10, Line 283).\n[Reviewer #3 Comment 6]: Line 271: “recent failure” is singular and time-specific, but\ndeforestation is pervasive and chronic. Tone down to “failures to limit deforestation”.\n[Response] Thank you for pointing this out. We have revised the phrase “recent failure” to\n“unchecked deforestation, agriculture expansion, urbanization, and expanding energy\ninfrastructure” in order to better reflect the ongoing intensification of ecosystem disturbances\n(Page 10, Lines 267-268).\nFinally, thank you again for reviewing our paper and providing valuable suggestions. We have\nacknowledged this in the paper: “We acknowledge the anonymous reviewers for detailed and\nhelpful comments on the original manuscript.” (Page 21, Lines 567-568)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Overall, the authors have done a good job of addressing the main concerns I raised in my original review. I think the reframing around LAI change (instead of greening vs browning) is much clearer. It’s also good to see that the authors have discussed some of the caveat with estimating LAI from satellite data and the comparison to the GEDI data is generally encouraging.", "reviewer_2": "I really appreciate that you made the effort to address my comments. Regarding my initial first comment, I may still not be convinced. But it could be due to my misunderstanding, I hope the authors can further make it clearer in the manuscript since if I misunderstood the method, likely other readers will too.\nI am sorry I may not understand the response \"it (R2) quantifies the proportion of variance explained by BVLS when decomposing LAI signal at the 0.05° x 0.05° pixels. However, this R2 value does not directly reflect the accuracy of the final reconstructed LAI\". Isn't that the definition of R2? And, I thought you aggregated the 0.05° pixels into 0.1° pixels (L360-362) and then fit the model, didn't you?\nTo my understanding, you fit a model w/ all the 0.1° pixels in a 3° grid. The R2 was calculated for the grid. Then, you used the model to predict a single LAI value for that grid with averaged predictor values. So, the model didn't change but the predictor values were averaged. How come the original model has an R2=0.5 yield reconstructed 16-year LAIs with a correlation coefficient (SI Fig. 6) of 0.99? Could it be that the reconstructed annual LAI for that 3° grid is just very close to the mean LAI of all 0.1° pixels? Then, the correlation coefficient between the 16-year \"mean\" values calculated by the models and the mean values from the observations got very high? In this sense, maybe the correlation coefficient at this 3° grid isn't a good metric for proving accuracy here? To me, the histogram in SI Fig. 4 makes more sense than SI Fig. 6. Why don't you map the R2 values in SI Fig.4 so readers would know where the model is more accurate?\nI am exaggerating the problem and I am sure the model works at some level. Still, my main point is that I hope the method validation can be clearer and the uncertainty be reasonably reported.\nBtw, a minor comment: L37-38. I understand the Abstract needs to be short, but do you think further explaining this sentence a bit (e.g., how the sensitivity is defined? Can you provide some quantitative values?) would be helpful for readers to understand?", "reviewer_3": "Thank you for addressing my minor comments from my positive earlier review."}, "author_response": "Response to the reviewers (NCOMMS-24-44410A)\nReviewer #1 (Remarks to the Author):\nReviewer #1 General Comments: Overall, the authors have done a good job of addressing\nthe main concerns I raised in my original review. I think the reframing around LAI change\n(instead of greening vs browning) is much clearer. It’s also good to see that the authors have\ndiscussed some of the caveat with estimating LAI from satellite data and the comparison to the\nGEDI data is generally encouraging.\n[Response] We sincerely appreciate your recognition of the revisions we made. Your comments\nwere instrumental in helping us improve the clarity and rigor of the manuscript, and we are\ngrateful for your constructive guidance throughout the review process.\nReviewer #2 (Remarks to the Author):\nReviewer #2 General Comments: I really appreciate that you made the effort to address my\ncomments. Regarding my initial first comment, I may still not be convinced. But it could be due\nto my misunderstanding, I hope the authors can further make it clearer in the manuscript since\nif I misunderstood the method, likely other readers will too.\n[Response] Thank you very much for acknowledging our efforts to address the concerns and\nfor taking the time to engage deeply with our work. We recognize that some aspects of our\nmethod may still not have been presented as clearly as intended. In response, we have further\nrevised the relevant section of the manuscript to enhance clarity and make our approach more\ntransparent and accessible to a broader audience.\nSpecific comments:\n[Reviewer #2 Comment 1]: I am sorry I may not understand the response “it (R2) quantifies\nthe proportion of variance explained by BVLS when decomposing LAI signal at the 0.05°×0.05°\npixels. However, this R2 value does not directly reflect the accuracy of the final reconstructed\nLAI”. Isn’t that the definition of R2? And, I thought you aggregated the 0.05° pixels into 0.1°\npixels (L360-362) and then fit the model, didn’t you?\nTo my understanding, you fit a model w/ all the 0.1° pixels in a 3° grid. The R2 was calculated\nfor the grid. Then, you used the model to predict a single LAI value for that grid with averaged\npredictor values. So, the model didn’t change but the predictor values were averaged. How\ncome the original model has an R2=0.5 yield reconstructed 16-year LAIs with a correlation\ncoefficient (SI Fig. 6) of 0.99? Could it be that the reconstructed annual LAI for that 3° grid is\njust very close to the mean LAI of all 0.1° pixels? Then, the correlation coefficient between the\n16-year ‘mean’ values calculated by the models and the mean values from the observations got\nvery high? In this sense, maybe the correlation coefficient at this 3° grid isn’t a good metric for\nproving accuracy here? To me, the histogram in SI Fig. 4 makes more sense than SI Fig. 6.\nWhy don’t you map the R2 values in SI Fig. 4 so readers would know where the model is more\naccurate?\nI am exaggerating the problem and I am sure the model works at some level. Still, my main\npoint is that I hope the method validation can be clearer and the uncertainty be reasonably\nreported.\n[Response] We sincerely thank the reviewer for the thoughtful follow-up and for highlighting\nthis important issue regarding our explanation of the Bounded-Variable Least Squares (BVLS)\nmodel and its validation. We recognize that our previous response lacked clarity, and we greatly\nappreciate the opportunity to clarify it further.\n“I thought you aggregated the 0.05° pixels into 0.1° pixels...”\nYes, thank you for pointing this out. In the disentanglement approach, we first aggregated\nthe 0.05° LAI and VCF data to 0.1°resolution by averaging valid values annually. These\n0.1° cells were then used as input to the BVLS regression within each 3° × 3° grid. We\nhave clarified this processing step in the Methods section (Page 13, Lines 363-369 in the\nclean version of revised manuscript). To ensure robustness, we also tested multiple grid\nspacings (Supplement Table 3). The confusion arose from an imprecise description in our\nprevious response, and we apologize for any misunderstanding it may have caused.\n“Isn’t that the definition of R²?”\nYou are absolutely right. R2 quantifies the proportion of variance explained by the model.\nOur earlier wording was unclear, and we regret the confusion. The R2 from the BVLS\nfitting reflects the model’s performance in capturing 0.1° pixel-level LAI variation\nwithin a 3° × 3° grid, using grid-averaged greenness parameters (𝛽 values). However,\nour final attribution analysis focuses on grid-level mean LAI values, not individual 0.1°\npixels, which may cause a potential disconnect between model fit (R2) and large-scale\npattern reconstruction (correlation with observed 3° LAI). We have now revised the\nMethods section to make this distinction explicit. (Page 14, Lines 393-396).\n“How come the model with R²=0.5 produces LAI with correlation = 0.99?”\nThis is an important point. The R2 value of ~0.5 reflects the model’s ability to predict 0.1°\nLAI values within a 3° grid using grid-level 𝛽 values, assuming 100% coverage of each\nvegetation type within the grid. It measures how well the 𝛽 values capture sub-grid\nvariability. Using these 𝛽 values, we reconstructed annual LAI for each 3° × 3° grid\nby multiplying them with the corresponding grid-averaged vegetation fractions. The\nreconstructed LAI closely matches satellite observations in magnitude and interannual\nvariability, with low bias and high temporal correlations at the grid level (Supplementary\nFigures 6-7). This indicates that while the pixel-level fit may be moderate, the model\nstill captures large-scale signals effectively in terms of magnitude and interannual\nvariability. We have clarified this distinction in the revised manuscript (Page 17, Lines\n476-480).\n“Maybe correlation at the 3° grid isn't a good accuracy metric?”\nWe fully agree with this suggestion. In response, we have added a global map of R2 values\n(i.e., grid-level model performance) at the 3° grid level in Supplementary Figure 4 (also\nFigure R1 below). This provides a more direct assessment of model goodness-of-fit across\nregions. This map shows that the model performs better (i.e., exhibits higher R2 values) in\nregions with LAI increase hotspots and/or significant vegetation cover change (e.g., East\nAsia, South Asia, Europe, Sahel, eastern South America, eastern Australia), while it shows\nlower R2 in areas with dense, uniform vegetation (e.g., the Amazon and Congo forests, and\nhigh-latitude tundra), likely due to limited variability in vegetation composition. We have\nadded this description in the revised manuscript (Page 15, Lines 404-408).\n“I hope the method validation can be clearer and the uncertainty be reasonably reported.”\nWe fully agree that clear validation and uncertainty reporting are essential. In the revised\nmanuscript, we have clarified the distinction between R2 (which reflects model fit at the\n0.1°-pixel level; Lines 393-394) and correlation coefficients (which assess agreement\nbetween reconstructed and observed annual LAI at the 3° grid-level; Lines 470-473). We\nalso revised the Methods section to specify that all decomposition and validation steps are\nperformed at the 3° grid level, and applying the estimated grid-level 𝛽 values to\ndisentangle and attribute LAI at finer scales (e.g., 0.1°) may introduce uncertainty. (Pages\n17-18, Lines 480-482).\nOnce again, we thank the reviewer for highlighting these important points, which has led to\nsubstantial improvement in the clarity of our methodological description and validation\npresentation.\nFig. R1. Spatial distribution of model goodness-of-fit (R2 values) for the improved BVLS\nmodel at 3° × 3° grid spacing. Grids with more than 16 years of statistically significant BVLS\nresults (F-statistic p-value <0.05) over the study period are defined as valid and are shown in\nthe map with their corresponding R2 values. The inset displays the probability density\ndistribution and the mean value of the multi-year average R2 across all valid grids.\n[Reviewer #2 Comment 2]: Btw, a minor comment: L37-38. I understand the Abstract needs\nto be short, but do you think further explaining this sentence a bit (e.g., how the sensitivity is\ndefined? Can you provide some quantitative values?) would be helpful for readers to\nunderstand?\n[Response] Thank you for this helpful suggestion. We have revised the sentence in both\nAbstract (Page 2, Lines 34-37) and main text (Page 7, Lines 192-197) to provide a clear\ndefinition of sensitivity and to report the values across continents. We have also reorganized\nthe abstract to improve clarity and expression while adhering to the word limit.\nReviewer #3 (Remarks to the Author):\nReviewer #3 General Comments: Thank you for addressing my minor comments from my\npositive earlier review.\n[Response] Thank you very much for your kind and supportive comments. We truly appreciate\nyour positive earlier review and the helpful suggestions you provided."}, {"round": 3, "reviewer_comments": {"reviewer_2": "Thank you for addressing my further comments. I think the manuscript is ready to be published, congrasts!"}, "author_response": ""}]

Decision letter