Enhanced understanding of dominant drivers of Water Yield change across China through the improved coupled carbon and water model

Rapid environmental changes, including climate change, escalating atmospheric CO₂ concentration ([CO₂]), and vegetation dynamics, have been significantly impacting hydrological processes. Yet it remains difficult to disentangle the contributions of climate, vegetation, and [CO₂] change to water yield (WY), especially clarifying [CO₂]-driven physiological effects. Therefore, this study improved the coupled carbon and water (CCW) model by integrating dynamic water use efficiency (WUE) to capture [CO₂]-physiological effects, and quantified the attribution of WY changes across China from 1982 to 2017 using a scenario analysis based on this model. The results showed that WY changes were attributed to the changes in climate, vegetation, and [CO₂]. Among them, climate change (especially precipitation change) emerged as the dominant driver, directly affecting over 70 % land area of China. The vegetation change was the second largest factor, especially in central China. The escalating [CO₂] had a relatively small effect. Spatially, the changes in vegetation and [CO₂] exerted greater influence within the 400–1600 mm precipitation range. Remarkably, under CMIP6 SSP585 projections, an accelerating [CO₂] rise will lead to a +1.29 % yr⁻¹ increase in annual WY by 2100 based on the elasticity estimate, which surpasses vegetation influence (−0.26 % yr⁻¹). By using dynamic WUE, this improved framework advances WY attribution and confirms that [CO₂] physiological regulation can partly offset vegetation-induced WY reductions. These findings deliver novel references for adaptive regional water resource management under climate change.

1 Introduction

The global environment has been undergoing rapid changes, impacting hydrological processes through climate change, escalating atmospheric CO₂ concentration [CO₂], and vegetation dynamics (Piao et al., 2007; Wei et al., 2024). Notably, China has experienced a visible greening trend in recent decades, prompting a heightened focus on ecological and water resource concerns (Chen et al., 2019). Investigating the influence of vegetation changes on runoff has thus emerged as a pivotal research area, aligning with China's increasing emphasis on environmental sustainability. Among hydrological metrics, water yield (WY) is especially relevant because it directly represents the amount of water remaining after evapotranspiration and is therefore closely linked to runoff generation and regional water availability. China's diverse climatic zones and pronounced greening make it an ideal natural laboratory for investigating these ecohydrological feedbacks, with implications for both China (Ogutu et al., 2021; Yang et al., 2019) and other semi-arid and monsoon-influenced regions such as the Sahel, South Asia, and the Mediterranean Basin (Nkiaka et al., 2025; Rahman et al., 2025; Serrano-Notivoli et al., 2022). Understanding how vegetation dynamics, climate change, and [CO₂] interact to regulate WY is therefore of both scientific and practical importance for water resource management and ecological restoration under accelerating environmental change.

Several methods have been employed to separate the effects of climate, vegetation, and [CO₂] change on runoff change, including paired catchment experiments, statistical methods, and modeling approaches (Zeng et al., 2020). Given that annual water yield (WY) equates to runoff through negligible soil water storage changes, these methodological evaluations directly inform WY attribution frameworks (Zhang et al., 2022b). The paired catchment experiment method, though classical, is limited to small-scale watersheds and is less applicable to larger regions (Peng et al., 2016). Statistical methods, while helpful in identifying correlations, lack a physical basis and are insufficient for explaining the underlying mechanisms of runoff changes (Chen et al., 2022). Modeling approaches for attribution fall into two broad classes: (i) process-based models that explicitly simulate coupled water–energy–carbon processes, and (ii) conceptual models that approximate these processes with parsimonious, physically interpretable relationships (Zhai and Tao, 2021). Process-based models can capture detailed mechanisms, but they require extensive inputs and many parameters, are sensitive to calibration and equifinality, and are computationally demanding – limitations that hinder basin-to-continental applications over long periods (Jiao et al., 2017; Ma et al., 2023). By contrast, conceptual models retain key ecohydrological mechanisms with far fewer parameters, scale well to large regions, and thus are well suited for large-scale attribution while preserving physical interpretability. Among these conceptual models, the Budyko framework, widely used to separate climate change effects on runoff, quantifies water balance through the aridity index (PET/precipitation) and incorporates a catchment-specific Budyko parameter (n) representing integrated land surface characteristics (e.g., vegetation, soil, topography) (Zhang et al., 2022a, 2016a). However, most Budyko-based applications primarily emphasize climate-driven attribution; vegetation and [CO₂] influences are typically introduced only indirectly – by assigning temporal changes in the Budyko parameter (n) to vegetation (Tan et al., 2024; Xue et al., 2022; Zhou et al., 2023) or correlating the Budyko parameter (n) with NDVI (Liu et al., 2024; Tan et al., 2023), and by embedding [CO₂] effects through PET adjustments (Liu et al., 2024). These practices conflate vegetation with other controls captured by the Budyko parameter (n) (e.g., soil, topography) and mix [CO₂]-physiological impacts with meteorological drivers in PET, making it difficult to isolate vegetation structural change from [CO₂]-induced stomatal adjustments and to ascribe mechanisms robustly (Gan et al., 2021). But accurately representing the physiological effects of elevated [CO₂] is critical for runoff attribution, because rising [CO₂] can reduce stomatal conductance, thereby suppressing vegetation transpiration and evapotranspiration (Lammertsma et al., 2011; Xu et al., 2016).

Meanwhile, the coupled carbon and water (CCW) model provides a suitable basis for mechanistic WY attribution, because its calculation framework explicitly couples carbon and water fluxes and links vegetation dynamics with evapotranspiration (Li et al., 2024b; Zhang et al., 2021b, 2022b). In contrast to empirical parameter n in the Budyko framework, the CCW model estimates actual evapotranspiration using gross primary productivity (GPP), vapor pressure deficit (VPD), and underlying water-use efficiency (UWUE), expressed as ET = GPP ⋅ VPD^0.5 $/$ UWUE. This carbon and water coupling enables ET to be directly derived from vegetation processes, making the attribution framework useful for separating vegetation effects from climate effects in WY. Nevertheless, the original CCW model, adopts a static UWUE and does not account for [CO₂]-induced physiological changes, specifically long-term enhancements in water use efficiency (WUE) resulting from elevated [CO₂], thereby limiting its capacity to isolate [CO₂] physiological effects from vegetation structural and climate influences (Adams et al., 2020; Li et al., 2023).

Therefore, our study prepares to enhance the CCW framework by incorporating dynamic WUE responses to [CO₂], allowing explicit attribution of runoff changes to three distinct drivers, namely climate change (e.g. precipitation, temperature, and so on), vegetation structural change (NDVI, and land use and land cover (LULC)), and [CO₂]-physiological effects (stomatal optimization). This extension advances beyond empirical or regression-based attribution, clarifies how [CO₂] modulates vegetation and hydrology interactions across large spatial scales, and provides policy-relevant evidence for sustainable water resource management and ecological restoration in China under accelerating environmental change.

2 Methods and Data

2.1 Data sources and processing

Four main datasets were employed in the improved CCW model: vegetation data (NDVI), climate data (precipitation, temperature, shortwave radiation, vapor pressure deficit, and atmospheric pressure), land use and land cover (LULC), and [CO₂]. The monthly NDVI dataset used in this study (Table 1) was derived from a daily 0.05° gap-free NDVI dataset in China (https://doi.org/10.6084/m9.figshare.c.7002225.v1, Li et al., 2024a), which was developed from the NOAA's daily NDVI dataset, applying effective data recognition and spatiotemporal gap-filling techniques. The dataset spans 1981–2023 and provides a spatial resolution of 0.05°, and we used bilinear interpolation to generate the dataset with a spatial resolution of 0.1°.

Climate data (Table 1), including precipitation, air temperature, surface downward shortwave radiation, relative humidity, and atmospheric pressure, were sourced from the China Meteorological Forcing Dataset (CMFD) at the National Tibetan Plateau Data Center (TPDC) of the Institute of Tibetan Plateau Research, Chinese Academy of Sciences (He et al., 2020). The dataset spans 1979–2018 and provides a spatial resolution of 0.1° and temporal resolutions at 3 h, daily, monthly, and annual scales. As the dataset did not provide vapor pressure deficit (VPD), we calculated VPD using the method from Howell and Dusek (1995), based on atmospheric pressure, temperature, and relative humidity.

LULC data (Table 1) were obtained from the Zhang et al. (2024) global dataset, which provides consistent multi-temporal global LULC maps at 30 m spatial resolution for 1985–2022. The dataset includes 35 fine-resolution LULC types. For the purposes of this study, and to facilitate LULC change analysis, we merged these 35 LULC types into 17 types using the IGBP classification, based on the method by Yang et al. (2017). Four primary LULC types – cropland, forest, grassland, and bare land – were determined following the method described by Mu et al. (2013). The data were resampled to the 0.1° spatial resolution, ensuring compatibility for modeling within the modified CCW framework.

[CO₂] data were sourced from the Mauna Loa Observatory (MLO), Hawaii (20° N, 156° W) (https://gml.noaa.gov/ccgg/trends/data.html, last access: 1 June 2026), with yearly observations used to represent national [CO₂] levels due to the minimal spatial variation in [CO₂] across China (Table 1). These datasets were then used to drive the improved CCW model.

In this study, the hydrological data for model validation from 145 hydrological stations (Fig. 1), each with at least 15 years of continuous data since 1982, was collected from the Hydrological Bureau of the Ministry of Water Resources of China (https://www.mwr.gov.cn/english/, last access: 18 June 2026). Annual runoff data were calculated from the daily runoff and the catchment area controlled by each hydrological station.

Table 1 Hydrology, climate, and vegetation data for the improved CCW model.

Figure 1 The geographic location and topography of the study area, where the black triangles mark the location of the hydrological gauging stations for model evaluation. Ten river basins considered in this study are: Songhua River basin (SHR), Liao River basin (LR), Hai River basin (HR), Huai River basin (HuR), Yangtze River basin (YZR), Yellow River basin (YR), Pearl River basin (PR), Southeast Rivers (SER), Southwest Rivers (SWR) and Northwest Rivers (NWR).

2.2 The improved CCW model

The Coupled Carbon and Water (CCW) model is a remote-sensing driven ecosystem model that mechanistically links terrestrial carbon uptake and water dynamics (Zhang et al., 2016b), and its detailed information is provided in Appendix A. In this study, we improve the CCW model by replacing the static underlying water use efficiency (UWUE) in the original CCW model with a dynamic ecosystem-scale WUE that responds to changing [CO₂], vapor pressure deficit (VPD), atmospheric pressure, and vegetation structure.

In the improved CCW model, evapotranspiration (ET) is calculated as the ratio of GPP to WUE:

$$\begin{array}{}\text{(1)}& \mathrm{ET}={\displaystyle \frac{\mathrm{GPP}}{\mathrm{WUE}}}\end{array}$$

where WUE is defined as the ratio of GPP (Gross primary productivity) to total ecosystem evapotranspiration. GPP is estimated as:

$$\begin{array}{}\text{(2)}& \mathrm{GPP}=\mathrm{APAR}\times \mathit{\epsilon }=\mathrm{PAR}\times \mathrm{FPAR}\times {\mathit{\epsilon }}_{\mathrm{pot}}\times R_{\mathrm{s}}\times T_{\mathrm{s}}\times W_{\mathrm{s}}\end{array}$$

where APAR is the absorbed photosynthetically active radiation (MJ m⁻²), which is calculated as the product of incident photosynthetically active radiation (PAR) and the fraction of PAR absorbed by vegetation (FPAR). PAR is typically assumed to be 45 % of the total shortwave radiation (Running et al., 2000), and FPAR is determined by the normalized difference vegetation index (NDVI) (Sims et al., 2005). ε is the realized light-use efficiency (g C MJ⁻¹), which is calculated by multiplying the potential light-use efficiency (ε_pot) and environmental scalars for diffuse radiation (R_s), temperature (T_s), and moisture stress (W_s).

WUE is estimated using the WEC (Water Efficiency and Carbon) equation proposed by Cheng et al. (2017):

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathrm{WUE}=& {\displaystyle \frac{C_{\mathrm{a}}\times P_{\mathrm{a}}}{\mathrm{1.6}\left(\mathrm{VPD}+g_{\mathrm{1}}\sqrt{\mathrm{VPD}}\right)}}\\ & \left[\mathrm{1}-\exp \left(-k\cdot \mathrm{LAI}\right)\right]\left(\mathrm{1}-f_{\mathrm{i}}\right)\end{array}\end{array}$$

where C_a is atmospheric CO₂ concentration (mol(CO₂) mol⁻¹(air)); P_a is atmospheric pressure (kPa); VPD is vapor pressure deficit (kPa); g₁ is an empirical parameter of the Ball stomatal conductance model; k is the radiation extinction coefficient, typically set at 0.6; LAI is the leaf area index; and f_i is fraction of interception to total evapotranspiration. In this study, the interception evaporation factor (f_i) was set to zero, since Cheng et al. (2017) reported that canopy interception account for a minor portion of total evapotranspiration at annual scales. Given that the improved CCW model focused on yearly water yield (WY) dynamics, neglecting interception loss reduces model complexity without substantially affecting WY estimation.

In order to ensure the consistency of NDVI and LAI trends, LAI is calculated using NDVI (Gutman and Ignatov, 1998):

$$\begin{array}{}\text{(4)}& \left\{\begin{array}{l}\mathrm{LAI}=-\mathrm{2}\ln \left(\mathrm{1}-f_{\mathrm{NDVI}}\right)\\ f_{\mathrm{NDVI}}=\frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{\mathrm{0}}}{{\mathrm{NDVI}}_{\mathrm{1}}-{\mathrm{NDVI}}_{\mathrm{0}}}\end{array}\right.\end{array}$$

where NDVI₀=0.04, NDVI₁=0.52

Under this framework, ET should be interpreted as an effective annual ecosystem-scale evapotranspiration suitable for long-term WY attribution, rather than a detailed partitioning of individual ET components. Compared with the original CCW model, this modification allows ET to vary with dynamic WUE and therefore enables the model to represent [CO₂]-induced physiological regulation in addition to vegetation structural and climatic effects.

Finally, the water yield (WY) is calculated as the difference between precipitation (P) and ET:

$$\begin{array}{}\text{(5)}& \mathrm{WY}=P-\mathrm{ET}\end{array}$$

On an annual scale, WY can be approximately equal to runoff, as changes in soil water storage over long periods (one year or longer) are considered negligible (Xiao et al., 2020; Zhang et al., 2021b). Thus, the attribution of WY can also be considered as the attribution of runoff. Accordingly, in this study, WY is used as the modelled output, while the term “runoff” is reserved for observed streamflow or literature values explicitly labelled as such. This approximation is most suitable for annual to multi-year analyses in regions without substantial long-term groundwater depletion or strong reservoir regulation.

2.3 Attribution analysis framework

In this study, we proposed an attribution analysis framework based on scenarios simulation to explore the combined and individual effects of climate, vegetation, and [CO₂] change on WY. Four scenarios were designed (Table 2). Scenario 1 (Actual) aimed to validate the improved CCW model and estimate the combined effects of climate, vegetation, and [CO₂] change on WY by allowing all factors to vary from 1982 to 2017. Scenario 2 (Vegetation Change) focused on estimating the direct effects of vegetation change on WY by allowing vegetation factors (NDVI and LULC) to vary while keeping climate and [CO₂] fixed at 1982 levels. In this case, the trend in WY obtained reflects the impact of vegetation change alone. Scenario 3 (Climate Change) aimed to estimate the direct effects of climate change on WY by allowing climate factors (precipitation, temperature, relative humidity, solar radiation, and atmospheric pressure) to change, while fixing vegetation and [CO₂] at 1982 levels. This scenario helps isolate the effects of climate change on WY. Scenario 4 ([CO₂] Change) was designed to estimate the direct effects of [CO₂] change on WY by varying [CO₂] levels from 1982 to 2017, while climate and vegetation factors were fixed at 1982 levels. The resulting WY trend reflects the impact of [CO₂] change alone. The resulting WY series under each scenario represents the direct impact of the corresponding driver under the assumption that the other drivers are fixed.

Table 2 Designed scenarios in the improved CCW model for WY attribution. LULC: Land use and land cover types; NDVI: Normalized difference vegetation index; TMP: Temperature; SRAD: Shortwave radiation; VPD: Vapor pressure deficit.

Note: The symbol “▴” denotes a changing input factor over time, whereas the symbol “△” represents a fixed input factor at the level of the initial year (1982).

For each scenario, the long-term “trend” in annual WY over 1982–2017 was quantified using the Theil–Sen slope of the annual WY series (Sen, 1968; Theil, 1992). The relative contributions of climate, vegetation, and [CO₂] to changes in WY were calculated using the following formula (Ma et al., 2023; Wang et al., 2022):

$$\begin{array}{}\text{(6)}& \left\{\begin{array}{l}{\mathrm{RC}}_{\mathrm{vegetation}}=\frac{{\mathrm{trend}}_{\mathrm{vegetation}}}{\left|{\mathrm{trend}}_{\mathrm{vegetation}}\right|+\left|{\mathrm{trend}}_{\mathrm{climate}}\right|+\left|{\mathrm{trend}}_{{\mathrm{CO}}_{\mathrm{2}}}\right|}\times \mathrm{100}\mathit{\%}\\ {\mathrm{RC}}_{\mathrm{climate}}=\frac{{\mathrm{trend}}_{\mathrm{climate}}}{\left|{\mathrm{trend}}_{\mathrm{vegetation}}\right|+\left|{\mathrm{trend}}_{\mathrm{climate}}\right|+\left|{\mathrm{trend}}_{{\mathrm{CO}}_{\mathrm{2}}}\right|}\times \mathrm{100}\mathit{\%}\\ {\mathrm{RC}}_{{\mathrm{CO}}_{\mathrm{2}}}=\frac{{\mathrm{trend}}_{{\mathrm{CO}}_{\mathrm{2}}}}{\left|{\mathrm{trend}}_{\mathrm{vegetation}}\right|+\left|{\mathrm{trend}}_{\mathrm{climate}}\right|+\left|{\mathrm{trend}}_{{\mathrm{CO}}_{\mathrm{2}}}\right|}\times \mathrm{100}\mathit{\%}\end{array}\right.\end{array}$$

where trend_vegetation, trend_climate, and ${\mathrm{trend}}_{{\mathrm{CO}}_{\mathrm{2}}}$ denote the Theil–Sen slopes of the annual WY series under vegetation, climate, and [CO₂] scenarios, respectively. These slopes represent the long-term rates of WY change attributable to each driver within the scenario framework; the relative contributions (RC_vegetation, RC_climate, and ${\mathrm{RC}}_{{\mathrm{CO}}_{\mathrm{2}}}$) are expressed as percentages, indicating the normalized magnitude of the scenario-derived direct effects.

At each grid, the dominant driver of WY change is identified as the factor with the largest absolute relative contribution among climate, vegetation, and [CO₂]. When the absolute difference in relative contributions between two factors do not exceed 5 %, these two factors are considered as joint significant contributors to the changes in WY at that grid (Jia et al., 2022). This criterion is used to distinguish grids dominated by a single factor from those where two drivers contributed comparably to WY change.

Furthermore, elasticity analysis is implemented to quantify the sensitivity of WY to individual drivers. Elasticity is defined as the relative change in WY caused by a 1 % perturbation in a given driver (Yang and Yang, 2011). We specifically focus on precipitation because, despite not always having the highest sensitivity, it is integral to the hydrological cycle and essential for assessing water yield (WY) under various climate change scenarios (Liu et al., 2017). Specifically, for each driver (x), including precipitation, NDVI, and [CO₂], the baseline S1 series is increased by 1 % over the entire 1982–2017 period while all other factors are held unchanged. Then the elasticity coefficient is estimated as:

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{\Delta }{R}_x}{R_x}}={\displaystyle \frac{{\mathrm{WY}}_{{\mathrm{mean}}_x}-{\mathrm{WY}}_{{\mathrm{mean}}_{\mathrm{S}\mathrm{1}}}}{{\mathrm{WY}}_{{\mathrm{mean}}_{\mathrm{S}\mathrm{1}}}}}\end{array}$$

where (x) represents precipitation, NDVI, or [CO₂]; ${\mathrm{WY}}_{{\mathrm{mean}}_{\mathrm{S}\mathrm{1}}}$ is the mean annual WY under the baseline scenario S1 during 1982–2017; and ${\mathrm{WY}}_{{\mathrm{mean}}_x}$ is the mean annual WY after increasing driver (x) by 1 % while keeping all other factors unchanged.

Mathematically, the elasticity coefficient is defined as the runoff change rate divided by 1 %, and the formula is as follows:

$$\begin{array}{}\text{(8)}& {\mathit{\epsilon }}_x={\displaystyle \frac{\frac{\mathrm{\Delta }{R}_x}{R_x}}{\frac{\mathrm{\Delta }x}{x}}}={\displaystyle \frac{\frac{\mathrm{\Delta }{R}_x}{R_x}}{\mathrm{1}\mathit{\%}}}\end{array}$$

Generally, while the scenario analysis above has identified which factors are most influential based on their relative contributions, the elasticity coefficients allow us to explain why these factors are critical by demonstrating their respective impacts on WY through sensitivity analysis. By separating the sensitivity of WY from the historical magnitude of driver changes, the elasticity analysis extends the scenario-based attribution and provides a basis for approximating future WY responses when projected changes in precipitation, vegetation, or [CO₂] are available.

3 Results

3.1 Changes in hydrometeorological factors

Figure 2 demonstrates the trends of annual precipitation, air temperature, relative humidity, atmospheric pressure, solar radiation, and NDVI across China during 1982–2017. Annual precipitation change exhibited a clear spatial distribution pattern, specifically decreases in central China, including the middle reaches of the Yellow River and the Yangtze River basins, and increases in the northwest and southeast. Air temperature exhibited a consistent warming trend across China. In contrast, relative humidity generally decreased across most China. Atmospheric pressure remained relatively stable. Regarding solar radiation, decreases were in northern China, while an increase was in southern regions. The decreasing solar radiation in northern China is likely due to increased aerosol concentrations (Liang et al., 2024). NDVI showed a significant increasing trend, which indicates an overall enhancement in vegetation growth across China. This trend was especially prominent in central and eastern regions, including the Yellow River Basin and the Yangtze River Basin. In these regions, LULC changes, such as afforestation and agricultural practices, likely contributed to the observed increases in NDVI (Chen et al., 2019).

Figure 2 Spatial patterns of trends in annual climatic and vegetation factors during 1982–2017. (a) Precipitation (mm yr⁻¹); (b) air temperature (°C yr⁻¹); (c) atmospheric pressure (Pa yr⁻¹); (d) shortwave radiation (W m⁻² yr⁻¹); (e) relative humidity (% yr⁻¹); (f) NDVI (yr⁻¹).

Significant changes in land use and land cover (LULC) occurred in China during 1982–2017, as illustrated in Fig. 3. Although the overall percentage distribution of major land cover types, namely grasslands, forests, croplands, and bare lands, remained relatively stable, these four categories dominated the landscape, with most changes concentrated within them. Notably, the transitions among these categories were characterized by mutual conversions, particularly from bare land to grasslands (Fig. 3). Spatially, the changes exhibited distinct regional patterns. In southern China, LULC changes were mainly characterized by the conversion of land to forests and grasslands. In contrast, the northeastern regions exhibited more complex transformations, with some areas shifting to bare land and croplands (Fig. 3).

Figure 3 Land use and land cover (LULC) changes from 1982–2017. (a) Spatial distribution of major conversion types among the four dominant LULC classes. F = forest, G = grassland, C = cropland, and B = bare land. The arrow indicates the direction of conversion from the original land-cover type to the converted type. (b) Percentage share of each major conversion type in the total converted area (%). (c) Chord diagram of the corresponding conversion flows, where the labels around the circle denote land-cover classes, and the surrounding axis values indicate converted area (km²).

3.2 Performance of the improved CCW model

As shown in Fig. 4a and b, the observed annual water yield (WY) and the simulated annual WY by the improved CCW model showed strong linear correlations (R² = 0.7), with the regression line slope being 1.45, R² being 0.7, and RMSE being 9.54 mm yr⁻¹. By contrast, the initial model without WUE showed weaker skill (slope = 1.45, R² = 0.68, RMSE = 9.62 mm yr⁻¹), indicating that explicitly representing [CO₂]-induced regulation of water-use efficiency measurably improves accuracy and reduces bias.

The estimated annual WY trends had distinct spatial patterns (Fig. 4c), which closely aligned with that of precipitation. Specifically, decrease trends in WY occurred in the central regions of the Yellow River Basin and the middle section of the Yangtze River Basin, while increase trends were found in other regions, with the southeast exhibiting the highest rate of increase.

Figure 4 (a) Validation of simulated WY trend using the improved CCW model; (b) validation of simulated WY trend using the original CCW model; (c) spatial distribution of WY trends under scenario S1 (actual situation) during 1982–2017.

3.3 Attribution analysis of annual WY changes

Figure 5 shows the distribution of WY changes caused by climate, vegetation, and [CO₂] changes, integrating both absolute magnitude (Fig. 5a, c, e) and relative dominance (Fig. 5b, d, f) of their contributions. Climate-driven WY changes exhibited marked spatial heterogeneity, with absolute increases exceeding 15 mm yr⁻¹ in southeastern China (Fig. 5a), corresponding to 60 %–90 % relative contributions (Fig. 5b). Central basins showed contrasting declines of 0–6 mm yr⁻¹ under climate forcing, while northeastern transitional zones displayed mixed positive/negative absolute changes (Fig. 5a) despite maintaining 40 %–70 % relative climate dominance (Fig. 5b). This spatial heterogeneity aligned with precipitation change patterns (Fig. 2a).

Figure 5 The absolute contributions of (a) climate, (c) vegetation, and (e) [CO₂], and the relative contributions of (b) climate, (d) vegetation, and (f) [CO₂] to changes in WY trends for 1982–2017.

Vegetation-mediated WY reductions reached 0–6 mm yr⁻¹ (Fig. 5c), accompanied by 0 %–60 % relative contributions (Fig. 5d). These effects originated from enhanced evapotranspiration through LULC changes and vegetation greening, particularly pronounced in central China. Specific regions in the Yangtze, Yellow, and northeastern rivers showed vegetation-driven relative contributions reaching 40 %–60 % (Fig. 5d). [CO₂] effects generated limited direct absolute impacts (< 5 mm yr⁻¹, Fig. 5e) but exerted 10 %–40 % relative influences (Fig. 5f) through stomatal closure mechanisms. This process partially counteracted vegetation-related WY decreases in transitional climates like northeastern China, where competing drivers created complex ecohydrological interactions (Fig. 5).

Figure 6 illustrated the spatial distribution of WY trend drivers over the past four decades. Climate change was the dominant factor of WY variation in more than 70 % regions, mainly in the Northwest, Southwest, Southeast, Pearl River basins, and other parts of the Yangtze and Yellow River basins. Vegetation changes ranked as the secondary control, dominating WY changes in parts of the Yangtze, Yellow, Songhua, Liao, and Hai Rivers. Remarkably, it was shown that the region where vegetation and [CO₂] had the dominant influence mainly distributes within precipitation ranges of 400–1600 mm. [CO₂]-induced effects were least influential at a national scale. This three-tiered hierarchy – climate changes as the primary forcing, vegetation changes as the secondary control, and [CO₂] effects as a localized modifier – reveals how hydrological regimes govern the spatial succession of dominant drivers across China's diverse ecohydrological gradients.

Figure 6 Spatial distributions of dominant factors controlling WY change. Driving factors include climate, vegetation, and [CO₂]. Climate: Areas where climate (e.g., precipitation, temperature) is the dominant factor influencing WY change; CO₂: Areas where [CO₂] is the primary driver of WY change; Vegetation: Areas where vegetation changes (e.g., NDVI, LULC) primarily drive WY changes. Climate-Vegetation: Areas where both climate and vegetation jointly influence WY; Climate-CO₂: Areas where both climate and [CO₂] jointly contribute to WY change; Vegetation-CO₂: Areas where vegetation changes and [CO₂] jointly control WY; Climate-Vegetation-CO₂: Areas where the combined effect of climate, vegetation, and [CO₂] jointly controls WY change. Additionally, the approximate isohyetal line shown in the figure were derived based on annual precipitation data from 1982–2017.

3.4 Elasticity of WY to main factors

The elasticity of WY to precipitation (εP), NDVI (εNDVI), and [CO₂] (εCO₂) exhibits distinct spatial patterns in Fig. 7. Nationally averaged elasticity coefficients showed that a 10 % increase in precipitation, [CO₂], and NDVI altered WY by 15.5 % (εP=1.55), 5.5 % (εCO₂=0.55), and −4.4 % (εNDVI = −0.44), respectively, indicating that, in terms of the sensitivity of runoff to changes in each factor, the ranking was precipitation > [CO₂] > NDVI. The historical relative changes in precipitation (0.48 % yr⁻¹), [CO₂] (0.49 % yr⁻¹), and NDVI (0.67 % yr⁻¹), suggest that precipitation increased WY by 27.2 %, [CO₂] increased WY by 9.97 %, and NDVI decreased WY by 10.7 % over 1982–2017. Therefore, although [CO₂] had a higher elasticity than NDVI, its historical contribution to WY change was slightly smaller than that of vegetation.

Figure 7 Spatial distribution of elasticity coefficients of WY relative to changes in hydrological factors such as (a) annual precipitation, (b) NDVI, and (c) [CO₂].

The elasticity coefficients of precipitation (εP), [CO₂] (εCO₂), and vegetation (|εNDVI|) all exhibited a coherent latitudinal decline across China's river basins, showing systematically higher sensitivity in northern regions than southern counterparts. Quantitatively, εP decreased from 2.09 in the Songhua River basin to 1.15 in the Southeastern Basin, accompanied by similar reductions in |εNDVI| (from 0.76 to 0.13) and εCO₂ (from 1.08 to 0.16) (Table 3).

A distinct abrupt transition zone in elasticity coefficients was identified around 33° N, closely aligning with China's traditional North-South physiographic divide. Around the zone, elasticity coefficients exhibited an abrupt decline from the Yellow River Basin to the Yangtze River Basin. Specifically, the Yellow River Basin showed higher sensitivities to precipitation (εP=1.87), [CO₂] (εCO₂=0.86), and NDVI (εNDVI = −0.53), which were approximately 1.4, 2.8, and 2.8 times greater, respectively, than those in the Yangtze River Basin (εP=1.31, εCO₂=0.31, εNDVI = −0.19).

Table 3 Elasticity Coefficients of Runoff to Precipitation, NDVI, and CO₂ in Different Watersheds.

Note: Some LULC types were excluded from the analysis. Due to many missing data points, the Liao River, Southwest, and Northwest river basins were also omitted.

Figure 8 Changes in different factors during 1982–2017: (a) relative change in annual precipitation, (b) relative change in NDVI, and (c) [CO₂] with its trend and relative annual change rate.

4 Discussion

4.1 Strength of the attribution analysis framework

This study developed an attribution framework based on the improved CCW model to separate the effects of climate change, vegetation structural change, and [CO₂]-induced physiological changes on WY changes. A key strength of this framework is that it explicitly links vegetation dynamics and [CO₂] effects to actual evapotranspiration at grid scale, rather than representing them indirectly through empirical catchment parameters or potential evapotranspiration adjustments.

First, the improved framework addresses a major limitation of Budyko approaches. Traditional Budyko approaches often attribute vegetation effects to temporal variations in the Budyko parameter n by either statistically regressing the Budyko parameter against vegetation proxies such as NDVI (Liu et al., 2024; Tan et al., 2023) or simplistically equating the Budyko parameter to vegetation effects (Li et al., 2020b; Zhou et al., 2023). In addition, some Budyko approaches account for [CO₂] effects by adjusting potential evapotranspiration (Liu et al., 2024). These treatments can conflate vegetation change and [CO₂] change with other catchment properties embedded in n (Zeng et al., 2020). Moreover, when [CO₂] effects are represented through PET adjustments, their influence on WY is transmitted through the Budyko function and further constrained by n, making it difficult to isolate [CO₂]-induced physiological effects from vegetation structural change and other land-surface controls. In contrast, our framework estimates actual ET through a carbon–water coupling structure and explicitly separates the effects of climate, LAI-related vegetation structural change, and [CO₂]-driven physiological changes, thereby providing a more direct basis for attributing WY changes than Budyko approaches.

Second, compared with the original CCW model using a static UWUE, our framework replaces it with a dynamic WUE formulation that accounts for long-term WUE changes associated with variations in [CO₂], VPD, and vegetation structure. Elevated [CO₂] reduces stomatal aperture, thereby lowering stomatal conductance and transpiration flux while only modestly increasing carbon assimilation, leading to an overall enhancement in WUE. This process is represented by the Medlyn-type stomatal conductance model (Medlyn et al., 2011), which links photosynthetic rate (A), transpiration (T), and vapor pressure deficit (D) as:

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{A}{T}}={\displaystyle \frac{C_{\mathrm{a}}{P}_{\mathrm{a}}}{\mathrm{1.6}\left(D+g_{\mathrm{1}}\sqrt{D}\right)}}\end{array}$$

where C_a is atmospheric CO₂ concentration, P_a is air pressure, D is vapor pressure deficit, and g₁ is an empirical slope parameter that quantifies plant sensitivity to [CO₂] and humidity. This formulation also indicates that WUE varies with [CO₂] and atmospheric water demand, rather than remaining constant over time. Incorporating this relationship into the CCW framework enables the model to distinguish NDVI-driven vegetation change effects on WY from [CO₂]-induced stomatal effects. Therefore, this improvement is not merely a statistical refinement, but also a process-based modification that helps explain why [CO₂] can partly offset vegetation-related WY reductions in some regions.

Third, the grid-scale implementation of the framework allows spatial heterogeneity in attribution results to be identified within basins. Many previous runoff attribution studies have been conducted at the basin scale (Liu et al., 2024, 2017; Yang et al., 2022), which is useful for summarizing overall basin responses but may mask substantial intra-basin differences. By conducting attribution at the grid scale, our framework captures spatial heterogeneity in WY controls within the same basin, distinguishing areas dominated by climatic variability from those more strongly affected by vegetation greening. This spatially explicit attribution provides more detailed and differentiated information for regional water resource assessment and ecological restoration planning within a basin.

4.2 New insights into attribution analysis

Our findings highlight climate change as the dominant driver of water yield (WY) changes (contributing > 70 %), which is consistent with previous assessments (Table 4). Climate impacts dominate in the Northwest and Southwest River Basins, as well as parts of the Yangtze, Yellow, Southeast, and Pearl River Basins. In contrast, vegetation and CO₂ dominate WY change in central China (portions of the Yangtze, Yellow, Songhua, Liao, and Hai River basin). The spatial pattern is slightly distinct from previous studies (Table 4), with noticeable discrepancies in some regions. These discrepancies may arise from differences in study periods, as vegetation effects have become more pronounced after 2000 in several regions but may be weakened when averaged over longer historical periods. This basin-scale contrast is also consistent with a broader hydroclimatic gradient, as vegetation and [CO₂] effects become relatively more important in the intermediate precipitation zone of 400–1600 mm. The relative importance of vegetation and [CO₂] is not randomly distributed but is concentrated along transitional hydroclimatic zones, where water availability is sufficient to support vegetation growth yet remains sensitive to changes in evapotranspiration.

Table 4 Comparative studies of the contribution of climate variability and vegetation to runoff changes.

Elasticity analysis (Sect. 3.4) revealed distinct sensitivities of WY to environmental drivers: precipitation exhibited the highest elasticity coefficient for the whole China (εP = 1.55), followed by [CO₂] (εCO₂=0.55) and NDVI (εNDVI = −0.44). However, spatial analysis showed that vegetation and [CO₂] collectively dominated WY changes in the 400–1600 mm precipitation zones, despite their lower sensitivity rankings. The contribution of each driver is determined by its elasticity multiplying the magnitude of its change. In the 400–1600 mm precipitation zones, NDVI displayed a larger relative temporal variation compared with precipitation (Fig. 8), which fluctuated within a narrower range. Therefore, high elasticity does not necessarily indicate high historical contribution. The larger relative change in NDVI enhanced the contribution of vegetation to WY change despite its lower elasticity, whereas the historical contribution of [CO₂] was limited by its slower accumulation rate (0.49 % yr⁻¹) despite its relatively high elasticity.

Although the historical contribution of [CO₂] was limited by its slow increase during 1982–2017, future [CO₂] may substantially amplify its effect on WY. Projections from the CMIP6 SSP585 scenario indicate [CO₂] will rise at 2.34 % yr⁻¹, which is nearly five times faster than the historical rate (Cheng et al., 2022). Using the historical elasticity coefficient as a simple sensitivity-based approximation, the increase of [CO₂] would drive a +1.29 % yr⁻¹ increase in annual WY by 2100. By comparison, future NDVI will increase at approximately 0.59 % yr⁻¹ under CMIP6 SSP585 (Li et al., 2025), its effect on WY would be only −0.26 % yr⁻¹. Therefore, hydrological effect of rising [CO₂] may be much larger than the direct effect of future vegetation greening and therefore should not be ignored in future WY assessments. However, this estimate assumes that the historical [CO₂] elasticity remains unchanged and does not account for possible counteracting effects, such as rising VPD, drought stress, or saturation of WUE under elevated [CO₂] (Adams et al., 2020; Li et al., 2023).

From a policy perspective, these spatial contrasts have distinct implications for regional water management. In regions primarily affected by vegetation change, such as the Yangtze and Yellow River basins, priority should be given to optimizing vegetation composition and preventing overgreening that may reduce runoff. Conversely, in regions primarily affected by climate change, such as Northwest and Southeast China, adaptive measures that emphasize precipitation variability, water storage capacity, and drought resilience are crucial. Tailoring water management strategies to these distinct regional conditions can enhance the effectiveness of both ecological restoration and climate adaptation programs across China.

4.3 Uncertainty in attribution analysis

This study used an improved CCW model to attribute annual WY changes in China to climate change, vegetation structural change, and [CO₂]-induced physiological effects, providing a more mechanistic understanding of vegetation and [CO₂] effects of hydrological responses. However, several limitations remain that need to be addressed in future work to improve the accuracy and robustness of the results.

First, the adopted WUE formulation is defined at the ecosystem scale, but the present implementation simplifies the interception term by assuming f_i=0. Soil evaporation is still implicitly represented through the LAI-based partitioning term, whereas canopy interception is not explicitly allowed to vary across space or time. This simplification may bias the absolute magnitude of ET and WY, particularly in humid and forested regions, although its influence on long-term WY trends is expected to be small. Nevertheless, given the small interannual variation of f_i (Zhao et al., 2022), its influence on runoff trends is negligible in our study (Cheng et al., 2017). Future work should prioritize its explicit calculation.

Second, the attribution framework treats climate, vegetation, and [CO₂] as separate drivers, whereas these factors interact with each other in real ecohydrological systems. Vegetation exhibits interactions and feedbacks with climate, making it difficult to strictly separate the impacts of climate change, vegetation dynamics, and [CO₂] on hydrological responses. Changes in vegetation, such as NDVI, reflect a combination of climate change, human activities (e.g., reforestation and irrigation), and natural vegetation growth. Additionally, vegetation greening in upwind regions can increase atmospheric moisture, potentially enhancing precipitation downwind (Zhang et al., 2021a), which may counteract some of the negative impacts of increased evapotranspiration on local WY. Although the climate data used in our model may implicitly capture some of these feedbacks, they cannot be explicitly separated in this analysis. This limitation also affects the simple elasticity-based SSP585 estimate of [CO₂] effects, because the estimate does not account for concurrent changes in VPD, drought stress, or possible saturation of WUE under elevated [CO₂]. Consequently, our results represent an attempt to estimate the direct first-order net impacts of climate, vegetation greening, and [CO₂] increase on WY (Zhang et al., 2021b). Future research should adopt more comprehensive models that consider soil-vegetation-atmosphere interactions to better differentiate the contributions of each driving factor to WY.

Third, the improved CCW model does not incorporate certain human activities, such as large-scale irrigation, groundwater pumping, and reservoir regulation, which also limits the assumption that WY = P − ET with negligible storage change. This assumption is mainly appropriate at the annual scale and becomes less reliable in regions where long-term groundwater depletion or reservoir operation substantially alters water storage and runoff routing. For instance, irrigation can sustain vegetation greening during dry seasons, potentially amplifying the vegetation–climate feedback on water yield. Incorporating such anthropogenic processes into the CCW framework through coupled irrigation and water management modules would enable more comprehensive attribution analyses in future studies. Our research also excludes water bodies and built-up land. While urbanization can increase flood risks due to the growing proportion of impervious surfaces (Wasko and Sharma, 2017), these land-use changes represent a small fraction of China's land area.

5 Conclusions

In this study, we improved the coupled carbon and water (CCW) model incorporating a dynamic water use efficiency (WUE) calculation to explicitly represent the [CO₂] physiological feedback on water yield (WY). This improvement enabled a comprehensive national-scale assessment quantifying the relative contributions of climate forcing, vegetation structural changes, and [CO₂]-driven stomatal changes to WY dynamics in China. The main conclusions are as follows:

1. The improved CCW model effectively simulated WY variations in most basins under elevated [CO₂] scenarios, demonstrating its applicability and reliability for modeling WY changes.

2. Climate change, particularly variations in precipitation, emerged as the primary driver influencing WY, with significant regional disparities in its effects. Vegetation change constituted the second most critical factor, predominantly resulting in WY reduction, notably in central China, whereas the effect of [CO₂]-induced stomatal closure on WY was comparatively minor. Spatial analysis aligned with isohyetal lines further revealed that vegetation change and [CO₂] exerted greater influence within the 400–1600 mm precipitation range.

3. The elasticity analysis of WY indicated that northern basins exhibit higher sensitivity to influencing factors, whereas southern basins demonstrate relatively lower elasticity. Specifically, the absolute elasticity coefficients for the whole China were ranked in descending order as follows: precipitation > [CO₂] > NDVI. Thus, under SSP585, rising [CO₂] may have a stronger effect on future WY than vegetation greening, with an estimated +1.29 % yr⁻¹ WY effect compared with −0.26 % yr⁻¹ from projected NDVI increases, although VPD changes and possible WUE saturation are not considered.

These insights provide a nuanced understanding of regional hydrological responses, which is essential for sustainable water resource management under changing environment.

Appendix A: Original CCW Model

The Coupled Carbon and Water (CCW) model is a data-driven, remote sensing-based model that estimates gross primary productivity (GPP) and evapotranspiration (ET) through a coupled carbon–water formulation (Zhang et al., 2016b). In the original CCW model, GPP is estimated using a light-use efficiency approach:

$$\begin{array}{}\text{(A1)}& \mathrm{GPP}=\mathrm{APAR}\times \mathit{\epsilon }=\mathrm{PAR}\times \mathrm{FPAR}\times {\mathit{\epsilon }}_{\mathrm{pot}}\times R_{\mathrm{s}}\times T_{\mathrm{s}}\times W_{\mathrm{s}}\end{array}$$

where APAR is the absorbed photosynthetically active radiation (MJ m⁻²), which is calculated as the product of incident photosynthetically active radiation (PAR) and the fraction of PAR absorbed by vegetation (FPAR). PAR is typically assumed to be 45 % of the total shortwave radiation (Running et al., 2000), and FPAR is determined by the normalized difference vegetation index (NDVI) (Sims et al., 2005). ε is the realized light-use efficiency (g C MJ⁻¹), which is calculated by multiplying the potential light-use efficiency (ε_pot) and environmental scalars for diffuse radiation (R_s), temperature (T_s), and moisture stress (W_s).

ET is then derived from GPP using the underlying water-use efficiency (UWUE):

$$\begin{array}{}\text{(A2)}& \mathrm{ET}={\displaystyle \frac{\mathrm{GPP}\times {\mathrm{VPD}}^{\mathrm{0.5}}}{\mathrm{UWUE}}}\end{array}$$

where VPD is vapor pressure deficit and UWUE is the underlying water-use efficiency. The UWUE parameter is prescribed by vegetation type and is calibrated using FLUXNET observations (Zhou et al., 2014; Zhang et al., 2016b). Therefore, UWUE is treated as static in the original CCW framework and does not explicitly vary with rising [CO₂], atmospheric pressure, vegetation structure, or long-term physiological adjustment.