Impacts of climate change and vegetation response on future aridity in a Mediterranean catchment

on projected water balance components and the resulting impact on aridity in a medium-sized catchment of Central Italy. We validate and couple a hydrological model with climate projections from five regional climate models and perform simulations considering the vegetation responses or not. Results show that their inclusion significantly affects potential evapotranspiration. The other water balance components, namely actual evapotranspiration, water yield, percolation, and irrigation, are also influenced but with less significant changes. Considering or not the CO 2 suppression effect on stomatal conductance, coupled with the uncertainty related to precipitation, highly affects the estimation of future aridity as the future climate classification ranges from “ humid ” to “ semi-arid ” depending on the simulation and climate model, even if model outputs need to be evaluated cautiously with CO 2 concentration higher than 660 ppm.


Introduction
The Mediterranean area is considered a hotspot for climate change since, compared to other regions, temperatures will rise 20% faster and precipitation will decrease 4% faster per degree warming than the global average (Lionello and Scarascia, 2018).Moreover, it will face increased extreme heat, heavy precipitation, and hydrological and agricultural droughts (Arias et al., 2021).Focusing on future precipitation in the Mediterranean region, many studies highlighted a clear North-South gradient, with the Southern areas facing the most severe impacts of climate change.Despite the great uncertainty, the zero-change line in precipitation is usually estimated to cross Northern and Central Italy (e. g., Coppola et al., 2021;Mariotti et al., 2015;Spano et al., 2020).
Understanding future precipitation trends and their spatial patterns over Italy is crucial given their large socio-economic and environmental impact (WHO, 2018).Indeed, according to the Organisation for Economic Co-operation and Development (OECD) classification, Italy is currently considered a medium-high water-stressed country since more than 30% of renewable water resources are used, with a large share dedicated to agriculture (PNACC, 2018).Furthermore, water availability and demand are very inhomogeneously distributed throughout the country.Past studies on precipitation variability showed that there is a negative trend in total annual precipitation for the whole Italian territory, which is most pronounced in the winter season (Caporali et al., 2021).Interestingly, these past trends oppose future precipitation projections that simulate a reduction in summer precipitation and a slight increase in winter (Spano et al., 2020).More specifically, over Central Italy, no considerable past trends were obtained for yearly precipitation, while positive trends in autumn and winter were identified.At the same time, meteorological drought analyses revealed an increasing trend in both severity and frequency for two other regions in Central Italy, Abruzzo and Umbria (Di Lena et al., 2014;Vergni and Todisco, 2011).
Past and future warming over land has a strong impact on the atmospheric evaporative demand (here referred to as potential evapotranspiration, PET), leading to a reduction in soil moisture and increased agricultural and ecological droughts (Douville et al., 2021;Seneviratne et al., 2021).Vicente-Serrano et al. (2022a) stressed the importance of PET changes in the observed increase in agricultural and ecological droughts.Temperature and humidity are strongly linked and the change in relative humidity is tightly related to the soil-moisture availability (Drobinski et al., 2020).The surface-drying effect due to increased temperatures will be reduced in the Northern Mediterranean due to the availability of sufficient soil moisture (Tramblay et al., 2020).Mariotti et al. (2015) inferred that in Europe future evapotranspiration changes over land were mainly linked to changes in projected precipitation.Surface insolation is projected to increase over the whole Mediterranean, mainly because of reduced cloudiness (Coppola et al., 2021).On the other hand, soil moisture is projected to decrease in many areas of the Mediterranean, but with high uncertainties; for example, no noticeable changes are projected over large parts of Italy (Mariotti et al., 2015).
Another major source of uncertainty when trying to assess the future aridity conditions is the direct effect of CO 2 concentration increase through changes in plant transpiration and growth (Manzoni et al., 2022;Vicente-Serrano et al., 2022b).CO 2 rising boosts crop growth through the CO 2 fertilization effect and reduces plant transpiration and PET through the CO 2 suppression effect on stomatal conductance (Zhang et al., 2022).Overall, the decrease in transpiration caused by the stomatal conductance reduction is compensated by the increase in transpiration caused by the CO 2 fertilization effect (Manzoni et al., 2022), especially in dry and semi-arid climates (Fatichi et al., 2016).Increasing CO 2 , therefore, has an indirect effect on runoff which is the main factor explaining the discrepancy between a projected increased future runoff as predicted by climate models and a drying trend that is projected from future drought and aridity estimations (Yang et al., 2019).Globally, the greening effect of CO 2 rising and climate change was demonstrated with evidence from the last ice age and the historical era (Scheff, 2018;Zhu et al., 2016).Assessing future aridity conditions using temperature-based indices without accounting for the CO 2 fertilization and stomatal suppression effects may result in an incomplete assessment (Scheff, 2018;Swann et al., 2016).Nevertheless, the current increase in CO 2 concentration is occurring at an unprecedented rate at which ecosystems might not be able to take advantage; also, nutrient availability might limit the positive effect of CO 2 fertilization (Scheff, 2018).Climate models already take the CO 2 fertilization mechanism into account, but this effect might be overestimated, and caution is suggested when directly using the outputs of climate models in the estimation of future drought (Vicente-Serrano et al., 2022a).
While the uncertainty related to projected precipitation is vastly explored in literature, the one related to PET, linked mainly to the CO 2 suppression effect on stomatal conductance, is rarely estimated.The objective of this study is to demonstrate the importance of considering the plant physiological responses to CO 2 , mainly stomatal conductance reduction, when calculating projected PET, quantifying the difference in simulated water balance components and aridity when considering or not the vegetation responses.The impacts of climate change on water resources such as the water balance components and river flows are frequently assessed by coupling climate and hydrological models (Tramblay et al., 2020).Among the various hydrological models, the Soil and Water Assessment Tool (SWAT) has been frequently used in the Mediterranean region, including for climate change analyses in Italy (Aloui et al., 2023).Here, we use as meteorological data five bias-corrected EURO-CORDEX Regional Climate Models (RCMs) and evaluate their projected trends for precipitation and temperature in the Ombrone catchment in Central Italy.These are then used to force the SWAT+ model (Arnold et al., 2018;Bieger et al., 2017) to simulate future water balance.We evaluate the model performance in predicting monthly streamflow using a multi-site calibration and validation process.We then simulate PET with constant and decreasing stomatal conductance and assess the effect of the plant physiological responses to CO 2 on aridity and other water balance components.In this study, we also focus on the SWAT+ approach to estimate PET upon exceedance of the 660 ppm threshold, which is considered the maximum value at which the equations used by the model are valid.Understanding the mechanisms behind and improving the quantification of the CO 2 suppression effect on stomatal conductance and CO 2 fertilization will allow the development of more robust aridity projections which, in turn, are required to optimally plan potential adaptation measures.

Study area: the Ombrone catchment
Coastal, small-to-medium-sized, temporary rivers prevail in the Mediterranean region, accounting for more than half of the total area (Ducrocq et al., 2016).The Ombrone catchment, located in Central and Southern Tuscany, is a typical example of a Mediterranean catchment, with an area of 3552 km 2 , a maximum elevation of 1738 m.a.s.l., and the river outlet in the Tyrrhenian Sea (Fig. 1).The Ombrone river is the second longest river of Tuscany with a length of 161 km and has several tributaries, including the Arbia, Merse, Farma and Orcia (Diodato et al., 2023).The Ombrone catchment is almost entirely included in the provinces of Siena and Grosseto.This part of Tuscany is considered substantially water-stressed due to the high concomitant water demand for agriculture and tourism, especially in the coastal areas (Villani et al., 2022).Southern Tuscany is also the area of the region that receives less precipitation and is experiencing the most pronounced increases in dry spell occurrence (Bartolini et al., 2022).Furthermore, the analysis conducted by Diodato and Bellocchi (2008) classified this area as prone to agricultural drought.According to the Köppen classification, the prevalent climate can be described as a hot-summer Mediterranean (Csa) climate (Beck et al., 2018), characterized by hot, dry summers and cool, humid winters, with an annual precipitation of 600-1100 mm according to the data used in this study.The climate of the internal areas of the Ombrone catchment is more continental compared to the coastal areas, with slightly higher precipitation and shorter summer (Diodato et al., 2023).Future projections indicate that this might shift towards hot and cold semi-arid (BSh and BSK) climates (Beck et al., 2018).The Ombrone catchment is characterized by hilly and mountainous areas with slopes of over 20% and torrential streams (Diodato et al., 2023).According to the Corine Land Cover of 2018 used in this study, the main land covers of the catchment are forest (39%) and herbaceous annual crops (46.87%), and the other types of vegetation are permanent pastures (4.28%), vineyards (3.84%), and olive groves (1.89%).The remaining parts of the catchment are covered with artificial surfaces, shrubland and bare land.
In the last sixty years, Tuscany has experienced considerable economic growth that led to a reduction in the number of farms and the abandonment of cultivated land, especially in the less productive areas (Napoli et al., 2017).The upland areas of the Ombrone catchment in the Siena province are mainly cultivated with non-irrigated crops such as cereals, olives, and grapevine.Vineyards and olive groves are widespread in the whole of Tuscany and are prevalent in rugged and sloped areas prone to erosion (Napoli et al., 2014;Napoli and Orlandini, 2015).In the coastal areas of the Grosseto province, irrigation is more widely used and orchards and horticultural crops are common.Water is thereby pumped mainly from the coastal aquifer of the Grosseto plain, which suffers from overexploitation, seawater intrusion, and pollution (Aldinucci et al., 2012;Zucaro and Tudini, 2008).

The SWAT+ model
The SWAT+ model is a restructured version of SWAT that offers better spatial discretization of the catchment, improved code maintenance and more flexibility in representing management practices compared to the previous versions (Arnold et al., 2018;Bieger et al., 2017).This dynamical model uses a daily time step, and represents the catchment with a semi-distributed approach, by dividing it into subbasins, landscape units, and Hydrological Response Units (HRUs).The HRUs each have homogeneous characteristics of soil, slope, and land use.A key improvement in SWAT+ upon SWAT is the inclusion of decision tables, which allow for an improved representation and modelling of complex rules related to water and agricultural management (Arnold et al., 2018).
The SWAT/SWAT+ modelling suite can be conveniently used to evaluate the impacts of the plant physiological responses to CO 2 , namely the CO 2 suppression effect on stomatal conductance and the CO 2 fertilization, on streamflow and water balance components (Wang et al., 2017;Wu et al., 2012).The modification of the Penman-Monteith approach to simulate the suppression effect on stomatal conductance included in SWAT is commonly applied (e.g., Lemaitre-Basset et al., 2022).In the traditional Penman-Monteith equation, stomatal resistance, which is the inverse of stomatal conductance, is assumed to remain constant and is therefore unrealistic (Lemaitre-Basset et al., 2022).Only when the Penman-Monteith method is used, SWAT+ has been adjusted to account for the stomatal suppression effect (Neitsch et al., 2011;Nkwasa et al., 2023).More specifically, in SWAT+ the stomatal resistance (r c ) is allowed to change with Eq. 1 (Easterling et al., 1992) based on an experiment that reached a CO 2 concentration of 660 ppm (Morison, 1987), which is much lower than the projected increase by the end of the century under Representative Concentration Pathway (RCP) 8.5.Moreover, the CO 2 fertilization effect is accounted for in SWAT+ by simulating increased radiation use efficiency (RUE) with Eq. 2, which affects daily biomass accumulation (Neitsch et al., 2011). (1) where r l is the minimum effective stomatal resistance of a single leaf (s m − 1 ), LAI is the Leaf Area Index of the canopy, r 1 and r 2 are the shape coefficients calculated by the model for each crop.
For the SWAT+ model setup, we used the EU-DEM (version 1.0) at 25 m resolution and the 2018 Corine Land Cover and Land Use map from the Copernicus Land Monitoring Service.As soil map, we used the Pedological Database of the Tuscany Region.In this database, soil texture and organic matter content are available as average in the whole soil profile, while hydraulic conductivity is reported for two layers.Other information, such as soil depth and salinity, are reported as categorical variables.To estimate soil properties, pedotransfer functions are typically used (Abbaspour et al., 2019).We estimated available water capacity with the widely used pedotransfer function of Saxton and Rawls (2006) as in Napoli et al. (2017), bulk density with the equation proposed by Manrique and Jones (1991) that performs well in Italy (Pellegrini et al., 2007), the soil erodibility factor with the method of Williams (1995), and the soil albedo with the equation introduced in Sugathan et al. (2014).Climate data were obtained from the Regional Hydrological Service (SIR) of the Tuscany Region.In SWAT+, weather stations are created based on climate input data and assigned to the HRUs.More details on climate input data are available in Table S1.
Since agricultural land occupies a large share of the Ombrone catchment, we defined a simplified representation of herbaceous cropland using the four main crops of the catchment to improve the representativeness of the model.For this, we used the land covers provided by L. Villani et al. the Tuscany region through the ARTEA agency, which are published yearly and contain detailed field-specific characteristics of crop type, crop variety, and crop management (ARTEA, 2018).Hence, in the model, we split the herbaceous cropland use into four classes: durum wheat as the rainfed winter crop (30%), sunflower as the rainfed spring crop (15%), maize as the irrigated spring crop (15%), and alfalfa as the forage crop (40%).We checked and slightly modified the default decision tables already available in the model for sowing and harvesting crops to match the typical sowing and harvesting dates.We included mouldboard and harrow tillage in the management schedule as well as representative fertilization schemes.We prepared a decision table for automatic sprinkler irrigation of 20 mm per event triggered by a water stress threshold of 0.241.By using these values, we obtained volumes of irrigation of the same order of magnitude as compared to the most updated data retrieved from the National Institute of Statistics (ISTAT, 2010).In SWAT+, the water stress is calculated comparing actual and potential plant transpiration (Neitsch et al., 2011).To confirm the consistency of the crop-management schedule, we compared it with reported crop schedules in the Tuscany region (Dalla Marta et al., 2010;Orlando et al., 2015;Giannini and Bagnoni, 2000;Tuscany region, 2010).The impacts of climate change on crop yield in the Ombrone catchment are evaluated and discussed in Villani et al. (2024).

SWAT+ setup and multi-site calibration and validation
We set up the SWAT+ model (revision 60.5.4) in QGIS for a period from 2010 until 2021.Then, we used the SWAT+ Toolbox to perform automatic sensitivity analysis, calibration, and validation for monthly streamflow.The parameters considered for the automatic sensitivity analysis were selected from literature, but we also included other parameters related to the soil map (Table S2).We considered two years of warm up for calibrations, validations and simulations.
Monthly streamflow data for three flow gauging stations were retrieved from SIR (Fig. 1).After the automatic sensitivity analysis performed with the SWAT+ Toolbox, we started calibrating the selected parameters in the upstream flow station of Buonconvento, using five years of monthly flow (2017)(2018)(2019)(2020)(2021).Then, for Sasso d'Ombrone we calibrated the parameters using the same time window and validated with monthly flows from 2012 to 2016, maintaining the calibrated parameters for the subbasins of Buonconvento.We finally repeated the same procedure for Istia, where we considered three years of data for both calibration (2019-2021) and validation (2013-2015).The periods for calibration and validation were selected according to the data availability.Hence, the calibration was performed for all three gauging stations, while the validation was only in two of them since for the Buonconvento gauging station we could rely only on five years of streamflow.Nevertheless, we don't consider this a major problem since Buonconvento is the most upstream gauging station.
To perform calibration and validation, the SWAT+ Toolbox allows the use of the Nash Sutcliffe Efficiency (NSE), the Mean Square Error (MSE) and the Root Mean Square Error (RMSE).As mostly applied with the SWAT/SWAT+ modelling suite, we selected as objective function the NSE during the automatic calibration for monthly streamflow.The percent bias (PBIAS) and the RMSE-observations standard deviation ratio (RSR) were calculated as additional statistics to evaluate the performance of the model.We used the equations and evaluated the model performances following the criteria of Moriasi et al. (2015) for NSE and PBIAS, while those of Moriasi et al. (2007) for RSR (Table S3).

Climate change scenarios
To estimate future climate changes, we used five EURO-CORDEX climate models (Jacob et al., 2014) (Table 1) and RCPs 4.5 and 8.5 (Moss et al., 2010).We considered two 30 year-long periods, 2041-2070 and 2071-2100, to evaluate medium-and long-term climate change impacts, comparing the projected values with those of the historical simulations (from 1976 to 2005) of each climate model.The criteria for selecting the climate models were (1) the availability of both RCPs, (2) the availability at daily frequency of precipitation and the climate variables needed to calculate PET with the Penman-Monteith method (maximum and minimum temperatures, relative humidity, solar radiation, and wind speed), (3) the use of a complete (Gregorian) calendar, (4) a horizontal RCM resolution of 0.11 • over the EURO-CORDEX domain.A bias correction of the climate-projection data is necessary as systematic biases are present in the meteorological data (Maraun and Widmann, 2018).Among the different methodologies that exist for bias correction of climate models, we adopted distribution mapping which is commonly used for climate and hydrological studies (Teutschbein and Seibert, 2012;Themeßl et al., 2011).
We used the CMHyd software to bias correct temperature and precipitation (Rathjens et al., 2016).The software reprojects the data and applies the selected bias-correction method using the station data provided by the user.The CMHyd outputs can be directly used in the SWAT+ model without further preprocessing.Since bias correction is performed by comparing the historical simulation of the climate models with observed values before 2005, we used a lower number of stations compared to the calibration and validation period.The stations used for bias correction were included in the Ombrone catchment or very close and had more than 10 years of data as indicated in Fung (2018).More details about the climate data used for bias correction can be found in Table S1.To process the other climate variables, solar radiation, relative humidity and wind speed, we used the Climate Data Operators (CDO) software, version 2.0.5 (Schulzweida et al., 2021).These bias-corrected climate data were then used to run historical and future simulations.In the result section, we detected the changes comparing future and historical simulations for each climate model and RCP.
In SWAT+, the crop cycle is defined with the number of days required to reach maturity, differently from the older version which used the number of heat units (Nkwasa et al., 2023).To account for the shortening of the crop cycles, at first we retrieved the heat units starting from the days to maturity used during the calibration and validation period.Then, we calculated the new crop cycle length considering the different temperatures in the historical and future periods.
After the calibration and validation of the SWAT+ model, we performed simulations to assess the impacts of climate change considering the different climate models, RCPs, and periods.At first, we evaluated the magnitude and sign of the climate change signal for future precipitation and average temperature in the Ombrone catchment, considering the basin scale outputs of the SWAT+ model.To quantify the impacts of the plant physiological responses to CO 2 on PET and other water balance components, we performed simulations with constant CO 2 at 400 ppm, as in the calibration and validation, and others considering the values as reported by Büchner and Reyer (2022), for RCPs 4.5 and 8.5, for the historical, near and far future periods.We used the 30-year-average CO 2 concentration since it is a fixed input parameter in the SWAT+ model.The average CO 2 concentration for RCP 4.5 (522 and 589 ppm for the near and far futures, respectively) and the near future for RCP 8.5 (611 ppm) fall between the limit of the Morison experiment (330-660 ppm), while the average CO 2 concentration for the far future of RCP 8.5 is much higher (939 ppm).
We evaluated the impacts of climate change and vegetation

Table 1
The five climate models used in the study with the General Circulation Model (GCM) and the Regional Climate Model.

GCM institute GCM model RCM institute RCM model
L. Villani et al. responses to CO 2 on future PET and estimated future dryness conditions with the Aridity Index (AI).AI is calculated as the ratio between precipitation and PET (Middleton and Thomas, 1997) and is useful to classify the climate following the UNEP climate classification.Despite the intrinsic simplifications of this method, AI is used for climate classification (Massari et al., 2022) and as a benchmark for aridity conditions at the regional scale (He et al., 2022).Hence, AI estimation is useful to evaluate the uncertainty related to RCPs and the inclusion of the stomatal conductance reduction caused by CO 2 concentration rising.Finally, we performed a similar analysis for the water balance components and other variables simulated with the SWAT+ model, namely actual evapotranspiration (ET), water yield (a direct output of the model defined as the sum of surface runoff, lateral soil flow, and tile flow), percolation, and irrigation.In the main analysis reported in the results section we refer to the water balance components considering average yearly values, but in the supplementary materials we reported other simulations outputs considering the 5th and 95th percentiles, summer and winter season values and additional variables such as streamflow.
We applied the Wilcoxon rank-sum non-parametric test to detect significant differences between historical and future periods and between simulations including and excluding the vegetation responses to CO 2 .We considered the yearly values of the variables considered in the analysis.The two samples used to test the significant differences were always selected from the same climate model.

Multi-site calibration of SWAT+
Out of the 15 parameters pre-selected for the sensitivity analysis, we considered only the most sensitive, specifically cn2, esco, epco, bd, and revap_co.More details about the calibrated parameters are available in the supplementary materials (Tables S2, S3, S4).Despite overestimation of peak flows, we obtained more than satisfactory performances according to the criteria considered in this study (Table 2, Fig. 2).More in detail, the model achieved very good performances in Buonconvento during the calibration period for the three statistics considered, while in Sasso d'Ombrone we obtained very good performances considering NSE and RSR and good for PBIAS.For the Istia gauging station, we obtained very good performances for NSE and RSR during validation and good for PBIAS, while during the calibration period satisfactory performances for NSE and PBIAS and good for RSR.The reduced performance for Istia during calibration can be mainly attributed to a discharge peak in February 2019 (see Fig. 2) which might be an error in the observed streamflow data since it is not present in the other sites.

Projected temperature and precipitation changes over the Ombrone catchment
A positive, significant climate change signal over the Ombrone catchment was obtained for annual average temperatures, for all future periods and RCPs considered (Fig. 3, Tables 3, S5).While the increase in temperature under RCP 4.5 was of similar magnitudes in the near and far futures, under RCP 8.5 it continued to rise in the far future.In the far future, the highest increases of 18% and 32% were found for NorESM1-M -REMO2015, under RCPs 4.5 and 8.5 respectively.The ensemblemean temperature increase at the end of the century was 2.1 • C and • C for RCPs 4.5 and 8.5, respectively (Table S5).Temperature increases were largest for the summer season and daily minimum temperature.For average, maximum, and minimum temperatures, the increases in summer were higher than those in winter, particularly for RCP 8.5.The ensemble-mean increases were 1.9 • C and 3.5 • C for winter average temperature, for RCPs 4.5 and 8.5, respectively, while for summer these were 2.3 • C and 4.8 • C (Table S5).
For annual precipitation, climate change projections were much more uncertain (Fig. 4, Tables 3, S5).The RCMs disagreed on the sign of change, with four RCMs predicting negligibly small changes or increases and one (NorESM1-M -REMO2015) a decrease (Table 3, Fig. S1).The ensemble-mean average annual precipitation increased at the end of the century by 70 and 32 mm for RCPs 4.5 and 8.5, respectively (Table S5).This difference is mainly caused by the significant decrease in precipitation by NorESM1-M -REMO2015 for RCP 8.5 (-21%) (Table 3).In RCP 4.5 the slight increase occurred by the end of the century, in contrast to RCP 8.5 in which the increase was found in the near future.The increases in precipitation were found mainly in winter, while the models indicated reduced increases or even decreases in spring and summer.More specifically, the ensemble-mean increases in winter average precipitation were 34 and 22 mm for RCPs 4.5 and 8.5 respectively, while the ensemble-mean changes for summer were 5 and − 7 mm for RCPs 4.5 and 8.5 respectively (Table S5).In the far future, 3 out of 5 models projected significant increases in precipitation under RCP 4.5.Under RCP 8.5, two predicted significant increases while NorESM1-M -REMO2015 significant decreases (Table 3).Under RCP 4.5 in the near future, CNRM-CM5-RACMO22E behaved differently as compared to the other climate models and was the only one showing significant increases of 27% in precipitation (Table 3).

Impacts of climate change and vegetation responses to CO 2 on future potential evapotranspiration and future aridity
Projected PET drastically changed, both in terms of magnitude and sign, considering RCPs 4.5 and 8.5 and due to the inclusion of the stomatal conductance reduction effect (Fig. 5a, Table 3).Consequently, the other components of the water balance such as ET, water yield, percolation, and irrigation were similarly influenced (Fig. 6, Table 4).As expected, the changes in PET in the far future were particularly high for RCP 8.5.When not including the suppression effect, PET increased in line with temperature.In that case, under RCP 4.5, the average PET increased up to 110 mm (Table S6) mainly in the near future, while in the far future the increases remained of a similar magnitude.This reflects the increase in temperature that occurred early in the near future and slowed down by the end of the century (Fig. 3).Under RCP 8.5, on the contrary, PET continued to increase until the end of the century, reaching an average increase of 225 mm (Table S7).The CO 2 suppression effect balanced the temperature-induced change in PET, quenching its increase under RCP 4.5 to 16 mm (Table S6).The CO 2 concentration used in the long-term future period under RCP 8.5, on the other hand, was 939 ppm, far above the upper limit of 660 ppm of the Morison experiment.This explains the unrealistic drop in PET in the period 2071-2100, with an ensemble-mean decrease of − 211 mm (Table S7).
The uncertainty in future precipitation as predicted by the five climate models considered (Fig. 4, S1, Table 3) and the one in future PET as predicted by the SWAT+ model forced with the climate projections (Fig. 5a) escalated when considering AI.While AI was simulated to be around 0.65 in the historical simulations, that is the threshold that divides the "dry sub-humid" and "humid" climates, AI drastically disperses depending on the RCPs, climate models, and whether or not the CO suppression effect on stomatal conductance is included (Fig. 5b).By the end of the century, under RCP 4.5, the deviation between the different models will be higher as compared to the historical simulations.When considering stomatal conductance reduction, the AI values ranged between 0.83 and 0.52, while they were slightly lower with constant stomatal conductance, between 0.78 and 0.48.Following the more uncertain predictions under RCP 8.5, the projected AI ranged between 0.72 and 0.36 when considering constant stomatal conductance and between 1.05 and 0.54 when analysing the simulations with stomatal conductance reduction.Therefore, according to the UNEP classification (Middleton and Thomas, 1997), the current humid/dry sub-humid climate is predicted to shift towards a much more humid climate for RCP 8.5 when considering the vegetation responses to CO 2 , while it shifts towards a much more arid climate when considering the same RCP but calculated with constant stomatal conductance.With the driest climate projections of the NorESM1-M -REMO2015 climate model, under RCP 8.5 and constant stomatal conductance, AI was predicted to be below 0.5, the threshold that divides "dry sub-humid" and "semi-arid" climates.

Impacts of climate change and vegetation responses to CO 2 on future water balance components
Since Eq. 1 used in the SWAT+ model was tested only until 660 ppm, we opted to analyze future impacts of climate change and vegetation responses to CO 2 on water balance components only in the near future (2041-2070), with CO 2 concentrations lower than the threshold.
Changes in water yield and percolation were strictly linked to precipitation.In the near future, the climate models which predicted precipitation increases were CNRM-CM5-RACMO22E and MPI-ESM-LR-RCA4 under RCP 4.5 and all except NorESM1-M-REMO2015 under RCP 8.5 (Table 3).If the precipitation increases were mostly not significant, for water yield and percolation the Wilcoxon test always resulted in p-values lower than 5% (Table 4).The percentage changes were much higher for water yield and percolation as compared to precipitation increases, reaching up to 105% and 73% increases for CNRM-CM5-RACMO22E under RCP 8.5 and with stomatal conductance suppression for water yield and percolation, respectively (Table 4).
The variables related to temperature -PET and ETwere reduced when considering the plant physiological responses to CO 2 , while those related to precipitationwater yield and percolationshowed an increase when the vegetation responses were included.If the changes caused by the inclusion/exclusion of the stomatal conductance suppression were significant for PET (Table 3), the changes found for the other balance components were not significant in most cases.We obtained significant changes only under RCP 8.5 for ET (3 models out of 5) and irrigation (2 models out of 5).Fig. 6 reports the absolute values of the ensemble water balances for the two cases considered in this study and their relative percentage difference.Differences in the historical period were minimal (1-2%) and they increased as the CO 2 concentration was higher.Under RCP 8.5, irrigation changed the most with a difference of 10.1%, followed by percolation (-8%), soil evaporation (7.1%), water yield (-5.5%) and transpiration (3.2%).Canopy evaporation was barely affected by the change in PET driven by different CO 2 concentrations.The magnitude of percentage changes caused by the vegetation responses to CO 2 ranged between − 4.8% and 4.1% for RCP 4.5, while between − 8% and 10.1% for RCP 8.5.

Climate models' uncertainty in the Northern Mediterranean area
Our study further confirmed that the projected increases in temperatures in the study region are remarkable, and especially high for RCP 8.5 and during summer, consistent with previous research (e.g., Spano et al., 2020).On the other hand, results revealed even more uncertainty regarding the future precipitation predicted by climate models for the Northern Mediterranean area.The ensemble mean of the five climate models considered in this study indicated an increase in precipitation, more accentuated during winter.Central and Northern Italy are in a transition zone between the arid North African and the humid Central European climate zones, and the zero-change line predicted in past studies usually crosses this area (e.g., Coppola et al., 2021;Mariotti et al., 2015;Spano et al., 2020).A deep analysis of the literature regarding only the Northern Mediterranean area showed very high uncertainty, and in particular when considering Central and Northern Italy.For example, Mariotti et al. (2015) analysed the outputs of CMIP5 experiments, using an ensemble of climate models for RCP 4.5, and found only a minor decrease of 0.2 mm/day during summer by the end of the 21st century, and no decrease when considering winter or annual precipitation over the Northern Mediterranean area.On the contrary, Lionello and Scarascia (2018) considered RCP 8.5 of CMIP5 experiments and found an overall reduction of precipitation except for winter and, to a lower extent, spring months in the Northern Mediterranean area.The analysis of the outputs of the first EURO-CORDEX RCMs ensemble showed that, for Italy, the climate change signal for precipitation was uncertain under RCP 4.5, while it was negative for Central and Southern Italy under RCP 8.5, again with a distinct gradient increasing southward (Jacob et al., 2014).More recently, Coppola et al. (2021) used a much larger ensemble of EURO-CORDEX RCMs and compared the outputs with those of CMIP5 and CMIP6 GCMs, considering RCP 8.5.The three ensembles agreed to show the precipitation zero-change line over the Northern Mediterranean during winter, while it shifted northward in summer, meaning that a decrease in precipitation is projected for summer months over Italy.Evin et al. (2021) also used a large ensemble of EURO-CORDEX RCMs considering RCPs 2.6, 4.5, and 8.5 to estimate future temperatures and precipitation and the relative uncertainty.For Italy, they projected slight increases in winter precipitation and noticeable decreases in summer precipitation, yet with very large uncertainties.It is interesting to underline that, in this study, Italy emerged as the country with the lowest reduction in precipitation among the Mediterranean countries.Focusing only on the Italian territory, CMCC carried out a climate change analysis considering an ensemble of EURO-CORDEX RCMs (Spano et al., 2020).The main results confirmed the North-South gradient for precipitation since a reduction was found for Southern and Central Italy, mainly in summer and, to a lower extent, in spring months, while an increase in winter precipitation over Northern Italy.Moreover, an increasing trend in maximum daily precipitation was found for summer and autumn.Studies summarized in the draft of the National Climate Change Adaptation Plan (PNACC, 2018) reported a reduction in total precipitation, more pronounced in the Southern areas, in the summer season and when considering RCP 8.5.A robust decrease in total precipitation with similar patterns was also found in the analysis carried out by Padulano et al. (2020) which also used an ensemble of EURO-CORDEX RCMs in Italy.The overall picture of increasing precipitation during winter and decreasing precipitation in summer is confirmed also by our analysis when considering RCP 8.5, but the winter increases were higher as compared to the spring and summer decreases, resulting in an overall, yet minor, increase.Instead, for RCP 4.5, minor increases were found also for the summer season.Notably, most precipitation changes in our study were not statistically significant (Table 3).
For water yield and percolation, the sign of change was in line with the predicted change in precipitation, consistent with other studies that coupled the SWAT model with climate models in Italy.For example, Fiseha et al. (2014) used climate variables from one GCM downscaled with three RCMs as input to simulate future precipitation and hydrological water balance components in the upper Tiber basin, in Central Italy, considering two different scenarios of future CO 2 concentrations.Except for one climate model in the lower emission scenario, their results showed a general decrease in precipitation and related variables, mostly during summer.Decreasing trends in precipitation and related water balance components such as water yield, groundwater recharge, and ET were also found for the Candelaro catchment in Southern Italy, considering three RCMs (De Girolamo et al., 2017).Pulighe et al. (2021) applied the SWAT+ model to simulate the future climate in the Sulcis catchment in Sardinia with two RCMs and RCPs 4.5 and 8.5.Their results showed a clear decrease in precipitation only for one climate model with RCP 4.5, while slight increases for the other simulations.PET was predicted to strongly increase and the other water balance components such as surface runoff and percolation decreased because of the increased water loss to the atmosphere.In a small catchment of the Po River delta in Northern Italy, Pesce et al. (2019) used 10 different GCM-RCM combinations and found an average decrease in future precipitation, although with an unclear tendency, especially for RCP 4.5 in the medium-term future (2041-2070).Future water flow was projected to increase in the wet season and decrease during spring and summer.Finally, Glavan et al. (2015) simulated climate change impacts with six different climate models in a small Slovenian catchment, very close to the Italian border.Their results showed that PET increased as well as ET if precipitation also increased.Also, precipitation was projected to increase with few exceptions by the end of the century, and stream flows showed consistent increases but higher in magnitude.This is in line with

Table 3
Precipitation, average temperature and potential evapotranspiration in the historical, near and far future periods for the five climate models considered in this study.For potential evapotranspiration, the values for both cases, constant and reduced stomatal conductance, are reported.Significance levels of the Wilcoxon test between future and historical periods and between the two cases considered are also included in the table.

Climate model
Precipitation (mm) Average temperature ( our study since we also observed that the percentage change in water yield was much higher as compared to the increase in precipitation.It is worth noting that in our study and the ones previously discussed, a subset of the available climate models was used, which might lead to under-representative estimates enhancing the uncertainty (Evin et al., 2021).

Impacts of vegetation responses to CO 2
The roles of the stomatal conductance reduction and the CO 2 fertilization effects caused by the CO 2 concentration rising are still unclear and debated.Experiments using earth system models showed that the CO 2 physiological response of vegetation on ET and long-term runoff had higher impacts compared to radiative or precipitation changes caused by CO 2 rising (Lemordant et al., 2018).Furthermore, plant physiological responses to CO 2 were also found to reduce future drought-stress predictions (Swann et al., 2016).Moreover, Skinner et al. (2018) demonstrated that the CO 2 -driven vegetation changes amplified the frequency and intensity of summer heat waves.GCMs predict higher temperatures compared to RCMs, and this was explained by the fact that the latter generally do not include the vegetation response to CO 2 rising since when including this effect, the temperature predicted by RCMs increased (Schwingshackl et al., 2019).However, the study of Taranu et al. (2022) showed that plant physiology had a limited effect and did not explain the large discrepancies observed between GCMs and RCMs.Finally, Vicente-Serrano et al. (2022a) argued that the physiological effect of vegetation as included in climate models might be overestimated and that their outputs should be used with caution when studying future droughts.
The CO 2 rising effect on PET and therefore on all the indices that use it to infer future drought and aridity conditions are remarkable, yet not completely understood (Scheff, 2018;Vicente-Serrano et al., 2022b).The effect of CO 2 rising on future mechanisms and processes relevant to the estimation of future aridity conditions was also analysed quantitatively in climate and hydrological studies, proving that the impacts are not negligible, and stressing the importance of understanding and quantifying them better to obtain reliable future projections of drought and aridity.For example, Greve et al. (2019) demonstrated that the future estimation of PET was largely influenced by the method used to calculate it and that this uncertainty also affects the validity of indexes such as AI.Zhou et al. (2022) found opposite trends for past conditions in China when calculating PET using modified and traditional Penman-Monteith equations by including or excluding stomatal conductance reduction.They concluded that ignoring this effect results in an important PET overestimation, especially in arid regions.However, ET in arid and semi-arid regions is mainly controlled by soil moisture and is not very sensitive to PET (Dakhlaoui et al., 2020), meaning that this overestimation might be a problem in humid regions.Our analysis confirmed that the increase in PET did not lead to a proportional increase in ET.More in detail, in our study the difference in annual ET between the two cases considered under RCPs 4.5 and 8.5 in the near future was around 5% (Fig. 6).This is consistent with the compensative effect caused by CO 2 fertilization (Manzoni et al., 2022) and in line with the magnitudes of changes of less than 8% reported by Fatichi et al. (2016).Lemaitre-Basset et al. (2022) showed that the stomatal conductance reduction effect has a strong impact on future runoff projections over France, while Boé (2021) reported that the decrease in ET caused by the physiological effect of CO 2 did not result in an increase in river flows and soil moisture due to reduced precipitation in summer over France.Yang et al. (2019) demonstrated that the outputs of climate models are similar to those of hydrological models when accounting for the suppression effect in the calculation of PET.Using SWAT, multiple studies showed that ET was reduced by the plant physiological responses to CO 2 , leading to substantial increases in runoff, recharge and discharge (Ficklin et al., 2009;Kishawi et al., 2022;Lee et al., 2018;Van Liew et al., 2012).Other studies evaluated the impact of the stomatal conductance reduction and CO 2 fertilization by modifying SWAT to include dynamic CO 2 concentration as input, finding similar results in terms of increased streamflow and reduced evaporation (Butcher et al., 2014;Wang et al., 2017;Wu et al., 2012).Notably, all these studies were conducted with the older SWAT model versions and run over the United States.
Our results are in line with previous studies, since the impact of vegetation responses to CO 2 on future PET, and therefore on future water fluxes and aridity, was high when considering RCP 8.5.

Table 4
Water yield, percolation, actual evapotranspiration and irrigation in the historical and near future periods for the five climate models and for both cases, constant and reduced stomatal conductance, considered in this study.Significance levels of the Wilcoxon test between future and historical periods and between the two cases considered are also included in the table.6, the magnitude of change in water yield, comparing simulations considering and ignoring the physiological responses, was approximately 5% under RCP 4.5, while it reached more than 20% under RCP 8.5 (Fig. S2).Butcher et al. (2014) reported increases ranging from 3% to 38% and a median of 11%.The contribution of the plant physiological responses to CO 2 estimated by Wu et al. (2012) amounted to 22% for streamflow by the end of this century, much higher than the effect quantified at 1-4% of the recent decades.Marginal increases in streamflow of approximately 1% were instead reported in the study conducted by Wang et al. (2017).
In SWAT+, both the CO 2 suppression effect on stomatal conductance and the CO 2 fertilization effect are considered.Nevertheless, they are calculated in different steps since the stomatal conductance is reduced when using the Penman-Monteith equation to retrieve PET and the fertilization effect when calculating daily accumulated biomass.This might cause some inconsistencies due to the leaf-and canopy-levels transpiration changes caused by these two effects (Manzoni et al., 2022), but the magnitudes of the reductions in ET and the other variables seem to confirm that the increase in biomass partially compensates the decrease in transpiration caused by reduced stomatal conductance.Furthermore, the approach of using climate inputs from GCMs or RCMs to force a hydrological model, used in our study and the papers previously discussed, has some limitations that need to be considered.With this one-way coupling, the interactions and feedback between climate and vegetation are mostly neglected (Wu et al., 2012).Coupling offline hydrological models which do not account for physiological responses of vegetation with climate models is questionable (Milly and Dunne, 2017), especially when the climate models consider these effects (Boé, 2021).As reported by Schwingshackl et al. (2019), most of the GCMs that they used in their study included the CO 2 vegetation response while none of the RCMs considered it.Indeed, the physiological effect of vegetation induces a larger decrease in precipitation which should be compensated by a decrease in ET (Boé, 2021).The opposite problem might occur if the offline hydrological model simulation of the physiological effect of CO 2 on ET is inconsistent with the strength of the vegetation response as simulated by the climate model (Boé, 2021).Furthermore, Swann et al. (2016) suggested that using outputs of earth system models in hydrological models may lead to overestimation of the future drought stress due to double counting of plant feedback on surface humidity, temperature and net radiation.

Conclusion
In this study, the SWAT+ model was forced with climate data from five EURO-CORDEX climate models to estimate future climate change impacts in a typical Mediterranean catchment, the Ombrone catchment in Central Italy.Future aridity conditions were also estimated considering constant and decreasing stomatal conductance.The model performed well after the multi-site calibration carried out for three gauging stations, considering monthly streamflow.
In contrast to temperature, high uncertainties exist in the future precipitation trends.Only one climate model predicted a clear decrease in future precipitation following RCP 8.5 while the others showed minor increases or constant values.The ensemble mean of winter precipitation increased while summer precipitation remained almost constant or slightly decreased, with an overall increase in annual average precipitation.
The impact of stomatal conductance suppression on future PET was significant and should be taken into account when performing future aridity or drought analyses and in particular when applying the SWAT/ SWAT+ modelling suite.Upon disregarding, high increases in PET were obtained, while minor PET increases or even decreasing values were found upon consideration of the suppression effect.Under RCP 8.5 in the far future period, the differences between disregarding or including the vegetation responses to CO 2 were nearly 50% in PET and ranged from 20 to 30% for the other water balance components (Fig. S2).
The SWAT+ model considers the CO 2 effect on future PET with a modification of the Penman-Monteith equation based on an experiment that was conducted for a range of CO 2 values between 330 and 660 ppm.RCP 8.5 predicts much higher CO 2 concentration values by the end of the century.For RCP 8.5, when considering stomatal conductance suppression, we found a dubious drop in PET of more than 200 mm.Further research is certainly needed, but the outputs of the SWAT+ model when excluding vegetation responses to CO 2 and when considering CO 2 concentrations much higher than 660 ppm are prone to large uncertainties.Nevertheless, the Penman-Monteith equation is recommended when using SWAT+ to assess future climate change impacts to account for the effect of reduced stomatal conductance.
The uncertainty in future precipitation and atmospheric evaporative demand patterns strongly increases when considering measures such as the Aridity Index and, consequently, this notably affects future climate classification.Unravelling the uncertainties related to future precipitation in transition zones, like the Northern Mediterranean area, and the plant physiological responses caused by rising CO 2 concentration on future atmospheric evaporative demand is crucial to better understand climate change impacts and plan more effective adaptation strategies.
interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.The Ombrone catchment with the three gauging stations, the subbasins, and the boundaries of the provinces of Siena and Grosseto.

Fig. 2 .
Fig. 2. Multi-site calibration and validation, with the location of the gauging stations of Buonconvento (a), Sasso d'Ombrone (b), and Istia (c), and the respective hydrographs with simulated and observed monthly average flows.

Fig. 3 .
Fig. 3. Lineplot and boxplots of future temperature ( • C).The lineplot (a) reports the absolute difference between the future years (2041-2100) and the yearly average of the historical period (1976-2005), for RCPs 4.5 and 8.5.The line represents the average of the five climate models, and the band is the confidence interval.The values of average temperature calculated with the same procedure considering annual, winter and summer values, for the near (b, 2041-2070) and far futures (c, 2071-2100), are used in the boxplots.

Fig. 4 .
Fig. 4. Lineplot and boxplots of future precipitation (mm).The lineplot (a) reports the absolute difference between the future years (2041-2100) and the yearly average of the historical period (1976-2005), for RCPs 4.5 and 8.5.The line represents the average of the five climate models, and the band is the confidence interval.The same values considering annual, winter and summer values, for the near (b, 2041-2070) and far futures (c, 2071-2100), are used in the boxplots.

Fig. 5 .
Fig. 5. Potential evapotranspiration (PET) and Aridity Index (AI) with constant or decreased stomatal conductance.(a) Lineplot of projected PET (mm), reported as the absolute difference between the future years (2041-2100) and the yearly average of the historical period (1976-2005), for RCPs 4.5 and 8.5, considering constant and decreasing stomatal conductance.The line represents the average of the five climate models, and the band is the confidence interval.The abrupt change occurring in 2070 is due to the fact that the CO 2 concentration value was averaged and considered constant for the two future periods (2041-2070 and 2071-2100).(b) The AI, calculated as the ratio between precipitation and PET for the historical (1976-2005), medium-term (2041-2070) and long-term (2071-2100) future periods considering RCPs 4.5 and 8.5, and constant and increasing CO 2 concentration values.The dots represent the average values and the uncertainty is reported considering the five climate models.Note the unrealistic change in PET, and consequently in AI, for the far future and RCP 8.5 when considering stomatal conductance suppression.

Fig. 6 .
Fig. 6.Differences in water balance components caused by inclusion of the vegetation responses to CO 2 .The bar plots report the water balance components of the various combinations of simulations considering the historical period, RCPs 4.5 and 8.5 and the different CO 2 concentrations.Only the near future period (2041-2070) is considered, so that both RCP simulations have CO 2 concentration values <660 ppm.The names and definitions of the water balance components correspond to those of the SWAT+ model and are divided in inputs (precipitation, snow and irrigation) and outputs (canopy evaporation, transpiration, soil evaporation, water yield and percolation).The water balance values used are the averages of the five climate models.For each couple of bar plots, in the adjacent tables are reported the absolute values of the water balance components (expressed in mm) for the cases of constant stomatal conductance (400 ppm) and modified CO 2 concentration.Furthermore, in the Δ column of the tables, the percentage change due to the vegetation responses to CO 2 is reported.

Table 2
Model performances for monthly streamflow during calibration and validation for the three gauging stations considered.For Buonconvento only the calibration was carried out.
a Very good; b Good; c Satisfactory.L.Villani et al.
Lemaitre-Basset et al. (2022)t and Ulbrich, 2003;Lee et al., 2018)stomatal conductance samples† †p-value < 0.01 for the Wilcoxon rank-sum test with constant and reduced stomatal conductance samples Particularly,Lemaitre-Basset et al. (2022)applied the two methods that we used to estimate PET, namely the non-modified Penman-Monteith equation and the modified Penman-Monteith as proposed byStockle et al. (1992), finding very similar future trends to those we identified.In our study, differences in PET estimation ranged from more than 200 mm increases when considering constant stomatal conductance to decreases of the same magnitude when considering stomatal conductance reduction under RCP 8.5.The magnitude and sign of changes when considering plant physiological responses caused by CO 2 concentration rising were consistent with other studies that applied SWAT.It is interesting to note that when past studies considered CO 2 concentration beyond 660 ppm they found high percentage decreases in future ET as compared to baseline periods.More specifically,Ficklin et al. (2009),Lee et al. (2018)andKishawi et al. (2022)used 970 ppm, 850 ppm, and 935 ppm, obtaining decreases of 40%, 30% and 32% respectively.It has been already hypothesized that the simulated reduction of ET caused by the plant physiological responses to CO 2 is overestimated by SWAT, but this was attributed to several simplifications in the equations used by the model and not by the invalidity of the equations used over 660 ppm(Butcher et al., 2014;Eckhardt and Ulbrich, 2003;Lee et al., 2018).Considering the findings ofLemaitre-Basset et al. (2022), our results confirm that the method included in the SWAT/SWAT+ modelling suite for CO 2 concentrations higher than the 660 ppm threshold is questionable, and simulation outcomes should be interpreted with caution.On the other hand, PET increases with constant stomatal conductance in our study amounted to 225 mm, which corresponds to more than 20% in relative change.Regarding water yield changes caused by the inclusion of the plant physiological responses to CO 2 , our findings agree with those of previous studies in the positive sign of change.As shown in Fig.
*p-value < 0.05 for the Wilcoxon rank-sum test with historical and future samples **p-value < 0.01 for the Wilcoxon rank-sum test with historical and future samples †