Watershed models are used worldwide to assist with water and nutrient management under conditions of changing climate, land use, and population. Of these models, the Soil and Water Assessment Tool (SWAT) and SWAT+ are the most widely used, although their performance in groundwater-driven watersheds can sometimes be poor due to a simplistic representation of groundwater processes. The purpose of this paper is to introduce a new physically-based spatially-distributed groundwater flow module called gwflow for the SWAT+ watershed model. The module is embedded in the SWAT+ modeling code and is intended to replace the current SWAT+ aquifer module. The model accounts for recharge from SWAT+ Hydrologic Response Units (HRUs), lateral flow within the aquifer, Evapotranspiration (ET) from shallow groundwater, groundwater pumping, groundwater–surface water interactions through the streambed, and saturation excess flow. Groundwater head and groundwater storage are solved throughout the watershed domain using a water balance equation for each grid cell. The modified SWAT+ modeling code is applied to the Little River Experimental Watershed (LREW) (327 km²) in southern Georgia, USA for demonstration purposes. Using the gwflow module for the LREW increased run-time by 20% compared to the original SWAT+ modeling code. Results from an uncalibrated model are compared against streamflow discharge and groundwater head time series. Although further calibration is required if the LREW model is to be used for scenario analysis, results highlight the capabilities of the new SWAT+ code to simulate both land surface and subsurface hydrological processes and represent the watershed-wide water balance. Using the modified SWAT+ model can provide physically realistic groundwater flow gradients, fluxes, and interactions with streams for modeling studies that assess water supply and conservation practices. This paper also serves as a tutorial on modeling groundwater flow for general watershed modelers.

Parameter sensitivity analysis plays a critical role in efficiently determining main parameters, enhancing the effectiveness of the estimation of parameters and uncertainty quantification in hydrologic modeling. In this paper, we demonstrate an uncertainty and sensitivity analysis technique for the holistic Soil and Water Assessment Tool (SWAT+) model coupled with new gwflow module, spatially distributed, physically based groundwater flow modeling. The main calculated groundwater inflows and outflows include boundary exchange, pumping, saturation excess flow, groundwater–surface water exchange, recharge, groundwater–lake exchange and tile drainage outflow. We present the method for four watersheds located in different areas of the United States for 16 years (2000–2015), emphasizing regions of extensive tile drainage (Winnebago River, Minnesota, Iowa), intensive surface–groundwater interactions (Nanticoke River, Delaware, Maryland), groundwater pumping for irrigation (Cache River, Missouri, Arkansas) and mountain snowmelt (Arkansas Headwaters, Colorado). The main parameters of the coupled SWAT+gwflow model are estimated utilizing the parameter estimation software PEST. The monthly streamflow of holistic SWAT+gwflow is evaluated based on the Nash–Sutcliffe efficiency index (NSE), percentage bias (PBIAS), determination coefficient (R²) and Kling–Gupta efficiency coefficient (KGE), whereas groundwater head is evaluated using mean absolute error (MAE). The Morris method is employed to identify the key parameters influencing hydrological fluxes. Furthermore, the iterative ensemble smoother (iES) is utilized as a technique for uncertainty quantification (UQ) and parameter estimation (PE) and to decrease the computational cost owing to the large number of parameters. Depending on the watershed, key identified selected parameters include aquifer specific yield, aquifer hydraulic conductivity, recharge delay, streambed thickness, streambed hydraulic conductivity, area of groundwater inflow to tile, depth of tiles below ground surface, hydraulic conductivity of the drain perimeter, river depth (for groundwater flow processes), runoff curve number (for surface runoff processes), plant uptake compensation factor, soil evaporation compensation factor (for potential and actual evapotranspiration processes), soil available water capacity and percolation coefficient (for soil water processes). The presence of gwflow parameters permits the recognition of all key parameters in the surface and/or subsurface flow processes, with results substantially differing if the base SWAT+ models are utilized.