Introduction

Forests provide a plethora of ecosystem services and are indispensable to the anthroposphere and the biosphere1,2,3,4,5,6. However, forests consume water for physiological functioning and growth, thus they might potentially deplete Earth’s water resources7,8,9,10. Conceptually, forest communities consume water resources via three routes. First, as living systems, forests absorb and assimilate water through transpiration11,12. Transpiration differs greatly among various plant species that compose the forests and among surrounding environments13. Second, rainfall can be intercepted partially by forest canopy and evaporated directly into the atmosphere without entering the forest hydrological system14. Third, as late participants of Earth’s terrestrial system, forests enhance the complexity of Earth’s surface system and thus affect soil evaporation15,16. The amount of water consumption due to increased forest complexity varies greatly with climatic conditions and watershed characteristics17,18,19.

However, compared to other vegetation covers or bare lands, forest communities are also equipped with advanced strategies to reduce water consumption. For example, forests can create a relatively closed physical environment under the canopy, in which air temperature and wind speed are significantly lower20,21 due to mutual conversion between water and vapor to expend energy inner the closed environment, and thus may reduce the evapotranspiration of the whole ecosystem. Such biophysical effects22,23, possibly together with biogeochemical effects24, of forest communities decrease the conductance from soil to canopy, reducing water loss by concentrating evapotranspiration to the outmost layer of the ecosystem25,26. Biophysical effects, determined by the structure of the ecosystem, are expected to be more apparent in forest ecosystems with a relatively higher complexity. Ecosystem complexity increases as succession progresses27, during which hydrological processes are also in continuous coevolution with other ecological processes—if free from external interferences. Besides, coevolution mechanisms progressively intensify water internal circulation toward inner forest ecosystems, promote the substitution of species with distinct water demands28, and thus maintain the balance between water supply and demand, or even reduce water use. Nevertheless, co-evolving hydrological and other functional processes may be asynchronous, delaying the emergence of water saving effect. The various processes that may negatively or positively affect water resources are restricted by energy and water supplies both inside and outside the forest ecosystems, which creates an uncertainty in the forest–water relationship29,30.

The forest–water relationship may also be complicated by inappropriately using the term “forest”. In the hydrological literature, the term “forest” represents the most extensive and conspicuous vegetation type on Earth, which comprises all the vegetation communities with trees as the constructive species regardless the number of trees on the land. The term “forest” may also refer to forests of different origins (planted forests or natural forests), various ages, and naturally or artificially promoted growth. These highly variable uses of the term “forest” combined with diversified climatic circumstances aggravate the uncertainty of the forest–water relation. Such an uncertainty also contributes to the debates of whether forests increase or decrease water yield, which dates back to the middle of the nineteenth century or earlier31,32, and continue to these days7,17,32,33,34. At present, the increasing global greening during the past decades35 once again induces concerns about declining water resources36,37.

This study aims to clarify the uncertainties of the forest–water relationship and to unearth the potential for more forests but with less water consumption. For this purpose, the responses of runoff coefficients (R/P, R-runoff, P-precipitation) to the forest changes of two causes (artificial interference events and forest natural growth) were investigated. Both artificial interference and natural growth could lead to forest changes, in which the former one is abrupt, ephemeral, and intermittent compared to the later one. Therefore, we quantify the impact of artificial interferences (i.e., disturbances) by comparing the R/Ps before and after disturbance events, while also evaluating the effect of forest natural growth by examining the temporal change of R/Ps. Thus, we compare the differences in the instantaneous responses of annual R/Ps to artificial interference events and identify the long-term changes of annual R/Ps under natural growth. We expect to expound the forest–water relationship to inform worldwide forest restoration projects and aid in achieving the UN Sustainable Development Goals (SDGs), in particular the SDG 15 for sustainable use of terrestrial ecosystems and sustainable forest management.

Results

Responses of runoff coefficient to treatment events

We first investigated if the ΔR/Ps induced by treatment events were significantly different from zero (Fig. 1). DECREASE treatment significantly enhanced the R/Ps (p < 0.05) in all combination categories except in the natural forests located in P/PET ≥ 1 region (Fig. 1a–c). In comparison, INCREASE treatment significantly reduced (p < 0.05) the R/Ps in three of the eight combinations (Fig. 1d–f). The three combinations of forests with the R/Ps significantly decreased are the categories of all planted forests, all forests in P/PET < 1 regions, and the planted forests in P/PET < 1 regions (Fig. 1d–f and Supplementary Table 1).

Fig. 1: Differences in the response of runoff coefficient change (ΔR/P) to treatments (DECREASE and INCREASE) among the eight combinations of forest and climatic region.
figure 1

a–c DECREASE treatment; Panels d–f INCREASE treatment. The number close to each column indicates the sample size (number of watersheds). Error bars indicate the standard error of mean. * indicates significant difference from zero (p < 0.05, one-sample t-test, see Supplementary Table 1 for detail). # indicates significant difference found between two columns connected by solid lines (p < 0.01, two-sample t-test, see Supplementary Table 2 for detail), while non-significant comparisons are linked by broken lines. DECREASE: felling or thinning treatment that directly reduce tree numbers; INCREASE: afforestation or reforestation treatment that favors increasing the number of trees.

Full size image

We then investigated ΔR/Ps between each pair of contrasting combinations (combinations with contrast origins or/and climate regions, Fig. 1). The results showed that both DECREASE (Fig. 1a–c and Supplementary Table 2) and INCREASE (Fig. 1d–f and Supplementary Table 2) treatments resulted into significant differences in ΔR/Ps between P/PET < 1 and P/PET ≥ 1 regions (Fig. 1a, d, p < 0.01), and between planted and natural forest watersheds (Fig. 1b, e, p < 0.01). Particularly, the ΔR/Ps were significantly higher in planted forests than in natural forests in P/PET < 1 regions for both INCREASE and DECREASE treatments (Fig. 1c, f, p < 0.01), but no significant difference was found between planted and natural forests in P/PET ≥ 1 regions in both INCREASE and DECREASE treatments (Fig. 1c, f, p > 0.05). Additionally, for planted forests, the ΔR/Ps were significantly higher for P/PET < 1 regions than in P/PET ≥ 1 regions under both DECREASE and INCREASE treatments (Fig. 1c, f, p < 0.01). In comparison, for natural forests, the ΔR/Ps were significantly higher for P/PET < 1 regions than in P/PET ≥ 1 regions under DECREASE treatment (Fig. 1c, p < 0.01), but no significant difference was found in INCREASE treatments (Fig. 1f, p > 0.05).

In addition to the impacts of forest origin and climate region, we also examined the impacts of pre-treatment forest coverages and forest ages on ΔR/P. Results revealed that ΔR/P decreased as pre-treatment forest coverages and forest ages increase under both DECREASE and INCREASE treatments (Fig. 2, p > 0.05).

Fig. 2: Absolute values of net change of runoff coefficient (│ΔR/P│) influenced by forest coverage and age before treatment events.
figure 2

a│ΔR/P│ with forest coverage in DECREASE treatment; b│ΔR/P│ with forest coverage in INCREASE treatment; c│ΔR/P│ with forest age in DECREASE treatment; d│ΔR/P│ with forest age in INCREASE treatment. Shaded area represents 95% confidence interval. DECREASE: felling and thinning treatments directly reduce tree numbers; INCREASE: afforestation and reforestation treatments that favors increasing the number of trees.

Full size image

Responses of runoff coefficient to forest natural growth

We examined the 278 watersheds under natural growth/recovery (i.e., free from any artificial interferences since the first measurement or after INCREASE/DECREASE treatment), which covers a span of 25 years. Results showed that, by pooling the 278 watersheds, the change of annual R/Ps was non-significant during the natural growth process (Fig. 3a, p > 0.05). Similarly, no significant trend was detected in the annual R/Ps for the treatment categories (i.e., DECREASE and INCREASE) and forest origin categories (Fig. 3b, d). However, it should be noted that mild decreases were detected for most temporal trend analyses for the first 10 data-years, which was followed by rapid increases in R/P (Fig. 3).

Fig. 3: Temporal trends of forest runoff coefficient (R/P) under natural growth.
figure 3

a All the 278 watersheds with 2–25 years of measurement records; b the watershed of DECREASE and INCREASE treatments; c P/PET < 1 and P/PET ≥ 1 regions; d planted and natural forests; e forest ages (pre-treatment); f forest coverage (pre-treatment). Error bars indicate the standard error of mean. Shaded area represents 95% confidence interval.

Full size image

Notably, the annual R/Ps significantly increased in both P/PET < 1 region (Fig. 3c, p < 0.05) and the region with pre-treatment forest coverage below 70% (Fig. 3f, p < 0.001). Likewise, when we divided all the 278 watersheds into two groups based on the pre-treatment forest age using the threshold of 20 years (which divided the samples into two groups with sample sizes almost equal), the annual R/Ps also significantly increased with the pre-treatment forest age (Fig. 3e, p < 0.001). Note that we attempted but failed to identify the age threshold of the shift in the R/Ps-age relationship due to lack of sufficient watersheds with old forest ages.

Discussion

Forest complexity increases hydrological resistance to disturbances

In general, natural forests, old forests, forests with high coverage, and forests located in low aridity regions (P/PET ≥ 1) are characterized by higher ecosystem complexity than planted forests, young forests, forests with low coverage, and forests located in arid regions (P/PET < 1)7,17,38. An ecosystem with a higher complexity is anticipated to be more resilient under the influence of external stresses and to maintain relatively stable ecosystem processes (including water cycles)39. For example, a decrease in transpiration and interception due to tree removal by DECREASE treatment can be partially compensated by the additional evapotranspiration from other routes, including more open forestland, released trees formerly suppressed, reduction of wind barriers, and increases in canopy conductance, etc40,41. In comparison, the increases in transpiration and interception from more trees introduced by INCREASE treatment can also be mitigated by the reduction of evapotranspiration resulted from more closed forestland, more suppressed trees, rising of wind barriers, and decreases in canopy conductance. These two distinct processes can be deemed as compensation and mitigation effects, respectively. Forests with higher ecosystem complexity possess stronger compensation or mitigation capabilities. The observed ΔR/Ps in PWE are, in fact, the results adjusted with compensation or mitigation effects, and thus the ΔR/Ps in PWE should be deemed as an “apparent effect” of forest changes on water yield.

The results presented in Figs. 1 and 2 confirmed the above perspective. For example, the absolute values of ΔR/P due to both DECREASE and INCREASE treatments decreased with the increases of forest coverage and age (Fig. 2), i.e., decreased with the increases of ecosystem complexity. Moreover, when comparing two contrasting combinations, we found that ΔR/Ps of the forests with lower ecosystem complexity were significantly more sensitive to both DECREASE and INCREASE treatments than the forests with higher ecosystem complexity (Fig. 1). This is consistent to a former study in which water yield was more sensitive to climate change in planted forests (lower complexity) than in natural forests (higher complexity)42. Similarly, evaluating whether ΔR/P differs significantly from zero, we found that DECREASE treatment significantly increased the R/Ps of seven combinations with relatively low ecosystem complexity, but this is not the case for the other one with relatively high ecosystem complexity (Fig. 1a–c). In parallel, INCREASE treatment significantly decreased the R/Ps of only three combinations with the lowest ecosystem complexity, but this did not apply to the other five ones with relatively high ecosystem complexity (Fig. 1d–f).

Generally, DECREASE treatment increases ΔR/Ps, whereas INCREASE treatment decreases ΔR/Ps (Fig. 1 and Supplementary Table 3). It is worth to note that R/P was more sensitive to DECREASE treatment than to INCREASE treatment (Fig. 1), which has been acknowledged since a long time ago31, but the underlying mechanism has not been clarified. The asymmetric effects are also closely related to ecosystem complexity changes caused by the two treatment events. Nevertheless, ecosystem complexity is a synthetic indicator depending on multiple factors and can barely be described by a single metric. Therefore, despite that DECREASE or INCREASE treatment is performed at the same magnitude, ecosystem complexity might not respond proportionally, which in turn exerts different hydrological effects on the system. This finding suggests that the impacts propagated through the chain from treatments, via the alteration of ecosystem complexity, to the changes of hydrological effect is highly nonlinear. Such a complexity, however, is also the primary cause of the high uncertainties of the forest–water relationship. A previous study has demonstrated theoretically that R/P changes nonlinearly with the variation of watershed characteristics17, but the nonlinear relationship between treatments and alteration of ecosystem complexity was not recognized.

Natural growth does not decrease forest R/Ps

Surprisingly, our results revealed no temporal downtrend for the mean R/Ps derived from the 278 watersheds under natural succession. This finding is robust as revealed by the comparisons of different combination categories, including the means of R/Ps derived from the two treatments, P/PET climate regions, forest origins, pre-treatment forest age categories and pre-treatment forest coverage categories (Fig. 3). Alternatively, significant increasing trends in the means of R/Ps were observed from the watersheds of lower forest coverages, various forest ages, and dryer regions (Fig. 3c, e, f). To sum up, our results suggest that natural growth of vegetation in forest ecosystems not only prevents water yield reduction in watersheds, but also increases the water yield to some extent depending on different climate conditions and underlying surface characteristics.

Ecological principles indicate that any surface on the Earth that is suitable for plant survival will be naturally covered by vegetation that is in harmony with the hydrothermal environment if given enough time28. Moreover, even if external disturbance happens, the vegetation complexity (including biomass, greenness, coverage, etc) will increase continuously and naturally after interference until reaching the climax stage28. This might explain the decrease and increase found in the periods before and after the 10 data-years (Fig. 3). Specifically, we suspect that forest natural succession process is artificially shortened by anthropogenic disturbance, such as tree-planting activities (e.g., INCREASE), but the ecosystem complexity has yet been fully developed in harmony with the environment (i.e., species composition, root structure, microbes, invertebrates) in the first decade after treatment. Thus, the water consumption increased due to asynchronous development of the disturbed ecosystem at first; however, as natural growth proceeds, co-evolving hydrological and other functional processes help reducing water consumption and increasing the runoff coefficient (Fig. 3). Therefore, in an ecosystem in progressive succession state over a long-term period, functional processes do not decline but are stable or increase28, especially key processes such as water cycling. Hence, our results revealed no reduced water yields for the 278 watersheds experiencing a progressive succession process (natural growth without anthropic disturbances over 25 years).

One may argue that such a pattern could be explained partially by elevated atmospheric CO2 concentration, which might promote stomatal closure43,44, reduce water consumption, and boost photosynthetic activities in trees40,45. Nonetheless, rising CO2 concentration alone is insufficient to explain why the increase of R/Ps was only found in the forest ecosystems with low forest coverages. In fact, a recent study reported that rising CO2 concentration would boost transpiration and reduce runoff in a region with limited vegetation cover since the effect of increases in leaf area index (LAI) on transpiration exceeds that of stomata closure in such ecosystems46. Thus, the enriching atmospheric CO2 concentration is unlikely the driver preventing forest ecosystem from water yield reduction during succession. Instead, a more plausible explanation may involve two contrasting processes that were sufficiently mutual-feedback during progressive succession. One process is the decreased water yield from enhanced transpiration and interception as vegetation complexity increases. The other process is reduced water demand for transpiration from coevolution of water cycling process and other ecological functional processes (see Supplementary Discussion). However, the enhancing effects of transpiration and interception by forest natural growth tend to diminish with the increasing “mitigation effect” as ecosystem succession progresses and ecosystem complexity develops. These mechanisms can also be explained by the changes of variables in potential evapotranspiration calculation using Penman equation47 (see Eq. (1)):

$${E}_{{{{{{\mathrm{p}}}}}}}=\frac{{\triangle \cdot R}_{n}+6.43\cdot \gamma \cdot (1+0.536\cdot {u}_{2})\cdot ({e}_{s}-{e}_{a})}{\lambda \cdot (\triangle +\gamma)}$$
(1)

where Ep is potential evapotranspiration, ∆ is the slope of the saturated vapor pressure curve, Rn is the net radiation, γ is the psychrometric constant, u2 is wind speed, es and ea are saturated and actual vapor pressure, respectively, and λ is the latent heat of vaporization. In Eq. (1), the values of u2 and es – ea will decrease definitely but that of Rn does not have determined change trends, which, thus, reduces Ep with forest ecosystem developing from pioneer stage to climax stage.

However, the coevolution among different ecological processes in ecosystems with progressive succession is omnipresent, even in the climax stage. With the coevolution of key functional processes such as water and carbon cycling processes, the species composition and population size, as well as seasonal and spatial community structures, are in adjustment to keep functional processes, including water cycling process, in the most effective state possible.

Perspective

Forests are not ephemeral and static entities, but exist for decades or centuries and grow constantly if no harsh external disturbances. Global PWE studies demonstrated well the effect of the abrupt forest change on water yields in an ephemeral time-span (Fig. 4). This study further highlights the effects of natural growth on water yields (Fig. 4). Our results indicate that global greening contributed from forest natural growth is unlike to decrease water resources. Even for the forestation in arid and semi-arid regions (P/PET < 1) or planting trees in the existing planted forests, the reduced water yield due to forestation will recover if no continuous artificial interference is exerted.

Fig. 4: Conceptual diagram of the approaches used to examine forest change impacts on water yield in this study.
figure 4

The classical paired-watershed experiments (PWE) focus on the effect of the abrupt forest change on water yields, whereas the effects of natural growth on water yields are highlighted in this study.

Full size image

The present study addressed highly variable results from PWE worldwide and the underlying mechanisms, in which a convergence was identified based on ecological principles. Our study newly revealed that natural succession of forest ecosystems along with natural growth does not decrease or even increases water resource. The findings provide a generalized understanding about the integrated forest–water interactions with a global perspective. It should be pointed out that our analyses of the natural growth impacts on water yield focus on a time-span of 25 data-years, while further studies are required if more, longer measurement data become available in the future. Despite such a limitation, these results encourage us to advocate that nature-based forest restoration should be highlighted in global greening projects, especially in regions vulnerable to water shortage. Our study will benefit the sustainable development and management for promotion of the global forests. Moreover, it also provides a roadmap to the forest managers to adjust the forest management measures and policies when evaluating the balance of water supply and forest ecosystem services.

Methods

Literature search

The peer-reviewed research literature was searched within the Web of Science and the Chinese Science Citation Database using the keyword combinations “forest hydrology” AND “runoff OR evapotranspiration OR water yield”. Compiling of the original database was performed from 2013 to 2015, including 234 and 1928 papers related to paired-watershed experiments (PWE) and temporal measurement/observation analyses, respectively17. Since 2015, ~800 additional papers were found to replenish the database using the same search terms. The nearly 3000 papers were filtered based on three inclusion criteria: (1) Watersheds information were either directly reported or could be found in other studies, especially details about the forest origins and treatment. Specifically, watersheds with more than 50% of trees artificially planted were defined as planted forest watersheds, whereas those with more than 50% of trees naturally originated were defined as natural forest watersheds. Felling and thinning were defined as DECREASE treatment (treatments directly reduced tree number), whereas planting trees was defined as INCREASE treatment (treatments favors increasing tree number). Forests free from any artificial intervention (DECREASE or INCREASE treatments) were classified as CONTROL. (2) The mean annual runoff data before treatments or at least the annual runoff data of the last year before treatments should be available. (3) No artificial interferences after treatments were exerted on watersheds but the hydrological changes were continuously monitored. After DECREASE or INCREASE treatments, there should be at least one year of annual runoff data or mean annual runoff data available, in which the measurement year should also be indicated. For watershed used for evaluating the impacts of forest natural growth on water yield, there should be at least 2 years of annual runoff data reported (see “Effects of natural growth on runoff coefficients” section below).

Following the application of the selection criteria, we identified 496 qualified watersheds scattered across 36 countries globally (Fig. 5 and Supplementary data 1). Almost all the data used were obtained directly from the corresponding research papers, except for a few watersheds in which annual precipitation was not clearly reported (see Supplementary Data 1 and 2). For watersheds with missing precipitation data, we extracted the information from the China Meteorological Data Service Centre (http://data.cma.cn/) and cross validated. We then calculated the yearly potential evapotranspiration (PET) and wetness index (P/PET) according to methods elaborated in previously published literatures48,49.

Fig. 5: Distribution of the 496 watersheds included in the study.
figure 5

DECREASE: felling or thinning treatment that directly reduces tree quantity; INCREASE: afforestation or reforestation treatment that favors increasing the number of trees; CONTROL: free from DECREASE, INCREASE treatments, and any other external disturbances during the studied time period.

Full size image

Runoff data

From the 496 watersheds, 216 watersheds are allocated to DECREASE treatment, 125 watersheds are allocated to INCREASE treatment, and the remaining 155 consist the CONTROL. Besides, 278 of the 496 watersheds have records of 2–25 years of annual runoff data after treatments, or they were continuously measured and the year is clearly indicated. For watersheds undergoing DECREASE or INCREASE treatment, the first data-year was defined as the year following the treatment. While for CONTROL watershed, the earliest year when the annual runoff data was available was defined as the first data-year.

Combinations of forests and climatic regions

To obtain robust conclusions from limited records, the 496 watersheds were grouped into eight contrast combination categories according to climate regions (P/PET < 1 and P/PET ≥ 1) and forest origins (planted and natural forests). The eight combinations are watershed categories of: (i) all planted forests, (ii) all natural forests, (iii) all forests in P/PET < 1, (iv) all forests in P/PET ≥ 1 regions, the planted forests in P/PET < 1 (v) or P/PET ≥ 1 (vi) regions, and the natural forests in P/PET < 1 (vii) or P/PET ≥ 1 (viii) regions. Each pair of the contrast combinations generally differed in ecosystem complexity (e.g., lower ecosystem complexity of natural forests is expected in P/PET < 1 region than in P/PET ≥ 1 region). We compared these combinations to explore the generalized relationship of ecosystem complexity and water consumption (i.e., ubiquity of the conclusions).

Effects of treatment events on runoff coefficients

Treatment effects were derived from DECREASE or INCREASE watersheds. The differences in mean annual runoff coefficients of pre- and post-treatments (DECREASE or INCREASE) were compared with zero (one-sample t-test) or between any two of the eight combinations aforementioned (two-sample t-test). Thus, eight one-sample t-tests and six two-sample t-tests were available, i.e., all planted forests vs. all natural forests, all forests in P/PET < 1 vs. P/PET ≥ 1 regions, the planted forests in P/PET < 1 vs. P/PET ≥ 1 regions, the natural forests in P/PET < 1 vs. P/PET ≥ 1 regions, the planted forests vs the natural forests in P/PET < 1 regions, and the planted forests vs the natural forests in P/PET ≥ 1 regions.

Effects of natural growth on runoff coefficients

The 278 watersheds with records of 2–25 years of annual runoff records were used to quantify the effects of natural growth on runoff coefficients (including 100 DECREASE, 100 INCREASE, and 78 CONTROL watersheds). The annual runoff coefficients in each combination category (i.e., all the 278 watersheds, two treatments, P/PET < 1 and P/PET ≥ 1 regions, planted and natural forests, two forest age groups, two forest coverage groups) were averaged for each data-year from the first data-year onward. Linear regression models were also applied to quantify the temporal trends of annual runoff coefficients under natural growth states.

Statistical analysis

One-sample t-test and two-sample t-test in SPSS v22.0 (IBM Corporation, USA) were adopted to examine whether the ΔR/Ps induced by treatment events are significantly different from zero and whether the difference in ΔR/Ps between two contrasting combinations (combinations with contrast origins or/and climate regions) is statistically significant. Pearson correlation analysis was applied to calculate the R2 and corresponding p-values between net changes of runoff coefficients (ΔR/P) and forest coverage or forest age, as well as the temporal trends of forest runoff coefficient (R/P), using OriginPro 2021 (OriginLab Corporation, USA).