Coupling of SUFI 2 and SWAT for Improving the Simulation of Streamflow in an Agricultural Watershed of South Dakota

Calibration and validation of process based hydrological models are two major processes while simulating the water balance components of watershed systems. However, these processes need a better understanding of the parameters which influence hydrologic processes within the system. In this study, we used SWAT model to simulate the stream flow for Skunk Creek (SK) watershed in South Dakota for the period from 1980-2000. Model calibration and validation were performed for both daily and monthly time periods using SUFI-2 within SWAT-CUP using 24 parameters selected from past available literature. Our calibration outputs for the period from 1987-1994 showed a good correlation between observed and model simulated values with NSE=0.56 and R2=0.70 for daily simulation. However, the model showed a better performance for monthly simulation with NSE and R2 values of 0.84 and 0.84 respectively. During validation period (1995-2000), the NSE and r2 values were 0.55 and 0.44, respectively for daily simulation and these statistical values were 0.76 and 0.77, respectively for monthly time step. Following calibration, the overall effect of each parameter used was ranked using global sensitivity function within SWAT-CUP. From the analysis, SOL_AWC was found to be the most sensitive parameter with absolute t-value of 17.50 and p-value of 0.00 to simulate the stream flow of the SK watershed. The CH_K2 was observed as the least sensitive parameter with t-statistic and p-value of 0.02 and 0.97, respectively. It was concluded from the study that coupling of the SWAT and SWAT-CUP made the calibration process quicker and reliable to simulate local hydrology within the watershed.


Introduction
Hydrologic models are widely used across the globe to simulate hydrologic processes including quality and quantity of stream flow in a basin [1]. It is highly labor, time, and cost intensive for maintaining gauging stations to collect water quality and quantity from a number of locations for long-term [2]. Therefore, hydrologic models play a significant role in simulating various hydrologic processes such as water quantity and quality, rainfall-runoff conceptualization, and sediment yield. There are ranges of models that are available for simulating longterm trends of hydrologic processes at small and larger watershed scales. Some of widely used hydrologic models include: Hydrological Simulation Program-Fortran (HSPF) [3], Systeme Hydrologique Europeen (SHE) [4], MIKE SHE [5], SHETRAN [6], Soil and Water Assessment Tool (SWAT) [7], Topographic Model (TOPMODEL) [8,9], and MOHID Land [10]. Of these, SWAT, a distributed process-based river basin continuous hydrologic model [7] is the most commonly and widely applied tool for simulating the management and climate change impacts on hydrologic processes at watershed scales.
SWAT model has been extensively used across the globe for assessing the long-term management [11], land use and climate change [1,[12][13][14][15][16] impacts on water quantity and quality at watershed scales. The development of SWAT is a continuation of United States Department of Agriculture-Agricultural Research Service (USDA-ARS) modeling experience that spans a period of over 30 years [17,18]. Accurate calibration of hydrologic models (including SWAT) is very important in assessing the water budget of a watershed for sustainable water resources management [19]. Improvement in calibration of models for better simulation of water quantity and quality has become a major topic of research among hydrologists. However, there are many uncertainties involved while working with hydrologic models that arise due to large spatial variability and multitude of input parameters [20]. These uncertainties over and under estimate the hydrologic processes, and can lead to wrong decisions [21]. Therefore, it is crucial to carefully carry the calibration and uncertainty analysis of hydrologic models for improved simulations [22][23][24].
To perform successful calibration and uncertainty analysis, there are various methods available that include such as parameter solution (PARASOL), Adaptive Clustering Covering (ACCO), general algorithm (GA), multi-start (M-Simplex), SWAT-CUP (including GLUE, SUFI-2, MCMC, PARASOL, and PSO), and uncertainty estimation based on local error and clustering (UNEEC). Among these, the Sequential Uncertainty Fitting 2 (SUFI-2) approach with SWAT Calibration Uncertainty Procedure (SWAT-CUP) is the most widely used approach to carry out parameterization, sensitivity analysis, uncertainty analysis, calibration, and validation of hydrologic parameters on both daily and monthly time-step [25]. The SUFI-2 is a semi-automated approach that makes the calibration process easier to carry within the realizable time bounds [26]. It has gained more interest among researchers for the reason that if it carried out manually, the incorporation of large number of parameters in a model can make the calibration process more complex and computationally extensive [27,28].
The first step in building any successful and reliable predictive hydrologic model is to carry Sensitivity Analysis (SA) which is helpful to identify and rank the parameters that have significant impact on specific model outputs of interest [29]. In SUFI-2, the techniques employed to perform the SA involves one-at-a time (OAT or Local Analysis) and Global Sensitivity Analysis methods. The OAT approach identifies the response from the output by sequentially varying each model parameter by a certain fraction while others are kept at their normal values [30][31][32]. This approach is less reliable as the parameter perturbation results in eccentricity from the nominal parameter value [33]. On the other hand, Global Sensitivity Analysis explores the entire range of the parameter values during model simulation process. In this method, all the parameters under consideration are simultaneously perturbed, allowing investigation of parameter's interaction and their influence on model outputs. This approach has the potential to capture the full range of model parameter values, and also to identify the interactions among parameters [34]. These two approaches yield different results [35], and therefore, it is always recommended to undergo comprehensive evaluation of parameter sensitivity [36].
There are number of factors involved in successful calibration of SWAT model such as parameter sensitivity, number of simulations, number of iterations, uncertainty associated among the parameters. Therefore, the specific objective of this study was to simulate the stream flow of an agricultural dominated watershed with minimum uncertainty among the parameters, and high reliability in predicting them using SWAT model coupled with SUFI-2 approach.

Study watershed
The present study was conducted for the Skunk Creek Watershed (SCW) located in the eastern part of South Dakota (SD) (Figure 1). The Skunk Creek is one of the main tributaries of the Big Sioux River. The Skunk Creek is 118 km permanent natural creek lying between 97.35-96.74° West and 43.45-44.13° North. The SCW covers an area of 1,606 km 2 and has its spread within six different counties of SD including Lake, Moody, Minnehaha, McCook, Lincoln, and Turner. The elevation of the watershed ranges between 351 and 574 m above mean sea level. The land use of the study watershed is dominated by cropland: corn (38%), soybean (26%), and pastures (14%). The c texture is dominated by silt clay loam with a few gravelly loams, loams, silt loams, loamy sand, and silt clay [37]. The average annual precipitation is about 590 mm, of which major proportion typically falls during the months from April to September. The average snowfall is 930 mm per year [37]. The average daily maximum and minimum temperature values are 15.5°C and 2.8°C, respectively.

SWAT model
The ArcGIS enabled SWAT model (version 2012) used for this study. The SWAT is an effective tool for evaluating the long-term impacts of management on hydrologic processes of diverse watershed systems [12,13,38,39]. In addition, this model also estimates the climate change impacts on plant growth, stream flow, and other responses such as total water yield, ET, snowmelt, and many others by taking into account the effects of increased atmospheric CO 2 concentrations on plant development and transpiration [40,41].
In SWAT, watershed hydrology is simulated in two major phases: land phase and routing phase [42]. The land phase controls the amount of water, sediment, nutrient, and pesticide loadings to the main channel in each sub-basin, whereas, the routing phase regulates the movement of these materials through the channel networks to the outlet of a watershed. The water balance equation that governs the land phase of hydrological cycle is given by the Equation 1 shown below as: where SW t is final soil water content (mm H 2 O), SW o is initial soil water content (mm H 2 O), t is time (days), R day is amount of precipitation on day i (mm H 2 O), Q surf is amount of surface runoff (mm H 2 O), E a is amount of evapotranspiration (mm H 2 O), W seep is amount of water entering the vadose zone from the soil profile (mm H 2 O), and Q gw is amount of return flow (mm H 2 O). The potential evapotranspiration in the model can be computed by three different approaches: Hargreaves [43], Priestley-Taylor [44], and Penman Monteith [45]. However, Penman-Monteith model was used for this study that incorporates energy and aerodynamic considerations [46]. To estimate the surface runoff, we used the Soil Conservation Service (SCS) Curve Number (CN) method [47] as shown below in Equation 2 which uses local land use, soil type, and antecedent moisture conditions.
where Q surf is accumulated runoff or rainfall excess (mm H 2 O), R day is rainfall depth for the day (mm H 2 O), I a is initial abstraction which includes surface storage, interception and infiltration prior to runoff (mm H 2 O), and S is the retention parameter (mm H 2 O). The I a was commonly approximated as 0.2S [47]. Therefore, the surface runoff will occur when R day >I a . The retention parameter S was computed as presented in the following Equation 3.
where CN is curve number for the day. The retention parameter which governs the value is spatially variable which is primarily based on local land use type and soil water content [42].
The SWAT has an embedded tool to reduce uncertainty through analyzing sensitivity of the parameters. It uses a method which is a dataset used by hydrologists across the globe for simulating hydrologic processes [13,51,52]. Temperature and precipitation data were derived from five weather stations located within the basin and named as Stations A, B, C, D, and E as shown in the Figure 1 and Table 1. The remaining meteorological inputs were automatically generated within the SWAT using daily precipitation and temperature data. Simulated stream discharge outputs from the generated model were compared with stream discharge taken from the USGS site no 06481500 located at Sioux Falls, SD for the period from 1987-2000. Detailed description of hydro-meteorological stations chosen for the study is presented in Table 1.

SWAT model set-up and calibration
After incorporating all the data inputs, the SWAT model was built with the total of 7 sub-basins and 364 Hydrological Response Units (HRUs). First, the model was set to run for 1980-2000 with the initial 7 years (1980-86) as model warm-up period that allows the model to stabilize for further simulations. Then, the consequent phases were established as calibration (1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994) and validation (1995)(1996)(1997)(1998)(1999)(2000) periods. This was done with care so as to ensure that both the periods should nearly have same water balance values. Calibration and validation of the model were done using the SUFI-2 approach within SWAT-CUP considering 24 key hydrologic parameters based on detailed previous literature sources [53,54]. Then, each parameter was set to a default lower and upper values as suggested by the SWAT expert group [42]. Finally, the best fitted parameter values obtained from SWAT-CUP were incorporated into the SWAT database for stream discharge simulations at both daily and monthly time steps. The model performance was evaluated using Nash-Sutcliffe efficiency (NSE) [55,56] coefficient of determination (R 2 ), and percentage of bias (PBIAS) [57] indices as presented in the following Equations 5-7.
where, Y obs is the measured data, Y sim is the model simulation output, and Y o mean and Y s mean is the mean of observed data and simulated data of stream flow, and m and s stand for measured and simulated, i is the i th measured or simulated data. Figure 2 demonstrates the schematic flow of processes those are involved while operating the SUFI-2. First, the objective function combination of Latin hypercube and one factor-at-a time sampling that allows global sensitivity analysis for various parameters with only limited number of model runs [21]. Therefore, it necessitates to run the model multiple times for better accuracy of the model simulation outputs. In our study, both the calibration and sensitivity analyses were carried out using semi-automated stochastic model, the SWAT-CUP. We performed the global sensitivity analysis where the parameter sensitivities were determined by calculating multiple regression system as presented in the Equation 4.

SUFI-2 algorithm
where g is the objective function, α and β are the variables and b i is parameter. A t-test was then used to identify the relative significance of each parameter b i through the application of inverse optimization approach.

Data source and analysis
Data on topography, land use, soil, climate, and stream discharge required for this study were compiled from different sources. The Digital Elevation Model (DEM) of 10 m by 10 m resolution was derived from Geospatial Data Gateway (GDG) to use as topographic data of the study basin. The DEM is required for delineating a watershed into sub-basins and then into smallest representative unit of the watershed, the Hydrologic Response Units (HRUs) based on specific land use, soil, and slope characteristic features. Land use map for this study was obtained in the form of Cropland Data Layer (CDL) which is created by USDA, National Agricultural Statistics Service (NASS). A CDL is a raster, geo-referenced, crop specific land cover dataset created annually for continental US using moderate resolution (30 m and 56 m) satellite imagery and extensive agricultural ground truth. The soil map was obtained from Soil Survey Geographic Data (SSURGO) collected by National Cooperative Soil Survey (NCSS), USDA, and Natural Resources Conservation Services (NRCS) with the scales ranging from 1: 12,000 to 1: 63,360 [48]. All these spatial datasets were set to the projection of WGS-1984 UTM Zone 14N using ArcGIS (version 10.0) for further simulations.
Similarly, the SWAT requires meteorological daily data inputs such as precipitation, maximum and minimum air temperature, wind speed, relative humidity, and solar radiation for hydrologic simulations. For the current study, the temperature and precipitation data were extracted from Daily Surface Weather and Climatological Summaries (DAYMET) Single Point Data Extraction (SPDE) [49,50] for five different spatial locations as Stations A, B, C, D, and E, and presented in the Figure 1 and Table 1. This dataset is available on a daily time scale with the resolution of 1 km × 1 km, and supported by the funding from National Aeronautics and Space Administration (NASA) through Earth Science Data and Information System (ESDIS) and terrestrial Ecosystem Program. This is commonly used meteorological (here to maximize the NSE) was defined. Then, the algorithm optimizes absolute maximum and minimum range of the parameters. It was assumed that all the parameters were uniformly distributed within the region bounded by these range of values. The SUFI-2 algorithm was designed in a manner that no automated optimization routine can replace the insight from physical understanding and knowledge of the effects of parameters on system response. Further, first round of Latin hypercube sampling [58] took place, where initial uncertainty ranges were next assigned to parameters, making parameter combinations. Following the Gaussian Newton method and neglecting higher order derivatives, equivalent of Hessian matrix was then calculated. Based on Cramer Rao theorem [59], an estimate of the lower bound of parameter covariance matrix was computed. It was observed that there was quite small correlation between any two parameters; as all parameters were allowed to change. This was because in SUFI-2, all the parameters are kept constant while a change is brought in one only at a time [60].
Global Sensitivity Analysis (GSA) was used to carry out sensitivity analysis of the parameters chosen for the calibration process. In this process, multiple regression system regress Latin hypercube generated parameters against objective function to determine sensitivity of the parameters. As statistical measurements, t-stat and p-value were used. A t-stat is the coefficient of a parameter divided by its standard error and is a measure of the precision with which the regression coefficient is measured. Therefore, the parameter is sensitive when the coefficient is larger than the standard error [60]. A p-value was determined from student's t-distribution table with the values obtained for t-stat for a parameter; where a lower p-value suggests higher sensitivity of the parameter, and vice-versa [60]. The overall uncertainty in the output was computed by 95 Percent Prediction Uncertainty (95 PPU) and dotty plots for each parameter. This helped in determining the new ranges and best fitted values that were applied for further iterations to maximize the objective function. The 95 PPU was calculated at 2.5 and 97.5% levels of the cumulative distribution of an output variable obtained through Latin hypercube sampling. Two indices were used to quantify the goodness of calibration/uncertainty performance: the p-factor, which is the percentage of data bracketed by the 95 PPU band (ideal value should approach closer to 1) and the r-factor, which is an average width of the band divided by standard deviation of the corresponding measured variable (ideal value should be close to 0). To minimize uncertainties and maximize the objective function, the number of sampling round was increased with the set of new parameter ranges. However, care was taken for using the physical acceptable parameter values for this process.

Climate and stream discharge data, and model performance
The climate data collected for 1987-2000 showed that the higher precipitation was received from April to September with the maximum value of 105 mm in June, and minimum value of 12.4 mm in December (Figure 4a). Nearly, 40% of total precipitation was received during summer season of the study period. The average monthly maximum stream flow (293 m 3 /s) was observed in the month of April. A total of about 49% of cumulative stream flow cumulated on seasonal basis was observed during the spring season (Figure 4b). The higher stream flow during spring season can be attributed to the fact that the snow melt water dominantly contributes during this season [61].
The model evaluation coefficients of the simulated daily stream discharge are presented in Table 2. The R 2 and NSE values for calibration (1987-94) were 0.70 and 0.56, respectively, whereas these values were 0.55 and 0.44 for the validation period (1995)(1996)(1997)(1998)(1999)(2000). The model showed better performance for the monthly simulations with R 2 and NSE values as 0.84 and 0.84, respectively, for the calibration, and 0.77 and 0.76 for the validation period. The measured and simulated stream discharge values were represented in the hydrographs shown in Figure 3. The model outputs revealed that the model is satisfactory in simulating the stream flow for the study watershed [55].

PPU plots
The hydrographs of 95 PPU plots derived from two different iterations (500 and 2000 simulations) are presented in the Figure 5. The p-factor and r-factor values were 0.57 and 0.90, respectively, for the 500 simulations. However, the factor values were 0.65 and 0.94 for the 2000 simulations obtained with the same set of input parameters. The combination of p-factor and r-factor together indicates the strength of model calibration and uncertainty assessment, as these are intimately linked [60]. The data showed that by keeping all the conditions same, the different number of model simulations produced different results for uncertainty indices, and therefore, affecting the goodness-of-fit. This approach of comparing two different simulations helped for better assessment of model uncertainties. The p-factor in either case was found to achieve desirable values [15,24,62]. The large r-factor values may be attributed to insufficient accounting of agricultural water use in the model [62] and inter-basin water transfer [63]. However, a reexamination of the conceptual model was needed if the p-factor value does not lie within the proposed acceptable range [64].

Parameter sensitivity
The most sensitive input hydrologic parameters identified on the basis of global sensitivity analysis and used for stream discharge simulations are presented in Table 3. We found different ranking of the same parameter with different simulation numbers, indicating the stochastic nature of the SWAT-CUP [65]. However, for both the simulations, the SOL_AWC and CN2 were observed as the most sensitive parameters to influence stream discharge of the basin. This was due to the reason that SOL_AWC represents the soil moisture characteristics, and plays an important role in evapotranspiration to influence the surface runoff [66]. Similar to these findings, Kannan et al. [67] also suggested that the soil water capacity correlated with various water balance components. In addition, SOL_AWC is considered to be directly proportional to the ability of soil to hold water affecting the stream flow. Similarly, the higher sensitivity of the CN2 is attributed to the higher influence of runoff generation in the basin [68]. However, the SMFMN which was moderately sensitive for 500 simulations was observed the least sensitive parameter for 2000 simulations. The ranking of most sensitive parameters observed in this study was also supported by the findings of Faramarzi et al. [62,67].
It was found that the parameter rankings impact the value of objective function when we altered the number of model simulations; though the change was small. However, we found larger changes when the simulation change was accompanied with the increasing iterations. In the latter case, we obtained the narrower parameter range with the maximum objective function values. The new fitted values were also observed beyond the limit of meaningful physical range of SWAT parameters. Therefore, these assessments demonstrated that awareness about meaningful physical range of hydrologic parameters is crucial while working with semi-automated stochastic calibration tool to monitor water balance components of a watershed system.

Dotty plots
The dotty plots mapped the model parameter and objective function values to help in identifying the relative sensitivity associated with each parameter influencing the objective function [69]. The plots obtained for 500 and 2000 simulations showed that there was trend followed by the points (corresponding to every simulation in an iteration) in case of SOL_AWC and CN2 that was influencing the objective function unlike with other parameters (Figures 6 and 7). Therefore, the SOL_AWC and CN2 were observed to be the most sensitive parameter in both the cases. This process also helped to define the new range of values for the parameters that further decided the best fitted values to work for next iteration level. With change in number of simulations and iterations, the range and best fitted values of the parameters were also changed   The dotty plots derived from SUFI-2 for 500 simulations are shown for the most and least sensitive parameters used in this study (Note: the above plots (a to d) are for the most sensitive and the below plots (e to h) are for the least sensitive parameters).

Statistics
Pre     (Table 4). For 500 simulations, the best fitted value for CN2 (where maximum NS approximately approached to 1) was observed to be approximately -0.1. However, as the simulations increased to higher number (2000 in our case), the value shifted to around 0.02. The relative sensitivity of each parameter was estimated by observing the impact on an objective function using the dotty plots. It was concluded that if the points on dotty plots are scattered or haphazard; the sensitivity is low for the parameter and if the points do follow a trend, the sensitivity is higher.

Superiority in using SUFI-2 over other algorithms
The SUFI-2 tool involves stochastic calibration, where the errors and uncertainties in model are recognized and expressed as ranges accounting for all driving variables, conceptual model, parameters and measured data [60]. The provision of inclusion of large number of parameters representing different processes in the objective function, in SUFI-2, helps to make the model result enveloping most of the observations well [70]. While simulating runoff and sediment modeling using SWAT in Gumera catchment in Ethiopia, used two approaches for calibration. They were fully automated Parameter Solution (ParaSol) and semi-automated Sequential Uncertainty Fitting 2 (SUFI-2) for the period from 2003 to 2006. They carried out calibration using 13 runoff producing parameters and concluded that SUFI-2 was more flexible to work with than Parasol and yielded higher values for coefficient of determination and NSE coefficient. Luo et al. [71] in their research on carrying calibration and uncertainty analysis of Japanese river Catchment using SWAT model using GLUE and SUFI-2 algorithm concluded that though the calibration results improved using GLUE approach but the processing time of GLUE approach is longer than SUFI2. In contrast to this, research carried by Singh et al. [54] showed that coefficient of determination and NSE was found to get reduced while using GLUE algorithm on daily time step. Even results inferred from this study also improved using SUFI-2, which can be seen in having more improved values for statistical evaluating parameters in case of calibration when compared with pre-calibration. This was possible because it reduced the uncertainty among the parameters and helped to build a highly reliable model with higher objective function of simulating the stream discharge relative close to observed values. This makes SUFI-2 to stand out from most of algorithms used for carrying calibration and uncertainty analysis of a hydrologic model because, in SUFI-2, uncertainty among the parameters accounts for various sources like uncertainties in input data, conceptual model, and among parameters itself because desegregation of the error into its source components is complex, especially in cases pertaining to hydrologic modeling where the models are nonlinear and different sources of error may interact to produce measured deviation [72,73].

Summary and Conclusions
The ArcGIS enabled SWAT (ArcSWAT) model was used to simulate stream flow from the Skunk Creek watershed using 24 different key hydrologic parameters. Parameterization, and uncertainty and sensitivity analyses were carried out using the SUFI-2 approach within  the SWAT-CUP. Our simulation results showed a reasonable accuracy between measured and model simulated stream flow values. The SWAT-CUP improved the stream flow simulations, and reducing uncertainty among the parameters. It was observed that due to the inclusion of larger confidential interval in less sensitive parameters, the uncertainty reduction among these parameters took more time than more sensitive parameters. Moreover, during parameterization process, awareness of physical meaningful range of parameters chosen for calibration led to better simulation results. It was also observed that in order to maximize the objective function, optimum number of iterations and simulations should be performed, else the best fitted value for the parameters may go beyond acceptable range. Finally, semi-automated stochastic model, the SWAT-CUP improved the SWAT simulations of stream flow with the meaningful physical acceptable range of the key hydrologic parameters and higher statistical evaluating parameters depicting more reliability of simulated results.