Quantifying Baseflow Changes Due to Irrigation Expansion Using SWAT+gwflow

1. Introduction

In hydrology, baseflow refer to the portion of streamflow sustained by groundwater seepage during dry periods [1], plays a critical role in maintaining ecological balance and supporting water-dependent activities. However, groundwater extractions for irrigated agriculture, can deplete aquifers and reduce baseflow, potentially impacting ecosystems and downstream users [2,3,4,5]. While modern irrigation systems improve efficiency, they may paradoxically encourages the expansion of irrigated areas, leading to unintended consequences for water [6,7,8,9,10,11]. Therefore, understanding the interactions between groundwater, irrigation, and baseflow is key for sustainable water management [12].

Climate variability is a major driver of yield fluctuations globally, accounting for around one-third of observed variation [13]. In this context, supplementary irrigation offers farmers a way to stabilize yields and reduce losses in dry years [14,15]. Irrigation schedules are typically adjusted based on crop needs, weather, and soil moisture conditions [16,17], meaning that water use fluctuates annually even when crop types and irrigated areas remains constant [18]. However, water permits are often granted based on fixed annual requirements, rarely accounting for seasonal or inter-annual variability in demand [19,20]. This disconnect limits the effectiveness of regulation, which is often poorly equipped to address the dynamic nature of real-world water use, potentially masking the cumulative impacts of variable irrigation on baseflows.

These challenges are especially pronounced in groundwater-fed irrigation systems, where dynamic feedbacks between surface water, groundwater, and irrigation demands are difficult to quantify [21]. Understanding these dynamics is essential for developing sustainable irrigation strategies [22]. Here, hydrological models provide valuable tools for exploring such interactions and supporting decision-making [23,24].

The SWAT model (Soil and Water Assessment Tool) [25] has been widely used to simulate agricultural catchments, including hydrological processes and water management [26,27,28,29,30]. Enhanced versions of SWAT have aimed to better represent groundwater and surface water exchanges (GW-SW). Among them, SWAT+gwflow [31] provides an improved, two-dimensional simulation of groundwater, capturing vertical (e.g. evapotranspiration, percolation, infiltration, exfiltration, pumping) and horizontal (e.g. later flows) water movements, and has been successfully applied in a wide range of studies [32,33,34,35,36]. Compared to SWAT+standalone [37], which represent vertical flow, SWAT+gwflow offers a better compromise between process realism and usability, making it well suited to the objectives of this work. Although SWAT-MODFLOW [38] enables full 3D simulation and multi-layer aquifer representation[39,40], its high data requirements and computational costs limit its applicability.

This study focuses on the San Antonio catchment, in northern Uruguay, an area dominated by intensive, groundwater-irrigated horticulture and citriculture [41,42]. Its high density of hydrometeorological data and sensitivity to El Niño–Southern Oscillation (ENSO)-driven climate variability [43], make it an ideal setting to investigate irrigation-groundwater-streamflow interactions. Increasing reliance on groundwater during dry years intensifies pressure on water resources [44], raising concerns about long-term impacts on baseflow.

This study aims to evaluate how irrigation expansion affects baseflows under inter-annual weather variability. It hypothesizes that increased summer groundwater extraction reduces baseflow during the irrigation season, with the magnitude of this impact varying spatially across the catchment according to the degree of hydraulic connectivity between the aquifer and the stream. These effects are further influenced by climatic conditions, with more pronounced reductions in baseflow during dry periods. While climate change is not explicitly modeled, climate variability is treated as a key driver of irrigation dynamics and constraints.

2. Materials and Methods

This study establishes a groundwater–surface water modeling framework for the San Antonio catchment, located in northern Uruguay, an area characterized by a humid subtropical climate and intensive agricultural activity. The methodological approach combines field observations, farmer interviews, and global datasets to support the calibration and validation of a SWAT+gwflow model, capable of explicitly representing surface and subsurface hydrological processes.

Section 2.1 describes the study area, the monitoring network, and the available datasets. Section 2.2 details the model setup, including subbasins, stream channels, and the aquifer. The model calibration was conducted in two phases: an initial calibration focused on streamflow using observed discharge data, followed by a second phase aimed to minimize errors in simulated groundwater heads and baseflow using groundwater level observations. Land use and agricultural management practices were characterized through field surveys and interviews to accurately incorporate typical crops, irrigation schedules, and water allocation to set up the model.

Section 2.3 outlines the assumptions made regarding groundwater pumping and irrigation expansion. The potential impacts of expanding irrigation in rainfed citriculture areas were evaluated over a 30-year historical period. Climate variability was accounted by classifying months as wet, normal, or dry based on precipitation percentiles, enabling a detailed analysis of water use dynamics under different climatic scenarios.

All the figures in this work were generated in R. The R scripts are available in the PDF file in the supplementary material.

2.1. Study Area and Dataset

The San Antonio Catchment, covers an area of 225 km², experiences a humid subtropical climate, as classified by the Köppen climate system [45]. The mean annual rainfall is 1,430 mm, with lower monthly precipitation during the winter months (June, July, and August). Average daily temperatures range from 10 to 15 °C in winter and from 20 to 30 °C in summer [41]. The monitoring network, maintained by the “Departamento del Agua - CENUR Litoral Norte” since 2018, includes eight observation wells for groundwater, one hydrometric station with reliable discharge computations for surface water (Figure 1c), and fourteen rain gauges distributed throughout the region.

The land use within the catchment (Figure 1a) is characterized by a significant proportion of open field horticulture (OFCP, 34.7%) and citriculture (ORAN, 20.3%), both of which rely on a combination of rainfed crops and irrigation [46]. Most of the irrigation zones are in the southwest of the catchment, while the northern area holds potential for developing supplementary irrigation, particularly for citrus crops (Figure 1b). In addition, the catchment features a large proportion of grassland (GRAS, 35.6%) and pastures (PAST, 0.6%), both utilized for grazing livestock. Other minor land uses include forestry plantations (EUCA, 1.6%), native forest (FRSE, 4.2%), greenhouse horticulture (GHCP, 0.8%), summer crops (AGRL, 2.1%), and urban areas (URBN, 0.1%). This diverse land use pattern reflects the multifunctional nature of the catchment.

The terrain is relatively flat, characterized by a rolling landscape with occasional small hills (Figure 1c). Groundwater heads at steady state range from 20 to 80 meters (Figure 1d) [44]. The upper catchment is dominated by silty clay-textured Haplic Vertisols (50–60% clay, 28–42% silt, 9–12% sand). These soils have high cation exchange capacity (CEC), elevated base saturation, and near-neutral pH, contributing to their natural fertility. Their fine texture and expansive clays enhance water and nutrient retention but also reduce infiltration under saturated conditions, increasing the risk of waterlogging. In contrast, the lower catchment is mainly composed of Luvic Phaeozems with a humic phase. These soils are acidic, with low base saturation, reduced CEC, and a sandy surface horizon that transitions abruptly to a clay-rich subsoil (6–36% clay, 7–11% silt, 56–87% sand). This profile limits moisture retention. While moderate organic matter may mitigate some constraints, these soils are prone to nutrient leaching and have low water-holding capacity [47].

Cross-slope tillage is predominant in OFCP. The use of residue cover under conventional tillage is common in OFCP and AGRL, while GRAS typically lacks residue cover. Heavy forest cover and shortgrass are typical for EUCA, FRSE, and ORAN. Geology comprises sedimentary deposits and fissured basalt rocks from the Cretaceous-Tertiary periods, which form part of the Salto-Arapey aquifer [42,48]. The GLobal HYdrogeology MaPS (GLHYMPS) of permeability dataset [49] identifies four zones with distinctive aquifer hydraulic conductivities (Figure 1e). Additionally, SoilGrids data [50] indicates an absolute depth to bedrock (aquifer thickness) ranging from 10 to 40 meters (Figure 1f).

2.2. Groundwater – Surface Water Model

SWAT+gwflow [31] is an extension of the SWAT+ model [37] that integrates a computationally efficient, 2D finite-difference groundwater flow model to enhance the simulation of GW-SW. This extension provides a more explicit representation of groundwater flow dynamics at a reasonable computation cost. In the model, surface water processes are represented through Hydrologic Response Units (HRUs), which are polygons defined by the combination of land use, soil type, and slope. HRUs serve as the primary elements for hydrological production functions, such as infiltration, exfiltration, and surface runoff. Their structure remains consistent with SWAT+standalone model. The model's transfer function—governing surface and subsurface routing as well as storage—is defined by channels and an aquifer grid. Each aquifer grid cell is linked to its corresponding HRU polygons through spatial overlap. The proportion of each HRU that falls within a grid cell, relative to its total area, is used to compute GW-SW and to calculate a partial recharge volume from a single HRU to the grid cell. As several HRUs can overlap a grid cell, the total recharge volume from HRUs is estimated by summing all partial recharge contributions. Stream cells are also identified by intersecting the stream network with the gwflow grid. This methodology follows the guidelines outlined in the gwflow tutorial [51]. These features allow for the representation of spatial variability in aquifer recharge (derived from SWAT+) and aquifer outflows (e.g., groundwater pumping or discharge to channels).

The change in groundwater storage (∆S) is calculated as the difference between all inflows and outflows (1).

∆S = rech + swgw – gwet - gwsw - satx - soil - ppag ± latl ± bndr

(1)

Inflows (positive values, indicating water added to an aquifer grid cell) include soil water percolating to groundwater (rech), stream seepage to groundwater (swgw). Outflows (negative values, indicating water removed from an aquifer grid cell) include groundwater lost through evapotranspiration (gwet), discharge to streams (gwsw), saturation excess flow (satx), upward transfer to the soil zone (soil), and pumping for agricultural irrigation (ppag). Bidirectional fluxes (positive or negative values, indicating water added to or removed from an aquifer grid cell) include lateral groundwater flow between cells (latl) and boundary fluxes at the watershed edge (bndr).

The San Antonio catchment is modelled with 411 HRUs, ranging in size from 0.1 to 2,200 hectares; 37 channels, varying in length from 70 to 9,000 meters; 24 subbasins, spanning from 0.22 to 31 km²; and an aquifer grid comprising 161 by 284 cells, each with a resolution of 100 by 100 meters. Precipitation from local rain gauges were pre-processed using inverse distance interpolation. This interpolation represents the average precipitation for each sub-catchment, which serves as the precipitation input for the model. Agricultural management practices, including crop types, planting and harvesting schedules, fertilization, and irrigation strategies, were identified through interviews with local farmers. This information was then incorporated into the model to reflect the dominant land uses across the catchment. Irrigation pumping rates were set based on the water demands of irrigated crops, ensuring that simulated water usage aligns with real-world agricultural requirements. Key aquifer properties, such as aquifer thickness and hydraulic conductivity, were defined using global datasets (Figure 1e and Figure 1f). Boundary conditions were assumed to have constant hydraulic heads, enabling groundwater exchange with adjacent catchments.

The model was calibrated during the period from 1 February 2019 to 1 August 2021. Results were validated for the periods from 8 August 2018 to 1 February 2019 and 1 September 2021 to 5 February 2021. These calibration and validation windows were selected to ensure the occurrence of dry and wet periods, with the calibration window comprising two-thirds of the total period of available data. For that purpose, a two-phase supervised random calibration process was used. This process involves a series of iterations, beginning with a uniform distribution across specified parameter ranges. In each iteration, the model runs 360 times, and each simulation is evaluated using the objective function specific to the corresponding calibration phase. The first phase prioritized achieving the best possible fit for streamflow simulations, assessed using the Kling-Gupta Efficiency (KGE) metric [52], as it has been proven to be a good criterion for model calibration [53]. The model parameters calibrated during this phase are listed in Table 1. For parameters marked as substitutive, values were assigned based on SWAT+ documentation [54], while for multiplicative changes, values were selected to limit the change to ±10% for the Curve Number and ±30% for the saturated hydraulic conductivity and the depth of the soil profile. The ±30% range was chosen to avoid large deviations from the expected values. In addition, Nash-Sutcliffe Efficiency (NSE) [55] and percentage bias (BIAS) were used solely for streamflow model validation in this phase.

The second phase aimed to minimize the normalized root mean square error (nRMSE) for groundwater heads and surface baseflow, further improving the model’s accuracy in simulating subsurface hydrological processes and baseflow dynamics. Baseflow separation is made by Lyne-Hollick filter [56] with the R package grwat [57]. The parameters adjusted during this phase are detailed in Table 2. For the parameters corresponding to the aquifer (gwflow.input), the multiplicative changes were set within a wider range than in Table 1, since the initial values were taken from GLHYMPS and SoilGrids, which may deviate significantly from local conditions. Substitutive values were assigned based on the SWAT+gwflow documentation.[51]. For lum.dtl (water stress thresholds), independent substitutive values to trigger irrigation were set between 0.5 and 1 (see section 2.3 for further details). Instead of using a partitioned period of the time series, six observation wells were used for calibration (gw46, gw56, gw57, gw67, gw75, gw84) and two for validation (gw52, gw83). This approach was chosen to address gaps in the records at certain sites.

Surface water prediction uncertainty was assessed by the streamflow logarithmic residuals (2).

L_res = log(Q_sim) – log(Q_obs)

(2)

where L_res are the logarithmic residuals, Q_obs the observed streamflow and Q_sim the simulated streamflow. Groundwater prediction uncertainty was spatially evaluated using the absolute error (3).

E_gw = H_sim - H_obs

(3)

where E_gw is the groundwater absolute error, H_sim is the simulated groundwater heads and H_obs the observed groundwater head at all locations.

2.3. Water Pumping and Irrigation Expansion Criteria

It is a challenging task to determine the amount of water extracted from groundwater in Uruguay, as the locations of all pumping wells are not well known, and the volume of water drawn from the aquifer by each well remains unknown. To address this issue, land use involving irrigated crops was identified through a field campaign and by referencing declared irrigation wells in the national database. Pumping wells used for other purposes, such as domestic and livestock water supply, were not considered, as the volume of water used for these purposes is assumed to be negligible compared to irrigation water in the catchment. Once the irrigated crops were identified, it was important to determine the amount of water extracted for irrigation. This volume was estimated based on crop water requirements. However, in practice, irrigation is applied in excess or deficit. For this reason, the water stress threshold that triggers irrigation in the model was calibrated (table 2). Water stress threshold in SWAT+ range from 0 (severe stress, no irrigation) to 1 (no stress, all water requirements are fully satisfied). Setting thresholds between 0.5 and 1 allows the model to simulate irrigation when water stress becomes critical (at 0.5) up to conditions where irrigation is triggered immediately when any water deficit occurs (close to 1). At this stage, citrus, greenhouse horticulture, and open-field horticulture were aggregated with different water stress thresholds. This approach represents the average condition for farmers and allows for differentiation by crop types without the over-parameterization that would result from setting individual thresholds.

The effect of irrigation expansion on baseflow was estimated using the calibrated SWAT+gwflow model over a 30-year period (1992–2021). To this end, it was assumed that irrigation expansion would occur only in areas currently dedicated to rainfed citriculture. This assumption was made on the basis that citriculture plays a crucial role in the local economy, providing employment and contributing to national exports, underscoring its potential for supplementary irrigation. As shown in Figure 1a and Figure 1b, irrigation expansion could take place in the northern part of the catchment, increasing the total irrigated area from 6,193 to 8,561 hectares, representing a rise from 30% to 41% of the catchment area. In addition, climate variability was classified based moisture conditions to analyze water use dynamics in greater detail. Wet moisture conditions were defined as the months in which precipitation exceeded the 66th percentile of the annual monthly precipitation distribution, while dry moisture conditions were those with precipitation below the 33rd percentile. Normal moisture conditions felt between these thresholds (33rd–66th percentiles). This classification provides a clearer understanding of how water use varies according to the combination of rainfall and groundwater irrigation needed to meet crop water requirements.

3. Results

3.1. Model Development

Figure 2 presents the calibration and validation results for the streamflow simulation. The box-and-whisker plot (Figure 2a), density function (Figure 2c), and flow duration curve (Figure 2e) demonstrate that the simulated values closely align with the observed data. The streamflow hydrograph shown in Figure 2g is displayed on a logarithmic scale to better highlight and identify the baseflow pattern. Overall, total streamflow is well represented, though the model occasionally underestimates baseflow, with less frequent instances of overestimation (Figure 2g). There are also intervals where baseflow is accurately simulated (e.g., mid-2019). The same types of plots were used for validation (Figure 2b, 2d, 2f) to allow for an easy comparison with calibration (Figure 2a, 2c, 2e). While a slight decline in performance is visually noticeable during the validation period, the model continues to produce acceptable results, with KGE values of 0.72 for calibration and 0.74 for validation, NSE values of 0.59 for calibration and 0.53 for validation, and BIAS values of 15.9 for calibration and 9.8 for validation. No seasonal biases were detected.

Streamflow logarithmic residuals (L_res) follow a piecewise relationship (Figure 3a). When simulated streamflow is greater than 0.2, L_res are normally distributed with a mean of 0.47 and a standard deviation of 0.82. For simulated streamflow values less than 0.2, L_res follow a power-law relationship (L_res=5.93Q_sim^0.85, see supplementary material for details on the segmented regression). The error of the piecewise regression follows a normal distribution with zero mean and a standard deviation of 0.83 (Figure 3b and 3c). This error model was used to determine the 95% prediction uncertainty shown in Figure 2g.

Figure 4 presents the hydrographs of groundwater levels at the observation wells. Inset tables compare statistics of observations and simulations. In general, the simulated groundwater levels exhibit less variability compared to the observed levels. Among the sites, the greatest variability is observed at the observation well gw67 (Figure 1c), which suggests that this well is more strongly connected to the stream, indicating a rapid response to streamflow changes. Additionally, the scatter plot demonstrates a strong correlation between the simulated and observed data, indicating the spatial reliability of the simulations in replicating the observed trends.

3.2. Irrigation Expansion and Aquifer Water Balance

A 30-year period (1992–2021) was simulated to assess long-term water balances under both current conditions and irrigation expansion. Table 3 presents the overall water balance of the aquifer. The supplementary material shows the daily groundwater balance error to verify that the time step is adequate.

The irrigation expansion will withdraw 19.5% more water from the aquifer. Table 4 presents the volume of irrigation water used by land use and the percentage of water allocation relative to the total water used in the catchment. The sector that consumes the most water is OFCP, followed by ORAN, with only a marginal water allocation for GHCP. This pattern aligns with the land use distribution, which follows a similar trend.

At the catchment scale, monthly average groundwater heads over the simulation period are shown in Figure 5. Under actual conditions, the mean groundwater head is 43.07 meters, while under the irrigation expansion scenario, it is 42.78 meters. The differences in mean groundwater heads range from -0.5 to 0 meters, indicating a modest but consistent decline in groundwater levels associated with increased irrigation. The difference is more pronounced during periods when the mean groundwater level is already low, suggesting that irrigation expansion may have greater impacts in drier months.

The amount of water used largely depends on moisture conditions, with dry periods requiring more supplementary irrigation. Figure 6 shows the monthly water volume used under wet (Figure 6a), normal (Figure 6c), and dry moisture conditions (Figure 6e). During wet moisture conditions, water demand is similar between the actual scenario and the irrigation expansion. Under normal moisture conditions, an increase in water demand is observed from January to April. In dry moisture conditions, water demand is significantly higher, extending from December to April (in this work, regarding boxplots, a significant difference refers to a clear change in the order of magnitude of the mean values and the width of the boxes). This variable water demand, driven by precipitation, also has a fluctuating effect on baseflow. Figure 6f shows that baseflow decreases under irrigation expansion for the location shown in Figure 7, with the effect lagging by four months after the irrigation season. However, under wet conditions (Figure 6b), no significant changes were detected, and under normal moisture conditions (Figure 6d), only a slight decrease in baseflow was observed from June to August.

At the local scale (gwflow grid), annual groundwater depletion of up to -1.2 meters is predicted (Figure 6b), with depletion zones closely matching areas of irrigation expansion (Figure 1b) and increasing mean annual pumping rates (Figure 6a). This groundwater depletion directly impacts GW-SW, with the predominant trend being a reduction in baseflows (Figure 6c). However, the spatial pattern of GW-SW indicates that baseflow could also increase in certain areas which are prone to greater contributions from the aquifer.

4. Discussion

4.1. Model Performance

The calibration and validation routines indicate that the model had a good performance for simulating streamflow. Streamflow performs similarly in terms of KGE and BIAS to the simulations obtained with the distributed HBV model, presented in a previous study in the same area [41]. The HBV model does not explicitly account for groundwater irrigation; thus, similarities in performance could be due to parameter compensation for the unaccounted irrigation [58]. In addition, over a relatively short period, with the goodness of fit statistic in the good to satisfactory range [59], indicating that the selected time frame may be sufficient for calibration and validation. This is relevant, as studies in other regions suggest that the optimal length for calibrating the SWAT model is five to seven years [60]. Factors such as climate variability, land use changes, model structure, length of the calibration/validation windows, and the uncertainty of forcing data are key considerations in defining model uncertainties [61,62,63,64]. Despite the limited length of the simulation period, the model successfully captures key hydrological processes. The two-phase calibration procedure focuses step by step on surface water and groundwater. This approach enables simplified calibration, reduces computational costs, and minimizes overfitting [65]. The use of logarithmic transformations in the streamflow hydrograph enables a more detailed assessment of model performance in estimating baseflow, particularly in identifying potential stream-aquifer interactions [66]. Additionally, flow duration curves, probability density functions, and box-and-whisker plots provide complementary insights, allowing for a comprehensive evaluation of how well the model matches the observations. These tools address limitations observed in studies that rely solely on time-series plots or aggregated metrics [67].

The streamflow error model indicates an underestimation of low flows and an overestimation of high flows. L_res follows a power-law relationship for low flows, meaning that the lower the streamflow, the lower the ratio of simulated to observed streamflow. In contrast, L_res associated with high flows are normally distributed around a mean of 0.47, indicating that the ratio of simulated to observed streamflow remains approximately constant. These opposing errors tend to compensate for each other in the overall streamflow prediction, contributing to results with relatively small biases over the simulation period. In addition, the error model shows a good fit and could be used to adjust the simulated streamflow accordingly.

The modelling approach produces groundwater head simulations with a level of uncertainty comparable to that reported in previous studies of the same area, which used the MODFLOW model in a fully 3D configuration, yielding a mean error of 0.67 and a correlation coefficient of 0.94 [44]. These results are consistent with those obtained with SWAT+gwflow, which yielded a mean error of 0.97 and a correlation coefficient of 0.93. This indicates only a slight increase in uncertainty. Given that SWAT+gwflow operates with a 100 m groundwater grid size and a global groundwater dataset, the accuracy remains reasonably [68,69]. As with streamflow, some parameter compensation may also occur in the groundwater simulations, particularly since MODFLOW represents steady-state conditions using a different spatial resolution and set of assumptions.

4.2. Assessing Irrigation Expansion

The model results indicate groundwater depletion, particularly in regions experiencing irrigation expansion and increasing groundwater extraction. These findings align with previous studies documented in other regions where groundwater declines in heavily irrigated agricultural areas [11,70]. The observed reductions in baseflows due to irrigation expansion further support established hydrological principles linking groundwater depletion to surface water reductions, where depletion can occur in two ways: an increased flux from streams to the aquifer and a reduced flux from the aquifer to streams. Numerous studies have emphasized that persistent groundwater withdrawals lead to diminished baseflows, reducing streamflow availability during dry periods [71,72]. However, the model also identifies regions where baseflows exhibit an increasing trend. This result suggests localized groundwater contributions influenced by irrigation return flows [73,74] or higher boundary inflow rates due to inconsistent model boundaries [75].

In regions with highly variable climates, such as Uruguay [43], seasonal fluctuations in precipitation significantly influence irrigation practices since the water requirements are filled with a combination of infiltrated water from precipitation and groundwater extractions for irrigation [18]. This characteristic leads to variability in pumping rates, resulting in differential impacts on baseflows that may affect both water quantity and quality on a seasonal basis. [76]. In wet years, when precipitation is abundant, aquifers experience natural recharge, leading to an overall increase in groundwater storage. As a result, the reliance on irrigation water decreases, which in turn reduces the stress on both groundwater reserves and surface water systems. This dynamic also minimizes potential negative impacts on baseflow, as less groundwater extraction translates to more stable streamflow. This characteristic is not considered in current Uruguayan regulations, which assume a constant monthly water supply for irrigation [77]. Conversely, during dry years, aquifers serve as a crucial water source for irrigation, often leading to localized depletion. The extent of this depletion depends on both the intensity of irrigation demands and the precipitation deficit, with the expanding use of groundwater pumping being the primary driver of groundwater depletion worldwide [78]. In some areas, prolonged dry periods can cause significant drawdowns in groundwater levels, potentially leading to hydrological shifts in nearby water bodies.

The lag between groundwater irrigation and its impact on baseflows arises from the fundamental differences in hydrological processes, velocities, and response times to external forcings such as precipitation variability and/or irrigation timings. Surface water typically responds within hours to days to such inputs, particularly in catchments dependent on surface water irrigation [73]. In contrast, groundwater moves through subsurface pathways at significantly lower velocities (ranging from centimeters to meters per day, depending on the aquifer properties), leading to a delayed response [79]. The impact of groundwater withdrawals on surface water can manifest over timescales from weeks to decades, depending on aquifer properties such as transmissivity, storage capacity, and connection to the stream network [80]. As a result, short-term increases in groundwater pumping may not immediately reduce streamflow, but prolonged pumping can lead to persistent baseflow declines, altering watershed hydrology [81]. A particularly concerning effect is the transition of some stream from perennial to ephemeral flow regimes due to an increase pressure on water resources [82], which can have significant ecological and hydrological consequences [83].

4.3. Model Limitations

Modelling GW-SW dynamics presents several challenges. First, the complex structure of fractured aquifers within the catchment often requires 3D modelling approaches to capture intricate flow patterns [84]. This was simplified with a 2D approach. Second, a major challenge arises from the limited knowledge of water extractions from the aquifer. This challenge was addressed through model calibration, validation and error estimation which helps to evaluate and adjust for these unknowns. Some authors have tackled this issue in larger catchments using satellite data to better estimate irrigation applications and/or soft-calibration techniques [85,86]. The main advantage of satellite data is that large areas can be easily estimated, avoiding the extensive time required for field surveys. However, these solutions are based upon further estimations that introduce uncertainty. Third, interactions with surrounding areas further complicate the modelling process. At this stage, transient boundary conditions governing how water enters, moves through, and exits the groundwater flow domain can introduce uncertainty [87,88]. Fourth, plain-dominated landscape can result in low hydraulic gradients, influencing groundwater movement in unexpected ways. Fifth, the small size of the catchment and the reliability of global datasets at such scales may require additional work to achieve the best fit with local datasets [89]. This could be particularly relevant for the aquifer thickness parameter, as partial audiomagnetotelluric scans have revealed that the bedrock of the aquifer follows a pattern, with the aquifer being thinner in the east than in the west of the catchment [42]. Sixth, the relatively short period for model calibration and validation. Usually, long periods are preferred because they may include a wide range of weather conditions, such as droughts or floods. The occurrence of severe droughts could lead to increased pressure on groundwater resources, which could impact the identification of model parameters [90,91], especially aquifer properties, and the quantification of irrigation. Thus, the model should be used with caution in climate change analyses or long memory studies [92,93]. Seventh, the scenario uncertainty itself. The assumed irrigation expansion is expected to occur in the future, but the potential adoption of new irrigation technologies that improve efficiency and reduce water use was neglected. This omission could impact groundwater extractions and return flows. Additionally, climate change was not considered [94], primarily due to the strong assumption that climate variability is the main driver of irrigation scheduling, having a more significant effect on groundwater depletion and baseflow changes than climate change itself. This is particularly relevant in regions with highly variable climates, where identifying the contribution of each component driving these changes remains both a challenge and an opportunity for future research [61].

4.4. Model Benefits

The model is a practical tool for simulating groundwater-surface water (GW-SW) interactions in the San Antonio catchment with a relatively low computation cost. This study demonstrates that the modelling approach produces comparable performance to that observed in previous studies [44], while providing additional benefits in terms of model flexibility and applicability. A key advantage of the model is its ability to incorporate a wide range of hydrological and management scenarios, making it particularly valuable for assessing climate variability and human-induced changes in water resources. Unlike traditional groundwater models, SWAT+gwflow seamlessly integrates surface water and groundwater dynamics while allowing for the inclusion of irrigation pumping schedules, nutrient transport processes in both surface and subsurface flows [31]. Furthermore, the model supports the evaluation of best management practices (BMPs) for environmental impact assessment, such as buffer zones, changes in land use, and optimized irrigation strategies. By linking hydrological and economic analyses, the model facilitates a joint assessment of irrigation expansion and its environmental and economic impacts [95]. This holistic approach enables stakeholders to evaluate trade-offs between agricultural productivity, water resource sustainability, and economic returns, aiding in the development of policies that promote sustainable water use.

5. Conclusion

This study examined the impact of irrigation expansion on baseflow and groundwater levels, with particular attention to seasonal dynamics. Results confirm that the aquifer and the stream are hydraulically connected, although the strength of this connection varies across different zones of the watershed. Increased groundwater extraction in the summer months led to a delayed effect, reducing baseflow in winter.

These developments have important implications for national water regulations in a country with a highly variable climate, such as Uruguay. It provides valuable information such as spatiotemporal dynamics and the hydraulic connectivity between the aquifer and the stream. Current regulations overlook these interactions due to the lack of integrated modeling approaches and fail to account for seasonal variability in irrigation water demand. Recognizing these dynamics is essential for developing more adaptative and science-based water allocation policies. Doing so would promote irrigation practices, optimize water use without depleting aquifers and remove key barriers for sustainable agriculture intensification.