Quoted from: Bailey, Ryan T., Tyler C. Wible, Mazdak Arabi, Rosemary M. Records, and Jeffrey Ditty. "Assessing regional‐scale spatio‐temporal patterns of groundwater–surface water interactions using a coupled SWAT‐MODFLOW model." Hydrological Processes 30, no. 23 (2016): 4420-4433. https://doi.org/10.1002/hyp.10933 

      The coupled SWAT‐MODFLOW framework developed in this study combines the updated SWAT 2012 model (Revision 591) with MODFLOW‐NWT. Within the framework, SWAT simulates land surface processes, crop growth, in‐stream processes and soil zone processes, whereas MODFLOW‐NWT simulates three‐dimensional groundwater flow and all associated sources and sinks (e.g. recharge, pumping, discharge to tile drains and interaction with stream network). Both modelling codes are combined into a single FORTRAN code that is compiled and run as a single executable.

      The basic process of linking the SWAT and MODFLOW models is to pass HRU‐calculated deep percolation (i.e. water that exits the bottom of the soil profile) as recharge to the grid cells of MODFLOW, and then pass MODFLOW‐calculated groundwater–surface water interaction fluxes to the stream channels of SWAT. With this approach, SWAT calculates the volume of overland flow and soil lateral flow to streams, MODFLOW calculates the volume of groundwater discharge to streams, and then SWAT routes the water through the stream network of the watershed. Groundwater–surface water interaction is simulated using the River package of MODFLOW, with Darcy's law used to calculate the volumetric flow of water Qleak [L3/T] through the cross‐sectional flow area between the aquifer and the stream channel

(1)

where Kbed is river bed hydraulic conductivity [L/T], Lstr is the length of the stream [L], Pstr is the wetted perimeter of the stream [L], hstr is river stage [L], hgw is the hydraulic head of groundwater [L], and zbed is the thickness of the river bed [L]. Qleak is negative if groundwater flows to the river (i.e. groundwater hydraulic head hgw is above the river stage hstr), and positive if river water seeps into the aquifer. These calculations are performed for any grid cell through which a stream passes. The SWAT‐MODFLOW model simulation and linking process is summarized in Figure 2. Upon reading input data for both the SWAT and MODFLOW models, the simulation runs through the repeated daily process of SWAT HRU calculations, passing data to MODFLOW, running MODFLOW, passing data to SWAT and routing water through the watershed's stream network.

      Data are passed between the models using “mapping” subroutines that relate HRUs to MODFLOW grid cells and MODFLOW River cells to SWAT stream channels (Figure 3). The main elements of this mapping scheme are: HRUs; Disaggregated HRUs (DHRUs), which divide each original HRU into individual, contiguous areas within a sub‐basin to allow HRU calculations to be geo‐located; MODFLOW grid cells; MODFLOW River cells; and SWAT stream channels. Calculated deep percolation (i.e. recharge) for HRUs are first mapped to each individual DHRU, and then mapped to each MODFLOW grid cell according to the percent area of the DHRU contained within the grid cell for use by the Recharge package. SWAT‐calculated channel depth from each sub‐basin is mapped to the group of River cells within the sub‐basin for use by the River package. MODFLOW then computes groundwater hydraulic head and groundwater–surface water interactions, which are passed to SWAT. Groundwater discharge volumes, computed on a cell‐by‐cell basis within MODFLOW, are summed and added to in‐stream flow for each of SWAT's sub‐basins. SWAT then completes the stream routing calculations for the day, with the daily loop continuing until the end of the simulation.

      For the possible scenario of a River cell intersecting more than one stream, the length of each stream within the cell is used to calculate a composite weighted value of channel depth for use by MODFLOW and to distribute the cell groundwater discharge volume to the associated sub‐basin main channels. Within this scheme, MODFLOW is called as a subroutine within the SWAT framework, providing a single compiled FORTRAN code.

      The spatial relationship between HRUs, DHRUs, MODFLOW grid cells and MODFLOW River cells is presented in Figure 3 and Table 1 for a hypothetical watershed with four original SWAT HRUs overlying a MODFLOW finite difference grid. SWAT HRUs #1–3 each are designated as a single DHRU each, because they are spatially contiguous, whereas SWAT HRU #4 is separated into three separate DHRUs (Table 1). For simplicity, DHRU‐to‐cell mapping is only shown for DHRU #6 Table 1, which is the northern‐most portion of HRU #4 (Figure 3). As seen in Table 1, 12 grid cells intersect DHRU #6 (the eastern‐most polygon of HRU #4 in Figure 3), with the portion of the DHRU overlapped by each cell shown in the “Portion of DHRU” column and the portion of the cell within the DHRU shown in the “Portion of Grid Cell” column. These relationships for each HRU, DHRU, and grid cells and also a list of the River cells contained within each sub‐basin, are read in at the beginning of a SWAT‐MODFLOW simulation and stored in memory for use at each time step (Figure 2).

      The linkage procedure allows for SWAT and MODFLOW models of different spatial extents, that is, one of the models can extend beyond the boundaries of the other. This was performed to facilitate linking of existing models, such as regional‐scale groundwater models (e.g. Rumman and Payne, 2003; Christenson et al., 2011; Paschke, 2011; Gannett et al., 2012; Mashburn et al., 2013). Beyond the overlap area, the original functionality of each model is retained. Stress period and time stepping information specified in an original MODFLOW model is used to determine calls to read stress data (e.g. pumping) from input files. If the MODFLOW time step is less than 1 day, the code calculates and executes the required number of time steps for the single day. If the time step is greater than 1 day, the code forces MODFLOW to run with a daily time step. The code allows for the MODFLOW subroutine to be called at any specified frequency, for example, every week, every month, etc…depending on the assumed hydraulic connection between the soil zone and the water table. However, this functionality was not used in this study. Also, the code accounts for and converts from time and length units specified in the original MODFLOW model.