*Article* **A New Physically-Based Spatially-Distributed Groundwater Flow Module for SWAT**+

#### **Ryan T. Bailey 1,\*, Katrin Bieger <sup>2</sup> , Je**ff**rey G. Arnold <sup>3</sup> and David D. Bosch <sup>4</sup>**


Received: 17 September 2020; Accepted: 6 October 2020; Published: 9 October 2020

**Abstract:** 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<sup>2</sup> ) 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.

**Keywords:** SWAT+; groundwater; modeling; groundwater–surface water interactions

### **1. Introduction**

The SWAT (Soil and Water Assessment Tool) watershed model [1] is used frequently worldwide to assess water supply, nutrient dynamics, and sediment dynamics under scenarios of climate change, water management practices, and conservation practices. More recently, the SWAT+ model [2] has been presented as an alternative modeling tool with an emphasis on connectivity between spatial features (hydrologic response units, aquifers, channels, reservoirs, ponds, canals, pumps, etc.) within the watershed. Both models, however, often are limited in groundwater-driven watersheds due to the use of simplified, steady-state flow equations to represent water table fluctuation, groundwater

storage, and groundwater discharge to streams [3–5]. Simulating groundwater head and flow using equations for unsteady flow in a heterogeneous aquifer system is important for estimating groundwater supply, groundwater–surface water interaction location and rates, and nutrient (nitrogen, phosphorus) loading to streams via subsurface flow for such watersheds and is essential for the correct assessment of scenarios.

To provide for more accurate simulation groundwater flow in watershed models, SWAT and SWAT+ have been linked to physically-based, spatially-distributed groundwater flow models, most notably MODFLOW [6], which simulates flow in heterogeneous three-dimensional (3D) aquifer systems. Examples for linked SWAT-MODFLOW models include [5,7–10], with [10] being used recently for many watersheds worldwide, e.g., [11–16]. To the authors' knowledge, only one study [17] has linked SWAT+ with MODFLOW.

Although these linked models have provided improved simulation capabilities for coupled stream-aquifer systems within watersheds, their general use is hampered by three principal limitations:


The first limitation has been overcome largely by the development of the graphical user interface QSWATMOD [26], a QGIS plugin tool that facilitates the linkage of SWAT and MODFLOW models and the running and visualization of SWAT-MODFLOW simulations. However, the second and third limitations are ongoing and unsolved.

The objective of this paper is to present a new physically-based spatially-distributed groundwater flow module for SWAT+ that addresses the three limitations listed above. Called *gwflow*, the module: (1) does not require the use of MODFLOW or other groundwater modeling codes, many of which are commercialized; (2) does not add significant run time to SWAT+ simulations; and (3) consists of only three new code subroutines, in keeping with the modular nature of the SWAT+ code. In addition, all data for the *gwflow* module can be prepared using a GIS (ArcMap or QGIS). Due to the detailed description of the *gwflow* module and the underlying solution algorithm presented in Section 2, this paper also serves as a primer for watershed modelers that desire to better understand groundwater flow modeling and how it relates to watershed hydrologic processes.

The capabilities of the *gwflow* module for SWAT+ are demonstrated through an application to the Little River Experimental Watershed (LREW) (327 km<sup>2</sup> ) in southern Georgia, wherein groundwater flow has been observed to be a significant portion of streamflow. In addition, global data sets for aquifer thickness and aquifer permeability are used to populate necessary data for the module to demonstrate how the model could be applied to data-scarce watersheds. Model results are compared against observed streamflow and groundwater levels at eight stream gages and eight groundwater monitoring sites, respectively, to demonstrate the general correctness of the module in simulating hydrological processes. As the purpose of this paper is to present the new *gwflow* module, detailed calibration is not performed. Hence, further calibration is needed for the LREW model if it is to be used for scenario analysis.

### **2. Materials and Methods**

This section presents an overview of the SWAT+ model, the theoretical foundation for the new *gwflow* module and its implementation into the SWAT+ modeling code, and an application of the new SWAT+ modeling code to the Little River Experimental Watershed.

### *2.1. SWAT*+

SWAT+ [2,27,28] is a restructured version of the SWAT watershed modeling code in which watershed hydrologic entities (hydrologic response units, aquifers, channels, reservoirs, ponds, point sources) are designated as spatial objects that can exchange water in any user-specified system. The watershed is divided into subbasins, which are then each divided into one or more landscape units, which lump hydrographs and route them to any spatial object. Within SWAT+, groundwater flow is simulated in the same manner as in the original SWAT model, with the following assumptions and limitations:


#### *2.2. Groundwater Flow Module (gwflow) for SWAT*+

#### 2.2.1. Overview of the gwflow Module

The *gwflow* module for SWAT+ presented in this section is based on a physically-based, spatially distributed groundwater flow model for unconfined aquifers that are hydraulically connected to streams. The aquifer is discretized laterally into a collection of square grid cells, with groundwater volume and groundwater head calculated for each cell on a daily time step. The module uses a single layer to represent the aquifer from the ground surface to the bedrock and is based on the Dupuit–Forchheimer assumption of horizontal flow in unconfined aquifers, and therefore ignores any vertical gradients in groundwater head. Groundwater stresses include recharge, lateral flow, groundwater evapotranspiration (ET), discharge to streams, stream seepage to the aquifer, pumping, and saturation excess flow. These stresses are listed in Table 1, along with whether each stress is a source (groundwater entering the aquifer) or a sink (groundwater leaving the aquifer) if the stress flux (i.e., volumetric flow rate m<sup>3</sup> /day) depends on groundwater head and/or groundwater storage, the required aquifer and streambed properties to compute the stress flux, and the general method for computation. Recharge is provided to cells by connecting SWAT+ HRUs, and cells that intersect SWAT+ channels can exchange water with these channels. The module is embedded into the SWAT+ code, as described later in this section, and replaces the current aquifer module of SWAT+. Methods for calculating groundwater volume and groundwater head are presented in Section 2.2.4.


**Table 1.** Groundwater stresses included in the *gwflow* module. Specific details regarding the method for computing each stress are presented in Section 2.2.4.

*K* = aquifer hydraulic conductivity (m/day). *S<sup>y</sup>* = aquifer specific yield (-). *Kbed* = streambed hydraulic conductivity (m/day). *dbed* = streambed thickness (m).

#### 2.2.2. Solution Strategy for the *gwflow* Module

This section describes the method to calculate groundwater volume and groundwater head throughout a heterogeneous unconfined aquifer system within a watershed domain. The hydrologic fluxes in a typical watershed domain are shown in Figure 1a: land surface and soil fluxes include rainfall, runoff, lateral flow, ET, and percolation; groundwater fluxes include lateral flow from each direction (*Qnorth*, *Qsouth*, *Qwest*, *Qeast*); flow across the watershed boundary, which is a component of lateral flow; recharge *Qrech*; groundwater ET *Qgwet*; pumping *Qpump*; groundwater discharge to streams *Qgwsw*; stream seepage to the aquifer *Qswgw*; and saturation excess flow *Qsatex*, which occurs typically during runoff events related to a rapid rise in the water table and subsequent runoff to streams [29–31]. For SWAT+, the land surface and soil hydrologic fluxes are simulated by original SWAT+ routines, whereas the groundwater fluxes are simulated by the new *gwflow* module. Figure 1a also shows the saturated thickness *s* of the aquifer (water table to the bedrock) and the head *h* of groundwater, measured from a datum, typically mean sea level. Within an unconfined flow system, the head *h* is generally equivalent to water table elevation.

As with all groundwater flow models (e.g., MODFLOW [6], SUTRA [32], CATHY [33], HydroGeoSphere [34]), a control volume approach based on mass conservation is used to establish a water balance equation for each grid cell in the model domain. These equations are then solved for groundwater volume and corresponding groundwater head for each time step of the simulation. Most models use an implicit approach to solve the system of equations, in which head values for all grid cells are updated concurrently (i.e., each head value depends on the head values of surrounding cells at the same time step), and therefore requires linear algebra algorithms or iteration methods. In contrast, the *gwflow* module uses an explicit approach, in which volume and head values for each grid cell are calculated using information (head values of surrounding grid cells; source and sink fluxes) only from the previous time step. Therefore, each grid cell water balance equation is independent of the equations for surrounding cells and can be solved using a simple loop within a computer program, negating the need for computationally intensive solution algorithms. The explicit solution method was chosen for the *gwflow* module due to (1) simplicity in coding; and (2) facility in explaining groundwater modeling within the context of watershed systems for watershed modelers. Readers interested in the implicit approach for groundwater modeling are referred to [35].

Although explicit methods have been available for modeling groundwater heads for many years (e.g., [36]), they often are not used due to the requirement of relatively small time steps for solution stability [36]. However, for linking with watershed models that often use daily time steps, the required time step is reasonable, as will be shown in the model application in Section 2.3.

**Groundwater** 

Discharge to streams

Saturation excess

generally equivalent to water table elevation.

**Stress Source/Sink** 

for computing each stress are presented in Section 2.2.4.

conductivity (m/day). *dbed* = streambed thickness (m).

saturated thickness *s* of the aquifer (water table to the bedrock) and the head *h* of groundwater, measured from a datum, typically mean sea level. Within an unconfined flow system, the head *h* is

**Table 1.** Groundwater stresses included in the *gwflow* module. Specific details regarding the method

Soil Recharge Source No No - SWAT+ HRUs Lateral Flow Mixed Yes Yes *K* Darcy's Law

Groundwater ET Sink Yes Yes - Linear

Stream seepage Source Yes No *Kbed, dbed* Darcy's Law Pumping Sink No Yes - User Specified

**Dependent on Groundwat er Storage** 

Sink Yes Yes *Kbed, dbed* Darcy's Law

**Required Aquifer/Strea m Properties**  **Method for Computing Flow Rate** 

Equation

Storage equation

**Dependent on Groundwat er Head** 

flow Sink Yes No *Sy*

**Figure 1.** Schematics showing the hydrologic fluxes in a typical watershed stream-aquifer system: (**a**) cross-section of stream-aquifer system showing main hydrologic elements and hydrologic fluxes; (**b**) general control volume from (**a**) showing all groundwater inputs/outputs and the groundwater volume *V*; (**c**) plan view of collection of connected control volumes (i.e., grid cells), showing lateral groundwater flow into a grid cell from the four surrounding grid cells. **Figure 1.** Schematics showing the hydrologic fluxes in a typical watershed stream-aquifer system: (**a**) cross-section of stream-aquifer system showing main hydrologic elements and hydrologic fluxes; (**b**) general control volume from (**a**) showing all groundwater inputs/outputs and the groundwater volume *V*; (**c**) plan view of collection of connected control volumes (i.e., grid cells), showing lateral groundwater flow into a grid cell from the four surrounding grid cells.

As with all groundwater flow models (e.g., MODFLOW [6], SUTRA [32], CATHY [33], HydroGeoSphere [34]), a control volume approach based on mass conservation is used to establish a water balance equation for each grid cell in the model domain. These equations are then solved for The derivation of the generic aquifer water balance equation is now presented. The aquifer domain within the watershed is discretized into individual control volume cells (i.e., 3D cubes), which in the *gwflow* module are required to be square on top, e.g., 100 m × 100 m or 200 m × 200, with the thickness of the cell equal to the thickness of the aquifer (ground surface to the bedrock). An example control volume cell and corresponding groundwater fluxes are shown in Figure 1b, with a plan view of the cells shown in Figure 1c. The collection of cells shown in Figure 1c is referred to as a "grid" and covers the area of the watershed. In Figure 1b, the water table is denoted by a dashed blue line, and the saturated thickness *s* of the cell is the vertical distance from the water table to the bedrock. As the cell is composed of both aquifer material and groundwater, total groundwater volume (m<sup>3</sup> ) in the cell is equal to (∆x·∆y·*s*·porosity). However, as we are concerned only with groundwater that can be added to/removed from storage, and not with that which is retained due to suction forces between aquifer sediment and groundwater, the groundwater volume *V* (m<sup>3</sup> ) of concern is:

$$V = \left(\Delta \mathbf{x} \cdot \Delta \mathbf{y} \cdot \mathbf{s}\right) \cdot \mathbf{S}\_y \tag{1}$$

where *S<sup>y</sup>* is specific yield (m<sup>3</sup> water / m<sup>3</sup> of bulk material), i.e., the volume of groundwater released from storage per unit surface area of aquifer per unit decline in groundwater head (i.e., water table) due to gravity.

A volumetric water balance is performed daily for each cell by tracking all groundwater fluxes into/out of the cell:

$$\frac{\Delta V}{\Delta t} = \sum Q\_{in} - \sum Q\_{out} \tag{2}$$

where *t* is time (day), *Qin* represents fluxes (m<sup>3</sup> /day) that add groundwater to the cell, and *Qout* (m<sup>3</sup> /day) represent fluxes that remove groundwater from the cell. Groundwater inputs and outputs can be further categorized into source fluxes, sink fluxes, and lateral fluxes:

$$\frac{\Delta V}{\Delta t} = \sum Q\_{\text{source}} - \sum Q\_{\text{sink}} \pm Q\_{\text{lateral}} \tag{3}$$

where source fluxes include recharge (*Qrech*) and stream seepage (*Qsw*→*gw*) and sink fluxes include groundwater discharge to streams (*Qgw*→*sw*), groundwater ET (*Qgwet*), pumping (*Qpump*), and saturation excess flow (*Qsatex*), as shown in Figure 1b. Lateral fluxes are gradient-driven fluxes that flow across the four sides of the cell. The grid in the *gwflow* module is required to be oriented along north-south and west-east axes, resulting in lateral fluxes in the north, south, west, and east directions (*Qnorth*, *Qsouth*, *Qwest*, *Qeast*), as shown in Figure 1c. The "±" sign indicates that flow can either enter the cell from the cell in that direction, or leave the cell towards the cell in that direction. Including these individual flux terms into Equation (2) yields:

$$\frac{\Delta V}{\Delta t} = \left(Q\_{\text{rech}} + Q\_{\text{sw}\to\text{gw}}\right) - \left(Q\_{\text{gwet}} + Q\_{\text{gw}\to\text{sw}} + Q\_{\text{stax}} + Q\_{\text{punup}}\right) \pm Q\_{\text{north}} \pm Q\_{\text{sath}} \pm Q\_{\text{usst}} \pm Q\_{\text{asst}} \tag{4}$$

where lateral fluxes are assumed to add groundwater to the cell. The actual direction of lateral fluxes, i.e., into or out of the cell, depends on groundwater head patterns, as described in Section 2.2.3.

Saturated thickness can be included in the water balance equation of (4) by inserting Equation (1):

$$\frac{\Delta\left[\left(\Delta\mathbf{x}\cdot\mathbf{l}\mathbf{s}\circ\mathbf{s}\right)\cdot\mathbf{s}\_{\mathbf{y}}\right]}{\Delta t} = \left(Q\_{\text{refl}} + Q\_{\text{SW}\to\text{gw}}\right) - \left(Q\_{\text{gwet}} + Q\_{\text{gw}\to\text{sw}} + Q\_{\text{satter}} + Q\_{\text{pump}}\right) \pm Q\_{\text{unvth}} \pm Q\_{\text{suntl}} \pm Q\_{\text{resl}} \pm Q\_{\text{loss}} \quad \text{(5)}$$

Since a change in saturated thickness is equal to a change in groundwater head, *h* can be substituted for *s* in Equation (5). Assuming *S<sup>y</sup>* of the aquifer does not change in time yields:

$$\frac{\Delta l}{\Delta l}(\Delta x \Delta y S\_{\mathcal{Y}}) = \left(Q\_{\text{rech}} + Q\_{\text{SW} \to \text{gw}}\right) - \left(Q\_{\text{gwet}} + Q\_{\text{gw} \to \text{sw}} + Q\_{\text{sletx}} + Q\_{\text{punp}}\right) \pm Q\_{\text{north}} \pm Q\_{\text{sourl}} \pm Q\_{\text{bresl}} \pm Q\_{\text{loss}} \quad \text{(6)}$$

Using Equation (3) to simplify notation, the change in head ∆*h* for the grid cell is computed by:

$$
\Delta h = \left(\sum Q\_{\text{source}} - \sum Q\_{\text{sink}} \pm Q\_{\text{lateral}}\right) \left(\frac{\Delta t}{S\_y \Delta x \Delta y}\right) \tag{7}
$$

which, if assessing a change in head between two points in time *n* and *n*+1, *h* at the next time *n*+1 can be computed by:

$$h^{n+1} = h^n + \left(\sum Q\_{\text{source}} - \sum Q\_{\text{sink}} \pm Q\_{\text{lateral}}\right) \left(\frac{\Delta t}{S\_y \Delta \mathbf{x} \Delta y}\right) \tag{8}$$

Equation (8) is solved for each cell in the grid using cell information (flux terms, *h*) from the previous time *n*. Each cell is identified using a row and column index using the *i* and *j* notation shown in Figure 1c, leading to:

$$h\_{i,j}^{\ n+1} = h\_{i,j}^{\ n} + \left(\sum Q\_{\text{source}\_{i,j}} - \sum Q\_{\text{sink}\_{i,j}} \pm Q\_{\text{lateral}\_{i,j}}\right) \left(\frac{\Delta t}{S\_{y\_{i,j}} \Delta x \Delta y}\right) \tag{9}$$

and now including all flux terms for the *gwflow* module:

$$h\_{i,j}^{n+1} = h\_{i,j}^n + \left[ \left( Q\_{\text{left}} + Q\_{\text{sur} \rightarrow \text{gw}} \right) - \left( Q\_{\text{gw}t} + Q\_{\text{gw}\to \text{sw}} + Q\_{\text{soft}} + Q\_{\text{punp}} \right) \pm Q\_{\text{unth}} \pm Q\_{\text{soft}} \pm Q\_{\text{aux}} \right] \left( \frac{\Lambda t}{S\_{\text{N}\_{ij}} \Delta x \Delta y} \right) \tag{10}$$

As the computation of *h* for each cell depends only on information from the previous time *n*, this solution strategy is an *explicit* numerical method. Inside the SWAT+ code (see Section 2.2.4 for an explanation of the code structure), the *gwflow* subroutine loops over the grid cells in the watershed domain, solving for head in each cell using Equation (10). As the head value *h <sup>n</sup>* must be known to solve for the head value at the next time step *h n*+1 , the model user must specify the initial head value for each grid cell so that the head values can be calculated on the first day of the SWAT+ simulation. The calculation of each source and sink flux is described in Section 2.2.3. The requirement for solution stability for this explicit approach has been reported [36] as:

$$\frac{\Delta t \cdot K \cdot s}{\Delta x \cdot \Delta y \cdot S\_y} \le \frac{1}{4} \tag{11}$$

Note that Equation (11) does not account for source and sink flux terms. Therefore, the largest ∆*t* that can be used without causing solution instabilities should be found by trial and error when running the model.

#### 2.2.3. Calculating Groundwater Stress Volumetric Fluxes

This section provides details for calculating the individual flux terms in Equation (10) within the context of the SWAT+ modeling code.

#### Soil Recharge

Recharge to the water table is assumed equivalent to the water percolating from the bottom of the soil profile in each HRU of the SWAT+ model. Similar to SWAT-MODFLOW [10], HRU soil percolation is mapped to grid cell recharge using intersected areas of the HRUs and grid cells. During each daily time step, the depth of recharge from an HRU is multiplied by the area of the HRU to provide a volumetric flow rate in m<sup>3</sup> /day. This volume is then distributed to grid cells that overlap the HRU, multiplying the volume by the fraction of the HRU that overlaps the grid cell. These fractions are calculated during model construction by intersecting the HRU shapefile and the grid cell shapefile within a GIS, and then placed in an input file *gwflow.hrucell*.

#### Lateral Flow

The lateral rate of groundwater flow across the north, south, west, and east sides of each grid cell (see Figure 1c) is calculated using Darcy's Law:

$$Q = A \cdot K \cdot \frac{\Delta h}{\Delta x} \tag{12}$$

where *A* is the cross-sectional area of groundwater flow (m<sup>2</sup> ), *K* is the hydraulic conductivity of the aquifer material (m/day), and ∆h/∆x is the hydraulic gradient, with groundwater flowing from areas of high head to areas of low head. Using the *i* and *j* notation from Figure 1c, lateral flux across the west side of a cell (*i*,*j*) is written as:

$$Q\_{\text{west}} = \left(\frac{\Delta x}{\frac{\Delta x/2}{K\_{i-1,j}} + \frac{\Delta x/2}{K\_{i,j}}}\right) \left[\left(\frac{s\_{i-1,j} + s\_{i,j}}{2}\right) \Delta y\right] \left(\frac{h\_{i-1,j} - h\_{i,j}}{\Delta x}\right) \tag{13}$$

where the first term is the harmonic mean *K* value of the cell (*i*,*j*) and the cell to the west (*i*−1,*j*), often used to simulate flow perpendicular to aquifer units; the second term uses the average saturated thickness *s* (head *h*–bedrock elevation) from the two cells to estimate *A*, and the third term uses the *h* values from the two cells to compute the hydraulic gradient. The third term is written so that higher head in the cell to the west results in a positive gradient, indicating an input of groundwater into the cell (*i*,*j*); conversely, lower head in the cell to the west results in a negative gradient, indicating removal of groundwater from the cell (*i*,*j*). Similar equations are written for the other three sides:

$$\begin{cases} Q\_{\text{cast}} = \left(\frac{\Delta x}{\frac{\Delta x/2}{k\_{i+1,j}} + \frac{\Delta x/2}{k\_{i,j}}}\right) \left[\left(\frac{s\_{i+1,j} + s\_{i,j}}{2}\right) \Delta y\right] \left(\frac{h\_{i+1,j} - h\_{i,j}}{\Delta x}\right) \\\ Q\_{\text{north}} = \left(\frac{\Delta x}{\frac{\Delta x/2}{k\_{i,j}} + \frac{\Delta x/2}{k\_{i,j}}}\right) \left[\left(\frac{s\_{i,j} - 1 + s\_{i,j}}{2}\right) \Delta x\right] \left(\frac{h\_{i,j-1} - h\_{i,j}}{\Delta y}\right) \\\ Q\_{\text{sound}} = \left(\frac{\Delta x}{\frac{\Delta x/2}{k\_{i,j} + 1} + \frac{\Delta x/2}{k\_{i,j}}}\right) \left[\left(\frac{s\_{i,j+1} + s\_{i,j}}{2}\right) \Delta x\right] \left(\frac{h\_{i,j+1} - h\_{i,j}}{\Delta y}\right) \end{cases} \tag{14}$$

Note that for homogeneous aquifer systems, the harmonic mean is simplified to the arithmetic mean of hydraulic conductivity.

#### Groundwater ET

*QSW*

package:

Pumping

ET from the saturated zone is simulated only if the water table in a cell is above a specified elevation *zbot* (m), calculated by subtracting a specified ET extinction depth *EXDP* (m) (i.e., the depth below which ET cannot occur) from the ground surface *zsurf* (Figure 2a). The maximum depth of ET that can be removed from the saturated zone is equal to the unsatisfied ET *ETremain* (mm), set equal to the difference between the potential ET (mm) and the actual ET (mm) simulated for each HRU for the day. The connection between HRUs and grid cells for mapping unsatisfied ET is the same as for mapping soil percolation to recharge. The depth of groundwater *ETGW* (mm) removed from the cell is calculated using the linear relationship from MODFLOW's EVT package [6]:

$$\begin{aligned} \mathbf{h}\_{\mathbf{i},j} &< z\_{\text{bot}} \to ET\_{\text{GW}} = \mathbf{0} \\ \mathbf{h}\_{\mathbf{i},j} &> z\_{\text{bot}} \to ET\_{\text{GW}} = ET\_{\text{remain}} \cdot \left(\frac{\mathbf{h}\_{\mathbf{i},j} - z\_{\text{bot}}}{z\_{\text{surf}} - z\_{\text{bot}}}\right) \end{aligned} \tag{15}$$

*Hydrology* **2020**, *6*, x FOR PEER REVIEW 9 of 26

**Figure 2.** Variable depictions and cross-section schematics for (**a**) groundwater Evapotranspiration (ET) calculations and (**b**) groundwater–surface water exchange. **Figure 2.** Variable depictions and cross-section schematics for (**a**) groundwater Evapotranspiration (ET) calculations and (**b**) groundwater–surface water exchange.

Discharge to Streams and Stream Seepage to the Aquifer For grid cells that contain streams, Darcy's Law is used to estimate flow rates between the aquifer and the stream. The direction of flow depends on the relationship between the groundwater head *h* and the stream stage *hstream*. Figure 2b shows the scenario of *h* > *hstream*, resulting in groundwater flow from the aquifer to the stream channel through the streambed with thickness *dbed* (m). Flow occurs along the entire length *L* (m) of the stream within the cell and along the entire bottom width *W* (m) of the channel. *QGWSW* (m3/day) for this scenario is estimated by: The depth of *ETGW* is multiplied by the horizontal spatial area of the HRU to provide a volumetric flow rate in m<sup>3</sup> /day, and then divided amongst the cells connected to the HRU, resulting in a *Qgwet* value for each cell. This groundwater volume is removed from the cell (see Figure 1b). Figure 2a shows the scenario of Equation (15) within the context of the variables. For the row crop (corn as an example), the condition on the left results in no groundwater ET, whereas the condition on the right would result in *ETGW* equal to approximately half of *ETremain*. For the tree, the condition would result in *ETGW* equal to more than half of *ETremain*. Groundwater ET can be significant for areas with high water tables and deep-reaching vegetation roots, e.g., in riparian areas of streams [37,38].

where *Kbed* (m/day) is the hydraulic conductivity of the streambed material and (*WL*) (m2) is the crosssection area of flow. If groundwater head *h* < *hstream*, then stream water recharges the aquifer and

*GW* is calculated. If *h* is greater than the streambed elevation *zbed*, the gradient uses the head difference between *hstream* and *h*; if *h* is lower than *zbed*, the gradient uses the channel depth. These scenarios are summarized in the following set of equations, similar to those in MODFLOW's River

0

*GW SW bed stream bed*

<sup>−</sup> <sup>=</sup> <sup>&</sup>gt;

*SW GW*

→

 = =

*Q*

*Q*

*Q*

*h h Q K WL h h <sup>d</sup>*

→

*GW SW*

→

→

*GW SW*

=

→

→

*bed stream stream SW GW bed*

*z hh h h Q K WL <sup>d</sup>*

*h z h z Q K WL <sup>d</sup>*

0

< < <sup>−</sup>

0

*bed stream bed SW GW bed*

<sup>&</sup>lt; <sup>−</sup>

*h h Q K WL <sup>d</sup>* <sup>→</sup>

*GW SW bed*

( ) *stream*

*bed*

( )

*stream*

*bed*

*bed*

( )

<sup>=</sup>

( )

<sup>=</sup>

(16)

(17)

Discharge to Streams and Stream Seepage to the Aquifer

For grid cells that contain streams, Darcy's Law is used to estimate flow rates between the aquifer and the stream. The direction of flow depends on the relationship between the groundwater head *h* and the stream stage *hstream*. Figure 2b shows the scenario of *h* > *hstream*, resulting in groundwater flow from the aquifer to the stream channel through the streambed with thickness *dbed* (m). Flow occurs along the entire length *L* (m) of the stream within the cell and along the entire bottom width *W* (m) of the channel. *<sup>Q</sup>GW*→*SW* (m<sup>3</sup> /day) for this scenario is estimated by:

$$Q\_{GW \to SW} = K\_{bed}(WL) \left(\frac{h - h\_{stream}}{d\_{bed}}\right) \tag{16}$$

where *Kbed* (m/day) is the hydraulic conductivity of the streambed material and (*WL*) (m<sup>2</sup> ) is the cross-section area of flow. If groundwater head *h* < *hstream*, then stream water recharges the aquifer and *QSW*→*GW* is calculated. If *h* is greater than the streambed elevation *zbed*, the gradient uses the head difference between *hstream* and *h*; if *h* is lower than *zbed*, the gradient uses the channel depth. These scenarios are summarized in the following set of equations, similar to those in MODFLOW's River package:

$$\begin{cases} h > h\_{\text{stream}} \\ z\_{\text{bed}} < h < h\_{\text{stream}} \\ h < z\_{\text{bed}} \end{cases} \quad \begin{cases} Q\_{\text{GW}\to\text{SW}} = K\_{\text{bed}}(\text{WL}) \Big( \frac{h - h\_{\text{stream}}}{d\_{\text{bed}}} \Big) \\ Q\_{\text{SW}\to\text{GW}} = 0 \\ Q\_{\text{GW}\to\text{SW}} = 0 \\ Q\_{\text{SW}\to\text{GW}} = K\_{\text{bed}}(\text{WL}) \Big( \frac{h\_{\text{stream}} - h}{d\_{\text{bed}}} \Big) \\ Q\_{\text{GW}\to\text{SW}} = 0 \\ Q\_{\text{SW}\to\text{GW}} = K\_{\text{bed}}(\text{WL}) \Big( \frac{h\_{\text{star}\text{max}} - z\_{\text{bed}}}{d\_{\text{bed}}} \Big) \end{cases} \tag{17}$$

Pumping

Pumping rates *Qpump* (m<sup>3</sup> /day) are specified by the model user. These can be specified for any grid cell, for any day of the simulation in the *gwflow.aqu* input file.

#### Saturation Excess Flow

Saturation excess flow occurs when groundwater head *h* rises above the ground surface, typically during rainfall events that rapidly increase the water table. This condition is tested for each cell during each daily time step. If *h* > ground surface elevation, the volumetric flux (m<sup>3</sup> /day) of groundwater excess flow is given by:

$$Q\_{\text{state}} = \left(h\_{\text{i,j}} - z\_{\text{surf}\_{\text{i,j}}}\right) (\Delta x \Delta y) S\_{y\_{\text{i,j}}} \tag{18}$$

The water is removed from the grid cell and added to the closest stream channel on that same day.

#### 2.2.4. SWAT+ Code Structure with the *gwflow* Module

The structure of the SWAT+ modeling code with the embedded *gwflow* module is presented in Figure 3. The *gwflow* module adds only two subroutines to the SWAT+ code: *gwflow\_read*, which reads data for all grid cells from the input file *gwflow.aqu*, and *gwflow\_simulate*, which performs the computations of groundwater fluxes and changes in groundwater head and groundwater storage. The subroutine *hyd\_read\_connect*, which reads connections between all spatial objects in the SWAT+ model, also reads the input file *gwflow.con*, which contains a list of connections between grid cells and SWAT+ channels to enable groundwater–surface water interactions (see Section 2.2.4). The inputs required in the *gwflow.aqu* and *gwflow.con* files will be discussed in Section 2.2.5.

**Figure 3.** Data flow of the Soil and Water Assessment Tool+ (SWAT+) modeling code. Subroutines for the new *gwflow* module are shaded in light red, input files for the *gwflow* module are highlighted in red text, and the calculations performed in the *gwflow* module are listed in bullet points for the **Figure 3.** Data flow of the Soil and Water Assessment Tool+ (SWAT+) modeling code. Subroutines for the new *gwflow* module are shaded in light red, input files for the *gwflow* module are highlighted in red text, and the calculations performed in the *gwflow* module are listed in bullet points for the *gwflow\_simulate* subroutine.

2.2.5. Required Inputs for the *gwflow* Module All data for the *gwflow* module are contained in three input files: *gwflow.con*, *gwflow.aqu*, and gwflow.*hrucell* (see Figure 3). *gwflow.con* contains a list of all grid cells connected to stream channels, with the ID number listed for each, so that *QSWGW* and *QGWSW* can be calculated for the correct cells. These cells are called "River Cells". The file *gwflow.aqu* contains general information for the module and a list of data values for each grid cell. A complete list of information in this file is shown in Table 2. Basic information includes the number of rows and columns in the grid, the cell size (m), flags for water table initiation (at the Following the reading of all input data, the *time\_control* subroutine is called, which controls the yearly and daily computation loops. For each daily time step of the simulation, the *command* subroutine is called, which loops through all objects (HRUs, Routing Units, channels) in the model, including the *gwflow* grid cells. When the *gwflow\_simulate* subroutine is called, all groundwater fluxes are calculated, whereupon the change in groundwater storage and groundwater head is calculated for each grid cell. Equations (see Sections 2.2.2 and 2.2.3) for each calculation are indicated in Figure 3. Head values, groundwater flux values for each grid cell, and a groundwater balance for the entire watershed model domain are then written to files. The groundwater balance is calculated and output for each day, each year, and then for the entire simulation (average annual).

#### beginning of the simulation) and the inclusion of saturation excess flow and groundwater ET, a recharge delay term (to transfer recharge from the bottom of the soil profile to the water table), and 2.2.5. Required Inputs for the *gwflow* Module

*gwflow\_simulate* subroutine.

the time step. The water table can be initiated by (1) ground surface minus a specified depth for each grid cell, (2) a fraction of the distance between the ground surface and the bedrock, or (3) a value specified for each cell. For *K* and *Sy*, the grid cells are divided into zones, with *K* and *Sy* specified for each zone. This facilitates calibration by changing only values for each zone rather than values for All data for the *gwflow* module are contained in three input files: *gwflow.con*, *gwflow.aqu*, and gwflow.*hrucell* (see Figure 3). *gwflow.con* contains a list of all grid cells connected to stream channels, with the ID number listed for each, so that *QSW*→*GW* and *QGW*→*SW* can be calculated for the correct cells. These cells are called "River Cells".

each grid cell. Groundwater head *h* for each cell can be output for any day of the simulation. Certain cells can also be designated as "observation cells", with *h* values for these cells output for each day of the simulation, resulting in time series that can be compared to data from groundwater monitoring wells. Many of these data can be obtained from national or global datasets. Sources for these data will be discussed in Section 2.3. The file *gwflow.aqu* contains general information for the module and a list of data values for each grid cell. A complete list of information in this file is shown in Table 2. Basic information includes the number of rows and columns in the grid, the cell size (m), flags for water table initiation (at the beginning of the simulation) and the inclusion of saturation excess flow and groundwater ET, a recharge delay term (to transfer recharge from the bottom of the soil profile to the water table), and the time

step. The water table can be initiated by (1) ground surface minus a specified depth for each grid cell, (2) a fraction of the distance between the ground surface and the bedrock, or (3) a value specified for each cell. For *K* and *Sy*, the grid cells are divided into zones, with *K* and *S<sup>y</sup>* specified for each zone. This facilitates calibration by changing only values for each zone rather than values for each grid cell. Groundwater head *h* for each cell can be output for any day of the simulation. Certain cells can also be designated as "observation cells", with *h* values for these cells output for each day of the simulation, resulting in time series that can be compared to data from groundwater monitoring wells. Many of these data can be obtained from national or global datasets. Sources for these data will be discussed in Section 2.3.


**Table 2.** Information required in the *gwflow.aqu* file for the *gwflow* module. Many of the parameters can be found from national or global datasets, as described in the model application in Section 2.3.

#### *2.3. Application to the Little River Experimental Watershed, Georgia*

#### 2.3.1. Overview of Study Region

The SWAT+ model using the new *gwflow* module is applied in this study to the 327 km<sup>2</sup> Little River Experimental Watershed (LREW), located within the Upper Suwannee River Basin in south-central Georgia (Figure 4). This watershed is ideal for testing the *gwflow* module due to the strong connection between the shallow unconfined aquifer and the many streams that make up the river system [31]. In addition, there are a number of streamflow gages and groundwater monitoring wells located throughout the watershed for model testing.

**Figure 4.** Geographic location and SWAT+ datasets and computational units, showing (**a**) digital elevation model (DEM) (m); (**b**) Hydrologic Response Units (HRUs), showing HRUs located in the upland areas and the floodplain areas; and (**c**) land use. **Figure 4.** Geographic location and SWAT+ datasets and computational units, showing (**a**) digital elevation model (DEM) (m); (**b**) Hydrologic Response Units (HRUs), showing HRUs located in the upland areas and the floodplain areas; and (**c**) land use.

The climate of the region is humid subtropical, with hot and humid summers and mild winters. Summers are characterized by high-intensity thunderstorms, leading to sharp increases in streamflow. The average annual rainfall is approximately 1200 mm, and the annual mean temperature is 19 °C. The average annual ET has been estimated to be approximately 70% of annual precipitation [37]. The watershed elevation ranges between 82 m and 148 m, with broad floodplains (Figure 4a). Land use (Figure 4c) consists of forest (50%), primarily pine trees in upland areas and hardwoods and evergreens in riparian areas; agriculture, primary cotton and peanut row-crops (41%); urban areas (7%), and open water (2%). Soils are loamy sands and sandy loams [39], which are underlain by the relatively impermeable Hawthorne formation which limits groundwater flow to deeper geologic layers [40]. The presence of the Hawthorne formation, referred to as the bedrock through the remainder of this paper, leads to significant lateral groundwater flow and associated groundwater discharge to streams [31,37,41]. Saturation excess flow conditions occur within the riparian areas [42]. The thickness of the alluvial aquifer is presented in Section 2.3.3. Many variations and applications of the SWAT and SWAT+ models have been applied to the LREW, including for model evaluation, calibration, and parameter sensitivity [43–46], analyzing the The climate of the region is humid subtropical, with hot and humid summers and mild winters. Summers are characterized by high-intensity thunderstorms, leading to sharp increases in streamflow. The average annual rainfall is approximately 1200 mm, and the annual mean temperature is 19 ◦C. The average annual ET has been estimated to be approximately 70% of annual precipitation [37]. The watershed elevation ranges between 82 m and 148 m, with broad floodplains (Figure 4a). Land use (Figure 4c) consists of forest (50%), primarily pine trees in upland areas and hardwoods and evergreens in riparian areas; agriculture, primary cotton and peanut row-crops (41%); urban areas (7%), and open water (2%). Soils are loamy sands and sandy loams [39], which are underlain by the relatively impermeable Hawthorne formation which limits groundwater flow to deeper geologic layers [40]. The presence of the Hawthorne formation, referred to as the bedrock through the remainder of this paper, leads to significant lateral groundwater flow and associated groundwater discharge to streams [31,37,41]. Saturation excess flow conditions occur within the riparian areas [42]. The thickness of the alluvial aquifer is presented in Section 2.3.3.

effect of conservation practices on water quality [47], landscape routing [48] and connectivity between upland areas and floodplains [28]. Indeed, the SWAT+ model was first published and introduced using an application to this watershed [2]. These applications, however, treated groundwater flow and groundwater–surface water interactions simplistically, based on steady flow assumptions and unconnected aquifer regions. The inclusion of the *gwflow* module is expected to simulate more realistic groundwater hydrologic fluxes, which can lead to improved assessments of riparian-stream hydrology and conservation practices within the study region. Many variations and applications of the SWAT and SWAT+ models have been applied to the LREW, including for model evaluation, calibration, and parameter sensitivity [43–46], analyzing the effect of conservation practices on water quality [47], landscape routing [48] and connectivity between upland areas and floodplains [28]. Indeed, the SWAT+ model was first published and introduced using an application to this watershed [2]. These applications, however, treated groundwater flow and groundwater–surface water interactions simplistically, based on steady flow assumptions and unconnected aquifer regions. The inclusion of the *gwflow* module is expected to simulate more realistic groundwater hydrologic fluxes, which can lead to improved assessments of riparian-stream hydrology and conservation practices within the study region.

### 2.3.2. SWAT+ Model Construction

The base SWAT+ model was constructed using QSWAT+, a QGIS interface for SWAT+ (https: //swat.tamu.edu/software/plus/). The interface guides the user through loading datasets (DEM, land use and crop data, soil type, and climate station data) and creating hydrologic response units (HRUs), routing units, and channels. The resolution and source for each data set are shown in Table 3. The interface then uses the SWAT+ Editor to write all necessary input files to run the SWAT+ modeling code. In this study, SWAT+ version 59.3 is used. The process resulted in 10796 HRUs and 343 channels. Figure 4b shows the HRUs divided into upland areas (yellow) and floodplains areas (blue). Channel delineation is shown in Section 2.3.3.

Before including the *gwflow* module, an initial calibration was performed for SWAT+ to provide hydrologic fluxes (surface runoff, lateral flow, groundwater discharge, ET) that are similar to those observed from field data assessment. This calibration was performed manually as well as using IPEAT+ (Integrated Parameter Estimation and Uncertainty Analysis Tool Plus) [49], an automated calibration tool developed for SWAT+, based on the Dynamically Dimensioned Search (DDS) algorithm [50], using the streamflow at the watershed outlet as testing data.



#### 2.3.3. Preparing the gwflow Module

Preparing the *gwflow* module for SWAT+ consists of populating the *gwflow.aqu*, *gwflow.hrucell*, and *gwflow.con* files. The following list provides details regarding the preparation of all required data for the *gwflow.aqu* file (see Table 2). All data can be prepared in a GIS (e.g., ArcMap, QGIS).


within the watershed boundary are "active" and assigned status = 1, except for cells that lie along the boundary and are designated as "boundary" cells and assigned status = 2. Boundary cells are simulated as constant-head cells, i.e., the groundwater head assigned to these cells at the beginning of the simulation remains constant throughout the simulation period.


This study uses 200 m grid cells for the *gwflow* module, resulting in 8647 active cells. The grid and associated *K* values are shown in Figure 5a. These *K* values (3 m/day, 5 m/day, 20 m/day) are the final values based on manual calibration (see Section 2.3.4), and were increased from the values (0.001 m/dy, 0.086 m/day, 1.369 m/day) provided by the global data set. Corresponding *S<sup>y</sup>* values are 0.10, 0.15, and 0.11. As shown in Figure 5a, there are only 3 zones of aquifer materials. *Kbed* and *dbed* were set to 0.001 m/day and 2 m, respectively, for all River Cells. The aquifer thickness (cm) (Figure 5b) ranged from 0.13 m to 88.2 m. The identified River Cells are shown in Figure 5c. *Hydrology* **2020**, *6*, x FOR PEER REVIEW 16 of 26

**Figure 5.** (**a**) gwflow grid cells showing SWAT+ channels and final hydraulic conductivity (m/day) values; (**b**) aquifer thickness (cm); and (**c**) grid cells designated as River Cells (i.e., cells that can **Figure 5.** (**a**) gwflow grid cells showing SWAT+ channels and final hydraulic conductivity (m/day) values; (**b**) aquifer thickness (cm); and (**c**) grid cells designated as River Cells (i.e., cells that can exchange water with SWAT+ channels).

#### 2.3.4. Overall Simulation

**3. Results and Discussion** 

*3.1. Water Balance* 

exchange water with SWAT+ channels).

The simulation is run for the 1980–2013 period, with the first two years used as a warm-up (i.e., the results from 1980 and 1981 are not compared to observed watershed data). Based on initial simulations, a time step ∆*t* of 0.25 days was required for the stability of the groundwater balance equation (Equation (10)). The run time of the model is only 20% higher than the original SWAT+ simulation due to the additional equations added to SWAT+ by the *gwflow* module. Simulated daily streamflow was compared to average daily streamflow at 8 gage sites, provided by the USDA-ARS

**Figure 6.** Map of the Little River Experimental Watershed (LREW) showing the location of the streamflow gages (B, O, N, F, I, J, K, M) and the U.S. Geological Survey (USGS) observation wells. The SWAT+ delineated channels also are shown, using color to indicate approximate channel width.

and its basic workings rather than prepare a model for scenario analysis.

Only manual calibration was performed in this study, to provide reasonable matches between simulated and observed streamflow and groundwater head. The Nash–Sutcliffe Efficiency (NSE) coefficient is used to quantify the performance of the model regarding streamflow simulation. For aquifer properties *K* and *Sy*, values were changed only according to the three identified zones (see Figure 5a) to remain consistent with the global data set, and therefore point-by-point calibration was not pursued to yield optimal matches between simulated and observed groundwater head at the 8 monitoring wells. Automated calibration could be performed to yield optimal matches for both groundwater head and streamflow. However, the intent of this paper is to present the *gwflow* module (Agricultural Research Service) Southeast Watershed Research Laboratory, and simulated daily groundwater head was compared to periodic groundwater levels measured at 8 monitoring wells, provided by the U.S. Geological Survey (USGS) (https://maps.waterdata.usgs.gov/mapper/index.html, accessed 15 February 2020). Figure 6 shows the location of the streamflow gages (green dots) and monitoring wells (red dots). **Figure 5.** (**a**) gwflow grid cells showing SWAT+ channels and final hydraulic conductivity (m/day) values; (**b**) aquifer thickness (cm); and (**c**) grid cells designated as River Cells (i.e., cells that can exchange water with SWAT+ channels).

*Hydrology* **2020**, *6*, x FOR PEER REVIEW 16 of 26

**Figure 6.** Map of the Little River Experimental Watershed (LREW) showing the location of the streamflow gages (B, O, N, F, I, J, K, M) and the U.S. Geological Survey (USGS) observation wells. The **Figure 6.** Map of the Little River Experimental Watershed (LREW) showing the location of the streamflow gages (B, O, N, F, I, J, K, M) and the U.S. Geological Survey (USGS) observation wells. The SWAT+ delineated channels also are shown, using color to indicate approximate channel width.

SWAT+ delineated channels also are shown, using color to indicate approximate channel width. Only manual calibration was performed in this study, to provide reasonable matches between simulated and observed streamflow and groundwater head. The Nash–Sutcliffe Efficiency (NSE) coefficient is used to quantify the performance of the model regarding streamflow simulation. For aquifer properties *K* and *Sy*, values were changed only according to the three identified zones (see Figure 5a) to remain consistent with the global data set, and therefore point-by-point calibration was not pursued to yield optimal matches between simulated and observed groundwater head at the 8 monitoring wells. Automated calibration could be performed to yield optimal matches for both groundwater head and streamflow. However, the intent of this paper is to present the *gwflow* module Only manual calibration was performed in this study, to provide reasonable matches between simulated and observed streamflow and groundwater head. The Nash–Sutcliffe Efficiency (NSE) coefficient is used to quantify the performance of the model regarding streamflow simulation. For aquifer properties *K* and *Sy*, values were changed only according to the three identified zones (see Figure 5a) to remain consistent with the global data set, and therefore point-by-point calibration was not pursued to yield optimal matches between simulated and observed groundwater head at the 8 monitoring wells. Automated calibration could be performed to yield optimal matches for both groundwater head and streamflow. However, the intent of this paper is to present the *gwflow* module and its basic workings rather than prepare a model for scenario analysis.

#### and its basic workings rather than prepare a model for scenario analysis. **3. Results and Discussion**

#### **3. Results and Discussion**  *3.1. Water Balance*

*3.1. Water Balance*  The overall water balance of the SWAT+ simulation is shown in Figure 7 using average annual flux values in mm (depth over the entire watershed). For simulations using the original SWAT+ code, fluxes would be limited to surface runoff (78 mm), lateral flow (57 mm), ET (843 mm), percolation and recharge (184 mm), and groundwater discharge. With the inclusion of the *gwflow* module, additional fluxes include boundary flow (0 mm), groundwater discharge (6 mm), stream seepage (1 mm), groundwater ET (11 mm), and saturation excess flow (161 mm), with state variables groundwater head (average of 107.2 m for the watershed) and saturated thickness (average of 19.2 m). Surface runoff (78 mm) and lateral flow (57 mm) account for 12% of annual precipitation, whereas groundwater accounts for 14% of annual precipitation. Water yield is 26% of annual precipitation, close to the observed value of

27% [31]. Groundwater entering streams occurs primarily via saturation excess flow (161 mm) rather than groundwater discharge (6 mm). water exchange, and groundwater ET. Saturation excess flow occurs primarily in the riparian areas, as does groundwater ET, with the latter due to riparian trees extracting water from the water table [54].

over the 1980–2013 simulation period for recharge, saturation excess flow, groundwater–surface

*Hydrology* **2020**, *6*, x FOR PEER REVIEW 17 of 26

saturation excess flow (161 mm) rather than groundwater discharge (6 mm).

shown. The percent error in the groundwater balance is < 0.00001%.

The overall water balance of the SWAT+ simulation is shown in Figure 7 using average annual flux values in mm (depth over the entire watershed). For simulations using the original SWAT+ code, fluxes would be limited to surface runoff (78 mm), lateral flow (57 mm), ET (843 mm), percolation and recharge (184 mm), and groundwater discharge. With the inclusion of the *gwflow* module, additional fluxes include boundary flow (0 mm), groundwater discharge (6 mm), stream seepage (1 mm), groundwater ET (11 mm), and saturation excess flow (161 mm), with state variables groundwater head (average of 107.2 m for the watershed) and saturated thickness (average of 19.2 m). Surface runoff (78 mm) and lateral flow (57 mm) account for 12% of annual precipitation, whereas groundwater accounts for 14% of annual precipitation. Water yield is 26% of annual precipitation, close to the observed value of 27% [31]. Groundwater entering streams occurs primarily via

Total water inputs include rainfall (1167 mm), total outputs equal 1155 mm, and the total change in groundwater storage and soil water storage is 15 mm and 1 mm, respectively, resulting in a water balance error of 0.4%. A time series of volumetric amounts (m3 × 106) for all hydrologic fluxes is shown in Figure 8a, with positive values indicating water entering the watershed, and negative values indicating water leaving the watershed. The year-by-year subsurface storage volume (groundwater + soil water) also is shown. A similar time series is shown in Figure 8b for the aquifer system, with positive values indicating water entering the aquifer (sources), and negative values indicating groundwater leaving the aquifer (sinks). The year-by-year groundwater storage volume also is

The aquifer system time series (Figure 8b) shows that recharge is largely balanced by saturation excess flow, i.e., recharge events during large storm events produce saturation excess flow near the

**Figure 7.** Average annual water balance for 1980-2013 for the LRW SWAT+ simulation, with fluxes shown in mm/yr. Fluxes groundwater ET, saturation excess flow, stream seepage, groundwater discharge, boundary flow, and recharge and state variables groundwater head (m) and saturated thickness (m) are simulated using the new *gwflow* module of SWAT+. **Figure 7.** Average annual water balance for 1980–2013 for the LRW SWAT+ simulation, with fluxes shown in mm/yr. Fluxes groundwater ET, saturation excess flow, stream seepage, groundwater discharge, boundary flow, and recharge and state variables groundwater head (m) and saturated thickness (m) are simulated using the new *gwflow* module of SWAT+.

Total water inputs include rainfall (1167 mm), total outputs equal 1155 mm, and the total change in groundwater storage and soil water storage is 15 mm and 1 mm, respectively, resulting in a water balance error of 0.4%. A time series of volumetric amounts (m<sup>3</sup> <sup>×</sup> <sup>10</sup><sup>6</sup> ) for all hydrologic fluxes is shown in Figure 8a, with positive values indicating water entering the watershed, and negative values indicating water leaving the watershed. The year-by-year subsurface storage volume (groundwater + soil water) also is shown. A similar time series is shown in Figure 8b for the aquifer system, with positive values indicating water entering the aquifer (sources), and negative values indicating groundwater leaving the aquifer (sinks). The year-by-year groundwater storage volume also is shown. The percent error in the groundwater balance is < 0.00001%. *Hydrology* **2020**, *6*, x FOR PEER REVIEW 18 of 26

**Figure 8.** Time series of hydrologic flux volumes for (**a**) the entire watershed and (**b**) the aquifer system, for 1980-2013. For the legends, "gw" refers to groundwater, "ET" refers to evapotranspiration, "sw" refers to surface water, "gwsw" refers to groundwater discharge to streams, "swgw" refers to stream seepage to the aquifer, "sat excess" refers to saturation excess flow, and "boundary" refers to groundwater flow across the watershed boundary. **Figure 8.** Time series of hydrologic flux volumes for (**a**) the entire watershed and (**b**) the aquifer system, for 1980-2013. For the legends, "gw" refers to groundwater, "ET" refers to evapotranspiration, "sw" refers to surface water, "gw→sw" refers to groundwater discharge to streams, "sw→gw" refers to stream seepage to the aquifer, "sat excess" refers to saturation excess flow, and "boundary" refers to groundwater flow across the watershed boundary.

#### *Hydrology* **2020**, *7*, 75

The aquifer system time series (Figure 8b) shows that recharge is largely balanced by saturation excess flow, i.e., recharge events during large storm events produce saturation excess flow near the streams. This is demonstrated further in the maps (Figure 9) of spatial varying volumetric fluxes (m<sup>3</sup> ) over the 1980–2013 simulation period for recharge, saturation excess flow, groundwater–surface water exchange, and groundwater ET. Saturation excess flow occurs primarily in the riparian areas, as does groundwater ET, with the latter due to riparian trees extracting water from the water table [54]. *Hydrology* **2020**, *6*, x FOR PEER REVIEW 19 of 26

**Figure 9.** Maps of total volume (m3) of groundwater stress flux for the 1980-2013 simulation period: (**a**) recharge, (**b**) lateral flow, (**c**) groundwater–surface water exchange, (**d**) groundwater ET, and (e) **Figure 9.** Maps of total volume (m<sup>3</sup> ) of groundwater stress flux for the 1980-2013 simulation period: (**a**) recharge, (**b**) lateral flow, (**c**) groundwater–surface water exchange, (**d**) groundwater ET.

#### saturation excess flow. *3.2. Streamflow Results*

*3.2. Streamflow Results*  Time series of daily simulated and observed streamflow (m3/sec) are shown in Figure 10 for four of the stream gage sites from the upstream to downstream end of the watershed. The NSE value for each of the eight sites is shown in the inset table. Although several sites show adequate values (J, I, F, B: ≥ 0.30), others show poor matches (K, O: 0.22 and 0.27) and two are below 0 (M, N: −0.24, −0.37). Note that sites K, O, M and N are along small tributaries of the river system with moderate (< 15 m3/s) or low flows (< 5 m3/s). However, visual comparisons (see time series charts for M and O in Figure 10) indicate that the model does track the timing of measured flow. This is encouraging, as these flows are controlled to a large degree by subsurface flows calculated by the *gwflow* module. Detailed calibration of subbasin parameters for these tributary regions would be necessary if the model is to be used for scenario analysis. As this model application is intended solely for *gwflow* module demonstration, the results are adequate. Time series of daily simulated and observed streamflow (m<sup>3</sup> /sec) are shown in Figure 10 for four of the stream gage sites from the upstream to downstream end of the watershed. The NSE value for each of the eight sites is shown in the inset table. Although several sites show adequate values (J, I, F, B: ≥ 0.30), others show poor matches (K, O: 0.22 and 0.27) and two are below 0 (M, N: −0.24, −0.37). Note that sites K, O, M and N are along small tributaries of the river system with moderate (< 15 m<sup>3</sup> /s) or low flows (< 5 m<sup>3</sup> /s). However, visual comparisons (see time series charts for M and O in Figure 10) indicate that the model does track the timing of measured flow. This is encouraging, as these flows are controlled to a large degree by subsurface flows calculated by the *gwflow* module. Detailed calibration of subbasin parameters for these tributary regions would be necessary if the model is to be used for scenario analysis. As this model application is intended solely for *gwflow* module demonstration, the results are adequate.

*Hydrology* **2020**, *6*, x FOR PEER REVIEW 20 of 26

**Figure 10.** Time series of daily simulated and observed streamflow (m3/s) for the eight stream gage locations in the LREW study region (see Figure 6). The circled gage *B* is the outlet of the watershed. **Figure 10.** Time series of daily simulated and observed streamflow (m<sup>3</sup> /s) for the eight stream gage locations in the LREW study region (see Figure 6). The circled gage *B* is the outlet of the watershed. The Nash–Sutcliffe Efficiency (NSE) coefficient is shown on the chart for each stream gage site.

The Nash–Sutcliffe Efficiency (NSE) coefficient is shown on the chart for each stream gage site.

The differences in simulated streamflow between the original SWAT+ model and the SWAT+ model with the *gwflow* module are shown in Figure 11 for the outlet (site B) for several years of the simulation. As seen by the time series, SWAT+ simulates near-zero flow between rainfall events, consistent with the measured streamflow, whereas SWAT+ with *gwflow* simulates low flows (1–3 m3/s) during these times. Groundwater discharge in SWAT+ is simulated in the "pulse" form, based on the magnitude of rainfall events and the resulting recharge. Using the *gwflow* module, however, groundwater can enter the streams at any time based on groundwater gradients, leading to small flows between rainfall events. The results of SWAT+ are more accurate in this instance in this regard, as observed flows also approach zero between rainfall events during certain times of the year. As such, further calibration is required for the *gwflow* module if these low-flow periods are of importance for the intended use of the model. However, the inclusion of the *gwflow* module allows groundwater discharge to streams, either via the streambed or via saturation excess flow, to be simulated in a more physically realistic manner for scenario analysis (e.g., climate change, water conservation, nutrient The differences in simulated streamflow between the original SWAT+ model and the SWAT+ model with the *gwflow* module are shown in Figure 11 for the outlet (site B) for several years of the simulation. As seen by the time series, SWAT+ simulates near-zero flow between rainfall events, consistent with the measured streamflow, whereas SWAT+ with *gwflow* simulates low flows (1–3 m<sup>3</sup> /s) during these times. Groundwater discharge in SWAT+ is simulated in the "pulse" form, based on the magnitude of rainfall events and the resulting recharge. Using the *gwflow* module, however, groundwater can enter the streams at any time based on groundwater gradients, leading to small flows between rainfall events. The results of SWAT+ are more accurate in this instance in this regard, as observed flows also approach zero between rainfall events during certain times of the year. As such, further calibration is required for the *gwflow* module if these low-flow periods are of importance for the intended use of the model. However, the inclusion of the *gwflow* module allows groundwater discharge to streams, either via the streambed or via saturation excess flow, to be simulated in a more physically realistic manner for scenario analysis (e.g., climate change, water conservation, nutrient management).

#### management). *3.3. Groundwater Head Results*

Groundwater head (m) at the end of the simulation (end of 2013) for each grid cell is shown in Figure 12. The groundwater head (i.e., water table elevation) mimics the ground surface elevation (see Figure 4a). Groundwater head fluctuations are also shown at the eight monitoring well locations, in comparison with the observed groundwater head values from the USGS data set. The dotted brown line on each time series chart represents the ground surface at that site. Observed rapid groundwater head fluctuations generally are captured by the model. The rapid changes occur due to recharge events and resulting lateral flow and/or saturation excess flow.

*Hydrology* **2020**, *6*, x FOR PEER REVIEW 21 of 26

**Figure 11.** Comparison of daily streamflow at the outlet of the watershed (site B in Figure 6) simulated for the original SWAT+ model and the SWAT+ model with the new *gwflow* module. A semi-log plot is used to show the differences between the model results for low flows. The measured streamflow is also shown. **Figure 11.** Comparison of daily streamflow at the outlet of the watershed (site B in Figure 6) simulated for the original SWAT+ model and the SWAT+ model with the new *gwflow* module. A semi-log plot is used to show the differences between the model results for low flows. The measured streamflow is also shown. in comparison with the observed groundwater head values from the USGS data set. The dotted brown line on each time series chart represents the ground surface at that site. Observed rapid groundwater head fluctuations generally are captured by the model. The rapid changes occur due to recharge events and resulting lateral flow and/or saturation excess flow.

**Figure 12.** Map of groundwater head (m) for each grid cell at the end of the simulation period, and time series of daily simulated groundwater head (m) and periodically observed groundwater head (m) for the eight USGS groundwater monitoring well locations (see Figure 6). The dotted brown line on each time series chart represents the ground surface at that site. MAE = Mean Absolute Error (m), i.e., the average residual between the observed and simulated groundwater head values. **Figure 12.** Map of groundwater head (m) for each grid cell at the end of the simulation period, and time series of daily simulated groundwater head (m) and periodically observed groundwater head (m) for the eight USGS groundwater monitoring well locations (see Figure 6). The dotted brown line on each time series chart represents the ground surface at that site. MAE = Mean Absolute Error (m), i.e., the average residual between the observed and simulated groundwater head values.

**Figure 12.** Map of groundwater head (m) for each grid cell at the end of the simulation period, and time series of daily simulated groundwater head (m) and periodically observed groundwater head (m) for the eight USGS groundwater monitoring well locations (see Figure 6). The dotted brown line on each time series chart represents the ground surface at that site. MAE = Mean Absolute Error (m), i.e., the average residual between the observed and simulated groundwater head values. The mean absolute error (MAE) of simulated head values that correspond in location and timing to measured head values is 1.7 m. All sites have an MAE of 3.0 m or less. Two locations (313144083335501, 313238083331901) show excellent tracking of the measured head by the model, while others show moderate agreement. For several locations, the model underpredicts head levels, The mean absolute error (MAE) of simulated head values that correspond in location and timing to measured head values is 1.7 m. All sites have an MAE of 3.0 m or less. Two locations (313144083335501, 313238083331901) show excellent tracking of the measured head by the model, while others show moderate agreement. For several locations, the model underpredicts head levels, The mean absolute error (MAE) of simulated head values that correspond in location and timing to measured head values is 1.7 m. All sites have an MAE of 3.0 m or less. Two locations (313144083335501, 313238083331901) show excellent tracking of the measured head by the model, while others show moderate agreement. For several locations, the model underpredicts head levels, and at two locations (313435083390101, 313630083385001) the model overestimates head levels. The underestimation and overestimation could be due to the misrepresentation of aquifer properties (*K*, *Sy*) or underestimation of localized recharge by HRUs. Likely the aquifer properties are not represented correctly for these localized areas. However, as with the streamflow results shown in Section 3.2, in this study aquifer property values are constrained by the zones established by the global permeability map (see Figure 5a). Point-by-point calibration could be performed to provide better matches between the observed and simulated head values; however, the purpose of this paper is to present the *gwflow* module and to highlight its capabilities. If the model is to be used for scenario analysis, localized aquifer conditions can be better represented using borehole lithology data. Such an approach would also use automated parameter estimation methods and employ ensemble-based uncertainty analysis and sensitivity analysis methods (e.g., [21–23]) to explore the effect of parameter values (e.g., hydraulic

conductivity: [24,25]) on system-response variables such as water table elevation, groundwater–surface water exchange rates, and streamflow. uncertainty analysis and sensitivity analysis methods (e.g., [21–23]) to explore the effect of parameter values (e.g., hydraulic conductivity: [24,25]) on system-response variables such as water table

an approach would also use automated parameter estimation methods and employ ensemble-based

*Hydrology* **2020**, *6*, x FOR PEER REVIEW 22 of 26

and at two locations (313435083390101, 313630083385001) the model overestimates head levels. The underestimation and overestimation could be due to the misrepresentation of aquifer properties (*K*, *Sy*) or underestimation of localized recharge by HRUs. Likely the aquifer properties are not represented correctly for these localized areas. However, as with the streamflow results shown in Section 3.2, in this study aquifer property values are constrained by the zones established by the global permeability map (see Figure 5a). Point-by-point calibration could be performed to provide better matches between the observed and simulated head values; however, the purpose of this paper

Finally, depth to water table (m) and saturated thickness (m) for each grid cell at the end of the simulation are shown in Figure 13. The map of depth to water table shows depths near 0 or at 0 for much of the riparian and floodplain areas, indicating near-saturated conditions, leading to groundwater ET (see Figure 9d) and saturation excess flow (see Figure 9b). The map of saturated thickness shows spatial patterns similar to the map of aquifer thickness (see Figure 5b). Such maps can assist in determining groundwater storage throughout the aquifer system and the strong interplay between the land surface and subsurface hydrologic processes. elevation, groundwater–surface water exchange rates, and streamflow. Finally, depth to water table (m) and saturated thickness (m) for each grid cell at the end of the simulation are shown in Figure 13. The map of depth to water table shows depths near 0 or at 0 for much of the riparian and floodplain areas, indicating near-saturated conditions, leading to groundwater ET (see Figure 9d) and saturation excess flow (see Figure 9b). The map of saturated thickness shows spatial patterns similar to the map of aquifer thickness (see Figure 5b). Such maps can assist in determining groundwater storage throughout the aquifer system and the strong interplay between the land surface and subsurface hydrologic processes.

**Figure 13.** (**a**) Depth to Water Table (m) and (**b**) Saturated Thickness (m) for each grid cell at the end of the simulation period. Depth to Water Table was calculated by subtracting groundwater head (see Figure 12 for map) from ground surface elevation for each grid cell. Saturated Thickness was calculated by subtracting bedrock elevation from groundwater head. **Figure 13.** (**a**) Depth to Water Table (m) and (**b**) Saturated Thickness (m) for each grid cell at the end of the simulation period. Depth to Water Table was calculated by subtracting groundwater head (see Figure 12 for map) from ground surface elevation for each grid cell. Saturated Thickness was calculated by subtracting bedrock elevation from groundwater head.

#### **4. Summary and Conclusions**

**4. Summary and Conclusions**  This paper presents a new physically-based spatially-distributed groundwater flow module (*gwflow*) for the SWAT+ watershed model. The module is embedded in the SWAT+ modeling code as with other hydrologic modules and is intended to replace the current SWAT+ aquifer module [2] for watersheds with strong stream-aquifer connections. The model accounts for recharge from SWAT+ HRUs, lateral flow within the aquifer, ET from shallow groundwater, pumping, groundwater– surface water interactions through the streambed, and saturation excess flow. Daily groundwater head and groundwater storage are solved using an explicit numerical method approach for the groundwater balance equation, with head and flow values for the current day based on head and This paper presents a new physically-based spatially-distributed groundwater flow module (*gwflow*) for the SWAT+ watershed model. The module is embedded in the SWAT+ modeling code as with other hydrologic modules and is intended to replace the current SWAT+ aquifer module [2] for watersheds with strong stream-aquifer connections. The model accounts for recharge from SWAT+ HRUs, lateral flow within the aquifer, ET from shallow groundwater, pumping, groundwater–surface water interactions through the streambed, and saturation excess flow. Daily groundwater head and groundwater storage are solved using an explicit numerical method approach for the groundwater balance equation, with head and flow values for the current day based on head and flow values from the previous day. The module outputs groundwater head and flux rates for each groundwater stress for analysis and mapping.

The modified SWAT+ model is applied to the 327 km<sup>2</sup> Little River Experimental Watershed in southern Georgia, USA, demonstrating its capabilities of simulating both surface and subsurface hydrological processes. Results from an uncalibrated model are compared against measured streamflow and groundwater levels at eight stream gage sites and eight groundwater monitoring locations, respectively. The model performed adequately regarding visual comparisons of time series and quantitative statistics but can be improved with additional calibration if it is to be used for scenario analysis. For this watershed, groundwater additions to the stream system are governed by saturation excess flow rather than discharge via the streambed, based on continual near-saturation conditions in

the riparian areas and measured low baseflow between rainfall events. For this watershed, inclusion of the *gwflow* module increases SWAT+ model run time by approximately 20%.

The modified SWAT+ modeling code with the *gwflow* module can be an important tool for modeling hydrological fluxes in watersheds that have a strong connection between the shallow unconfined aquifer and the stream system. Using this code can provide physically realistic groundwater flow gradients, flux values, and interactions with streams, which are important for assessing the effects of conservation practices on hydrology and nutrient transport.

**Author Contributions:** Conceptualization, R.T.B. and J.G.A.; methodology, R.T.B.; software, R.T.B. and J.G.A.; validation, R.T.B., K.B., J.G.A., and D.D.B.; formal analysis, R.T.B.; investigation, R.T.B. and K.B.; resources, K.B. and D.D.B.; data curation, K.B. and D.D.B.; writing—original draft preparation, R.T.B.; writing—review and editing, K.B., J.G.A., and D.D.; visualization, R.T.B. and K.B.; supervision, R.T.B.; project administration, R.T.B. and J.G.A.; funding acquisition, R.T.B. and J.G.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the United State Department of Agriculture–Agricultural Research Service through Cooperative Agreement number 59-3098-8-002 and a grant from the Agriculture and Food Research Initiative of the USDA National Institute of Food and Agriculture, grant number #2017-67019-26332.

**Acknowledgments:** The authors wish to thank several anonymous reviewers for their suggestions which greatly improved the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
