Comparing SWAT with SWAT-MODFLOW hydrological simulations when assessing 1 the impacts of groundwater abstractions for irrigation and drinking water 2

ion scenarios simulated by SWAT and SWAT-MODFLOW showed different signals in 744 streamflow change. The simulations by both models indicated that drinking water abstraction caused 745 streamflow depletion and that irrigation abstraction caused a slight total flow increase (but decreased 746 the soil or aquifer water storage, which may influence the hydrology outside the catchment). However, 747 the impact of drinking water abstraction on streamflow depletion by SWAT was minimal and 748 underestimated, and the streamflow increase caused by irrigation abstraction was exaggerated 749 compared with the SWAT-MODFLOW simulation, which produced more realistic results. 750 Overall, the new SWAT-MODFLOW model calibrated by PEST, which included the Drain Package 751 and a new auto-irrigation routine, presented a better hydrological simulation, wider possibilities for 752 groundwater analysis, and more realistic assessments of the impact of groundwater abstractions (for 753 either irrigation or drinking water purposes) on streamflow compared with SWAT. Thus, SWAT754 MODFLOW can be used as a tool for managing water resources in groundwater-affected catchments, 755 taking into account its higher computational demand and more time consumption. 756 https://doi.org/10.5194/hess-2019-232 Preprint. Discussion started: 25 June 2019 c © Author(s) 2019. CC BY 4.0 License.


26
Being able to account for temporal patterns of streamflow, the distribution of groundwater resources, 27 as well as the interactions between surface water and groundwater is imperative for informed water 28 resources management. We hypothesize that, when assessing the impacts of water abstractions on 29 streamflow patterns, the benefits of applying a coupled catchment model relative to a lumped semi- groundwater-related outputs, such as spatial-temporal patterns of water table elevation. In the 43 abstraction scenarios analysis, both models indicated that abstraction for drinking water caused some 44 degree of streamflow depletion, while abstraction for auto-irrigation led to a slight total flow increase 45 (but a decrease of soil or aquifer water storages, which may influence the hydrology outside the 46 catchment). In general, the simulated signals of SWAT-MODFLOW appeared more plausible than 47 1. Introduction 58 The interaction between groundwater and surface water is an important aspect of the water cycle, and 59 the management or use of one often impacts the availability and temporal patterns of the other. 60 Improper management and over-exploitation of these water resource components influence the 61 sustainability of both the water resource itself and also the ecosystems that it supports. Groundwater Garcia, 2016). However, interactions between groundwater and surface water are difficult to observe 67 and measure, and it is, therefore, difficult to determine how much of the reduced streamflow recorded 68 in some rivers is due to abstractions and how much is due to natural weather-induced variability in 69 water table elevation. 70 For quantitative assessment of the impacts of pumping wells on streamflow, a hierarchy of modeling 71 tools has been developed, ranging from analytical models based on simple water balance equations to 72 regional, three-dimensional numerical models, depending on the complexity and available data source . Nevertheless, as they do not simulate many of the physical processes and ignore the real-78 world complexity, they may render unrealistic results. In contrast, numerical, process-based models 79 consider the entire complexity and heterogeneity of river-aquifer systems. Such models can simulate 80 the regional groundwater dynamics as well as the interactions between groundwater and surface water. 81 They are therefore part of local water management applications including estimation of streamflow 82 depletion, although they are generally more time-consuming and costly to set up, calibrate, test, and 83 apply. 84 MODFLOW is a physically-based, fully-distributed, and three-dimensional (3D) finite-difference 85 groundwater model, and it is considered a state-of-the-art international standard for simulating and the national government has endeavored to regulate the abstraction of surface and groundwater to a 163 level preventing negative impacts on in-stream biota. Gradually, direct abstraction from surface waters 164 has been prohibited and groundwater abstraction is regulated to secure a certain minimum flow in all 165 Danish rivers, mainly by moving the abstraction wells away from riverbanks and wetlands and 166 implementing a groundwater abstraction permit authority system. However, there still remains some 167 areas where groundwater exploitation is above the sustainable yield and causes streamflow depletion 168 according to the national water resource model (Henriksen et al., 2008). 169 To better understand how abstraction wells used for drinking water or irrigation may influence nearby 170 streamflow, we applied both SWAT and SWAT-MODFLOW to a lowland catchment in Northern 171 Denmark -the Uggerby River Catchment. We hypothesize that, when assessing impacts of water 172 abstractions on streamflow patterns, the benefits of applying SWAT-MODFLOW relative to SWAT 173 outweigh the costs of additional data requirements and computational resources. We compared the 174 performance of the two models and assessed the simulated signals of streamflow in a range of 175 groundwater abstraction scenarios with real wells and abstraction rates for either drinking water or 176 irrigation with both models. The SWAT-MODFLOW complex used in this study was further 177 developed based on the publically available version (https://swat.tamu.edu/software/swat-modflow/) 178 to enable application of the Drain Package of MODFLOW and to allow auto-irrigation. In addition, an 179 approach based on PEST (Doherty, 2018)  The Uggerby River Catchment lies between latitude 57°17′10"-57°35′25" N and longitude 9°58′47"-187 10°19′55" E. It covers an area of 357 km 2 and is located in the Municipality of Hjørring, which is 188 situated in the northern part of Jutland, Denmark (Fig. 1). The Uggerby River originates from the 189 southern part of Hjørring in Sterup and winds its way through the area of Hjørring, Sindal, Mosbjerg,190 Bindslev, and Uggerby and then discharges into the bay Tannisbugten at the coast of the North Sea.

191
The study area has a typical Atlantic climate, which is temperate with an average annual temperature   208 We used the QSWAT 1.5 interface (George, 2017), which works with the latest SWAT Editor version 209 2012.12.19 and is integrated into a QGIS 2.8.1 interface. The input data for the SWAT model in this 210 study include topography, land use, soil, climate, agricultural management, wells, and wastewater 211 discharge as point sources.

212
The catchment was divided into 19 subbasins ( Fig. 1) based on the 32 m pixel size Digital Elevation 213 Model (DEM), which has been resampled from a 1.6 m LIDAR DEM (Knudsen and Olsen, 2008). For   214   the creation of HRUs, we used the land use map based on the Danish Area Information System (NERI,   215 2000) and the soil map based on a national three-layer soil map with a 250 m grid resolution (Greve 216 et al., 2007), and surface slope type was classified into three classes (<2%, 2-6%, >6%). To reduce the 217 number of HRUs and facilitate the posterior model linkage process, land use for range-grasses and 218 range-brush, which covered only 1.3% and 1.9% of the total catchment area, respectively, were merged 219 into pasture, and water (0.9%) was merged into wetland areas. In order to represent the agricultural 220 management practices in detail, the agricultural area was split into three equally sized types with 221 different five-year crop rotation schedules (Table 1) (Table 1). We do not know the specific tile drain distribution within the entire  We assumed that irrigation only occurs in the HRUs where irrigation pumping wells exist (based on a         The Sequential Uncertainty Fitting Algorithm (SUFI2), which is implemented in the SWAT-CUP 302 software (Abbaspour, 2015), was used to calibrate discharge performance in SWAT. The latest 303 SWAT-CUP version (5.1.6.2) was used. Calibration was performed based on daily discharge records  In the study area, two hydrologically connected monitoring stations are found, located at the outlet of 310 subbasin 13 (station A) and subbasin 18 (station B), respectively (Fig. 1). The two stations represent a 311 small (average discharge 1.95 m 3 s -1 ) and relatively large (average discharge 4.56 m 3 s -1 ) stream in 312 Denmark, and both were used for calibration and validation in this study. Station A is located upstream 313 from station B and its flow therefore has an influence on station B. Thus, the simulated discharge of 314 station A was preliminarily calibrated first (initial range of related parameters are shown in Table 2 were calibrated using the same initial parameter range as for station A (Table 2), while the basin-wide 322 parameter ranges from the final calibration step for station A were used as initial ranges. By this 323 approach, we attempted to make the basin-level parameters representative for both upstream and 324 downstream areas. Afterwards, the water stress threshold was calibrated manually to ensure proper 325 simulation of the annual irrigation amount, which ranges from 80 to 120 mm yr -1 and occurs in the 326 period April to October (Aslyng, 1983). Once the calibration was completed and the parameters were 327 fixed, we validated the model by running one simulation from 1 Jan. 1997 to 31 Dec. 2015 using the 328 first 12-years as a warm-up period.

329
To analyze parameter sensitivity and make the sensitivity analysis comparable with SWAT-330 MODFLOW, an additional iteration with 500 simulations was run for the calibration period. In this 331 iteration, the ranges of basin level parameters and subbasin level parameters for the area upstream folder, can be used to adjust SWAT parameters within the PEST routine.

345
The framework using PEST to calibrate SWAT-MODFLOW was firstly introduced by (Park, 2018). 346 We applied this framework to this study as well but with BEOPEST (Doherty, 2018) instead of PEST 347 as the PEST executable file. Figure Table 3. 385 The observed streamflow used for calibrating SWAT-MODFLOW was identical to that used for 386 calibrating SWAT. Relatively continuous observations of the groundwater table were available at the 387 location of two grid cells, and these were used for calibrating the variation of the groundwater table 388 simulated by SWAT-MODFLOW. Because station A is located upstream from station B and its flow 389 thus has an influence on station B, the weight for deriving the objective function for station A, which 390 represents a small stream, was set to 2, and the weight for station B was set to 1. The weights for the 391 two grid cells were set to 1.    438 The SWAT and SWAT-MODFLOW models both represented well the streamflow hydrographs during 439 the calibration period, while during the validation period, one high peak flow event occurred in the 440 SWAT and SWAT-MODFLOW simulations but not in the observations (Fig. 7). The baseflow was 441 generally reproduced well by both models, but the SWAT-MODFLOW visibly performed better.  (Table 5) suggested "very good" performance of both models during the calibration 446 period based on percent bias (PBIAS). During the validation period, the models performed "good" at 447 station A and "satisfactory" at station B. For NSE values, the performance was "very good" for SWAT-448 MODFLOW calibration at station B, "good" for SWAT-MODFLOW calibration at station A, 449 "satisfactory" for SWAT calibration and SWAT-MODFLOW validation at both stations and SWAT 450 validation at station A, but "unsatisfactory" for SWAT validation at station B. For R 2 values, the 451 performance was "good" for SWAT-MODFLOW calibration, "satisfactory" for SWAT calibration 452 and SWAT-MODFLOW validation, but "unsatisfactory" for SWAT validation. 453 The statistical performances of SWAT-MODFLOW with and without PEST calibration were 454 compared. After calibration by PEST, the summary statistics of SWAT-MODFLOW were improved, 455 especially for the validation period at station B where the performance increased from "unsatisfactory" 456 to "satisfactory" according to NSE values (Table 5). In addition, the weighted residuals between 457 simulation and observation were reduced after calibration by PEST, with the reduced residuals mainly 458 coming from streamflow simulation (Table 6). 459 Table 6. 460 In SWAT, almost all the top 12 sensitive parameters (Fig. 8) were surface process parameters (Table   461 2) except for the groundwater parameter GW_DELAY. In contrast, for SWAT-MODFLOW (Table 3) For the water balance, the evaporation simulated by SWAT-MODFLOW was a little higher (13 mm 474 yr -1 ) than that simulated by SWAT, while the total water yields (total stream flow) simulated by SWAT 475 and SWAT-MODFLOW were almost equal ( Table 7). The water balance components, however,  Table 7. The annual abstrations by drinking water wells or irrigation wells set up in the two models were 486 approximately equivalent (Table 8). In the SWAT simulations, compared with the no-wells scenario 487 (scenario 1), a decrease in the average annual stream flow was observed in scenario 2 (only drinking 488 water wells), while an increase was recorded in scenario 3 (only irrigation wells) and scenario 4 (both 489 drinking water and irrigation wells). In the SWAT-MODFLOW simulations, the average annual 490 streamflow decreased not only in scenario 2, but also in scenario 4, and at subbasin 18 outlet in 491 scenario 3, while a slight increase occurred at subbasin 13 outlet in scenario 3. The decrease in scenario 492 2 simulated by SWAT-MODFLOW was much larger than that by SWAT and also closer to the 493 abstracted amount, and the increase at subbasin 13 outlet in scenario 3 simulated by SWAT-494 MODFLOW was apparently lower than that simulated by SWAT (Table 8). Table 8. 496 In SWAT, the decrease of average annual total flow in scenario 2 was minimal as a result of a tiny 497 decrease in the groundwater return flow (Fig. 11a). In scenario 3 and scenario 4, with unchanged tile 498 flow, all the other flow components rose, especially groundwater and lateral soil discharge. In SWAT-499 MODFLOW, the decrease of average annual total flow in scenario 2 also resulted from a decreased 500 groundwater return flow, but the decrease was much larger than that simulated by SWAT. In scenario  Figure 11. 509 When comparing the temporal patterns of streamflow with the no-wells scenario (scenario 1), we found 510 the daily discharge difference in scenario 2 (only drinking water wells) to be almost always negative 511 (sometimes zero), while in scenario 3 (only irrigation wells) and scenario 4 (both drinking water and 512 irrigation wells) it fluctuated around zero in simulations by both SWAT and SWAT-MODFLOW (Fig.   513   12). Thus, the daily flow in the scenario with drinking water wells was almost always lower than the 514 scenario without drinking water wells, and the daily flow in the scenario with only irrigation wells or 515 the scenario with both irrigation and drinking water wells could be higher or lower than the scenario 516 without wells. The daily discharge difference between scenario 2 and the no-wells scenario simulated 517 by SWAT-MODFLOW was obvious, but the difference using SWAT was minimal. In the comparison   Accordingly, the corresponding summary performance statistics were also better for SWAT-539 MODFLOW. The simulated peak flow on 16 October 2014 by both models was much higher than the 540 observed data (Fig. 7). This discrepancy may be attributed to a high record of precipitation on that day 541 based on a 10 by 10 km grid, which may not be representative for the wider catchment. Additionally, 542 it is also likely that the observed streamflow was underestimated as it is calculated from the Q-h 543 relation, which typically does not adequately cover peak flow events (Poulsen, 2013).

544
In the parameter sensitivity analysis, the surface process parameters of the two models shared the same 545 ranges, while the models had different groundwater modules and parameters. While the SWAT-546 MODFLOW calibration was based on an objective function that took into account not only streamflow 547 but also groundwater heads at the location of two wells, the calibration by PEST mainly improved the 548 streamflow simulation performance (Table 4). According to the parameter sensitivity ranking, the  showed that the impact of drinking water abstractions on streamflow in the SWAT simulation was 572 negligible. In the SWAT-MODFLOW set-up, because the aquifer in the Uggerby River Catchment is 573 connected to and interactive with an area outside of the topographical catchment (Fig. 3), the 574 abstraction from an aquifer located in the Uggerby River Catchment not only impacts the hydrology 575 inside but potentially also outside the catchment. According to the water balance, we expected that the  (Table 8), and a slight total flow increase within the catchment (Fig.   594 11b). The subbasin aquifers in the SWAT set-up are closed and have no interaction with areas outside 595 a subbasin. Meanwhile, the abstracted amount of water from aquifers for irrigation is larger than the 596 amount of returning aquifer recharge from irrigated water, and we would therefore expect decreasing 597 groundwater discharge to streamflow in SWAT simulations. However, the SWAT simulations also 598 showed that irrigation led to enhanced streamflow (Table 8, Fig. 11a), which apparently was even 599 higher than the increase simulated by SWAT-MODFLOW. This supports the point mentioned above 600 that SWAT underestimates the abstraction effect on streamflow depletion. SWAT simulations can, 601 therefore, lead to incorrect assessments of the impacts of groundwater abstractions for irrigation on 602 streamflow, while SWAT-MODFLOW provided more realistic assessments.

603
Upon inspecting the SWAT source code, it appears that the groundwater discharge calculation discharge on the current day is highly related to the groundwater discharge on the previous day, and 609 the increase of the groundwater discharge resulting from each irrigation application could then lead to 610 enhanced groundwater discharge for several days in a row. This may explain why the increased 611 discharge following an irrigation event descended more smoothly in SWAT than in SWAT-612 MODFLOW (Fig. 12).

613
In the SWAT-MODFLOW model, the exchange rate between groundwater and surface water is based 614 on the head difference between the river stage (or drain cell stage) and the head of its surrounding   648 In previous studies, after coupling a calibrated SWAT and calibrated MODFLOW model, the SWAT-  (Table 6).   Fortunately, the groundwater heads of the study area did not change much during the study period ( Fig.   687   9, Fig. 10) and the difference inherently exists between the observed and simulated heads, indicating 688 that the error between the ideal simulated initial heads and the actually used simulated initial heads 689 was small. observed head, but this limitation is anticipated to be minor in our study as the groundwater head did 696 not change much in our simulations and the change mainly followed the variation of recharge with 697 precipitation as its source.

698
The average annual streamflow difference and the regular pattern of daily streamflow difference 699 between the abstraction scenarios and the no-wells scenario were generally explained well, but, 700 surprisingly and unexpectedly, the streamflow difference between the scenario with only drinking 701 water wells and the no-wells scenario on 24 March, 2010, simulated by SWAT-MODFLOW at two 702 stations, were positive, being 1.54 and 0.55 m 3 s -1 , respectively (Fig. 12). The streamflow difference 703 between the scenario with only irrigation wells and the no-wells scenario at station B on the extreme 704 peak flow day (16 October, 2014) simulated by SWAT was -5.2 m 3 s -1 but then became positive next 705 day, which cannot be explained well to our best of knowledge so far. However, we found that the 706 general results of this study were not influenced when modifying the value of these two unexpected 707 points to be expected.

708
Both the SWAT and SWAT-MODFLOW simulations were based on the "best" parameter combination 709 achieved through calibration, which was deemed to be satisfactory for the purpose of this study.

710
However, complex models such as SWAT and SWAT-MODFLOW are subject to non-uniqueness (i.e. 711 more than one parameter combination may yield satisfactory results), so future studies may need to 712 consider the uncertainty due to, for example, parameter uncertainty. The calibration tool SWAT-CUP 713 has already been able to evaluate SWAT parameter uncertainty, whereas the new approach based on 714 PEST to calibrate SWAT-MODFLOW needs to be further explored to adapt for model uncertainty 715 analysis.

716
Our results support our original hypothesis that SWAT-MODFLOW can produce more reliable results 717 in the simulation of the effects of groundwater abstraction for either drinking water or irrigation on 718 streamflow patterns. In addition, SWAT-MODFLOW can produce more outputs than SWAT.

719
However, SWAT-MODFLOW also requires more effort and data to be set up and calibrated, and