Uneven decline in the hydrological efficiency of China's natural and plantation forests

The vegetation transpiration fraction (TF) links terrestrial water and carbon cycles, yet the sensitivity of TF to changes in leaf area index (LAI) (θ) and its hydroclimatic controls remain poorly constrained, particularly in China's contrasting natural forests (NF) and plantation forests (PF). Using forest-type maps to identify natural forests (NF) and plantation forests (PF), together with GLEAM-derived evapotranspiration variables and ERA5-Land hydroclimatic data from 1990–2020, we quantified the spatiotemporal patterns of θ and evaluated the relative roles of soil moisture (SM) and vapor pressure deficit (VPD). We found that θ increased spatially from humid to semi-arid regions and was consistently higher in NF than in PF, but declined significantly over time in both forest types, with a larger decrease in PF. Along the joint SM–VPD gradient, θ exhibited a nonlinear response surface, with higher values concentrated under moderate SM and intermediate-to-high VPD conditions. Although VPD still explained a slightly larger share of the present spatial pattern of θ, the influence of SM strengthened over time, indicating increasing soil moisture limitation. Overall, China's forests appear to be shifting toward a more conservative water-use regime, in which the capacity of additional leaf area to enhance TF has weakened. These results provide new insights into how contrasting forest types regulate water use under concurrent warming and greening, with implications for climate-adaptive forest and water-resource management.

1 Introduction

Terrestrial ecosystem evapotranspiration (E) is the primary pathway by which land surface moisture returns to the atmosphere and thus plays a central role in the global water cycle (Liu et al., 2023; Sun et al., 2022). E comprises three components: soil evaporation (E_s), vegetation canopy interception (E_i), and plant transpiration (T) (Niu et al., 2020; Wei et al., 2017). Among these, vegetation transpiration, by which water is released to the atmosphere through plant stomata, accounts for over 60 % of total terrestrial E (Li et al., 2024; Wei et al., 2017). As transpiration is closely linked to carbon assimilation and energy exchange, the vegetation transpiration fraction (TF, defined as $T/E$) quantifies the contribution of vegetation to land–atmosphere water–vapor flux and is a key indicator of vegetation–climate coupling strength (Schlesinger and Jasechko, 2014; Wei et al., 2017). Understanding changes in TF is crucial for revealing ecohydrological mechanisms and accurately predicting climate–change impacts.

Persistent climate warming and rising atmospheric CO₂ over the past few decades have driven global vegetation greening and altered ecosystem water balances (Denissen et al., 2022; Hu et al., 2023). Increases in Leaf area index (LAI) directly enhance canopy interception and transpiration potential; global plant transpiration is estimated to have increased by about 6 % from 1990–2020, primarily as a consequence of increased LAI (Chen et al., 2024). However, the positive effect of LAI on transpiration depends on water availability: when soil moisture (SM) is scarce or atmospheric drought, often represented by vapor pressure deficit (VPD), intensifies, plants close their stomata to suppress water loss, thereby reducing both transpiration and photosynthesis (Liu et al., 2020b; Zahra et al., 2023). SM represents the supply-side constraint, whereas VPD represents the atmospheric demand-side pull; these two factors often co-vary and can together impose compound drought stress on ecosystems (Song et al., 2024). High VPD induces partial stomatal closure to prevent excessive water loss and hydraulic failure, causing transpiration rates per unit leaf area to saturate or even decline under very high VPD conditions (Grossiord et al., 2020; Novick et al., 2016). Therefore, even as LAI increases, the additional leaf area struggles to further increase T under extreme atmospheric drought (Xu et al., 2023). Conversely, under moderate SM and appreciable atmospheric demand, vegetation can maintain both water supply and evaporative demand at levels that support transpiration, allowing increases in LAI to exert a stronger positive effect on TF (Liu et al., 2020a). During extreme soil drought, increased surface sensible heat further dries the near-surface atmosphere, creating positive SM–VPD feedback that exacerbates drought conditions (Qing et al., 2022; Zhou et al., 2019). Accurately characterizing the response of TF to changes in LAI therefore requires a unified framework that accounts for the synergistic and nonlinear effects of both soil water supply and atmospheric demand (Koehler et al., 2023). In recent years, discrepancies have emerged regarding the relative roles of SM and VPD: some studies have emphasized that atmospheric drought imposes greater limits on ecosystem water–carbon cycles (Novick et al., 2016), whereas others, after decoupling the coupled effects of VPD and SM, have found that soil moisture is the dominant factor, especially in semi-arid regions, where SM typically imposes a stronger limitation on productivity (Liu et al., 2020a). Because SM and VPD are often strongly correlated, quantifying their relative roles under coupled hydroclimatic conditions remains a central challenge in ecohydrological research.

Since the 1970s, China has implemented some of the world's largest-scale afforestation and ecological restoration programs, increasing forest cover from approximately 12 % in the 1970s to over 22 % in recent years (Cheng et al., 2025). China's forests now comprise extensive natural and planted forests. Natural forests are concentrated mainly in the mountainous regions of northeastern and southwestern China, whereas planted forests are more widespread in eastern, central, and southern China (Cheng et al., 2024b). Compared with natural forests, planted forests in China are more often younger, structurally simpler, and dominated by single-species or even-aged stands under more intensive management (Cheng et al., 2024a; Farooq et al., 2021). These differences in origin, structure, and water-use strategies may underlie substantial disparities between NF and PF in soil water acquisition, stomatal regulation, and hydraulic safety margins, thereby leading to different responses of TF to changes in LAI (hereafter, LAI–TF sensitivity, denoted as θ). Exploring these disparities between forest types can deepen the understanding of vegetation water-use mechanisms and improve regional water-resource management, afforestation benefit assessment, and climate-adaptive forestry strategies. However, systematic research on this issue is currently lacking. Existing studies have mostly examined the sensitivity of transpiration or productivity to drought stress at global or broad regional scales (Liu et al., 2020a; Novick et al., 2016), or have focused on ecohydrological processes in specific arid regions, without a comparative assessment of θ across forest types along a unified climatic gradient. This raises four interrelated questions: (1) Are there systematic differences in θ between NF and PF along the climatic gradient from humid to semi-arid conditions? (2) How have these differences changed over the last 30 years? (3) Under coupled changes in SM and VPD, what are the relative roles of these two factors in shaping θ, and do these roles vary across climatic backgrounds or over time? (4) Through which hydrological pathways (SM and/or VPD) do changes in macroclimate, such as radiation, precipitation, temperature, and wind speed, indirectly affect θ? Addressing these questions will help fill important research gaps and deepen our understanding of forest–water relations under climate change.

The objectives of this study were to quantify the spatiotemporal patterns of transpiration sensitivity (θ) in China's natural forests (NF) and plantation forests (PF) from 1990 to 2020 and to identify the hydroclimatic controls underlying these patterns. To this end, we integrated multi-source remote-sensing and reanalysis data to derive growing-season LAI, evapotranspiration components, and meteorological variables for forest pixels at 0.1° resolution. We first used the aridity index (AI) framework to characterize the climatic background across four climatic zones and then focused the statistical analyses on the humid, semi-humid, and semi-arid zones to ensure comparability. We then applied complementary analytical approaches at different statistical scales: pooled spatial binning was used to characterize the response surface of θ along the joint SM–VPD gradient and to compare spatial contrasts under comparable hydroclimatic backgrounds, whereas sliding-window partial correlation and ridge regression were used to quantify the independent temporal associations of SM and VPD with θ and their evolution over time. Finally, pathway analysis was used to decompose the indirect effects of macroclimate change on θ through local hydrological pathways. All trends were estimated using the Theil–Sen slope and evaluated with the Mann–Kendall test. We further assessed robustness through cross-product comparison and examined whether atmospheric CO₂ and stand age altered the main interpretation. Together, this framework was designed to reveal both the spatial heterogeneity and temporal evolution of hydroclimatic controls on θ in China's contrasting forest types.

2 Materials and methods

2.1 Data

We used datasets describing forest type, canopy structure, evapotranspiration components, hydroclimate, and several auxiliary background factors (Table 1). Forest type data were obtained from the 1990–2020 spatial distribution dataset of natural and plantation forests in China (Cheng et al., 2024b). This dataset was derived from Landsat imagery and extensive forest inventory samples using machine-learning classification and was resampled to 0.1° resolution to match the hydroclimatic datasets. Pixels with 100 % coverage of either natural forest (NF) or plantation forest (PF) were treated as pure pixels, and only these pixels were included in all major analyses. To maintain temporal consistency between forest masks and climate data, we adopted a period-based dynamic matching strategy. Specifically, forest maps for 1990, 1995, 2000, 2005, 2010, 2015, and 2020 were used to represent the periods 1990–1994, 1995–1999, 2000–2004, 2005–2009, 2010–2014, 2015–2019, and 2020, respectively.

Table 1 Data overview.

LAI, Leaf area index; E, Evapotranspiration; T, Transpiration; SM, Root-zone soil moisture; T_a, Air temperature; T_d, Dewpoint temperature; P_a, Air pressure; P, Precipitation; PET, Potential evapotranspiration; R_n, Net radiation; WS, wind speed.

LAI was derived from the GIMMS LAI4g dataset (Cao et al., 2023). The main evapotranspiration dataset used in this study was GLEAM v4.2a, from which transpiration (T), evapotranspiration (E), and soil moisture (SM) were extracted. GLEAM separately estimates the major components of terrestrial evaporation and provides both surface and root-zone soil moisture. In this study, SM refers to the root-zone soil moisture product, rather than surface soil moisture, because it more closely represents plant-available water relevant to transpiration regulation at the seasonal to interannual scales considered here. In GLEAM, root-zone soil moisture is represented as a vegetation-accessible multi-layer soil profile constrained by assimilated surface observations, rather than as a single shallow soil layer. Near-surface meteorological variables, including air pressure (P_a), relative humidity (RH), air temperature (T_a), precipitation (P), dewpoint temperature (T_d), wind speed (WS), net radiation (R_n), and potential evapotranspiration (PET), were obtained from ERA5-Land. All environmental variables were aggregated over the growing season (April–October) for 1990–2020, using only growing-season data and excluding non-growing-season observations. Means were used for state variables (e.g., LAI, SM, T_a, and VPD), whereas cumulative values were used for flux variables (e.g., E, T, and P).

To evaluate the robustness of the estimated θ patterns and trends to data-product choice, we additionally used two alternative gridded products for cross-product comparison in supplementary analyses. One was the Simple Terrestrial Hydrosphere v2 (SiTHv2) product, which provides independent estimates of evapotranspiration and transpiration (Zhang et al., 2024). The other was the China terrestrial ecosystem transpiration fraction dataset (Niu et al., 2020), from which an alternative transpiration-fraction-based θ estimate was derived for comparison. These two datasets were used only for robustness assessment of θ, whereas all main calculations in the study were based on GLEAM.

To assess whether long-term changes in θ may also be modulated by non-hydroclimatic background factors, we further compiled several auxiliary datasets for supplementary analyses, including a global 1 km atmospheric carbon dioxide concentration dataset (Wang, 2026), a long-term reconstructed forest age dataset for China (Xia et al., 2023), and a 2020 forest age spatial distribution dataset for China (Cheng et al., 2024a). These auxiliary datasets were resampled to 0.1° resolution and used only to assess whether CO₂ and stand age materially altered the observed θ patterns and trends. They were not included in the main analytical framework, in which the focus was on hydroclimatic controls.

2.2 Study region

Our study area covers terrestrial China, spanning temperate, subtropical, and tropical climate zones. The analysis focuses on China's NF and PF. NF are mainly concentrated in the mountains of Northeast and Southwest China, whereas PF are widely distributed across the plains and hills of East, Central, and South China (Fig. 1a). The vast majority of NF and PF are located in humid and semi-humid zones (Fig. 2).

Figure 1 Spatiotemporal characteristics of China's natural forests (NF) and plantation forests (PF) from 1990 to 2020. (a) Average spatial distribution of NF (green) and PF (orange) in China from 1990 to 2020. (b) Dynamic changes in the area (km²) of NF and PF from 1990 to 2020. (c) Rate of change (%) in NF and PF area for each five-year period from 1990–2020.

Figure 2 Spatiotemporal distribution and composition of forests in different climatic zones in China from 1990 to 2020. (a) Spatial distribution of the four climatic zones (humid, semi-humid, semi-arid, and arid) in China; the pie charts show the total forest area within each zone and its composition of natural forest (green) and planted forest (orange). (b) Dynamic changes in forest area in each climatic zone from 1990 to 2020. For each year, the left solid stacked bar represents natural forest (NF), and the right hatched stacked bar represents plantation forest (PF). Bar colors indicate the forest area contributed by different climatic zones. The value at the top of each year indicates the total national forest area, and labels with “K” denote 10³ km².

We used the aridity index (AI) to characterize climatic background and to organize the spatial analyses within a unified framework. AI was defined as the ratio of annual potential evapotranspiration to annual precipitation ($\text{AI}=\text{PET}/P$). Based on the 1990–2020 multi-year mean AI, China was divided into four climatic zones: humid (AI<1), semi-humid ($\mathrm{1}\le \text{AI}<\mathrm{1.5}$), semi-arid ($\mathrm{1.5}\le \text{AI}<\mathrm{4}$), and arid (AI≥4). This four-zone classification was used as the geographic framework for mapping and descriptive comparisons. Because forest samples in the arid zone were extremely sparse, subsequent statistical analyses were restricted to the humid, semi-humid, and semi-arid zones.

2.3 Methods

2.3.1 VPD calculation

Vapor pressure deficit (VPD, hPa) is expressed, as follows, as the difference between saturation vapor pressure (P_s, hPa) and actual vapor pressure (P_w, hPa):

$$\begin{array}{}\text{(1)}& \text{VPD}=P_{\mathrm{s}}-P_{\mathrm{w}}\end{array}$$

To calculate VPD, we used directly observed meteorological parameters, including air temperature (T_a, °C), relative humidity (RH, %), and air pressure (P_a, hPa). P_s was calculated using the improved Magnus equation (Yuan et al., 2019):

$$\begin{array}{}\text{(2)}& P_{\mathrm{s}}=\mathrm{6.112}\times f\times e^{\frac{\mathrm{17.67}\times T_{\mathrm{a}}}{T_{\mathrm{a}}+\mathrm{243.5}}}\end{array}$$

Here, f is the atmospheric pressure enhancement factor, which corrects for the effect of pressure on saturation vapor pressure in moist air (Buck, 1981). This coefficient is calculated using P_a as follows:

$$\begin{array}{}\text{(3)}& f=\mathrm{1}+\mathrm{7}\times {\mathrm{10}}^{-\mathrm{4}}+\mathrm{3.46}\times {\mathrm{10}}^{-\mathrm{6}}\times P_{\mathrm{a}}\end{array}$$

P_w is then derived from P_s and RH as follows:

$$\begin{array}{}\text{(4)}& P_{\mathrm{w}}=P_{\mathrm{s}}\times \text{RH}\times {\displaystyle \frac{\mathrm{1}}{\mathrm{100}\mathit{\%}}}\end{array}$$

2.3.2 Sensitivity of TF to LAI

We define the transpiration fraction as

$$\begin{array}{}\text{(5)}& \text{TF}={\displaystyle \frac{T}{E}}\end{array}$$

In this study, T and E were first aggregated from daily data to annual growing-season totals (April–October), and LAI was represented by the corresponding annual growing-season mean. TF was then calculated from these annual growing-season quantities, and all subsequent θ analyses were therefore conducted at the annual growing-season scale rather than at the daily or monthly scale.

Note that $\text{TF}\in (\mathrm{0},\mathrm{1})$ is bounded. To explicitly account for the inherent nonlinearity and saturation in the TF–LAI relationship, we estimated the sensitivity θ as a marginal response under a bounded nonlinear framework.

Within the full period (1990–2020) and within each 11-year moving window, we fit a quadratic model in logit space:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}z=& \text{logit}\left(\text{TF}\right)=\ln \left({\displaystyle \frac{\text{TF}}{\mathrm{1}-\text{TF}}}\right)\\ =& a+b\left(\text{LAI}-{\text{LAI}}_{\text{med}}\right)+c{\left(\text{LAI}-{\text{LAI}}_{\text{med}}\right)}^{\mathrm{2}}\end{array}\end{array}$$

Here, LAI_med is the pixel-specific median LAI within the same period (median-centering was used to improve numerical stability). Prior to the logit transform, TF was clipped to $[\mathit{\epsilon },\mathrm{1}-\mathit{\epsilon }](\mathit{\epsilon }=\mathrm{0.005})$ to avoid numerical singularities.

We then recovered $\widehat{\text{TF}}=\mathrm{1}/(\mathrm{1}+e^{-z})$ and computed the marginal sensitivity as

$$\begin{array}{}\text{(7)}& \mathit{\theta }={\displaystyle \frac{\partial \text{TF}}{\partial \text{LAI}}}=\widehat{\text{TF}}\left(\mathrm{1}-\widehat{\text{TF}}\right)[b+\mathrm{2}c\left(\text{LAI}-{\text{LAI}}_{\text{med}}\right)]\end{array}$$

For ease of interpretation, Eq. (7) can be written as $\mathit{\theta }=\underset{\text{damping}}{\underbrace{\widehat{\text{TF}}(\mathrm{1}-\widehat{\text{TF}})}}\times \underset{\mathit{\eta }}{\underbrace{[b+\mathrm{2}c(\text{LAI}-{\text{LAI}}_{\text{med}}\left)\right]}}$, in which the damping term reflects boundary effects and η represents the structural sensitivity in logit space.

To separate the position-dependent saturation effect from changes in the shape of the fitted TF-LAI relationship, we evaluated (i) θ_total at LAI = LAI_med, representing the sensitivity under the actual state of each pixel and (ii) θ_clim at a fixed reference LAI_ref (global median LAI), which removes the influence of shifting along the TF–LAI curve. Additionally, the boundary component was quantified as

$$\begin{array}{}\text{(8)}& {\mathit{\theta }}_{\text{bound}}={\mathit{\theta }}_{\text{total}}-{\mathit{\theta }}_{\text{clim}}\end{array}$$

We calculated θ at two time scales. (1) The full-period mean θ characterized the average state during 1990–2020. (2) The θ time series characterized decadal changes, generated using an 11-year moving window recorded for the central year (1995–2015), a data processing step that was also used for subsequent trend and attribution analyses. Pixels with insufficient valid years and negligible LAI variability were excluded to ensure robustness.

We further used a fixed-curve counterfactual decomposition to separate the temporal change in θ_total into a position/saturation effect and a shape effect. The position effect describes changes in θ caused by LAI moving along an otherwise fixed TF–LAI curve, and therefore reflects the bounded nature of TF and the saturation compression associated with TF(1−TF). In contrast, the shape effect describes changes in θ caused by temporal changes in the fitted TF–LAI curve itself, as reflected by changes in the logit-quadratic coefficients.

For each forest type, the first 11-year window was used as the reference curve. The fixed-curve prediction was calculated by holding the regression coefficients of this reference curve constant while allowing LAI to vary through time. The position effect was then defined as the change predicted by this fixed-curve model, whereas the shape effect was defined as the residual between the observed θ_total trajectory and the fixed-curve prediction. Thus,

$$\begin{array}{}\text{(9)}& \mathrm{\Delta }{\mathit{\theta }}_{\text{total}}=\mathrm{\Delta }{\mathit{\theta }}_{\text{position}}+\mathrm{\Delta }{\mathit{\theta }}_{\text{shape}}\end{array}$$

Where Δθ_position represents movement along a fixed bounded response curve and Δθ_shape represents changes in the response-curve shape.

2.3.3 Spatial association under the joint SM–VPD gradient: binning analysis

To assess the relative roles of SM and VPD in the spatial differentiation of θ, we applied a pooled spatial-binning analysis. Using all eligible pure NF and pure PF pixels as samples, we extracted the full-period mean θ, SM, and VPD values and z-score standardized them. To evaluate the relative effect of VPD under comparable SM backgrounds, denoted as Δθ (VPD|SM), pixels were first grouped into bins according to SM. Within each SM bin, the difference in mean θ between high- and low-VPD groups was calculated, and these differences were then averaged across bins. The relative effect of SM, denoted as Δθ(SM|VPD), was evaluated in the same way by binning pixels according to VPD and then comparing high- and low-SM groups within each bin. This approach thus summarized spatial contrasts under comparable hydroclimatic backgrounds using a control-variable logic. These analyses are complementary but not identical to the subsequent temporal attribution analyses; specifically, the binning analysis characterized cross-site contrasts in θ across the long-term SM–VPD climate space, whereas the partial-correlation and ridge analyses quantified within-pixel temporal associations of θ with SM and VPD.

2.3.4 Dynamic temporal attribution: partial correlation and ridge regression

To investigate the pixel-scale temporal relationships of SM and VPD with θ, we used partial correlation and ridge regression based on the 1995–2015 11-year sliding-window series. This analysis was conducted at two levels. First, for each pixel, we calculated an overall partial correlation coefficient and ridge-regression coefficient based on the full 1995–2015 series, representing the average independent association of SM and VPD with θ throughout that period after accounting for their covariation. Second, to evaluate change over time, we performed a secondary sliding-window analysis on the same 1995–2015 series and generated time series of the coefficients themselves. Trends in these coefficient series were then used to infer whether the effects of SM and VPD became stronger or weaker over time. Thus, the first-level analysis described the mean temporal attribution pattern, whereas the second-level analysis characterized its change over time.

Partial correlation analysis: We calculated two pixel-wise partial correlation coefficients: $r(\mathit{\theta },\text{SM}|\text{VPD})$ and $r(\mathit{\theta },\text{VPD}|\text{SM})$. The former represents the correlation between θ and SM after controlling for interannual VPD fluctuations, whereas the latter represents the correlation between θ and VPD after controlling for SM variability.

Ridge regression analysis: Given the collinearity between SM and VPD, we further employed a ridge regression model, $\mathit{\theta }={\mathit{\beta }}_{\text{SM}}\cdot \text{SM}+{\mathit{\beta }}_{\text{VPD}}\cdot \text{VPD}$, to more robustly quantify their relative contributions to θ. By introducing an L2 regularization term, ridge regression stabilizes coefficient estimates under multicollinearity. The regression was performed on standardized time series, and the resulting coefficients (β_SM and β_VPD) directly reflect the relative strengths of the SM and VPD effects. We further calculated the relative contribution (RC) of SM and VPD as ${\text{RC}}_{\text{SM}}=\left|{\mathit{\beta }}_{\text{SM}}\right|/\left(\right|{\mathit{\beta }}_{\text{SM}}|+|{\mathit{\beta }}_{\text{VPD}}\left|\right)$ and ${\text{RC}}_{\text{VPD}}=\left|{\mathit{\beta }}_{\text{VPD}}\right|/\left(\right|{\mathit{\beta }}_{\text{SM}}|+|{\mathit{\beta }}_{\text{VPD}}\left|\right)$, to determine the dominant hydrological driver at each pixel. To further diagnose multicollinearity between SM and VPD, we also calculated the variance inflation factor (VIF).

To quantify temporal changes in θ and in the coefficient series derived from the 11-year moving-window analysis, we estimated trends using the Theil–Sen median slope estimator. Trend significance was assessed using the Mann–Kendall (MK) test (two-tailed, α=0.05). To reduce potential bias caused by temporal autocorrelation, the MK statistic was variance-corrected. The same trend-analysis framework was applied consistently to all moving-window-derived time series used in the temporal attribution analysis.

2.3.5 Pathway analysis of macroclimate effects on θ

To analyze the pathways through which macroclimate change may affect θ via local hydrological processes, we constructed a pathway model. The independent variable X served as the climate-factor trend (including P, T_a, R_n, and WS), the mediator variable M corresponded to the local water supply–demand status (SM and VPD trends), and the dependent variable Y represented the θ trend. RH was not included separately because its effect on atmospheric water demand is already represented by VPD. Pathway analysis was used to partition the associations between macroclimate trends and θ into direct effects and indirect effects mediated through SM and VPD. The statistical significance of all pathway effects was tested by bootstrap resampling (1000 iterations). Because SM and VPD may remain correlated under coupled hydroclimatic conditions, the pathway coefficients were interpreted as complementary association pathways rather than as a strict causal separation of two fully independent mediators. This interpretation was further evaluated using a supplementary mediator-specification comparison based on SM-only, VPD-only, and joint pathway models (Fig. S19 in the Supplement). The trends of P, T_a, R_n, WS, SM, VPD, and θ used in the pathway analysis were all estimated using the same Theil–Sen and variance-corrected Mann–Kendall framework described above.

3 Results

3.1 Spatiotemporal patterns and long-term trends of θ

Mean θ showed clear spatial heterogeneity and was consistently higher in NF than in PF (Fig. 3a, b, and e). High-θ areas in NF were concentrated mainly within the semi-arid to semi-humid transition belt, whereas PF showed a distinctly more fragmented pattern. During 1990–2020, growing-season θ declined widely in both forest types, with a stronger mean decline in PF than in NF (Fig. 3c and d). Significant negative trends were concentrated mainly in semi-humid and semi-arid transition regions, while humid regions showed weaker negative trends. The zone-averaged time series further confirmed that θ decreased significantly across climatic zones in both forest types, with the strongest declines occurring under drier climatic conditions (Fig. 3f and g). Overall, the marginal enhancement of TF by increasing LAI weakened over the past three decades, especially in PF and in relatively dry regions.

Figure 3 Spatiotemporal patterns and long-term trends of LAI–TF sensitivity (θ) in China's natural forests (NF) and plantation forests (PF) during the 1990–2020 growing seasons. (a) and (b) show the multi-year mean spatial distribution of θ for NF and PF, respectively, derived from the overall logit–quadratic model using TF_E. (c) and (d) show the corresponding long-term trends in θ estimated from the 11-year moving-window time series using the Theil–Sen slope; black dots indicate significant trends based on the Mann–Kendall test (p<0.05). (e) shows the distribution of θ across climatic zones (humid, semi-humid, and semi-arid) for NF and PF. (f) and (g) show the temporal evolution of zone-averaged θ for NF and PF, respectively, with Theil–Sen trend lines overlaid.

The decomposition of θ_total provided additional insight into the mechanisms behind the observed decline (Fig. 4). Under the main TF_E definition, the long-term decrease in θ_total reflected both position-dependent saturation and changes in the fitted TF–LAI curve shape, but their relative roles differed between forest types. In natural forests, the position/saturation component made a negative contribution to θ_total, whereas the shape component differed in sign and partly offset the decline. This contrast indicates that the total θ trend would be difficult to interpret from θ_total alone. In plantation forests, the position/saturation component also contributed to the weakening of θ_total, while the shape component indicated additional changes in the fitted TF–LAI sensitivity. These results show that the weakening of the LAI enhancement effect on TF is not merely a bounded geometric consequence of higher LAI or TF approaching saturation; it also reflects changes in the biophysical TF–LAI relationship under evolving hydroclimatic conditions.

Figure 4 Decomposition of temporal changes in LAI–TF sensitivity (θ) into position/saturation and shape/climate-related components. (a, b) Time series of spatially averaged θ_total, θ_clim, and θ_bound for natural forests (NF) and plantation forests (PF), respectively, under the main TF_E definition. (c, d) Counterfactual decomposition of Δθ_total into Δθ_position and Δθ_shape for NF and PF, respectively. The fixed-curve null model holds the logit-quadratic regression coefficients at their first-window values while allowing LAI to vary, so Δθ_position reflects movement along a fixed bounded TF–LAI curve. Δθ_shape is the residual between the observed trajectory and the fixed-curve prediction and reflects changes in the fitted response-curve shape. Lines indicate Theil–Sen trends; significance was evaluated using the Mann–Kendall test.

3.2 Hydroclimatic controls of θ: SM and VPD

3.2.1 Response of θ to the joint SM–VPD gradient

Pooled binning analysis showed that θ varied nonlinearly along the joint SM–VPD gradient in both forest types, with relatively high values occurring under intermediate SM and moderate-to-high VPD conditions, and lower values toward both the wet/low-VPD and dry/high-VPD ends of the gradient (Fig. 5a and b). The contrast analysis further revealed clear differences in the relative effects of SM and VPD (Fig. 5c). At the national scale, NF showed a much stronger SM-related contrast than VPD-related contrast [$\mathrm{\Delta }\mathit{\theta }\left(\text{SM}\right|\text{VPD})=-\mathrm{0.396}$ versus $\mathrm{\Delta }\mathit{\theta }\left(\text{VPD}\right|\text{SM})=-\mathrm{0.065}$], whereas PF showed the opposite pattern, with a much stronger VPD-related contrast [$\mathrm{\Delta }\mathit{\theta }\left(\text{VPD}\right|\text{SM})=\mathrm{0.228}$ versus $\mathrm{\Delta }\mathit{\theta }\left(\text{SM}\right|\text{VPD})=\mathrm{0.009}$]. Across climatic zones, NF exhibited stronger VPD-related contrasts in the humid and semi-arid zones but stronger SM-related contrasts in the semi-humid zone. In contrast, PF was mainly VPD-dominated in the humid and semi-arid zones, while both contrasts were relatively weak in the semi-humid zone. Within the pooled spatial-binning framework, these results indicate that the relative roles of soil water supply and atmospheric demand differed markedly between NF and PF. Zone-specific pooled SM–VPD binning patterns for NF and PF are shown in Figs. S14 and S15, respectively.

Figure 5 Pooled spatial binning analysis of LAI–TF sensitivity (θ) under the joint soil moisture (SM)–vapor pressure deficit (VPD) gradient during the 1990–2020 growing seasons. (a) and (b) show the mean standardized θ values in 10×10 SM–VPD bins for natural forests (NF) and plantation forests (PF), respectively; only bins with sample size ≥10 are shown. (c) shows the relative effects of VPD and SM as indicated by Δθ(VPD|SM) and Δθ(SM|VPD), calculated from the top 20 % versus bottom 20 % contrasts within each bin and then averaged across bins. Circles denote climatic zones, while squares denote the national total. θ was z-score standardized in this analysis.

3.2.2 Dynamic evolution of the independent effects of SM and VPD

Partial-correlation analysis revealed that the independent associations of θ with both SM and VPD were predominantly negative and were generally stronger in PF than in NF (Fig. 6). At the national scale, the mean values of $r(\mathit{\theta },\text{SM}|\text{VPD})$ and $r(\mathit{\theta },\text{VPD}|\text{SM})$ were −0.0725 and −0.0992 in NF, compared with −0.1517 and −0.1630 in PF, respectively, indicating slightly stronger VPD-related associations in both forest types and overall stronger hydroclimatic constraints in PF. Across climatic zones, the strongest negative partial correlations in NF occurred in the semi-arid zone, especially for VPD, whereas PF showed consistently negative correlations across all zones. The temporal analyses further revealed spatially heterogeneous changes in these relationships (Fig. 7). In NF, the SM-related partial correlation became less negative in the humid and semi-arid zones, while the VPD-related partial correlation also became less negative in the semi-humid zone. In PF, the SM-related partial correlation became more negative in the semi-humid and semi-arid zones, whereas the VPD-related partial correlation became less negative in the semi-humid zone. Overall, both SM and VPD constrained θ, but the constraints were generally stronger in PF, and their temporal changes were regionally differentiated rather than spatially uniform.

Figure 6 Spatial distribution and zonal means of partial correlations between forest LAI–TF sensitivity (θ) and hydroclimatic factors during the growing season from 1995 to 2015. (a) and (b) show the spatial patterns of $r(\mathit{\theta },\text{SM}|\text{VPD})$ for natural forests (NF) and plantation forests (PF), respectively; (c) and (d) show the corresponding patterns of $r(\mathit{\theta },\text{VPD}|\text{SM})$. (e) and (f) show the mean partial correlations across climatic zones. All variables were z-score standardized before analysis, and pixels with fewer than 10 years of valid data were excluded.

Figure 7 Temporal evolution and spatial trends of partial correlations between forest LAI–TF sensitivity (θ) and hydroclimatic factors during the growing season from 1995 to 2015. (a) and (d) show the temporal evolution of zone-averaged $r(\mathit{\theta },\text{SM}|\text{VPD})$ for NF and PF, respectively; (b) and (c) show the corresponding spatial trends. (e) and (h) show the temporal evolution of zone-averaged $r(\mathit{\theta },\text{VPD}|\text{SM})$ for NF and PF, respectively; (f) and (g) show the corresponding spatial trends. Trend slopes were estimated using the Theil–Sen estimator, and significant pixels according to the variance-corrected Mann–Kendall test (p<0.05) are indicated by black dots.

3.2.3 Quantification and spatiotemporal heterogeneity of dominant hydrological drivers

Ridge regression showed that both SM and VPD were generally negatively associated with θ, but the magnitude was usually stronger for VPD (Fig. 8). At the national scale, both NF and PF remained slightly more VPD-dominated in terms of the present spatial pattern, with mean relative contributions of VPD of 0.541 and 0.534 in NF and PF, respectively. Consistently, VPD-dominated pixels accounted for 57.1 % of all NF pixels and 57.0 % of all PF pixels, compared with 42.9 % and 43.0 %, respectively, for SM-dominated pixels. The relative-contribution and dominance statistics therefore consistently indicated a modest present-day advantage of VPD over SM, particularly in the semi-arid zone. The trend analysis further indicated a general strengthening of the SM effect and a weakening of the VPD effect in several regions, especially across semi-humid to semi-arid areas (Fig. 9). Within the pixel-wise temporal analysis, these results suggest that VPD explained a slightly larger share of the spatial pattern of θ, whereas the relative role of SM increased over time. RC and VIF supported the robustness of this interpretation (Fig. S16).

Figure 8 Spatial distribution of ridge-regression standardized coefficients and dominant hydrological driver for forest θ. For natural forests (NF), (a) and (b) show the spatial distributions of β_SM and β_VPD, respectively, and (c) shows the corresponding dominance pattern. (d)–(f) show the same results for plantation forests (PF). (g) and (i) show boxplots of pixel-level standardized β_SM and β_VPD across climatic zones for NF and PF, respectively. The central line denotes the median, boxes indicate the interquartile range, whiskers extend to 1.5 times the interquartile range, and points outside the whiskers are shown as outliers. (h) and (j) show the fractions of SM- and VPD-dominated pixels across climatic zones. Warm colors indicate positive coefficients and cool colors indicate negative coefficients; blue and red in dominance maps denote SM-dominated and VPD-dominated pixels, respectively.

Figure 9 Temporal evolution and spatial trends of ridge-regression coefficients for forest θ. (a) and (e) show the temporal evolution of zone-averaged β_SM for NF and PF, respectively; (b) and (f) show the corresponding spatial trends in β_SM. (d) and (h) show the temporal evolution of zone-averaged β_VPD for NF and PF, respectively; (c) and (g) show the corresponding spatial trends in β_VPD. Trend slopes were estimated using the Theil–Sen estimator, and significant pixels based on the Mann–Kendall test (p<0.05) are indicated by black dots. Solid trend lines indicate significant trends, and dashed lines indicate non-significant trends.

3.3 Pathway attribution of macro-climate effects on θ

Pathway analysis indicated that the effects of macroclimate change on θ varied strongly among drivers, climatic zones, and forest types (Fig. 10). In natural forests (NF), precipitation generally showed negative direct effects, especially in the semi-humid and semi-arid zones, while its indirect effects differed among pathways, including positive VPD-mediated effects in humid and semi-arid zones. Temperature was the most consistent negative direct driver in NF, with particularly strong negative effects in humid and semi-arid regions, while its indirect effects differed according to pathway, including positive VPD-mediated effects in humid and semi-arid zones. In plantation forests (PF), the pathway structure was more heterogeneous; specifically, precipitation showed a significant positive direct effect in the semi-arid zone, temperature showed opposite indirect effects through SM and VPD in the humid zone, and net radiation and wind speed mainly exerted pathway-specific effects in humid and semi-humid regions. Overall, the pathway coefficients indicated strong dependence on hydroclimatic context and forest type in how macroclimate trends were associated with changes in θ, with no single pathway dominating across all climate drivers and zones. A supplementary mediator-specification comparison showed that the direction of VPD-mediated indirect effects was generally more robust across single-mediator and joint models than that of SM-mediated effects, whereas SM-mediated pathways were more sensitive to model specification under SM–VPD coupling (Fig. S19).

Figure 10 Heatmap of pathway coefficients across climate zones for natural and plantation forests. Standardized path coefficients show the decomposition of the effects of four macroclimate drivers (precipitation, temperature, net radiation, and wind speed) on θ_total trends into direct effects and indirect effects mediated by SM and VPD for natural forests (a) and plantation forests (b). Rows correspond to the four macro-climate drivers, while columns are grouped by climate zone (humid, semi-humid, and semi-arid), each subdivided into direct effects, indirect effects via SM, and indirect effects via VPD. Cell values are standardized coefficients from parallel dual-mediation pathway analysis, and asterisks indicate significance at the 95 % level (i.e., bootstrap confidence interval excludes zero). Warm colors indicate positive effects, while cool colors indicate negative effects.

4 Discussion

4.1 Nonlinear hydroclimatic regulation of θ

The joint SM–VPD analysis shows that θ was highest under intermediate SM and moderate-to-high VPD, corresponding to relatively strong atmospheric demand, and declined toward both wetter/low-demand and drier/high-demand conditions (Figs. 5, S14 and S15). This pattern indicates that increased leaf area enhances TF most effectively when water supply is sufficient but not excessive and evaporative demand is strong enough to sustain canopy transpiration. Once soil water becomes limiting or atmospheric demand becomes too strong, stomatal regulation increasingly constrains transpiration to avoid hydraulic damage, so the marginal gain in transpiration per unit increase in LAI declines (Novick et al., 2016; Grossiord et al., 2020; McDowell et al., 2022). Under wetter and weak-demand conditions, the capacity of additional leaf area to increase TF is also reduced because energy limitations, weaker canopy–atmosphere coupling, and larger non-transpirational evaporation fractions all weaken the translation of additional leaf area into transpiration (Konings et al., 2017; Stoy et al., 2019). The θ decomposition and supplementary diagnostics further suggest that this nonlinearity cannot be explained solely as a bounded geometric effect of TF; part of the long-term change also reflects shifts in the LAI–TF relationship itself (Figs. 4 and S6–S13).

4.2 Long-term weakening of the LAI enhancement effect on TF

The widespread decline in θ, especially across semi-humid to semi-arid transition zones, indicates that the capacity of increasing LAI to enhance TF has weakened over the last three decades (Fig. 3). This trend is consistent with a background of warming, rising VPD, and increasing frequency of soil moisture limitation, all of which reduce the effective transpiration response of forests to additional canopy development (Yuan et al., 2019; Lian et al., 2020; Denissen et al., 2022). Rising atmospheric CO₂ may also have contributed to this decline by lowering stomatal conductance and increasing intrinsic water-use efficiency, thereby reducing transpiration demand per unit leaf area even where greening continued (Keenan et al., 2013; Keenan et al., 2016; Lavergne et al., 2019; Liang et al., 2023). At the same time, the supplementary analyses are more consistent with interpreting CO₂ as a temporally coherent background forcing rather than a spatially heterogeneous driver of θ across China's forests (Fig. S17 and Table S1 in the Supplement). Stand development may likewise have modulated this trajectory. Age-related changes in stand structure, rooting depth, canopy roughness, and interception can shift the partitioning of evapotranspiration between transpiration and non-transpirational components, and therefore alter the LAI–TF relationship even without large changes in total leaf area (Fan et al., 2017; Forrester, 2015). However, the age-stratified results show that the higher θ in NF than in PF persisted across age groups, and including stand age did not overturn the broad interpretation of the relative roles of SM and VPD, although the magnitude of change was heterogeneous across forest types and climatic zones (Fig. S18 and Table S2). Together, these results suggest that CO₂ rise and stand development may act as background modifiers of θ, while the dominant large-scale patterns remain primarily associated with hydroclimatic variability.

4.3 Temporal changes in hydroclimatic effects on θ

Partial correlation, ridge regression, and pathway analysis consistently indicate that the hydroclimatic controls on θ have changed over time, with the relative role of SM strengthening and that of VPD weakening in many regions (Figs. 6–10). This does not mean that atmospheric demand has ceased to matter; VPD still explains a slightly larger share of the present spatial pattern of θ, as also supported by the dominance statistics and multicollinearity diagnostics (Figs. 8 and S16). The more important point is that the sensitivity of θ to soil moisture limitation has intensified under recent climatic conditions. This interpretation is physically consistent with the increasing prevalence of water-limited ecosystem behavior under warming, in which rising evaporative demand is increasingly translated into ecological stress through soil-water depletion rather than through atmospheric forcing alone (Berg et al., 2016; Denissen et al., 2022). Once ecosystems approach critical soil-moisture thresholds more frequently, stomatal and hydraulic regulation become more tightly constrained by water supply, and the explanatory role of SM rises accordingly (Fu et al., 2022; Liu et al., 2025). The pathway results add the same message at a broader scale: temperature exerted generally negative effects on θ, whereas the effects of precipitation and radiation varied by forest type and climatic zone, implying that macroclimate trends influence θ through multiple hydroclimatic pathways rather than through a single dominant mechanism.

4.4 Contrasting ecohydrological responses of natural and plantation forests

The contrast between natural and plantation forests remained one of the clearest features of the analysis. NF consistently showed higher θ, and high-θ values in NF were maintained across a wider range of SM–VPD backgrounds, whereas PF showed lower θ and a more restricted range of hydroclimatic conditions under which high θ occurred. This pattern likely reflects differences between long-term spatial hydroclimatic contrasts and within-pixel temporal variability. Spatial contrasts integrate longer-term differences in stand structure, species composition, hydraulic diversity, and belowground water access, whereas temporal anomalies are expressed more directly through canopy–atmosphere coupling and stomatal–hydraulic regulation during drought (Martínez-Vilalta and García-Forner, 2017; Grossiord et al., 2020; Novick et al., 2024; Bachofen et al., 2024). Under this interpretation, the broader hydroclimatic response of NF may reflect the greater structural and functional heterogeneity often associated with natural forests, while the narrower response of PF may reflect the relative structural simplification typical of many plantation stands. This interpretation is also consistent with evidence that higher hydraulic diversity can buffer ecosystem drought responses (Anderegg et al., 2018). The age-stratified supplementary results are consistent with this interpretation, because the higher θ in NF than in PF persisted across age classes and this pattern is therefore unlikely to be explained solely by stand age.

4.5 Implications and limitations

These results imply that further greening will not necessarily translate into proportionally higher transpiration fractions under continued warming and drying. In regions where hydroclimatic water stress is increasing, especially across semi-humid to semi-arid transition zones, the ecohydrological consequences of increasing canopy leaf area depend increasingly on soil water supply, and this dependence appears more evident in plantation forests than in natural forests. This has direct management implications: plantation design may benefit from greater rooting-depth and trait diversity and from stand structures that reduce rapid soil-water depletion, whereas natural forests may benefit more from maintaining multilayer canopy structure and hydraulic diversity under rising atmospheric demand. Several limitations should remain explicit. A nationwide in situ benchmark for θ trends is not yet available because long-term observations that simultaneously resolve transpiration partitioning, LAI dynamics, and stand development remain too sparse across China's forests. We therefore evaluated robustness through cross-product consistency rather than direct site-based validation. In addition, because θ, SM, and VPD were analyzed from growing-season aggregates and moving-window series, the results characterize seasonal-to-interannual hydroclimatic controls rather than instantaneous stomatal responses to short-lived atmospheric fluctuations. Even with these constraints, the convergence among the main analyses, the supplementary robustness checks, and the CO₂- and age-sensitivity diagnostics supports the central inference that the enhancement of TF by increasing LAI has weakened across China's forests and is increasingly shaped by soil moisture limitation.

5 Conclusion

Across China's forests, the enhancement of TF caused by increased LAI has weakened during 1990–2020. VPD still explains a slightly larger share of the present spatial pattern of θ, but the role of SM has strengthened over this timespan, indicating there has been an increase in SM limitation. Differences between NF and PF were obvious: NF maintained higher θ values, whereas PF showed a greater long-term decline in θ. Collectively, these results indicate that continued greening under warming will not necessarily result in proportionally higher TF, particularly in semi-humid to semi-arid transition regions.