Quantifying the streamflow response to groundwater abstractions for irrigation or drinking water at catchment scale using SWAT and SWAT–MODFLOW

Groundwater abstraction can cause a decline in the water table, and thereby affects surface streamflow connected to the aquifer, which may impair the sustainability of both the water resource itself and the ecosystem that it supports. To quantify the streamflow response to groundwater abstractions for either irrigation or drinking water at catchment scale and compared the performance of the widely used semi-distributed hydrological model SWAT and an recently integrated surface–subsurface model SWAT–MODFLOW, we applied both SWAT and SWAT–MODFLOW to a groundwater-dominated catchment in Denmark and tested a range of groundwater abstraction scenarios. To accommodate the study area characteristics, the SWAT–MODFLOW model complex was further developed to enable the Drain package and an auto-irrigation routine to be used. A PEST (parameter estimation by sequential testing)-based approach which enables simultaneous calibration of SWAT and MODFLOW parameters was developed to calibrate SWAT–MODFLOW. Both models demonstrated generally good statistical performance for the temporal pattern of streamflow, with better R2 and NSE (Nash–Sutcliffe efficiency) for SWAT–MODFLOW but slightly better PBIAS (percent bias) for SWAT. Both models indicated that drinking water abstractions caused some degree of streamflow depletion, while abstractions for returned irrigation led to a slight total flow increase, but may influence the hydrology outside the catchment. However, the streamflow decrease caused by drinking water abstractions simulated by SWAT was unrealistically low, and the streamflow increase caused by irrigation abstractions was exaggerated compared with SWAT–MODFLOW. We conclude that the SWAT–MODFLOW model produces much more realistic signals relative to the SWAT model when quantifying the streamflow response to groundwater abstractions for irrigation or drinking water; hence, it has great potential to be a useful tool in the management of water resources in groundwater-dominated catchments. With further development of SWAT–MODFLOW and the PEST-based approach developed for its calibration, this study would broaden the SWAT–MODFLOW application and benefit catchment managers.

Background The interaction between groundwater and surface water is an important aspect of the water cycle, and the management or use of one often impacts the availability and temporal patterns of the other. Improper management and over-exploitation of these water resource components influence the sustainability of both the water resource itself and also the ecosystems that it supports. Groundwater abstraction can cause a decline of the water table, and it thereby directly affects surface water bodies connected to the aquifer [1][2][3]. For rivers in which a considerable portion of the streamflow is base flow, this can have a strong influence on the general flow and deteriorate the function of river ecosystems [4,5]. However, interactions between groundwater and surface water are difficult to observe and measure, and it is, therefore, difficult to determine how much of the reduced streamflow recorded in some rivers is due to abstraction and how much is due to natural weather-induced variability in water table elevation. In addition, the specific field techniques (e.g., flow analysis, permeameter tests, thermal regime tests and tracer tests) used to estimate patterns of groundwater/surface water interaction are typically performed at small spatial scales and over a short time period [6][7][8][9][10]. In this regard, surface-subsurface hydrological models can overcome the above limitations to some extent because of their ability to simulate longterm groundwater-surface water interactions through a holistic approach and also enable scenario analysis (e.g., climate change, groundwater abstraction and land use planning, etc.).
The Soil and Water Assessment Tool (SWAT) [11,12] is a semi-distributed catchment-scale hydrological model that has been widely used to predict the influence of project performance, management practices or climate change on water quantity and quality at different geographical locations and scales [13][14][15][16]. In SWAT, the basin is divided into subbasins through a topographybased delineation, each subbasin containing a tributary of the river. Each subbasin is further divided into Hydrologic Response Units (HRUs), which are unique combinations of land use, soil type, and surface slope. HRUs are modelled as lumped and non-geo-located within each subbasin [17], which makes SWAT computationally efficient for long-term simulation. SWAT has also been used to simulate and quantify groundwater resources [18][19][20] or the effects of drinking water or irrigation pumping on streamflow [21,22]. However, the SWAT model has traditionally emphasized surface processes as the model only includes a relatively simple representation of groundwater dynamics, and its output does not give any spatially explicit information on the groundwater table. In the most recent version of SWAT (v. 670), groundwater is represented by a lumped module in individual subbasins divided into a shallow and a deep aquifer. Both the shallow and the deep aquifer may contribute to streamflow as baseflow through a linear reservoir approximation, ignoring distributed parameters such as hydraulic conductivity and storage coefficients [23]. With this simplified implementation of groundwater dynamics in SWAT, the model can mislead evaluation of groundwater resources or perform rather poorly in catchments where the streamflow is strongly dependent on groundwater discharge [24].
To the best of our knowledge, there are two main approaches for making SWAT perform better in groundwater-dominated catchments. One approach is to modify the SWAT groundwater module code itself. For example, Zhang et al. [25] modified the subroutines in the SWAT source code by converting the shallow aquifer water storage change into water table fluctuation with three groundwater parameters added, namely specific yield, the bottom bed burial depth, and shallow aquifer porosity. The modified SWAT could then simulate both water table fluctuations and water storage of the shallow aquifer in time and space. However, it still applied a lumped, linear reservoir approach to simulate groundwater storage and derive the water table at HRU level, which could give rise to errors as the HRUs are not spatially explicit within a subbasin. Pfannerstill et al. [26] implemented a three-storage concept in the groundwater module by splitting the shallow aquifer into a fast and a slow contributing aquifer. Nguyen and Dietrich [27] replaced the deep aquifer in the original SWAT model with the multicell aquifer model. In both of these studies, the modified SWAT model achieved a better prediction of baseflow than the original SWAT model. However, both models only improved a part of aquifer system simulation, either the shallow aquifer or the deep aquifer. In addition, they maintained the semi-distributed approach.
The other approach for improving the performance of SWAT in groundwater-dominated catchments is to couple SWAT with a physically based, spatially distributed numerical groundwater model, such as MODFLOW (modular finite-difference flow model). A number of studies have applied MODFLOW to assess the impact of groundwater abstraction on surface water resources [3,20,28,29]. However, MODFLOW does not simulate surface processes such as land-atmosphere interactions, agricultural management practices, and surface runoff [30,31]. To obtain spatial-temporal varying recharge rates, MODFLOW is therefore often linked with landsurface models such as the precipitation-runoff modelling system [32,33] and the soil and water assessment tool (SWAT) [34,35].
There are a few studies that have integrated SWAT and MODFLOW code into one model complex [17,23,[36][37][38]. The coupled SWAT-MODFLOW version developed by Bailey et al. [37] has several advantages over others: an efficient HRU-grid cell mapping scheme (including generation of geographically explicit HRUs), the ability to use SWAT and MODFLOW models of different spatial domains, public availability of codes, and a graphical user interface that has been recently developed for its application [39,40]. Recently, the current published SWAT-MODFLOW code (Version 2 on the SWAT website) has been applied to several catchments of varying sizes, for water resources assessment or management, such as in the USA [37,41,42], Canada [43], Denmark [44], Iran [45], Japan [46] and India [47]. It has also been further developed for application in large-scale mixed agro-urban river basins [48] and further coupled with the solute transport model reactive transport in 3D (RT3D) [35]. Within the coupled SWAT-MODFLOW framework, SWAT simulates surface hydrological processes, whereas MODFLOW-NWT (a Newton-Raphson formulation for MODFLOW-2005 [49], which improves the solution of unconfined groundwater-flow problems) simulates groundwater flow processes and all associated sources and sinks on a daily time step. Only the domain covered by both SWAT and MODFLOW was coupled, and the original functionality of MODFLOW or SWAT was retained beyond the common domain. At first, the HRUs of SWAT are disaggregated to make them spatially explicit. Then the HRU-calculated deep percolation from SWAT is passed to the grid cells of MODFLOW as recharge, and MODFLOW-calculated groundwater-surface water interaction fluxes are passed to the stream channels of SWAT. Hence, the model complex accounts for two-way interactions between groundwater and surface waters considering both river losses and gains, and allows full distribution of the groundwater domain, thereby enables a potentially much better representation and thus understanding of the spatial-temporal patterns of groundwater-surface water interactions, which are of key importance to catchment management in groundwater-dominated catchments.
In Denmark, approximately 800 million m 3 of water are abstracted annually and used for irrigation (175-259 million m 3 during 1989-2017) or drinking water [50,51], making the country highly dependent on groundwater. Since the very dry summers in 1975 and 1976 led to drying out of many watercourses around some cities in Denmark, the national government has endeavored to regulate the abstraction of surface and groundwater to a level preventing negative impacts on in-stream biota. Gradually, direct abstraction from surface waters has been prohibited and groundwater abstraction is regulated to secure a certain minimum flow in all Danish rivers, mainly by moving the abstraction wells away from riverbanks and wetlands and implementing a groundwater abstraction permit authority system. However, there still remains some areas where groundwater exploitation is above the sustainable yield and causes streamflow depletion according to the national water resource model [52].
To better understand how abstraction wells used for drinking water or irrigation may influence nearby streamflow, we applied both SWAT and SWAT-MODFLOW to a groundwater-dominated catchment in Northern Denmark-the Uggerby River catchment. Even though SWAT-MODFLOW was designed to improve the SWAT performance, the performance of SWAT and SWAT-MODFLOW, especially when assessing the impacts of groundwater abstractions (for either irrigation or drinking water) on streamflow patterns, has rarely been compared. In this study, we compared the performance of the two models and assessed the simulated signals of streamflow in a range of groundwater abstraction scenarios with real wells and abstraction rates for either drinking water or irrigation with both models. To accommodate the study area, the SWAT-MODFLOW complex used in this study was further developed based on the coupling framework developed by Bailey et al. [37] to enable application of the Drain package of MODFLOW and to allow auto-irrigation. In addition, an approach based on PEST (parameter estimation by sequential testing) [53] was developed to calibrate the coupled SWAT-MODFLOW by adjusting SWAT and MODFLOW parameters simultaneously against the observations of both streamflow and groundwater table.

Study area
The Uggerby River catchment lies between latitude 57°17ʹ 10ʹʹ − 57° 35ʹ 25ʹʹ N and longitude 9° 58ʹ 47ʹʹ− 10° 19′ 55" E. It covers an area of 357 km 2 and is located in the Municipality of Hjørring, which is situated in the northern part of Jutland, Denmark (Fig. 1). The Uggerby River originates from the southern part of Hjørring and discharges into the coast of the North Sea. The study area has a typical Atlantic climate, which is temperate with an average annual temperature around 8 °C, being warmest in August (17 °C average) and coldest in January (0.5 °C average). The average annual precipitation during the study period 2002-2015 was approximately 933 mm with no obvious distinctions among seasons.
The mean catchment elevation is 34.5 m a.s.l. and ranges from 0 to 108 m. Land cover in the catchment is dominated by arable agricultural land, and the other land use types include evergreen forest, pasture, wetland, and urban areas. The soil types are loamy sand, sandy loam, and sand. The main crops grown in the area include winter wheat, winter rape, barley, corn, and grass. Artificial tile drains have been installed in parts of the agricultural land in the catchment, although the precise drainage locations are somewhat uncertain [54]. According to an investigation carried out by Hjørring Municipality in 2009, 101 drinking water pumping wells and 57 irrigation pumping wells placed on pasture and agricultural land were registered within the Uggerby River catchment, and another 256 wells exist outside the catchment but inside Hjørring Municipality (Fig. 1). Generally, irrigation only occurs from April to October. The average annual irrigation amount varies from 80 to 200 mm depending on the types of crop and soil conditions [55].

Model set-up and coupling SWAT model set-up
We used the QSWAT 1.5 interface [56], which works with the latest SWAT Editor version 2012.12.19 and is integrated into a QGIS 2.8.1 interface. The input data for the SWAT model in this study include topography, land use, soil, climate, agricultural management, wells, and wastewater discharge as point sources.
The catchment was divided into 19 subbasins (Fig. 1) based on the 32 m pixel size digital elevation model (DEM), which has been resampled from a 1.6 m LIDAR DEM [57]. For the creation of HRUs, we used the land use map based on the Danish Area Information System [58] and the soil map based on a national three-layer soil map with a 250 m grid resolution [59], and surface slope type was classified into three classes (< 2%, 2-6%, > 6%). To reduce the number of HRUs and facilitate the posterior model linkage process, land use for range-grasses and range-brush, which covered only 1.3% and 1.9% of the total catchment area, respectively, were merged into pasture, and water (0.9%) was merged into wetland areas. In order to represent the agricultural management practices in detail, the agricultural area was split into three equally sized types with different 5-years crop rotation schedules (Table 1) based on the real contour of agricultural field plots and the land use map. Similar to land use, soil types covering a minor part of the catchment (1% or less) were merged into similar soil types. The distribution and proportion of each land use, soil type, and slope band after reclassification are shown in Fig. 2. Based on the combination of land use, soils, and slope, the catchment was discretized into 2620 HRUs.
Climate data used in the model comprised the 10-km grid national daily precipitation data (six stations inside the catchment), 20-km grid daily solar radiation and wind speed data (five stations inside or near the catchment), gauged-level daily maximum and minimum temperatures, and relative humidity data (one station, 27 km from the catchment) during 1997-2015 from the Danish Meteorological Institute [60].
Farm type and manure/mineral fertilizer application of each agricultural rotation as well as dates of sowing, harvesting, and tillage were assigned based on reported statistics for 2005 available from [61] (Table 1). We do not know the specific tile drain distribution within the entire catchment. In general, loamy soils in relatively flat areas are known to be tile drained in Denmark (Olesen [54]. To represent this situation, tile drains were set up in agricultural land with a slope less than 2% and for soil types with a clay content above 8% [62], representing 27% of the agricultural land in the catchment. We assumed that irrigation only occurs in the HRUs where irrigation pumping wells exist (their locations were obtained from a MODFLOW model created by NIRAS A/S). It is difficult to know the exact dates and water amount used for irrigation. Thus, to simulate the irrigation, auto-irrigation management was set up based on heat unit scheduling for the HRUs containing irrigation pumping wells. For the auto-irrigation of crops, the water resource used for pumping was defined as the shallow aquifer, and the soil water content, commonly used as an indicator in actual field irrigation [63], was selected as the water stress identifier with 70 mm as the initial water stress threshold. With the number and location of pumping wells as well as their pumping rates obtained from the Well package in the MODFLOW model, the water abstraction amounts from drinking water wells were added up in each subbasin and set as the water use pumped from the shallow aquifer in SWAT.
The only significant point source of the study area is the discharge from the wastewater treatment plant in Sindal located in subbasin 16. With a few other minor sources aggregated to a total discharge from the wastewater treatment plant, a total of 2768.8 m 3 of water was discharged into the stream per day (data are based on an average for the period 2007-2010).

MODFLOW-NWT model set-up
A steady-state version of the MODFLOW-NWT model has previously been set up for the entire Hjørring Municipality, covering an area of 930 km 2 , in which the Uggerby River catchment is situated (Fig. 1). The model set-up was first established in 2011 and then updated in 2016 by the consultant company NIRAS A/S and Hjørring Water Supply Company, and has been applied for water resources management in the Hjørring Municipality.
In the model set-up, the geology is represented by 5 hydro-stratigraphic layers, discretized into 183,112 grids (376 rows and 487 columns) with a discretization of 100 × 100 m. The first, third and fifth layers are dominated by sand with relatively large hydraulic conductivities, while the second and fourth layers are dominated by clay with lower hydraulic conductivities (Fig. 3) [53]. Twenty-two different hydraulic conductivity values exist in the originally calibrated MODFLOW model (Fig. 3). In order to facilitate the posterior SWAT-MODFLOW calibration, we reclassified and grouped the specific hydraulic conductivities into five groups. The grouping was made for grid cells of similar specific hydraulic conductivities, representing the sedimentary materials of clay, silt, silty sand, mixture of silty sand and clean sand, and clean sand, respectively. Each group was assigned a unique specific hydraulic conductivity, which could be targeted for calibration.
For the SWAT-MODFLOW set-up, we converted the modified calibrated steady-state model into a transient model by assigning values to the specific yield (only for the unconfined layer) and specific storage of each cell according to the type of sedimentary materials of the cell and representative values of storage coefficients. The simulated heads generated by the steady-state model were used as the initial head conditions for the transient model.

SWAT-MODFLOW coupling
SWAT and MODFLOW were combined using the coupling framework developed by [37] and following the procedures described in the instructions available from the SWAT website [65].
For this study, the following changes were made to the original SWAT-MODFLOW code: (1) the grid cells in the Drain package were linked with SWAT subbasins. The Drain package for SWAT-MODFLOW are used for removal of groundwater via subsurface drains and the groundwater from drainage in the watershed that are not included in SWAT's subbasin channel network; and (2) groundwater pumping in agricultural areas or pastures is dictated by irrigation applied to HRUs through SWAT's auto-irrigation routines. For the latter, this is achieved by calculating the daily volume of applied irrigation water (irrigation depth * HRU area) and then extracting this volume from the underlying grid cells using MOD-FLOW's Well package (Fig. 4). In this study, the irrigation pumping source was defined as the third layer. When applying the Drain package of MODFLOW, the original tile drain routine in SWAT was disabled.  [44]. After being coupled with MODFLOW, the overland part of SWAT model remains semi-distributed, while the HRU-calculated percolation from SWAT model is explicitly spatial The steps in the coupling procedure included: (1) disaggregation of HRUs to disaggregated Hydrologic Response Units (DHRUs) through GIS processing to make the model spatially explicit; and (2) creation of six linking text files (HRUs to DHRUs, DHRUs to MODFLOW grids, MODFLOW grids to DHRUs, MODFLOW river cells to SWAT subbasin rivers, MODFLOW drain cells to subbasin rivers, irrigation pumping wells in HRUs to MODFLOW grids) through GIS processing. All related files (MODFLOW input files, original SWAT model files, linkage files) were stored in one working directory for SWAT-MODFLOW execution.

Model calibration SWAT calibration
The Sequential Uncertainty Fitting Algorithm (SUFI2), which is implemented in the SWAT-CUP software [66], was used to calibrate discharge performance in SWAT. The latest SWAT-CUP version (5.1.6.2) was used. Calibration was performed based on daily discharge records from 1 Jan. 2002 to 31 Dec. 2008, with a previous 5-years model warm-up period and using Nash-Sutcliffe efficiency (NSE) as the objective function. Five parameters at basin-wide level and 34 parameters at subbasin level related to streamflow were selected and assigned initial calibration value ranges based on expert judgement and previous SWAT applications in Danish catchments [67,68].
There are two hydrologically connected monitoring stations in the study area, located at the outlet of subbasin 13 (station A) and subbasin 18 (station B), respectively ( Fig. 1). The two stations represent a small (average discharge 1.95 m 3 s −1 ) and relatively large (average discharge 4.56 m 3 s −1 ) stream in Denmark, and both were used for calibration and validation in this study. Station A is located upstream from station B and its flow therefore has an influence on station B. Thus, the simulated discharge of station A was preliminarily calibrated first (initial range of related parameters is shown in Table 2), running 3 iterations with 500 simulations each. After the final iteration for station A, the subbasin level parameters for the area upstream station A were fixed, while the final ranges of the basin-wide parameters were used in the subsequent calibration of station B. As the basinwide parameter values can impact the hydrology of the entire catchment, for the calibration of station B, discharge data from both station A and B were included in the objective function. An additional three iterations with 500 simulations were run, where the subbasin level parameters for the remaining area upstream station B were calibrated using the same initial parameter range as for station A (Table 2), while the basin-wide parameter ranges from the final calibration step for station A were used as initial ranges. By this approach, we attempted to make the basin-level parameters representative for both upstream and downstream areas. Afterwards, the water stress threshold was calibrated manually to ensure proper simulation of the annual irrigation amount, which ranges from 80 to 120 mm years −1 and occurs in the period April to October [55]. Once the calibration was completed and the parameters were fixed, we validated the model by running one simulation from 1 Jan. 1997 to 31 Dec. 2015 using the first 12-years as a warm-up period.
To analyze parameter sensitivity and make the sensitivity analysis comparable with SWAT-MODFLOW, an additional iteration with 500 simulations was run for the calibration period. In this iteration, the ranges of basinlevel parameters and subbasin-level parameters for the area upstream station A were the same as those in the final calibration step for station A, while the ranges of subbasin level parameters for the area upstream station B were identical with the final calibration step for station B. Model accuracy during calibration and validation was evaluated using three performance metrics including coefficient of determination (R 2 ), Nash-Sutcliffe efficiency coefficient (NSE [69]), and percent bias (P BIAS ) based on the best behavioral solution.

SWAT-MODFLOW calibration
After model coupling, the SWAT-MODFLOW was calibrated by adjusting SWAT and MODFLOW parameters simultaneously against the observations of both streamflow and groundwater table through a combination of manual calibration and auto-calibration by the widely used PEST approach [53]. The periods used for model warm-up, calibration, and validation were identical to those used for SWAT. Since SWAT-CUP can only be used to calibrate SWAT parameters, the PEST approach was developed and utilized to adjust SWAT and MODFLOW parameters simultaneously. However, SWAT-MOD-FLOW can also be run through SWAT-CUP, whereby the summary statistics of model performance can be derived and directly compared between SWAT and SWAT-MODFLOW. In addition, model.in and Swat_Edit.exe, which are included in the creation of the SWAT-CUP project folder, were used to adjust SWAT parameters within the PEST routine.
The framework for using PEST to calibrate SWAT-MODFLOW was firstly introduced by [70]. We applied the same framework to this study as well but using the distributed, parallel implementation of PEST, BEOPEST [53] instead of PEST as the PEST-executable file, thereby shortening the calibration time considerably (Fig. 5). Five types of files are required to run PEST: PEST control file, PEST-executable file, model batch file, model input template files, and model output . The extracted simulated data are read by PEST using information from the model output instruction file and then compared against the corresponding observations. Each iteration includes a number of model runs according to the control variable set in the PEST control file to allow adjustment of parameter values. After each iteration, the objective function and a Jacobian matrix are calculated, based on which the PEST will make its decision for the next iteration until one of its stopping criteria, specified in the PEST control file, is met. More detailed information about the optimization process and principles of PEST can be found in [71] and the PEST manual [53]. As shown in Table 3, 26 parameters from SWAT related to surface hydrological processes and 13 parameters from MODFLOW were selected and calibrated through PEST. For SWAT parameters, with the parameters related to tile drains and groundwater excluded, the final calibrated parameter values used in SWAT were applied as the initial values in PEST, and the parameter ranges used in the iteration for SWAT parameter sensitivity analysis were employed as the parameter ranges in PEST. By manually adjusting MODFLOW parameter values to test their impact on model outputs, storage coefficients (SY and SS), horizontal hydraulic conductivity (HK), and two drain conductance (COND) were deemed as the potential sensitive parameters, with the value of HANI (the ratio of hydraulic conductivity along columns to hydraulic conductivity along rows) always being 1 and the values of VKA (the ratio of horizontal to vertical hydraulic conductivity) fixed being 3, 5, or 10 from the original MODFLOW set-up. For MODFLOW parameters, the originally calibrated and modified parameter values in the steady-state MODFLOW version were used as the initial parameter values in PEST, and a small range around the initial values was assigned as the parameter range according to the experience from manual calibration and representative values (derived from [72]).
The observed streamflow used for calibrating SWAT-MODFLOW was identical to that used for calibrating SWAT. Relatively continuous observations of the groundwater table were available at the location of two grid cells, and these were used for calibrating the variation of the groundwater table simulated by SWAT-MODFLOW. Because station A is located upstream from station B and its flow thus has an influence on station B, the weight for deriving the objective function for station A, which represents a small stream, was set to 2, and the weight for station B was set to 1. The weights for the two grid cells were set to 1.
In order to establish template files and facilitate the process of modifying parameter values (HK, SS, SY) in the UPW package while running PEST, the parameter value file (PVAL) and Zone file [73] were first established based on the original UPW package through running a code file in FORTRAN. Ten iterations were specified as the stop criteria in the PEST control file. To shorten the calibration time, 11 BEOPEST slaves were created on three computers with BEOPEST as the PEST-executable file so that 11 simulations could be run simultaneously. A total of 638 simulations were run before the stop criteria was achieved. With the calibrated parameters fixed, the water stress threshold was calibrated manually to ensure proper simulation of the annual irrigation amount (ranging from 80 to 120 mm years −1 , occurring in the period between April to October) and make the simulated average annual irrigation amount in the irrigated HRUs (mm years −1 ) comparative with that in the calibrated SWAT model. Finally, the SWAT-MODFLOW model performance was validated following a procedure equivalent to that used for SWAT.

Groundwater abstraction scenarios
In order to evaluate the impacts of both irrigation and drinking water abstractions on streamflow for streams of difference sizes, four abstraction scenarios were designed and applied to the Uggerby River catchment using both models: (1) the no-wells scenario, where all abstractions are terminated; (2) the irrigation-wellsstop scenario, where all abstractions in irrigation wells are terminated, while abstractions in drinking water wells remain; (3) the drinking-wells-stop scenario, where all abstractions in drinking water wells are terminated, while abstractions in irrigation wells remain; and (4) the baseline scenario, where abstractions in all wells are included, which represents the current level of abstraction. We assumed that the point source discharge to the stream in subbasin 16 would remain the same in all scenarios. Once the scenarios were simulated, their impacts on streamflow were analyzed by assessing the average annual runoff amount, the contribution of water balance components, and the temporal dynamics of streamflow. The simulated signals of SWAT and SWAT-MODFLOW in the abstraction scenarios were then compared.

Steady-state MODFLOW performance
Visualization of the proximity of simulated and observed heads (Fig. 6) was used to evaluate how well the modified calibrated MODFLOW model performed at the steady state and three summary statistics were used as Table 3 Initial values, ranges, and calibrated values of the selected parameters for SWAT-MODFLOW calibration using PEST a Means that the parameter applies to the upstream areas, including subbasins: 4, 5, 7-13, while b applies to downstream areas, including subbasins 1, 3, 6, 14-19. "-" indicates that the corresponding parameters can be found in Table 2 Parameter   indicators for goodness of model fit ( Table 4). The simulated heads and summary statistics have changed little compared with the original calibrated MODFLOW setup. Thus, the modified calibrated MODFLOW model was satisfactory and suitable as a basis for coupling to SWAT in transient mode.

SWAT and SWAT-MODFLOW transient model performance
The SWAT and SWAT-MODFLOW models both represented well the streamflow hydrographs during the calibration period, while during the validation period, one high peak flow event occurred in the SWAT and SWAT-MODFLOW simulations but not in the observations (Fig. 7). The baseflow was generally reproduced well by both models, but the SWAT-MODFLOW visibly performed better. Compared with the recommended evaluation criteria by [74], the statistical performance percent bias (P BIAS ) of both models during the calibration period were "very good" (Table 5). During the validation period, the P BIAS of both models were "good" at station A and "satisfactory" at station B. For NSE values, the performance was "very good" for SWAT-MODFLOW calibration at station B, "good" for SWAT-MODFLOW calibration at station A, "satisfactory" for SWAT calibration and SWAT-MOD-FLOW validation at both stations and SWAT validation at station A, but "unsatisfactory" for SWAT validation at station B. For R 2 values, the performance was "good" for SWAT-MODFLOW calibration, "satisfactory" for SWAT calibration and SWAT-MODFLOW validation, but "unsatisfactory" for SWAT validation. The statistics NSE and R 2 were better for SWAT-MOD-FLOW than SWAT during both calibration and validation periods, however, the P BIAS in validation period for SWAT was slightly better than SWAT-MODFLOW during validation period ( Table 5). The statistical performances of SWAT-MODFLOW with and without PEST calibration were compared. After calibration by PEST, the summary statistics NSE and R 2 for SWAT-MODFLOW performance were improved, especially for the validation period at station B where the performance increased from "unsatisfactory" to "satisfactory" according to NSE values (Table 5). In addition, the weighted residuals between simulation and observation were reduced after calibration by PEST, with the reduced residuals mainly coming from streamflow simulation (Table 6). However, the P bias for SWAT-MODFLOW performance in the validation period became slightly worse (Table 5).
In SWAT, almost all the top 12 sensitive parameters (Fig. 8) were surface process parameters ( Table 2) except for the groundwater parameter GW_DELAY. In contrast, for SWAT-MODFLOW (Table 3), all the top 12 sensitive parameters were groundwater parameters with the exclusion of only one surface process parameter OV_N.
Compared with SWAT, the SWAT-MODFLOW model not only produced output for streamflow but also for the groundwater table in each cell on any given day. There was generally a good agreement between the groundwater head level and dynamics simulated by SWAT-MODFLOW and that recorded at the two observation wells within the catchment, though the seasonal well drawdowns in Well A did not always occur in the observations (Fig. 9).
For the water balance, the evaporation simulated by SWAT-MODFLOW was 2.5% higher (13 mm years −1 ) than that simulated by SWAT, while the total water yield (total stream flow) simulated by SWAT-MOD-FLOW was 1% (4 mm years −1 ) lower than SWAT ( Table 7). The water balance components, however, differed substantially. Compared with SWAT, the surface runoff simulated by SWAT-MODFLOW was 36.4% higher (8 mm years −1 ), while the lateral subsurface flow was 29% lower (25 mm years −1 ). In SWAT-MOD-FLOW, the largest contributor to streamflow was the drain flow simulated by the Drain package (constituting Table 5   70% of the streamflow). Conceptually, however, this can also be viewed as a surface-near groundwater contribution. Hence, when lumping the contribution from drains and groundwater, these are clearly the dominant sources for streamflow in both the SWAT and SWAT-MODFLOW model.

Groundwater abstraction scenarios simulation
The annual abstractions by drinking water wells or irrigation wells set up in the two models were approximately equivalent ( The composite parameter senstivity Parameters b SWAT-MODFLOW More details regarding composite parameter sensitivity can be found in [53]. The "a" indicates that the parameter applies to the upstream areas, including subbasins: 4, 5, 7-13, and "b" indicates that the parameter applies to downstream areas, including subbasins 1, 3, 6, 14-19 2 (only drinking water wells), while an increase was recorded in scenario 3 (only irrigation wells) and scenario 4 (both drinking water and irrigation wells). In the SWAT-MODFLOW simulations, the average annual streamflow decreased not only in scenario 2, but also in scenario 4, and at subbasin 18 (Table 8).
In SWAT, the decrease of average annual total flow in scenario 2 was minimal as a result of a tiny decrease (0.3 mm years −1 ) in the groundwater return flow (Fig. 10a). In scenario 3 and scenario 4, with unchanged tile flow, all the other flow components rose, especially groundwater and lateral soil discharge. In SWAT-MOD-FLOW, the decrease of average annual total flow in scenario 2 also resulted from a decreased groundwater return flow, but the decrease (5.5 mm years −1 ) was much larger than that (0.3 mm years −1 ) simulated by SWAT. In scenario 3, the lateral soil runoff and drain flow increased in SWAT-MODFLOW similar to SWAT, while in scenario 4, reduced drain flow was recorded (Fig. 10b). Compared with the no-wells scenario, the amount of evapotranspiration remained unchanged in scenario 2, whereas it increased by 5 mm years −1 in the scenarios with irrigation wells in both the SWAT and SWAT-MODFLOW simulations. In the scenario with only irrigation, evapotranspiration and total flow increased in both the SWAT and SWAT-MODFLOW simulations, but the increase of soil or aquifer water storage (ΔS) decreased according to the water balance. a b Fig. 9 Hydrograph of daily simulated and observed groundwater heads (m a.s.l) of the two wells located in layer 1 used for calibrating the variation of groundwater heads simulated by SWAT-MODFLOW where relatively continuous observed data are available. Also shown are summary performance statistics When comparing the temporal patterns of streamflow with the no-wells scenario (scenario 1), we found the daily discharge difference in scenario 2 (only drinking water wells) to be almost always negative (sometimes zero), while in scenario 3 (only irrigation wells) and scenario 4 (both drinking water and irrigation wells) daily discharge difference fluctuated around zero in simulations by both SWAT and SWAT-MODFLOW (Fig. 11). Thus, the daily flow in the scenario with drinking water wells was almost always lower than the scenario without drinking water wells, and the daily flow in the scenario with only irrigation wells or the scenario with both irrigation and drinking water wells could be higher or lower than the scenario without wells. The daily discharge difference between scenario 2 and the no-wells scenario simulated by SWAT-MODFLOW was obvious, but using SWAT alone did not predict this difference. In the comparison of scenario 3 with the no-wells scenario, when the discharge difference was positive after an irrigation event, it descended smoothly in the SWAT simulation and more sharply in the SWAT-MOD-FLOW simulations.
In the SWAT-MODFLOW set-up, the water exchange between aquifer and streams occurs between each MODFLOW river/drain cell and its surrounding cells. The newly developed SWAT-MODFLOW model complex can output the daily rate of water exchange between aquifer and streams for each subbasin. When the water exchange is positive, it is indicative of water flow from the aquifer to the stream. The temporal pattern of groundwater discharge was the same as for the stream flow, and the temporal patterns of the differences in groundwater discharge between the abstraction scenarios and the no-wells scenario were similar to the differences in streamflow, except for some peak flow days (Figs. 11, 12), which indicates that the abstractioninduced streamflow change followed the groundwater discharge change.

Performance and parameter sensitivity of SWAT and SWAT-MODFLOW
Both the SWAT and SWAT-MODFLOW model simulated the temporal patterns of streamflow generally well at the two hydrological stations during the calibration and validation periods according to evaluation criteria recommended by Moriasi et al. [74]. However, visually SWAT-MODFLOW performed better, especially during recession curves and low flow periods, suggesting a better simulation of the interaction between surface water and groundwater. The peak flow on 16 October 2014 by both models was much higher than the observed data (Fig. 7). This discrepancy may be attributed to a high record of precipitation on that day based on a 10 by 10 km grid, which may not be representative for the wider catchment. Additionally, it is also likely that the observed streamflow was underestimated as it is calculated from the Q-h relation, which typically does not adequately cover peak flow events [75].
In the parameter sensitivity analysis, the surface process parameters of the two models shared the same ranges, while the models had different groundwater modules and parameters. While the SWAT-MODFLOW calibration was based on an objective function that took into account not only streamflow but also groundwater heads at the location of two wells, the calibration by PEST mainly improved the streamflow simulation performance (Table 4). According to the parameter sensitivity ranking, the parameters regarding groundwater processes in SWAT-MODFLOW played an important role in the streamflow simulation performance, while in SWAT, the impact of groundwater module parameters on streamflow simulation was generally insignificant. This reflects the shortcoming of the concept for SWAT groundwater module, which ignores the variability in distributed parameters such as hydraulic conductivity and storage coefficients, represents groundwater by a lumped module in individual subbasins, and contributes to the stream network as baseflow based on a linear reservoir approximation. With this simplified implementation of groundwater dynamics and water exchange between surface water and groundwater in SWAT, the discharge simulated by SWAT cannot be optimized to the same extent as that simulated by SWAT-MODFLOW.
At the first sight of Fig. 9, Well A poorly represents the dynamic and Well B seems to have a systematic bias. However, compared with the original MODFLOW-NWT set-up, which has a mean absolute error of 2.22 m between observed and simulated heads ( b SWAT -MODFLOW Station B -subbasin 18 outlet Fig. 12 The hydrograph of simulated daily groundwater discharge to the stream network in the no-wells scenario and daily groundwater discharge differences between the abstraction scenarios (only drinking water wells, scenario 2; only irrigation wells, scenario 3; both drinking water and irrigation wells, scenario 4) and the no-wells scenario (scenario 1) in the upstream area of station A (a) and upstream area of station B (b), respectively, during the entire study period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015), based on SWAT-MODFLOW errors between the observed and simulated head shown in Fig. 9 are much smaller. It is worth noting that the original MODFLOW-NWT set-up has been used in the management of water resources in the Hørring Municipality since 2009. The availability of spatial-temporal patterns of the groundwater head in SWAT-MODFLOW could significantly benefit groundwater resources management and provide spatial explicitly water resources dynamics within a catchment. The outputs of SWAT-MODFLOW in this study showed that the model performed well, not only in streamflow simulations but also with respect to the spatial-temporal patterns of the simulated groundwater head. In contrast, since no information of groundwater table output is provided by SWAT, its goodness in streamflow simulation may potentially be based on an improper groundwater simulation where its performance on groundwater simulation is unknown.

Models ability to simulate effects of groundwater abstractions on streamflow
In scenario 2 where only drinking water wells are active according to the water balance where there is no change in evaporation compared with the no-wells scenario, we expected that the streamflow depletion simulated by SWAT would be approximately equivalent to the abstracted water volume, taking into account a possible small change in the aquifer or soil storage. However, results in this study showed that the impact of drinking water abstractions on streamflow in the SWAT simulation was negligible. In the SWAT-MODFLOW set-up, because the aquifer in the Uggerby River catchment is connected to and interactive with an area outside of the topographical catchment (Fig. 1), the abstraction from an aquifer located in the Uggerby River catchment not only impacts the hydrology inside but potentially also outside the catchment. According to the water balance, we expected that the SWAT-MODFLOW simulated streamflow depletion in the catchment would be lower at a level somewhat equivalent to the abstracted water volume. With equivalent abstraction for drinking water, the annual flow decrease simulated by SWAT-MODFLOW was much larger than that by SWAT and closer to the abstracted volume. Therefore, we conclude that SWAT simulations underestimate the impacts of groundwater abstraction for drinking water on streamflow depletion, while SWAT-MODFLOW provided more realistic assessments.
The simulated irrigation operation abstracts water from an aquifer and then applies the water onto the surface of agricultural land or pasture. Most of the water infiltrates back into the soil and is then utilized by the vegetation and partly lost through evapotranspiration or infiltrates deeper to the aquifer, and a small part of the water might flow to streams directly as a small increase in surface runoff. Though the abstraction causes groundwater depletion, the recharge from the irrigated water can partly refill the aquifer and produce groundwater discharge. Since in the SWAT-MODFLOW set-up the aquifer in the Uggerby River catchment was connected and interactive with an outside area, after each event of groundwater abstraction for irrigation, the aquifer storage would be recharged not only from the irrigated land area, but also by the groundwater flowing from the outside area. If the recharge rate is larger than the abstracted water amount, the groundwater discharge to the stream will presumably increase. Hence, the irrigation events also brought about a slight increase of average annual stream flow at the subbasin 13 outlet (Table 8), and a slight total flow increase within the catchment (Fig. 10b). Another possible reason for the streamflow increase is that irrigation transfer water from less pervious to more pervious formations, which are more directly connected to the stream.
The subbasin aquifers in the SWAT set-up are closed and have no interaction with areas outside a subbasin. Meanwhile, the abstracted amount of water from aquifers for irrigation is larger than the amount of returning aquifer recharge from irrigated water, and we would therefore expect a decrease or a lower increase than SWAT-MOD-FLOW in groundwater discharge to streamflow in SWAT simulations. However, the SWAT simulations also showed that irrigation led to enhanced streamflow ( Table 8, Fig. 10a), which apparently was even higher than the increase simulated by SWAT-MODFLOW. This supports the point mentioned above that SWAT underestimates the abstraction effect on streamflow depletion. SWAT simulations can, therefore, lead to incorrect assessments of the impacts of groundwater abstractions for irrigation on streamflow, while SWAT-MODFLOW provided more realistic assessments. " Upon inspecting the SWAT source code, it appears that the groundwater discharge calculation equation used in SWAT does not take into account the impact of water abstraction from shallow aquifers on water table fluctuations. Thus, the groundwater removal by abstractions in the SWAT simulation does not have a direct effect on the groundwater discharge, which may explain the somewhat surprising simulation signals of SWAT. In addition, in the equation, the groundwater discharge on the current day is highly related to the groundwater discharge on the previous day, and the increase of the groundwater discharge resulting from each irrigation application could then lead to enhanced groundwater discharge for several days in a row. This may explain why the increased discharge following an irrigation event descended more smoothly in SWAT than in SWAT-MODFLOW (Fig. 11).
In the SWAT-MODFLOW model, the exchange rate between groundwater and surface water is based on the head difference between the river stage (or drain cell stage) and the head of its surrounding groundwater grid cells. This can reflect the temporally dynamic hydrological processes and also the impacts from all the external stressors (e.g., temporally and spatially varying recharge and groundwater abstractions) on water table fluctuations. Naturally, this should also allow SWAT-MOD-FLOW to provide more realistic assessments of the impacts of groundwater abstractions on streamflow in comparison with SWAT.
While setting up the drinking water abstractions in SWAT, three limitations were identified, also reported in [44]. The first is that SWAT only allows one decimal point for abstraction numerical inputs with a unit of 10 4 m 3 days −1 for each month. This means that pumping rate variations within 1 m cannot be simulated by SWAT and that the accuracy of abstraction dynamics thus cannot be guaranteed. As a result of this limitation, the abstraction amount in SWAT and SWAT-MODFLOW was not completely identical. The second limitation is that the abstraction from deep aquifer did not result in any streamflow change. Therefore, all the abstraction sources had to be defined as the shallow aquifer in SWAT to achieve a signal in streamflow despite that we had at least three wells receiving water from a deep aquifer (the fifth layer according to the MODFLOW-NWT set-up). The last limitation is that the abstraction rates of all wells in each subbasin in SWAT have to be summed up to one input value, thereby ignoring the specific location of wells within individual subbasins.
SWAT-MODFLOW overcomes the limitations in SWAT by exploiting the spatial explicitness of MOD-FLOW where groundwater abstraction can be simulated using the Well package, which allows many decimal points for abstraction inputs as well as user-defined units, pumping rates at potentially daily intervals, and wells located in any vertical layer and any grid cell within a subbasin. In addition to the outputs from SWAT, SWAT-MODFLOW also provides fully distributed groundwater-related outputs such as spatial-temporal patterns of water table elevation, distributed aquifer recharge, and groundwater-surface water exchange rates at a cell level, permitting detailed analysis of groundwater and its interaction with surface water. This may be an important input to groundwater resources management (e.g., groundwater abstraction) and the solving of surface water rights issues. These capabilities demonstrate the advantage of SWAT-MODFLOW over modifying the SWAT groundwater module codes to improve groundwater flow simulation [25][26][27], which remains a semi-distributed way to simulate subsurface hydrologic processes and does not generate detailed groundwater outputs. This point supports the findings about the advantages of SWAT-MODFLOW over SWAT in [44] but using a much more complex set-up.

Performance of SWAT-MODFLOW and SWAT relative to other recent studies
In previous studies, after coupling a calibrated SWAT and calibrated MODFLOW model, the SWAT-MOD-FLOW model complex was applied without further calibration [37,43], with calibration against only streamflow observations [44], with separated calibration for streamflow and groundwater head [17], or with simple manual calibration by graphically comparing the simulated and observed streamflow and groundwater head [46]. However, the coupling of a calibrated SWAT and a calibrated MODFLOW cannot guarantee a proper or sufficiently optimized parameter set for the integrated SWAT-MODFLOW model. Furthermore, because groundwater and surface water interact with each other, calibrating the simulation of one part does not guarantee proper simulation of the other part. Application of a combined calibration approach based on PEST allowed us to calibrate the SWAT-MODFLOW model by adjusting simultaneously SWAT and MODFLOW parameters and using observations of both streamflow and groundwater table when deriving the objective function, though a more ideal calibration involves more data, such as feedback fluxes from the groundwater domain, infiltration fluxes from irrigated areas. The calibration results demonstrated that the summary statistics of the SWAT-MODFLOW performance were generally improved by this approach ( Table 6).
The ability of SWAT and SWAT-MODFLOW to evaluate the impacts of groundwater abstraction on streamflow or groundwater-surface water interactions has been tested in previous studies [17,43,44] for example, also found that the SWAT model showed almost no impact of groundwater abstraction on streamflow depletion. Besides the simple representation of groundwater dynamics, the other cause of this, we believe, is that same as suggested above, that the impact of groundwater water removal by abstractions on water table fluctuations is currently not accounted for in the groundwater discharge calculation in the SWAT source code. Our findings are generally consistent with those of these previous studies, although all of the studies tested the effects of groundwater abstraction only by drinking water without considering irrigation and based on assumed drinking water pumping wells. In addition, in all the previous studies using the SWAT-MODFLOW developed by Bailey et al. [37], the River package in the MODFLOW model was the only package used for simulating groundwater-surface water interaction, ignoring the potential drain flow processes. The SWAT-MODFLOW complex used in our study was further developed to allow application of the Drain package and to allow also an auto-irrigation routine to extract water from groundwater grid cells; in this way the impacts of groundwater abstraction for both drinking water and irrigation could be assessed.

Limitations and future research
Several limitations to this study need to be acknowledged. The simulated head generated by the steady-state model was used as the initial head conditions for the transient model, as also suggested in other studies [76,77]. The ideal simulated initial heads should be calibrated with the observed initial heads. However, we did not have enough observed heads at the beginning of the simulation period (1997), so we used the observed heads covering the period 1996-2010 for calibrating the original steady-state MODFLOW-NWT to obtain the simulated initial heads. Fortunately, the groundwater heads of the study area did not change much during the study period ( Fig. 9) and the difference inherently exists between the observed and simulated heads, indicating that the error between the ideal simulated initial heads and the actually used simulated initial heads was small.
An approach based on PEST was utilized to calibrate streamflow and groundwater table variation simultaneously in our SWAT-MODFLOW simulation, which improved the model performance and enabled parameter sensitivity analysis for the model. However, only two wells with relatively continuous time series of observed groundwater head were available and used to calibrate the groundwater variation. Ideally, calibration would involve more wells with continuous time series of observed head, but this limitation is anticipated to be minor in our study as the groundwater head did not change much in our simulations and the change mainly followed the variation of recharge with precipitation as its source.
The average annual streamflow difference and the regular pattern of daily streamflow difference between the abstraction scenarios and the no-wells scenario were generally explained well, but, surprisingly and unexpectedly, the streamflow difference between the scenario with only drinking water wells and the nowells scenario on 24 March, 2010, simulated by SWAT-MODFLOW at two stations, were positive, being 1.54 and 0.55 m 3 s −1 , respectively (Fig. 11). The streamflow difference between the scenario with only irrigation wells and the no-wells scenario at station B on the extreme peak flow day (16 October, 2014) simulated by SWAT was − 5.2 m 3 s −1 but then became positive next day, which to date we have been unable to explain. However, we found that the general results of this study were not influenced when fixing the value of these two unexpected points.
The groundwater abstraction scenario simulations by both the SWAT and SWAT-MODFLOW were based on the "best" parameter combination achieved through calibration, which was deemed to be satisfactory for the purpose of this study. However, complex models such as SWAT and SWAT-MODFLOW are subject to nonuniqueness (i.e., more than one parameter combination may yield satisfactory results), so future studies may need to consider the uncertainty due to, for example, parameter uncertainty. The calibration tool SWAT-CUP has already been able to evaluate SWAT parameter uncertainty, whereas the new approach based on PEST to calibrate SWAT-MODFLOW needs to be further explored and adapted to enable uncertainty analysis.
It was proved that SWAT-MODFLOW can produce more reliable results in the simulation of the effects of groundwater abstraction for either drinking water or irrigation on streamflow patterns. In addition, SWAT-MODFLOW can produce more outputs than SWAT. However, SWAT-MODFLOW also requires more effort and data to be set up and calibrated, and longer time to run (around 6 h for a 19-years simulation in SWAT-MODFLOW by a desktop with an Intel ® Core ™ Processor i7-6700 CPU and 16 GB installed RAM versus 6 min for a SWAT simulation). Therefore, the balance between scientific accuracy and the computational burden should be defined relative to the study goal when choosing between SWAT and SWAT-MODFLOW in a future study. But clearly, if the purpose of a study is to investigate effects of groundwater abstraction on streams, the efforts should be focused on setting up and applying a fully distributed model in the groundwater domain, such as SWAT-MODFLOW. A graphical user interface has also been developed to couple SWAT and MODFLOW based on the publically available version of the SWAT-MODFLOW complex [40]. Since the SWAT-MODFLOW complex used in this study was newly developed and allowed use of the Drain package and auto-irrigation, a new graphical user interface based on the new SWAT-MODFLOW complex could ensure that a study such as that presented here is repeated with less effort and technical challenges.

Conclusions
Generally both models simulated well the temporal patterns of streamflow at the two hydrological stations during the calibration and validation periods, with better R 2 and NSE for SWAT-MODFLOW but slightly better P BIAS for SWAT during validation period. Both models indicated that drinking water abstractions caused streamflow depletion and that irrigation abstractions caused a slight total flow increase. However, abstraction scenarios simulated by SWAT and SWAT-MODFLOW showed different signals in streamflow change. The impact of drinking water abstractions on streamflow depletion by SWAT was minimal and underestimated, and the total flow increase caused by irrigation abstractions was exaggerated, compared with the SWAT-MODFLOW simulation, which produced more realistic results. Overall, the SWAT-MODFLOW model produced more realistic assessments of the impact of groundwater abstractions (for either irrigation or drinking water purposes) on streamflow compared with SWAT. Thus, SWAT-MODFLOW has great potential to be a useful tool for managing water resources in groundwater-dominated catchments. The further development of SWAT-MOD-FLOW for use of Drain package and auto-irrigation routine and the developed PEST-based calibration approach would broaden the SWAT-MODFLOW application.
Abbreviations SWAT : Soil and water assessment tool; MODFLOW: Modular finite-difference flow model; SWAT-MODFLOW: A model coupling the soil and water assessment tool and the modular finite-difference flow model; PEST: Parameter estimation by sequential testing; NSE: Nash-Sutcliffe efficiency; P bias : Percent bias; DHRU: Disaggregated Hydrologic Response Units (DHRUs); HRUs: Hydrologic Response Units.