The aim of this paper is to demonstrate the extent of disparity between reference evapotranspiration (ET_o) estimated with and without incorporating the influence of CO₂ on contiguous India during three timeframes: the near-term (2021–2030), mid-term (2051–2060) and the long-term (2091–2100) periods. To achieve this, data from five Global Climate Models (GCMs) (see Supplementary Table S1) were obtained from the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) (Warszawski et al., 2014; Hempel et al., 2013). These five models were chosen because their GCM projections were bias-corrected for the systematic deviation from observations and made freely accessible through the ISIMIP portal of the Potsdam Institute of Climate Impact Research.¹ Furthermore, ISIMIP data effectively capture the uncertainties in global temperature change projections (Ito et al., 2020), which is essential since temperature is an important variable in estimating ET_o.

We focus on two contrasting climate change scenarios: SSP5-8.5 and SSP1-2.6 (see Supplementary Figure S1), which correspond to projected global CO₂ concentrations of approximately 1,130 ppm and 474 ppm, respectively (Büchner and Reyer, 2022). These scenarios are associated with a projected mean temperature rise in India of 4.0-to-4.4 °C under SSP5-8.5 and 1.2-to-1.8 °C under SSP1-2.6 by 2100. We selected SSP5-8.5 because it is commonly used in climate change impact assessment under the worst-case, high-emissions scenario (Jaiswal et al., 2017; Pielke, 2021; Climatedata.ca, 2024), which closely followed observed CO₂ emission trends until recent years (Fuss et al., 2014). In contrast, SSP1-2.6 represents a low-emissions, sustainable development pathway, serving as a benchmark for the most optimistic future with aggressive mitigation.

We employed two approaches to estimate ET_o, one without considering the effect of atmospheric CO₂ concentration (Equation 1) and the second after incorporating atmospheric CO₂ concentration (Equation 6). The first approach is based on the FAO-PM equation (Allen et al., 1998), which combines the aerodynamic component with the energy component and is idealized for a hypothetical reference crop (Allen et al., 1998). The second approach modifies FAO-PM equation by considering stomatal conductance as a function of atmospheric CO₂ (Equation 5) (Li et al., 2019) instead of a fixed value of surface resistance of 70 s m⁻¹ as used in the original FAO-PM equation (Equation 1) (Allen et al., 1998).

According to the original FAO-PM equation (Allen et al., 1998),

where, ET o original is the reference evapotranspiration (mm day⁻¹), R_n is the net radiation (MJ m⁻² day⁻¹) at the canopy surface, G is the soil heat flux density (MJ m⁻² day⁻¹) (G is assumed negligible and hence equals zero), γ is the psychrometric constant (kPa °C⁻¹), T is mean daily air temperature (°C) at 2 m height, u₂ is the wind speed at 2 m height (m s⁻¹), e_s is the saturation vapor pressure (kPa), e_a is the actual vapor pressure (kPa), e_s - e_a is vapor pressure deficit (kPa), and Δ is the slope of the saturated vapor pressure curve (kPa °C⁻¹).

During the formulation of Equation 1, the term ( 1 + r s r a ) from the Penman-Monteith (PM) model (Monteith, 1965) is substituted with r_s = 70 s m⁻¹ and r a = 208 u 2 , to obtain (1 + 0.34u₂). As a first step to incorporating CO₂ in FAO-PM equation, we modify Equation 1 in the following way (Allen et al., 1998; Li et al., 2019; Jarvis et al., 1997):

where, r_s is the bulk surface resistance (s m⁻¹), r_a is the aerodynamic resistance (s m⁻¹), g_c is the canopy conductance (m s⁻¹), g_s is the leaf stomatal conductance (m s⁻¹), LAI_active is effective leaf area index (m² m⁻²), h is the hypothetical crop height (equals 0.12 assumed in Equation 1 in m).

Substituting Equations 2, 3 in Equation 1, we get the modified FAO-PM model (Equation 4)

In the above equation, stomatal conductance g_s (m s⁻¹) appears on the right-hand side of Equation 4. Li et al. (2019) developed a modified hyperbolic model that express g_s as a function of atmospheric CO₂ as shown in Equation 5.

We replaced g_s (m s⁻¹) from Equation 5 to Equation 4 to obtain CO₂ dependent value of reference evapotranspiration ( ET o CO 2 ) (in mm day⁻¹). Our modified version of the FAO-PM equation, presented in Equation 6, provides a simplified representation of ET o CO 2 as a function of atmospheric CO₂, and all the other meteorological inputs used in the original FAO-PM equations (Allen et al., 1998).

To validate our proposed modified FAO-PM equation, we used measured data from six independent sites included in the AmeriFlux network (Novick et al., 2018), which is a part of the global FLUXNET (Pastorello et al., 2020) network of eddy covariance towers.² The six sites were selected based on details provided about vegetation and/or crop cover, availability of crop coefficients, and availability of the planting and harvest dates (see Supplementary Table S2; Nass, 2010). They comprise of four agricultural and natural vegetation types including alfalfa cultivation (Twitchell and Bouldin Islands), managed pastures (Medford hay pasture), irrigated croplands (continuous maize at Mead) and native prairie ecosystems (Rosemount Prairie and Konza prairie; see Supplementary Table S2). All of the sites also span across different types of climates. Twitchell alfalfa and Bouldin islands come under the Csa Koppen climate classification (Mediterranean) with mild winters and dry hot summers. Medford hay pasture and Konza Prairie exhibits Cfa climate (Humid subtropical) with mild winters, hot summers and year-round rainfall. Rosemount prairie and Mead’s maize site experience a Dfa climate (Humid Continental) with severe winters, hot summers, and no dry season. The downloaded data from the AmeriFlux site consisted of the following variables; air temperature, vapor pressure deficit, net radiation, wind speed, atmospheric pressure, atmospheric CO₂ and latent heat flux. Additionally, weather data from two representative stations located in the southern Indian state of Kerala were used to compare model performance: one at Trivandrum (8.544°N, 76.913°E) and the other at Palakkad (10.807°N, 76.7258°E). These datasets were used to compare the performance between the modified FAO-PM equation and two other methods: the original FAO-PM equation (Allen et al., 1998) and the widely used Priestley-Taylor equation (Priestley and Taylor, 1972).

We converted the daily observed latent heat flux data into actual evapotranspiration by dividing latent heat flux with latent heat of vaporization [λ = 2.45 MJ kg⁻¹ (Allen et al., 1998)] (see Supplementary Equation S1). The actual evapotranspiration was then divided by crop coefficients (see Supplementary Table S2) corresponding to vegetation type (see Supplementary Equation S2) to estimate reference evapotranspiration ( ET o obs ). We used ET o obs to validate our predictions of ET o CO 2 made using modified FAO-PM equation (Equation 6) for the six sites.

To compare our model estimations against commonly used approaches for selected sites in India, we used the weather data from specific stations in the South Indian state of Kerala (see Supplementary Table S3) to estimate reference evapotranspiration using both the original and modified FAO-PM models. Additionally, the ET o CO 2 values were further compared with the estimated values by the Priestley-Taylor equation, based on the same weather station data.

The statistics used for the validation and comparison were Root Mean Square Error (RMSE in mm day⁻¹), coefficient of determination (R²), and the correlation coefficient (r) (see Supplementary Equations S3–S5).

The climate data to estimate ET_o across India was obtained from the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP), ISIMIP3b protocol (Hempel et al., 2013; Warszawski et al., 2014). The ISIMIP project provides bias-corrected gridded global climate projection data from 1800 to 2100 at daily (and coarser) time steps with a spatial resolution of 0.5° × 0.5°. We used projected climate data for contiguous India for the period from 2021–2030 (near-term), 2051–2060 (mid-term), and 2091–2100 (long-term) corresponding to the five GCMs, namely GFDL-ESM4, IPSL-CM6A-LR, MPI-ESM1-2-HR, MRI-ESM2-0, and UKESM1-0-LL under SSP1-2.6 and SSP5–8.5 scenarios (Lange and Büchner, 2021; see Supplementary Table S1). We downloaded six variables, namely near surface relative humidity, surface air pressure, surface downwelling shortwave radiation, near surface wind speed, daily maximum near surface temperature, and minimum near surface air temperature (Lange and Büchner, 2021). Additionally, ISIMIP3b atmospheric composition input data for annual mean CO₂ concentrations under SSP5-8.5 and SSP1-2.6 (Büchner and Reyer, 2022) was also downloaded (see Supplementary Figure S1) to calculate ET o CO 2 using modified FAO-PM equation (Equation 6). The soil heat flux density (G) was assumed to be negligible in our calculations (Allen et al., 1998; Varghese and Mitra, 2024).

We estimated the daily ET o original (Equation 1) and ET o CO 2 (Equation 6) for each year from 2021–2030 (near-term), 2051–2060 (mid-term) and 2091–2100 (long-term) across contiguous India using data from each of the five GCMs for both the scenarios (SSP1-2.6 and SSP5-8.5). From this point forward, any reference to India refers to contiguous India.

Maps were classified using manually determined class breaks after identifying the minimum and maximum values projected by all GCMs for each step (ESRI, 2024).

To quantify relative impact of incorporating CO₂ on national average evapotranspiration we have calculated ET o , near − term original , ET o , mid − term original and ET o , long − term original , which represent spatio-temporal averages (for entire country over a period of 10 years) of ET_o estimated using the original FAO-PM (Equation 1) for near-term period (2021–2030), mid-term (2051–2060) and long-term periods (2091–2100), respectively. Similarly, we calculated ET o , near − term CO 2 , ET o , mid − term CO 2 and ET o , long − term CO 2 , which represent the spatio-temporal averages of ET o CO 2 estimated using the modified FAO-PM (Equation 6). We calculated percentage change in ET o original , as we move from near-term to mid-term ( 100 × ( ET o , near − term original – ET o , mid − term original ) ÷ ET o , near − term original ) , and long-term ( 100 × ( ET o , near − term original – ET o , long − term original ) ÷ ET o , near − term original ) with the assumption that original FAO-PM equation would be continued to be used. Similarly, we also calculated percentage changes in the estimated ET_o as a result of incorporating CO₂ for near-term ( 100 × ( ET o , near − term original – ET o , near − term CO 2 ) ÷ ET o , near − term original ), mid-term ( 100 × ( ET o , mid − term original – ET o , mid − term CO 2 ) ÷ ET o , mid − term original ), and long-term periods ( 100 × ( ET o , long − term original – ET o , long − term CO 2 ) ÷ ET o , long − term original ) using ET o original as base value.

After the spatio-temporal analyses, we identified two locations in India which exhibit drastically different response of rising CO₂ and temperature on ET_o. These two locations were used to conduct a sensitivity analysis of ET o CO 2 with respect to temperature and atmospheric CO₂. To perform the sensitivity analysis, the CO₂ concentration was varied from 400 to 1,200 ppm while keeping all the remaining input parameters unchanged. Similarly, temperature sensitivity analysis was done by varying temperature from 18 to 34 °C while keeping all the remaining input parameters unchanged. These ranges were determined based on meteorological data for the three decadal periods for these two specific locations. The values of all the other variables during sensitivity analysis were kept constant based on their average values and are given in Supplementary Table S4. All the data produced during this analyses is available publicly as an archive (Surendran et al., 2025).

We developed the modified FAO-PM (Equation 6) which effectively predicted the daily reference evapotranspiration ( ET o CO 2 ) at the six AmeriFlux (Novick et al., 2018) sites (Figure 2) with atmospheric CO₂ concentrations ranging from 370 ppm in 2001 to 418 ppm in 2021 (see Supplementary Table S2). The correlation coefficient (r) indicates moderate to strong linear relationships between the observed ( ET o obs – see Supplementary Equation S2) and estimated ET_o using Equation 6 ( ET o CO 2 ), with the site US-A32 showing the highest correlation (r = 0.72) and the site US-Ne1 showing the lowest correlation (r = 0.606). The slopes of the regression lines through the origin ranged from 0.8 to 1.04, demonstrating close alignment with the 1:1 line across all six sites (Figure 2). RMSE values ranged from 1.152 to 1.566 mm day⁻¹, with the site US-Ro4 exhibiting the lowest RMSE and the site US-Bi1 showing the highest.

To compare the estimation of the modified FAO PM equation (Equation 6) against commonly used methods for selected Indian sites, we used daily weather data from two stations located in the southern Indian state of Kerala: Trivandrum (8.544°N, 76.913°E) and Palakkad (10.807°N, 76.7258°E). Using this data, reference evapotranspiration was estimated using the original FAO-PM model (ETₒ FAO–PM), the modified CO₂-sensitive FAO-PM model ( ET o CO 2 ), and the Priestley-Taylor model (ET_PT). Figure 3 compares ET o CO 2 with ET_o FAO–PM and ET_PT for both locations. The modified FAO-PM model exhibited a strong agreement with the original FAO-PM model at both sites, with high R² values (0.999 for Trivandrum and 0.988 for Palakkad), strong correlation coefficients (r = 0.996 and 0.992, respectively), and low RMSE values (0.084 and 0.329-mm day⁻¹, respectively). Comparisons with the Priestley-Taylor model showed relatively lower agreement, with higher RMSE values (0.486 and 0.887 mm day⁻¹, respectively) and underprediction tendencies (slopes of 0.82 and 0.76, respectively).

Analyses of ET o original and ET o CO 2 under SSP1-2.6 and SSP5-8.5 scenario utilizing five GCMs (see Supplementary Table S1; Supplementary Figure S1) during near-term (2021–2030), mid-term (2051–2060) and long-term (2091–2100) periods, and all spatial locations in India, showed that the former exceeded the latter consistently in all cases (Figure 4) due to the impact of rising atmospheric CO₂ concentration that was not included in the original FAO-PM model (Equation 1). In general, both ET o original and ET o CO 2 , are decreasing as we move from west to east for all cases (Figure 4). The highest values of ET o original and ET o CO 2 (Figure 4) and the difference between ET o original and ET o CO 2 (Figure 5) were observed covering parts of desert regions in Rajasthan by all the five GCMs for all the time periods and scenarios.

While the differences between ET o original and ET o CO 2 were similar for both SSP1-2.6 and SSP5-8.5 in the near term, they became significantly larger under SSP5-8.5 during the mid-term and long-term time periods, indicating divergent temporal trajectories across the two scenarios (Figure 5). The regions predicted to have lowest values of ET o original and ET o CO 2 (Figure 4) and their difference (Figure 5) included high-altitude deserts of Ladakh and northeastern states.

The spatially averaged values of ET o original and ET o CO 2 across India (Figure 6) indicated a small variation in the time-series of estimated values within a period of 10-years corresponding to the near-term, mid-term and long-term time periods for all the GCMs and under both SSP1-2.6 and SSP5-8.5 scenarios. The predicted values of ET o original was consistently higher than ET o CO 2 for all the scenarios and time periods. The differences between the two were largest for the long-term (2091–2100) period under SSP5-8.5 when CO₂ concentrations are projected to reach 1,130 ppm (see Supplementary Figure S1).

Spatial variations in the magnitude of the differences between ET o original and ET o CO 2 for the four seasons are seen in all three time periods under both the scenarios (Figure 7). The differences between ET o original and ET o CO 2 was the largest during the pre-monsoon season, in most parts of India except the far-northern and north-eastern states. The lowest differences were observed during the post-monsoon and winter season (Figure 7). Spatially, the magnitude of these differences for the winter and post-monsoon seasons were the least in the northern and north-eastern regions and highest in major parts of Rajasthan, Gujarat, Maharashtra and southern India for all time periods and scenarios. The differences were higher for the long-term period than for the near-term period and mid-term period only under scenario SSP5-8.5 for all the locations and seasons. Under the SSP1-2.6 scenario, the differences were smaller and ranged from 0.05 to 2.5 mm day⁻¹ across all seasons and time periods.

The absolute spatio-temporal average values of the seasonal (monsoon, post-monsoon, winter, and pre-monsoon) variations of ET o original and ET o CO 2 for the near-term, mid-term and long-term periods under SSP1-2.6 and SSP5-8.5 are shown in Table S5 and S6, respectively (see Supplementary material) for the five GCMs. The seasonal patterns are nearly consistent among GCMs with higher values of ET o original and ET o CO 2 observed during pre-monsoon and lower values observed during the post-monsoon and winter season for all the three simulation periods and scenarios. All seasons show a decrease in the predicted ET_o when incorporating the effect of CO₂ for SSP5-8.5 (see Supplementary Tables S5, S6).

The sensitivity analysis of ET o CO 2 with respect to temperature and CO₂ concentration (Figure 8) revealed distinct trends at both the locations selected. Location 1 (26.75°N, 70.25°E), situated in the arid desert region of Rajasthan with a hot desert climate, exhibited ET o CO 2 values ranging from 3.22 to 9.65 mm day⁻¹ as CO₂ concentration increased from 400 to 1,200 ppm and temperature was varied from 26 to 33 °C. This location showed the greatest difference in ET_o​ estimated using original and modified FAO-PM equations. Location 2 (23.75°N, 93.25°E), situated in Arunachal Pradesh with a tropical rainforest climate, showed ET o CO 2 values ranging from 0.6 to 1.35 mm day⁻¹ for the same CO₂ range and temperature variation from 18 to 24 °C. This location exhibited the lowest difference in ET_o estimated using original and modified FAO-PM equations.

Spatio-temporal averages across India (Figure 9) shows that incorporating CO₂ concentrations in calculating ET_o results in a reduction of 12.4, 13.9, and 13.0% for near-term, mid-term, and long-term time periods respectively under SSP1-2.6. Under SSP5-8.5, the percentage difference (12.6%) was more-or-less similar to the ones observed under SSP1-2.6 for the near-term period but were much larger for the mid-term and long-term periods (18.0 and 29%, respectively).

If one continues to use original FAO-PM equation (Equation 1), which does not consider the effect of rising CO_2, then the increase in the ET_o is expected to be 4.0 and 4.6% as we move from near-term (2021–2030) to mid-term (2051–2060), and near-term to long-term (2091–2100) periods, respectively, under scenario SSP1-2.6. These numbers change to 6.0 and 16.0%, respectively under the scenario SSP5-8.5.

Our model performed satisfactorily against the AmeriFlux network data (Novick et al., 2018), with R² values ranging from 0.75 to 0.90 (Figure 2). These values are comparable to those reported by Li et al. (2019), who found R² values of 0.76 to 0.83 while validating ET o CO 2 against water balance-based estimates for maize under controlled conditions. In contrast, our validation included six sites with diverse vegetation types and management practices (see Supplementary Table S2), where site selection was guided by the availability of crop coefficient (Basketfield, 1985; Allen et al., 1998; Nass, 2010; Pereira et al., 2023), which is essential for converting actual evapotranspiration (ET_actual) (see Supplementary Equation S1) from eddy covariance towers to ET o obs using crop coefficients (see Supplementary Equations S1, S2). The greater scatter observed in our validation plot (Figure 2) can be attributed to uncertainties in crop coefficients (Peng et al., 2019), as well as site-specific factors such as spatial heterogeneity, variations in planting and harvest dates, irrigation practices, environmental stresses, and other management practices.

The comparison of the modified FAO-PM model using daily weather data from two locations in Kerala, Palakkad and Trivandrum, showed strong agreement with the original FAO-PM formulation, yielding high R² values (0.98 to 0.99) (Figure 3). These sites represent typical humid tropical environments characterized by high moisture availability and seasonal variability, making them suitable test references for model evaluation in such climates. Palakkad includes agriculturally important areas, further supporting the relevance of these results. The high agreement with the original FAO-PM method suggests that the modified model preserves the core structure and reliability of the standard formulation while integrating the physiological response of vegetation to elevated CO₂. Previous studies, such as Nandagiri and Kovoor (2006) have shown that radiation-based models like Priestley-Taylor (PT) perform reasonably well in humid regions. Supporting this, our comparison of ET o CO 2 with ET_PT at both Kerala locations showed strong statistical relationships (R² = 0.936-to-0.966; r = 0.919-to-0.988), indicating that PT captures the temporal patterns of evapotranspiration well. However, the slopes of the lines (0.76 and 0.82) and relatively higher RMSE values (0.486 to 0.887 mm day⁻¹) point to a consistent underestimation by the PT method compared to the modified FAO-PM model (Figure 3). This underprediction highlights the advantage of including physiological and aerodynamic controls, as well as CO₂ sensitivity, which are absent in simpler radiation-based models.

Together, these results show that the modified FAO-PM model not only performs reliably under controlled or semi-controlled conditions but also maintains robustness across complex, real-world scenarios. This consistency across different climates, vegetation types, and data sources supports the model’s potential for large-scale applications in climate impact studies and agricultural water management.

The CO₂ range in the validation data for the modified FAO-PM equation (Equation 6) was relatively narrow (370–418 ppm) (see Supplementary Table S2), reflecting past natural environmental conditions. However, this modified FAO-PM equation (Equation 6) has been previously validated under controlled conditions with CO₂ levels up to 900 ppm (Li et al., 2019). This higher level is comparable to the atmospheric CO₂ concentrations projected by the SSP5-8.5 scenario in the long-term period (2091–2100) (see Supplementary Figure S1), considered in our study to evaluate the spatial variation of the impact of CO₂ on ET_o across India (Figures 4b, 5b).

We predicted a consistent trend of ET o CO 2 being less than ET o original (Figures 4, 5) due to the CO₂ impacts on stomatal closure, which is similar to reported trends in previous studies (Li et al., 2019; Yang et al., 2019; Varghese and Mitra, 2024). In our knowledge, the absolute values of ET o CO 2 and ET o original covering whole India have not been reported previously. The predicted ET o original for the near-term period ranged from 1.94-to-7.04 mm day⁻¹ under both SSP1-2.6 and SSP5-8.5 (Figure 4), aligning with recent estimates done for smaller regions within India (Jhajharia et al., 2009; Nag et al., 2014; Pandey et al., 2016; Das et al., 2023). The increase in ET o original from the near-term to long-term (Figure 4b) can be attributed to rising temperatures under SSP5-8.5, as studies have shown a positive correlation between temperature and ET_o (Wang et al., 2022; Zhou et al., 2022) with temperature contributing up to 45% of ET_o variation (Varghese and Mitra, 2024). However, this increase in ET_o is moderated ( ET o CO 2 < ET o original in Figure 4) as a consequence of incorporating the effect of CO₂ in our calculations. Estimations using all five GCMs (Figures 4, 5) are consistent in predicting the ET o CO 2 to be less than the ET o original but they differ in the magnitude as well as spatial variations of the differences. The spatial averages (Figure 6) reveal no distinct temporal trends within 10-year intervals for either the near-term or long-term periods across all GCMs. However, transitioning from the near-term to long-term period highlights the dominant influence of CO₂ on ET_o, with reductions in ET_o due to elevated CO₂ levels outweighing increases driven by rising temperatures under SSP5-8.5 (Figure 6b). Consequently, long-term ET o CO 2 is projected to be lower than the current estimates, under the SSP5-8.5 scenario, contrary to several previous studies (Liu et al., 2020; Zhai et al., 2020) that did not account for the impact of CO₂.

The spatial variability of differences between ET o CO 2 and ET o original (Figure 5) is influenced by climate inputs beyond CO₂ as shown in the sensitivity analyses of the modified FAO-PM (Equation 6) for two distinct locations (Figure 8). Location 1 (26.75°N, 70.25°E), characterized by the hot desert climate in the northwest, exhibits a pronounced response to rising CO₂, with ET o CO 2 ranging from 3.22-to-9.65 mm day⁻¹, as CO₂ concentration increases from 400 to 1,200 ppm and temperature varies from 26-to-33 °C. In contrast, Location 2 (23.75°N, 93.25°E), characterized by a tropical rainforest climate of the northeast, shows minimal response, with ET o CO 2 values varying from 0.6-to-1.35 mm day⁻¹ for the same CO₂ range and temperature variation from 18-to-33 °C. This non-linear, location-specific interaction between CO₂, temperature, and ET_o may help explain deviations from the typically positive correlation between ET_o and temperature (Wang et al., 2022; Zhou et al., 2022), a phenomenon often referred to as the “evapotranspiration paradox” (Rao and Wani, 2011; Varghese and Mitra, 2024). The “evapotranspiration paradox” may result from overly simplistic vegetation representation in hydrological models, despite the fact that leaf stomatal transpiration can account for over 80% of evapotranspiration (Nelson et al., 2020; Yu et al., 2024). Similar behavior was observed by Vremec et al. (2024) in the Austrian Alps. Consequently, the effects of climatic factors such as CO₂, temperature, humidity, wind speed, and radiation on leaf stomatal behavior are often unappreciated (Ainsworth and Long, 2021) and continue to remain a challenge in the field of hydrology (Blöschl et al., 2019). However, substantial uncertainties also remain in predicting these variables (temperature, humidity, wind speed, and radiation), making it essential to select GCMs based on performance indicators tailored to specific regions (Raju and Kumar, 2020). Unfortunately, most studies evaluating the suitability of GCMs have focused on specific regions within India (Song et al., 2023; Verma et al., 2023). Panjwani et al. (2019) have covered all of India but could not identify a single GCM capable of reliably predicting all variables required for evapotranspiration calculations. Consequently, it is challenging to determine which of these models (see Supplementary Table S1) is best for ET predictions across India. However, multimodal ensemble methods are often preferred over single-model predictions (Khan et al., 2018), as demonstrated in several studies on issues related to water resources under climate change (Haddeland et al., 2011; Davie et al., 2013).

Incorporating changing CO₂ concentrations in the mid-term period (2051–2060) allows us to assess not only the long-term implications of elevated CO₂ but also the potential transitional effects that may influence water demand and crop planning strategies over the coming few decades. Although the mid-term atmospheric CO₂ levels (~550-to-650 ppm under SSP5-8.5) are lower than those projected for the long-term period, they still represent a significant increase compared to the near-term period (2021–2030). Our results show that even at these intermediate concentrations, ET o CO 2 is consistently lower than ET o original across most regions of India (Figures 4, 5), suggesting that stomatal closure effects begin to noticeably influence evapotranspiration well before the end of the century (Figure 9). This mid-term reduction in ET_o has critical implications for regional irrigation scheduling and water resource allocation, especially in semi-arid and arid zones where small changes in evaporative demand can significantly alter water availability (Konapala et al., 2020). Furthermore, by capturing the gradual onset of CO₂-driven feedbacks on evapotranspiration, our study emphasizes the importance of accounting for dynamic CO₂ trajectories even in near- and mid-term projections, an aspect often overlooked in traditional ET_o estimation frameworks.

The outcomes of this study suggest that atmospheric CO₂ can greatly impact India’s annual water budget. Rainfall in India is expected to rise by 6 to 14% under various climate scenarios (Chaturvedi et al., 2012; Kumar et al., 2013) by the end of the century. The combination of reduced ET_o due to incorporation of CO₂ and seasonal variations may exacerbate the difference between water demand and supply both spatially and temporally, potentially leading to water scarcity during peak agricultural demand and flooding during the monsoon season when water demand is minimal. Evapotranspiration, accounting for approximately 40% of India’s water budget (Narasimhan, 2008), is anticipated to be over-estimated by approximately 29% over the long term (Figure 9) under SSP5-8.5 if the effects of CO₂ are disregarded. However, under SSP1-2.6, the contribution of CO₂ is not significant because of limited increase in atmospheric CO₂ concentration and temperature. This underscores the critical necessity to incorporate CO₂ in models that forecast water demand and supply, especially when considering business-as-usual scenario such as SSP5-8.5.

The differences between ET o original and ET o CO 2 averaged across five GCMs using the multimodal ensemble method, shows significant spatial variations across India for the four seasons (monsoon, post-monsoon, winter, and pre-monsoon; Figure 7). These variations in the ET_o caused by atmospheric CO₂, often ignored in climate change impact assessments on agricultural water demand (Sreeshna et al., 2024) and drought (Aadhar and Mishra, 2020; George and Athira, 2025; Varghese and Mitra, 2025), could play a critical role in future water resources planning in India (Varghese and Mitra, 2024). Generally, a reduction in ET_o corresponds to a decrease in agricultural water demand, aligning with field observations (Ainsworth and Long, 2021). However, this effect is often overlooked or inadequately represented in climate change impact assessments on water resources (Döll et al., 2015; Athira et al., 2023) resulting in poor predictions of water availability and demand in the agricultural sector. Our projections indicate that the impact of CO₂ on ET will remain moderate from the near-term to mid-term and long-term time periods under SSP1-2.6, but will intensify significantly under SSP5-8.5, potentially leading to severe consequences for various sectors intricately linked to climate change. For example, the major grain producing states in India (Uttar Pradesh, Madhya Pradesh; Ministry of Finance, 2023) are expected to experience a decrease of 1.1-to-2.8 mm day⁻¹ in ET o CO 2 (Figure 5b) by the end of the century, which appears to be beguiling in terms of reduced agricultural water demand, but could pose serious challenges for the management of water resources and extreme hydrological events, if CO₂ effects are not considered. In this context itself, we must also consider seasonal variations (Kingra et al., 2024; Figure 7b) as water consumption may differ throughout the seasons. The seasonal variation (see Supplementary Table S6; Figure 7b) has major consequences for extreme events such as flooding (Döll et al., 2015) and heatwaves (Ford and Schoof, 2017). While runoff is estimated to be more responsive to variations in precipitation than ET_o (Bharat and Mishra, 2021), the role of ET_o is likely to become more important in a future with higher levels of CO₂ (Davie et al., 2013; Meng et al., 2016). Flood-prone regions in India (Chakraborty and Joshi, 2016) may likely experience a decrease in ET_o of up to 2.1 mm day⁻¹ during the monsoon season (Figure 7b), potentially worsening the flood conditions. The increasing severity of heatwaves attributed to climate change in Rajasthan, Bihar, West Bengal, specific areas of Kerala, and northeastern India may be underestimated, as prior work on impact of climate change on heatwaves (Dubey and Kumar, 2023; Ravindra et al., 2024) did not account for the influence of rising CO₂ levels on ET_o and, subsequently, on heatwaves. Additionally, neglecting CO₂ in seasonal ET_o estimates can influence prediction of flash droughts which are closely linked to evapotranspiration rates (Mahto and Mishra, 2020; Wang et al., 2016; Pendergrass et al., 2020). Accurately representing the role of CO₂ in estimating ET is crucial for hydroclimatic forecasting, as it helps explain contradictory phenomena like the observed greening of the earth despite continental drying (Milly and Dunne, 2016; Chen et al., 2023) and inconsistent runoff estimations (Kooperman et al., 2018; Zhou et al., 2023; Lesk et al., 2024). The implementation and scaling of large-scale land-based climate solutions (Jaiswal et al., 2025) that rely on plants will also be influenced by the accurate representation of evapotranspiration, particularly in the context of water demand and supply.

In addition to CO₂, other factors such as vegetation, temperature, rainfall (Lovelli et al., 2010), and vapor pressure deficit (VPD) (Ort and Long, 2014) will also impact ET. Assessing the complex interactions among these variables including teleconnections between climate variables is challenging (He et al., 2022; Sidhan and Singh, 2025). Both ET_o (Soni and Syed, 2021) and climate teleconnections (Sharma et al., 2020; Sahu et al., 2025) can influence the water budgets of India’s 20 major river basins, which collectively provide an average of 1,914 billion cubic meters of replenishable water resources (Bassi et al., 2020). India’s multiple river basins are not only hydrologically fragmented but also affected by the non-uniform distribution of rainfall, both of which pose major challenges to achieving nationwide water security. In response, the national river-linking project was proposed to redistribute water from flood-prone regions to water-scarce areas and to manage rainfall variability according to regional demand. However, the original design of this initiative did not account for the impacts of climate change. This omission is particularly critical, as rising atmospheric CO₂ levels can significantly influence ET_o, especially in large river basins where hydrological responses are highly sensitive to land use/land cover changes (Das et al., 2018) and climatic variability (Sarker, 2022). To ensure long-term sustainability and effectiveness, the river-linking project must incorporate climate change considerations — particularly the effects of elevated CO₂ on ET_o — in evaluating impacts on river networks (Abed-Elmdoust et al., 2016), the role of critical hydrological monitoring nodes (Singhal et al., 2024), the maintenance of river network integrity (Sarker et al., 2019), and the strategic placement of dams (Gao et al., 2022).

Our approach is based on empirical data that does not make distinction between C3 and C4 crops (Li et al., 2019). It is well recognized that these plant types respond differently to elevated CO₂ (Leakey et al., 2019) but semi-empirical approaches such as FAO-PM are not suitable to incorporate such details. Making ET_o predictions after accounting for the photosynthetic pathway (C3 or C4 types) would require using biophysical approaches (Lochocki et al., 2022) that model behavior of stomatal conductance to CO₂ concentration (Ball et al., 1987) while accounting for leaf biochemical characteristics. Model coupling tools (Surendran and Jaiswal, 2023) can also be a simpler way to enhance existing models neglecting CO₂ concentration to make reliable predictions under rising CO₂ concentrations. The sensitivity of evapotranspiration water losses may also vary (Lockwood, 1999) across different landcover types (such as grasslands, slow growing tall canopies) and these factors are also not included in our modified approach.

Our analysis aimed to provide a broad, country-level perspective on the impact of CO₂ on reference evapotranspiration using ISIMIP data at a spatial resolution of 0.5° × 0.5°. For studies focused on finer spatial scales or specific regions within India (Mishra et al., 2020), incorporating spatial downscaling or using regional climate models (RCMs) would enhance the resolution and enable more localized insights. Such approaches can be particularly valuable for translating large-scale climate projections into actionable information at the state or district level.