Quantifying Baseflow Changes Due to Irrigation Expansion Using SWAT+gwflow

1. Introduction

In hydrology, baseflow refers to the portion of streamflow that originates from groundwater seepage into a river or stream [1]. It represents the sustained flow of water in a stream during periods without significant rainfall or snowmelt, serving as the primary source of streamflow during dry weather. Understanding the interactions between baseflow, groundwater, and irrigation is crucial for effective water resource management [2]. These interactions are especially significant in catchments where groundwater is the primary source for irrigation. Excessive extraction of groundwater can deplete aquifers, diminish baseflows, and disrupt the delicate balance between groundwater and baseflow. Ensuring sustainable water use and maintaining long-term water availability requires quantifying these interactions effectively.

Irrigation practices play a critical role in shaping these interactions, with different types of irrigation having varying impacts on water resources [3]. Surface water irrigation, which relies on rivers, lakes, or reservoirs, can alter streamflow dynamics by modifying natural flow regimes, potentially reducing environmental flows needed to sustain aquatic ecosystems [4]. Additionally, to add complexity, some exceptions to these patterns have been identified, where excess irrigation can, in certain cases, contribute to increased baseflows [5]. In contrast, groundwater irrigation depends on aquifer withdrawals, which can lower water tables when extraction rates exceed natural recharge [6]. These impacts are complex to evaluate, as they depend on factors such as recharge rates, soil properties, climate conditions, and irrigation efficiency. For this reason, coupled models that integrate surface water and groundwater dynamics provide a valuable tool for assessing these interactions, allowing for more accurate predictions of water availability and the effects of different irrigation strategies on hydrological processes [7].

The SWAT+gwflow model represents a significant advancement in hydrological modelling. By coupling a physically based, spatially distributed groundwater flow module [8] with the SWAT+ model [9], this tool is designed to simulate the complex exchanges between groundwater - surface water (GW-SW) more effectively. The combined model takes advantage of the widely used SWAT+ to represent surface hydrological processes and agricultural practices [10,11,12,13,14], while adding a two-dimensional approach to simulate groundwater fluxes. The SWAT+gwflow model has demonstrated improved performance in representing baseflows compared to the standalone SWAT+ [15,16,17,18]. This makes it a valuable tool for assessing the impacts of irrigation expansion.

The San Antonio catchment, characterized by intensive horticultural and citrus crops irrigated by pumping wells, serves as an ideal experimental area for the characterization of GW-SW [19,20]. This catchment is particularly relevant due to its dependence on groundwater resources and the availability of a relatively rich hydrometeorological network. Furthermore, the area has potential for further irrigation expansion, which could significantly affect water availability and the balance between groundwater and surface water [20]. Such an expansion would increase the demand for groundwater resources and potentially alter baseflows, requiring careful evaluation to ensure sustainable water management. In this context, the objective of this work is to evaluate the impacts of irrigation expansion on baseflows and to develop a predictive tool for assessing the effects of future agricultural practices on surface and groundwater resources, ensuring sustainable water management in groundwater-dependent catchments.

2. Materials and Methods

2.1. Study area and dataset

The San Antonio Catchment, covering an area of 225 km², is in north-western Uruguay and experiences a humid subtropical climate, as classified by the Köppen climate system [21]. 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 [19]. The monitoring network, maintained by the “Departamento del Agua - CENUR Litoral Norte” since 2018, consists of 8 observation wells for groundwater, 1 hydrometric station with reliable discharge computations for surface water (Figure 1c), and 14 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 supplemental irrigation, depending on the practices and resources available to each farmer [22]. Most of the irrigation zones are in the southwest of the basin, 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 utilised 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). The soils are predominantly Argiudolls and Hapluderts, with textures ranging from silty clay to silty clay loam [23]. Geology comprises sedimentary deposits and fissured basalt rocks from the Cretaceous-Tertiary periods, which form part of the Salto-Arapey aquifer [24,25]. The Global Permeability (GLHYMPS) dataset [26] identifies four zones with distinctive aquifer hydraulic conductivities (Figure 1e). Additionally, SoilGrids data [27] indicates an absolute depth to bedrock (aquifer thickness) ranging from 10 to 40 meters (Figure 1f).

Figure 1. (a) Landuse, (b) irrigation scenarios, (c) hypsometry, stations for observation of surface water (sw) and groundwater (gw), (d) steady-state hydraulic head [20], (e) hydraulic conductivity zones [26], (f) aquifer thickness [27].

Figure 1. (a) Landuse, (b) irrigation scenarios, (c) hypsometry, stations for observation of surface water (sw) and groundwater (gw), (d) steady-state hydraulic head [20], (e) hydraulic conductivity zones [26], (f) aquifer thickness [27].

2.2. Groundwater – Surface water model

SWAT+gwflow [8] is an extension of the SWAT+ model [9] 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. 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 the standalone SWAT+ 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 HRUs polygons. This feature allows for the representation of spatial variability in aquifer recharge (derived from SWAT+) and aquifer outflows (e.g. groundwater pumping or discharge to channels).

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 data from local rain gauges were pre-processed using inverse distance interpolation to calculate the average precipitation for each sub-catchment, ensuring accurate spatial distribution of rainfall inputs. 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 optimized during the period from 1 February 2019 to 1 August 2021. Results were verified for the periods from 8 August 2018 to 1 February 2019 and 1 September 2021 to 5 February 2021. These optimization and verification windows were selected to ensure the occurrence of dry and wet periods, with the optimization window comprising two-thirds of the total period of available data. For that purpose, a two-phase supervised random optimization process is 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 optimization phase. The first phase prioritized achieving the best possible fit for streamflow simulations, assessed using the Kling-Gupta Efficiency (KGE) metric [28]. The model parameters optimized during this phase are listed in Table 1.

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 [29] with the R package grwat [30]. The parameters adjusted during this phase are detailed in Table 2. Instead of using a partitioned period of the time series, six observation wells were used for optimization (gw46, gw56, gw57, gw67, gw75, gw84) and two for verification (gw52, gw83). This approach was chosen to address gaps in the records at certain sites.

2.3. Assumptions for water pumping in irrigation and irrigation expansion

Determining the amount of water extracted from groundwater in Uruguay is a challenging task, 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 could be estimated based on crop water requirements; however, in practice, farmers often lack formal irrigation guidelines. Instead, they rely on empirical knowledge and expert judgement regarding crop needs. As a result, irrigation may sometimes be applied in excess or deficit. For this reason, the water stress threshold that triggers irrigation in the model was optimized (table 2). At this stage, citrus, greenhouse horticulture, and open-field horticulture were modelled with three different water stress thresholds, each representing the average condition for its respective crop type.

The effect of irrigation expansion on baseflow is estimated using the SWAT+gwflow model. To this end, it is assumed that irrigation expansion will occur only in areas currently dedicated to rainfed citriculture. This assumption is made because 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.

3. Results

Figure 2 presents the optimization and verification results for the streamflow simulation. The box-and-whisker plot (Figure 2a), density function (Figure 2c), and flow duration curve (Figure 2d) demonstrate that the simulated values closely align with the observed data. The streamflow hydrograph (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. There are also intervals where baseflow is accurately simulated (e.g., mid-2019). The same types of plots were used for verification (Figure 2b,d,f). While a slight decline in performance is visually noticeable during the verification period, the goodness-of-fit metrics—KGE, Nash-Sutcliffe Efficiency (NSE), and percentage bias (BIAS)—indicate that the model continues to produce acceptable results.

Figure 3 presents the hydrographs of groundwater levels at the observation wells. The data has been derived by extracting the groundwater head elevation from each cell corresponding to an observation well. In general, the simulated groundwater levels exhibit less variability compared to the observed levels. Among the sites, the greatest variability is observed at site gw67, which is the closest observation well to the stream. 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.

The model predicts mean annual groundwater depletion of up to -1.2 meters (Figure 4a), with depletion zones closely matching areas of irrigation expansion (Figure 1b) and increasing mean annual pumping rates (Figure 4c). This groundwater depletion directly impacts baseflows, where the predominant trend is a reduction in baseflows (Figure 4b). However, some channels indicate that baseflow could also increase in certain areas. This finding underscores the complexity of assessing GW-SW at the catchment scale, as trends may vary in different directions. In addition, the model provides mean annual GW–SW outputs, including a gridded representation from the gwflow module. The mean annual GW–SW absolute difference (Figure 4d) exhibits a similar spatial pattern to the mean annual baseflow difference (Figure 4b), which is a coarser, semi-distributed output derived from the model’s channel network via the SWAT+ component rather than the gridded gwflow module. Furthermore, flow duration curves for two channels provide additional insight. Channel 28 (Figure 4f) shows a reduction in lower streamflow due to irrigation expansion, whereas Channel 37 (Figure 4e) experiences an increase. In both cases, high streamflow is less affected, suggesting that irrigation primarily influences low-flow conditions while having a limited impact on peak flows.

4. Discussion

The optimization and verification procedures highlight the robustness of the model in simulating streamflow over a relatively short period, demonstrating that the selected time frame is sufficient for effective optimization and verification. This finding is particularly relevant, as studies in other regions suggest that the optimal length for calibrating the SWAT model is five to seven years [31]. Factors such as climate variability, land use changes, model structure, and the uncertainty of forcing data are key considerations in defining model uncertainties [32,33,34]. Despite the limited length of the simulation period, the model successfully captures key hydrological processes. 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 [35]. Additionally, flow duration curves, probability density functions, and box-and-whisker plots provide complementary insights, allowing for a comprehensive evaluation of the model fit to observations. These tools address limitations observed in studies that rely solely on time-series plots or aggregated metrics [36]. They are particularly valuable in hydrological modelling, as they effectively highlight the alignment between simulated and observed values, as well as the distribution of streamflow magnitudes, including extremes and baseflow conditions.

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 [37]. Second, a major challenge arises from the limited knowledge of water extractions from the aquifer. This uncertainty has been addressed through model optimization, which helps estimate and adjust for these unknowns. In contrast, other authors have tackled this issue using satellite data to better estimate irrigation applications and/or soft-optimization techniques [38,39]. Third, interactions with surrounding areas further complicate the modelling process. At this stage, assuming that constant head boundary conditions govern how water enters, moves through, and exits the groundwater flow domain can be a source of uncertainty. Fourth, flat relief 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 [40]. 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 [25]. Despite these limitations, the results show a similar level of uncertainty in simulated waterheads to that observed in previous studies using 3D MODFLOW [20]. This indicates that the model’s performance is consistent with earlier findings while offering additional flexibility for GW-SW, as the new predictive tool can integrate agricultural practices.

The model results indicate groundwater depletion, particularly in regions experiencing irrigation expansion and increasing groundwater extraction. These findings align with previous studies documenting groundwater declines in heavily irrigated agricultural areas [6,41]. The observed reductions in baseflows due to irrigation expansion further support established hydrological principles linking groundwater depletion to surface water reductions. Numerous studies have emphasized that persistent groundwater withdrawals lead to diminished baseflows, reducing streamflow availability during dry periods [42,43]. However, the model also identifies regions where baseflows exhibit an increasing trend, which is an exception rather than the norm. This result suggests localized groundwater contributions influenced by irrigation return flows [5,44] or higher boundary inflow rates due to inconsistent model boundaries [45].

The gridded GW–SW interaction outputs reinforce the spatial patterns observed in mean annual baseflow differences, highlighting consistency between the distributed gwflow module and the semi-distributed SWAT+ channel network representation. This agreement suggests that the model effectively captures the dominant hydrological processes driving groundwater and surface water exchanges. Studies employing similar modelling frameworks have reported comparable findings, demonstrating that integrated hydrologic models can delineate regional trends but may require further calibration to resolve finer-scale heterogeneities [46].

Overall, these findings underscore the complexity of groundwater–surface water interactions, where depletion trends and baseflow responses vary spatially due to hydrological, geological, and anthropogenic factors. The consistency of the model’s predictions with previous research validates its application for assessing catchment-scale water balance changes. Furthermore, the model's capabilities can be expanded in future work to estimate water quality impacts associated with groundwater depletion and altered baseflows. By integrating nutrient transport and contaminant fate modelling, this approach could enhance understanding of how groundwater extraction influences pollutant concentrations and ecosystem health. However, further field validation and sensitivity analyses are recommended to refine predictions and enhance understanding of localized hydrological responses to groundwater extraction.

5. Conclusions

This research has studied the impact of irrigation expansion on baseflow and groundwater levels using the SWAT+gwflow model, marking its first application in Latin America. The model results indicate a reduction in baseflow and groundwater depletion due to irrigation expansion. The findings suggest a minor impact on both baseflow and groundwater levels from increased irrigation pumping. The benefits of this approach include the identification of spatial changes in baseflow, surface water-groundwater exchanges, and groundwater levels, which could contribute to a deeper analysis and understanding of implications at a sub-catchment scale.

The insights from this study reinforce the well-established hydrological principle that groundwater depletion leads to reductions in baseflows, demonstrating a clear link between irrigation expansion and decreased streamflow. However, exceptions were observed, where localized groundwater contributions resulted in increasing baseflows. These anomalies may be attributed to irrigation return flows or inconsistencies in model boundary conditions leading to higher boundary inflow rates.

Despite challenges associated with fractured aquifers, flat topography, and uncertainties in groundwater extractions, the model successfully captures key hydrological processes and GW–SW interactions. The spatial agreement between model outputs and observed patterns supports its reliability as a tool for assessing GW–SW interactions at the catchment and sub-catchment scales.

Future research should focus on incorporating water quality assessments to evaluate the impact of groundwater depletion on pollutant transport. Expanding the model’s predictive capabilities will contribute to more informed water resource management strategies, ensuring a sustainable balance between agricultural water use and long-term groundwater and surface water availability.