Assessment of Surface Water Resources in the Big Sunflower River Watershed Using Coupled SWAT – MODFLOW Model

The groundwater level in the Big Sunflower River Watershed (BSRW) in the U.S. has declined significantly in the past 30 years. Therefore, it is imperative to assess surface water resources (SWR) availability in BSRW to mitigate groundwater use for irrigation. This research applied the coupled Soil and Water Assessment Tool–Modular Groundwater Flow model (SWAT–MODFLOW) to assess SWR in BSRW. This study aimed at: (1) Assessing the reliability of SWAT–MODFLOW in BSRW, (2) analyzing temporal and spatial variations of SWR, and (3) assessing the potential availability of SWR in BSRW. Calibration and validation results showed that SWAT–MODFLOW can well simulate streamflow and groundwater levels in BSRW. Our results showed that BSRW had lower average monthly total stream resources (MSR = 8.8 × 107 m3) in growing seasons than in non-growing seasons (MSR = 11.0 × 107 m3), and monthly pond resources (MPR from 30,418 to 30,494 m3) varied less than stream resources. The proportion of sub-basins in BSRW with stream water resources greater than 700 mm was 21% in dry years (229 to 994 mm), while this increased to 35% in normal years (296 to 1141 mm) and 57% in wet years (554 to 991 mm). The Water Stress Index (WSI) ranged from 0.4 to 2.1, revealing that most of the sub-basins in BSRW have net SWR available for irrigation. Our results suggested that surface water resources might be supplementary irrigation sources to mitigate the water resources scarcity in this region.


Introduction
Groundwater levels are decreasing, and this is a serious problem in many regions of the United States, including Mississippi, where the predominant groundwater consumption is for agriculture irrigation.The Mississippi River alluvial plain in northwestern Mississippi (the Mississippi Delta) is one of the most important agricultural regions in the United States, with nearly 81% of the area being crops, including soybean, corn, cotton, and rice [1].Over the past 30 years, this region has been impacts on agro-ecosystems in the Fort Cobb Reservoir experimental watershed of southwestern Oklahoma, USA.Ni and Parajuli [39] used the SWAT-MODFLOW to determine the effects of management practices on surface and groundwater resources in the BSRW.Chunn et al. [40] also used the SWAT-MODFLOW to study the influence of climate change on the interactions between surface and groundwater.SWAT-MODFLOW was adopted in our study for the following two reasons: (1) SWAT-MODFLOW integrates the surface water hydrology and groundwater hydrology with the consideration of pond resources, and (2) SWAT or MODFLOW have been individually applied in the Mississippi Delta, which can provide some available information for the model calibration and validation [2,3,41].
Limited information is available regarding the temporal and spatial variations of SWR in the Mississippi Delta, especially in BSRW (Figure 1).The objectives of this study were to: (1) Assess the reliability of the SWAT-MODFLOW model in BSRW based on calibration and validation; (2) analyze the temporal and spatial variations of stream and pond water resources in BSRW using the coupled SWAT-MODFLOW model, and (3) estimate the potential availability of SWR for agriculture irrigation in BSRW.

Study Area
The BSRW is located in north-western Mississippi, USA, and lies between latitude 32° to 35° N and longitude 91° to 90° E, with an area of 10,488 km 2 (Figure 1).It is relatively flat, and slope ranges from 5.7 to 7 cm km −1 in BSRW [42,43].The BSRW is a highly productive agricultural region in the Mississippi Delta, owing to fertile soils and long growing seasons.Nearly 80% of the total land use in BSRW is used for crops, including soybean, corn, rice, and cotton [44].Over 90% of the soil types in this region are Alligator, Sharkey, Dowling, Forestdale, and Dundee, which have high clay and silt content [41].The climate of BSRW is subtropical, with a mean annual precipitation of 1371 mm and

Study Area
The BSRW is located in north-western Mississippi, USA, and lies between latitude 32 • to 35 • N and longitude 91 • to 90 • E, with an area of 10,488 km 2 (Figure 1).It is relatively flat, and slope ranges from 5.7 to 7 cm km −1 in BSRW [42,43].The BSRW is a highly productive agricultural region in the Mississippi Delta, owing to fertile soils and long growing seasons.Nearly 80% of the total land use in BSRW is used for crops, including soybean, corn, rice, and cotton [44].Over 90% of the soil types in this region are Alligator, Sharkey, Dowling, Forestdale, and Dundee, which have high clay and silt content [41].The climate of BSRW is subtropical, with a mean annual precipitation of 1371 mm and annual average temperature of 18 • C between 2000 and 2016 [32].Only 30% of annual rainfall occurs during the growing seasons.As a result, irrigation is needed for soybean, corn, cotton, and rice in this region.Carvajal et al. [45] reported that impounded water in on-farm ponds is one of the available resources for supplying water for the irrigation required by those crops.Ouyang et al. [23] also demonstrated that farm ponds can provide multiple benefits, such as recharging groundwater, storing rain water, and collecting runoff.Due to the increasing demand of surface and groundwater resources in the Mississippi Delta, many farm ponds have been recently built for agriculture irrigation in this region [23].

Definition of Stream Water Resource, Pond Water Resource, and SWR
SWR in this study is defined as the volume or depth of stream and pond water in each sub-basin [46], and the unit of SWR is given as m 3 or mm [8].The volume or depth of stream and pond water was calculated by the following equations [46]: Depth pond = V pond SA subba sin (4) Depth SWR = Depth stream + Depth pond (6) where V stream is the volume of stream water resources in each sub-basin (m 3 ), A stream is the cross-sectional area of stream in each sub-basin (m 2 ), L stream is the stream length (km) in each sub-basin, Depth stream is the depth of stream water resources in each sub-basin for a given time step, Z stream is the inverse of the stream side slope, W btm is the bottom width of the stream in each sub-basin (m), V pond is the volume of pond water resources at the end of the day (m 3 ), V stored is the volume of water stored in the pond water body at the beginning of the day (m 3 ), V flowin is the volume of pond water entering the water body during the day (m 3 ), V flowout is the volume of pond water flowing out of the water body during the day (m 3 ), V pcp is the volume of precipitation falling on the pond water body during the day, V evap is the volume of water removed from the pond water body by evaporation during the day (m 3 ), V seep is the volume of water lost from the pond water body by seepage (m 3 ), Depth pond is the depth of pond resources during the day in each sub-basin, and SA sub-basin is the surface area of each sub-basin.
V SWR and Depth SWR are the volume and depth of SWR in each sub-basin, respectively.There are only three gauge stations in BSRW.Therefore, we applied SWAT-MODFLOW to estimate availability of stream and pond water resources on a sub-basin scale.system (GIS) interface for preparing the input data and running the model [47].The input data for the SWAT model consist of a digital elevation model (DEM), land use, soil chemicals, physical and hydrological properties, weather data (daily precipitation, maximum and minimum air temperature, solar radiation, average daily wind speed, and daily relative humidity), and land management and operation information.SWAT can be applied to simulate the volume and depth of stream and pond water [47] by dividing the watershed into sub-basins, which can be further categorized into a number of hydrological response units (HRU) with same land use, soil type, and slope.Groundwater recharge in SWAT can be calculated as: where w rchrg,i is the amount of recharge entering the aquifers on day i (mm), δ gw (days) is the delay time or drainage time of the overlying geologic formations, w seep is the total amount of water exiting the bottom of the soil profile on day i (mm), and w rchrg,i−1 is the amount of recharge entering the aquifers on day i−1 (mm).

MODFLOW
MODFLOW is a finite-difference groundwater model based on the physical and mathematical concepts of groundwater flow theory [48].It can be used to simulate groundwater recharge, vadose zone percolation, evapotranspiration, pumping, and interactions between surface and groundwater.The groundwater aquifer system in MODFLOW is discretized into finite-difference cells, with a linear equation written for each cell and solved for groundwater head at each time step [48].The simulation of MODFLOW is based on modular packages.For example, the river package in MODFLOW is used to simulate the exchange of water between the aquifer and river through the river bed.The well package is designed to simulate wells that withdraw water from or add water to the aquifer at a constant rate during a stress period.Groundwater level and discharge can be simulated on a cell scale as follows: where QRIV n is the flow between the river and the aquifer (mm), taken as positive if it is directed into the aquifer; HRIV n is the water level in the river (m); CRIV n is the hydraulic conductance of the river-aquifer interconnection (cm day −1 ); and h i,j,k is the groundwater level at the node in the cell underlying the river reach.MODFLOW can represent multiple aquifer conditions, including confined, unconfined, leaky, delayed yield, and variably confined/unconfined conditions.Both steady state and transient conditions can be simulated [49].

SWAT-MODFLOW
There are limitations if only SWAT or MODFLOW is applied.Therefore, we applied the coupled SWAT-MODFLOW, which was developed by [36].Compared to the previous SWAT-MODFLOW model, this model has a wide range of advantages.It has geographically located HRU and an efficient HRU grid mapping scheme.In addition, the unconfined groundwater-flow equation in the model was improved by using the Newton formulation of MODFLOW-2005 [36].The key procedure of the coupled SWAT-MODFLOW model simulation is: SWAT passes its deep percolation to MODFLOW as recharge in each grid cell, and MODFLOW provides groundwater discharge to SWAT in each sub-basin.There are three major procedures for linking SWAT and MODFLOW.Firstly, HRUs are spatially disaggregated in preprocessing GIS routines to provide a geographic location for each HRU.Secondly, these disaggregated HRUs (DHRUs) are then intersected with MODFLOW grid cells to specify HRU-cell relationships.Thirdly, the MODFLOW river cells are grouped by SWAT sub-basin so that MODFLOW-calculated groundwater/surface water exchange rates can be passed to the correct Water 2019, 11, 528 6 of 21 sub-basin channel in preparation for streamflow routing.Full details regarding SWAT-MODFLOW are provided in Bailey et al. [36].

SWAT Input
Topography data were obtained in the form of a digital elevation model (DEM) at 10 m resolution from the United States Geological Survey (USGS) (http://viewer.nationalmap.gov/basic/?basemap=b1&category=ned,nedsrc&title=3DEP%20View).The DEM data are used to delineate our study area into 23 sub-basins and stream networks.The landuse data in 2016 were obtained from the United States Department of Agriculture-National Agriculture Statistics Service (NASS) (https://nassgeodata. gmu.edu/CropScape/) at 30 m resolution, and then was resampled to 10 m resolution.The soil data were downloaded in the form of a 30 × 30 m raster layer from the soil survey geographic (SSURGO) database from the Natural Resources Conservation Service (NRCS), and also resampled into 10 m resolution.The daily maximum and minimum air temperatures, solar radiation, relative humidity, and wind speed of 20 weather stations in our study region were downloaded from the NRCS (http://wcc.sc.egov.usda.gov/nwcc/site?sitenum=2064)during the periods of 2000 to 2016 (Figure 1).The missing data for these 20 stations were generated using the weather generator model in SWAT.Table 1 shows the monthly irrigation amounts for soybean, corn, cotton, and rice in our study region.Irrigation data were obtained from Dakhlalla et al [41], which were calculated from the Mississippi Delta annual report summarized by the Yazoo Mississippi Delta Joint Water Management District (http://www.ymd.org/).

MODFLOW Input
The MODFLOW model domain is the entire Mississippi Delta region, an area of 17,866 km 2 .The observed groundwater level data from 2000 to 2016 were obtained from the Yazoo Mississippi Delta Joint Water Management District (http://www.ymd.org).Groundwater levels in each year were recorded in the beginning of April and October.The locations of 46 wells were displayed in Figure 1.The river stage data used for the river package were obtained from the U.S. Corps of Engineers website (www.rivergages.com).DEM and landuse data in SWAT were also taken as the input data for simulating groundwater recharge.
The left boundary of the BSRW was set as the specific head boundary, while the right boundary of the BSRW was set as general head boundary.The finite-difference grid in BSRW consists of 297 rows and 112 columns, with cell grid size of 1 × 1 km and with daily time steps and monthly stress periods.The parameters required by MODFLOW were hydraulic conductivity in the horizontal direction, vertical anisotropy, specific storage, specific yield, recharge multipliers, streambed conductance, and horizontal flow barriers.The specific values of these parameters can be found in the literature reported by Clark et al. [2].The aquifer in BSRW consists of Quaternary-age sand and gravel deposits overlying an erosional Tertiary-age surface.Although BSRW has multiple sources of recharge to the aquifer, including rainfall, leakage from the Mississippi river, and interflow from Bluff Hills, the recharge rate is 64 mm year −1 , which is 5% of the mean annual rainfall in this region [50].

Model Calibration and Validation Method
The model was calibrated from 2000 to 2008 and validated from 2009 to 2016 using USGS daily and monthly observed streamflow data at 3 gauge stations in our study region (Figure 1): Merigold (07288280), Sunflower (07288500), and Leland (07288500).Since streamflow is an important component of water balance, we calibrated and validated streamflow to assess the performance of the SWAT-MODFLOW model in BSRW.The model was also calibrated and validated using groundwater level data from 2000 and 2016.Each sub-basin has two wells, of which one was used for calibration and the other for validation.
The calibration and uncertainty analysis tool (SWAT-CUP) with the sequential uncertainty fitting (SUFI-2) algorithm was employed to perform calibration and validation for streamflow and groundwater levels in this study.Arnold et al. [47] reported that input parameters in SWAT are process-based and should be set in a reasonable uncertainty range.Based on the previous studies conducted by Dakhlalla et al. [41] and Jayakody et al. [51] in BSRW, five hydrological parameters were selected for calibration.They were the Manning's roughness coefficient for the main channel (CH_N2), soil conservation service curve number (CN2), surface runoff lag time (SURLAG), saturated hydraulic conductivity of the soil layers (SOL_K), and the threshold depth of water in the shallow aquifer (GWQMN).

Evaluation Criteria
Statistical parameters of the coefficient of determination (R 2 ), Nash-Sutcliffe efficiency (NSE), Root Mean Square Error (RMSE), and percentage bias (PBIAS) were used to evaluate the performance of the SWAT-MODFLOW model.The R 2 is a measure of strength between simulated and observed values.NSE is a measure to illustrate the relationship between observed and simulated values along a 1:1 line [52].RMSE is also a measure for assessing agreement between observed and simulated values.PBIAS measures whether the simulated results are systematically overestimated (negative PBIAS) or underestimated (positive PBIAS).Arnold et al. [47] and Coffey et al. [53] suggested that there were 20 statistic methods for evaluating model performance, such as R 2 , NSE, RMSE, nonparametric tests, t-test, objective functions, autocorrelation, and cross-correlation.By far, the most widely used statistics reported for calibration and validation are R 2 and NSE.Moriasi et al. [54] suggested that the model performance would be acceptable when NSE and R 2 > 0.5, and Tuppad et al. [55] reported that the model performance was rated as very good when NSE and R 2 > 0.75.Moriasi et al. [54] also reported that the model performance of streamflow can be rated as very good when PBIAS ≤ 10%, and rated as satisfactory when ±15% ≤ PBIAS ≤ ±25% [56].

Water Stress Index and Potential Availability of Surface Water Resources
In order to assess SWR potential for agriculture irrigation, the WSI (Water Stress Index) was calculated as follows: PSWR SWR − Environmental stream f low depth (10) where SWR is the sum of stream and pond water resources and PSWR is the potential availability of surface water resources.Some scientists suggested that stream water resources can be utilized in a more sustainable way if environmental streamflow regulation is established to prevent stream water resources from over-withdrawals [32,57].Ouyang [32] calculated the environmental streamflow of 0.8 m 3 s −1 , using the frequent-low flow and the conventional 7Q10 method in BSRW.We converted the environmental streamflow value to environmental streamflow depth for each sub-basin based on their areas, which ranged from 3.3 to 123.4 mm.Annual surface water demand is the irrigation amount required by the major crops in the Mississippi Delta (Table 1).In order to clearly show the whole process of the model application, the flowchart of methodology was presented in Figure 2.

Calibrated and Validated Parameters
The final ranges and reasonable values of the five parameters we selected for calibration and validation are listed in Table 2.The final calibrated CH_N2 is identical to the value of 0.3 reported by Dakhlalla et al. [41], which indicates that the surfaces of the channels in the Big Sunflower River Watershed are relatively rough [58].The fitted value of CN2 was −0.015, agreeing with the result of −0.017 reported in the literature [41].While the fitted value of SURLAG in the study of Dakhlalla et al. [41] was 1.4, our results had a higher value of 9.5, indicating that more surface runoff in BSRW should be released to the main channel due to the lower infiltration rate caused by the overlying clay and fine-grained material in the upper part of the aquifer [22].Compared to the lower value of 929 mm [41], the fitted value of GWQMN in our study was 1075 mm, which may have been caused by an absence of rainfall or irrigation return flow during the dry summer months [22].

Calibrated and Validated Parameters
The final ranges and reasonable values of the five parameters we selected for calibration and validation are listed in Table 2.The final calibrated CH_N2 is identical to the value of 0.3 reported by Dakhlalla et al. [41], which indicates that the surfaces of the channels in the Big Sunflower River Watershed are relatively rough [58].The fitted value of CN2 was −0.015, agreeing with the result of −0.017 reported in the literature [41].While the fitted value of SURLAG in the study of Dakhlalla et al. [41] was 1.4, our results had a higher value of 9.5, indicating that more surface runoff in BSRW should be released to the main channel due to the lower infiltration rate caused by the overlying clay and fine-grained material in the upper part of the aquifer [22].Compared to the lower value of 929 mm [41], the fitted value of GWQMN in our study was 1075 mm, which may have been caused by an absence of rainfall or irrigation return flow during the dry summer months [22].Simulated and observed monthly streamflow at three gauges were compared for calibration and validation in Figure 3.This indicated that the calibrated and validated results were acceptable in simulating streamflow in the study region.Tuppad et al. [55] reported that the simulated results in BSRW were in very good agreement with the observed values at Merigold during both calibration and validation periods, and at Sunflower during calibration (both R 2 and NSE > 0.75).Although not all the simulated streamflow values were close to the observed values, the R 2 and NSE values at the three gauge stations were > 0.5, indicating all the simulated results were acceptable during both calibration and validation periods [54].RMSE ranged from 10.19 to 16.82 m 3 s −1 , suggesting the simulated results were acceptable within our study region.PBIAS at the three gauge stations ranged from −22% to 25% during the calibration and validation periods, which suggests the least satisfactory model performance [54].
Water 2019, 11, x FOR PEER REVIEW 9 of 21 Simulated and observed monthly streamflow at three gauges were compared for calibration and validation in Figure 3.This indicated that the calibrated and validated results were acceptable in simulating streamflow in the study region.Tuppad et al. [55] reported that the simulated results in BSRW were in very good agreement with the observed values at Merigold during both calibration and validation periods, and at Sunflower during calibration (both R 2 and NSE > 0.75).Although not all the simulated streamflow values were close to the observed values, the R 2 and NSE values at the three gauge stations were > 0.5, indicating all the simulated results were acceptable during both calibration and validation periods [54].RMSE ranged from 10.19 to 16.82 m 3 s −1 , suggesting the simulated results were acceptable within our study region.PBIAS at the three gauge stations ranged from −22% to 25% during the calibration and validation periods, which suggests the least satisfactory model performance [54].The comparison statistics values of NSE, R 2 , RMSE, and PBIAS between observed and simulated daily streamflow values by SWAT and SWAT-MODFLOW were also displayed in Table 3.As compared with the performance of SWAT, SWAT-MODFLOW had higher R 2 and NSE but lower RMSE during both calibration and validation periods.SWAT-MODFLOW significantly improved at Merigold and slightly improved at Sunflower and Leland.This might have been due to inadequate spatial coverage of the precipitation inputs [47], as there was only one weather station for Leland, no weather station for Sunflower, and two weather stations for Merigold.Gupta et al. [59] reported that PBIAS values for streamflow might vary greatly among different hydrologic models and various auto-calibration methods, but this would not influence model evaluation if PBIAS values are within the reasonable range.Although PBIAS changed from negative to positive values for the Sunflower station in our study, the performance of SWAT-MODFLOW was rated as satisfactory (±15% ≤ PBIAS ≤ ±25%) [56].The comparison statistics values of NSE, R 2 , RMSE, and PBIAS between observed and simulated daily streamflow values by SWAT and SWAT-MODFLOW were also displayed in Table 3.As compared with the performance of SWAT, SWAT-MODFLOW had higher R 2 and NSE but lower RMSE during both calibration and validation periods.SWAT-MODFLOW significantly improved at Merigold and slightly improved at Sunflower and Leland.This might have been due to inadequate spatial coverage of the precipitation inputs [47], as there was only one weather station for Leland, no weather station for Sunflower, and two weather stations for Merigold.Gupta et al. [59] reported that PBIAS values for streamflow might vary greatly among different hydrologic models and various auto-calibration methods, but this would not influence model evaluation if PBIAS values are within the reasonable range.Although PBIAS changed from negative to positive values for the Sunflower station in our study, the performance of SWAT-MODFLOW was rated as satisfactory (±15% ≤ PBIAS ≤ ±25%) [56].The comparison of observed and simulated monthly streamflow at three gauges by SWAT and SWAT-MODFLOW were displayed in Figure 4.Although Dakhlalla et al. [41] reported that SWAT performed well in BSRW, the high degree of uncertainty has existed in simulating streamflow due to the overestimation of low streamflow.Figure 4 showed that SWAT overestimated low streamflow values (10 to 28 m 3 s −1 ) by as high as 93%, while SWAT-MODFLOW overestimated low streamflow by as low as 12%.Bejranon et al. [30] reported that the low streamflow may be overestimated due to the higher groundwater discharge simulated by the SWAT model.Barlow and Leake [50] also reported that even lower groundwater discharge would have significant impacts on streamflow.Therefore, MODFLOW was incorporated into SWAT, and the SWAT-MODFLOW was further calibrated and validated in our study, instead of re-calibrating and re-validating SWAT individually.Dakhlalla et al. [41] calibrated and validated SWAT in BSRW, but they did not simulate groundwater discharge in their study.We employed SWAT [41] and the SWAT-MODFLOW model to simulate groundwater discharge, and then made comparisons (Figure 5). Figure 5 displays comparisons of simulated monthly groundwater discharge between SWAT and SWAT-MODFLOW in BSRW.The monthly groundwater discharges simulated by SWAT were higher (1.9 to 4.5 mm) than those simulated by SWAT-MODFLOW (0.4 to 2.3 mm).Compared to SWAT, groundwater discharge can be well simulated by SWAT-MODFLOW, with consideration of groundwater pumping.As a result, the discrepancies between simulated and observed streamflow values were reduced [50].Barlow and Clark [22] reported that the average monthly groundwater recharge in BSRW from the Mississippi river and unsaturated zone were only 0.01 mm, due to the overlying clay and fine-grained material in the upper part of the aquifer.Therefore, the impact of groundwater recharge on groundwater discharge were neglected in BSRW.The comparison of observed and simulated monthly streamflow at three gauges by SWAT and SWAT-MODFLOW were displayed in Figure 4.Although Dakhlalla et al. [41] reported that SWAT performed well in BSRW, the high degree of uncertainty has existed in simulating streamflow due to the overestimation of low streamflow.Figure 4 showed that SWAT overestimated low streamflow values (10 to 28 m 3 s −1 ) by as high as 93%, while SWAT-MODFLOW overestimated low streamflow by as low as 12%.Bejranon et al. [30] reported that the low streamflow may be overestimated due to the higher groundwater discharge simulated by the SWAT model.Barlow and Leake [50] also reported that even lower groundwater discharge would have significant impacts on streamflow.Therefore, MODFLOW was incorporated into SWAT, and the SWAT-MODFLOW was further calibrated and validated in our study, instead of re-calibrating and re-validating SWAT individually.Dakhlalla et al. [41] calibrated and validated SWAT in BSRW, but they did not simulate groundwater discharge in their study.We employed SWAT [41] and the SWAT-MODFLOW model to simulate groundwater discharge, and then made comparisons (Figure 5). Figure 5 displays comparisons of simulated monthly groundwater discharge between SWAT and SWAT-MODFLOW in BSRW.The monthly groundwater discharges simulated by SWAT were higher (1.9 to 4.5 mm) than those simulated by SWAT-MODFLOW (0.4 to 2.3 mm).Compared to SWAT, groundwater discharge can be well simulated by SWAT-MODFLOW, with consideration of groundwater pumping.As a result, the discrepancies between simulated and observed streamflow values were reduced [50].Barlow and Clark [22] reported that the average monthly groundwater recharge in BSRW from the Mississippi river and unsaturated zone were only 0.01 mm, due to the overlying clay and fine-grained material in the upper part of the aquifer.Therefore, the impact of groundwater recharge on groundwater discharge were neglected in BSRW.The comparison of observed and simulated monthly streamflow at three gauges by SWAT and SWAT-MODFLOW were displayed in Figure 4.Although Dakhlalla et al. [41] reported that SWAT performed well in BSRW, the high degree of uncertainty has existed in simulating streamflow due to the overestimation of low streamflow.Figure 4 showed that SWAT overestimated low streamflow values (10 to 28 m 3 s −1 ) by as high as 93%, while SWAT-MODFLOW overestimated low streamflow by as low as 12%.Bejranon et al. [30] reported that the low streamflow may be overestimated due to the higher groundwater discharge simulated by the SWAT model.Barlow and Leake [50] also reported that even lower groundwater discharge would have significant impacts on streamflow.Therefore, MODFLOW was incorporated into SWAT, and the SWAT-MODFLOW was further calibrated and validated in our study, instead of re-calibrating and re-validating SWAT individually.Dakhlalla et al. [41] calibrated and validated SWAT in BSRW, but they did not simulate groundwater discharge in their study.We employed SWAT [41] and the SWAT-MODFLOW model to simulate groundwater discharge, and then made comparisons (Figure 5). Figure 5 displays comparisons of simulated monthly groundwater discharge between SWAT and SWAT-MODFLOW in BSRW.The monthly groundwater discharges simulated by SWAT were higher (1.9 to 4.5 mm) than those simulated by SWAT-MODFLOW (0.4 to 2.3 mm).Compared to SWAT, groundwater discharge can be well simulated by SWAT-MODFLOW, with consideration of groundwater pumping.As a result, the discrepancies between simulated and observed streamflow values were reduced [50].Barlow and Clark [22] reported that the average monthly groundwater recharge in BSRW from the Mississippi river and unsaturated zone were only 0.01 mm, due to the overlying clay and fine-grained material in the upper part of the aquifer.Therefore, the impact of groundwater recharge on groundwater discharge were neglected in BSRW.

Groundwater Level
The calibration and validation results of groundwater levels within the 46 wells are presented in Figure 6, indicating that groundwater levels in our study region were well simulated by SWAT-MODFLOW based on the 1:1 line fitting diagram.R 2 was 0.95 for calibration and 0.88 for validation, and NSE values were 0.99 and 0.93 for calibration and validation, respectively.The RMSE values for the calibration and validation periods was 0.51 and 0.59 m, demonstrating a good relationship between simulated and observed groundwater levels.The observed and simulated spatial distributions of groundwater levels averaged from April 2000 to October 2014 are displayed in Figure 7.Both the observed and simulated results showed a similar spatial distribution of groundwater levels in BSRW, ranging from 14 to 36 m.The groundwater levels were higher in the north of BSRW than in the south, since extensive irrigation amounts were used in the south where the groundwater storage loss has reached up to 1.8 × 10 8 m 3 per year [22,41].

Groundwater Level
The calibration and validation results of groundwater levels within the 46 wells are presented in Figure 6, indicating that groundwater levels in our study region were well simulated by SWAT-MODFLOW based on the 1:1 line fitting diagram.R 2 was 0.95 for calibration and 0.88 for validation, and NSE values were 0.99 and 0.93 for calibration and validation, respectively.The RMSE values for the calibration and validation periods was 0.51 and 0.59 m, demonstrating a good relationship between simulated and observed groundwater levels.The observed and simulated spatial distributions of groundwater levels averaged from April 2000 to October 2014 are displayed in Figure 7.Both the observed and simulated results showed a similar spatial distribution of groundwater levels in BSRW, ranging from 14 to 36 m.The groundwater levels were higher in the north of BSRW than in the south, since extensive irrigation amounts were used in the south where the groundwater storage loss has reached up to 1.8 × 10 8 m 3 per year [22,41].

Temporal and Spatial Variability of Stream and Pond Resources
Tables 4 and 5 show the average monthly total stream resources (MSR), average monthly pond resources (MPR), annual mean stream water resources (ASR), and pond water resources (APR) in BSRW during the periods from 2000 to 2016.Our results showed that BSRW had lower MSR values

Temporal and Spatial Variability of Stream and Pond Resources
Tables 4 and 5 show the average monthly total stream resources (MSR), average monthly pond resources (MPR), annual mean stream water resources (ASR), and pond water resources (APR) in BSRW during the periods from 2000 to 2016.Our results showed that BSRW had lower MSR values in growing seasons (MSR = 8.8 × 10 7 m 3 ) than in non-growing seasons (MSR = 11.0 × 10 7 m 3 ).This might be caused by the lower average monthly rainfall (91 mm) and greater evapotranspiration (ET) (71 mm) during growing seasons, and relatively greater rainfall (111 mm) and less ET (40 mm) during off-growing seasons in BSRW.Our results also showed that pond resources vary less than stream resources, since pond resources in our study region were less disturbed by human activities such as irrigation [24].In general, Vorosmarty et al. [60] and Stonestrom et al. [61] demonstrated that SWR was affected by precipitation and landuse.Sun and Ren [8] applied the simple linear regression analysis between SWR and rainfall, and found that SWR was more influenced by rainfall, since lower rainfall mainly affects crop water use, rather than SWR.They also indicated that SWR may be underestimated without considering the effect of landuse changes over time.However, simulation of streamflow resources still had uncertainties.Wang et al. [10] conducted similar modeling research in a semi-arid region of China, where the simulated average annual total SWR was 0.86 × 10 9 m 3 , which is lower than our results.Meanwhile, the annual average streamflow resources simulated by Sun and Ren [8] in a semi-arid region had a higher value of 4.3 × 10 9 m 3 compared to our results.This indicates that streamflow resources are not determined by rainfall only.As a result, the simulated stream resources in our study may also be underestimated, as Tao et al. [62] reported that Mississippi has experienced intensive landuse changes in recent years.Therefore, landuse changes should be taken into consideration when simulating the SWR. Figure 8 shows the spatial distributions of ASR in dry, normal, and wet years in BSRW.The empirical frequency analysis method was used to categorize 16 years into dry, normal, and wet years.In general, ASR in the entire BSRW ranged from 365 to 972 mm.ASR of 23 sub-basins in dry years, normal years, and wet years were 481, 605, and 726 mm, respectively.The proportion of sub-basins in BSRW with stream resources greater than 700 mm was 21% in dry years (ranging from 229 to 994 mm), while this increased to 35% in normal years (ranging from 296 to 1141 mm) and 57% in wet years (ranging from 554 to 991 mm).We focused more on values of 700 mm because only sub-basins 6, 8, 10, and 11 had stream resources >700 mm in dry years.The spatial distributions of pond resources for the 23 sub-basins in the BSRW are also displayed in Figure 9.The sub-basins with higher pond resources were found in the northwest of this region (>7.3 mm), since most of the on-farm ponds were built in this region [21].Although many studies demonstrated that stream water is the main source for agriculture irrigation [63][64][65], stream water might not be sustainably usable for irrigation in BSRW in the future, due to the stream water depletion caused by climate change and human activities [23].Instead, on-farm ponds should be constructed in BSRW as an alternative to collect rain and runoff.Although the annual and monthly pond volumes in BSRW were relatively small compared to stream water, pond water can play a significant role in irrigation, since Ouyang et al. [23] reported that 10,000 m 2 soybean fields could save approximately 542 m 3 groundwater each year if the ratio of pond size to irrigated soybean land was 1:18.Furthermore, irrigation water from streams are conveyed and disturbed by complex aqueduct systems, since most fields that need irrigation are often not close to streams.It is expensive to construct such water delivery systems.Ponds can be less expensive than aqueduct systems to construct adjacent to fields.irrigated soybean land was 1:18.Furthermore, irrigation water from streams are conveyed an sturbed by complex aqueduct systems, since most fields that need irrigation are often not close eams.It is expensive to construct such water delivery systems.Ponds can be less expensive tha ueduct systems to construct adjacent to fields.

Assessment of Surface Water Resources
Figure 10 displayed averaged annual WSI and PSWR in each sub-basin in dry, normal, and wet years.WSI in most of the sub-basins was less than 1.0, ranging from 0.6 to 2.1 for dry years (Figure 10b), 0.5 to 1.6 for normal years (Figure 10d), and 0.4 to 1.0 for wet years (Figure 10f), indicating that annual PSWR was higher than annual surface water demand.Therefore, SWR has the potential to be used for irrigation without adverse impacts on the environment and ecosystems.Girolimetto and Venturini [66] also estimated WSI in Oklahoma, USA, which has a similar climate as BSRW, and reported that the WSI in normal years ranged from 0.5 to 0.7.Kar and Kumar [67] reported that WSI was between 0.61 and 0.63 before application of irrigation.Our results are comparable with those published results, indicating that WSI can be used in our study area for the assessment of surface water resources.Figure 10a,c,e showed that PSWR was higher than irrigation demand in most of the sub-basins in dry, normal, and wet years.The values of PSWR-irrigation ranged from −170 to 651 mm in dry years, −13 to 712 mm in normal years, and 196 to 651 mm in wet years, indicating that there was net SWR available for irrigation in BSRW.Groundwater scarcity has been a serious problem in BSRW, and agriculture irrigation in this region cannot depend only on groundwater resources for a long time.However, surface water resources might be a supplementary irrigation source to mitigate the water resources scarcity in this region.

Assessment of Surface Water Resources
Figure 10 displayed averaged annual WSI and PSWR in each sub-basin in dry, normal, and wet years.WSI in most of the sub-basins was less than 1.0, ranging from 0.6 to 2.1 for dry years (Figure 10b), 0.5 to 1.6 for normal years (Figure 10d), and 0.4 to 1.0 for wet years (Figure 10f), indicating that annual PSWR was higher than annual surface water demand.Therefore, SWR has the potential to be used for irrigation without adverse impacts on the environment and ecosystems.Girolimetto and Venturini [66] also estimated WSI in Oklahoma, USA, which has a similar climate as BSRW, and reported that the WSI in normal years ranged from 0.5 to 0.7.Kar and Kumar [67] reported that WSI was between 0.61 and 0.63 before application of irrigation.Our results are comparable with those published results, indicating that WSI can be used in our study area for the assessment of surface water resources.Figure 10a,c,e showed that PSWR was higher than irrigation demand in most of the sub-basins in dry, normal, and wet years.The values of PSWR-irrigation ranged from −170 to 651 mm in dry years, −13 to 712 mm in normal years, and 196 to 651 mm in wet years, indicating that there was net SWR available for irrigation in BSRW.Groundwater scarcity has been a serious problem in BSRW, and agriculture irrigation in this region cannot depend only on groundwater resources for a long time.However, surface water resources might be a supplementary irrigation source to mitigate the water resources scarcity in this region.The coupled SWAT-MODFLOW model has advantages in estimating SWR.Firstly, SWAT-MODFLOW can be used to estimate SWR, even in an ungauged sub-basin.The Bureau of Hydrology [68] reported that SWR can be calculated by a subentry investigation method, but multiple gauge stations are required for this method to estimate stream water resources in each sub-basin, based on water use data and water balance equations.However, there are 23 sub-basins but only three different gauge stations in BSRW (Figure 1), and thus it is impossible to obtain reliable results for all the subbasins using the subentry investigation method.While SWAT-MODFLOW can estimate stream water resources in each sub-basin, even in an ungauged sub-basin [40], the subentry investigation method can only be applied in small anthropogenic-affected watersheds (Bureau of Hydrology 1986).However, BSRW was greatly influenced by human activities, such as heavy groundwater pumping.The calibrated and validated SWAT-MODFLOW can simulate SWR, even in regions similar to The coupled SWAT-MODFLOW model has advantages in estimating SWR.Firstly, SWAT-MODFLOW can be used to estimate SWR, even in an ungauged sub-basin.The Bureau of Hydrology [68] reported that SWR can be calculated by a subentry investigation method, but multiple gauge stations are required for this method to estimate stream water resources in each sub-basin, based on water use data and water balance equations.However, there are 23 sub-basins but only three different gauge stations in BSRW (Figure 1), and thus it is impossible to obtain reliable results for all the sub-basins using the subentry investigation method.While SWAT-MODFLOW can estimate stream water resources in each sub-basin, even in an ungauged sub-basin [40], the subentry investigation method can only be applied in small anthropogenic-affected watersheds (Bureau of Hydrology 1986).However, BSRW was greatly influenced by human activities, such as heavy groundwater pumping.The calibrated and validated SWAT-MODFLOW can simulate SWR, even in regions similar to BSRW, by taking groundwater discharge into consideration.Thus, it is more reliable to assess SWR using SWAT-MODFLOW.

Figure 1 .
Figure 1.Location of the study area, sub-basins, United States Geological Survey (USGS) gauge stations, weather stations, and observation wells.

Figure 1 .
Figure 1.Location of the study area, sub-basins, United States Geological Survey (USGS) gauge stations, weather stations, and observation wells.

Figure 2 .
Figure 2. Flowchart of methodology for assessing the surface water resources in the Big Sunflower River Watershed.

Figure 2 .
Figure 2. Flowchart of methodology for assessing the surface water resources in the Big Sunflower River Watershed.

Figure 3 .
Figure 3. Observed and simulated (by Soil and Water Assessment Tool-Modular Groundwater Flow model (SWAT-MODFLOW)) monthly streamflow values for calibration and validation at three United States Geological Survey (USGS) gauge stations (Merigold 07288280, Sunflower 07288500, and Leland 07288500) in the Big Sunflower River Watershed (07288280, 07288500, and 07288500 represent the gauge station identification numbers).

Figure 3 .
Figure 3. Observed and simulated (by Soil and Water Assessment Tool-Modular Groundwater Flow model (SWAT-MODFLOW)) monthly streamflow values for calibration and validation at three United States Geological Survey (USGS) gauge stations (Merigold 07288280, Sunflower 07288500, and Leland 07288500) in the Big Sunflower River Watershed (07288280, 07288500, and 07288500 represent the gauge station identification numbers).

Figure 4 .
Figure 4. Comparison of simulated (by SWAT and SWAT-MODFLOW) and observed monthly averages of daily streamflow at three United States Geological Survey (USGS) gauge stations (Merigold 07288280, Sunflower 07288500, and Leland 07288500) in the Big Sunflower River Watershed (07288280, 07288500, and 07288500 represent the gauge station identification numbers).

Figure 4 .
Figure 4. Comparison of simulated (by SWAT and SWAT-MODFLOW) and observed monthly averages of daily streamflow at three United States Geological Survey (USGS) gauge stations (Merigold 07288280, Sunflower 07288500, and Leland 07288500) in the Big Sunflower River Watershed (07288280, 07288500, and 07288500 represent the gauge station identification numbers).

Figure 4 .
Figure 4. Comparison of simulated (by SWAT and SWAT-MODFLOW) and observed monthly averages of daily streamflow at three United States Geological Survey (USGS) gauge stations (Merigold 07288280, Sunflower 07288500, and Leland 07288500) in the Big Sunflower River Watershed (07288280, 07288500, and 07288500 represent the gauge station identification numbers).

Figure 5 .
Figure 5. Comparisons of simulated monthly groundwater discharge to river between the SWAT model and the coupled SWAT-MODFLOW model for the entire Big Sunflower River Watershed.

Water 2019 , 21 Figure 5 .
Figure 5. Comparisons of simulated monthly groundwater discharge to river between the SWAT model and the coupled SWAT-MODFLOW model for the entire Big Sunflower River Watershed.

Figure 6 .
Figure 6.Comparisons of observed and simulated groundwater levels for the selected 46 observation wells from 2000 to 2016.

Figure 6 . 21 Figure 7 .
Figure 6.Comparisons of observed and simulated groundwater levels for the selected 46 observation wells from 2000 to 2016.Water 2019, 11, x FOR PEER REVIEW 2 of 21

Figure 7 .
Figure 7. Observed and simulated groundwater levels across the Big Sunflower River Watershed.

Figure 8 .
Figure 8. Spatial distributions of annual stream water resources in dry, normal, and wet years from 2000 to 2016 in the Big Sunflower River Watershed.

Figure 8 .
Figure 8. Spatial distributions of annual stream water resources in dry, normal, and wet years from 2000 to 2016 in the Big Sunflower River Watershed.

Figure 9 .
Figure 9. Spatial distributions of average annual pond water resources in the Big Sunflower River Watershed.

Figure 9 .
Figure 9. Spatial distributions of average annual pond water resources in the Big Sunflower River Watershed.

Figure 10 .
Figure 10.Simulated surface water stress index and potential availability of surface water resources in each sub-basin during the dry years, normal years, and wet years.

Figure 10 .
Figure 10.Simulated surface water stress index and potential availability of surface water resources in each sub-basin during the dry years, normal years, and wet years.

Table 1 .
Irrigation amounts for soybean, corn, cotton, and rice production in the Big Sunflower River Watershed.

Table 2 .
Calibrated results of the parameters based on the objective functions of streamflow and groundwater levels in the Big Sunflower River Watershed.
a CH_N2 is the Manning's roughness coefficient for the main channel; b CN2 is curve number; c SURLAG is the surface runoff lag time; d SOL_K is the saturated hydraulic conductivity of the soil layers; and e GWQMN is the threshold depth of water in the shallow aquifer for baseflow to occur.3.1.2.Streamflow

Table 2 .
Calibrated results of the parameters based on the objective functions of streamflow and groundwater levels in the Big Sunflower River Watershed.
a CH_N2 is the Manning's roughness coefficient for the main channel; b CN2 is curve number; c SURLAG is the surface runoff lag time; d SOL_K is the saturated hydraulic conductivity of the soil layers; and e GWQMN is the threshold depth of water in the shallow aquifer for baseflow to occur.

Table 3 .
Statistical comparison of observed and simulated (by both SWAT and SWAT-MODFLOW models) daily streamflow values at three United States Geological Survey (USGS) gauge stations (Merigold 07288280, Sunflower 07288500, and Leland 07288500) in the Big Sunflower River Watershed (07288280, 07288500, and 07288500 represent the gauge station identification numbers).
[41]lalla et al. (2016)t of determination between observed and simulated daily streamflow; b NSE is the Nash-Sutcliffe efficiency between observed and simulated daily streamflow; c RMSE is the root mean square error between observed and simulated daily streamflow; d PBIAS is the percentage bias between observed and simulated daily streamflow; and e The Statistical values for SWAT were obtained fromDakhlalla et al. (2016)[41].

Table 4 .
Simulated monthly average stream resources (MSR), monthly average pond water resources (MPR), and measured rainfall from 2000 to 2016 in the Big Sunflower River Watershed.

Table 5 .
Simulated annual stream resources (ASR), annual pond water resources, (APR) and measured rainfall in the Big Sunflower River Watershed.