Quantifying the Hydrological Impacts of Irrigation on a Mediterranean Agricultural Context Through Explicit Satellite‐Derived Irrigation Estimates

Water budget over anthropized basins is profoundly altered by human interventions, among which irrigation prevails. This study investigates how such water use affects hydrological flux and state variables by comparing two SURFEX/ISBA (SURFace EXternalisée/Interaction Sol Biosphère Atmosphère) Land Surface Model simulations performed over a portion of the Ebro basin (Spain) that includes a heavily irrigated area. The first simulation is forced by atmospheric input only, while the second one is fed with atmospheric forcing plus remote sensing‐derived irrigation amounts added to the liquid precipitation. Four main hydrological output variables are compared under the two configurations: surface soil moisture (SM), total evaporative flux (i.e., the sum of components from bare soil and vegetation) (E), drainage, and runoff. Results indicate that SM and E are the variables mostly impacted during the main watering season. Considering mean daily rates per each month, a maximum increase of +30% and of +220% is found over the irrigation districts in July for SM and E, respectively. The reliability of the results is supported by comparisons with satellite evapotranspiration data from MODIS (MODerate‐resolution Imaging Spectroradiometer). The inclusion of irrigation reduces the RMSE between monthly anomalies of modeled E and of MODIS reference data by 41% during the main watering season. Further validation through wavelet coherence analysis is also proposed. This study sheds light on the potential of explicit satellite‐derived irrigation water use data that, when integrated into modeling platforms, allow us to track and quantify anthropogenic impacts on water resources.


Introduction
Irrigation is the predominant human intervention on the hydrological cycle.Agricultural water represents the primary source of anthropogenic water use, largely prevailing on domestic and industrial components (Abbot et al., 2019;Campbell et al., 2017;Foley et al., 2011).The supply of irrigation systems, if not rationally managed, is known to put pressure on water resources.In fact, several studies highlight irrigation-induced threats on freshwater availability (Ferguson & Maxwell, 2012;Hoang et al., 2019;Khazaei et al., 2019;Taye et al., 2021;Yin et al., 2021) and quality (Foster et al., 2018;Merchán et al., 2013;Park et al., 2018), as well as implications of irrigation practices on climate (de Vrese et al., 2016;Mishra et al., 2020;Rosa et al., 2020;Thiery et al., 2020).Nevertheless, the general absence of explicit spatio-temporal data about irrigation water use makes the quantitative analysis of irrigation impacts on hydrological flux and state variables a challenging task.
So far, studies addressing this research topic through different models rely on irrigation data from their own schemes or from statistical surveys (Dorigo et al., 2021).Lo and Famiglietti (2013) highlighted an anthropogenic loop in the water cycle over California's Central Valley triggered by irrigation-induced evapotranspiration increases and resulting in rising summer precipitation rates and streamflow.They used the Community Atmosphere Model (Gent et al., 2009) in combination with the Community Land Model (CLM) (Oleson et al., 2008).A similar finding was obtained by Wei et al. (2013), by exploiting reanalysis data from MERRA (Modern-Era Retrospective Analysis for Research and Applications) (Rienecker et al., 2011).They showed that regions where irrigation is largely practised can be characterized by evapotranspiration rates double than elsewhere, with implications on precipitation patterns.Leng et al. (2017) developed an irrigation module for the ACME (Accelerated Climate Modeling for Energy) Land Model (ALM) and provided evidence of soil moisture and evapotranspiration increases due to irrigation; an overview of different implications of several parameterized irrigation techniqueswas also given.Wang et al. (2022) carried out SWAT (Soil and Water Assessment Tool) (Arnold et al., 1998) simulations with and without irrigation taken into account, which indicated runoff generation as the mostly impacted mechanism by such practices.Druel et al. (2022) recently implemented an irrigation scheme in the SURFEX/ISBA (SURFace EXternalisée/Interaction Sol Biosphère Atmosphère) land surface model (LSM) (Boone et al., 1999;Habets et al., 2003;Noilhan & Mahfouf, 1996) showing that simulated irrigation amounts alter the magnitudes of some output variables, such as soil moisture and evapotranspiration.Similarly, irrigation schemes have been implemented into other LSMs as well, as the Community Land Model (CLM; Pokhrel et al., 2012), the ORCHIDEE (Organizing Carbon and Hydrology in Dynamic Ecosystems) model (De Rosnay et al., 2003), and the Noah-MP (Niu et al., 2011).Such irrigation schemes proved to be useful tools to investigate irrigation impacts on both water and energy balances at the regional scale (Modanesi et al., 2022).Recent exploitation of remotely sensed data to estimate irrigation amounts has led to the development of the first spatially distributed irrigation water data sets at different scales (Dari et al., 2020(Dari et al., , 2023;;Zhang et al., 2022).More in detail, Dari et al. (2020) developed a data set of irrigation estimates at 1 km spatial resolution over a heavily irrigated area in Northeastern Spain through the SM-based (Soil Moisture-based) inversion approach, firstly proposed by Brocca et al. (2018).The data set, relying on DISPATCH (DISaggregation based on Physical And Theoretical scale CHange algorithm) (Merlin et al., 2013) SMOS (Soil Moisture and Ocean Salinity) and SMAP (Soil Moisture Active Passive) soil moisture, covers an almost-7-year time span.Later on, Dari et al. (2023) implemented the methodology to develop high-resolution and regional-scale irrigation data sets for three major basins: the Ebro basin (in Spain), the Po river valley (in Italy), and the Murray-Darling basin (in Australia).In this case, the SM-based inversion approach was forced with Sentinel-1 soil moisture data obtained through the RT1 (firstorder Radiative Transfer) model (Quast et al., 2019(Quast et al., , 2023) ) for the European sites and with CYGNSS (Cyclone Global Navigation Satellite System) soil moisture (Freeman et al., 2020) for the Australian pilot area.A similar framework was adopted by Zhang et al. (2022) for producing global irrigation estimates relying on coarse resolution soil moisture data from multiple sensors.In a nutshell, such products rely on observations and thus they allow to overcome uncertainties in estimating irrigation impacts through modeling-only approaches (Dorigo et al., 2021), since real irrigation practices can strongly differ from theoretical assumptions adopted in modeling platforms.
In this work, the impacts of irrigation on hydrological state and flux variables over an agricultural portion of the Ebro basin are quantitatively assessed by comparing two SURFEX/ISBA simulations, a reference one forced only by atmospheric input, that is defined as the "natural" simulation, and another one having as an input atmospheric forcing plus explicit irrigation estimates from Dari et al. (2020), being defined as "human-altered" simulation.The physical robustness of the proposed framework, and thus the reliability of the obtained results, has been assessed by comparison with reference satellite evapotranspiration from MODIS (MODerate-resolution Imaging Spectroradiometer) and through an indirect validation carried out by applying the WCA (Si, 2008).

Study Area
The analysis presented in this study has been carried out over a portion of the Ebro basin, in Spain, including six sub-basins enclosed between the Cinca (at the West side) and the Segre (at the East side) catchments.This domain is henceforth defined as the "Ebro-CS" catchment.In general, the Ebro basin in its entirety is characterized by a varying orography, which is reflected on an uneven rainfall distribution across the catchment, with annual rates ranging from less than 500 mm/year in the central valley to more than 1,800 mm/year in the upper mountain areas (Dari et al., 2023).According to the Köppen-Geiger climate classification (Beck et al., 2018) the focus area, that is, the Ebro-CS basin, is mainly characterized by a cold semiarid climate with oceanic climate influences in the upper portion.
Figure 1a provides an overview of the Ebro-CS basin, while the Digital Elevation Model at 25 m spatial resolution derived from the Copernicus EU-DEM v1.1 is shown in Figure 1b.The Ebro-CS basin has a total area of 19,056.4km 2 .It involves the largest irrigated area of the whole Ebro basin, organized in four districts: Urgell, Algerri Balaguer, Pinyana, and Catalan and Aragonese.Different irrigation strategies and techniques are adopted in each district.In the Urgell there is an old system designed for flood irrigation that is still working; in the other districts, sprinkler irrigation is generally adopted for herbaceous crops, while drip irrigation is used for fruit trees (Dari et al., 2021;Paolini et al., 2022).

SURFEX/ISBA Simulations
The impacts of irrigation practices on the natural hydrological cycle have been evaluated by comparing the outputs of two SURFEX/ISBA simulations, one forced with atmospheric data only (henceforth "SURFEX/ ISBA") and another one forced with atmospheric input plus remote-sensing-based irrigation estimates from Dari et al. (2020) added to the liquid precipitation (henceforth "SURFEX/ISBA_RSI").A conceptualization of the experiment carried out in this study is provided in Figure 2.
A 6-year time span ranging from the hydrological year 2011/12 to 2016/2017 has been considered.Both experiments have been carried out through the version 8.1 of the SURFEX modeling platform, obtained from the Center National de Recherches Météorologiques.The ISBA scheme (Boone et al., 1999;Habets et al., 2003;Noilhan & Mahfouf, 1996) for natural surfaces has been used; more in detail, the performed simulations have been carried out through the ISBA-DIF version, in which the soil is discretized in 14 vertical layers.Several validation studies showing the reliability of the ISBA scheme (see, e.g., Calvet et al., 1998;Quintana Seguí et al., 2009;Decharme et al., 2016;Albergel et al., 2018;Le Moigne et al., 2020), of which some focused on the Iberian Peninsula (Cenobio-Cruz et al., 2023;Quintana Seguí et al., 2017), are available in literature.For both experiments, ERA-5 reanalysis data (Hersbach et al., 2020) have been used as atmospheric forcing.The required input variables are hourly 2 m air temperature and relative humidity, wind speed, precipitation, and downward visible and infrared radiation.A regular grid with spatial resolution of 1 km has been adopted.Hence, the coarser resolution atmospheric data (0.25°× 0.25°) has been resampled to this finer resolution grid.Moreover, possible mismatches between the ERA-5 relief and GTOPO30, implemented in SURFEX, have been prevented by applying a lapse rate correction to relative humidity and temperature (Dari et al., 2021).For both simulations, inputs have been provided at hourly resolution and ECOCLIMAP II (Faroux et al., 2013) has been used as a land cover map.In the SURFEX/ ISBA_RSI simulation, spatially-distributed irrigation estimates at 1 km spatial resolution obtained from remotely sensed soil moisture through the SM-based inversion approach (Brocca et al., 2018;Dari et al., 2020Dari et al., , 2022b) ) have been added to the precipitation as a model forcing.The irrigation estimates data developed in Dari et al. (2020) have been used for this aim.The data set consists of irrigation water use covering a period ranging from January 2011 to September 2017 relying on 1 km DISPATCH SMOS (2011)(2012)(2013)(2014)(2015) and SMAP (2016-September 2017) soil moisture.Validation carried out in Dari et al. (2020) highlighted satisfactory performances in reproducing districtaggregated benchmark irrigation amounts; referring to the portion of the data set relying on DISPATCH SMAP, a mean RMSE across districts equal to 0.95 mm/day was found.Similar performances were found for the validation of data relying on DISPATCH SMOS.For more details about the development of the irrigation water use data set, the reader is referred to the Supporting Information S1.
In order to resemble the sprinkler method, the irrigation amounts have been added to the hourly liquid precipitation during a 4-hr time span in the morning (from 07:00 to 11:00) of each irrigation day within the period of interest.Thus, daily irrigation amounts have been evenly distributed over the 4 hr in the morning.The maps provided in Figure 3 show, for each year, the spatial patterns of mean daily irrigation rates during the highestintensity watering season (i.e., from May to August).The maps focus on the irrigated domain.
Four output hydrological variables have been compared in two configurations: surface soil moisture (henceforth SM, obtained as the weighted average of the first two ISBA-DIF layers, whose depth is up to 1 and 4 cm), total evaporative flux (henceforth E, indicating the sum of evaporative contributions from bare soil and vegetation), runoff, and drainage (see Figure 2).As an additional analysis, comparisons in terms of root zone soil moisture (up to 1.5 m) also are provided in the Supporting Information S1.The overall idea is to assess the impact of irrigation practices on the main components of the hydrological cycle by comparing the outputs of the human-altered simulation with the natural ones.2020) are added to the atmospheric forcing as an input to SURFEX/ISBA Land Surface Model.

Comparison With Satellite-Derived Evapotranspiration Data
Simulated E rates with and without irrigation have been compared with reference actual evapotranspiration obtained from MODIS observations at 500 m spatial resolution, which is consistent with the scale of this study; 8day aggregated data belonging to the MOD16A2 product have been extracted for the area of interest through the Google Earth Engine platform (Gorelick et al., 2017).The rationale of this comparison is that if a higher agreement with MODIS evapotranspiration is found when including irrigation amounts in the simulation, then the physical reliability of the experiment carried out is proven.Moreover, Ramos et al. ( 2009) pointed out the sensitivity of remote sensing techniques to crop conditions.Giving the different intrinsic ranges of variability of satellite-retrieved and modeled evaporative fluxes, the comparison has been carried out in terms of anomalies normalized by the standard deviation, En(σ).Thus, monthly time series of both modeled and satellite-retrieved evapotranspiration rates have been scaled by subtracting their temporal mean and dividing by their standard deviation, σ.The choice of considering remotely sensed data as a benchmark is due to the absence of records of flux tower stations or in-situ soil moisture observations over the pilot area during the study period.Nevertheless, a validation based on satellite (hence, potentially tracking irrigation) soil moisture reference data has not been performed because of the willing of avoiding unfair comparisons, as the irrigation amounts injected in the SURFEX/ISBA_RSI simulation rely on remote sensing soil moisture observations (Dari et al., 2020).Nevertheless, as described in the following section, an additional indirect validation has been performed as well.

Wavelet Coherence Analysis
In order to further check the physical robustness of the experiment carried out, the continuous WCA (Si, 2008) has been applied.The WCA allows an in-depth correlation analysis between irrigation input signals and differences in the output hydrological variables examined, as it quantifies the strength of the correlation (i.e., the coherence) and the phase shift (delay) between the studied signals for each point in the time and frequency domains (Rahmati et al., 2020).The possible range for coherence in WCA varies between 0 (with lower correlation between two examined signals in that given time and frequency domains) and 1 (with higher correlation between two examined signals in that given time and frequency domains).The phase shift, also referred as delay, is determined by applying the phase angle function to complex numbers of the cross-wavelet spectrum and then using the following equation to convert it from radians to time resolution (in this case, weeks) (Rahmati et al., 2020): where n represents the period (1/frequency) of the signals provided as one of the outputs in WCA.When performing WCA, one of the signals is provided as the base signal (i.e., the irrigation signal) and the other is provided as the second signal (i.e., the differences in the selected outputs induced by irrigation), and therefore both the phase shift (in weeks) and the phase angle (in radians) refer to the difference in time between two consecutive maximum values of the base and second signals.The possible range for the phase shift varies between ½ n and +½ n, where a negative phase shift means that the base signal lags behind the second signal, while a positive phase shift means the opposite (Rahmati et al., 2020).Lag can be quantified from the phase shift using the following equations, which can help to determine which signal controls the other (Rahmati et al., 2020): Positive lag values mean that the base signal controls the second signal, while negative lag values mean the opposite.According to the above equation, the |lag| value quantifies the length of the time window in which the current state of the base signal/second signal is controlled (depending on the sign of the lag) by the previous states of the other signal that fall within this time window.Results of the WCA analysis can be summarized in coherence-delay maps; they consist of heatmaps quantifying the coherence with arrows whose direction provides information on the correlation signal (right-aligned arrows indicate positive correlation, while left-aligned ones denote negative correlation).For further clarification on the graphical interpretation of outcomes of the WCA analysis, as well as for additional details on the method, the reader is referred to the Supporting Information S1.

Irrigation Impacts on Hydrological Output Variables
Changes induced by irrigation in the spatio-temporal dynamics of the selected components of the water cycle have been assessed.Figure 4 provides, for each hydrological variable considered, the weekly time series produced by the simulation with (SURFEX/ISBA_RSI, indicated by the red line) and without (SURFEX/ISBA, indicated by the black line) irrigation, which is indicated by the green shaded area.For each panel, the differences between the two configurations are quantified by the upper horizontal bar.All the magnitudes are spatially averaged over the irrigated domain.Among the diagnostic variables, SM and E are the two more responsive to the irrigation signal.In fact, the Pearson correlation coefficient, r, between satellite-retrieved irrigation estimates and differences detectable in SM (ΔSM) is equal to 0.75; the r value against differences in E (ΔE) is even higher, equal to 0.86.As a consequence, the largest discrepancies in these two components occur during the highest-intensity irrigation periods, that is, from May to August.A different behavior is observed for the runoff.Irrigation produces an increase in this component which is more evenly distributed in time with respect to SM and E. The differences between the two simulations, ΔRUNOFF, are maximized when runoff peaks occur.The r between ΔRUNOFF and the input irrigation rates is equal to 0.43 that denotes a positive correlation of minor magnitude because the runoff differences between the two simulations increase when the irrigation rates are lower.The dynamics of irrigation impacts on the drainage component is opposite to that for SM and E. Irrigation-induced discrepancies in drainage reproduced by the two ISBA simulations carried out are maximum during the winter season.In fact, inverse correlation between irrigation and its effects on this component, that is, ΔDRAINAGE, does exist (r = 0.49); during the highest-intensity irrigation periods, almost no differences are detected.It is important to point out that all the above-mentioned correlation coefficients are statistically significant at a level of 1%.
Figure 5 provides, for each investigated variable, the climatologies (i.e., the mean daily rates per each month of the year) referring to the two ISBA simulations carried out; results referring to SURFEX/ISBA_RSI and to SURFEX/ISBA are indicated by the red and the black line, respectively.In each plot, the differences between the two climatologies are quantified by the color of the lower circles.As already discussed, SM and E are largely impacted during the highest-intensity irrigation season, which is highlighted by the light gray shaded area in the plots.Conversely, the highest discrepancies in drainage are detected during winter months.The effects of irrigation on the runoff term do not differ significantly across seasons.On average, the month in which SM is altered the most by irrigation is July, with a difference of +0.05 m 3 /m 3 .Likewise, for E also the largest differences are detected in July.In this case, an average E increase of +2.13 mm is found.August is the month in which the highest changes of runoff are observed, quantified in an average increase of +0.36 mm.Nevertheless, changes in this variable referring to the other months are comparable.Finally, Figure 5 shows how low irrigation rates occurring in the study area during the winter season lead to the highest discrepancies in the drainage term, which occur in November and are quantified in +1.15 mm.In relative terms, the maximum increase in average SM and E during July is equal to +30% and +220%, respectively.
The spatial patterns of the differences induced by irrigation on each diagnostic variable during the highestintensity watering season are provided in Figures 6-9, referring to ΔSM, ΔE, ΔRUNOFF, and ΔDRAINAGE, respectively.Each map provides, for each year, the differences between the temporal mean of the considered variable during the period May-August referring to the ISBA simulation with irrigation added to precipitation and its analogous referring to the natural experiment.Results are focused over the irrigated area, as it is the only portion experiencing irrigation-induced changes.SM and E variations between the two configurations investigated (see Figures 6 and 7, respectively) appear more evenly distributed in space rather than irrigation-induced discrepancies in runoff and drainage (Figures 8 and 9, respectively).Referring to ΔSM, a more pronounced pattern can be observed in the East side of the irrigated domain during summer 2016 and 2017, likely due to the change of the source soil moisture information used for estimating irrigation (from SMOS to SMAP, see Section 2.2).
Higher discrepancies in the Eastern portion of the irrigated land are detected in most of the maps showing ΔRUNOFF (Figure 8) and especially ΔDRAINAGE (Figure 9).Such an area corresponds to the Urgell district, over which irrigation estimates are generally higher than over the other districts (Dari et al., 2020).Hence, this surplus of water increases the differences induced in the runoff generation and drainage percolation over the Urgell more than over the rest of the irrigated land, which is reasonable considering that the district is mainly irrigated with flood irrigation.This circumstance opens up the interesting prospect of a potential application of the proposed framework as a tool to assess whether over-irrigation is occurring, which is expected to produce excess runoff and drainage flows.

Comparison Against MODIS Actual Evapotranspiration
Figure 10 provides a comparison between monthly En(σ) calculated for E rates output from the simulations with (red) and without (black) irrigation against corresponding values calculated for actual evapotranspiration amounts from the MOD16A2 product; all the magnitudes are averaged over the irrigated area.The left panel refers to the whole study period (September 2011-September 2017), while the right one provides a focus on the main watering seasons (May-August).In both cases, the inclusion of the satellite-derived irrigation amounts brings to En(σ) values better matching MODIS ones.Considering the whole study period, the RMSE is reduced from 0.83 (referring to the SURFEX/ISBA simulation) to 0.52 (referring to SURFEX/ISBA_RSI).Considering the main watering seasons only, a higher performance increase is observed.In fact, in this case the RMSE referring to the simulation with and without irrigation is equal to 0.63 and to 1.06, respectively.The comparison corroborates the reliability of E fluxes obtained by the simulation with the injected satellite-derived irrigation estimates, thus supporting, in turn, the reliability of the results referring to irrigation-driven alterations in the monitored output variables.The annual coherence of irrigation with ΔSM, ΔE, and ΔRUNOFF shown by right-aligned arrows in the coherence-delay maps (panel a of Figures 11-13) is positive but negative with ΔDRAINAGE, as suggested by left-aligned arrows in the coherence-delay map (panel a of Figure 14).Furthermore, the results show that there is a near-perfect correlation between irrigation and ΔSM between 2012 and 2015, with lag values fluctuating around zero and then transitioning to a lagged correlation where lag values increase up to 2 weeks (Figure 11c).The positive lag values between irrigation and ΔSM mean that the first signal controls the second one from 2015 onward.This means that positive anomalies for soil moisture occur whenever irrigation is applied, which can confirm the physically sound inclusion of irrigation in the simulations.
In contrast to soil moisture, the coherence between irrigation and ΔE is almost lagged in all years studied, with negative lag values before 2015 and positive lag values after (Figure 12e).This simply means that before 2015, E was driving irrigation, while after 2015, irrigation was driving E. This can be interpreted to mean that the study area was an energy-limited system or a less water-limited system before 2015, while it appears to have changed completely to a water-limited system after 2015.This is because prior to 2015, any increase in E required more water, which appears to have been provided by irrigation, but whenever E remained lower, less irrigation was also used.This is also confirmed by the results obtained for irrigation and ΔSM, where there is a near perfect and nonlagged correlation between irrigation and soil moisture before 2015.However, for the period after 2015, the picture is completely different.During this period, E seems to have always been constrained by lower soil water availability, and whenever soil moisture availability was increased by irrigation, the system was able to transpire more water.
The relationship between irrigation and runoff seems to be even more interesting.As shown in Figure 13c, the lag values between irrigation and ΔRUNOFF are mostly positive until 2016 (with a phase shift peak of 4 weeks in 2015), while they turn negative after 2016.This therefore implies that in all the years studied (except for 2016 and 2017), irrigation controls the change in runoff, which means that whenever irrigation was applied, more runoff was generated.
In contrast to all other difference signals examined, the coherence between irrigation and ΔDRAINAGE is almost perfectly anti-correlated (see left-aligned arrows of Figure 14a).The delay values fluctuating around zero (Figure 14c) also confirm this perfect anti-correlation.This simply means that these two signals negatively but instantly respond to each other, which means that theoretically any increase in one signal immediately leads to a decrease in the other.

Discussion
Results showcase that during the main watering season irrigation impacts the most SM and E. This is an expected outcome, coherent with the physical implications induced by irrigation, whose scope is to provide crops with enough water to make evapotranspiration occur at the potential rate.The maximum monthly-averaged increase for E found in this work (equal to +220%, see Figure 5) is reasonable if compared to previous studies (Wei et al., 2013).It is important to highlight that the other two considered fluxes (i.e., runoff and drainage) experience even higher relative increments, but the magnitude of E series and of their range of variation strongly prevail in absolute terms.In fact, the range of E variations is up to 4 times higher than those referring to runoff and drainage during the summer season, as it can be deduced by Figure 7 with respect to Figures 8 and 9.It is noteworthy that lower ratios between E and runoff or drainage would imply that the way irrigation is applied in the model is mostly inefficient.The maximum differences in drainage fluxes detected in winter can be explained by irrigation occurring outside the main watering season to feed fruit trees and greenhouses (Dari et al., 2020).Such irrigation amounts, even though much lower than those delivered to the crops during summer (Dari et al., 2023), produce a surplus of water when the consumption due to evapotranspiration is lower and the soil moisture is closer to saturation.This condition reasonably leads to higher groundwater recharge rates (Jasechko et al., 2014).Moreover, the phenomenon is accentuated by the orography of the pilot agricultural area, which is mostly flat, thus favoring the infiltration process (Flammini et al., 2023).Similarly, adding irrigation to liquid precipitation when there is a lower loss rate due to evapotranspiration and the soil is generally wetter can emphasize runoff peaks (see Figure 4).Nevertheless, winter irrigation amounts are not the only responsible for the behavior observed in terms of drainage dynamics.During the main watering season, the water stored in the soil is greatly depleted, resulting in an almost-absent drainage rate (see Figure 4); irrigation contributes in progressively replenishing the soil profile and increasing root zone soil moisture, which is expected to be maximum right after the end of the irrigation season thus favoring drainage in early winter.Analysis of the response of root zone soil moisture to the injected irrigation, provided in the Supporting Information S1, fully supports the dynamics described above.
In terms of spatial patterns (see Table 1), ΔSM shows the highest correspondence with the irrigation amounts applied.Such a result corroborates the assumption of using soil moisture as a proxy of irrigation at the basis of several studies aimed at detecting, mapping, and quantifying irrigation (Dari et al., 2023;Dari, Brocca, et al., 2022;Elwan et al., 2022;Filippucci et al., 2020;Le Page et al., 2020, 2023;Zappa et al., 2022, to cite a few).
The values of R between spatial patterns of irrigation and ΔE are lower than those referring to ΔSM, circumstance that could be explained by the fulfillment of potential conditions, determining a situation in which ΔE is less responsive to water applied than in a water-limited regime.It is interesting to highlight that for ΔSM and ΔE both minimum and maximum R values occur in the same seasons; the same happens for ΔRUNOFF and ΔDRAIN-AGE also, but in different seasons.Results provided in Table 1 show how surplus of water applied for irrigation does not impact the selected output variables with the same hierarchy everywhere, thus meaning that some of them are strongly influenced by other conditions as well; for instance, soil texture is expected to play an important role for drainage dynamics.
As a general comment, it is important to highlight that the way in which irrigation has been taken into account in the SURFEX/ISBA_RSI simulation (i.e., addition to liquid precipitation) is assimilable to the assumption of having sprinkler irrigation over the whole irrigated domain.Actually, over the pilot districts drip and flood irrigation are practised as well and in the real world different techniques impact on the hydrological processes with different dynamics and magnitudes (Leng et al., 2017).As an example, in case of drip irrigation, there is no water interception from the leaves (Druel et al., 2022).Evaluating the hydrological response to different irrigation techniques is surely among next steps foreseen for this research, but explicit and spatially-distributed information on irrigation type are not easy to obtain.Nevertheless, spatial data sets of irrigation techniques adopted over the Ebro basin and novel physiographical data considering irrigation related processes are currently under development and future simulations will benefit of this valuable information (Barella-Ortiz et al., 2023).Similarly, the integration with spatial information on crop dynamics as for instance land cover maps from the Copernicus Global Land Service or the WorldCereal data set (Van Tricht et al., 2023) would surely help in reducing the gap between LSM simulations and real world irrigation practices.Another potential limitation of the proposed work is the use of a climatological plant phenology; as highlighted by Etchanchu et al. (2017), the exploitation of satellite data to derive the actual phenology can bring benefits in the study of surface fluxes and it represents a possible step forward for this research.
Validation carried out through reference actual evapotranspiration from MODIS shows a better agreement when the most impacted flux among those investigated, namely E, is simulated with the addition of irrigation amounts (SURFEX/ISBA_RSI configuration).The WCA offers the possibility of further corroborating the physical robustness of the experiment carried out, disclosing in-depth relationships between irrigation and hydrologic output variables.For instance, the transition to a more water-limited regime from 2015 onwards highlighted by the WCA for ΔE finds a correspondence in the 2015 European drought, which actually started in The Iberian Peninsula during the spring (European Commission -Joint Research Centre, 2015) and then strongly affected almost whole Europe (Ionita et al., 2017;Laaha et al., 2017;Van Lanen et al., 2016).The negative lag values between irrigation and runoff in 2016 and 2017 mean that the runoff signal controlled the irrigation signal during this period.While this may seem odd at first glance, it most likely highlights a scenario in which surface water availability dropped below critical levels after the 2015 drought year, thus determining a regime characterized by lower resources for irrigation practices, which can, in turn, aggravate streamflow droughts (Van Loon et al., 2022).Another interesting aspect of the relationship between irrigation and runoff in 2016-2017 is that the coherence of these signals appears to strengthen in seasonal cycles (see Figure 13a) during this period, supporting our hypothesis about the role of water resource availability for irrigation.Nevertheless, it is important to remind that 2016 is the year in which the source of remotely sensed soil moisture used to estimate irrigation amounts changed from SMOS to SMAP (see Section 2.2) to overcome issues linked to Radio-Frequency Interference problems affecting SMOS-derived retrievals over the Eastern portion of irrigated lands (Dari et al., 2021).Such a circumstance may have an influence on the observed change in the lag between irrigation and ΔRUNOFF, as well as on the relationships between irrigation and changes in the other variables considered.Referring to the anti-correlation between irrigation and ΔDRAIN-AGE, this cannot be interpreted to mean that when irrigation is higher, drainage is also lower.It might be reasonable to interpret this to mean that whenever drainage was high, irrigation was lower, which indirectly means that whenever moisture in the soil layer was higher, resulting in drainage (e.g., during rainy periods), less irrigation occurred.However, it is worth noting that the anti-correlation is only observed in the annual cycle.While the areas with higher frequency coherence all show a positive correlation, albeit localized at different time windows.These areas are mostly between fall and spring, which means that watering an already well-filled soil has a direct effect on drainage, which is physically realistic.

Conclusions
Water used for agricultural purposes determines imbalances affecting all the components of the water cycle.
This study aims at quantifying such effects by analyzing the differences in SM, E, runoff and drainage output from two SURFEX/ISBA simulations, namely one forced with atmospheric inputs only and another one forced with atmospheric data plus explicit and remote-sensing-based irrigation estimates.By analyzing temporal dynamics, results highlight that SM and E are the two components most responsive to irrigation, as larger discrepancies occur during the highest-intensity irrigation season (see Figures 4 and 5).Referring to the average day of each month, SM and E experience an increase reaching a maximum up to +30% and +220%, respectively, in July.It is important to remark that the impacts quantified in this study refer to the irrigated portion of the pilot basin; obviously, the analogous basin-scale magnitudes are expected to be lower, since the agricultural districts roughly cover 12% of the total area.The injected irrigation determines an increase in the runoff term which is quite evenly distributed in time, while the largest impacts on the drainage process are detected in winter and they are likely due to winter irrigation which, even modest, introduces into the system more water when the consumption due to evapotranspiration is lower.From a spatial point of view, the patterns of irrigation during the main watering season (May-August) are faithfully reproduced by those referring to ΔSM (average R equal to 0.93).A good correspondence is found for ΔRUNOFF (average R equal to 0.86) and ΔE (average R equal to 0.72) as well, while ΔDRAINAGE spatial patterns appear to be strongly driven by ancillary information also (e.g., soil texture).The physical reliability of the experiment carried out has been assessed through comparison between E fluxes simulated with and without irrigation against reference data from MODIS.In terms of monthly anomalies, the inclusion of satellite-derived irrigation amounts leads to a decrease in the RMSE against MODIS actual evapotranspiration of 37% and 41% considering the whole study period and the main watering seasons only, respectively (see Figure 10).The WCA highlights a strong correlation at the annual cycle between injected irrigation and all the studied signals (ΔSM, ΔE, ΔRUNOFF, and ΔDRAINAGE).Moreover, the relationships between irrigation and ΔSM, ΔE, and ΔRUNOFF signals highlight a transition to a more water-limited regime after 2015, when part of Europe experienced an important drought whose track is likely contained in the remotely sensed data used to estimate the irrigation water amounts.This study represents a first step toward the assessment of the hydrological impacts of irrigation practices, which is now possible thanks to spatially distributed irrigation amounts estimated through satellite data (Brocca et al., 2023); the evaluation of irrigation impacts on the water balance at larger scale (e.g., at the regional scale over wide anthropized basins) is a foreseen next step of this research, which will require the coupling with river routing systems for being able to evaluate the impacts on river discharge as well.Another planned activity is the comparison with hydrological impacts determined by irrigation volumes simulated through the ISBA own irrigation scheme (Druel et al., 2022).The final scope of this research is to combine the latest satellite retrievals with modeling platforms to provide insights useful for a rational management of water resources in a context of adaptation to climate change and human exploitation.

Data Availability Statement
The simulations presented in this study have been carried out through the version 8.1 of the SURFEX modelling platform, available at: https://www.umr-cnrm.fr/surfex/spip.php?article387.The atmosphering forcings come from the ERA-5 reanalysis data set (Hersbach et al., 2020), which can be downloaded from the Copernicus Climate Data Store (CDS): https://cds.climate.copernicus.eu/#!/home.The MOD16A2 product is available on the Earth Engine Data Catalog (https://developers.google.com/earth-engine/datasets/catalog/MODIS_061_MOD16A2).
For the satellite-derived irrigation estimates, the data set presented in Dari et al. (2020) has been used.

Figure 1 .
Figure 1.The Ebro-CS domain: (a) overview of the area and sub-basins that constitute the Ebro-CS domain and (b) Digital Elevation Model at 25 m spatial resolution derived by the Copernicus EU-DEM v1.1 with irrigated areas outlined in white.

Figure 2 .
Figure 2. Conceptualization of the experiment carried out in this study.Panel (a) shows the setting of the "natural" simulation, while panel (b) provides the configuration of the "human-altered" simulation, in which irrigation estimates developed in Dari et al. (2020) are added to the atmospheric forcing as an input to SURFEX/ISBA Land Surface Model.

Figure 3 .
Figure 3. Spatial distribution of the mean daily irrigation rates during the period May-August of each considered year.The maps refer to the estimates obtained in Dari et al. (2020).

Figure 4 .
Figure 4. Weekly time series of the considered hydrological variables referring to the configuration with (SURFEX/ISBA_RSI, red line) and without (SURFEX/ISBA, black line) irrigation.In each panel, the weekly differences between the two configurations are quantified by the color of the horizontal upper bars.The considered irrigation estimates are represented by the green shaded area.The light gray shaded areas indicate the highest-intensity watering seasons.All the data shown in this plot are spatially averaged over the irrigated domain.

Figure 5 .
Figure 5. Mean daily surface soil moisture (SM), total evaporative flux (E), runoff, and drainage per each month with (SUREX/ISBA_RSI, red line) and without (SURFEX/ISBA, black line) irrigation.The differences can be inferred by the colors of the circles in the lower part of each panel.The light gray shaded area indicates the highest-intensity watering season.

Figure 6 .
Figure 6.Maps showing the difference between mean daily surface soil moisture (ΔSM) during the period May-August of each considered year output from the simulation with irrigation (SURFEX/ISBA_RSI) and the analogous referring to the simulation without irrigation (SURFEX/ISBA).

Figure 9 .
Figure 9. Same as Figure 6 but referring to ΔDRAINAGE.

Figure 10 .
Figure 10.Comparison of monthly En(σ) values calculated for E rates output from the SURFEX/ISBA (black) and the SURFEX/ISBA_RSI (red) simulation against corresponding values calculated for actual evapotranspiration amounts from the MOD16A2 product.Left panel refers to the whole study period (September 2011-September 2017), while the right one refers to the main watering seasons (May-August).

Figure 11 .
Figure 11.Coherence-delay plots between irrigation and ΔSM signals.ΔSM is defined as the district-aggregated difference between soil moisture series from simulation with (SUREX/ISBA_RSI) and without (SURFEX/ISBA) irrigation.Panel (a) shows the coherence-delay map, panel (b) the scale-averaged coherence, panel (c) the scale-averaged phase shift, panel (d) the time-averaged phase shift, and panel (e) the time-averaged coherence.

Figure 12 .
Figure 12.Same as Figure 11 but referring to ΔE.

Figure 13 .
Figure 13.Same as Figure 11 but referring to ΔRUNOFF.

Figure 14 .
Figure 14.Same as Figure 11 but referring to ΔDRAINAGE.

Table 1
Spatial Correlation, R, Between Mean Irrigation Rates, M.I., Shown in Figure3and Changes inMean SM, E, Runoff, and