Emergent constraints on the hydrological impacts of land use and land cover change

Domain: chemistry
Published: 2026-02-18
Editorial outcome: major_revision

Abstract
Land use and land cover changes have substantial effects on the terrestrial water cycle, but their sign and magnitude remain elusive at large scales. State-of-the-art Earth system models disagree on how these changes affect terrestrial evapotranspiration. Here we use the observation-based transpiration-specific Bowen ratio to correct modelled evapotranspiration changes induced by land use and land cover changes globally and regionally within a hierarchical emergent constraint framework. We show that the constraint reverses the sign of the original model estimates at the global scale and over Central and South America, and narrows the inter-model spread. The misrepresentation of transpiration-specific Bowen ratio and its variations across plant functional types in models is the main source of this bias. Applying an analogous constraint framework to a future afforestation scenario, the constrained simulations project stronger evapotranspiration enhancements and weaker decreases in terrestrial water availability compared to the original simulations, particularly in tropics and subtropics.

Introduction
About three-quarters of the Earth’s land surface has been altered by human activities like afforestation, deforestation for agricultural expansion, grazing and urbanization1,2. The effects of such anthropogenic land use and land cover change (LULCC) on the terrestrial water cycle are considered of comparable magnitude to those originating from climate change3–7. LULCC directly influences the scale and timing of terrestrial evapotranspiration (ET) due to the inherent differences between plant functional types (PFTs) in growing-season length, root development, leaf traits, and water-use efficiency4,8,9. In addition, it contributes to indirectly controlling precipitation (P) patterns locally and remotely in downwind regions through moisture recycling5,10. Such LULCC-driven processes may lead to shifts in soil moisture and total terrestrial water availability (WA, defined as P minus ET)5,9 and trigger cascading detrimental or beneficial impacts on socio-ecological systems and services crucial for the human well-being, such as biodiversity conservation11, carbon uptake and storage12,13, food security and timber production14. Furthermore, afforestation, one of the most obvious manifestations of LULCC, has a strong capacity to absorb atmospheric CO2 and is considered a key solution to meet climate targets15,16. But the scalability of this land use management option has recently been questioned due to its potential threat to water supply for societies17,18. Therefore, a robust and accurate assessment of the long-term effects of LULCC scenarios on the terrestrial water cycle is of paramount importance for developing sustainable water management and more integrated and realistic land-based climate mitigation strategies.

Numerous studies based on in-situ hydrological measurements have explored the ongoing changes in the terrestrial water cycle in response to LULCC but their results are elusive and conclusions are still largely controversial19. Previous regional and catchment-scale studies have led to contradictory findings, i.e., that forest gain has either a negative20, negligible21, or positive effect22,23 on water yield. These results suggest that the hydrological response to forest change may vary greatly from one catchment to another24 due to different background climate, geomorphological characteristics and soil properties25. Furthermore, it is unclear to what extent the results of such analyses based on ground-based hydrological observations can be extrapolated to larger areas26,27. Repeated and consistent satellite observations can provide a more systematic characterization in space and time of the terrestrial water cycle at the regional to global scales28. However, such observational products can only partially capture the non-local biophysical effects driven by changes in large-scale circulation and the associated impacts on cloud cover and P typically associated to large-scale LULCC29,30. Earth system models (ESMs) as well as land surface models (LSMs, i.e., the land component of ESMs) can provide a comprehensive assessment of the hydrological response to LULCC at the global scale by integrating local and non-local effects31,32. Nevertheless, their simulations lack consensus on the direction and magnitude of the response of key hydrological variables (e.g., ET, P, and runoff) to LULCC, because of the highly heterogeneous representation of key biological and physical processes across ESMs/LSMs33,34. Results from a recent modelling study focusing on the hydrological response to a global deforestation experiment emphasize the substantial uncertainty in hydrological projections under LULCC27. This is particularly critical considering the strong coupling between terrestrial water and carbon cycles12,35,36. Indeed, the resulting uncertainties in projections of hydrological response to LULCC largely hamper a proper quantification of the associated impacts on the terrestrial carbon budget and may ultimately limit the identification and development of effective land-based climate mitigation policies.

Here we use multi-source observations to enhance the representation of the modelled LULCC effect on terrestrial evapotranspiration (δETLULCC) within a hierarchical emergent constraint framework37. Factorial simulations of an ensemble of ESMs participating in the Coupled Model Intercomparison Project Phase 6 (CMIP6)38 are used to disentangle the δETLULCC signal during the 1982-2014 at the global scale (Supplementary Tables 1–3). We apply such an emergent constraint framework to future scenarios to derive more accurate projections of afforestation impacts on terrestrial ET (δETAFF) and WA (δWAAFF) over 2015-2099. We extend the constraint analyses over the reference regions39 used in the Intergovernmental Panel on Climate Change Sixth Assessment Report (IPCC AR6), to provide more spatially detailed information about the LULCC-driven hydrological changes and better assist water management strategies and land-based mitigation plans. Finally, we investigate the dependence of the modelled LULCC impact on the surface energy partitioning across PFTs to elucidate the underlying mechanisms responsible for the inter-model spread.

Results
Large disagreement among ESMs on the LULCC hydrological effects

The CMIP6 model ensemble – under prescribed land use and land cover forcing from the LUH2 dataset40,41 – shows that global forest area has decreased significantly (p<0.01, Mann-Kendall test) by 0.018±0.008 million km2 yr-1 (or 0.016±0.007% yr-1 in area fraction) during the period 1982-2014 (Fig. 1a). Over the same period, global cropland area has increased significantly (p<0.01) by 0.017±0.002 million km2 yr-1 (or 0.014±0.002% yr-1), whereas global grassland area has experienced an increase (0.012±0.010 million km2 yr-1 or 0.010±0.009 % yr-1, p<0.01) and then a decrease (-0.019±0.010 million km2 yr-1 or -0.016±0.010% yr-1, p<0.01) with a change point around year 2000 (Fig. 1a). Overall, global land has experienced a massive conversion from forests to croplands during the past three decades. The inter-model spread in the trends of forest, cropland and grassland area (here expressed in terms of standard deviation) can largely be attributed to the difference in the assumptions embedded within the land-use translator of the ESMs (the program that translates the original LUH2 dataset into cover fractions of trees, crops, and grasses)40. In addition, ESMs such as CESM2 and UKESM1-0-LL explicitly consider the competition between natural vegetation types and these effects may further amplify their difference in terms of forest and grassland areas with respect to other ESMs42,43.

According to an ensemble of historical simulations from twelve CMIP6 models, the LULCC – generally characterized by the conversion from forests to croplands over the period 1982-2014 – has increased terrestrial ET globally (δETLULCC) by an average rate of 0.057±0.155 mm yr-2 (Fig. 1b). Such positive effects on ET are counter-intuitive with regard to the widespread reduction in forest cover during this period (Fig. 1a), raising further concerns about the capacity of ESMs in properly capturing ET dynamics in response to LULCC27. Moreover, the standard deviation of the ensemble is 2.5 times the mean value and further underscores the substantial uncertainty and the large discrepancy amongst single ESM runs. Only seven out of twelve individual ESMs agree on the overall positive sign of δETLULCC, and the spread in magnitude is considerable (global average signal ranging from 0.005 mm yr-2 to 0.380 mm yr-2, Fig. 1b). Similar patterns emerge in terms of the relative contribution of LULCC on historical change in terrestrial ET (13.3±33.4%, Supplementary Fig. 1). These results highlight the lack of consensus in direction and magnitude of δETLULCC across state-of-the-art ESM simulations when the signal is aggregated at the global scale.

The spatially aggregated estimates described above largely mask the variability in the LULCC effect across different regions of the Planet. The ensemble mean δETLULCC is positive in Eastern and Southern Africa, and India, whereas it is negative in Australia, Southern Siberia, and South America (Fig. 1c, d). We point out that the ensemble standard deviation of δETLULCC is more than 5 times the ensemble mean value over 44.7% of the global vegetated land. In contrast, it appears lower than the ensemble mean only in about 0.6% of the domain (Fig. 1e). Such remarkable inter-model spread is more pronounced in high-latitude regions of the Northern Hemisphere (50°N-65°N) and over the tropics of the Southern Hemisphere (5°S-10°S) (Fig. 1e, f). Estimates of mean and standard deviation of δETLULCC obtained from different model ensembles and different formulations (full ESM ensemble vs. LUMIP ESM ensemble, details in Methods) exhibit consistent spatial patterns (Fig. 1d, f). Altogether, these results unequivocally show substantial uncertainties in ESM-based assessment of the hydrological response to land use and land management practices and call for additional analysis to improve predictions.

Emergent constraints on modelled global LULCC effects over the historical period

To minimize the uncertainty and further improve the robustness of ESM-based assessment on the hydrological effect of LULCC, we propose the use of a hierarchical emergent constraint37. To this aim, we combine the CMIP6 modelled global δETLULCC during the historical period, for which no direct observations exist, with the corresponding global-averaged mean annual ratio between sensible heat and latent heat through transpiration (Tr) (represented by a Tr-specific Bowen ratio (Bts), details in Methods), for which observational estimates are available. Observational estimates of Bts are obtained based on the combination of four observation-based Tr/ET products and two reanalysis climate datasets (details in Methods). Specifically, the observation-based Tr/ET estimates used here are derived from the LAI-based upscaling approach44, and remote sensing models45–47, in combination with 20CR v348 and ERA5-Land climate datasets49. This emergent constraint approach grounds on the tight linear positive relationship across ESMs (r=0.931, p<0.01) between the modelled global δETLULCC and the natural logarithm of modelled global Bts (i.e., ln(Bts)) derived for the observational period (1982-2014) (Fig. 2a). Results derived from one-thousand random sampling further confirm the robustness and significance of the inter-model relationship between δETLULCC and ln(Bts) (Supplementary Fig. 2a, b; details in Supplementary Text 1). Similar to the model-based estimates of δETLULCC (Fig. 1), the global mean ln(Bts) also exhibits substantial difference across ESMs, ranging from 0.352 to 1.360 (Fig. 2a, c). Such emergent relationship shown above reveals that ESMs with high ln(Bts) values generally have a large LULCC-induced ET increase, while ESMs with low ln(Bts) values show a weak increase or even decrease in ET under same LULCC conditions (all CMIP6 ESMs adopt the same land-use forcing dataset LUH241) (Fig. 2a). Within the proposed modelling framework, Bts is transformed to its natural logarithm (1) to better reflect the inherent non-linear relationship between Bts and δETLULCC (Supplementary Figs. 3a and 4; details in Supplementary Text 2), and (2) to meet the linear assumption of emergent constraints50. This emergent relationship provides a robust and effective constraint to the wide range of modelled δETLULCC with the global ln(Bts) values from an eight-member ensemble of observation-based datasets (four Tr/ET products × two climate datasets).

Emergent constraint results show that historical LULCC has led to a decrease in ET with a global average rate of -0.031 mm yr-2 during the period 1982-2014, which is opposite in sign to the original CMIP6 ensemble mean (0.057 mm yr-2) (Fig. 2b). After applying the emergent constraint, the overall transition from positive to negative δETLULCC during the past three decades results from the combination of (1) the positive relationship between the effect of LULCC (characterized as forest-cropland conversion globally) and ln(Bts) and (2) the overall overestimation of Bts in current-version ESMs (modelled ln(Bts): 0.633±0.314; observed ln(Bts): 0.404±0.137) (Fig. 2a, c). The latter bias is largely attributable to the substantial ESM underestimation of the ratio of Tr to total terrestrial ET (i.e., Tr/ET in Supplementary Fig. 5c). The negative value of δETLULCC obtained after constraint is highly consistent with previous findings of observation-based assessments on global ET change4,51, and in accordance with the previously reported decrease in latent heat induced by the conversion of forests to croplands52. The emerging negative signal of δETLULCC over the observational period can be partly attributed to the massive deforestation in the tropics and the induced decrease in Tr and canopy interception that exceed any potential increase in soil evaporation8,53. Furthermore, the constrained distribution shows a standard deviation of 0.081 mm yr-2, which corresponds to a significant (p<0.05, F-test) reduction in model spread by 48.0% compared to the unconstrained (original) CMIP6 ensemble (Fig. 2b).

A comprehensive set of experiments performed on the selection of the model ensemble (Supplementary Table 4) and on the source of observational data including the subsets from the full ensemble of observation-based datasets (Supplementary Table 5) and the eddy covariance measurements at FLUXNET sites54,55 (Supplementary Fig. 6 and Data 1 and 2), further confirms the robustness of our results and the validity of our methodology (details in Supplementary Text 1). Moreover, additional analyses based on the classic emergent constraint56,57 show a consistent sign switch in δETLULCC after constraint and a comparable magnitude of the constrained δETLULCC (-0.048±0.089 mm yr-2), aligning with the estimates derived from the hierarchical emergent constraint described above (Fig. 2 and Supplementary Fig. 7). Such high consistency demonstrates the substantial independence of our results on the type of emergent constraint method applied and further reinforces the robustness of our results (details in Supplementary Text 3).

Emergent constraints at the regional scale

Global-scale assessments illustrated in the previous section clearly show the effectiveness of the emergent constraint approach in correcting the LULCC-related signal derived from the original CMIP6 ESM simulations and in reducing its uncertainty (Fig. 2). Nevertheless, when emergent constraints are applied at the global scale, they are potentially subject to spatial compensatory effects that may originate from areas with diverse land use changes or opposite sensitivities in the land response to environmental drivers. To better elucidate such potential spatial heterogeneity in the response mechanism, we further apply the proposed emergent constraint framework at the regional scale (Fig. 3). To this end, we estimate the observationally constrained hydrological response to historical LULCC over forty-four IPCC AR6 reference regions (Supplementary Fig. 8, source from Iturbide et al.39) (details in Methods). Unconstrained results based on the original simulations from CMIP6 model ensemble generated for the historical period suggest that the LULCC-related signal is much stronger at the regional scale (magnitude ranging from -1.376±1.480 mm yr-2 to 0.886±0.967 mm yr-2 across reference regions, Fig. 3 and Supplementary Data 3) compared to the global mean level (0.057±0.155 mm yr-2, Fig. 1b). According to such simulations, LULCC contributes more than 25% of the historical change in local ET over 49.5% of global land. In hotspot regions experiencing dramatic forest loss like Southeastern South America (SES), Madagascar (MDG), and Southeast Asia (SEA), LULCC plays a dominant role in local ET decrease (relative contribution higher than 50%) (Supplementary Data 3). These results highlight the key relevance of land use and land cover change in long-term variations of the regional water cycle.

We find that applying the emergent constraint always reduces the uncertainty in estimates of ET response (i.e., δETLULCC), although the degree of reduction varies across regions (Fig. 3 and Supplementary Data 3). The inter-model correlation between δETLULCC and ln(Bts) is statistically significant (p<0.05) for twenty-one reference regions corresponding to 50.2% of the global vegetated land (hexagons with the diagonal lines in Fig. 3). These regions exhibit an average reduction in standard deviation after constraint by 21.7% during the historical period, more than 10-fold greater than the analogous estimates obtained for the remaining regions without significant correlations (1.8%) (Supplementary Data 3). Correction direction for estimates of δETLULCC varies spatially, primarily because of the distinct types of land cover conversion occurring across regions during the period 1982-2014 and of their interplay with the region-specific background climate. Over 50.2% of global vegetated land where the emergent relationship holds, areas characterized by an increase in forest area (17.0% of global vegetated land), generally shows a higher δETLULCC after applying the emergent constraint compared to the unconstrained estimate, which is opposite to the global mean results (Fig. 3). By contrast, regions that have experienced a decrease in forest area generally show a constrained estimate of δETLULCC lower than that one obtained from the unconstrained outputs, consistently with the global mean results (Figs. 2 and 3). In this respect, δETLULCC over 13.2% of global land turns from positive to negative after applying the emergent constraint (Fig. 3 and Supplementary Data 3). Unlike the counter-intuitive positive effect in these regions from the original ESM simulations, the constrained results of the negative δETLULCC occurring in regions such as Southern Central America (SCA), Northern South America (NSA) and SES are highly consistent with findings reported in previous literature8,58. Overall, these results demonstrate that original ESM outputs tend to misrepresent the direction and magnitude of historical LULCC effect on regional ET (Fig. 3), and moreover, confirm the effectiveness of emergent constraint approach in correcting the regional signal of LULCC effect derived from the original ESM simulations.

Improved hydrological impact assessment under future afforestation scenarios

Afforestation, one of the most important climate-related LULCC activities, is considered a nature-based and cost-effective strategy for enhancing terrestrial carbon sequestration and therefore is a key element in the portfolio of mitigation actions of many nations15,16. Despite the great potential for carbon neutrality and climate mitigation, afforestation may also pose additional pressures on WA for environment and human society17,18. It is therefore of critical importance to provide robust estimates of the potential hydrological impact that may originate from future afforestation scenarios. Global forest area is projected to increase significantly (p<0.01, Mann-Kendall test) during 2015-2099 under the SSP1-2.6 scenario and to decline under the SSP3-7.0 scenario (Fig. 4a). Accordingly, forest area in SSP1-2.6 is higher than that in SSP3-7.0 over most of global land (Fig. 4b), reflecting the former’s afforestation-oriented sustainability policies versus the latter’s limited climate governance41,59 (details in Methods). To disentangle the potential hydrological impact of afforestation scenarios, we quantify global δETAFF from original model projections as the difference in ET trend generated under two future scenarios driven respectively by SSP1-2.6 and SSP3-7.0 land use and land cover forcing (Eq. (3), details in Methods) and analyze its inter-model relationship with ln(Bts). We find that δETAFF during 2015-2099 exhibits a significantly negative inter-model relationship with modelled historical (1982-2014) ln(Bts) at the global scale (r=-0.878, p<0.01), demonstrating the applicability of the emergent constraint approach also for a future afforestation scenario (Supplementary Fig. 9a). ESMs with low ln(Bts) values generally have a large ET increase induced by future afforestation, while ESMs with high ln(Bts) values show a weak increase or even decrease in ET under the same afforestation scenario. The opposite relationship of ln(Bts) with future δETAFF and historical δETLULCC is because LULCC during the period 1982-2014 is characterized by a progressive forest loss, therefore with an obvious opposite effect on ET with respect to that one originating from the afforestation scenario considered here (Figs. 1a, 2a and 4a and Supplementary Fig. 9a).

The original CMIP6 model ensemble shows an overall positive effect of afforestation on terrestrial ET during 2015-2099, whose global average rate increases from 0.024±0.048 mm yr-2 to 0.051±0.031 mm yr-2 after applying the emergent constraint (Supplementary Fig. 9b) (details in Supplementary Text 4). Similar to results derived for the historical period, more pronounced signals emerge at regional scale under the future afforestation scenario. For example, compared to the global mean level (0.024±0.048 mm yr-2, Supplementary Fig. 9b), unconstrained projections suggest a much higher δETAFF in Central Africa (CAF) whose magnitude reaches 0.420±0.515 mm yr-2, indicating that afforestation could explain 74.1% of future change in local ET (CAF in Fig. 5a and Supplementary Data 4). Under the simulated future afforestation scenario, conversions from croplands to forests are expected to occur in most regions (Supplementary Data 4, exception cases are explained in Methods), generally leading to a higher constrained estimate of δETAFF than the analogous metric derived from the original model ensemble (Fig. 5a). The constrained value of δETAFF differs greatly across afforested regions, ranging from -0.158±0.357 mm yr-2 to 0.447±0.171 mm yr-2 (Fig. 5a). Occurrence of such spatial difference in δETAFF may arise from (1) the varying magnitude of increase in forest cover (Supplementary Data 4), and (2) the divergent hydrological sensitivity to forest gain across regions which is largely controlled by the local climate and landscape conditions25. Previous evidence suggested that ET is generally more sensitive to forest cover change in drier regions, highlighting the importance of accessing deeper soil moisture for ecosystems in water-limited environments60. Consistent with the global mean level (0.051±0.031 mm yr-2, Supplementary Fig. 9b), constrained δETAFF is positive over most (58.8%) part of the afforested land, suggesting the widespread positive effect of afforestation on ET as reported in previous studies8,61. Regions with the most notable afforestation-driven increase in ET manifest in Sub-Saharan Africa and South Asia (Fig. 5a and Supplementary Data 4). On the other hand, the negative value of δETAFF still exists after constraint in some regions such as Northeastern South America (NES) and SEA. This may occur because the shading by increased foliage cover decreases the amount of radiation reaching the ground ultimately leading to a reduction in soil evaporation which offsets the increase in Tr due to afforestation62.

For twenty-two regions (covering 51.1% of global vegetated land) where emergent relationship and correlation between δETAFF and future afforestation impact on precipitation (δPAFF) across ESMs are both statistically significant, local δPAFF is of comparable magnitude to δETAFF after constraint, particularly in tropics and subtropics (Fig. 5a, b and Supplementary Data 4 and 5). For example, δPAFF in North Eastern Africa (NEAF) increases from 0.008±0.659 mm yr-2 to 0.352±0.470 mm yr-2, along with the increase in δETAFF from 0.091±0.354 mm yr-2 to 0.312±0.253 mm yr-2 after applying the emergent constraint (Supplementary Data 4). The ultimate effect of future afforestation on WA (i.e., δWAAFF) in the constrained projection is derived as the difference between constrained δPAFF and δETAFF (Eq. (11), details in Methods). We estimate that sixteen out of twenty-two reference regions which cover 32.9% of global vegetated land, exhibit an increase in δWAAFF after constraint, corresponding to a weakened negative effect of afforestation on WA and even a transition to positive effect in the constrained projection (Fig. 5c and Supplementary Data 4). Such weakening of afforestation-induced negative effect on WA after constraint is particularly evident in Central and South America, Sub-Saharan Africa and South Asia, where constrained δPAFF is substantially higher than unconstrained one (Fig. 5b, c). It can be seen that although agreeing on the negative impact of afforestation on WA, the constrained results show a less pronounced afforestation-driven reduction in WA compared to the original projections both at the global scale (unconstrained δWAAFF: -0.035 mm yr-2; constrained δWAAFF: -0.029 mm yr-2) and over broad regions (Fig. 5c and Supplementary Fig. 9). This suggests the general overestimation of original ESM simulations on the negative impact of afforestation on WA. Furthermore, the existence of regions with higher magnitude of positive δPAFF than that of δETAFF (i.e., constrained δWAAFF > 0) suggests that the increase in ET associated with afforestation could be outweighed by more pronounced moisture recycling in these regions. Such patterns emerge in regions such as NEAF, Western Africa (WAF), South Asia (SAS), and Eastern Australia (EAU) (Fig. 5), and appear broadly consistent with observational evidence documenting the buffering effect of forest-climate feedback on local water depletion18,63.

Altogether these results demonstrate that original ESM simulations overestimate the water security induced by afforestation scenario, particularly in tropics and subtropics (Fig. 5c), and support the value of bias-corrected simulations for a more robust assessment of land-based climate mitigation plans both at the global and regional levels.

Methods
CMIP6 simulations

To explore the effect of historical land use and land cover change (LULCC) on terrestrial water cycle, we use outputs from an ensemble of eleven Earth system models (ESMs) participating in the Land Use Model Intercomparison Project (LUMIP)40 within the framework of the Coupled Model Intercomparison Project Phase 6 (CMIP6) (https://aims2.llnl.gov/search/cmip6/): BCC-CSM2-MR, CanESM5, CESM2, CMCC-ESM2, EC-Earth3-CC, EC-Earth3-Veg, GFDL-ESM4, IPSL-CM6A-LR, MIROC-ES2L, NorESM2-LM, and UKESM1-0-LL (Supplementary Table 1). These models are selected because they provide all the diagnostic variables required for this study under different factorial experiments. The first available realization (r1) for each model is used here to ensure equal weight to each model in the multi-model ensemble mean. For each model, one experiment with all forcings time-varying and one experiment with land use and land cover kept constant (“historical” and “hist-noLu” in the CMIP6 terminology) during the historical period (1982-2014) are analyzed within a factorial simulation framework (Supplementary Table 2). Except land use and land cover, all other forcings such as greenhouse gases, aerosols and nitrogen deposition are identical for these two set of experiments, which are run in fully coupled mode.

To project the hydrological impact of future afforestation (AFF), one SSP3-7.0 scenario experiment with corresponding land cover and land use and one SSP3-7.0 scenario experiment with land cover and land use from SSP1-2.6 (“ssp370” and “ssp370-ssp126Lu” in the CMIP6 terminology) for the period 2015-2099 are analyzed in this study. According to the LUMIP protocol40, SSP1-2.6 represents the afforestation scenario while SSP3-7.0 represents the deforestation scenario, as clearly reflected in Fig. 4. SSP1-2.6 represents the scenario where the combination of high economic growth and technological improvements and climate policies results in substantial reduction in total agricultural land and increase in forest area. Important climate policies from the land use perspective involved in SSP1-2.6 include the increased bioenergy use, and programs of avoided deforestation, reforestation and afforestation59. SSP3-7.0 represents the scenario where land use change is hardly regulated and strong expansion of cropland and pasture to cope with food shortage results in large-scale deforestation41. Therefore, the difference of outputs between “ssp370-ssp126Lu” experiment and “ssp370” experiment can be used to disentangle the impact of future afforestation scenario (details in section “Assessing LULCC effects from model outputs”). We stress that despite the overall afforestation at the global extent, forest cover in SSP1-2.6 is not always higher than that in SSP3-7.0 for all regions, as shown in Fig. 4b and Supplementary Data 4. Previous studies show possible scattered deforested regions (e.g., Eastern Europe and Southern China) under this future afforestation scenario9. Such cases could be related to the anticipated socioeconomic trends and region-specific climate policy considered in future scenarios41,59,84. Our study period ends in 2099 instead of 2100, mainly because of the limitation from time coverage of NorESM2-LM simulations in “ssp370-ssp126Lu” experiment. It is noteworthy that GFDL-ESM4 is excluded from historical analyses, since this model shows an unrealistic decreasing trend in ET over 1982-2014 in the “historical” experiment in contrast with simulations from all other ESMs (Supplementary Fig. 11a) as well as previous reports8,85. Nevertheless, given the good agreement in future trend in ET between GFDL-ESM4 and the other ESMs and to maximize the sample size when establishing the emergent relationship related to hydrological effect of future afforestation, GFDL-ESM4 is included in subsequent analyses (Supplementary Fig. 11b). Additional results derived by leaving GFDL-ESM4 out of our analyses further demonstrate that including/excluding GFDL-ESM4 would have limited impact on our findings (Supplementary Figs. 9a, b and 12). EC-Earth3-CC and EC-Earth3-Veg are employed solely for historical analyses since their future projections under SSP3-7.0 are not available. Details of experiments have been provided for each model in Supplementary Table 3. In addition to the aforementioned experiments run in fully-coupled mode, we use outputs from two sets of experiments run in ocean-uncoupled mode conducted within the Aerosol Chemistry Model Intercomparison Project (AerChemMIP)86 (Supplementary Table 2). For these simulations, sea surface temperature (SST) and sea ice concentration (SIC) are prescribed. Differently from ESMs in LUMIP that are constructed and operated with the full coupling between modules, ocean module and its connection with other modules have been switched off within ESMs in AerChemMIP. In view of this, a pair of experiments in AerChemMIP, i.e., “histSST” and “histSST-noLu”, is considered equivalent to the “historical” and “hist-noLu” in LUMIP86. Two ESMs, namely GISS-E2-1-G and MRI-ESM2-0, participating in these two sets of experiments provide all variables required for this study, and therefore are employed for our analyses (Supplementary Tables 1 and 3).

Land use and land cover forcing in both historical and future scenario periods for CMIP6 (i.e., in all experiments mentioned above) is prescribed by Land-Use Harmonization 2 (LUH2) dataset41. LUH2 is generated based on the History database of the Global Environment (HYDE3.2) and the Food and Agriculture Organization of the United Nations (FAO) wood harvest data, and constrained by Landsat remote sensing observations. As one of the key inputs of LUH2 dataset, the HYDE 3.2 involves the dynamics of cropland and grazing land, the distinction between rain-fed and irrigated cropland, the distinction between pasture and rangeland according to the intensity of use, and the delineation of protected areas87. Furthermore, the inclusion of FAO wood harvest data implies that LUH2 can also reflect the spatio-temporal patterns of wood harvesting.

A set of variables simulated by the ESMs and representative of the mixture of multiple plant functional types (PFTs) within each grid cell are selected for our analyses. These variables include: monthly-scale evapotranspiration (ET), transpiration (Tr), latent heat (LE), sensible heat (H), leaf area index (LAI), precipitation (P), and mean air temperature (T). Full time series of cover fractions of trees, crops and grasses during both historical and future scenario periods are available for five ESMs (CanESM5, CESM2, CMCC-ESM2, IPSL-CM6A-LR, and UKESM1-0-LL), which are also used here to explore the land cover condition and its dynamics globally and regionally (Figs. 1a and 4 and Supplementary Data 3 and 4). We point out that original LUH2 data cluster the global land into five main categories including cropland, grazing land, primary land, secondary land, and urban land41. ESMs are based on different assumptions when translating LUH2 raw information into datasets of tree, crop, and grass covers, thereby leading to the slight model spread in estimates of areas of global forests, croplands, and grasslands (Figs. 1a and 4a and Supplementary Data 3 and 4). In addition, the simplification and aggregation of numerous vegetation types within these ESMs may result in some inconsistency across models and between individual model and observed patterns in terms of cover fraction of PFTs. CMIP6 outputs are first resampled to a common 0.5°×0.5° global grid-cell using the bilinear method of interpolation, and then aggregated at the annual temporal resolution (i.e., accumulated value of monthly ET, Tr, and P in a given year; average value of monthly LE, H, LAI, T, and fractional covers of vegetation types in a given year). Following Forzieri et al.68, grid-cells with multi-year (1982-2014) average LAI < 0.15 m2 m-2 – thus characterized by low vegetation density – are excluded from the analyses. To keep consistent the spatial domain of simulations, the multi-year averaged LAI derived from the multi-model ensemble mean is used as common mask for all ESMs.

Furthermore, output from three additional ESMs (i.e., CESM2, NorESM2-LM, and UKESM1-0-LL) provided at subgrid-scale are exploited for the historical period. These simulations are used to quantify the difference in surface energy partitioning between PFTs (details in section “Identifying the difference in energy partitioning between PFTs”). Subgrid data provide information of LE, H and LAI distributed across four main land use categories within each grid-cell40: (1) primary and secondary lands (i.e., forests, grasslands, and bare ground), (2) croplands, (3) pasture lands, and (4) urban lands. These subgrid-scale categories are generally consistent with the original ones in LUH2 dataset, and primary and secondary lands are not further divided into PFTs like forests and grasslands, as performed at the scale of entire grid-cell.

Observation-based products

The observational ET and Tr are derived from four global gridded datasets, including the LAI-based upscaling44, the Global Land Evaporation and Amsterdam Model (GLEAM v3.8a)45,88 (https://www.gleam.eu/), the Gerrits’ model46 (https://data.4tu.nl/articles/_/12718811/1), and the Penman-Monteith-Leuning Model (PML-V2)47 (https://doi.org/10.6084/m9.figshare.14185739.v6). Further details of these products are provided in Supplementary Table 7. The above-mentioned datasets are selected since they cover the entire reference period 1982-2014 and have been well validated. In addition, the terms of surface energy partitioning, including H and LE are required in this study. Observational LE is derived by combining latent heat of vaporization (2.45 MJ kg-1) and ET estimates derived from the four different observation-based datasets mentioned above. H is obtained from the closure of the energy balance by subtracting LE from the surface net radiation (Rn), under the assumption of no changes in soil heat at the mean annual level34,68. Rn is retrieved for the period 1982-2014 from two different reanalysis datasets, namely 20CR v348 (https://psl.noaa.gov/data/20thC_Rean/) and ERA5-Land49 (https://cds.climate.copernicus.eu/). All observational products have been resampled to the common 0.5°×0.5° spatial resolution and used to derive the long-term annual average of variables mentioned above (i.e., Tr, ET, LE, Rn, and H). Based on such an approach, we derive an eight-member ensemble of H estimates, and an eight-member ensemble of Bowen ratio (H/LE) estimates (Supplementary Fig. 5b), each one derived from a different combination of observation-based ET product and climate reanalysis dataset. Sensible fluxes derived from FLUXCOM Model Tree Ensemble89 (https://www.bgc-jena.mpg.de/geodb/projects/Home.php) are used to verify the consistency of H estimates derived from the closure of the energy balance (Supplementary Table 8). To cover the whole period, FLUXCOM dataset used here is the ensemble mean that encompasses members of different machine learning methods and energy balance corrections within the “RS+METEO” setup.

Eddy covariance measurements of water and energy fluxes derived from 132 FLUXNET sites are also used in this study to further validate the robustness of our results derived from the observation-based gridded products mentioned above (Supplementary Fig. 6a and Data 1 and 2). Site-level datasets of H and LE (ET) are collected directly from the FLUXNET2015 database55 (https://fluxnet.org/data/fluxnet2015-dataset/). Tr dataset54 for these 132 sites is derived from the combination of FLUXNET2015 eddy covariance data with three different Tr estimation methods (https://zenodo.org/records/3978408). The method is based on the Transpiration Estimation Algorithm (TEA) method92. Here we use the mean of Tr estimates obtained from the aforementioned methods as benchmark, and the standard deviation of the three sets of estimates as a measure of uncertainty in observed Tr (Supplementary Data 2).

Spatial domain

In this study, we conduct analyses both at the global level and over reference regions used in the Intergovernmental Panel on Climate Change Sixth Assessment Report (IPCC AR6). Analyses focusing on the global level provide a clear overview of land use change effect on terrestrial water cycle during both historical and future scenario periods (Fig. 2 and Supplementary Fig. 9). This help to resolve the long-standing argument across model-based assessments about the direction and magnitude of changes in ET and terrestrial water availability (WA) in response to LULCC and afforestation actions27 and elucidate the underlying mechanism behind such argument (Fig. 6). Given the fact that global mean results may partially mask emerging patterns due to possible spatial compensatory effects (e.g., opposite LULCC effects amongst different areas of the Planet, Fig. 1c), we perform additional experiments to investigate these patterns at the regional scale and better capture the relevance of LULCC in long-term variations in regional water cycle (Figs. 3 and 5). Original IPCC AR6 reference regions include forty-six land regions and fifteen ocean regions with distinct regional climatic features39. Here, all ocean regions and two land regions located in Antarctica are excluded from our analyses, since we focus exclusively on the vegetated land (Supplementary Fig. 8). To this end, all estimation approaches described below have been applied both at the global scale and over forty-four IPCC AR6 reference regions.

Assessing LULCC effects from model outputs

Following similar approaches reported in literature93–95, the effect of LULCC on annual terrestrial ET (expressed as δETLULCC) during the period 1982-2014 is derived from factorial simulations of multiple CMIP6 experiments by calculating the difference between the trend in annual ET generated in the “historical” experiment and that one derived in the “hist-noLu” experiment:

ΔETLULCC = ΔEThist - ΔEThist-noLu (1)

where EThist and EThist-noLu are the trends in annual ET in the experiment with all forcings time-varying (i.e., “historical”) and in the experiment with land use and land use kept temporally constant (i.e., “hist-noLu”), respectively. The statistical significance of EThist and EThist-noLu is evaluated using the nonparametric Mann-Kendall test (Supplementary Fig. 11), and the significance of the resulting δETLULCC is assessed through t-test (Fig. 1b). The term δETLULCC excludes the effects of all other forcing agents (e.g., greenhouse gases and aerosols) and their associated changes in climate on ET, as these components have been removed from the factorial simulations (Eq. (1)). For both global and regional scale δETLULCC estimates, the ET terms reported in Eq. (1) are obtained by spatial average weighting each grid-cell value based on its land area (Figs. 1b and 2a). Land area for each grid-cell is derived by the product of the grid-cell area and the percentage of the grid-cell occupied by land (i.e., variable “sftlf” in CMIP6 archive). The same methodology is applied consistently for all global-scale and regional-scale aggregated metrics described in the following sections.

For GISS-E2-1-G and MRI-ESM2-0, δETLULCC can be derived from the pair simulations run with ocean-uncoupled mode, by replacing EThist and EThist-noLu in Eq. (1) with trends in annual ET in the “histSST” experiment and in the “histSST-noLu” experiment, respectively:

ΔETLULCC = ΔEThistSST - ΔEThistSST-noLu (2)

Our results demonstrate that including/excluding these two ESMs and their associated estimates does not affect our findings obtained from the full model ensemble (Fig. 1d, f). However, the inclusion of these two additional ESMs help to expand the sample size and to establish a more robust emergent relationship. For this reason, they are adopted in this study.

Moreover, the impact of future afforestation on terrestrial ET (expressed as δETAFF) is quantified through factorial simulations run under the future scenario period (2015-2099) as follows:

ΔETAFF = ΔETssp370-ssp126Lu - ΔETssp370 (3)

where ETssp370-ssp126Lu and ETssp370 are the trends in annual ET in the “ssp370-ssp126Lu” experiment and the “ssp370” experiment, respectively. The only difference between these two sets of experiments is the land use and land cover forcing, which ensures that Eq. (3) could effectively factor out the effects of all other forcings (e.g., greenhouse gases like CO2, and aerosols) and the secondary radiative forcing on terrestrial ET. As mentioned in the previous section, SSP1-2.6 is an afforestation scenario while SSP3-7.0 is a deforestation scenario. The large difference in forest cover and its dynamics between these two scenarios (Fig. 4) allow to capture the afforestation-induced signal of ET variation.

To improve the accessibility, here we use “LULCC” and “AFF” to represent the land cover conversion during the historical period (1982-2014) and the future afforestation scenario period (2015-2099), respectively.

Constraining LULCC effect on evapotranspiration based on observations

The emergent constraint approach has gained prominence in recent years to correct present-day simulations and constrain future projections of unobserved quantities of interest in ESMs or their single components (i.e., land surface models and climate models)50,96,97. This approach help to reduce uncertainty in ESM simulations through the combination of an ensemble of original model simulations with contemporary observations. The core concept relies on the fact that, despite substantial differences across ESMs, relationships between two elements (x and y, respectively) are implicit within ESM solutions of complex equations and associated parameterizations50. Estimates of variable x and variable y may both exhibit the large spread across models, but there may exist strong and robust relationship between these two variables, i.e., y=f(x)+ε, where ε is a relatively small departure from f. If x is a quantity for which observation-based data are readily available, then the relationship f may place a useful constraint on y, provided the observation uncertainty in x is small compared to the range of simulated values. Given the fact that such functional relationship cannot be diagnosed from a single ESM but rather becomes apparent through the analysis of a suitably large ESM ensemble, the approach mentioned above is therefore termed ‘emergent’98. In this respect, when applying the emergent constraint framework, bias correction and uncertainty reduction are not performed for one specific ESM as evaluated in studies focusing on model development31,99,100. By contrast, simulation from each ESM is treated as one sample, and all samples are clustered together to identify an overall functional relationship across ESMs (i.e., f).

The hierarchical emergent constraint proposed by Bowman et al.37 is applied in this study. Compared with the classic emergent constraint, the hierarchical emergent constraint accounts for the uncertainty in the observational datasets and it is, therefore, believed to yield more accurate and precise estimates101. The emergent constraint framework identified in this study links the historical δETLULCC and future δETAFF to the contemporary mean transpiration-specific Bowen ratio (Bts). Compared with the standard Bowen ratio (B=H/LE, the ratio of sensible heat to the latent heat) that represent the overall partition of available energy at the surface, Bts proposed here reflect explicitly the ratio between H and LE used for Tr. As Tr is more sensitive to LULCC than ET, the effect of LULCC on the transpiration-specific Bowen ratio is more direct (Fig. 3b). The analytical formulation of Bts is as follows:

Bts = H / (LE * Tr / ET) (4)

where H and LE are given in W m-2, Tr and ET are in mm yr-1. All variables in Eq. (4) correspond to mean annual values during the period 1982-2014. In line with B mentioned previously, the combination of four observation-based Tr/ET products and two climate reanalysis datasets produces an eight-member ensemble of observed Bts (Supplementary Figs. 5a, b and 13).

Since the natural logarithmic model outperforms in capturing the underlying relationship between Bts and the LULCC effect compared to the alternative linear and quadratic models (Supplementary Fig. 3a), predictor x is set to the natural logarithm of Bts during the period 1982-2014 (i.e., ln(Bts)). The target variable y included in such a framework refers to δETLULCC over 1982-2014 or δETAFF over 2015-2099, for both constraint applications at the global level (Fig. 2 and Supplementary Fig. 9a, b) and at the regional scale (Figs. 3 and 5a). This non-linearity is further validated through random sampling over global grid-cells with one-thousand replications (Supplementary Fig. 4, details in Supplementary Text 2). In addition, sensitivity analyses of the above-mentioned modelling choices related to the predictor variables have been performed and are elaborated at the end of this section.

Following related literature37,101–103, the mean value of constrained predictand (yc) can be estimated by the following equation:

yc = y + kc (x0 - x) (5)

where x and y are the multi-model ensemble mean of x (ln(Bts)) and y (δETLULCC or δETAFF). x0 is the mean of observed predictor (x0), which is derived from the mean of the eight-member ensemble of observation-based ln(Bts) estimates (Fig. 2c). The slope kc is calculated as follows:

kc = r(x,y) * σ(y) / σ(x) * (1 + SNR-1)-1 = k * (σ2(x) / (σ2(x) + σ2(x0))) (6)

where k is the slope of the linear equation relating the target y with the predictor x, and it is estimated by least square method; σ(*) and σ2(*) denote the standard deviation and variance of * (i.e., x or y), respectively. r(x,y) refers to the correlation coefficient between x and y. kc is regarded as the corrected k, by explicitly accounting for the signal-to-noise ratio (SNR) in the observation104. SNR is derived as the ratio between the variance across models (σ2(x)) and the variance across observational datasets (σ2(x0)):

SNR = σ2(x) / σ2(x0) (7)

The constrained standard deviation of predictand (i.e., the uncertainty) is estimated as below37,101:

σ(yc) = σ(y) * sqrt(1 - r2(x,y) * (1 + SNR-1)-1) = σ(y) * sqrt(1 - r2(x,y) * (σ2(x) / (σ2(x) + σ2(x0)))) (8)

where r2(x,y) is the square of the correlation coefficient between x and y.

The unconstrained and constrained ranges of historical δETLULCC and future δETAFF are estimated by assuming a Gaussian distribution37 (Fig. 2b and Supplementary Figs. 6c, 9b and 12b). The relative reduction in standard deviation (RRσ) is calculated as follows:

RRσ = (1 - σ(yc) / σ(y)) * 100% (9)

Clearly, the difference between y and yc is referred to the correction of magnitude in the unconstrained result. Furthermore, the correction magnitude is proportional to the correlation strength between x and y. In view of this, the correlation strength largely determines the correction performance after applying the emergent constraint. Therefore, for analysis at the regional scale, regions where correlation between δETLULCC (δETAFF) and ln(Bts) is non-significant or small, the correction to δETLULCC (δETAFF) by the emergent constraint approach generally has little effect (Figs. 3 and 5a).

We perform an additional set of experiments to test the sensitivity of our constraint results on the log-transformation of Bts. To this aim, we propose a new emergent constraint framework that incorporates Bts in place of ln(Bts) (details in Supplementary Text 2). Outcomes of these analyses suggest that emergent constraint results based on Bts are generally consistent with original ones based on ln(Bts) both at the global and regional scales, but exhibit a weaker reduction in uncertainty (reflected by RRσ) for some regions due to the relatively weaker correlation of ET response with Bts compared to that with ln(Bts) (Figs. 2, 3 and 5 and Supplementary Figs. 2a–d, 9, 14 and 15). The high agreement between the two sets of estimates further demonstrates the robustness of our emergent constraint results. On the other hand, to maximize the performance of emergent constraint in reducing model uncertainty, we adopt the emergent constraint framework built upon ln(Bts) throughout our study.

Moreover, we perform additional analyses replicated by using the classic emergent constraint56,57, and compare the associated estimates against analogous estimates based on the hierarchical emergent constraint described above (details in Supplementary Text 3). Results show a high consistency between the two sets of estimates, further demonstrating the robustness of our constrained results (Fig. 2 and Supplementary Figs. 7 and 9a, b).

Deriving constrained responses of precipitation and water availability

The original model projection of afforestation impact on P (δPAFF) during the period 2015-2099 is derived by applying Eq. (3) to P (in place of ET). We find a strong correlation between afforestation-induced P change and ET change across ESMs at the global scale (r=0.707, p<0.05) (Supplementary Fig. 9c). Additional experiments performed through random sampling over global grid-cells (Supplementary Fig. 2e, f) further confirms the robustness of our finding about the significant inter-model relationship between δPAFF and δETAFF (details in Supplementary Text 4). Similar results have also been reported in previous literature94. The strong inter-model correlation between δPAFF and δETAFF also exists widely over reference regions (Fig. 5b and Supplementary Data 5). Such relationship enables to constrain δPAFF by using observationally constrained estimate of δETAFF (i.e., yc in Eq. (5)) both at the global and regional scales. Hence, following methodological approaches documented in previous studies101,105, the constrained δPAFF is derived as follows:

constrained ΔPAFF = a * constrained ΔETAFF + b (10)

where a denotes the regression coefficient between raw estimates of δPAFF and δETAFF, and b the intercept. As shown in Supplementary Fig. 9c, for analysis at the global mean level focusing on the period 2015-2099, a=1.2358, and b=-0.0408. Here, we use a constrained estimate of δETAFF to derive δPAFF, mainly because P changes show a closer relationship with long-term ET changes compared with ln(Bts). This is likely because some key P determinants, such as afforestation’s regulation to atmospheric circulation and vertical ascent, are better represented by ET change than Bts dynamics106.

Terrestrial water availability (WA), i.e., the net flux of water at the land surface, is here defined as P minus ET9,69,107. Thus, WA change during a given period can be computed as ΔWA=ΔP - ΔET. Correspondingly, the observationally constrained net effect of future (2015-2099) afforestation on WA can be expressed as:

constrained ΔWAAFF = constrained ΔPAFF - constrained ΔETAFF (11)

Identifying the difference in energy partitioning across PFTs

To elucidate the potential mechanisms behind the large spread in the estimates of LULCC effect across ESMs, we perform an additional analysis based on subgrid-scale model outputs, and quantify the difference in ratio of H to LE used for Tr between croplands and forests (i.e., ΔBts forest-crop = Bts crop – Bts forest) during the period 1982-2014. As mentioned previously, forests are not an individual category in LUMIP subgrid data, but are classified as primary and secondary lands along with grasslands and bare ground (details in section “CMIIP6 simulations”). In this study, primary and secondary lands are considered as forests, and therefore, subgrid-scale outputs for category of primary and secondary lands (i.e., Hpsl, LEpsl, and LAIpsl, where psl is the abbreviation of primary and secondary lands) under the “historical” experiment are used to identify the energy partitioning for forests. To minimize the potential effect of this rough PFT classification on our results, only grid-cells dominated by forests (cover fraction of trees ≥ 60%, Forzieri et al.34) are considered in our subgrid-scale analyses (Supplementary Fig. 16). Cover fraction of trees for each grid-cell is directly provided by ESMs, as reported in section “CMIIP6 simulations” and grid-cells with mean annual LAI < 0.15 m2 m-2 are excluded from our subgrid-scale analyses as also performed in previous experiments (Fig. 6c and Supplementary Fig. 10). Since Tr is not available for different PFTs in subgrid-scale data, Tr/ET for forests within each grid-cell at the multi-year (1982-2014) mean level is approximated by Beer’s Law in combination with the corresponding mean annual value of LAIpsl108:

Tr/ETforest = 1 - exp(-kRn * LAIpsl) (12)

where kRn is the radiation extinction coefficient, and is generally assumed as 0.6108. Beer’s Law can provide reasonable and accurate ET partitioning into Tr at the long-time scales (i.e., monthly to annual time scales)109,110, and therefore is applied here for the estimation at the mean annual level. To further validate the performance of Beer’s Law in representing mean annual ET partitioning, Tr/ET estimate based on Beer’s Law is compared against that one directly derived from model outputs at the grid-cell scale (Supplementary Fig. 17). To preserve consistency, we derive the estimate based on Beer’s law by exploiting the LAI value for the mixture of PFTs for each grid-cell, and use model outputs of mixture Tr and ET under the “historical” experiment. Results show the general agreement between these two sets of estimates, demonstrating the efficiency of Beer’s law in deriving multi-year mean level of Tr/ET (Supplementary Fig. 17).

On the basis of Eq. (12), mean annual Bts for forests within each grid-cell during 1982-2014 can be estimated by the following equation:

Bts forest = Hpsl / (LEpsl * (1 - exp(-kRn * LAIpsl))) (13)

Similarly to the estimations of Tr/ET and Bts for forests shown above, Tr/ET and Bts for croplands within each grid-cell can be estimated by replacing variables in Eqs. (12) and (13) with those obtained based on the category of croplands (i.e., LAIcrop, Hcrop, and LEcrop).

Considering the limited numbers (three) of ESMs providing subgrid-scale outputs, the abovementioned analysis applied for the estimation of ΔBts forest-crop is complemented by a more general statistical framework extendible to the full ESM ensemble, thereby further increasing the reliability of our findings about the mechanism responsible for the large inter-model spread in δETLULCC estimate (Fig. 6a, b). Similar to approaches applied in previous studies focusing on biophysical climate impact of LULCC30,70, this framework assumes that, for a given grid-cell, the mean difference in Bts between two periods (ΔBts) under the realistic simulation (i.e., the “historical” experiment) is equal to the sum of Bts variation induced by LULCC (ΔBts LULCC) plus the residual signal (ΔBts res) due to the rising atmosphere CO2 concentration and climate change. ΔBts LULCC reflects the change in energy partitioning associated to LULCC-induced changes in vegetation structure and physiology. For grid-cells experiencing significant conversions of forests to croplands, ΔBts LULCC ≈ ΔBts forest-crop in general, and its magnitude can be estimated as follows:

ΔBts forest-crop ≈ ΔBts LULCC = ΔBts - ΔBts res (14)

where the Δ term represents the mean difference between periods 2000-2014 and 1982-1996. These two 15-year periods are selected since they are temporally independent and represent energy conditions in early 21st century and late 20th century, respectively. Here, grid-cells experiencing significant conversions of forests to croplands are defined as those with both significant (p<0.05, t-test) decrease in cover fraction of trees and significant increase in cover fraction of crops between the two periods (2000-2014 versus 1982-1996). In other words, for each grid-cell included in the land mask, t-test is performed twice: once between the two time-series (2000-2014 and 1982-1996) for cover fraction of trees, and once between the two time-series for cover fraction of crops. ΔBts res is estimated from the signal of Bts variation captured in the adjacent eight grid-cells that surround the target grid-cell. The identified reference grid-cells need to meet three criterions: (1) similar initial land cover and land use condition, defined as the mean differences in cover fractions of trees and crops between reference grid-cells and target grid-cell during 1982-1996 within the range (-10%~10%); (2) stable land cover and land use, namely non-significant (t-test, p>0.05) changes in both cover fractions of trees and crops between the two periods; and (3) limited difference in elevation with respect to the target grid-cell (-100 m~100 m). Elevation information is provided by CMIP6 ESMs, which is independent of time (variable “orog” in the CMIP6 terminology). For these reference grid-cells, we could assume ΔBts LULCC ≈ 0 and consequently ΔBts ≈ ΔBts res. The aforementioned processing enables us to minimize the potential impact of spatial heterogeneity in environmental conditions, such as background vegetation characteristics, climate111 and elevation112 on our results. For cases where two or more adjacent reference grid-cells are identified, we use an inverse distance weighting to estimate ΔBts res from the n (n≥2) reference grid-cells, according to the following equation:

ΔBts res = (Σ i=1 to n ΔBts i / di) / (Σ i=1 to n 1 / di) (15)

where di represents the distance from the target grid-cell to the reference grid-cell i. In such statistical framework, ΔBts forest-crop can be derived as the difference between Bts variation in target grid-cell (i.e., ΔBts in Eq. (14)) and that in adjacent grid-cells (i.e., ΔBts res in Eqs. (14) and (15)). Given that this statistical framework is primarily based on the pairing of neighboring grid-cells that meet specific criterions, we therefore name this statistical framework here “the paired grid-cell method”.

Overall, two different methods are applied here to estimate ΔBts forest-crop, an indicator to quantify the difference in energy partitioning between croplands and forests assumed by various ESMs. The first method relies on the subgrid-cell outputs of three ESMs (i.e., CESM2, NorESM2-LM, UKESM1-0-LL), while the second is based on the paired comparison at the entire grid-cell scale and is applicable to all ESMs. Results based on the aforementioned methods suggest the substantial disagreement between ESMs in terms of the magnitude and even direction of ΔBts forest-crop, and such disagreement largely contributes to the large inter-model spread in estimate of hydrological effect of LULCC (Fig. 6).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Review of \"Emergent constraints on the impacts of land use and land cover change on the terrestrial water cycle\" by Chen et al. The study examines the effects of land use and land cover changes on evapotranspiration in CMIP6 model simulations and reports a global emergent constraint using a set of observational datasets. The constraint differs from the model mean, and the uncertainty also decreases. The authors analyze existing model output and publicly available observational datasets; the novelty is rooted in the type of analysis based on emergent constraints. The manuscript is very well written. However, I do not see a great added value in this study, because (a) some of the emergent constraint relationships do not appear robust, and (b) the reported global signals are very low and thus insignificant in the global water cycle. In my opinion, the drawn conclusions are overstated. I recommend restructuring the study to focus more on hot-spots of land use and land cover change (LULCC) and ET responses, and understand different ET responses on a regional scale, as this is where these water cycle responses matter. Therefore, the manuscript in its current state, with its focus on global constraints, does not fit the outlet of Nature Communications. General Comments: 1. What do we really learn from a global change estimate in ET? Evapotranspiration patterns and changes are diverse, often balance out on continental scale, and LULCC are also quite diverse effects; Please argue why such a global estimate really matters. I think regional constraints or biome-specific constraints are more useful. 2. How does this compare to the overall terrestrial evapotranspiration? Global evapotranspiration is around 600 mm yr-1 (https://www.sciencedirect.com/science/article/abs/pii/S0022169423011666). Thus, an average change of ~0.06 mm yr-2 corresponds to a change of 0.01 percent. Even the highest outlier estimate in your study of 0.38 mm yr-2 translates into a 0.06% change per year; Given all uncertainties in the LULCC reconstructions and the observational datasets you use, these signal are negligible. 3. Robustness of fits: Given the twelve data points in Supp. Fig 2a, it is not convincing to me to fit non-linear functions to such a small point cloud. Why do you expect the relationship to be inherently non-linear across models? Within one model, I could see that, but across the ensemble, I do not see why this would be the case. How does the constraint change if you do not apply the log-transform? 4. Relative change vs. absolute change: You write, e.g \"Results show that the constrained δETLULCC is -0.031 mm yr-2 at the global scale, which is opposite in sign and 155.4% lower than that of the original CMIP6 ensemble mean (0.057 mm yr-2)\"; One can put small numbers in relation to each other so they appear high. This is misleading. Your CMIP6 mean indicates 0.01% yr-2 increase in global evapotranspiration, and your constraint suggests 0.005% yr-1 less global evapotranspiration (assuming 600 mm yr-1). These numbers are negligible. If you compare 2 drops of water to 1 drop of water, you see a 100% increase, but you still have only 2 drops. Specific Comments: L23: Is this not an “enhanced” estimate? You report a negative rate in evapotranspiration. L30: I could not see the increase in \"credibility in projections\" aspect emerging from your analyses. Could you specify? L112: Could this also be related to sensor fusion and the onset of the MODIS time series? L114: A verb is missing; please rephrase to \"can largely be attributed to\". L145: Yes, but perhaps not a global estimate. L156-158: Maybe add a sentence to provide more details on the products, e.g., upscaling products or remote sensing, to improve the flow of the manuscript for the reader. L195: Please use the present tense consistently throughout the manuscript; change \"estimated\" to \"estimate\". L221-226: Again, this is misleading, as both numbers are very small. L299: \"Greatly\" reads a bit strangely here; better to use \"strong\". L382-385: Given the low values, CO2 fertilization and climate change globally have considerably stronger impacts, and the LULCC signal likely disappears completely. Fig 3c: I do not see an inter-model spread emerging here; these fits are not convincing. The linear fit in Fig 3e is a stretch, and there is no clear pattern in the cloud of a few points.", "reviewer_2": "This is an interesting and compelling manuscript. I raise a number of questions regarding both the methodology, as well as some of the discussion. I generally think that some re-writing is necessary in order to make the manuscript more compelling and more accessible. In particular, I find that some of the discussion through lines 311-324 or so provides a lot of insight into why some of the findings turn out they way they do herein. I would recommend bringing some of this information forward, so that the Reader is better prepared for where this manuscript is headed. Further, I find that the points raised in this section and others suggest there is a very strong latitudinal feature to this analysis and that presumably drives some of the results, in particular any decline in ET with LULCC... This suggests that perhaps a more latitudinal analysis, or perhaps a more regionally-based analysis that compares ACROSS regions might be more appropriate...? Or at least that latitude should figure far more prominently as an important independent variable in the analysis...? Further, I am not entirely convinced by the methodology is entirely appropriate. Perhaps it could be a bit better defended...? There is no clear reason why a multi-model mean is likely to be better than some of the independent models integrated into it...? Further, why not simply base the analysis on Bts instead of the modelled representations of T...? Why only use this output to constrain already \"strained\" estimates...? In general, I do think this methodological orientation could be better defended and explained. The findings, however, are uniquely interesting. Thus I look forward to reading the revised version of this submission. Please see the attached submission pdf for additional comments and questions."}, "author_response": "First of all, we would like to thank the two reviewers for their insightful and constructive comments. In the revision of the manuscript, we tried to address all their comments and suggestions to improve the clarity of the methodology and the robustness of the analysis. In the following, we respond to each reviewer’s comment by referring to line numbers of the revised tracked version, when not differently indicated. References cited along with this response letter are reported at the bottom of the document. Reviewer #1: Questions and our responses The study examines the effects of land use and land cover changes on evapotranspiration in CMIP6 model simulations and reports a global emergent constraint using a set of observational datasets. The constraint differs from the model mean, and the uncertainty also decreases. The authors analyze existing model output and publicly available observational datasets; the novelty is rooted in the type of analysis based on emergent constraints. The manuscript is very well written. However, I do not see a great added value in this study, because (a) some of the emergent constraint relationships do not appear robust, and (b) the reported global signals are very low and thus insignificant in the global water cycle. In my opinion, the drawn conclusions are overstated. I recommend restructuring the study to focus more on hot-spots of land use and land cover change (LULCC) and ET responses, and understand different ET responses on a regional scale, as this is where these water cycle responses matter. Therefore, the manuscript in its current state, with its focus on global constraints, does not fit the outlet of Nature Communications. We thank the reviewer for her/his comments. We have substantially revised the manuscript according to the reviewer’s comments by further assessing the validity of our methods and the significance of the emerging signals. In addition, we have conducted new experiments focusing on hotspots of land use and land cover change. Such new analyses are described in response to specific comments below. We sincerely appreciate the important comments raised by the reviewer, particularly the inclusion of regional scale analysis and the increase in confidence of robustness of emergent relationships. Addressing these comments have indeed substantially improved the quality of our analysis and further confirmed the relevance of our work and made the overall narrative more comprehensive. We would like to stress here the novelty of our study, which – in our opinion – the reviewer may have not fully understood. Our study is not a mere numerical exercise based on the emergent constraint framework. On the contrary, by applying such a framework our contribution aims to provide a clearer understanding of the impact of historical land use and land cover change (LULCC) on terrestrial evapotranspiration (ET), both at the global extent and over IPCC AR6 reference regions. With the analogous emergent constraint, we further quantify the potential effect of future afforestation on ET and terrestrial water availability, all critical aspects that remain elusive and not yet comprehensively analyzed in previous studies. In brief, our study provides a robust and comprehensive data-driven assessment of the LULCC impact on terrestrial water cycle during the historical period and under future afforestation scenarios, and identify hotspot regions where water availability is heavily influenced by land use change. We believe our model-data integration approach and associated findings add significant value in the following ways: 1) We investigate the effect of historical LULCC on terrestrial ET both at the global and regional scales, through the integration of multi-source observations and Earth system model (ESM) simulations. In previous studies based on original simulations of ESMs and land surface models (LSMs, i.e., the land component of ESMs), the key LULCC-related signal is largely obscured by the uncertainty and bias of models themselves (Wang and Zeng, 2024). Such model limitations have so far hampered a comprehensive understanding and a robust assessment of the ongoing variations in the water cycle under land surface changes. 2) We project the impact of future afforestation on terrestrial water cycle and identify key hotspot areas where afforestation is expected to play a critical control on local water availability. Although some efforts have been recently devoted into related topics (Zan et al., 2024; Tang et al., 2025), our bias-corrected assessment with the help of emergent constraint framework provides a more robust projection of the potential threat of afforestation to water supply for societies, facilitating the development of effective land-based climate mitigation strategies. 3) We provide novel evidence that the inaccurate assessment of the hydrological effect of large-scale LULCC in CMIP6 ESM simulations is attributable to the model bias in representing the ratio between sensible heat and latent heat through transpiration and its variations across plant functional types (PFTs). Our study is the first one that has revealed the cause of the substantial discrepancy between state-of-the-art ESMs in representing LULCC effects on terrestrial water cycle. This finding highlights the urgent need to address the uncertainty of ESMs driven by the surface energy balance and provides new and important technical reference to the development of enhanced ESMs. Action taken: We have clarified the novelty of our contribution in the revised version of the manuscript (see lines 360-363, 490-494). Furthermore, changes made to address other comments (e.g., regional scale analysis and additional validation tests) have been specifically described in the following responses. General comments: 1. What do we really learn from a global change estimate in ET? Evapotranspiration patterns and changes are diverse, often balance out on continental scale, and LULCC are also quite diverse effects; Please argue why such a global estimate really matters. I think regional constraints or biome-specific constraints are more useful. Response: We perform the emergent constraint analysis at the global scale, primarily to provide a clearer overview of land use change effect on terrestrial water cycle during both historical and future scenario periods. The resulting bias-corrected global scale estimate helps to resolve the long-standing argument across model-based assessments about the direction and magnitude of changes in ET in response to LULCC and afforestation actions (Wang and Zeng, 2024) and elucidate the underlying mechanism. We agree with the reviewer that regional assessments are more useful, especially given that global-scale results are potentially subject to spatial compensatory effects that may originate from areas with diverse LULCC and opposite ET response. In fact, in the originally submitted version, local-scale analyses were already considered but were initially performed in the climate- instead of the spatial-domain and reported at the climate bin level (Fig. 4 in the previous version). Now, considering the reviewer’s suggestion, we have further enhanced our focus on hotspot regions and elaborated the emerging LULCC-related signals at the regional scale in more details. To this end, we have performed a series of new experiments over forty-four reference regions used in the Intergovernmental Panel on Climate Change Sixth Assessment Report (IPCC AR6), to provide specific constraints at regional level. Original IPCC AR6 reference regions include forty-six land regions and fifteen ocean regions, each with consistent regional climatic features (Iturbide et al., 2020). Here, all ocean regions and two land regions located in Antarctica are excluded from our analyses, since we focus exclusively on the vegetated land (Supplementary Fig. 7). To calculate the emergent constraints at the regional scale, we first extract the gridded datasets from observation-based products and model outputs for each reference region based on the map provided by Iturbide et al. (2020). We then estimate the regional averaged effect of historical LULCC on annual ET (δETLULCC) and that of future afforestation on annual ET (δETAFF) for each model (Eqs. (1)-(3)), the regional mean levels of observed and modelled transpiration-specific Bowen ratio (B , Eq. (4)) and their natural logarithm (ln(B )). Consistent with our global-scale aggregated metrics, these regional-scale estimates are also obtained by spatial average weighting each grid-cell value within the region based on its land area. Unconstrained results based on the original simulations from CMIP6 model ensemble generated for the historical period suggest that the LULCC-related signal is much stronger at the regional scale (magnitude ranging from -1.376±1.480 mm yr-2 to 0.886±0.967 mm yr-2 across reference regions, Fig. 4 and Supplementary Table 9) compared to the global mean level (0.057±0.155 mm yr-2, Fig. 1b). This finding is in line well with the reviewer’s intuition. Over 26.0% of global vegetated land, LULCC is believed to dominate historical change in local ET (relative contribution higher than 50%), according to the original CMIP6 model simulations (Supplementary Table 9). Similarly, more pronounced signals emerge at regional scale under the future afforestation scenario. For example, compared to the global mean level (0.024±0.048 mm yr-2, Fig. 3d), unconstrained projections suggest a much higher δETAFF in Central Africa (CAF) whose magnitude reaches 0.420±0.515 mm yr-2, indicating that afforestation could explain 74.1% of future change in local ET (CAF in Fig. 5a and Supplementary Table 10). These results highlight the key relevance of land use and land cover change in long-term variations of the regional water cycle. We find that applying the emergent constraint always reduces the uncertainty in estimates of ET response (i.e., δETLULCC and δETAFF), although the degree of reduction varies across regions (Figs. 4 and 5a and Supplementary Tables 9 and 11). The inter-model correlation between δETLULCC and ln(B ) is statistically (p<0.05) significant for twenty-one reference regions corresponding to 50.2% of the global vegetated land (hexagons with the diagonal lines in Fig. 4). Slightly more widespread patterns are found for the interplay between δETAFF and ln(B ) with twenty-five reference regions – covering 55.5% of the global vegetated land – showing statistically significant (p<0.05) correlations (hexagons with the diagonal lines in Fig. 5a). These regions exhibit an average reduction in standard deviation after constraint by 21.7% and 24.9% during the historical and future scenario periods, respectively, more than 10-fold greater than the analogous estimates obtained for the remaining regions without significant correlations (1.8% and 1.3%) (Supplementary Tables 9 and 11). Correction direction for estimates of δETLULCC varies spatially, primarily because of the distinct types of land cover conversion occurring across regions during the period 1982-2014 and of their interplay with the specific background climate of the region. Over 50.2% of global vegetated land where the emergent relationship holds, areas characterized by an increase in forest area (17.0% of global vegetated land), generally shows a higher δETLULCC after applying the emergent constraint compared to the unconstrained estimate, which is opposite to the global mean results (Fig. 4). By contrast, more regions have experienced a decrease in forest area, and generally show a constrained estimate of δETLULCC lower than that one obtained from the unconstrained outputs, consistently with the global mean results (Figs. 2 and 4). In this respect, δETLULCC over 13.2% of global land turns from positive to negative after applying the emergent constraint (Fig. 4 and Supplementary Table 9). Under the simulated future afforestation scenario, conversions from croplands to forests are expected to occur in most regions (Supplementary Table 10, exception cases are explained in Methods), generally leading to a higher constrained estimate of δETAFF than the analogous metric derived from the original model ensemble (Fig. 5a). The constrained value of δETAFF differs greatly across afforested regions, ranging from -0.158±0.357 mm yr-2 to 0.447±0.171 mm yr-2 (Fig. 5a). Occurrence of such spatial difference in δETAFF may arise from (1) the varying magnitude of increase in forest cover (Supplementary Table 10), and (2) the divergent hydrological sensitivity to forest gain across regions which is largely controlled by the local climate and landscape conditions (Zhang et al., 2017). Previous evidence revealed that ET is generally more sensitive to forest cover change in drier regions, highlighting the importance of the ability to access deeper soil moisture in water-limited environments (Ning et al., 2020). Consistent with the global mean level (0.051±0.031 mm yr-2, Fig. 3d), constrained δETAFF is positive over most (58.8%) part of the afforested land, suggesting the widespread positive effect of afforestation on ET as reported in previous studies (Teuling et al., 2019; Yang et al., 2023). On the other hand, the negative value of constrained δETAFF still exists in some regions after constraint, possibly because the shading by increased foliage cover could decrease the amount of radiation reaching the ground and further lead to the reduction in soil evaporation, which counteracts the increase in Tr due to afforestation (Ukkola et al., 2016). For twenty-two regions (covering 51.1% of global vegetated land) where emergent relationship and correlation between δETAFF and δPAFF across ESMs are both statistically significant, local δPAFF is of comparable magnitude to δETAFF after constraint, particularly in tropics and subtropics (Fig. 5a,b). For example, δPAFF in North Eastern Africa (NEAF) increases from 0.008±0.659 mm yr-2 to 0.352±0.470 mm yr-2, along with the increase in δETAFF from 0.091±0.354 mm yr-2 to 0.312±0.253 mm yr-2 after applying the emergent constraint (Supplementary Table 10). We estimate that sixteen out of twenty-two reference regions which cover 32.9% of global vegetated land, exhibit an increase in δWAAFF after constraint, corresponding to a weakened negative effect of afforestation on WA and even a transition to positive effect in the constrained projection (Fig. 5c and Supplementary Table 10). Such weakening of negative afforestation effect on WA after constraint is particularly evident in Central and South America, Sub-Saharan Africa and South Asia, where constrained δPAFF is substantially higher than unconstrained one (Fig. 5b,c). The existence of regions with higher magnitude of positive δPAFF than that of δETAFF (i.e., constrained δWAAFF > 0) suggests that the increase in ET associated with afforestation could be outweighed by more pronounced moisture recycling in these regions. Such patterns emerge in regions such as NEAF, Western Africa (WAF), South Asia (SAS), and Eastern Australia (EAU) (Fig. 5), and appear broadly consistent with observational evidence documenting the buffering effect of forest-climate feedback on local water depletion (O’Connor et al., 2021; Hoek van Dijke et al., 2022). Overall, these results demonstrate that original ESM outputs tend to misrepresent the direction and magnitude of historical LULCC effect on regional ET (Fig. 4), and more importantly, overestimate the water scarcity induced by afforestation scenario, particularly in tropics and subtropics (Fig. 5c). These results ultimately support the value of bias-corrected simulations for a more robust assessment of land-based climate mitigation plans at the regional level. Action taken: We have added the aforementioned results focusing on IPCC AR6 reference regions in the Results section (“Emergent constraints at the regional scale”) (see lines 261-357), and correspondingly, added the new Figs. 4 and 5 and Supplementary Fig. 7 and Tables 9-11 in the revised version. Meanwhile, we have added related discussion on findings obtained at the regional scale (see lines 389-393, 398-415, 472-484), and described the methodological details related to the experiments performed at the regional scale in the Methods section (“Spatial domain”) (see lines 700-718). To keep concise, Supplementary Tables 9-11 are not presented below. 2. How does this compare to the overall terrestrial evapotranspiration? Global evapotranspiration is around 600 mm yr-1 (https://www.sciencedirect.com/science/article/abs/pii/S0022169423011666). Thus, an average change of ~0.06 mm yr-2 corresponds to a change of 0.01 percent. Even the highest outlier estimate in your study of 0.38 mm yr-2 translates into a 0.06% change per year; Given all uncertainties in the LULCC reconstructions and the observational datasets you use, these signals are negligible. Response: We believe that the reviewer may have misunderstood the analysis we have conducted and the relevance of our results in terms of magnitude. Although the global mean ET is around 600 mm yr-1, as correctly reported by the reviewer, its overall trend during the period 1982-2014 derived from the ensemble of twelve CMIP6 models is 0.532±0.122 mm yr-2. A similar increase rate is estimated from an ensemble of observation-based ET products used in this study (0.412±0.195 mm yr-2), and is further corroborated by previous studies using different datasets (e.g., Yan et al. (2013): 0.46 mm yr-2; Pan et al. (2020): 0.23±0.28 mm yr-2; Ma et al. (2021): 0.31 mm yr-2; Yang et al. (2023): 0.66±0.38 mm yr-2). The value of δETLULCC estimated in our study represents the trend in ET induced solely by land use and land cover change during 1982-2014, and therefore its magnitude should be confronted to the overall trend in ET (0.532±0.122 mm yr-2), rather than the amount of annual ET (600 mm yr-1). The above-mentioned approach is typically used in related literature to quantify the marginal effect of individual factor/process on the long-term trend in target variable (e.g., Zhu et al., 2017; Bastos et al., 2019; Zhang et al., 2022; Liu et al., 2025). The global mean value of δETLULCC derived from the original CMIP6 model ensemble is 0.057±0.155 mm yr-2, corresponding to a relative contribution of LULCC on historical trend in ET of 13.3±33.4% (Fig. 1b and Supplementary Fig. 1). Such overall high relative contribution indicates that LULCC-induced signal should not be neglected, instead, deserves to be analyzed and constrained. We point out that the direct effect of changes in climate factors (e.g., precipitation, temperature and solar radiation) has a comparable contribution (22.5%) to global ET rise over a similar period (Chen et al., 2023). As for elevated atmospheric CO concentration (eCO ), its physiological effect is believed to have a small negative contribution due that the negative effect from eCO -induced reduction in stomatal conductance could overweigh the positive effect from increase in biomass (i.e., CO fertilization) (Lemordant et al., 2018; Liu et al., 2021; Yang et al., 2023). Overall, the LULCC effect on ET trend (the focus of our work) is of comparable magnitude to the effects on trend originating from climate change and eCO . Such comparison with other key drivers of ET clearly highlights the relevance of LULCC effect in the long-term environmental influence on the terrestrial water cycle. The new set of experiments designed to estimated δETLULCC at the regional scale (see our response to your previous comment) further corroborates the considerable effect of LULCC on water fluxes. Indeed, the emerging signal in δETLULCC estimated at the global level is subject to compensatory effects originating from different regions. When δETLULCC is estimated at regional scale, we find that the ET response to historical LULCC and future afforestation can be much stronger compared to analogous estimates obtained at the global level. We estimate that in 78.4% (58.1%) of the global vegetated land, the absolute magnitude of constrained δETLULCC is more than 2 times (5 times) the global mean level (-0.031±0.081 mm yr-2); and meanwhile, in 44.3% (17.6%) of the reference regions, the absolute magnitude of δETAFF is more than 2 times (5 times) the global mean level (0.051±0.031 mm yr-2), as derived from the observationally constrained CMIP6 model simulations (Figs. 2-5). Correspondingly, LULCC shows a higher contribution to ET variation at the regional scale (Supplementary Tables 9 and 10). We therefore believe that the additional analyses conducted at the regional scale suggested by the reviewer further clarify the importance of this phenomenon on the ET trend and the necessity of a better understanding of the LULCC impact on terrestrial water cycle. Action taken: We have added these contents into the Discussion section, to highlight the relevance of LULCC in variations of terrestrial water cycle (see lines 398-415). 3. Robustness of fits: Given the twelve data points in Supp. Fig 2a, it is not convincing to me to fit non-linear functions to such a small point cloud. Why do you expect the relationship to be inherently non-linear across models? Within one model, I could see that, but across the ensemble, I do not see why this would be the case. How does the constraint change if you do not apply the log-transform? Response: We apply the log-transformation to the global averaged value of B mainly because the natural logarithm regression has a better performance in fitting the inter-model relationship between B and δETLULCC over 1982-2014 compared with the linear regression, as reflected by a higher determination coefficient (r2) for the former fitting function (Supplementary Fig. 3a). This regression model shows the p-value lower than 0.01, and therefore its result is statistically significant even if it uses only twelve data points. To further test whether the relationship between B and δETLULCC is inherently non-linear across models, we perform a statistical analysis by applying a random sampling technique with one-thousand replications. We therefore generate one-thousand data sets by randomly extracting samples from 60% of vegetated grid-cells globally without replacement. Each data set includes twelve δETLULCC estimates and twelve associated B values, each pair of which corresponds to one CMIP6 ESM. Subsequently, we fit linear, natural logarithmic, and quadratic regressions to the inter-model relationship between δETLULCC and B for each data set, and use r2 and root mean squared error (RMSE) to identify the regression model that provides the best fit in each data set. We find that mean r2 of natural logarithmic regression fitting reaches 0.878 and is higher than that of linear and quadratic ones (0.851 and 0.876) (Supplementary Fig. 4a). In addition, mean RMSE of natural logarithmic regression fitting is 0.061 mm yr-2, which is lower when compared to the fitting by linear and quadratic regressions (0.067 mm yr-2 and 0.063 mm yr-2) (Supplementary Fig. 4b). As natural logarithmic regression has the best performance in fitting the inter-model relationship between B and δETLULCC, ln(B ) is used in our emergent constraint framework as reference. To further increase the confidence in our approach, we also compute and compare the inter-model correlation between δETLULCC and ln(B ), against that between δETLULCC and B , during one-thousand replications of random sampling similar to that described above. The Pearson’s correlation coefficient between δETLULCC and ln(B ) derived based on these one-thousand data sets reaches 0.936±0.005, higher than that one between δETLULCC and B (0.922±0.006) (Supplementary Fig. 2a). Similarly, δETAFF also exhibits a stronger correlation with ln(B ) across models compared to that one obtained with B (-0.878±0.006 vs. -0.866±0.006) (Supplementary Fig. 2c). Given the proportional relation between the correlation strength and the performance of emergent constraint, these additional analyses further confirm the utility of applying a log-transformation to B in emergent relationship. Moreover, we conduct additional analysis to examine whether applying the log-transformation or not would produce substantially different emergent constraint results. To this end, we replicate the analyses by replacing ln(B ) with B , and compare the constrained estimates derived based on B with the original ones derived based on ln(B ). Results based on B show that the constrained δETLULCC reaches -0.030±0.072 mm yr-2 at the global scale during the period 1982-2014, highly consistent with the estimate derived based on ln(B ) (-0.031±0.081 mm yr-2) (Fig. 2 and Supplementary Fig. 12a,b). As for δETAFF δPAFF and δWAAFF during the period 2015-2099, their global-mean values after the B -based constraint are 0.052 mm yr-2, 0.024 mm yr-2 and -0.028 mm yr-2, respectively, which are also in accordance with those after the ln(B )-based constraint (0.051 mm yr-2, 0.022 mm yr-2 and -0.029 mm yr-2) (Fig. 3 and Supplementary Fig. 12c,d). We further compare the constrained estimates derived based on B against those derived based on ln(B ) over IPCC AR6 reference regions (Supplementary Fig. 13). Results show the general agreement between these two sets of estimates, as reflected by the fact that the absolute differences in their magnitudes are very small in all regions (Figs.4 and 5a and Supplementary Fig. 13). Such high consistency of results derived based on B and those based on ln(B ) indicates that not applying the log-transformation to B has a negligible impact on our results and further demonstrates the robustness of our constraints on hydrological impact assessments. Altogether, these additional experiments demonstrate the robustness of our methodology and clarify the utility of applying the log-transformation to B to better capture the inherent relationship between LULCC effect and surface energy partitioning across models. Action taken: We have added these contents into the new Supplementary Text 2 (“Constrained results derived based on B ” section) and briefly recalled them in the main text (see lines 815-819, 852-863). Correspondingly, we have added the new Supplementary Figs. 2, 4, 12, and 13 in the revised version. Therein, Supplementary Fig. 2 has been presented in our response to the Specific Comment #11. 4. Relative change vs. absolute change: You write, e.g “Results show that the constrained δETLULCC is -0.031 mm yr-2 at the global scale, which is opposite in sign and 155.4% lower than that of the original CMIP6 ensemble mean (0.057 mm yr-2)”; One can put small numbers in relation to each other so they appear high. This is misleading. Your CMIP6 mean indicates 0.01% yr-2 increase in global evapotranspiration, and your constraint suggests 0.005% yr-1 less global evapotranspiration (assuming 600 mm yr-1). These numbers are negligible. If you compare 2 drops of water to 1 drop of water, you see a 100% increase, but you still have only 2 drops. Response and action taken: According to the reviewer’s comment, we have changed the expression of relative term to the absolute term when comparing the difference in magnitude between constrained estimates and unconstrained ones. Nevertheless, we still would like to use relative reduction in standard deviation (RR ) to reflect the reduction in uncertainty after applying the emergent constraint, since RR as well as relative reduction in variance (RRV) is a key indicator for the emergent constraint performance evaluation, and has been widely used in related studies (e.g., Shiogama et al., 2022; Dai et al., 2024; Chai et al., 2025). Furthermore, as illustrated in our response to the previous comment (General Comment #2), the LULCC effect on ET trend is substantial and is of comparable magnitude to the effects on the trend from other environmental factors, such as climate change and elevated CO . On the basis of original simulations from CMIP6 model ensemble, we find that compared to the global mean level (0.057±0.155 mm yr-2), such LULCC-related signal is much stronger at the regional scale, whose magnitude ranges from -1.376±1.480 mm yr-2 to 0.886±0.967 mm yr-2 over IPCC AR6 reference regions (Fig. 4 and Supplementary Table 9). Results obtained at the regional scale further demonstrate the strong interplay between changes in land surface and terrestrial water cycle. Specific comments: 1. L23: Is this not an “enhanced” estimate? You report a negative rate in evapotranspiration. Response and action taken: “Enhanced” here represents the improvement in estimation robustness after bias correction and reduction in uncertainty, instead of the positive rate in evapotranspiration. To avoid misunderstanding, we have changed “enhanced” to “more robust” (see lines 23-24). 2. L30: I could not see the increase in “credibility in projections” aspect emerging from your analyses. Could you specify? Response: Our study shows that the original simulations of ESMs usually report highly uncertain and counterintuitive results about the ET response to historical LULCC. By contrast, when integrating model simulations with the proposed emergent constraint framework, the derived constrained results could well reproduce the historical patterns of LULCC effect on ET, that is a negative δETLULCC at the global scale. How well the model or forecast method reproduces historical data and known patterns is a crucial indicator of credibility (Brekke et al., 2008). In this respect, we stressed in previous version that our model-data integration indeed could enhance the credibility of hydrological projections under LULCC. However, we realize that this may not be a widely accepted statement, so we have avoided using the word “credibility” in the revised version of the manuscript (see line 32). 3. L112: Could this also be related to sensor fusion and the onset of the MODIS time series? Response: Land use and land cover forcing employed in CMIP6 ESMs is derived from the Land-Use Harmonization product (LUH2) (Hurtt et al., 2020). LUH2 product for the historical period (up to 2014) is generated based on the History database of the Global Environment (HYDE 3.2) and the Food and Agriculture Organization of the United Nations (FAO) agricultural land use data, and constrained by Landsat remote sensing observations. Information from MODIS satellite imagery are not exploited in the LUH2 product. Therefore, the contrasting patterns of change in global grassland area before and after 2000 cannot be attributed to the sensor fusion and the onset of the MODIS time series. 4. L114: A verb is missing; please rephrase to “can largely be attributed to”. Response and action taken: We have modified the text according to the reviewer’s comment (see line 113). 5. L145: Yes, but perhaps not a global estimate. Response: Following the reviewer’s suggestion (General Comment #1), we have reshaped the narrative, focusing prominently on the emergent constraint results at the regional scale. We believe that our response and actions (a series of new experiments performed at regional scale) detailed in the answers to the previous comments have addressed this specific issue. 6. L156-158: Maybe add a sentence to provide more details on the products, e.g., upscaling products or remote sensing, to improve the flow of the manuscript for the reader. Response and action taken: Following the reviewer’s suggestion, we have added a sentence to further clarify the observation-based products used here. Details are as follow: “Specifically, observation-based Tr/ET estimates used here are derived from LAI-based upscaling approach44, and remote sensing models45-47, in combination with 20CR v348 and ERA5-Land climate datasets49.” (see lines 159-162). 7. L195: Please use the present tense consistently throughout the manuscript; change “estimated” to “estimate”. Response and action taken: We have carefully checked the manuscript and corrected the text – where appropriate – using the present tense consistently, as suggested by the reviewer. 8. L221-226: Again, this is misleading, as both numbers are very small. Response: We believe that our response and associated actions provided in response to the reviewer’s General Comment #4 have already addressed this specific issue. 9. L299: “Greatly” reads a bit strangely here; better to use “strong”. Response and action taken: We have changed “greatly” to “strong”, as proposed by the reviewer (see line 376). 10. L382-385: Given the low values, CO fertilization and climate change globally have considerably stronger impacts, and the LULCC signal likely disappears completely. Response: We stress that LULCC effect on global ET trend is of comparable magnitude to the effects on trend originating from climate change and elevated atmospheric CO , as shown in our response to the reviewer’s General Comment #2. Furthermore, results of our new analyses performed at the regional scale show that LULCC plays a far larger and even dominant role in influencing terrestrial ET and water availability in some regions across the globe. For example, δETAFF derived from the unconstrained model outputs is projected to be 0.420±0.515 mm yr-2 in Central Africa (CAF), suggesting that afforestation could explain 74.1% of future change of local ET (Fig. 5 and Supplementary Table 10). Even after applying the emergent constraint, the constrained estimate of δETAFF could reach 0.433±0.363 mm yr-2, and corresponding δWAAFF is expected to be -0.195 mm yr-2, which is ~6.7 times of the global mean level (-0.029 mm yr-2) (Fig. 5 and Supplementary Table 10). These results reinforce the importance of accounting for the LULCC to properly understand the driving mechanism behind the hydrological change, particularly at the regional scale. In this respect, we believe that the addition of a series of regional-scale analyses mentioned previously in response to the reviewer’s General Comment #1 has addressed this specific issue. Action taken. We have added the new results derived at the regional scale in the revised version of the manuscript (see lines 261-357, 389-393, 398-415, 472-484), and correspondingly, added the new Figs. 4 and 5 and Supplementary Fig. 7 and Tables 9-11. 11. Fig 3c: I do not see an inter-model spread emerging here; these fits are not convincing. The linear fit in Fig 3e is a stretch, and there is no clear pattern in the cloud of a few points. Response: The inter-model correlation between δETAFF and ln(B ) reaches 0.878, and that one between δETAFF and δPAFF is 0.707 (Fig. 3c,e). Both these two fits are statistically significant, with the p-value lower than 0.01 and 0.05, respectively, even if using only nine points. Therefore, we firmly believe that these analyses are robust, and the model fittings are significant. Nine points are considered here because only nine ESMs provide all required diagnostic variables under the future afforestation scenarios within CMIP6 project (details in Methods). To further enhance the confidence on the robustness of the estimated inter-model relationships between δETAFF and ln(B ) (and δPAFF), we apply a random sampling technique with one-thousand replications, and generate one-thousand data sets by randomly extracting samples from 60% of vegetated grid-cells globally without replacement. The additional validation analyses are similar to those presented in response to the reviewer’s General Comment #3 to test the potential effect of log-transformation to B . We compute the inter-model correlation between ln(B ) and δETAFF, and that between δETAFF and δPAFF for each data set. Results show that Pearson’s correlation coefficient between ln(B ) and δETAFF derived based on these one-thousand data sets reaches -0.878±0.006, with the p-value of 18.6×10-4±3.0×10-4 (Supplementary Figure 2c,d). The Pearson’s correlation coefficient between δETAFF and δPAFF reaches 0.719±0.013, with the p-value of 0.029±0.004 (Supplementary Figure 2e,f). Such consistent and significant inter-model correlations (all p-value lower than 0.05) demonstrate the robustness of the estimated relationships between ln(B ) and δETAFF, and between δETAFF and δPAFF, regardless of the sampling set. More importantly, the performance of the emergent constraint in terms of bias correction and reduction in uncertainty is proportionally dependent on the correlation strength, as shown and discussed in the manuscript (Figs.4 and 5, see lines 288-300, 845-851). The consistent strong correlations revealed here reinforce the efficiency of our emergent constraint framework in improving model-based assessment on hydrological effect of future afforestation. Moreover, to further verify the robustness of our emergent constraint results derived under future afforestation scenario, we conduct nine sets of additional experiments, each of which leaves one specific ESM out of the analysis of the emergent relationship. Results show high inter-model correlations between ln(B ) and δETAFF (all p-value lower than 0.01) and between δETAFF and δPAFF (p-value lower than 0.05 for seven out of nine, and lower than 0.10 for the remaining two) across all additional experiments (Supplementary Table 8). Furthermore, the constrained values of δETAFF derived in all experiments distribute in a narrow range, from 0.047±0.032 mm yr-2 to 0.055±0.033 mm yr-2, and it is generally close to that obtained from the full ESM ensemble (0.051±0.031 mm yr-2) (Fig. 3d). Similar results can also be found in the constrained projections of δPAFF and the ultimate effect on terrestrial water availability (δWAAFF = δPAFF − δETAFF) (Supplementary Table 8). Such high consistency between results from these additional experiments performed on reduced ensemble members and those derived based on the full ESM ensemble, further demonstrates that the strong inter-model correlations consistently exist and our emergent constraints are robust and not conditioned by any specific ESM. Action taken: We have described the aforementioned experiments and associated results in the new Supplementary Text 1 (“Test on the robustness of emergent constraint approach” section) and briefly recalled them in the main text (see lines 228-232, 869-872) with references to the new Supplementary Fig. 2c-f and Table 8. We believe that these new analyses could help to (1) improve the robustness of our results about the strong inter-model correlation between ln(B ) and δETAFF, and that between δETAFF and δPAFF, and (2) support the validity of the use of emergent constraint framework to improve the model projection of the hydrological response to future afforestation. Moreover, to further improve the robustness of our emergent constraint results derived for the historical period, we have also applied the random sampling technique to the inter-model relationship between ln(B ) and δETLULCC (Supplementary Fig. 2a,b), and briefly recalled the associated results in the main text (see lines 165-168)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Review of revised manuscript \"Emergent constraints on the impacts of land use and land cover change on the terrestrial water cycle\" by Chen et al. I thank the authors for the in‑depth and comprehensive revisions, and for addressing all my major and minor concerns. Especially, the inclusion of the regional assessment in the form of to‑the‑point visuals makes this study considerably stronger and clearer. Overall, I still think that the global constraint the authors place at the centre of the manuscript is a) IMO not correctly interpreted and b) does not add much to our understanding of the global effects on the hydrological cycle induced by LULCC. In their rebuttal, the authors claim that I did not understand the analysis; however, they repeat the same rationale as in the manuscript: To demonstrate that the effect is non‑negligible, they present the numbers as a relative change rather than an absolute change. As I wrote in my first review report, one can put small numbers in relation to each other so they appear higher; this is misleading. Your global constraint is 0.031 ± 0.081 mm yr⁻²; again, a) this is a very small number in the global hydrological cycle and it does not become more impactful when expressed as a relative contribution to the overall ET trend, and b) the interpretation that your constraint is qualitatively different from the CMIP6 mean is a huge over‑interpretation. Yes, it switches the sign in the mean, but your uncertainty estimate is considerably larger than the effect size. Is the sign switch still present if you use a bootstrap approach when deriving the emergent relationship in the model ensemble (randomly sample models, conduct your EC analysis, and test whether the sign switch is robust, e.g., EC in Winkler et al., 2019)? Returning to the small effect size, looking at your PDFs in Figure 2, my interpretation is that neither the unconstrained nor the constrained estimate is robustly non‑zero (if you leave out GISS‑E2‑1‑G in your multi‑model mean, the unconstrained mean estimate is likely spot-on zero). Thus, to back your central claim—“The constraint reverses the sign of the original ESM estimates obtained over the observational period (0.057 ± 0.155 mm yr⁻², mean ± s.d.) and narrows the inter‑model spread, ultimately leading to a more robust global estimate of the LULCC‑driven ET change of –0.031 ± 0.081 mm yr⁻²”—you must a) test that the sign switch is robust, and b) show that your effect is statistically different from zero. Bottom line, I stick to my recommendation to rearrange the manuscript to focus on your new analysis on hot‑spots of land‑use and land‑cover change (LULCC) and ET responses, and give the different regional ET responses the spotlight, as this is where these water‑cycle responses matter. Therefore, the manuscript should be revised accordingly before being considered for publication in Nature Communications.", "reviewer_2": "I appreciate the effort the authors have put into revising the manuscript. I find it is much improved. I likewise appreciate the effort that the Authors have put into looking at and assessing such broad range of datasets and resources. The attempts at producing more realistic estimates of ET, in particular, are both compelling and fascinating. I do still have many questions. But I think. it is appropriate to move forward to eventual publication of this submission. I have still included some comments here and there in the pdf submission. Mostly minor things. But there are still some questions that continue to nag at me... I do think, for example, that some of the discussion about the relative variation in estimations of the Bowen Ratio (Bts), for example, that are available in the Supplementary information, could be brought forward a bit. And I am still a bit concerned about our ability to estimate ET based on tree type, since there is SO MUCH variation in the real world empirical output that is not any where near adequately captured in the model simplifications... That said, I do think the Authors do an exceptional job of putting this information together and assessing it in meaningful and important ways... Thus, though I would ask the Authors to consider my remaining comments, I do think this submission is more or less ready for publication."}, "author_response": "First of all, we would like to thank the two reviewers for their insightful and constructive comments. We tried to address in this revision their remaining comments and suggestions to improve the robustness of the analysis and the quality of the text. In the following, we respond to each reviewer’s comment by referring to line numbers of the revised tracked version, when not differently indicated. References cited along with this response letter are reported at the bottom of the document. Reviewer #1: Questions and our responses I thank the authors for the in‑depth and comprehensive revisions, and for addressing all my major and minor concerns. Especially, the inclusion of the regional assessment in the form of to‑the‑point visuals makes this study considerably stronger and clearer. Overall, I still think that the global constraint the authors place at the centre of the manuscript is a) IMO not correctly interpreted and b) does not add much to our understanding of the global effects on the hydrological cycle induced by LULCC. In their rebuttal, the authors claim that I did not understand the analysis; however, they repeat the same rationale as in the manuscript: To demonstrate that the effect is non‑negligible, they present the numbers as a relative change rather than an absolute change. As I wrote in my first review report, one can put small numbers in relation to each other so they appear higher; this is misleading. Your global constraint is 0.031 ± 0.081 mm yr⁻²; again, a) this is a very small number in the global hydrological cycle and it does not become more impactful when expressed as a relative contribution to the overall ET trend, and b) the interpretation that your constraint is qualitatively different from the CMIP6 mean is a huge over‑interpretation. Yes, it switches the sign in the mean, but your uncertainty estimate is considerably larger than the effect size. Is the sign switch still present if you use a bootstrap approach when deriving the emergent relationship in the model ensemble (randomly sample models, conduct your EC analysis, and test whether the sign switch is robust, e.g., EC in Winkler et al., 2019)? Returning to the small effect size, looking at your PDFs in Figure 2, my interpretation is that neither the unconstrained nor the constrained estimate is robustly non‑zero (if you leave out GISS‑E2‑1‑G in your multi‑model mean, the unconstrained mean estimate is likely spot-on zero). Thus, to back your central claim — “The constraint reverses the sign of the original ESM estimates obtained over the observational period (0.057 ± 0.155 mm yr⁻², mean ± s.d.) and narrows the inter‑model spread, ultimately leading to a more robust global estimate of the LULCC‑driven ET change of –0.031 ± 0.081 mm yr⁻²” — you must a) test that the sign switch is robust, and b) show that your effect is statistically different from zero. Bottom line, I stick to my recommendation to rearrange the manuscript to focus on your new analysis on hot‑spots of land‑use and land‑cover change (LULCC) and ET responses, and give the different regional ET responses the spotlight, as this is where these water‑cycle responses matter. Therefore, the manuscript should be revised accordingly before being considered for publication in Nature Communications. We thank the reviewer for her/his comments. We would like to stress that the emergent constraint approach applied in Winkler et al. (2019), which is typically referred as the classic emergent constraint originally proposed by Cox et al. (2013), presents some key differences from the hierarchical emergent constraint approach applied in our study. This classic emergent constraint approach involves convolving the prediction error implied by the regression fit of the scatter plot to the inter-model relationship between x and y, with the observational uncertainty in x, to produce a constrained distribution for the y-axis variable (Cox et al., 2018). The prediction error in classic emergent constraint framework is reflected by the probability contours around the best-fit straight line, which are derived from bootstrap resampling of the model pairs (i.e., (x, y), where i indexes each model). In other words, the classic emergent constraint tempts to simply use the regression line itself to project the observational estimate of y, and the prediction error for fit can only be approximated empirically via bootstrap. In contrast, the hierarchical emergent constraint relates model-simulated x and y, and observations of x through conditional-probability distributions derived from Bayes’ theorem, with all sources of uncertainty explicitly represented as probability distributions (Bowman et al., 2018). The prediction error therefore can be analytically derived from the posterior distribution within this framework, avoiding the need for approximation by bootstrap resampling. In this respect, there is no necessity to use a bootstrap approach when deriving the emergent relationship within the hierarchical emergent constraint framework. To further assess the validity of our approach and better compare its robustness with methods suggested by Rev. 1, we have performed a set of additional experiments based on the classic emergent constraint approach (Cox et al., 2013; Winkler et al., 2019). The obtained estimates are then compared against analogous values based on the hierarchical emergent constraint approach currently shown in the text. Such assessment allows to build confidence on the robustness of our results about the sign switch of historical land use and land cover changes effect on global evapotranspiration (δETLULCC) after constraint. Following the approach reported in Winkler et al. (2019), we have applied bootstrapping to estimate the 68% confidence of the emergent linear relationship between δETLULCC and natural logarithm of transpiration-specific Bowen ratio (ln(B )) across CMIP6 ESMs. To this aim, we have randomly resampled the data with replacement keeping the resample size equal to the size of the original sample (i.e., twelve). We then have computed the least-squares linear best fit for the resampled data, and repeated this procedure one-thousand times, and further derived the 68% confidence contours of equal probability based on the set of one-thousand random regression lines. By combining this 68% confidence contours estimated by bootstrapping with ln(B ) observations and associated uncertainty, we have finally derived the probability density function of the classic emergent constraint on δETLULCC. Results based on the classic emergent constraint show that historical LULCC has led to a decrease in ET with a global average rate of -0.048 mm yr-2 during the period 1982-2014, which is generally consistent with the current estimate derived based on the hierarchical emergent constraint (-0.031 mm yr-2) and shown in our manuscript. Moreover, both sets of results consistently show a reversal in the sign of original ESM estimates (0.057 mm yr-2). The high consistency emerging among the two sets of results demonstrates the substantial independence of our results on the selection of emergent constraint type. In additional to the new experiment described above, the previous sets of sensitivity experiments performed on the selection of the model ensemble (by iteratively leaving one specific model out of the emergent relationship) and on the source of observational data also show a consistently sign switch of δETLULCC after applying the emergent constraint. We believe that all these additional experiments together demonstrate the robustness of our constrained results about the negative value of global δETLULCC. Furthermore, the inter-model spread after applying the classic emergent constraint ranges between -0.136 and 0.042 mm yr-2 expressed by the 68% interval assuming a Gaussian distribution and corresponds to a standard deviation of 0.089 mm yr-2. The wider spread obtained when replacing the hierarchical emergent constraint by the classic one (0.081 mm yr-2 vs. 0.089 mm yr-2) is consistent with findings reported in previous literature (Bowman et al., 2018). Such better performance of hierarchical emergent constraint arises from its analytical Bayesian conditioning and the absence of bootstrap-induced resampling noise. To minimize the uncertainty and further derive a more robust ESM-based assessment on the hydrological effect of LULCC, we therefore adopt the hierarchical emergent constraint instead of the classic one. To further illustrate the efficacy of the hierarchical emergent constraint in reducing inter-model spread, we have additionally performed the F-test in the revised version to evaluate whether the post-constraint variance in δETLULCC estimate is significantly smaller than that of the original CMIP6 model ensemble. Test results show a p-value of 0.020, demonstrating that applying hierarchical emergent constraint could significantly reduce the uncertainty in δETLULCC estimate originated from CMIP6 ESMs. We point out that the LULCC effect at the global mean level is not statistically different from zero (one-sample t-test, p>0.05), even after applying the emergent constraint. This situation occurs largely due to the fact that the contribution of LULCC on ET is relatively small at the global mean level, which is subject to spatial compensatory effects that my originate from areas with diverse LULCC and opposite ET response. In light of these considerations, we agree with the reviewer that the work should prominently focus on a regional-scale assessment. To this end, in the new round of revision, we have restructured the manuscript to focus more on the hotspots of LULCC and ET response. A detailed explanation has been given below. Action taken: We have added the aforementioned results based on the classic emergent constraint in Supplementary Text 3 (“Results obtained within the classic emergent constraint framework” section), and briefly recalled them in the main text (see lines 205-212 and 867-872). To further improve the robustness of our emergent constraint results derived for the future scenario period, we have also used the classic emergent constraint on the projections of afforestation impact on ET (δETAFF) over 2015-2099. Correspondingly, we have added the new Supplementary Fig. 7 in the revised version. We have also added the F-test results in the main text (see lines 197-199). Furthermore, to make the manuscript more concise and readable and better focused on the core topic, we have streamlined the global-scale analyses in the main text and particularly moved most of the analyses related to the global impact of future afforestation into Supplementary Text 4 (“Details in constrained projections of global impact of future afforestation” section). Meanwhile, we have added more details about the results derived at the regional scale into the main text. The new contents have been structured in two Results sections to present the hotspot regions of LULCC and hydrological responses and to clarify the differences across regions in hydrological responses (“Emergent constraints at the regional scale” section and “Improved hydrological impact assessment under future afforestation scenarios” section)."}, {"round": 3, "reviewer_comments": {"reviewer_1": "I thank the authors for the detailed additional analyses and the restructuring/revisions of their manuscript. I appreciate that the authors took my concerns about presenting a global constraint that is not sig. different from zero as the main finding into account and agreed that the work should prominently focus on a regional-scale assessment. Thank you for the constructive discussion here. After three rounds of review and given the rebuttals to Reviewer 2 also, this work is ready for publication in Nature Communications.", "reviewer_2": "I am impressed at the lengths the Authors have gone to, to respond to all the very detailed comments of both Reviewers (myself included). I think this is an impressive piece of work and look forward to seeing it in print. I do find that I am repeatedly inspired, in reading this text, to think about further issues that could be addressed and discussed. But this is perhaps the best reflection of the fact that this is an important piece of work and should further inspire others to take these ideas and points still further in future endeavours. Thus, rather than go into additional detailed comments, I would like to encourage publication at this time and hope that others will find this publication as thought provoking as I have. I did notice here and there that some very minor editing of the manuscript is still possible... The attached provides some indication of possible improvements. Please note that I have made no attempt to be exhaustive. Thus more editing is likely possible... However, these editing issues are minor and not worthy of further discussion."}, "author_response": "We thank the reviewer for her/his overall positive comments. Specific comments and associated responses have been summarized below, according to the annotations shown in the pdf document uploaded. Specific comments: 1. Line 18: Replace “effect” with “effects” Response and action taken: To meet the journal’s requirement for Abstract section, we have modified this sentence completely, and the word “effect” is no longer used. 2. Line 20: Replace “effect” with “effects” Response and action taken: To comply with the journal’s word limit for Abstract section, we have deleted this sentence. 3. Line 46: If anthropogenic change includes re/afforestation, there can also be positive effects... This may suggest only negative effects are the likely outcome... Response and action taken: We agree with the reviewer that LULCC-driven processes can generate not only negative but also positive impacts on socio-ecological systems and services. To avoid this one-sided interpretation, we have modified this sentence to present a more balanced perspective. Details are as follows: “Such LULCC-driven processes may lead to shifts in soil moisture and total terrestrial water availability (WA, defined as P minus ET)5,9 and trigger cascading detrimental or beneficial impacts on socio-ecological systems and services crucial for the human well-being, such as biodiversity conservation11, carbon uptake and storage12,13, food security and timber production14.” (see lines 34-39). 4. Line 51: New sentence Response and action taken: According to the reviewer’s suggestion, we have split the original sentence into two separate sentences (see lines 39-43). 5. Line 122: Replace “effect” with “effects” and replace “is” with “are” Response and action taken: We have modified the text according to the reviewer’s suggestion (see line 112). 6. Line 126: Replace “underlies” with “underscores” Response and action taken: We have modified the text according to the reviewer’s suggestion (see line 116). 7. Line 150: Replace “effect” with “effects” Response and action taken: We have modified the text according to the reviewer’s suggestion (see line 140). 8. Line 295: Add “the” before “original CMIP6 model ensemble” Response and action taken: We have added “the” in this sentence (see line 285). 9. Line 429: Replace “that” with “the” Response and action taken: We have modified the text according to the reviewer’s suggestion (see line 419)."}]

Decision letter