Enhancing evapotranspiration estimates under climate change: the role of CO₂ physiological feedback and CMIP6 scenarios

The future state of global evapotranspiration (ET) estimation under climate change remains uncertain. Current formulations primarily developed based on the high emission CMIP5 scenario, have been widely used to represent conditions under elevated greenhouse gas pathways. However, these formulations may not adequately capture the enhanced vegetation–climate interactions projected under the lower-emission scenarios of CMIP6. Without updates to account for evolving plant physiological responses to rising CO₂, projections may overlook critical feedbacks between atmospheric CO₂ concentrations, vegetation behavior, and hydrological processes.

To address this, developing CMIP6-specific formulations is essential to leverage its improved datasets and reduce uncertainties in future ET simulations. In this study, we update the Penman-Monteith evapotranspiration (PM-ET) model by incorporating the CO₂-vegetation coupling effect. This is achieved using outputs from four Coupled Model Intercomparison Project Phase 6 (CMIP6) global climate models (GCMs) under four Shared Socioeconomic Pathways (SSP1-2.6, SSP2-4.5, SSP3-7.0, SSP5-8.5).

Results indicate a sustained historical increase in potential evapotranspiration (Ep). The inclusion of CO₂ physiological effects reduces the deviation in projected ET trends by approximately 15 %–20 % compared to CMIP5-based frameworks, accounting for the increase in stomatal resistance driven by CO₂ concentrations rising from ∼284 to ∼ 935 ppm. Furthermore, our model predicts an increasing dependence of ET projections on emission scenario, highlighting the growing influence of pathway-specific feedbacks.

Overall, our approach demonstrates greater compatibility with CMIP6 simulations, allowing for more accurate representation of ET responses to future CO₂ increases. These findings provide valuable insights for advancing the analysis of nonlinear vegetation-atmosphere interactions and hydrological uncertainty under climate and physiological forcings.

1 Introduction

Vegetation-climate interactions under rising atmospheric CO₂ concentrations drive complex biogeochemical feedbacks that shape the global carbon-water cycle across diverse biomes. As a result, the relationship between vegetation dynamics and climate change has gained increasing prominence in Earth system science. The bidirectional feedbacks are fundamental to understand changes in global hydrological and carbon cycles (Yang et al., 2019; Xu et al., 2024).

On one hand, elevated CO₂ levels trigger complex plant physiological responses, which influence the ecosystem hydrology (Sun et al., 2023; Li, 2024). For example, studies have shown that photosynthetic rates can increase by 12 %–25 % for every 100 ppm rise in CO₂, while stomatal conductance decreases by 20 %–40 %. This reduction limits evaporative cooling, exacerbating extreme temperatures by 0.8–1.3 °C during drought periods (Zarakas et al., 2020; Li et al., 2024; Wu et al., 2024). Moreover, biome-specific responses create climatic asymmetries (Yang et al., 2023). Temperate forests exhibit adaptive resilience through phenological shifts, extending the growing season by 15.6±3.2 d to mitigate respiratory carbon losses during heat extremes (Ainsworth and Rogers, 2007; Keenan and Williams, 2018). In contrast, tropical ecosystems experience accelerated biomass loss due to compound warming and precipitation anomalies (Betts et al., 2007; Gimeno-Sotelo et al., 2024a, b). Soil carbon-climate interactions add further complexity, with elevated CO₂ increased mineralization rates ...(Kong et al., 2023). Critical thresholds in plant hydraulic strategies also contribute to these differences (Medlyn et al., 2011; Wu et al., 2019).

On the other hand, the strength of CO₂-driven vegetation feedbacks is constrained by water availability. Since 2000, compound drought-heatwave events have reduced global terrestrial water storage by 12 %–18 % and vegetation productivity by 9 %–15 % (Piao et al., 2007; Jones et al., 2016; Zhou et al., 2016; Lu et al., 2025). These phenomena are attributed to sustained evapotranspiration pressure (Xu et al., 2024). Although CO₂ fertilization initially enhances vegetation growth in water-limited regions, subsequent hydrological constraints often offset these benefits (Jasechko et al., 2013; Young et al., 2022). Therefore, accurately capturing the interplay between biome-specific responses and hydrological limitations is essential for understanding the cascading effects of climate and physiological drivers on the carbon–water cycle.

Progress in evapotranspiration (ET) modelling is challenged by the need to account for CO₂-climate coupling and regional hydrological variability. Widely used ET estimation methods, such as the Penman-Monteith Reference Crop (PM-RC) model, lack explicit consideration of CO₂ effects (Stocker et al., 2018; Wu et al., 2021). As a result, traditional PM-RC models tend to overestimate future ET trends by 68 %–100 % due to inadequate representation of CO₂-induced stomatal closure (Yang et al., 2016). To address this, Yang et al. (2019) developed a modified framework incorporating CO₂-dependent stomatal resistance, with improved ET projections under high CO₂ scenarios (Luo et al., 2018a, b). While high-resolution ET datasets have enhanced monitoring of extreme events (Pereira et al., 2015; Wu et al., 2021), alignment between climate model outputs and regional hydrological dynamics remains a critical challenge (Gimeno-Sotelo et al., 2024a). Consequently, the development of next-generation ET models that explicitly couple CO₂-climate feedbacks with regional hydrological dynamics is urgently needed.

The Coupled Model Intercomparison Project Phase 6 (CMIP6) provides an improved framework for addressing these challenges. Its long-term, multi-scenario structure and inclusion of dynamic vegetation modules offer greater fidelity in representing vegetation–climate coupling (Eyring et al., 2016; O'Neill et al., 2016). The Shared Socioeconomic Pathways (SSPs) embedded within CMIP6 allow for systematic exploration of ET trajectories under divergent emission scenarios (Zeng et al., 2016; Jones et al., 2016; Wu et al., 2019). However, PM-RC-CO₂ model proposed by Yang et al. (2019) still relies on formulas derived from CMIP5-era models, that the original coefficient ($\mathrm{2.4}\times {\mathrm{10}}^{-\mathrm{4}}$) may introduce inconsistencies when applied to CMIP6 scenarios. This can reduce confidence in cross-generational model comparisons and potentially exaggerate ET sensitivity to CO₂ rise. To fully leverage CMIP6's enhanced vegetation–climate framework, it is crucial to develop updated, scenario-specific formulations.

In this study, we advance ET modelling by integrating the CMIP6 climate projection system with the PM-RC-CO₂ model proposed by Yang et al. (2019). We use outputs from four global climate model (GCMs) from CMIP6, across four shared socio-economic pathways (from scenarioMIP, SSP1-2.6, SSP2-4.5, SSP3-7.0, SSP5-8.5) for the period 2015–2100. Meanwhile, we used monthly outputs of precipitation, actual evapotranspiration, short and long wave radiation, and near-surface air temperature, air pressure, wind speed and specific humidity, to derive an updated PM-RC-CO₂ Ep formulation. This framework improves the representation of interactions between stomatal resistance and aerodynamic drivers based on the characteristics of CMIP6 data. To assess performance, we compare the updated model against earlier formulations, evaluating Ep trends, scenario dependencies, and model outputs. The proposed framework provides a more robust characterization of Ep dynamics under the CMIP6 multi-scenario structure, thereby improving the simulation of future hydrological changes.

2 Data and methods

2.1 Data and model performance evaluation

To comprehensively assess the effects of CO₂ concentration on ET changes, we used outputs from four selected CMIP6 GCMs: IPSL-CM6A-LR, GFDL-ESM4, CNRM-CM6-1, and MPI-ESM1-2-HR. These models were obtained from the CMIP6 data portal (https://esgf-node.llnl.gov/search/cmip6/, last access: 30 August 2024) and include simulations for both the historical period (1850–2014) and four future emission scenarios (2015–2100): SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5.

Each model provides essential variables required for ET estimation, including monthly land surface data (runoff, precipitation, evapotranspiration, shortwave and longwave radiation) and near-surface atmospheric parameters (temperature, pressure, wind speed, and specific humidity). To ensure consistency across datasets, all outputs were resampled to a uniform spatial resolution of 0.25° × 0.25° by bilinear interpolation.

To evaluate the performance and reliability of the climate model simulations, we validated key variables against observational data from the W5E5 v2.0 dataset (https://data.isimip.org/search/simulation_round/ISIMIP3a/product/SecondaryInputData/climate_forcing/w5e5v2.0/, last access: 21 December 2024, 1979–2020). This comparison enabled the assessment of model skill in replicating observed climate and hydrological conditions, serving as a benchmark for subsequent ET modeling and scenario analyses.

The land use type data we used in this paper are from GLASS-GLC (http://www.glass.umd.edu/Overview.html, last access: 16 November 2024, period 1982–2015) (Liang et al., 2013, 2025).Compared with other land cover data sets, GLASS data has significant advantages in long-term time series and cross-year change monitoring, and can be widely used in ecological environment assessment, surface process simulation and global change research.

2.2 Determining r_s over non-water-limited regions and months

To quantify the sensitivity of stomatal resistance (r_s) under rising CO₂ concentrations, we applied the widely used non-water limited screening method (Milly and Dunne, 2016; Yang et al., 2019). This method systematically integrates hydrological constraints and eliminates the influence of frozen water through temperature constraints, so as to effectively screen out the water limited areas and periods.

The procedure is outlined as follows:

1. The period of analysis (1861–2100) was divided into eight 30-year climatological periods.

2. For each period, monthly ET and P outputs were fitted to a parabolic function. Grid cells with a maximum slope of the ET–P curve less than or equal to 0.05 were retained, indicating minimal hydrological limitation.

3. Based on existing research, we further filtered the data by selecting only those grid–month combinations where the ET/P ratio was less than 2.0, excluding regions under strong evaporative demand relative to precipitation.

4. The intersection of non-water-limited domains across all eight 30-year periods was retained to ensure consistent spatiotemporal coverage.

5. To eliminate the influence of frozen water, we exclude grid cells and months with average temperatures below 10 °C.

This filtering process isolates vegetated regions and time periods with minimal hydrological constraints, thereby enabling a more accurate assessment of stomatal resistance responses to elevated CO₂ concentrations.

2.3 Adjustment of the PM-RC-CO₂ model

2.3.1 The Penman–Monteith model

The Penman-Monteith (PM) equation provides a robust framework for estimating reference evapotranspiration (ET0), synthesizing surface energy balance and aerodynamic transfer principles (Monteith, 1977; Monteith and Unsworth, 2013; Milly and Dunne, 2016). The models calculate evapotranspiration (E) as:

$$\begin{array}{}\text{(1)}& \mathit{\lambda }E={\displaystyle \frac{s{R}_{\mathrm{n}}^{\ast }+{\mathit{\rho }}_{\mathrm{a}}{C}_{p}D/r_{\mathrm{a}}}{s+\mathit{\gamma }\left(\mathrm{1}+r_{\mathrm{s}}/r_{\mathrm{a}}\right)}}\end{array}$$

where s represents the slope of the saturation vapor pressure-temperature relationship (Pa K⁻¹), γ denotes the psychrometric constant (Pa K⁻¹), and ρ_a corresponds to air density (kg m⁻³). The specific heat at constant pressure (C_p, J kg⁻¹ K⁻¹) quantifies energy storage capacity. The vapor pressure deficit (D, Pa), drives evaporative demand, modulated by the temperature-dependent latent heat of vaporization (λ, J kg⁻¹). Aerodynamic resistance (r_a, s m⁻¹) is derived from logarithmic wind profile theory.

2.3.2 PM-RC Ep model

The standardized PM-RC model adopts fixed biophysical parameters representing typical C₃ crop physiology under non-drought conditions, including surface stomatal resistance (r_s=70 s m⁻¹), canopy height (0.12 m), and shortwave albedo (α=0.23), among others (Allen et al., 1998). Ep (mm d⁻¹) is computed as:

$$\begin{array}{}\text{(2)}& \text{Ep}={\displaystyle \frac{\mathrm{0.408}s{R}_{\mathrm{n}}^{\ast }+\mathit{\gamma }\frac{\mathrm{900}}{T+\mathrm{273}}uD}{s+\mathit{\gamma }\left(\mathrm{1}+\mathrm{0.34}u\right)}}\end{array}$$

2.3.3 PM-RC Ep model modified to account for atmospheric [CO₂](PM-RC-CO₂)

Yang et al. (2019) introduced a CO₂ responsive modification to the PM model by parameterizing stomatal resistance. This modification is based on empirically derived significant regression coefficients from controlled experiments, enabling quantification of vegetation-atmosphere feedbacks under high CO₂ conditions. This approach not only retains the biophysical basis of the original framework but also relates with the relationship between CO₂ concentration and stomatal dynamics. The equation is:

$$\begin{array}{}\text{(3)}& \text{Ep}={\displaystyle \frac{\mathrm{0.408}s{R}_{\mathrm{n}}^{\ast }+\mathit{\gamma }\frac{\mathrm{900}}{T+\mathrm{273}}uD}{s+\mathit{\gamma }\left\{\mathrm{1}+u\left[\mathrm{0.34}+\mathrm{2.4}\times {\mathrm{10}}^{-\mathrm{4}}\left(\left[{\mathrm{CO}}_{\mathrm{2}}\right]-\mathrm{300}\right)\right]\right\}}}\end{array}$$

Here, the term $\mathrm{2.4}\times {\mathrm{10}}^{-\mathrm{4}}$ ([CO₂] − 300) reflects the empirical modulation of stomatal resistance by CO₂ concentration, where 300 ppm represents the preindustrial baseline (1860–1960). This coefficient was obtained from nonlinear regression analysis of plant physiological responses in controlled environments, and encapsulates vegetation feedbacks aligned with CMIP6 scenario-driven CO₂ sensitivities.

2.3.4 Updated PM-RC-CO₂ model

The original PM-RC-CO₂ formulation was based on CMIP5-driven parameterization. Applying it directly to CMIP6 scenarios may introduce uncertainty due to differences in CO₂–climate feedback representation. Therefore, we recalibrated the stomatal resistance–CO₂ relationship using the CMIP6 multi-model ensemble, integrating vegetation physiological response and CO₂-forcing within a physically consistent parameter space. The updated formulation is:

$$\begin{array}{}\text{(4)}& \text{Ep}={\displaystyle \frac{\mathrm{0.408}s{R}_{\mathrm{n}}^{\ast }+\mathit{\gamma }\frac{\mathrm{900}}{T+\mathrm{273}}uD}{s+\mathit{\gamma }\left\{\mathrm{1}+u\left[\mathrm{0.34}+\mathrm{1.9}\times {\mathrm{10}}^{-\mathrm{4}}\left(\left[{\mathrm{CO}}_{\mathrm{2}}\right]-\mathrm{300}\right)\right]\right\}}}\end{array}$$

This version introduces a revised CMIP6-constrained empirical coefficient $\mathrm{1.9}\times {\mathrm{10}}^{-\mathrm{4}}$ [CO₂] − 300). Derived from a collection of four global climate models (GCMs), this parameter better captures the CO₂-induced stomatal resistance (r_s) effect and enhances model performance under the CMIP6 climate scenario.

2.4 Statistical Methods

To study changes in potential evapotranspiration (Ep) trends, we use the Bayesian Estimator of Abrupt Change, Seasonality, and Trend (BEAST) method (Zhao et al., 2019), which decomposes time-series into trend, seasonality, and residuals components and effectively handles nonlinear shifts (Li et al., 2022).

M-K trend test is a non parametric statistical test method (Mann, 1945). It can clearly detect the rise, fall or no trend in the time series data, and can more accurately capture the long-term trend without being disturbed by short-term fluctuations. It has been widely used in the field of climate change (Hamed, 2008; Machiwal et al., 2022). This study uses this method to analyze the global Ep results calculated by the original formula and the updated formula from 1850 to 2100.

We used Wilcoxon rank test to compare and analyze the results of different GCMS in different models at different periods. Wilcoxon signed rank test is a nonparametric statistical test used to compare the differences between two dependent variable samples (Cuzick, 1985). Like other nonparametric tests, this test does not assume a specific distribution of the analyzed data. The parameter test equivalent to Wilcoxon signed rank test is the dependent variable sample t test (or paired t test). If p<0.05, there is a significant difference between the two groups.

3 Results

3.1 Model performance evaluation and CO₂ Driven Surface Resistance Sensitivity

To assess the availability and reliability of climate model data, we compared four key climate variables (relative humidity, downward longwave radiation, downward shortwave radiation, and temperature) from the four global GCMs with observational data (Kling et al., 2012; Dahri et al., 2021; Zhang et al., 2024). Results indicate a strong agreement between model outputs and observations. For relative humidity (hurs), all models show high correlation coefficients (>0.85) and standard deviations within 10 % of observed values, indicating accurate humidity simulation (Fig. 1a). For the downward longwave radiation (rlds), model variability fits well with observations, with correlation coefficients above 0.75 and standard deviations mostly within 15 % (Fig. 1b). The downward shortwave radiation (rsds) shows slightly larger discrepancies, but still acceptable model performance (Fig. 1c), with correlation coefficients above 0.65 and standard deviations within 20 %. The simulation of temperature (tas) is particularly robust, with correlation coefficients exceeding 0.9 and standard deviations within 5 % (Fig. 1d). These results confirm that the selected GCMs are capable of accurately reproducing historical climatic conditions at the global scale.

Figure 1 Taylor diagrams of simulated climate variables compared to observations across the globe for the period 1979–2019. Panels show performance of four GCMs in simulating: (a) relative humidity (hurs), (b) downward longwave radiation (rlds), (c) downward shortwave radiation (rsds), and (d) near-surface air temperature (tas). The dashed green lines represent root mean square deviation (RMSD).

Using these validated models under multiple SSP scenarios, we evaluated the sensitivity of surface resistance (Δr_s) to changes in atmospheric CO₂concentration (Δ[CO₂]). A robust linear dependence was found between the two variables (Fig. 2), consistent across models. The inter model variability of the slope coefficient within ±15 % of the mean, demonstrating robustness. In historical simulations (1850–2014), the baseline CO₂ concentration is ∼284 ppm, with surface resistance around ∼52 s m⁻¹. Under the high emission SSP5-8.5 scenario (2071–2100), r_s rises to ∼78 s m⁻¹, while CO₂ increases to ∼935 ppm. This implies a ∼50 % increase in r_s for a ∼229 % increase in CO₂, revealing a stronger ET model sensitivity compared to earlier CMIP5-based projections.

Additionally, the relative sensitivity parameter Sr${}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$ remains spatially and temporally stable, ranging from 0.08 % ppm⁻¹ to 0.11 % ppm⁻¹ (Fig. 3). This confirms that CMIP6 models retain the CO₂ r_s coupling dynamics previously observed in CMIP5, but with increased sensitivity to CO₂.

Figure 2 Ensemble mean relationship between changes in stomatal resistance (rs) and atmospheric CO₂ concentration ([CO₂]) across four ESMs, relative to the 1851–1950 baseline, under different SSP scenarios. Panels (a)–(d) show IPSL-CM6A-LR results for SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5; (e–h) GFDL-ESM4 results for the same scenarios; (i–l) CNRM-CM6-1 projections under equivalent SSPs; and (m–p) MPI-ESM1-2-HR simulations following the SSP scenario hierarchy. All analyses use a consistent baseline climatology (1851–1950) and apply the same radiative scaling conventions.

Figure 3 Sensitivity of climate models to atmospheric CO₂ concentrations under different SSP scenarios. The horizontal axis shows the sensitivity metric $S_{\mathrm{rs}-\left[{\mathrm{co}}_{\mathrm{2}}\right]}$ (% ppm⁻¹), while the vertical axis lists the climate models. Colored bars represent different SSP scenarios. The dashed red line indicates the multi-model mean value of 0.1086.

Figure 4 Normalized annual mean Ep values across historical and future SSP scenarios. The black line represents the historical mean (1860–2014), while colored lines correspond to SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5. The shaded gray area denotes ±1 standard deviation.

3.2 Annual Changes in Ep from the Updated PM-RC-CO₂ Model

The updated PM-RC-CO₂ model, calibrated using CMIP6 data, shows a persistent upward trend in potential evapotranspiration (Ep) during historical simulations (1860–2014). By more accurately capturing nonlinear vegetation–climate interactions, the updated formulation improves the representation of CO₂-induced physiological feedbacks. This improvement is particularly evident under high-emission scenarios. From 2015 to 2100, the rates of Ep intensification differ across scenarios: under SSP5-8.5, the average decadal growth rate is approximately 2.1 %, whereas under SSP1-2.6, it stands at 1.2 %. Under SSP5-8.5, standardized Ep uncertainty (±1σ) rises from ±0.05 (2020–2040) to ±0.12 (2081–2100), indicating growing climate variability with increased CO₂ concentrations.

Figure 5 Time series decomposition of climate variables from 1850 to 2100 under different SSP scenarios: (a) SSP1-2.6, (b) SSP2-4.5, (c) SSP3-7.0, (d) SSP5-8.5. Each panel illustrates long-term trends, seasonal patterns, and residual variations within the scenario.

In order to further explore the characteristics of scenario related trend changes, we decompose the standardized annual Ep trend from 1850 to 2100. The low-emission SSP1-2.6 scenario shows gradual Ep increase of 58.3±14 mm yr⁻¹. The high-emission scenario (SSP5-8.5) exhibits a steeper Ep increase of 167.9±36 mm yr⁻¹. Seasonal decomposition shows ongoing cyclical patterns, but under SSP5-8.5, the amplitude of seasonal signals weakens after 2050, reflecting rising climate instability. Despite uncertainties, BEAST's probabilistic framework confirms robust stratification under SSPs, with SSP5-8.5 trends exceeding SSP1-2.6 by 187 % by 2100.

Figure 6 Global Ep trend patterns derived from the Mann-Kendall (MK) method using the original ET-RC-CO₂ model under (a) SSP1-2.6, (b) SSP2-4.5, (c) SSP3-7.0, and (d) SSP5-8.5. Associated histograms display the frequency distribution of Ep trend magnitudes across spatial domains.

3.3 Global Trends in Ep Calculated from Original and Updated PM-RC-CO₂ Models

The original ET-RC-CO₂ model Ep trend under different scenarios is shown in Fig. 6. Under the SSP1-2.6 scenario, the Ep trends show a relatively modest increase. This indicates that in this scenario, the changes in Ep are relatively stable although significant. However, when it comes to the SSP5-8.5 scenario, there are more pronounced increases in Ep, especially in tropical and mid-latitude regions. These areas experience a more substantial upward trend in Ep compared to other regions. The histograms in Fig. 6 further demonstrate that as scenarios become more extreme, there are higher frequencies of strong positive trends in Ep. This suggests that under the influence of different socioeconomic and climate scenarios, Ep exhibits diverse trends, with more extreme scenarios leading to more intense and frequent positive Ep trends.

Compared with the original model, the updated model using CMIP6 data (Fig. 7) has a similar spatial pattern, but there are significant differences in the trend amplitude. In SSP2-4.5 and SSP3-7.0, the trend magnitudes are consistently larger, indicating a more pronounced increase in Ep. This suggests that the updated model has an enhanced sensitivity to climatic shifts, as it captures more significant changes in Ep under these scenarios. The histograms also reflect this, showing a shift towards higher frequencies of strong positive trends in the more extreme scenarios.

Figure 8 compares the original and updated models. The spatial distribution of trend differences shows that the updated model predicts higher Ep trends by 2–3 mm yr⁻¹ in many regions under SSP3-7.0 and SSP5-8.5. This indicates an enhanced sensitivity of the updated model to climatic shifts, particularly in these more extreme scenarios. The differences are more pronounced in certain areas, suggesting that the updated model may better capture the regional variations in Ep trends.

Figure 7 Same as Fig. 6, but for the updated ET-RC-CO₂ model.

Figure 8 Differences in global Ep trends between the updated and original ET-RC-CO₂ models using CMIP6 data for (a) SSP1-2.6, (b) SSP2-4.5, (c) SSP3-7.0, and (d) SSP5-8.5. Associated histograms display the frequency distribution of Ep trend differences across spatial domains.

Regardless of the original or updated models, among the land use types such as Forest, Grassland and Tundra, the proportion of the area with MK trend greater than 0 (Slope >0) is generally very high (more than 95 %). These ecosystems are driven by climate change, and the Ep has a significant upward trend. The trend differentiation between Cropland and Shrubland is relatively obvious. For example, in the high emission scenario (SSP5-8.5), the proportion of farmland slope <0 in the updated model (1.8 %) is lower than that in the original model (5.3 %), and the proportion of shrub slope <0 in the updated model (3.9 %–9.5 %) is also different from that in the original model (8.3 %–18.6 %).

At the same time, updating the model effectively reduced the bias in trend estimation. The updated model also enhances the ability to capture CO₂ physiological feedbacks and scenario-specific responses, thereby reducing biases in trend estimation. Under different scenarios, the proportion of Ep rising trend of various land uses under high emission scenarios (such as SSP5-8.5) is more prominent than that under medium and low emission scenarios (such as SSP1-2.6). This indicates that emission intensity exerts an amplifying effect on the differentiation of evapotranspiration trends across different land use types.

Figure 9 Ep trend persentage across different land uses derived from the Mann-Kendall (MK) method using the original ET-RC-CO₂ model under (a) SSP1-2.6, (b) SSP2-4.5, (c) SSP3-7.0, and (d) SSP5-8.5. The pink part represents MK trend greater than 0, the blue and purple part represents less than zero, and the light green part represents equal to zero.

Figure 10 Same as Fig. 9, but for the updated ET-RC-CO₂ model.

3.4 Comparison Between Scenarios and Time Periods

We quantified Ep differences between original and updated models across scenarios (Fig. 11). Results show that the original model consistently underestimates Ep. Under SSP1-2.6, the updated formula's average Ep is 1.4 % higher than the original in the base period and 2.9 % higher in the late period. In SSP2-4.5, the difference widens from 1.3 % in the base period to 3.7 % in the late period. For SSP3-7.0, the late period difference reaches 5.0 %. In the high-emission SSP5-8.5, the late period difference is the largest at 6.3 %.

Moreover, we conducted a comparative analysis of the distributional results produced by the original and updated formulations across multiple models and time periods. Projections based on the original formulation (Fig. 12) indicate a marked increase in Ep severity from the historical baseline (1850–2014) to the late 21st century (2071–2100), with the biggest increases occurring after 2050. This upward trend persists under the updated formulation (Fig. 13), but the revised method exacerbates the projected Ep severity by an additional 12 %–18 %, particularly during the mid-century (2051–2070) and end-of-century under high-emission scenarios.

The comparative analysis (Figs. 12, 13) underscores the substantial sensitivity of Ep projections to the choice of model parameterization. Statistically significant differences (p<0.05) between the two methods were identified in 78 % of the late-century simulations. Notably, the divergence between original and updated formulations becomes increasingly pronounced with rising CO₂ concentrations, particularly under high-forcing scenarios.

In the late 21st century under SSP3-7.0, the updated formulation systematically yields higher Ep means compared to the original: 5.1 % increase (599.9 vs. 578.7 mm yr⁻¹) for IPSL-CM6A-LR, and 4.0 % (569.4 vs. 547.5 mm yr⁻¹) for GFDL-ESM4. Under SSP5-8.5, MPI-ESM1-2-HR exhibits a 6.8 % increase (919.8 vs. 861.4 mm yr⁻¹) by the end of a century, with its late-period Ep difference being 3.2 times larger than in the base period. Similarly, in CNRM-CM6-1, the Ep difference at the end of the century under SSP5-8.5 (758.2 vs. 715.4 mm yr⁻¹) doubles that of SSP1-2.6 (773.8 vs. 751.9 mm yr⁻¹) in the same period.

Overall, the updated formulation consistently produces higher Ep projections across all emission scenarios, with differences becoming more pronounced over time and with increasing emission intensity.

Figure 11 Total annual Ep severity using both the original and updated ET-RC-CO₂ formulas across different SSP scenarios and time periods: base (1850–2014), early (2015–2050), middle (2051–2070), and late (2071–2100). Panels: (a) SSP1-2.6, (b) SSP2-4.5, (c) SSP3-7.0, (d) SSP5-8.5. Significance levels: *** for p≤0.001, ** for p≤0.01, and * for p≤0.05.

Figure 12 Total annual Ep severity using the original ET-RC-CO₂ formula across SSP scenarios and time periods: base (1850–2014), early (2015–2050), middle (2051–2070), and late (2071–2100). Panels: (a) SSP1-2.6, (b) SSP2-4.5, (c) SSP3-7.0, (d) SSP5-8.5. Significance levels: *** for p≤0.001, ** for p≤0.01, and * for p≤0.05. The upper green sign indicates p values between two GCMs; the lower sign indicates p value between all GCM pairs, except the lower green one.

Figure 13 Same as Fig. 12 but for the updated ET-RC-CO₂ formula.

4 Discussion

The observed increase in Ep severity is consistent across multiple analytical methods, supporting earlier findings by Yang et al. (2023). Nevertheless, our results reveal that numerical estimates remain highly sensitive to model structure and parameterization schemes. Notably, evapotranspiration projections display heightened sensitivity during the latter half of the 21st century, especially under high-emission scenarios. This emphasizes the added value of incorporating CMIP6 data, which improves model responsiveness to evolving climatic drivers (Ma et al., 2018; Liu et al., 2022; Bai et al., 2025). These findings underscore the need for comprehensive, well-calibrated models in future climate risk assessments. A key insight from this study is the significant role of CO₂-induced plant physiological responses in modulating evapotranspiration. To improve the robustness and predictive capability of future hydrological models, better integration of dynamic vegetation processes is essential.

4.1 Uncertainties Across Emission Scenarios

Substantial uncertainty persists in the current understanding and modelling of the ET process (Pan et al., 2020). As shown in Fig. 11, low-emission scenarios (such as SSP1-2.6) exhibit a more dispersed distribution across models than high-emission scenarios, leading to slightly elevated average Ep values under low-emission conditions.

To explore these uncertainties, we performed comparative analysis using ET models driven by the same forcing data. Distributions from both the original and updated formulations were evaluated across models and time periods (Figs. 10 and 11), helping to pinpoint sources of model uncertainty and guiding future improvements in ET estimates (Warszawski et al., 2014; Miralles et al., 2016).

Our results indicate that low- and medium-emission scenarios show greater inter-model variability. Moreover, the original model consistently underestimates long-term climate responses under medium- to high-emission scenarios. This underestimation may be attributed to a parameterization process that heavily relied on historical high-emission conditions, limiting it's ability to capture future feedback mechanisms under intensifying radiative forcings (Knowling et al., 2019; Cui et al., 2021).

Although the updated formulation demonstrates improved sensitivity to CO₂ variations, there is still room for improvement – particularly under low-emission scenarios. These trends are shaped by varying driving factors, model architectures, and differences in parameterization approaches across existing ET datasets (Badgley et al., 2015; Michel et al., 2016). Future work should focus on reducing systemic bias in projections under low-emission scenarios.

4.2 Impacts, Limitations, and Prospects

This study relied on unadjusted CMIP6 outputs, which may introduce systematic errors when applied to the improved ET formulation. Future work could benefit from using bias-corrected CMIP6 datasets and accounting for systematic discrepancies in atmospheric forcing variables. Previous studies have shown that calibrated CMIP6 products can significantly enhance the representation of precipitation seasonality and CO₂ trajectories (Raziei and Pereira, 2013a, b). Such improvements are especially critical in semi-arid regions, where more refined modeling can incorporate nonlinear feedbacks between stomatal conductance and atmospheric CO₂, as well as biochemical parameterization for C₃/C₄ species (Ding et al., 2013; Potkay et al., 2025; Wu et al., 2025).

The current model assumes fixed vegetation responses to CO₂, which is appropriate for large-scale analyses but less suitable for regional applications (Bao et al., 2021; Cui et al., 2023). Given the high specificity of vegetation types at regional scales, fixed-response models struggle to capture physiological variations among species. For finer-scale assessments, model accuracy can be enhanced by integrating detailed representations of vegetation physiological responses, such as dynamic correlation data between photosynthetic rate and CO₂concentration (Luo et al., 2018a, b). For example, dynamic stomatal conductivity models that respond to rapid CO₂ and light fluctuations can increase transpiration estimates by 22 %–30 % under extreme conditions (Lawson and Vialet-Chabrand, 2019; Poyatos et al., 2016; Poyatos et al., 2021). The use of canopy conductance algorithms that account for photosynthetic pathway differences is also critical for accurately simulating ET in C₄-dominated dryland ecosystems (Croft et al., 2017; Wei et al., 2019).

Advanced parameterization, such as linking foliar nitrogen and phosphorus content to photosynthetic efficiency, enables more ecosystem-specific transpiration estimates (Cernusak et al., 2010). When combined with vegetation indices like Normalized Differences Vegetation Index (NDVI) (Suarez et al., 2008; Sayago et al., 2017; Ballester et al., 2018), these refinements can reduce simulation error by up to 0.38 mm d⁻¹ compared to FAO-56 methods (Allen et al., 2007; Pereira et al., 2021).

Integrating structural vegetation indices such as Leaf Area Index (LAI) and canopy height gradient, further improves the soil vegetation water coupling. Recent advances, such as the use of GEDI-derived canopy height products, have shown strong potential in separating soil evaporation from canopy transpiration via aerodynamic resistance calibration – particularly valuable in arid regions (Bao et al., 2022; Chen et al., 2022; Wu et al., 2024). Additionally, using the functional trait diversity threshold for grid selection has reduced uncertainty in carbon water coupling prediction in forest ecosystems by 12 %–15 % (Li et al., 2013; Xu et al., 2015; Zhang et al., 2018; Joswig et al., 2022; Wang et al., 2022).

These developments align with the improved Budyko framework, which integrates vegetation-mediated runoff elasticity and helps mitigate scale-dependent biases in regional ET partitioning (Roderick et al., 2014; Zeng et al., 2016; Mianabadi et al., 2019; Yang and Roderick, 2019; Scheff et al., 2022). Finally, improved spatial screening through moisture-limited grid selection criteria has been shown to reduce spatial heterogeneity in ET projections (Talsma et al., 2018; Lian et al., 2018; Lian et al., 2021).

5 Conclusion

This study integrates CMIP6 multi-scenario projections and CO₂-vegetation coupling effects into the Penman-Monteith evapotranspiration (PM-ET) model, enhancing the representation of climate-vegetation interactions for improved hydrological projections. Methodological advances include the dynamic parameterization of stomatal resistance and scenario-dependent sensitivity analysis. The main findings are as follows:

1. The updated PM-ET model reduces CO₂-induced evapotranspiration (ET) bias by 15 %–20 % compared to earlier approaches, showing improved consistency with CMIP6 simulations.

2. Evapotranspiration potential (Ep) exhibits a consistent upward trend – particularly under high-emission scenarios such as SSP5-8.5, where increases reach up to 635.1 mm yr⁻¹ – driven by CO₂-climate synergies.

3. Model parameterization plays a critical role in capturing CO₂-physiological feedbacks; however, regional uncertainties remain due to heterogeneity in vegetation responses and methodological sensitivities – for instance, the distinct differences in responses between croplands and shrublands. Additionally, regional uncertainties persist due to heterogeneous vegetation responses and methodological sensitivities..

These findings highlight the need to further refine hydrological models by incorporating CMIP6-specific mechanisms, such as dynamic vegetation modules and biome-specific feedbacks. By improving the accuracy of scenario-based projections, this work contributes to more robust assessment of water resource risks under climate change and provides valuable insights for adaptation planning in both ecological and agricultural systems.

Appendix A: Variable declaration

λ	temperature-dependent latent heat of vaporization (J kg−1)
s	gradient of the saturation vapour pressure with respect to temperature(Pa K−1)
$R_{\mathrm{n}}^{\ast }$	available energy (MJ m−2 d−1)
ρa	air density (kg m−3)
Cp	specific heat at constant pressure(J kg−1 K−1)
rs	surface stomatal resistance (s m−1)
ra	Aerodynamic resistance (s m−1)
D	vapor pressure deficit (Pa)
γ	psychrometric constant (Pa K−1)
u	wind speed (m s−1)
T	air temperature (°C)
α	shortwave albedo