Regional impacts of irrigation on the atmospheric and terrestrial water cycle of the Iberian Peninsula in a climate model

This study presents regional simulations over the Iberian Peninsula between 2010 and 2022 with the atmospheric (ICOLMDZ) and land surface (ORCHIDEE) components of the IPSL climate model in a new limited area model configuration (25 km resolution). It uses a recently developed river routing and irrigation scheme based on a water-conservative supply-and-demand approach. Two simulations, with and without irrigation, are compared to isolate the impacts of simulated irrigation on land-atmosphere interactions and the water cycle. First, an evaluation of the simulations is conducted to characterize existing model biases in river discharge, precipitation, evapotranspiration (ET) and surface soil moisture (SSM), and assess whether they can be improved by simulating irrigation. The simulated irrigation is too low in southern Spain because of a lack of available water in the reservoirs, and likely because of the absence of representation of river dams. In northern regions such as the Ebro Valley, the simulated irrigation is more realistic and reduces the biases of river discharge and ET in summer and autumn. In general, SSM is not strongly impacted by irrigation as most additional water is evaporated. Second, atmospheric changes induced by irrigation are studied in summer (JJA). Large atmospheric responses are found over intensely irrigated areas, mainly consisting of a shift in energy partitioning between the turbulent fluxes (increase in latent heat flux and decrease in sensible heat flux, up to 50 W m⁻²), and a lowering of the atmospheric boundary layer (−100 m) and of the lifting condensation level (−250 m). Increases in precipitation are statistically significant only over the mountainous areas surrounding the Ebro Valley, and are closely linked to increases in convective available potential energy. Finally, atmospheric moisture recycling over the Iberian Peninsula is identified by showing that the increase in ET in the presence of irrigation exceeds the amount of water added by irrigation. This is made possible by an increase in precipitation over land, although most of this increase is located in lightly irrigated areas rather than in intensively irrigated areas. These results point to remote atmospheric effects of irrigation and motivate further investigation into land-atmosphere coupling processes in the presence of irrigation in the IPSL model.

1 Introduction

Physical processes at the interface between the soil surface and the lower atmosphere can influence meteorological and hydrological variables at various spatial and temporal scales, contributing to complex feedback loops. In areas where the transition regime described by Budyko (1956); Budyko and Miller (1974) is most frequent, soil moisture (SM) plays a central role in these processes as a strong driver of evapotranspiration (ET), which conditions the partitioning of energy at the interface between land and atmosphere (Betts and Ball, 1995; Seneviratne et al., 2010). The GLACE experiment used general circulation models (GCMs) to identify regions of strong land-atmosphere coupling between SM and precipitation, which were mostly found in semiarid transition regions (Koster et al., 2004). This was confirmed by other modelling studies that also identified various mechanisms through which land surface conditions can impact the atmosphere in these coupling hotspots, and metrics to quantify them (Dirmeyer, 2011; Santanello et al., 2018; Zou et al., 2023; Hay-Chapman and Dirmeyer, 2023). The GLACE-CMIP5 experiments (Seneviratne et al., 2013) extended these conclusions by highlighting the importance of this coupling in the response of climate and the water cycle to global warming in these hotspots (Berg et al., 2015). Moreover, some warm biases of the CMIP5 models (Christensen and Boberg, 2012; Mueller and Seneviratne, 2014) were linked to underestimations of SM (Al-Yaari et al., 2019) and partly attributed to the representation of land-atmosphere coupling processes (Cheruy et al., 2013, 2014; Williams et al., 2016; Ma et al., 2018). This investigation was pursued in the CMIP6 framework with the Land Surface, Snow and Soil moisture Model Intercomparison Project (LS3MIP, van den Hurk et al., 2016). The higher sensitivity of precipitation to changes in soil moisture in transition zones between wet and dry climates was confirmed, as well as the more intense response of 2 m temperature and precipitation to moisture changes under a strong climate change scenario (Catalano et al., 2021). In CMIP6 models, regional analyses found improvements in the representation of precipitation and soil moisture over the USA (Srivastava et al., 2020; Yuan et al., 2021), but excessive land-atmosphere coupling strength in East and Southern Africa (Mwanthi et al., 2024). Substantial model uncertainty on land surface conditions remains in CMIP6 models (Yuan et al., 2021) and the mean climate simulated by current climate models is still highly dependent on the representation of land processes (May, 2023; Zarakas et al., 2024), and of the atmospheric feedbacks they induce (Laguë et al., 2019). The importance of land surface description and its interactions with the atmosphere to properly represent energy partitioning and precipitation has been renewed for this generation of climate models (Devanand et al., 2020; Singh et al., 2024), justifying ongoing research on their understanding and improvement within models.

A great source of complexity in land-atmosphere processes comes from the spatial heterogeneity of SM, which may be derived from the diversity of vegetation, soil types, orographic features, and anthropogenic processes. In particular, in semiarid climates, artificial water inputs in irrigated fields create strong contrasts with the surrounding environment. Observational studies have shown that it increases the latent heat flux and decreases the sensible heat flux in irrigated areas (Rappin et al., 2022; Boone et al., 2025), leading to a moister and cooler atmosphere near the surface (Bonfils and Lobell, 2007; McDermid et al., 2023). In the American Midwest, Nocco et al. (2019) showed that such effects could impact the regional climate and even mask the rise in temperature induced by global warming. Regarding the impacts of irrigation on precipitation, two opposite processes are at play. In general, increases in SM and ET have been associated with increases in precipitation in both modelling and observational studies (Koster et al., 2003; Guo et al., 2006; Wei and Dirmeyer, 2012; Findell et al., 2011), constituting a positive feedback loop (moisture recycling, as presented in Eltahir and Bras, 1996). However, they can also lead to a stabilization of the atmospheric boundary layer (ABL), inhibiting vertical development and convection (Findell and Eltahir, 2003; Ek and Holtslag, 2004). In heterogeneous soil moisture conditions, this can lead to a negative feedback loop where convective rainfall is more likely to occur over drier soil patches. This negative spatial coupling was frequently noticed in observations (Taylor et al., 2012; Klein and Taylor, 2020), although Guillod et al. (2015) emphasized the importance of temporal variability, showing that precipitation over drier soils occurs more often when soils are moister than usual at the regional scale. Specific responses of precipitation to irrigation can be hard to disentangle from other signals in observations, but decreases in irrigated areas associated with concurrent increases in downwind regions have been observed in different contexts (DeAngelis et al., 2010; Alter et al., 2015; Szilagyi and Franz, 2020).

To analyse the atmospheric processes involved and isolate the impacts of irrigation, observational campaigns are often complemented by mesoscale modelling studies. The Great Plains irrigation experiment (GRAINEX, Rappin et al., 2021) measurements and simulations with the Weather Research and Forecasting (WRF) model revealed a lower ABL over irrigated areas (Lachenmeier et al., 2024), although land-atmosphere coupling metrics indicated that irrigated fields were more favourable to convection than non-irrigated ones during the day (Whitesel et al., 2024). A reduction of existing mesoscale slope-induced circulations in the presence of irrigation was identified (Rappin et al., 2022; Phillips et al., 2022), and Lawston-Parker et al. (2023) emphasized the need to account for soil moisture heterogeneities induced by irrigation in modelling studies. Held in the Ebro Valley (northern Spain) in July 2021, the Land surface Interactions with the Atmosphere over the Iberian Semiarid Environment campaign (LIAISE Boone, 2019) highlighted strongly contrasting ABL properties between an irrigated and a rain-fed site, separated by 14 km (Boone et al., 2025). They were partly attributed to a high impact of the contrasts induced by irrigation on the development of the ABL through the morning transition (Brooke et al., 2023). Mesoscale simulations without irrigation over the LIAISE study area exhibited strong limitations over irrigated areas (Jiménez et al., 2025), and near-surface conditions and vertical profiles over the area were greatly improved by representing irrigation (Lunel et al., 2024a; Udina et al., 2024). A lower, more stable ABL was identified over irrigated and surrounding areas (Udina et al., 2024), as well as a weakening of the regional sea-breeze regime due to irrigation (Lunel et al., 2024b). However, irrigation was not associated with a better representation of precipitation during the campaign which was mostly driven by large-scale processes (Udina et al., 2024). In various other contexts, mesoscale simulations with the WRF model showed that the representation of irrigation improves model performance, especially limiting warm biases (Qian et al., 2013; Yang et al., 2017; Qian et al., 2020; Liu et al., 2021). Over the American Great Plains, they identified a lowering of the boundary layer and lifting condensation level (Qian et al., 2013), and a reduction of the summertime precipitation deficit by increasing the frequency of mesoscale convective systems (Qian et al., 2020). In the Colorado River basin, Yang et al. (2017) showed an increase in precipitation linked to irrigation due to moisture recycling over the Sierra Nevada mountains, and similar processes of regional precipitation recycling were identified by Zhang et al. (2025) in the US Corn Belt region. Conversely, a negative feedback loop on precipitation in Saudi Arabia was described in Lo et al. (2021), with moisture convergence leading to increased precipitation in a remote area west of the irrigated region.

However, these mesoscale modelling studies are mostly limited to timescales smaller than a year or focused on a specific field campaign. Moreover, they rely on idealized representations of irrigation, usually forcing soil moisture to remain at field capacity in irrigated areas (as done with Meso-NH in Lunel et al., 2024a). In contrast, this work analyses regional impacts of irrigation on a multiyear timescale and relies on a water-conservative parametrization of irrigation, which allows a study of the full water cycle. The impacts of irrigation on global climate have been studied with such approaches in land-atmosphere coupled simulations, with regional effects mostly visible and analysed in irrigation hotspots such as India, Eastern China and the United States of America (USA) (Sacks et al., 2009; Puma and Cook, 2010; Pokhrel et al., 2012; Cook et al., 2015, 2020; Arboleda-Obando et al., 2025). However, in the CMIP6 framework, only three models included a representation of irrigation, and Al-Yaari et al. (2022) showed that they better captured the trends of several climate variables over irrigated regions. Compared to the mesoscale studies mentioned above, climate models use coarser resolutions and can not explicitly capture the surface heterogeneities created by irrigation. They use parametrizations to account for sub-grid heterogeneities, which have generally been identified as a major challenge for future climate models (Fisher and Koven, 2020). Such modelling choices can greatly influence the response of the atmosphere to changes in surface conditions (Alexander et al., 2022) and partly explain why the impacts of irrigation on climate are still not well constrained (de Vrese and Hagemann, 2018). Additional choices in the parametrization of irrigation itself have also been shown to affect model performance compared to observed surface fluxes (Yao et al., 2022), and to generate contrasting responses in groundwater levels and river flows (Leng et al., 2017). Therefore, efforts are currently being pursued to better account for the impacts of irrigation in climate models, as illustrated by the Irrigation Model Intercomparison Project (IRRMIP), which aims to compare the responses of models to various representations of irrigation (Yao et al., 2023).

This study investigates the regional impacts of irrigation on land-atmosphere coupling variables and the water cycle of the Iberian Peninsula over a decadal timescale. The impact of irrigation on climate has rarely been looked into in this region, even though it has been identified as very dependent on irrigation and subject to future climate disruptions, like the rest of the Mediterranean basin (Cook et al., 2020; Arboleda-Obando et al., 2025; IPCC, 2021). This work is based on two regional simulations (with and without irrigation) run between 2010 and 2022. It leverages a new limited area model configuration developed for the Institut Pierre-Simon Laplace climate model (IPSL-CM) to perform the first regional climate simulations with the land surface and atmospheric components of the global model, ORCHIDEE (Krinner et al., 2005; Cheruy et al., 2020) and ICOLMDZ (Dubos et al., 2015; Hourdin et al., 2020). This setup provides insight on the behaviour of the IPSL-CM components while running with a higher resolution (25 km) and lower computational costs than in global applications, unlocking the possibility to study coupling processes and regional specificities in detail on a multiyear timescale. It also paves the way for future comparisons of the model to measurements of the LIAISE campaign to analyse coupling processes and evaluate the IPSL-CM in irrigated areas.

Section 2 provides a detailed description of the modelling design and reference products used to evaluate the simulations. Section 3 presents an evaluation of irrigation volumes, river discharge, precipitation, ET, SM, and a comparison of the two simulations to isolate the effects of irrigation and analyse continental moisture recycling. Finally, Sect. 4 presents the conclusions, as well as a discussion of the limitations and future developments of this work.

2 Methods

2.1 ICOLMDZOR limited area model

This study uses the atmosphere and land surface components of the IPSL-CM, which has been a regular participant in CMIP, including CMIP6 (Boucher et al., 2020). The simulations are run with prescribed sea surface temperature and sea ice content from the Atmospheric Model Intercomparison Project (AMIP) dataset.

2.1.1 ICOLMDZ atmospheric model

The atmospheric component of the model is the association of the dynamical core DYNAMICO (Dubos et al., 2015), which uses an icosahedral grid, and the LMDZ6A physics version NPv6.2 used for CMIP6, with 79 vertical levels (Hourdin et al., 2020). The physics of the model are run every 15 min and include the following parametrizations:

• a surface layer description based on Monin-Obukhov similarity theory (Monin and Obukhov, 1954). Turbulent fluxes are expressed using the specific humidity and temperature gradients between the surface and the lowest atmospheric level and the aerodynamic resistance. A neutral drag coefficient is computed from the surface roughness, and the formulations of Louis (1979) and King et al. (2001) are used to account for stability conditions;

• an eddy-diffusivity mass flux (EDMF) scheme of boundary layer vertical transfer composed of a turbulent diffusion scheme based on Yamada (1983) with recent improvements described in Vignon et al. (2018), and a thermal plume model for shallow convection (Rio and Hourdin, 2008; Hourdin et al., 2019);

• a mass-flux scheme for deep convection based on Emanuel's scheme (Emanuel, 1991; Grandpeix et al., 2004; Rio et al., 2013), with stochastic triggering (Rochetin et al., 2014a, b);

• a parametrization of the cold pools created below cumulonimbus clouds by reevaporation of convective rainfall (Grandpeix et al., 2010; Grandpeix and Lafore, 2010);

• a large scale condensation scheme based on a statistical distribution of the subgrid total water content, from which the cloud fraction and water contents are derived (Madeleine et al., 2020);

• the radiative transfer model RRTM (Mlawer et al., 1997).

Regional simulations are run using a limited area model (LAM) configuration, which was first used and described in Raillard et al. (2024). Lateral boundary conditions for the LAM are read at each time step of the dynamics, and are taken from ERA5 reanalysis hourly values at 0.25° resolution (Hersbach et al., 2020). The LAM domain comprises 3 zones: a raw forcing zone that contains values directly given by the forcing, a transition zone where the model is nudged towards the forcing with decreasing strength, and a free zone at the centre of the domain where there is no direct influence of the lateral forcing. The outputs of the model are natively on a hexagonal grid but are interpolated to a more traditional longitude-latitude grid of similar resolution, to simplify posttreatment and comparisons to evaluation products.

2.1.2 ORCHIDEE land surface model

Water and energy budgets

The ICOLMDZ LAM is coupled to the ORCHIDEE v2.2 land surface model (LSM) (Krinner et al., 2005; Cheruy et al., 2020). The spatial grid of the LSM is the same as that of the LAM and the timestep is dictated by that of the atmospheric physics (15 mn). ORCHIDEE represents a 2 m soil column discretized along 11 vertical layers of increasing thickness. It computes the coupled water and energy budgets at the surface by simulating, among other processes, latent and sensible heat fluxes, surface runoff and water infiltration in the soil column, with a free drainage condition at the bottom. Each grid cell is assigned the dominant United States Department of Agriculture (USDA) soil texture according to the 5 arcmin resolution map from Reynolds et al. (2000). Vegetation is described using 15 plant functional types (PFTs) from a 0.1° resolution input map based on the Land Use Harmonization 2 (LUHv2) dataset (Hurtt et al., 2020; Lurton et al., 2020). In each grid cell, PFTs are clustered into 3 soil tiles, for bare soil, forests, and low vegetation (including C₃ and C₄ crops). The model computes a separate water budget for each soil tile but only one energy budget per grid cell, using a composite approach with aggregated parameters (roughness length, albedo, aridity factor) to compute surface temperature and turbulent fluxes. The variable for surface soil moisture (SSM) gives water content in the first 10 cm of soil for the entire grid cell. Snow is represented with three layers of varying density and thermal conductivity (Wang et al., 2013). The energy budget of the snowpack is computed, and when a layer's simulated temperature exceeds freezing, it is automatically reset to 0 °C in the snow, and the excess energy is used to melt snow. The resulting liquid water may then percolate through the snowpack or be refrozen. If it reaches the soil surface, it is integrated into the surface water budget, contributing to infiltration and runoff in the same manner as rainwater.

Routing scheme

River discharge is simulated using a routing scheme based on the one described in Ngo-Duc et al. (2007). This scheme solves horizontal water transfers on the digital elevation model (DEM) grid, here upscaled at a 2 km resolution from the MERIT DEM (Yamazaki et al., 2019). Each DEM grid cell contains three linear reservoirs that provide overland flow, groundwater flow and river flow. For each reservoir, the characteristic residence time of water in a grid cell depends on a fixed transfer coefficient and on the slope provided by the DEM. The surface runoff and drainage computed for each ORCHIDEE soil column are interpolated to the routing grid to feed the overland and groundwater reservoirs, respectively. Water from all three reservoirs flows into the river reservoir of the downstream grid cell. In this scheme, river reservoirs constitute the only means of grid-to-grid horizontal water transfer.

Irrigation

Irrigation is modelled using the scheme extensively described in Arboleda-Obando et al. (2024), which is based on a water-conservative supply-and-demand approach. It computes a moisture deficit by comparing SM in the upper layers of the ORCHIDEE soil column (64 cm, corresponding to the root zone) to a target SM described as a fraction of the SM at field capacity. In the default version of the global model, this target is set to 90 % of the SM at field capacity, but this parameter was calibrated on global average to reflect a wide variety of irrigation practices, including flooding in rice paddies. In the Iberian Peninsula, the irrigation methods are less water-intensive, and after calibrating the routing and irrigation schemes, the target value was adjusted to 60 % of the SM at field capacity (see Appendix A). This value had already been identified in Arboleda-Obando et al. (2024) as a more suitable value for irrigated areas with other methods than paddy irrigation, confirming its relevance for the study area. To avoid computing irrigation requirements on grid cells without plants, the SM deficit is set to zero if the leaf area index (LAI) is below a given threshold LAI_min=0.1.

At each time step, the SM deficit is weighted by the irrigated fraction of each ORCHIDEE grid cell, based on the historical irrigation dataset (HID, Siebert and Döll, 2010) at 5 arcmin resolution, to define the irrigation requirement. This demand is then compared with the available water in the three reservoirs of the routing module. Only 90 % of the volume can be withdrawn from each reservoir to maintain a minimum environmental flow. This parameter value was selected after a sensitivity analysis in Arboleda-Obando et al. (2024) and found to have a much lower impact on withdrawals than the target SM parameter, so it was not tuned again for this study. Water is withdrawn preferentially from either surface water (overland and rivers) or groundwater depending on the nature of irrigation equipment, defined by the map of Siebert et al. (2010). It must be noted that although the original irrigation scheme from Arboleda-Obando et al. (2024) included the possibility to withdraw water from neighbouring grid cells to represent adduction systems, this option is not yet compatible with the new version of the routing scheme and was not used in this study. As a final step, the amount of water withdrawn from the reservoirs is added at the top of the ORCHIDEE soil column at the next time step for infiltration.

Figure 1 Altitude (m) over the simulation domain on the native hexagonal grid.

Figure 2 Stations used for river discharge evaluation and major rivers of the Iberian Peninsula, with river basin district borders from the WISE WFD Reference Spatial Datasets reported under Water Framework Directive. Grey diamonds show river dams from FAO – Aquastat (2009).

2.2 Simulation setup

The Iberian Peninsula was chosen as the study region since it is a hotspot for irrigation in Europe, particularly in the Ebro Valley (north-eastern Spain, see Fig. 2). In addition, it is clearly bounded by the Atlantic Ocean, the Mediterranean Sea, and the Pyrenees mountain range in the north, which ensures that all watersheds are fully included in the simulation domain, making the regional study of the continental water cycle rather straightforward. Moreover, although it was not identified by Koster et al. (2004) as a region of strong coupling, probably owing to its small size relative to the resolution of the model used and to other identified hotspots, its semiarid climate should make it very sensitive to the influence of land-atmosphere coupling processes. Therefore, modifications of the SM heterogeneity patterns induced by irrigation are expected to have a clear influence on the local and regional climates.

To study this region, the LAM simulation domain is a hexagon (centred at 40.4° N, −3.7° E) with a radius of 1500 km (Fig. 1). The radius is composed of 60 grid cells, so the diameter of each cell is 25 km. The simulations are run for 13 years from 2010 to 2022, which enables capturing some interannual variability of the current climate, making the averages less sensitive to anomalies or biases from any single year or extreme event. Before this 13 year simulation, 3 years were run as a spin-up to allow for the hydrological and vegetation variables to reach a satisfactory equilibrium (see Appendix B for more details). Two simulations are run with the same setup except for the inclusion of irrigation. They are referred to as irr (with irrigation activated in ORCHIDEE) and no_irr (without irrigation). Considering the focus on land-atmosphere coupling, the analysis is restricted to the Iberian Peninsula, thus excluding ocean grid cells and the land grid cells north of the Pyrenees, which flow to France. To prevent mismatches between the simulations and the evaluation datasets because of different coastlines, the analysis is further reduced to grid cells where the continental fraction is greater than 95 % .

Rustemeier et al. (2020)Martens et al. (2017)Miralles et al. (2025)Dorigo et al. (2017)Gruber et al. (2019)Dari et al. (2023)Fekete (2003)

Table 1 Reference datasets used for evaluation.

2.3 Evaluation datasets

The simulations are evaluated against the monthly mean values of the reference gridded products listed in Table 1. Precipitation data from the Global Precipitation Climatology Centre (GPCC) Full Data Monthly Product Version 2020 (Rustemeier et al., 2020) is used. This reanalysis product provides precipitation data until 2019 over land on a 0.25° × 0.25° grid using in situ rain gauges. For ET, the Global Land Evaporation Amsterdam Model (GLEAM) dataset is used, in its fourth version (Miralles et al., 2025). This product computes ET using a large set of input variables obtained from reanalyses as well as in situ and satellite observations. Monthly values at 0.25° resolution are used, initially given in mm per month but converted to mm d⁻¹. GLEAM4 is available until 2022 but since precipitation and ET are evaluated jointly, it is only used over the availability period of GPCC data (2010–2019). The evaluation of surface soil moisture uses the European Space Agency Climate Change Initiative (ESA CCI, Dorigo et al., 2017; Gruber et al., 2019) COMBINED product, which merges multiple active and passive satellite remote sensing instruments in the microwave domain (Preimesberger et al., 2021) to provide daily SSM (in m³ m⁻³) at 0.25° resolution. This product contains low-quality flags (e.g. in the presence of snow, dense vegetation, or radio-frequency interference in the measurements) and an estimated uncertainty of the measurement. A screening of the data for the Iberian Peninsula over the period 2010–2022 was conducted, following the procedure described in Mizuochi et al. (2021). Data records with quality flags other than zero were excluded (representing 19 % of daily data records), as well as those with an uncertainty larger than 0.6 m³ m⁻³ (representing only 0.1 % of the remaining data records). The remaining daily values were then averaged to monthly time steps. Direct comparison of the CCI SSM product with the ORCHIDEE model is limited by the inherent differences in soil representation between ORCHIDEE and the LSM used to scale the satellite measurements (GLDAS-Noah), as analysed in Raoult et al. (2018). Therefore, a spatio-temporal normalization was applied following Mizuochi et al. (2021); Polcher et al. (2016) to compare statistically normalized values rather than absolute values:

$$\begin{array}{}\text{(1)}& {\mathrm{SSM}}_{\mathrm{norm}}={\displaystyle \frac{\text{SSM}-\stackrel{\mathrm{‾}}{\text{SSM}}}{{\mathit{\sigma }}_{\mathrm{SSM}}}},\end{array}$$

where SSM_norm is the normalized SSM, and $\stackrel{\mathrm{‾}}{\text{SSM}}$ and σ_SSM are the mean and standard deviation, respectively, of the available SSM values over the Iberian Peninsula for all the considered monthly time steps. Similarly to ET, although the product is available until 2022, the evaluation was restricted to 2010–2019 to match the evaluation of precipitation, and the normalization is conducted over this period.

Irrigation is evaluated in the Ebro Valley region using a high-resolution remote sensing product from the European Space Agency (ESA) Irrigation+ project (Dari et al., 2023). This product estimates irrigation with a soil moisture based approach using satellite measurements from Sentinel-1, and provides data for three intensely irrigated areas: the Ebro Basin in Spain, the Po Valley in Italy and the Murray-Darling Basin in Australia. From 2016 to 2020 in the Ebro Basin, the median values of the RMSE, Pearson correlation coefficient r and bias are 12.4 mm over 14 d, 0.66, and −4.62 mm over 14 d, respectively. Weekly values are aggregated to a monthly average in mm d⁻¹. The simulated river discharge is evaluated against monthly observation data from discharge stations of the Global Runoff Data Center (GRDC, https://grdc.bafg.de, last access: 14 May 2026, Fekete, 2003). Stations were positioned on the MERIT DEM grid with tools presented in Polcher et al. (2023), which use the GPS position of the stations as well as the upstream catchment area to find the most appropriate grid cell for comparison with the observations. The 18 selected stations with available data over the simulation period and an adequate position on the DEM grid are described in Table 2 and shown in Fig. 2. Most stations have available data from January 2010 to September 2017, and river discharge was therefore evaluated using the first eight years of simulation (2010–2017).

Table 2 Characteristics of river discharge stations used for evaluation. Stations marked with ^* are the largest of the five major basins of the Peninsula, and are shown in Fig. 5.

2.4 Statistical significance

To assess the influence of irrigation on the model, the two simulations (irr and no_irr) are compared and a statistical significance test is used to filter out differences that may be the result of natural variability only. For the maps of changes induced by irrigation in Sect. 3.4, a Student t test is used to assess for each grid cell whether the mean difference between the two simulations (irr − no_irr) significantly differs from 0, with a p value of 0.05 as the limit to reject the null hypothesis. Grid cells with nonsignificant changes are partly hidden with hatches.

Figure 3 Simulated irrigation and its drivers. Input maps of (a) grid cell irrigated fraction (%, derived from Hurtt et al., 2020) and (b) the share of irrigation equipment for surface withdrawals, as opposed to groundwater withdrawals (%, derived from Siebert et al., 2010). Annual means (2010–2022) of (c) simulated irrigation requirement (mm d⁻¹), (d) irrigation (mm d⁻¹), and relative changes (irr − no_irr, %) in water volumes in (e) groundwater and (f) river reservoirs.

3 Results

3.1 Simulated irrigation

The computed irrigation demand (Fig. 3c) is highly dependent on the irrigated fraction (Fig. 3a) and much greater than the applied irrigation (Fig. 3d). This shows that irrigation is often constrained by water availability, with clear regional differences. Indeed, irrigation is much greater in the northern regions (Ebro and Douro river basins) than in southern regions (Guadiana and Guadalquivir basins) even though these regions have similar levels of irrigation demand (Fig. 3c, d). As shown in Fig. 3b, southern regions (Guadalquivir Basin, upper Guadiana Basin) are more dependent on groundwater equipment for irrigation water withdrawals than the Ebro Basin, where withdrawals are taken mainly from surface water (overland and river reservoirs in the model). Considering that the groundwater reservoir is much more depleted in the presence of irrigation than the river reservoir is (Fig. 3e, f), it is not surprising that the irrigation requirement cannot be met in these regions as much as it is in the north. This depletion can be explained by the fact that the groundwater reservoir can only be filled with drainage in the grid cell, whereas the river reservoir can be fed from upstream grid cells. It is also important to note that ORCHIDEE does not model deep groundwater storage, which is an important source of water for irrigation in southeastern Spain (Custodio et al., 2016).

Figure 4 Evaluation of simulated irrigation from January 2016 to July 2020 over the Ebro Valley. (a) Mean seasonal cycle of irrigation for the irr simulation and the remote sensing product (Ebro_estimate, mm d⁻¹). Annual mean values are shown in the legend. (b) Mean bias compared with the product (mm d⁻¹). The simulation outputs are interpolated on the grid of the remote sensing product.

Figure 5 Impacts of irrigation on river discharge. Mean seasonal cycle of river discharge (m³ s⁻¹) in observations (black) and the no_irr (red) and irr (blue) simulations at five stations: (a) Tortosa (Ebro), (b) Tore (Douro), (c) Talavera (Tagus), (d) Azud de Badajoz (Guadiana), and (e) Mengibar (Guadalquivir). Annual mean values are shown in the legend. A mask is applied to the simulations to filter out months without corresponding observation data.

In the Ebro Basin, simulated irrigation is evaluated using the irrigation remote sensing product from Dari et al. (2023) with good overall performance, particularly in summer (Fig. 4a). In winter, the model simulates almost no irrigation since it requires a minimum LAI to be activated, whereas the product shows irrigation all year long, which can be explained by the presence of winter crops not represented in the model. A delay of the summer peak in the model compared with the product is noticeable, but the model might not be very far from actual irrigation since this product was found to be slightly ahead of actual irrigation based on benchmark volumes in some districts of the Ebro Valley (Fig. 5 in Dari et al., 2023). Spatially, the remote sensing product shows greater irrigation on the hillslopes than in the thalwegs, whereas the model simulates the opposite, with more intense irrigation next to the large rivers (Ebro, Segre, Cinca). The resulting bias pattern (Fig. 4b) can be explained by the fact that in the model, water is mainly withdrawn from the river reservoir in this region, which is much greater in grid cells holding a large river than in upper areas of the valley. In reality, infrastructures such as the Canal d'Urgell in the Segre basin (Farran, 2024), enable gravity irrigation of hillslopes by diverting water from large rivers to neighbouring croplands. Including a representation of water adduction in the irrigation scheme by enabling withdrawal from adjacent grid cells could be a way to improve this bias. Overall, the spatial biases of the simulated irrigation offset each other relatively well. Averaged over the subdomain where the satellite product provides values, the simulated irrigation is 0.20 mm d⁻¹ while the product estimates it at 0.23 mm d⁻¹.

3.2 Impacts of irrigation on river discharge

Figure 5 shows the average seasonal cycle of river discharge for the five largest rivers of the Peninsula (Ebro, Douro, Tagus, Guadiana, and Guadalquivir), for the two simulations and the GRDC observation data. The station with the greatest river discharge was selected for each river, to reflect integrated impacts of irrigation over the basin. A similar figure with all eighteen stations is presented in (Fig. C2). In most cases, the model shows a slight delay and a large overestimation of river discharge compared with observations, particularly in winter and spring. Snow melt does not seem to play a large role in the delay since it occurs rather concurrently with snowfall (Fig. F1a). These errors are most likely related to the precipitation biases of the simulations, described in the next section, and to the lack of river dams in the model, which have a strong impact on actual discharge, given their high density in the Iberian Peninsula (Fig. 2, Sabater et al., 2022; Morán-Tejeda et al., 2012; Lobera et al., 2015). In particular, the model overestimates the winter and spring discharge, a period when water is stored in dam reservoirs, which reduces the actual river flow. The presence of dams also leads to unnatural seasonal cycles in the observations if river discharge is artificially increased in summer by the release of stored water (Fig. 5e).

Table 3 Station mean discharge (obs, m³ s⁻¹), discharge bias (for the no_irr and irr simulations, in %), and change in evaluated metrics (irr − no_irr) for the RMSE (relative change in %), Pearson correlation coefficient r, Nash-Sutcliffe efficiency and Kling-Gupta efficiency. Model performance improvements when using irrigation are shown in bold. Stations marked with ^* are the largest of the five major basins of the Peninsula, and are shown in Fig. 5.

The simulation of irrigation used cannot improve either of these aspects since it does not include a specific reservoir to store water. Its impact becomes noticeable in spring, with water withdrawals resulting in lower discharge during summer and autumn, generally leading to a much better match with observations (Fig. 5). As shown in Table 3, for all fifteen stations where the mean bias is positive in the no_irr simulation, this bias is reduced in the presence of irrigation (in three cases it becomes negative but the absolute value is reduced). However, for the three stations where the no_irr bias is negative (no. 10, 13, 17), it is worsened in the irr simulation. On average, the irr simulation exhibits clear improvements of the mean bias (−41.8 %) and root mean square error (RMSE, −7.12 %). The Pearson correlation coefficient is 0.56 on average in no_irr and is slightly improved (+0.02), with mostly small changes except for stations 6 and 8 (+0.09 and +0.08). Improvements are also observed for Nash-Sutcliffe efficiency (NSE, +0.67) and Kling-Gupta efficiency (KGE, +0.2), mostly as a consequence of improvements in the mean bias. However, only eight out of eighteen stations have a positive value of NSE and KGE values in the no_irr simulation (not shown), limiting the relevance of this average increase. In particular, the average NSE value is strongly influenced by a few stations (no. 5, 8, 15) with initial NSE values below −10.

Overall, the performance is clearly improved for 12 stations but partly degraded for 6 stations, although three of them (no. 10, 13, 17) have a very small average discharge. This may explain why small changes in the model can lead to large changes in performance and limits their relevance compared with larger stations. If only stations with an average annual discharge greater than 10 m³ s⁻¹ are considered, irrigation improves performance in nine out of twelve cases, and degrades it for three stations. One station (no. 9) exhibits an unexpected increase in the spring discharge peak which originates mostly from a single year (2011, see Fig. C1) and worsens an already-existing 250 % bias in this season. The other two (no. 4 and 5) are close to the Pyrenees Mountains and present very strong biases in both simulations (300 % overestimation in spring, unexpected double peak in March and June, Fig. C2). These discrepancies are likely related to biases in precipitation (discussed hereafter) and were not positively impacted by irrigation, apart from the mean bias. It must also be mentioned here that in the Pyrenees, snowfall represents up to 35 % of annual precipitation in the simulation (Fig. F1b) and largely contributes to river discharge in these two stations via snow melt.

Figure 6 Evaluation of simulated precipitation (P) and evapotranspiration (ET) over the Iberian Peninsula continental subdomain, from 2010 to 2019. (a) Mean seasonal cycle of P (mm d⁻¹) for the two simulations and GPCC product (annual mean values are shown in the legend), annual mean P bias of the (b) no_irr and (c) irr simulations relative to the GPCC product (mm d⁻¹). (d) Mean seasonal cycle of ET (mm d⁻¹) for the two simulations and GLEAM4 product (annual mean values are shown in the legend), annual mean ET bias of the (e) no_irr and (f) irr simulations relative to the GLEAM4 product (mm d⁻¹). (g) Mean seasonal cycle of normalized SSM for the two simulations and CCI product, annual mean normalized SSM bias of the (e) no_irr and (f) irr simulations relative to the CCI product.

3.3 Evaluation and influence of irrigation on precipitation, ET, and SSM

The simulated precipitation, ET, and normalized SSM are evaluated from 2010 to 2019 using the GPCC, GLEAM4 and CCI products, respectively (Fig. 6). On average over the domain, the two simulations present very similar seasonal cycles of precipitation. The model is in good agreement with GPCC until June (Fig. 6a), but presents a strong underestimation of precipitation in summer, followed by a delayed and overestimated peak in autumn. This peak is mostly the consequence of precipitation in mountainous regions such as the Pyrenees and the northern coast (see Fig. D1 for seasonal comparisons to GPCC) and is dominated by rainfall rather than snowfall (Fig. F1a). The delay in precipitation likely contributes to the biases of river discharge winter peaks visible in Fig. 5. This seasonal cycle is largely representative of the whole peninsula, although some spatial disparities persist. The two simulations exhibit very similar spatial patterns of annual mean precipitation, with a strong overestimation in elevated areas (Fig. 6b, c), which is a known bias of climate models (Arjdal et al., 2024; Adhikari et al., 2024). This is partly compensated by smaller underestimates of precipitation over large neighbouring areas, as seen in the Ebro Valley.

Both simulations match the GLEAM4 ET product well from January to May but underestimate ET for the rest of the year, particularly in summer (Fig. 6d). As expected, ET increases when irrigation is accounted for, particularly from May to September, which is the period where vegetation is the most developed and irrigation is the greatest. This partially alleviates the dry bias, but ET remains underestimated, even in the irr simulation. No similar patterns of biases between ET and incoming radiative fluxes were identified, and the remaining ET bias can be related to the underestimation of precipitation in the southwestern part of the Peninsula and in plains, such as the northern Ebro Valley (Fig. 6b, e). The biases in ET and precipitation are not consistent in the north-west of the Peninsula, which is explained by seasonal differences. Indeed, the overestimation of precipitation mostly occurs in winter and spring and does not strongly affect ET since soils are already wet in this season, whereas the underestimation of ET is dominated by summer (see appendix Figs. D1 and D2). Along large rivers, the ET underestimation almost disappears in the irr simulation (Fig. 6f). In contrast, the ET bias remains significant in lightly irrigated grid cells such as hillslopes, which is consistent with the limits of simulated irrigation described in Sect. 3.1 (Fig. 4).

The seasonal cycle of the normalized SSM matches the CCI product in both simulations (Fig. 6g). It is very similar in the two simulations, although in absolute value, SSM is slightly higher in summer in irr, and annual increases up to 11 % are visible in the most irrigated grid cells (Fig. E1a–b). Compared to the increase in ET, this indicates that the water added by irrigation is mostly evaporated, as will be further illustrated in Sect. 3.5 and Fig. 8. This confirms that the transition regime described by Budyko (1956) is dominant in the region. Spatial biases in normalized SSM compared with the CCI product are very similar in the two simulations and match the precipitation underestimation in the Ebro valley, Guadalquivir river mouth, and overestimation in the Pyrenees, but show a contrasting pattern in the West compared to precipitation and ET, due to a positive bias in summer (Fig. D3).

Figure 7 Summer irrigation and its impacts (JJA, 2010–2022). (a) Irrigation (mm d⁻¹) and 10 m wind (irr simulation). Mean changes in the presence of irrigation (irr − no_irr): (b) evaporative fraction, (c) 2 m temperature (K), (d) atmospheric boundary layer height (m), (e) lifting condensation level (m), (f) precipitation (mm d⁻¹), (g) convective available potential energy (J kg⁻¹), (h) moisture convergence (mm d⁻¹). Hatching indicates areas where the change is not statistically significant.

3.4 Atmospheric impacts of irrigation in summer

Here, the impacts of irrigation on atmospheric variables are studied in summer (JJA) since it is the season with the highest levels of simulated irrigation, with seasonal mean values of up to 1.5 mm d⁻¹ in the most intensely irrigated grid cells (Fig. 7a). In the presence of irrigation, the simulated latent heat flux (LE) increases across the entire Iberian Peninsula, by up to 50 W m⁻² in the Ebro Valley. As expected from the surface energy partitioning, this is compensated by a decrease in the sensible heat flux (H), which is almost equivalent in irrigated areas and leads to large increases in the evaporative fraction ($\text{EF}=\frac{\mathit{LE}}{\mathit{LE}+H}$) shown in Fig. 7b. When the sensible heat flux decreases, less energy is transmitted from the surface to the air, leading to a decrease in the 2 m air temperature which spatially matches the increase in EF. The order of magnitude remains low over most of the Peninsula, with the most important changes reaching −0.35 K in the Ebro Valley (Fig. 7c). The decreases in sensible heat flux and temperature also lead to a more stable boundary layer over most of the peninsula, but mostly in intensely irrigated areas where it is lowered by 100 m (Fig. 7d). Moreover, the presence of irrigation results in a moister lower atmosphere, with an average specific humidity over the Peninsula increasing by $\mathrm{2.8}\times {\mathrm{10}}^{-\mathrm{4}}$ kg kg⁻¹ in summer (+3.4 %) and maximal local increases in the Ebro Valley of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$ kg kg⁻¹ (+10 %). Since air temperature changes in the atmospheric column are rather small, the lowering of the lifting condensation level (LCL) reflects this atmospheric moistening very well. It is most marked in the Ebro Valley, where the LCL is lowered by 250 m (−13 %) in the most intensely irrigated grid cells, and remains significant even in areas where irrigation is low (Fig. 7e).

The lowering of the ABL and LCL theoretically favour opposite effects on precipitation. On the one hand, a lower and more stable ABL inhibits vertical mixing and convection, reducing the likelihood of cloud formation and deep convection initiation. On the other hand, if the LCL is lower, air parcels do not need to be lifted as high to condense, which increases the likelihood of cloud formation. Over the most intensely irrigated areas, ABL stabilization seems to dominate and inhibit convective development since no significant change in precipitation is observed. However, mountainous areas surrounding the Ebro Valley show significant increases in precipitation (Fig. 7f). This can be understood because ABL stabilization remains weak in these zones whereas humidity can still be increased if moisture is advected (Fig. 7d, e). In particular, the dominant wind patterns in the Ebro Valley (Fig. 7a) indicate that the additional atmospheric moisture from irrigated areas is driven towards the valley slopes, which is consistent with the increases in moisture convergence (Fig. 7h) and precipitation over the Pyrenees. The competing interactions of ABL stabilization and atmospheric moistening are reflected by the increases in convective available potential energy (CAPE) which are most important in elevated areas around the valley (Fig. 7g), where increases in precipitation are significant.

In other seasons (see Appendix G for equivalent figures), almost no significant impacts are visible in winter (DJF) since irrigation is very low and no clear pattern emerges in other variables. In spring (MAM) and autumn (SON), the impacts resemble those of JJA with lower values for all variables and less statistical significance over the domain, and a greater amplitude in SON than MAM. Precipitation increases in the Pyrenees are also present in these seasons, but the area where changes are significant is reduced, and no clear signal is found in the southern rim of the Ebro Valley. The only difference in spatial pattern is found for CAPE in the spring, with increases at the centre of the Peninsula rather than the mountainous areas, although the amplitude of the change remains almost ten times lower than in summer.

Figure 8 Domain-averaged influence of irrigation on monthly changes in ET and P (2010–2022). Each data point corresponds to the average value over the Iberian Peninsula continental domain for a single month of simulation (156 data points for 12 months over 13 years). The average amount of water added by irrigation in the irr simulation (mm d⁻¹) is plotted against the average change (irr − no_irr) in (a) ET (mm d⁻¹) and (b) P (mm d⁻¹). The data points for the winter months are all concentrated around (0,0) for both figures because of very small irrigation volumes and changes in ET and P during this season.

3.5 Atmospheric moisture recycling over the Iberian Peninsula

On average over the continental domain, the monthly change in ET in the presence of irrigation is well correlated with the amount of water added by irrigation and even exceeds it, particularly in summer months (the orange JJA data points in Fig. 8a are all on or above the 1:1 line). In the simulation, ET is constrained by available water, and almost all the water added by irrigation is evaporated or transpired, meaning that this additional increase in ET comes from an additional input of water into the soil. This is explained by a systematic increase in precipitation over the domain (all the data points are on or above the x axis on Fig. 8b). This increase is also roughly proportional to the amount of applied irrigation, although the correlation is weaker than that for the increase in ET, and its values remain lower than the amount of water added by irrigation. Therefore, it appears that irrigation contributes to an increase in atmospheric moisture, and that a part of this moisture is recycled as continental precipitation, which can then be reevaporated.

Figure 9 Changes in the atmospheric moisture budget for subdomains with different irrigation intensities. Each bar plot shows the annual mean irrigation in the irr simulation (2010–2022, mm d⁻¹) alongside the changes (irr − no_irr) in ET (mm d⁻¹) and P (mm d⁻¹) averaged over distinct subsets of the domain: (a) low irrigation, grid cells with an annual average irrigation lower than 0.06 mm d⁻¹, (b) medium irrigation, those where it is between 0.06 and 0.12 mm d⁻¹, (c) high irrigation those where it is higher than 0.12 mm d⁻¹, and (d) the Iberian Peninsula includes all 3 subsets. The three subdomains are shown in (e).

Table 4 Subdomains of different irrigation intensity.

To look further into this recycling, three subdomains were defined, namely, low, medium and high irrigation areas, on the basis of the mean irrigation thresholds given in Table 4. In the map of simulated irrigation (Fig. 3d), the low irrigation domain corresponds to the first colour bin (yellow), the medium irrigation domain to the second bin (light green), and the high irrigation domain to the eight other bins. The three subdomains are also shown distinctly in Fig. 9e. On average, the increase in ET is slightly superior to irrigation for each subdomain (Fig. 9). However, the increase in precipitation is more than twice as large for the low irrigation subdomain than for the medium and high irrigation subdomains. Since irrigated areas are mostly in plains and valleys, this result is consistent with the increase in precipitation already described over mountainous areas in summer (Fig. 7f). It points towards a nonlocal moisture recycling, with atmospheric moisture transfer from intensely irrigated areas to neighbouring lightly irrigated areas, meaning that a significant part of the additional rainfall does not occur on irrigated crops. Over the entire Iberian Peninsula, the increase in precipitation represents 25 % of the irrigated volume, whereas the increase in ET amounts to 112 % of irrigation.

4 Discussion and conclusions

This study analyses the regional impacts of simulated irrigation on land-atmosphere coupling variables and the water cycle over the Iberian Peninsula. It uses a regional model at 25 km resolution, with the same physics as the global IPSL-CM, to better understand the strengths and limits of its parametrizations in representing the impacts of irrigation.

It first shows that the ORCHIDEE irrigation scheme simulates realistic values from April to September in areas where surface water withdrawals are most important, such as the Ebro Valley, although it partly relies on the compensation between an overestimation of irrigation in the thalwegs and an underestimation in the hillslopes. However, the irrigation scheme cannot represent winter irrigation, or satisfy irrigation demand in southern regions, where actual irrigation is more dependent on groundwater pumping and river dams, due to low available volumes in rivers and groundwater routing reservoirs. Ongoing developments to add river dams into the ORCHIDEE routing scheme (Baratgin et al., 2024) could very likely improve this aspect by representing interseasonal water storage, making more water available in summer for irrigation. Explicit dam representation could generally improve the fidelity of the model by limiting the winter and spring overestimation of river discharge in anthropized areas (since water would be stored in the dam reservoirs during this season instead of flowing in the rivers), and by accounting for environmental flow regulations that increase discharge in summer using water stored in dam reservoirs (Sadki et al., 2023). This would disentangle the impacts of irrigation on river discharge from those of dam operation, and likely require a new parameter tuning of the irrigation scheme. In its current state, the irrigation parametrization reduces river discharge and enables better agreement with observations, but since it is only active when the LAI is above a defined threshold, these impacts are mostly visible in summer and autumn. Future work with a looser activation threshold for irrigation could help to represent winter crop irrigation, although it is not expected to have as significant an impact on discharge as an explicit dam representation since simulated irrigation demand would still remain low in winter. Nevertheless, biases induced by precipitation and the melting of snow in mountainous regions are very likely to remain major drivers of discharge biases, largely independent of irrigation or dam representation.

The simulation of precipitation and ET over the Iberian Peninsula is satisfactory in winter and spring, but this study highlights underestimates in summer and contrasted spatial patterns with positive precipitation biases in elevated regions and negative biases in plains. ET underestimation is partly improved by simulated irrigation, but remains present on average and over most of the domain. These linked biases might be improved with a different simulation setup, particularly in the lateral forcing. Preliminary analyses (not shown) revealed an abnormal behaviour of the model in the transition zone between the ERA5 forcing zone and the central free zone, which was attributed to discrepancies between the physics used in the model and in the reanalysis. This resulted in precipitation underestimations throughout the entire simulation domain, which were largely improved by using a larger domain for the simulations presented here. A good lead for future works would be to use lateral forcing from global simulations of the ICOLMDZ model or nested LAM simulations rather than a reanalysis, but these options are not yet technically available. The precipitation biases can also be due to structural flaws of the IPSL-CM parametrizations, as mentioned for mountain precipitation in Arjdal et al. (2024), and improving them might require more work in the modelling of radiative processes, shallow and deep convection (whose tuning often focuses on tropical regions), or surface processes (roughness, albedo, components of ET). This highlights the fact that the results of this study are necessarily limited by the modelling choices, uncertainties, and biases of the IPSL-CM, and therefore remain largely model-specific.

The atmospheric impacts of irrigation are analysed in more detail in summer, since it is the season with the largest irrigation values and the most significant response for all variables of interest, although it is the driest season, with very little precipitation. In JJA, the strong response of turbulent fluxes to irrigation leads to cooling and moistening of the lower atmosphere and significantly affects its structure (LCL and ABL height), with stronger effects on intensely irrigated regions, which is consistent with the findings of Rappin et al. (2022). In contrast, significant increases in precipitation are mostly detected in lightly irrigated mountainous areas surrounding the highly irrigated Ebro Valley. This points to a dominant effect of ABL stabilization, described by Findell and Eltahir (2003); Ek and Holtslag (2004), in intensely irrigated areas, and remote effects of atmospheric moistening as in DeAngelis et al. (2010), Lo and Famiglietti (2013), Yang et al. (2017). An improved representation of winter and spring irrigation could either allow to generalize the following results or to identify different responses to irrigation under moister atmospheric conditions. Furthermore, over the Iberian Peninsula, increases in ET are proportional to applied irrigation and actually exceed it for almost every simulation month. This is made possible by small but systematic increases in average precipitation over the domain, forming evidence of continental moisture recycling over the Iberian Peninsula. The precipitation increases are of lower magnitude than those of ET and occur much more in lightly irrigated regions than in intensely irrigated regions, confirming that the recycling is partial and mostly nonlocal.

These findings call for an analysis of land-atmosphere coupling processes in the presence of irrigation at the diurnal scale to better describe the impacts on the ABL structure in both irrigated areas and neighbouring regions. In particular, it would be relevant to compare the model to field observations from the LIAISE campaign held in the Ebro Valley in July 2021 (Boone et al., 2025). High resolution modelling experiments using irrigation parametrizations have shown large improvements of performance relative to LIAISE observations for turbulent fluxes, air temperature and humidity (Lunel et al., 2024a; Udina et al., 2024); stressed the importance of the convection parametrization for the response of precipitation to irrigation (Udina et al., 2024); and identified interactions of irrigation-induced heterogeneities with regional breeze circulations (Lunel et al., 2024b). Conducting similar analyses with the simulation setup used in this study should provide insights into the ability of a climate model to reproduce the complex structure of these heterogeneities (Mangan et al., 2023) and their impacts on the ABL and atmospheric water cycle.

Appendix A: Offline parameter tuning of irrigation scheme

The irrigation scheme used in this study was introduced in Arboleda-Obando et al. (2024), who also conducted a calibration and sensitivity analysis of the parameters introduced. The soil moisture target parameter (named β) was found to be the dominant parameter that drives the irrigation withdrawals, and a default value of 0.9 was selected for the global model (meaning the target soil moisture is 90 % of field capacity soil moisture). Although the parametrization uses a single value for the entire simulation domain, Arboleda-Obando et al. (2024) suggested that the value obtained on global average was mainly the result of two main irrigation practices, flooding in rice paddies and other types of irrigation that are less water-intensive:

The results roughly show at least two classes for β, with the first with values of 1.2 and 1.4 (for instance, China and north India) and the second with values of 0.6. Using at least two β values is not enough to reduce the irrigation bias at a global scale, but it has an important effect on the spatial distribution of the irrigation bias in Southern Asia, the region with the largest paddy rice area. These results suggest that the β parameter should have at least two values, namely 1.3 in areas with paddy rice and 0.6 in the rest of the irrigated areas.

Figure A1 Mean seasonal cycle of river discharge at four stations, for the no_irr, irr_default_beta0.9 and irr_reduced_beta0.6 offline simulations and GRDC observations (1983–2010).

Therefore, for this study of the Iberian Peninsula, preliminary offline simulations were used to experiment with different values of the β parameter. Three simulations are presented here:

• without irrigation (no_irr),

• with irrigation using the default value β=0.9 (irr_default_beta0.9),

• with irrigation using the reduced value β=0.6 (irr_reduced_beta0.6)

They were run using the WATCH Forcing Data ERA-Interim (WFDEI Weedon et al., 2014) forcing dataset from 1980 to 2010. Four GRDC discharge stations (see Sect. 2.3 for more details) were selected on the largest rivers of the Iberian Peninsula (Ebro, Tagus, Douro and Guadiana), to represent integrated impacts of irrigation over large river basins. The first 3 years were considered as a spin-up for vegetation and hydrological variables and the simulations were analysed over the remaining period (1983–2010). They revealed very large differences in simulated discharge between no_irr and irr_default_beta0.9 with one mostly overestimating river discharge and the other underestimating it, especially in summer, the season where irrigation is the largest. The irr_reduced_beta0.6 shows an intermediate behaviour with lower discharge than no_irr, but avoiding complete depletion of the river reservoir in summer. These preliminary results in offline simulations justified the choice of a target soil moisture at 60 % of field capacity SM for the coupled ICOLMDZOR LAM simulations.

Appendix B: Three-year spin-up

The simulations started after a 3-year spin-up run to enable the stabilization of the vegetation and hydrological variables, in particular irrigation and ET. The ORCHIDEE model starts with a saturated soil and no vegetation and therefore the first year of the spin-up shows no irrigation on average over the Iberian Peninsula, and a very small mean value of LAI (Fig. B1g–h). After three years, variations of total runoff, ET, total soil moisture in the ORCHIDEE 2-metre soil column, and routing reservoirs are largely driven by the interannual variability of precipitation, as illustrated by the correlation coefficients of Tables B1 and B2. It must be noted that over the final years, the simulations exhibit a decrease in precipitation and therefore in other hydrological variables. This trend in precipitation is also present in ERA5 and mostly dictated by the lateral boundary conditions of the simulation, which is why it was not considered as an indicator of an incomplete spin-up process.

Figure B1 Annual mean values over the Iberian Peninsula in the 3-year spin-up and irr simulation for (a) precipitation, (b) total runoff, (c) evapotranspiration, (d) soil moisture in the 2 m soil column, (e) groundwater reservoir, (f) river reservoir, (g) irrigation, (h) mean Leaf Area Index. Annual values are very similar for the no_irr simulation but irr was used to include irrigation.

Table B1 Correlation matrix (Pearson's r coefficient) for annual values of hydrological and vegetation variables in the irr simulation (2010–2022). The level of significance of the correlation coefficient is indicated by ^* if p value <0.05, ${}^{**}$ if p value <0.01, ${}^{***}$ if p value <0.001.

Table B2 Correlation matrix (Pearson's r coefficient) for annual values of hydrological and vegetation variables in the no_irr simulation (2010–2022). The level of significance of the correlation coefficient is indicated by ^* if p value <0.05, ${}^{**}$ if p value <0.01, ${}^{***}$ if p value <0.001.

Appendix C: River discharge

Figure C1 Time series of river discharge for the irr and no_irr simulations and GRDC observations.

Figure C2 Mean seasonal cycle of river discharge for the irr and no_irr simulations and GRDC observations. A mask is applied to the simulations to filter out months without corresponding observation data.

Appendix D: Seasonal comparison to evaluation products

Figure D1 Seasonal mean biases of precipitation (mm d⁻¹) in the no_irr simulation compared with the GPCC product (2010–2019).

Figure D2 Seasonal mean biases of ET (mm d⁻¹) in the no_irr simulation compared with the GLEAM4 product (2010–2019).

Figure D3 Seasonal mean biases of normalized SSM (no unit) in the no_irr simulation compared with the CCI product (2010–2019).

Appendix E: Influence of irrigation on surface soil moisture

Figure E1 (a) Mean seasonal cycle of surface soil moisture (mm contained in the top 10 cm) in the no_irr and irr simulations (2010–2022). (b) Annual mean relative difference (%) in surface soil moisture between the two simulations (irr − no_irr, 2010–2022).

Appendix F: Influence of snow

Figure F1 Influence of snow in the no_irr simulation (2010–2022). (a) Mean seasonal cycle of rainfall, snowfall, snow melt and total precipitation. (b) Annual mean share of snowfall in total precipitation (%).

Appendix G: Impacts of irrigation on land-atmosphere interactions in other seasons

As Fig. 7 only presented results in summer, equivalent figures are shown for winter, spring and autumn, with a different set of colour bar boundaries.

Figure G1 Winter irrigation and its impacts (DJF, 2010–2022). (a) Irrigation (mm d⁻¹) and 10 m wind (irr simulation). Mean changes in the presence of irrigation (irr − no_irr): (b) evaporative fraction, (c) 2 m temperature (K), (d) atmospheric boundary layer height (m), (e) lifting condensation level (m), (f) precipitation (mm d⁻¹), (g) convective available potential energy (J kg⁻¹), (h) moisture convergence (mm d⁻¹). Hatching indicates areas where the change is not statistically significant.

Figure G2 Spring irrigation and its impacts (MAM, 2010–2022). (a) Irrigation (mm d⁻¹) and 10 m wind (irr simulation). Mean changes in the presence of irrigation (irr − no_irr): (b) evaporative fraction, (c) 2 m temperature (K), (d) atmospheric boundary layer height (m), (e) lifting condensation level (m), (f) precipitation (mm d⁻¹), (g) convective available potential energy (J kg⁻¹), (h) moisture convergence (mm d⁻¹). Hatching indicates areas where the change is not statistically significant.

Figure G3 Autumn irrigation and its impacts (SON, 2010–2022). (a) Irrigation (mm d⁻¹) and 10 m wind (irr simulation). Mean changes in the presence of irrigation (irr − no_irr): (b) evaporative fraction, (c) 2 m temperature (K), (d) atmospheric boundary layer height (m), (e) lifting condensation level (m), (f) precipitation (mm d⁻¹), (g) convective available potential energy (J kg⁻¹), (h) moisture convergence (mm d⁻¹). Hatching indicates areas where the change is not statistically significant.