Future response of ecosystem water use efficiency to CO₂ effects in the Yellow River Basin, China

Ecosystem water use efficiency (WUE) is pivotal for understanding carbon–water cycle interplay. Current research seldom addresses how WUE might change under future elevated CO₂ concentrations, limiting our understanding of regional ecohydrological effects. We present a land–atmosphere attribution framework for WUE in the Yellow River basin (YRB), integrating the Budyko model with global climate models (GCMs) to quantify the impacts of climate and underlying surface changes induced by CO₂. Additionally, we further quantitatively decoupled the direct and secondary impacts of CO₂ radiative and biogeochemical effects. Attribution results indicate that WUE in the YRB is projected to increase by 0.36–0.84 gC kg⁻¹H₂O in the future, with climate change being the predominant factor (relative contribution rate of 77.9 %–101.4 %). However, as carbon emissions intensify, the relative importance of land surface changes becomes increasingly important (respective contribution rates of −1.4 %, 14.9 %, 16.9 %, and 22.1 % in SSP126, SSP245, SSP370, and SSP585). Typically, WUE is considered a reflection of an ecosystem's adaptability to water stress. Thus, we analyzed the response of WUE under different scenarios and periods and various drought conditions. The results show a distinct “two-stage” response pattern of WUE to drought in the YRB, where WUE increases under moderate–severe drought conditions but decreases as drought intensifies across most areas. Furthermore, GCM projections suggest that plant adaptability to water stress may improve under higher-carbon-emission scenarios. Our findings enhance the understanding of regional ecohydrological processes and provide insights for future predictions of drought impacts on terrestrial ecosystems.

1 Introduction

Ecosystem water use efficiency (WUE) is commonly defined as the ratio of the ecosystem's gross primary productivity (GPP) to water evapotranspiration (ET), i.e., WUE = $\mathrm{GPP}/\mathrm{ET}$, which reflects the carbon gain per unit of water lost (Keenan et al., 2013; Li et al., 2023a; Naeem et al., 2023). This metric couples the water and carbon cycle processes, serving as a key characteristic indicator for ecosystem function (Liu et al., 2020). Investigating the spatiotemporal characteristics and driving factors of WUE holds substantial practical importance for exploring the response mechanisms of ecosystems under changing environmental conditions (Huang et al., 2017; Tan et al., 2023; Zhou et al., 2017).

The increase in atmospheric CO₂ concentrations brings profound impacts on regional WUE, GPP, and ET through radiative effects and biogeochemical effects (Naeem et al., 2023; Yang et al., 2022). The interactions between land and atmosphere further complicate the changes associated with rising CO₂ levels (De Kauwe et al., 2021; He et al., 2023; Zhan et al., 2022). Regarding radiative effects, elevated CO₂ directly affects hydro-atmospheric processes, altering global and regional precipitation and evaporation patterns (Bintanja and Andry, 2017; L. Gu et al., 2023; Yin et al., 2018). The biogeochemical effects are primarily manifested through the impacts of CO₂ fertilization and changes in plant stomatal conductance on the underlying vegetation structure. Fertilization effect refers to the increased CO₂ concentrations enhancing photosynthesis rates, thereby increasing vegetation productivity and promoting plant growth (Chen et al., 2024; Sun et al., 2018). He et al. (2023) demonstrated that CO₂ was the dominant driver of the increase in forest carbon sinks over the past few decades. The growth in plant biomass can lead to an increase in evapotranspiration (Xie et al., 2020). However, elevated CO₂ concentrations can also contribute to a reduction in leaf stomatal conductance, potentially resulting in a decrease in transpiration (De Kauwe et al., 2021; Zhu et al., 2011). Beyond these direct climatic changes through radiative and direct biogeochemical effects on the land surface, the climate and land surface are also subject to secondary effects. For example, climate change can further impact vegetation cover (Berg et al., 2017); changes in the structure of surface vegetation will alter the water–energy exchange process, thereby affecting precipitation and other processes (Zhou et al., 2022). Few studies have focused on the relative contributions of CO₂'s direct and secondary effects on WUE. However, understanding these contributions is crucial for exploring the mechanisms of WUE changes and accurately predicting future changes in WUE.

Conducting a joint analysis of hydrological variables using conceptual models and global climate models (GCMs) serves as an effective method for attributing land–atmosphere processes (Zhou et al., 2023). Since its inception, the Budyko model (Budyko, 1974) has been extensively employed to evaluate the relationships between water–energy exchanges and changes in surface characteristics under different climate scenarios (Fang et al., 2020; Fathi et al., 2019; Roderick and Farquhar, 2011). As a conceptual framework, the Budyko model is characterized by its simplicity and excellent performance, elucidating the relationship between $\mathrm{ET}/P$ ($\mathrm{evapotranspiration}/\mathrm{precipitation}$) and $\mathrm{PET}/P$ ($\mathrm{potentialevapotranspiration}/\mathrm{precipitation}$) (Cheng et al., 2011; Choudhury, 1999; Xu et al., 2022). Yang et al. (2015) applied the Budyko framework to the carbon cycle, revealing the relationship between GPP and energy exchange. Fang et al. (2020) explored the relationship between GPP, WUE, and water–energy relations in semi-arid basins through the Budyko model, investigating the correlations between WUE and underlying surface parameters. Despite previous efforts to extend the scope of the Budyko model from water–energy relations to water–carbon–energy relations, there have been few descriptions of how these relationships may change in the future along with the relative contribution of elements to the studied variables.

The complexity of drought events can have significant impacts on the regional economy and environment (Yin et al., 2023; Yuan et al., 2023). Variations in WUE under drought conditions have been recognized as an indicator of ecosystems' adaptability to drought stress (Huang et al., 2017; Ponce-Campos et al., 2013). Ponce-Campos et al. (2013) suggested that ecosystems increased their WUE in response to water stress. However, observations have shown considerable variability in WUE responses to drought events across different regions (Huang et al., 2017; Lu and Zhuang, 2010; M. Wang et al., 2021; Xie et al., 2016; Yang et al., 2021). This indicates the need for more in-depth research into the WUE–drought relationship. According to Pokhrel et al. (2021), in the context of climate change, the risk of drought in most regions worldwide is expected to continue increasing, which raises the question of how WUE will respond to drought in the future. This study aims to explore the future response mechanisms of WUE to drought to better assess the resilience of ecosystems under climate change. The drought severity index (TWSA-DSI) based on the terrestrial water storage anomaly (TWSA) has been widely used in drought monitoring and analysis (Yin et al., 2022, 2023; Zhao et al., 2017). This index effectively monitors regional drought characteristics and development paths (Yin et al., 2022). Compared to other traditional drought indicators such as the Palmer drought severity index (PDSI), soil moisture drought index (SMI), standardized precipitation index (SPI), standardized runoff index (SRI), etc., TWSA-DSI is better able to detect overall regional drought conditions and to explore the effects of drought on hydrological systems and vegetation growth (A et al., 2017; Du et al., 2019; Zhao et al., 2017). This study chooses to use TWSA-DSI to identify regional drought events, allowing for a more nuanced understanding of how WUE may adapt or respond to these conditions in various ecosystems.

The Yellow River basin (YRB) is an important basin in China, situated in the arid and semi-arid regions. The ecosystem within the basin is complex and vulnerable to drought (Huang et al., 2015). Changes in the ecohydrological processes of the YRB in the future will profoundly impact China's socioeconomic development. This paper selects the YRB as the study area, trying to deepen the understanding of the spatiotemporal variation mechanisms of WUE in this region. This study aims to answer the following three research questions over both historical and future timescales:

• a.

  What changes will occur in the spatiotemporal distribution of WUE within the basin?

• b.

  How can the impact of elevated atmospheric CO₂ concentrations on WUE changes be quantitatively explained?

• c.

  How does the response mechanism of WUE to drought within the basin evolve under future scenarios?

The reminder of this paper is organized as follows. Section 2 describes the study area and the data used. Section 3 defines the methodology. Section 4 and Sect. 5 give the results and discussion of the study, respectively. Finally, Sect. 6 concludes the paper.

2 Study area and data

2.1 Study area

The Yellow River is the second-longest river in China and the fifth-longest river in the world, with a total length of about 5464 km. The YRB (32–41° N, 95–119° E) is the second-largest drainage basin in China (Fig. 1), with an approximate area of 795 000 km² (including a 42 000 km² inland flow area). The YRB spans the Tibetan Plateau, the Inner Mongolian Plateau, the Loess Plateau, and the North China Plain. The average altitude of the basin exceeds 3000 m, and the terrain becomes progressively lower from west to east. The sub-basins are characterized by a wide range of climatic and geomorphologic variations. Overall, the basin has a largely semi-humid climate in the southeast and an arid/semi-arid climate in the center and the northwest. The YRB is one of the earliest origins of human civilization, and the Yellow River is also known as the “mother river” of China. The YRB has a population of more than 120 million people, a gross national product of about 8 trillion RMB (≈ 1.1 trillion USD), and a cultivated agricultural area of about 200 000 km². The YRB has been regarded as one of the most important basins in China because of its important role in socioeconomic development, ecological resource protection, and agricultural food production.

Figure 1 Location and topography of the YRB.

2.2 Data collection

2.2.1 Climate and GPP data

The monthly precipitation, potential evapotranspiration, and evapotranspiration and GPP datasets during 1997–2016 were both retrieved from National Tibetan Plateau Data Center (https://data.tpdc.ac.cn/, last access: 1 February 2024). The precipitation and potential evapotranspiration dataset are based on the Climatic Research Unit (CRU) time series (TS) v4.02 dataset and WorldClim v2.0 data to generate data with high spatial resolution and accuracy (Ding and Peng, 2020, 2021). The evapotranspiration dataset employed a calibration-free nonlinear complementary relationship model, and the evaluation suggests that the evapotranspiration dataset has a good performance (Ma et al., 2019a; Ma and Szilagyi, 2019). The GPP dataset was developed based on the satellite-based near-infrared reflectance of vegetation (NIRv) (S. Wang et al., 2021). The accuracy and reliability have been validated by many studies (Cai et al., 2023; Feng et al., 2023; P. Li et al., 2023; Peng et al., 2023; Xing et al., 2023). For consistency, the above data were bilinearly interpolated to 0.25°×0.25°. To avoid uncertainties introduced by extremely small values, this study calculated only grid-point monthly WUE with GPP > 10 gC m² and ET > 10 mm following the approach used by Naeem et al. (2023).

2.2.2 Vegetation indicator

In this study, the normalized difference vegetation index (NDVI) and the leaf area index (LAI) were selected as vegetation indicators to characterize the underlying surface conditions. The GIMMS3g NDVI and GIMMS4g LAI data, acquired from the Advanced Very High Resolution Radiometer (AVHRR) sensor from NOAA's series of weather satellites, provide high-spatiotemporal-resolution and high-accuracy observed vegetation indicators (Pinzon and Tucker, 2014; Cao et al., 2023a). These data are freely available on NASA's website (https://daac.ornl.gov/VEGETATION/guides/Global_Veg_Greenness_GIMMS_3G.html, last access: 1 February 2024; https://doi.org/10.5281/zenodo.8281930, Cao et al., 2023b). To match the other data, the vegetation data were processed into monthly 0.25° × 0.25° gridded data.

2.2.3 TWSA datasets

To assess long-term changes in drought conditions in the basin, monthly reconstructed TWSA datasets from 1997–2016 were used in this study (Humphrey and Gudmundsson, 2019a). This dataset is based on two different Gravity Recovery and Climate Experiment (GRACE)/Gravity Recovery and Climate Experiment Follow-On (GRACE-FO) satellite products and three different meteorological forcing datasets, producing six reconstructed TWSA datasets with 100 ensemble members each. Its performance is effectively validated (H. Gu et al., 2023; Yin et al., 2023; Zhong et al., 2023). In this study, we calculated the mean results of the six datasets and thus bilinearly interpolated them to 0.25° × 0.25°.

2.2.4 GCM output data

In order to project future climate scenarios, we selected seven GCMs in CMIP6 (Table 1). These models were chosen because they encompass most or even all of the variables required for our analysis. Furthermore, they incorporate the experimental settings needed for our research and integrate the dynamic response of vegetation to climate change. In comparison to CMIP5, CMIP6 uses the matrix framework of the Shared Socioeconomic Pathway (SSP) and the Representative Concentration Pathway (RCP). In this study, we utilized simulated monthly data for GPP (gpp), evapotranspiration (evspsbl), precipitation (pr), and terrestrial water storage (mrtws), where the terms in parentheses represent the variable names within CMIP6. These data encompassed the historical period from 1985–2014 and projected future scenarios, including SSP126 (sustainable development; 2015–2100), SSP245 (moderate development; 2015–2100), SSP370 (regional rivalry development; 2015–2100), and SSP585 (conventional development; 2015–2100). These scenarios represent a range of possible futures from sustainable, low-emission trajectories to high-emission pathways based on varying degrees of socioeconomic development and environmental policy implementation. To assess the impact of rising CO₂ concentration on WUE, we also used data from three CO₂ sensitivity experiments – namely, 1ptCO2, 1ptCO2-bgc, and 1ptCO2-rad. In 1ptCO2, the concentration of CO₂ in the atmosphere increases by 1 % per year. In 1ptCO2-bgc and 1ptCO2-rad, the growth rate of CO₂ is the same as that in 1ptCO2, but the ones ending in -rad is only coupled with the atmospheric part (CO₂ radiative forcing), while that ending in -bgc is only coupled with the land part (CO₂ biogeochemically forcing), and the other part maintained the fixed CO₂ concentration level of pre-industrialization. It is notable that in all three experiments, the CO₂ concentration rises at the same rate as in SSP585.

Table 1 CMIP6 GCMs used in this study. The eighth column shows the model outputs for the three CO₂ sensitivity experiments (1ptCO2, 1ptCO2-bgc, and 1ptCO2-rad). GPP, ET, PET, P, and TWSA are represented by red, blue, yellow, green, and black crosses, respectively.

To obtain the TWSA series from CMIP6, the common practice is to select the same baseline period (2004–2009) as GRACE/GRACE-FO satellite data and deduct the series mean value during this period to obtain the TWSA series (Pokhrel et al., 2021; Yin et al., 2022). In this study, we also used this method. Since there is no direct output of potential evapotranspiration from CMIP6, we used near-surface air temperature (Tas), near-surface wind speed (sfcWind), relative humidity (hur), surface pressure (ps), heat flux, net surface radiation, and other data to calculate potential evapotranspiration with the Penman–Monteith equation (Eq. 1), which is the only standard method proposed by the Food and Agriculture Organization of the United Nations (FAO) (Allan et al., 1998). All GCM outputs were bilinearly interpolated to a common spatial resolution of 0.25° × 0.25°.

$$\begin{array}{}\text{(1)}& \mathrm{PET}={\displaystyle \frac{\mathrm{0.408}\mathrm{\Delta }\left(R_{\mathrm{n}}-G\right)+\mathit{\gamma }\frac{\mathrm{900}}{T+\mathrm{273}}{\mathit{\mu }}_{\mathrm{2}}\left(e_s-e_a\right)}{\mathrm{\Delta }+\mathit{\gamma }\left(\mathrm{1}+\mathrm{0.34}{\mathit{\mu }}_{\mathrm{2}}\right)}},\end{array}$$

where PET is the potential evapotranspiration (mm), Δ is the slope of the saturation vapor pressure curve at temperature T (kPa °C⁻¹), R_n is net surface radiation (MJ m⁻²), G is the soil heat flux (MJ m⁻²), γ is the psychrometric constant (kPa °C⁻¹), T is the near-surface air temperature (°C), μ₂ is the near-surface wind speed (m s⁻¹), and (e_s−e_a) is the vapor pressure deficit (kPa).

3 Methodology

3.1 Trend-preserving bias correction

Bias correction is the adjustment of GCM simulation data to reduce their systematic deviation as compared to observations. Compared to traditional bias correction methods, the trend-preserving bias correction method enables more robust bias adjustments for extremes and more accurately maintains trends across quantiles. The main steps of the trend-preserving bias correction are as below (Lange, 2019).

• Step 1. Detrend observation data ($x_{\mathrm{his}}^{\mathrm{obs}})$, GCM outputs during historical period ($x_{\mathrm{his}}^{\mathrm{sim}})$, GCM outputs during the future period or CO₂ experiment period ($x_{\mathrm{fut}}^{\mathrm{sim}})$ to obtain three new series, $x_{\mathrm{his}}^{\mathrm{Dobs}}$, $x_{\mathrm{his}}^{\mathrm{Dsim}}$ and $x_{\mathrm{fut}}^{\mathrm{Dsim}}$, respectively.

• Step 2. Transfer the change signals between the simulation series ($x_{\mathrm{his}}^{\mathrm{Dsim}}$ and $x_{\mathrm{fut}}^{\mathrm{Dsim}})$ to the observation series ($x_{\mathrm{his}}^{\mathrm{Dobs}})$ to obtain the pseudo-observation series, $x_{\mathrm{fut}}^{\mathrm{Dobs}}$.

• Step 3. Adjust the distribution of $x_{\mathrm{fut}}^{\mathrm{Dsim}}$ based on $x_{\mathrm{fut}}^{\mathrm{Dobs}}$ by quantile mapping. The formula is as follows:

  $$x_{\mathrm{fut}}^{\mathrm{Dobs}}=x_{\mathrm{his}}^{\mathrm{Dobs}}+(x_{\mathrm{fut}}^{\mathrm{Dsim}}-x_{\mathrm{his}}^{\mathrm{Dsim}}).$$

• Step 4. Add the results from Step 3 to the trend values previously subtracted to get the corrected results.

This study took into account the variations in the frequency distribution functions of TWSA across different periods, thereby implementing corrections for individual months separately. The results indicate that this month-specific correction approach yields better outcomes than a correction approach that does not differentiate by month.

3.2 WUE attribution framework

Many studies applied the Budyko model to isolate the effects of climate change and land surface properties on evapotranspiration or runoff at the basin/catchment scale (Roderick and Farquhar, 2011; Yang and Yang, 2011; Xu et al., 2022). The general Budyko model contemplates the interconnections and feedbacks between water and energy cycles while also considering the impact of basin characteristics. The expression of the model is as follows:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{E}{P}}=f\left({\displaystyle \frac{\mathrm{PET}}{P}},c\right),\end{array}$$

where E and P, respectively, represent evapotranspiration (mm) and precipitation (mm) and c is the basin underlying surface parameter.

In this study, we analyzed the relationship between hydrological elements and vegetation structure. After a selection process, we employed the linear Budyko function structure for WUE as proposed in previous studies (Cheng et al., 2011; Fang et al., 2020). The formula is expressed as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{WUE}}{P}}=f\left({\displaystyle \frac{\mathrm{PET}}{P}},m,n\right)=m{\displaystyle \frac{\mathrm{PET}}{P}}+n,\end{array}$$

where WUE is the water use efficiency (g C kg⁻¹ H₂O) and m and n represent the underlying surface conditions of the basin. The Spearman correlation coefficient and its significance were used to study the relationships between variables in Eq. (3). The model regression results and the coefficient of determination (R²) were derived from data analysis to evaluate the accuracy of the model.

To avoid the effects of climate variability and seasonal changes in the basin underlying parameters in the Budyko model, this study focused on a timescale of about 30 years (Ning et al., 2019).

In traditional attribution analysis based on the Budyko model, the total differentiation method is widely used. Therefore, based on Eq. (3), it can be formulated as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle }\begin{array}{rl}\mathrm{dWUE}& ={\displaystyle \frac{\partial \mathrm{WUE}}{\partial P}}\mathrm{d}P+{\displaystyle \frac{\partial \mathrm{WUE}}{\partial \mathrm{PET}}}\mathrm{dPET}+{\displaystyle \frac{\partial \mathrm{WUE}}{\partial m}}\mathrm{d}m\\ & +{\displaystyle \frac{\partial \mathrm{WUE}}{\partial n}}\mathrm{d}n,\end{array}\text{(5)}& {\displaystyle \frac{\partial \mathrm{WUE}}{\partial P}}=f\left({\displaystyle \frac{\mathrm{PET}}{P}},m,n\right)-{\displaystyle \frac{\mathrm{PET}}{P}}{\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}},\text{(6)}& {\displaystyle \frac{\partial \mathrm{WUE}}{\partial \mathrm{PET}}}={\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}},\text{(7)}& {\displaystyle \frac{\partial \mathrm{WUE}}{\partial m}}+{\displaystyle \frac{\partial \mathrm{WUE}}{\partial n}}=P{\displaystyle \frac{\partial f}{\partial m}}+P{\displaystyle \frac{\partial f}{\partial n}}.\end{array}$$

In Eqs. (5) and (6), the terms on the left side of the equality represent the impact of climate change on WUE. Meanwhile, in Eq. (7), the term on the left represents the influence of land surface changes. In accordance with Zhou et al. (2016), we assume that there is such an equation that holds: $\frac{\partial \mathrm{WUE}}{\partial m}+\frac{\partial \mathrm{WUE}}{\partial n}=P\mathrm{d}\left(\frac{\partial \mathrm{WUE}}{\partial P}\right)+\mathrm{PETd}\left(\frac{\partial \mathrm{WUE}}{\partial \mathrm{PET}}\right)$. In other words, the impact of land surface changes can be represented by the calculated result of $P\mathrm{d}\left(\frac{\partial \mathrm{WUE}}{\partial P}\right)+\mathrm{PETd}\left(\frac{\partial \mathrm{WUE}}{\partial \mathrm{PET}}\right)$. The proof proceeds as follows:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& P\mathrm{d}{\displaystyle \frac{\partial \mathrm{WUE}}{\partial P}}+\mathrm{PETd}{\displaystyle \frac{\partial \mathrm{WUE}}{\partial \mathrm{PET}}}=P\mathrm{d}\left[f-{\displaystyle \frac{\mathrm{PET}}{P}}{\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}}\right]\\ & +\mathrm{PETd}\left[{\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}}\right]\\ & =P\left({\displaystyle \frac{\partial f}{\partial \mathrm{PET}}}\mathrm{dPET}+{\displaystyle \frac{\partial f}{\partial P}}\mathrm{d}P+{\displaystyle \frac{\partial f}{\partial m}}\mathrm{d}m+{\displaystyle \frac{\partial f}{\partial n}}\mathrm{d}n\right)\\ & -P\mathrm{d}\left[{\displaystyle \frac{\mathrm{PET}}{P}}{\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}}\right]+\mathrm{PETd}\left[{\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}}\right]\\ & =P\left[{\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}}{\displaystyle \frac{\partial \left(\frac{\mathrm{PET}}{P}\right)}{\partial \mathrm{PET}}}\mathrm{dPET}\right.{\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}}{\displaystyle \frac{\partial \left(\frac{\mathrm{PET}}{P}\right)}{\partial P}}\mathrm{d}P+\\ & \left.+{\displaystyle \frac{\partial f}{\partial m}}\mathrm{d}m+{\displaystyle \frac{\partial f}{\partial n}}\mathrm{d}n\right]-P\mathrm{d}\left[{\displaystyle \frac{\mathrm{PET}}{P}}{\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}}\right]+\mathrm{PETd}\left[{\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}}\right]\\ & =P\left[{\displaystyle \frac{\mathrm{1}}{P}}{\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}}\mathrm{dPET}-{\displaystyle \frac{\mathrm{PET}}{P^{\mathrm{2}}}}{\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}}\mathrm{d}P\right]-P\mathrm{d}\left[{\displaystyle \frac{\mathrm{PET}}{P}}{\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}}\right]\\ & +\mathrm{PETd}\left[{\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}}\right]+P{\displaystyle \frac{\partial f}{\partial m}}\mathrm{d}m+P{\displaystyle \frac{\partial f}{\partial n}}\mathrm{d}n\\ & ={\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}}\mathrm{dPET}+\mathrm{PETd}\left[{\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}}\right]-{\displaystyle \frac{\mathrm{PET}}{P}}{\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}}\mathrm{d}P\\ & -P\mathrm{d}\left[{\displaystyle \frac{\mathrm{PET}}{P}}{\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}}\right]+P{\displaystyle \frac{\partial f}{\partial m}}\mathrm{d}m+P{\displaystyle \frac{\partial f}{\partial n}}\mathrm{d}n\\ & =\mathrm{dPET}{\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}}-\mathrm{d}P{\displaystyle \frac{\mathrm{PET}}{P}}{\displaystyle \frac{\partial f}{\partial \left(\frac{\mathrm{PET}}{P}\right)}}+P{\displaystyle \frac{\partial f}{\partial m}}\mathrm{d}m+P{\displaystyle \frac{\partial f}{\partial n}}\mathrm{d}n\\ & =P{\displaystyle \frac{\partial f}{\partial m}}\mathrm{d}m+P{\displaystyle \frac{\partial f}{\partial n}}\mathrm{d}n.\end{array}\end{array}$$

In practical applications, Eq. (4) is commonly expressed as the first-order Taylor approximation for attribution analysis:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}\mathrm{\Delta }\mathrm{WUE}& ={\displaystyle \frac{\partial \mathrm{WUE}}{\partial P}}\mathrm{\Delta }P+{\displaystyle \frac{\partial \mathrm{WUE}}{\partial \mathrm{PET}}}\mathrm{\Delta }\mathrm{PET}+{\displaystyle \frac{\partial \mathrm{WUE}}{\partial m}}\mathrm{\Delta }m\\ & +{\displaystyle \frac{\partial \mathrm{WUE}}{\partial n}}\mathrm{\Delta }n+O,\end{array}\end{array}$$

where ΔWUE, ΔP, ΔPET, Δm, and Δn, respectively, represent the actual changes in WUE; precipitation; potential evapotranspiration; and two underlying surface parameters, i.e., m and n. It is noteworthy that in Eq. (9), there is consistently a residual term, O, which represents the higher-order terms in the Taylor expansion. The value of O will affect the accuracy of the attribution results.

Zhou et al. (2015) have proven that if the function structure between a dependent variable Z and N independent variables (X₁, X₂, …, X_N) is

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{Z}{X_{\mathrm{1}}}}=f\left({\displaystyle \frac{X_{\mathrm{2}}}{X_{\mathrm{1}}}}{\displaystyle \frac{X_{\mathrm{3}}}{X_{\mathrm{1}}}}{\displaystyle \frac{X_N}{X_{\mathrm{1}}}}\right),\end{array}$$

the sum of the elasticity coefficients is then unity, i.e. the following:

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\partial Z/Z}{\partial X_{\mathrm{1}}/X_{\mathrm{1}}}}+{\displaystyle \frac{\partial Z/Z}{\partial X_{\mathrm{2}}/X_{\mathrm{2}}}}+\mathrm{\dots }+{\displaystyle \frac{\partial Z/Z}{\partial X_N/X_N}}=\mathrm{1}.\end{array}$$

Demonstrating that a similar complementary relationship exists in $\frac{\mathrm{WUE}}{P}=f\left(\frac{\mathrm{PET}}{P}mn\right)$ is straightforward; namely, the following applies:

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\partial \mathrm{WUE}/\mathrm{WUE}}{\partial P/P}}+{\displaystyle \frac{\partial \mathrm{WUE}/\mathrm{WUE}}{\partial \mathrm{PET}/\mathrm{PET}}}=\mathrm{1}.\end{array}$$

Therefore, according to previous studies (Zhou et al., 2016; Yang et al., 2023; Zhang et al., 2023), if such a complementary relationship (Eq. 12) exists, the difference in dependent variable Z (ΔWUE) can be expressed in the following form without any residual term:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}\mathrm{\Delta }\mathrm{WUE}& =\mathit{\alpha }\left[{\left({\displaystyle \frac{\partial \mathrm{WUE}}{\partial P}}\right)}_{\mathrm{1}}\mathrm{\Delta }P+{\left({\displaystyle \frac{\partial \mathrm{WUE}}{\partial \mathrm{PET}}}\right)}_{\mathrm{1}}\mathrm{\Delta }\mathrm{PET}\right.\\ & \left.+P_{\mathrm{2}}\mathrm{\Delta }\left({\displaystyle \frac{\partial \mathrm{WUE}}{\partial P}}\right)+{\mathrm{PET}}_{\mathrm{2}}\mathrm{\Delta }\left({\displaystyle \frac{\partial \mathrm{WUE}}{\partial \mathrm{PET}}}\right)\right]\\ & +\left(\mathrm{1}-\mathit{\alpha }\right)\left[{\left({\displaystyle \frac{\partial \mathrm{WUE}}{\partial P}}\right)}_{\mathrm{2}}\mathrm{\Delta }P+{\left({\displaystyle \frac{\partial \mathrm{WUE}}{\partial \mathrm{PET}}}\right)}_{\mathrm{2}}\mathrm{\Delta }\mathrm{PET}\right.\\ & \left.+P_{\mathrm{1}}\mathrm{\Delta }\left({\displaystyle \frac{\partial \mathrm{WUE}}{\partial P}}\right)+{\mathrm{PET}}_{\mathrm{1}}\mathrm{\Delta }\left({\displaystyle \frac{\partial \mathrm{WUE}}{\partial \mathrm{PET}}}\right)\right],\end{array}\end{array}$$

where subscripts 1 and 2, respectively, represent the initial- and final-state quantities; Δ denotes the difference between the two states; and α is the weighting factor, with a value range of 0 to 1.

Considering Eqs. (4) and (8), ΔWUE can be further written as changes due to the climate (ΔWUE_c) and due to the underlying surface (ΔWUE_u):

$$\begin{array}{}\text{(14)}& {\displaystyle }\mathrm{\Delta }\mathrm{WUE}=\mathrm{\Delta }{\mathrm{WUE}}_{\mathrm{c}}+\mathrm{\Delta }{\mathrm{WUE}}_{\mathrm{u}},\text{(15)}& {\displaystyle }\begin{array}{rl}\mathrm{\Delta }{\mathrm{WUE}}_{\mathrm{c}}& =\mathit{\alpha }\left[{\left({\displaystyle \frac{\partial \mathrm{WUE}}{\partial P}}\right)}_{\mathrm{1}}\mathrm{\Delta }P+{\left({\displaystyle \frac{\partial \mathrm{WUE}}{\partial \mathrm{PET}}}\right)}_{\mathrm{1}}\mathrm{\Delta }\mathrm{PET}\right]\\ & +\left(\mathrm{1}-\mathit{\alpha }\right)\left[{\left({\displaystyle \frac{\partial \mathrm{WUE}}{\partial P}}\right)}_{\mathrm{2}}\mathrm{\Delta }P+{\left({\displaystyle \frac{\partial \mathrm{WUE}}{\partial \mathrm{PET}}}\right)}_{\mathrm{2}}\mathrm{\Delta }\mathrm{PET}\right],\end{array}\text{(16)}& {\displaystyle }\begin{array}{rl}\mathrm{\Delta }{\mathrm{WUE}}_{\mathrm{u}}& =\mathit{\alpha }\left[P_{\mathrm{2}}\mathrm{\Delta }\left({\displaystyle \frac{\partial \mathrm{WUE}}{\partial P}}\right)+{\mathrm{PET}}_{\mathrm{2}}\mathrm{\Delta }\left({\displaystyle \frac{\partial \mathrm{WUE}}{\partial \mathrm{PET}}}\right)\right]\\ & +\left(\mathrm{1}-\mathit{\alpha }\right)\left[P_{\mathrm{1}}\mathrm{\Delta }\left({\displaystyle \frac{\partial \mathrm{WUE}}{\partial P}}\right)+{\mathrm{PET}}_{\mathrm{1}}\mathrm{\Delta }\left({\displaystyle \frac{\partial \mathrm{WUE}}{\partial \mathrm{PET}}}\right)\right].\end{array}\end{array}$$

From Eqs. (14)–(16), the impacts of climate and underlying surface changes on WUE can be precisely decomposed within the Budyko model. However, attribution solely through the Budyko model is insufficient to decouple the effects of CO₂ on land and atmosphere. The CO₂ sensitivity experiments of GCMs mentioned in Sect. 2.2.4 had been applied to explore the atmosphere–surface interactions caused by elevated CO₂ (Piao et al., 2007; Fowler et al., 2019). In the 1ptCO2-rad experiment, the increase in CO₂ concentration has a direct radiative effect on the climate and an indirect effect on vegetation physiology; conversely, in the 1ptCO2-bgc experiment, the increase in CO₂ concentration has a direct biogeochemical impact on vegetation and then indirectly affects the climate. The difference between the combined results of the -rad and -bgc experiments and the results of the 1ptCO2 experiment are attributed to land–atmosphere interactions. By applying the Budyko attribution method in the three CO₂ experiments (Zhou et al., 2023), the impact of land–atmosphere coupling on WUE can be quantitatively interpreted.

3.3 TWSA-DSI

The method proposed by Zhao et al. (2017) was employed to identify terrestrial drought conditions using the TWSA-DSI index. TWSA-DSI is capable of comparing drought characteristics across regions and periods, and it can also sensitively capture the available water amount affecting vegetation growth. The calculation is as follows:

$$\begin{array}{}\text{(17)}& \mathrm{TWSA}-{\mathrm{DSI}}_{i,j}=({\mathrm{TWSA}}_{i,j}-{\stackrel{\mathrm{‾}}{\mathrm{TWSA}}}_j)/{\mathit{\sigma }}_j,\end{array}$$

where $\mathrm{TWSA}-{\mathrm{DSI}}_{i,j}$ and TWSA_i,j represent the TWSA-DSI and TWSA data for the ith year and jth month and ${\stackrel{\mathrm{‾}}{\mathrm{TWSA}}}_j$ and σ_j, respectively, denote the average and standard deviation of TWSA for the jth month within the study period.

In this study, the TWSA-DSI index was calculated for each grid point within the basin for each study period (experiment). Based on these results, different levels of terrestrial wetness/dryness were classified as shown in Table 2 (Yin et al., 2022). Three main categories of terrestrial wet and dry conditions were investigated in this study: non-drought (TWSA-DSI > −0.8), moderate–severe drought (−1.6 < TWSA-DSI ≤ −0.8), and extreme–exceptional drought (TWSA-DSI ≤ −1.6).

Table 2 Classification standards for wetness (W) or dryness (D) levels based on TWSA-DSI.

3.4 Conditional probability

The impact of different drought levels on WUE can be assessed by considering the conditional probability distribution of the drought index. Referring to previous studies (Feng et al., 2019; Wu and Jiang, 2022), WUE is transformed into a standardized form using a meta-Gaussian model specifically as the standardized water use efficiency index (SWI). The differences in the distribution of SWI under different drought conditions (with TWSA-DSI as the drought index) were used to evaluate the response of WUE to different drought severity levels. As TWSA-DSI is also a standardized normal random variable (X), the bivariate conditional distribution of SWI (Y) given X can be straightforwardly expressed as follows:

$$\begin{array}{}\text{(18)}& Y\mathrm{|}X\sim N\left({\mathit{\mu }}_{Y\mathrm{|}X},{\mathrm{\Sigma }}_{Y\mathrm{|}X}\right),\end{array}$$

where μ_Y|X is the conditional mean and Σ_Y|X is the conditional variance.

4 Results

4.1 Assessment of bias correction performance

We applied bias correction to the monthly GCM data using the method described in Sect. 3.1 and compared the correction results for the overlapping period (1997–2014) with observation data. It is noteworthy that different GCMs include a varying number of variables. As a result, even though we have endeavored to select the same models, the ensemble GCM members for different variables are not exactly the same. The performance of correction was presented in the form of Taylor diagrams (Fig. 2). The correction results for GPP and ET were substantially better than those for P, PET, and TWSA. The correlation coefficients (CCs) for GPP and ET after correction were all above 0.97, with multi-model average root mean square errors (RMSEs) of 9.90 gC m⁻² and 7.35 mm, respectively. The inter-model differences were also small. Thus, it can be assumed that the WUE calculated from the corrected output data was credible. For P, PET, and TWSA, the models with the best correction effects were IPSL-CM6A-LR, CMCC-ESM2, and CESM2, respectively. Their CC reaches 0.89, 0.99, and 0.82, with RMSE values of 18.13, 8.35, and 10.16 mm. In comparison, CMCC-ESM2 performed relatively poor for P with a CC of 0.74. CNRM-ESM2-1 exhibited suboptimal correction performance for PET and TWSA, with CC values of 0.78 and 0.70, respectively. Overall, the outputs of all GCMs showed good performance after bias correction (with CCs all greater than or equal to 0.70). The trend-preserving bias correction method effectively eliminated systematic errors between GCMs and observation data. The corrected data can be used for the assessment and attribution of future climate scenarios.

Figure 2 Taylor diagrams for corrected GCM data of (a) GPP, (b) ET, (c) P, (d) PET, and (e) TWSA.

4.2 Changes in WUE over YRB

Figure 3 illustrated the time series of multi-model annual averages for WUE, GPP, and ET under historical and different future scenarios. During the historical period, WUE and GPP in the YRB showed an increasing trend, while ET exhibited a slight rise with fluctuations. From 2015 to 2100, the changes in each indicator varied in different scenarios, but the trend patterns of WUE and GPP are relatively similar (Fig. 3a and b). Under the low-emission pathway scenario (SSP126), WUE and GPP first experienced a slight increase in the future (peaking around 1.68 gC kg⁻¹H₂O and 675 gC m⁻², respectively), which was then followed by a decline. In the moderate-emission scenario (SSP245), WUE and GPP increased before stabilizing at the end of the century (2070–2099), with values of 1.55 gC kg⁻¹H₂O and 771 gC m⁻², respectively. Both WUE and GPP exhibited a continuous upward trend in the high-emission scenarios (SSP370 and SSP585), and SSP585 showed a higher growth rate compared to SSP370. For ET (Fig. 3c), there was a subtle increase under the low-emission scenario (SSP126) and a slight decrease towards the end of the century. Guo et al. (2022) observed that temperature and precipitation demonstrate similar developmental trends in the YRB under SSP126 as the low-emission pathway aims to limit global warming. Given that temperature and precipitation are crucial drivers of evapotranspiration, it is logical to see ET following this pattern of development. In other scenarios, ET showed an overall upward trend, and with increasing carbon emissions, ET in the YRB increased more rapidly. Additionally, it can be observed from the box plots that the outputs for WUE, GPP, and ET from the multi-model ensemble were more concentrated during the historical period, whereas they were more dispersed in the future periods. Furthermore, the degree of dispersion increased with the rise in carbon emission concentrations.

Figure 3 Annual averaged time series and box plots during historical-period (1985–2014) and future-period (2070–2099) (a) WUE, (b) GPP, and (c) ET from the CMIP6 multi-model ensemble in the YRB.

The first column in Fig. 4 displays the spatial distribution of annual WUE in the YRB during the historical period (1985–2014) and the future period (2070–2099) under different SSPs. The multi-model average results indicated that during the historical period, the multi-year average WUE in the YRB predominantly exhibited a pattern of higher values in the east and lower values in the west, with higher values in the south and lower values in the north. High WUE values were observed in the upstream areas and downstream coastal areas, while low values appeared in the northwest part of the basin, specifically the Inner Mongolia Plateau. Under different scenarios, the future distribution characteristics of WUE remained similar. As shown in the first column of Fig. 5, there was a predominant upward trend in WUE in the middle-downstream region in the future, and the relative growth increased with increasing carbon emissions. In some areas, the growth exceeded 200 %. However, the Hetao region in the northwest part of the Yellow River consistently experienced a downward trend in WUE under different scenarios, with the greater decline as carbon emissions increased. In the SSP585 scenario, this region experienced a negative growth of close to 20 %. The upstream areas of the YRB showed relatively consistent changes under different scenarios except for a noticeable decrease in the southwestern of the basin, i.e., the Qinghai–Tibet Plateau region, while changes in other regions were less evident.

To further investigate the spatial distribution characteristics of WUE, we plotted the spatial distribution and relative changes in GPP and ET (different columns in Figs. 4 and 5). The spatial distribution characteristics of GPP closely resembled those of WUE (the first and second column in Fig. 4), with the locations of high and low values aligning well. Under different scenarios, the entire YRB experienced an increase in GPP, though the degree varies (the second column in Fig. 5). Overall, there was relatively little change in GPP in the source area of the YRB, while significant growth was observed in the Loess Plateau and the intra-basin areas. The growth in GPP increased sequentially from SSP126 to SSP585. In SSP585, most regions in the basin, including the Loess Plateau and Ningxia Plain, exhibited a 260 % growth in GPP. Extremely low values ET were observed in the Hetao area, with ET increasing in a radial pattern, indicating higher ET values at locations increasingly distant from this region (the third column in Fig. 4). Throughout different periods and scenarios, the minimum values consistently appeared in the northwest part of the basin, specifically the Inner Mongolia Plateau, all below 60 mm. The maximum values occurred along the Qinling (which also serves as the boundary between northern and southern China), where ET exceeds 780 mm. The relative changes in ET in the YRB exhibited the opposite trend, increasing from southeast to northwest (the third column in Fig. 5). In the low-emission scenario (SSP126), ET showed minimal variation within the basin. However, with increasing carbon emissions, the growth in ET became more drastic, especially in the Inner Mongolia Plateau. In SSP585, the growth in this area and its vicinity reached 52 %. In comparison, ET changes in other regions are not significant except for in the source area, where ET also increased.

Figure 4 Spatial distribution of multi-model average results in the YRB under different scenarios. The three columns represent the multi-year average of WUE (the first column), GPP (the second column), and ET (the third column) during the historical (the first row) period (1985–2014) and future periods (2070–2099) under the SSP126 (the second row), SSP245 (the third row), SSP370 (the fourth row), and SSP585 (the fifth row) scenarios.

Figure 5 Spatial distribution of the relative changes in WUE (the first column), GPP (the second column), and ET (the third column) compared to the historical period under different scenarios.

4.3 Attribution of WUE changes over the YRB

We constructed a Budyko-type model of water–energy exchange for WUE based on the method outlined in Sect. 3.2. Figure 6a illustrates the Spearman correlation coefficients between the underlying surface parameters of the Budyko model, WUE, and vegetation indices based on observation data. WUE showed a strong positive correlation with LAI and NDVI, with correlation coefficients reaching 0.97 and 0.95, respectively. The correlation with the underlying land surface parameter, m, was also significant (p<0.01), with a coefficient of 0.56. Although the negative correlation between WUE and n was relatively weaker (−0.26), it was still statistically significant (p<0.01). The basin underlying surface parameters (m and n) exhibited good correlations with vegetation indices (LAI and NDVI) at 0.53, 0.65, −0.30, and −0.33, respectively. In general, an increase in LAI and NDVI will cause an increase in m and a decrease in n, consequently resulting in an increase in WUE.

Figure 6 Performance of the Budyko model. (a) Cross-correlation coefficients among WUE, Budyko model parameters, and vegetation indices based on observation data (^* indicates a significance of p<0.01). (b) Budyko model performance for WUE in the historical period (1985–2014) and under different SSPs (2070–2099). The shading represents the 95 % confidence interval.

Furthermore, we established Budyko models for WUE using the multi-model average data under different periods and scenarios to assess the combined impacts of climate and vegetation patterns on WUE. As depicted in Fig. 6b, the selected linear Budyko models demonstrated good performance across various periods and scenarios. It exhibited the best fits in SSP370 and historical scenarios, with R² values of 0.73 and 0.71, respectively, while performing relatively less well in SSP245, with an R² of 0.60. Overall, the models' performance had been satisfactorily validated in the YRB (R² consistently exceeding 0.60), providing confidence for subsequent attribution analyses.

Based on the established Budyko models, we calculated the changes in WUE and the attribution results for the YRB (Fig. 7a). The results indicated an increase trend in the multi-year average WUE for the basin in the future period (2070–2099) compared to the historical period (1985–2014). Moreover, with intensifying carbon emissions, the increase in WUE became more substantial, with average increases in WUE sequentially from SSP126 to SSP585 by 0.36, 0.54, 0.75, and 0.84 gC kg⁻¹ H₂O. SSP585 exhibited the highest growth in WUE, although we observed a slowdown in the increase from SSP370 to SSP585. According to the attribution results, climate change consistently remained the predominant factor influencing WUE changes, especially in SSP126, where almost the entire variation in WUE was attributable to climate change (101.4 %). In SSP126, SSP245, SSP370, and SSP585, the proportions dominated by underlying change were −1.4 %, 14.9 %, 16.9 %, and 22.1 %, respectively. However, we also observed that with increasing carbon emissions, the impacts of underlying surface changes on WUE became more significant. In the low-emission scenario (SSP126), the influence of land surface changes in the basin was minimal, only playing a slightly negative role (−1.4 %). But in SSP245, SSP370, and SSP585, the promoting effects of land surface factors on WUE became more pronounced. In SSP585, nearly a quarter of the WUE growth was attributed to changes in the underlying surface. These phenomena could be explained by the fertilization effect of carbon dioxide. In the SSP126 scenario, atmospheric CO₂ concentrations for the future period (2070–2099) were projected to be lower than in the historical period (1985–2014), which could be detrimental to plant growth and especially GPP. Therefore, under the low-emission scenario, changes in the land surface might have a slight impact on WUE. With the increase in carbon emissions, the growth of vegetation was further promoted, and the changes in vegetation productivity due to different carbon emission scenarios might be more pronounced compared to the impacts of climate change. Thus, the role of underlying land surface changes became increasingly important in the variation in WUE.

Figure 7 Attribution results for WUE changes. (a) Changes in WUE in the YRB under different SSPs along with the relative contributions of climate factors and underlying surface factors. The stacked chart represents the relative contribution rates of climate change and underlying changes to WUE variations under different scenarios. The line graph depicts the average WUE changes (during 1985–2014 and 2070–2099) across the basin under various scenarios. (b) Land–atmosphere decoupling attribution results based on SSP585 and CO₂ experiments. The percentages represents the relative contribution rates.

Due to the consistent increase in CO₂ concentration in the three CO₂ sensitivity experiments at the same rate as in SSP585, we combined the CO₂ experiments with the attribution in SSP585 to decouple the direct and secondary effects of CO₂-induced radiative and biogeochemical impacts. To ensure the validity of the results, we selected the first 30 years and the subsequent 30 years of the CO₂ experiment, with a time span consistent with SSP585. As shown in Fig. 7b, compared to the radiative effect (−1.3 %), the biogeochemical effect of CO₂ overwhelmingly dominated the change in WUE (100.6 %). This was understandable as the majority of the effects from the biogeochemical experiment directly acted on land surface (94.5 %) and specifically on the vegetation structure. Climate-induced feedback on changes in vegetation structure also played a positive role in WUE growth (6.1 %). Analyzing the attribution results from the radiative experiment, we found that although the radiative effect had a minimal negative impact on WUE (−1.3 %), this outcome was due to the significant weakening or even surpassing of the indirect response of the land surface vegetation structure to the climate (−56.3 %) compared to the direct radiative forcing effect of CO₂ on the climate (55.0 %). The mutual offsetting of these two results led to the radiative effect having a very small negative impact on WUE changes. However, the direct forcing of the radiative effect on climate (55.0 %), the indirect feedback in the biogeochemical experiment (6.1 %), and the interaction between the two experiments (16.8 %) ultimately made climate change the dominant factor (77.9 %) in the growth of the YRB.

4.4 WUE response to drought

In order to investigate how WUE responds to different levels of drought under various scenarios, we calculated the anomaly of WUE at each pixel at different drought levels and averaged the anomaly results across different GCMs. To obtain a sufficient number of event samples at the grid scale and derive more universally applicable conclusions, we focused only on the moderate–severe drought (−1.6 < TWSA-DSI ≤ −0.8) and extreme–exceptional drought (TWSA-DSI ≤ −1.6) categories by treating all other conditions as non-drought situations (TWSA-DSI > −0.8) (Yin et al., 2022). In the analysis of future drought responses, particular attention was given to the potential influence of long-term trends on the calculated values of WUE and TWSA. To mitigate this effect, we calculated the mean and standard deviation of WUE and TWSA for the same time spans separately to get the WUE anomaly and TWSA-DSI, which ensures that our analysis focuses on the specific variations and conditions within each period. As shown in Fig. 8 (the first and second columns), in most regions of YRB, the response of WUE to drought can be divided into two stages: an increase in WUE during moderate–severe drought and a decrease in WUE during extreme–exceptional drought. The area percentages of the two-stage model in different periods and scenarios are 98.0 % in historical, 96.8 % in SSP126, 97.6 % in SSP245, 97.0 % in SSP370, and 99.2 % in SSP585. The spatial characteristics of WUE response patterns were also similar across different periods and scenarios.

Figure 8 Responses of WUE to different levels of drought. The first column shows the anomalies of WUE during moderate–severe drought. The second column shows the anomalies of WUE during extreme–exceptional drought. The third column illustrates the conditional probability function of SWI under different drought levels.

In the historical period (1985–2014), WUE generally increased during moderate–severe drought, with a widespread increase, and only a few areas (specifically, two pixels in the Ningxia Plain) showed a decrease. The middle reaches of the YRB showed a relatively small increase, while more significant increases were observed in the upstream source area and downstream estuary. As drought intensified to extreme–exceptional levels, there was a large-scale decrease in WUE across the basin. The areas with increased WUE compared to the previous stage also experienced a decline in this stage. Similar to moderate–severe drought, the middle reaches of the YRB saw a significant decrease in WUE, while the source and downstream areas showed a slight increase.

Looking into the future (2070–2099), WUE's positive response to drought became more prominent under all SSPs. This was evident in the significantly higher increase in WUE during moderate–severe drought compared to the historical period. Moreover, the areas with a decrease in WUE during extreme–exceptional drought became smaller. With increasing carbon emissions, the positive response of WUE to drought became more pronounced. This was due to the fact that CO₂ not only promoted the growth of the plant, but also reduced stomatal conductance, allowing water to be better retained in the body during dry periods. In SSP126, the WUE response pattern to drought was similar to that of the historical period. However, in SSP585, almost all regions within the basin showed a positive response of WUE to moderate–severe drought, with WUE anomalies reaching around 1.0 gC kg⁻¹H₂O. Even during extreme–exceptional drought, there were very few areas where WUE decreased across the basin.

Figure 8 (the third column) quantitatively described the conditional probability distribution of SWI under different dry and wet conditions. Across different periods and SSPs, the mean (μ) of SWI under extreme–exceptional drought conditions was lower than that under moderate–severe drought conditions. Combining the previously mentioned two-stage response characteristics, this is understandable because as drought intensified, the WUE response in the basin tended to shift towards negative values. However, we also observe that μ under non-drought conditions was generally higher than that under drought conditions. A reasonable explanation was that, compared to the wet season, the basin tended to have lower WUE during seasons or months prone to drought. Therefore, even with a positive response during drought, the value of μ remained lower than in non-drought periods. Additionally, compared to the historical period, the future μ under different scenarios showed varying degrees of growth. This growth was less pronounced in the low-emission scenario (SSP126) but became more significant with increasing carbon emissions, reaching its maximum in SSP585. This was consistent with the results in Sect. 4.2.

Furthermore, we observed that with the increase in carbon emissions, the gap between μ under non-drought conditions and that under drought conditions became more apparent. The difference in μ (Δμ) represents the difference in the mean values of WUE across varying intensities of drought within the same period. This measure, utilized to analyze the differential impacts of drought levels, also became larger (from SSP126 to SSP585, Δμ is 0.01, 0.01, 0.02, and 0.12, respectively). This implied that a high-carbon-emission scenario would exacerbate the numerical differences in WUE between dry and wet seasons in the YRB while also amplifying the response differences in WUE to different drought levels. Moreover, in SSP585, these differences were significantly greater than in other scenarios.

5 Discussion

5.1 Spatiotemporal variation characteristics of WUE in the YRB

Unlike many studies that focus solely on the historical trends of WUE (Kim et al., 2021; Huang et al., 2017; Lu and Zhuang, 2010), this research systematically projected the development trends of WUE, GPP, and ET in the YRB for both historical and future periods. GPP, as a key variable in calculating WUE, exhibits significant differences among various GPP datasets, such as NIRv, GIMMS, FLUXCOM (based on upscaled eddy covariance flux tower measurements), LUE (based on the light-use efficiency model), and TRENDY (a recent model intercomparison project) (S. Wang et al., 2021). Compared to NIRv GPP data, GIMMS ignores the impact of CO₂ on LUE (De Kauwe et al., 2016), and FLUXCOM GPP lacks consideration of the CO₂ fertilization effect (Anav et al., 2015), thus underestimating the interannual sensitivity of GPP to climate change. A recent study has improved the traditional LUE method using eddy covariance, which produced a revised global GPP product (EC-LUE GPP) (Zheng et al., 2020). However, satellite-based APAR (absorbed photosynthetic radiation) data used may lead to inaccuracies in some regions due to saturation (Zhang et al., 2020). Furthermore, the GPP trends in TRENDY vary significantly between different models and offer lower spatial resolution (Zheng et al., 2020). Considering all these factors, we opted for the advanced GPP dataset based on NIRv with high temporal and spatial resolutions. Our results also indicate that these data are reasonable and accurate. During the historical period, the average results from multiple GCMs closely align with previous findings (H. Li et al., 2023; Sun et al., 2022; Zhao et al., 2022b), which indirectly demonstrated the excellent performance of bias correction done in this study. According to Fig. 3, WUE, GPP, and ET were projected to increase to varying degrees in the future. The growth rate was higher under high-emission scenarios, while under the SSP126 scenario, all three showed an initial increase followed by a decline. F. Li et al. (2023) suggested that global WUE had been approaching saturation in recent years and may maintain the saturated state in the future. According to our results, under SSP245 (which is also most consistent with the Chinese context), WUE in the YRB is expected to reach a relatively stable “saturated” state by the end of this century, with a value of around 1.55 gC kg⁻¹H₂O. The spatial distribution of WUE in the YRB can be summarized as “higher in the south and lower in the north, higher in the east and lower in the west” (Fig. 4). This spatial pattern and the regions where extreme values occur align well with the study of Sun et al. (2022). We also observed that the future spatial distribution of WUE was very similar to GPP (Fig. 4) (Liu et al., 2020). Combining these results with those in Fig. 2, despite the general increase in ET, WUE in most regions of the YRB is still on the rise. Therefore, we conclude that in most areas of the YRB, GPP's relative contribution to WUE change is higher than that of ET. In other words, GPP dominates the future changes in WUE in the basin, which is in good agreement with other scholars' findings (Naeem et al., 2023; Tan et al., 2023). However, ET remains a primary driving factor for WUE in arid regions (Yang et al., 2016), a point especially evident in the Yellow River source region. In the upstream source region of the YRB, ET shows an increasing trend under different SSPs, while GPP's growth is relatively slow. This asymmetrical growth between GPP and ET has led to a decline in WUE in the source region.

5.2 Interpretation of Budyko model parameters and attribution results

In this study, we employed the linear Budyko model to describe the relationship between WUE and the water–energy nexus in the YRB. In Eq. (3), m and n are the underlying surface parameters of the basin. To validate the rationality of these parameters, we conducted a correlation analysis between m, n, and vegetation indices as shown in Fig. 6a. The results indicate good correlations between underlying surface parameters and LAI, NDVI. This suggests that the chosen parameters can effectively reflect the underlying surface vegetation conditions closely related to WUE, further validating the physical significance of the established Budyko model. Additionally, we demonstrated a strong positive correlation between WUE and both LAI and NDVI due to the significant enhancement of the region's GPP by increased vegetation factors. In the mostly arid areas of the YRB, the impacts of these vegetation changes on ET is minimal (Liu et al., 2020; Yang et al., 2022). Moreover, our results also support the viewpoint that there is a positive correlation between WUE and m (Fang et al., 2020). Previous studies have shown that WUE is closely related to precipitation, vapor pressure deficit, temperature, wind speed, humidity, et cetera (F. Li et al., 2023; Lin et al., 2020; Liu et al., 2020; Yang et al., 2022). This study employed the Penman–Monteith equation to calculate PET, incorporating the aforementioned factors; hence, P and the calculated PET can be selected as the climatic attribution factors for WUE.

Figure 9 Spatial distribution of the relative changes in WUE (the first column), GPP (the second column), and ET (the third column) across different CO₂ sensitivity experiments.

To further elucidate the attribution results shown in Fig. 7b, we plotted the relative changes in WUE, GPP, and ET across different CO₂ sensitivity experiments (Fig. 9). Over the same time span, the spatial distribution pattern of relative changes in WUE, GPP, and ET within the 1ptCO2 scenario aligns substantially with that of the SSP585 scenario. In the 1ptCO2-rad experiment, WUE showed a declining trend in the peripheral areas of the YRB but an increase in the central region, thus presenting a slight negative contribution overall (−1.3 %) as depicted in Fig. 7b. Conversely, in the 1ptCO2-bgc experiment, WUE exhibited a comprehensive upward trend across the basin, serving as the primary driving factor for the increase in WUE observed in the comprehensive experiment (1ptCO2) (100.6 %). To explain the difference in WUE performance, we further analyzed the relative changes in GPP and ET across different CO₂ sensitivity experiments. In the 1ptCO2-rad experiment, due to the absence of CO₂'s direct biogeochemical effect on plants, the change in GPP was not as sharp as that in the 1ptCO2-bgc. Moreover, the radiative effects of CO₂ increase atmospheric evaporative demand, i.e., PET (Milly and Dunne, 2016, 2017), which is the most likely reason for the widespread increase in ET and the broad decline in WUE within the YRB. However, combining the analysis with Fig. 7b, such changes encompass CO₂'s direct impact on climate and the indirect feedback from the underlying surface, with the latter possibly exerting a greater influence on ET and GPP variations. However, due to the minimal numerical changes in variables in 1ptCO2-rad, the uncertainty is also higher. Under the 1ptCO2-bgc experiment conditions, owing to the fertilization effect of CO₂, plant growth was stimulated (He et al., 2023; Zhao et al., 2022a), resulting in a significant widespread increase in GPP across the basin. The increased plant biomass also led to a rise in ET (Mankin et al., 2019; Piao et al., 2007; Zhan et al., 2022), while the elevated CO₂ concentrations further reduced plant stomatal conductance, potentially diminishing transpiration (Guerrieri et al., 2019; F. Li et al., 2023; Mathias and Thomas, 2021; Zhang et al., 2022). We found that with the combined effects of these factors, changes in ET within the YRB are minimal, generally ranging between −5 % and 5 % in 1ptCO2-bgc experiment. Consequently, in 1ptCO2-bgc, WUE across various regions of the basin exhibits varying degrees of increase.

Table 3 Model parameter settings for simulating stomatal conductance in the GCMs.

CLASS-CTEM: coupled Canadian Land Surface Scheme and the Canadian Terrestrial Ecosystem Model. CLM5: Community Land Model version 5. CLM4.5: Community Land Model version 4.5. SURFEXv8.0: Surface Externalisée version 8.0. ORCHIDEE: Organizing Carbon and Hydrology in Dynamic Ecosystems Land Surface Model. g_s is the stomatal conductance; g₀ is the minimum stomatal conductance as assimilation rate reaches zero; g₁, D₀, b, m, n, and β are the empirical coefficients; A is the net photosynthesis rate; D is the humidity deficit; c_s is the CO₂ concentration at the leaf surface; τ is the CO₂ compensation point; c_a is the atmospheric CO₂ concentration; h_r is the relative air humidity; P_atm is the atmospheric pressure; and c_i is the CO₂ concentration in the intercellular spaces.

5.3 Two-stage response pattern of WUE to drought

WUE can reflect the impacts of local water availability on the ecosystem's carbon sequestration capacity to a certain extent and can serve as one of the reference indicators for measuring ecosystem resilience (Ponce-Campos et al., 2013). Therefore, we considered the response of WUE in the YRB to drought stress. Many scholars have suggested that WUE will continue to increase or decrease under drought conditions in most regions of the world (Liu et al., 2023; Ponce-Campos et al., 2013; Xie et al., 2016). However, these studies have not considered the differential responses of WUE to various drought intensities. Our results (Fig. 8) illustrate that WUE exhibits a two-stage response pattern during drought; i.e., it increases in moderate–severe drought but begins to decrease as the drought intensifies. This finding is similar to the phenomenon observed by Lu and Zhuang (2010) over most of the US. This means that under the moderate–severe drought condition, plants physiologically adapt to water stress, thereby maintaining a certain level of GPP under limited water conditions (Vicente-Serrano et al., 2010). However, as drought intensifies, plant physiology is damaged, thereby reducing WUE, which is consistent with findings by Yang et al. (2021) in semi-arid/semi-humid regions of the world. As mentioned before, the increase in atmospheric CO₂ concentration directly leads to a reduction in plant leaf stomatal conductance. In Table 3, it is evident that although different GCMs employ various models for stomatal conductance, the physical implications across these models consistently align with this principle. It is the change in stomatal conductance that is suggested to help alleviate the stress of drought on plants (De Kauwe et al., 2021; Leakey et al., 2006; Swann, 2018). Additionally, we note that with the intensification of carbon emission scenarios (from SSP126 to SSP585), the future response of WUE in the basin to drought becomes more positive, indicating an enhanced adaptability of plants to water stress. The results also indirectly prove that plants are more likely to benefit from drought stress in the future (De Kauwe et al., 2021).

6 Conclusion

WUE, due to its unique connotation, has become an important indicator for ecosystem health, reflecting the tradeoff between regional carbon assimilation and water loss. In this study, after applying trend-preserving bias correction to CMIP6 data informed by satellite-monitored reanalysis datasets, we investigated spatiotemporal variations in WUE, GPP, and ET spanning from 1985 to 2100. Our construction of a Budyko model, underpinned by GCMs, offers a robust framework for dissecting the intricacies of WUE evolution, revealing a predominantly climate-driven narrative, albeit with an escalating influence from underlying surface modifications under severe carbon emission trajectories.

The spatial dynamics of WUE across the YRB have historically been characterized by higher values in the east and the south contrasted with lower values in the west and the north. Looking forward, our analysis projects a universal uptick in WUE under varying future climate scenarios, with increases ranging from 0.36 to 0.84 gC kg⁻¹ H₂O across SSP126 to SSP585. The constructed Budyko model displays strong performance in different periods and under different scenarios (R²= 0.60–0.73). The land–atmosphere attribution framework underscores the pivotal role of climate change as the primary driver of WUE alterations, overshadowing underlying surface changes. Nevertheless, the impact of biogeochemical CO₂ effects emerges as a dominant force, illustrating the complex interplay of carbon dioxide in shaping WUE trajectories. Besides, our findings delineate a distinctive two-stage response pattern of WUE to drought conditions in the YRB and forecast enhanced ecosystem resilience to drought stress in future scenarios, particularly under heightened carbon emissions.

While this study advances our understanding of ecohydrological processes in the YRB, it is not without its constraints. The reliance on reanalysis datasets, despite their robust performance, introduces a layer of uncertainty when compared against actual measurements. Moreover, the complexity of vegetation responses to elevated CO₂ levels coupled with the variability in dynamic global vegetation models (DGVMs) across different GCMs calls for a deeper exploration into the quantification and assessment of uncertainties surrounding CO₂'s impact on WUE.

By charting the interdependencies of CO₂ dynamics, climate change, and land surface alterations, this study concerns the future of ecohydrological sustainability in the YRB, inviting further inquiry into the resilience of terrestrial ecosystems under anthropogenic stress.