1. Introduction
A critical repercussion of agroecosystem growth is the deterioration of water quality from nonpoint source pollution derived from locations that cannot be ascribed spatially in the landscape [
1]. Transport of nitrogen (N) and phosphorous (P) [
2,
3,
4], sediment [
5] and pesticides [
6] is responsible for a range of adverse impacts on human and environmental conditions [
7], constituting the primary source for impairment of rivers and streams in the United States [
8]. Given the difficulty in identifying and controlling nonpoint source pollution, traditional monitoring and regulatory strategies, such as point-based monitoring, are difficult. Further, attempts to test possible management solutions are encumbered by cost, uncertainty and potential lack of participation.
Hydrologic models provide the capability to represent and predict natural processes, especially where observations are limited, thereby aiding in scenario analysis of system response to the implementation of nonpoint source control measures. Accordingly, there is a need to understand the accuracy and fidelity of these models in representing known processes to improve the predictive certainty for a range of intended applications. Uncertainty in hydrologic models is attributed to a combination of measurement error, systematic error, natural variation, model uncertainty, subjective judgement and inherent randomness [
9,
10,
11]. Therefore, efforts to improve the capability of hydrologic models to effectively capture and mathematically represent natural processes in a physically meaningful way necessitates continuing adaptation and refinement of methodology.
River basins and watersheds are common boundaries of most hydrologic models focused on nonpoint source pollution, as they provide an integrated response of natural and anthropogenic processes [
12]. Reviews of several hydrologic models underscore differences in the choice of discretization and model complexity to capture system response [
13,
14,
15]. Mathematical representation of physical processes depends on the level of discretization used and should reflect the scale of interest of the intended application. Discretization can be defined either temporally or spatially and seeks to relate the modeled scale to observation and process scales [
16]. Changes in scale may require development of new mathematical models and approaches. Based on data availability, hydrologic models generally adopt a daily time step, although effects of temporal scale are not well understood [
17]. Simple, lumped approaches often are appropriate spatially and can produce comparable results to more detailed, complex approaches depending upon the application [
18,
19]. However, many simple hydrologic models lack necessary physical processes and management scenarios essential for modeling the associated transformation and transport of nonpoint source pollution [
12]. Balancing the need for explicitly modeling spatial processes with data availability and computational cost is thus an ongoing research issue [
20]. Calibration of the model is of equal importance, as it may contribute to improved results or may mask deficiencies in simulations [
21].
The Soil and Water Assessment Tool (SWAT) is a continuous-time, distributed-parameter model, capable of simulating hydrology and water quality at the watershed scale [
22]. According to Wellen et al. [
15], SWAT is the most commonly-used model for distributed watershed nutrient modeling and has been used in a wide range of applications [
23]. SWAT predominantly relies upon discretizing landscapes based on common soil, land use and slope characteristics, known as hydrologic response units (HRUs). While the HRU approach provides a simple, computationally-efficient framework, processes modeled on HRUs are lumped and therefore spatially disconnected, as they are routed directly to subbasin outlets. This was identified as a key weakness of the model by Gassman [
23] and Douglas-Mankin et al. [
24], among others. This lack of definition of landscape position makes implementation of spatially-targeted management difficult to incorporate in the model.
Although readily accessible in the SWAT model framework, few studies have evaluated modified discretization or routing beyond standard HRU routing through subbasin networks [
20,
25,
26,
27,
28]. To overcome the spatial limitations of the HRU approach, a grid-based version of the SWAT model, SWATgrid [
29], was developed to perform landscape simulation on a regularized grid [
30] and employ a modified landscape routing algorithm developed by Volk et al. [
31] and Arnold et al. [
20] to spatially connect grids during routing. While other gridded hydrologic models exist, few incorporate detailed management, plant growth and water quality algorithms intrinsic to SWAT.
Rathjens et al. [
29] found the approach was able to simulate streamflow reasonably well. However, SWATgrid remains largely untested, with little understanding of the impact of user-defined model spatial resolution. Previous studies have been conducted at spatial resolutions of 100 m or greater and have not considered watersheds with artificially drained by subsurface tiles ubiquitous to poorly-drained agricultural landscape, such as the U.S. Midwestern agricultural systems [
29,
30,
32]. Unfortunately, computational overhead is significant for SWATgrid, effectively precluding direct calibration and necessitating strategies such as parameter transfer from HRU-based SWAT models [
33]. Parameters have been shown to exhibit scale dependency [
34], so correct characterization of their distributions at various resolutions is critical prior to transfer.
Use of subbasins is the most commonly-employed spatial discretization approach in watershed modeling [
15]; however, input data are generally initially available in gridded raster format prior to reconceptualization. Typical spatial input data include a digital elevation model (DEM), land use and soils. For hydrologic modeling, modifying the resolution of a DEM is considered the most impactful as it affects watershed delineation, slopes and channel length, while land use and soil classification primarily modify flow partitioning [
35]. Several studies have assessed the impact of DEM resolution on hydrology and water quality predictions with SWAT. Runoff, sediment and nutrients were found to decrease with coarser resolution by Cho and Lee [
36], Chaubey et al. [
37] and Ghaffari et al. [
38], due to smoothing of average watershed slope. Conversely, Chaplot [
39], Lin et al. [
40], Shen et al. [
41] and Zhang et al. [
42] identified minimal change in flow while water quality variables remained impacted.
Research focused on the impacts of resolution on SWAT has not yet extended beyond the standard HRU-based approach. The overall goal of this study was to test the efficacy of the SWATgrid approach while gaining insight into its sensitivity to user-defined spatial input resolution. Accordingly, three objectives were considered: (1) identifying a target resolution range for SWATgrid applications to maximize predictive accuracy while minimizing computation time; (2) comparing SWATgrid simulations to HRU simulations; and (3) identifying differences in calibrated parameter ranges for varying resolutions. Future SWAT applications, such as identifying critical source areas or testing best management practices, may require greater spatial detail [
43] than afforded by subbasin discretization. The grid-based approach offers one possibility to capture such spatial details in order to develop effective land management options for reducing losses of nonpoint source pollutants.
2. Materials and Methods
To accomplish the objectives of this study, several models were created using multiple input data resolutions for both the SWAT and SWATgrid. SWAT HRU models used the standard SWAT routing methods, while SWATgrid models (GRID) implemented both standard and landscape routing. This enabled a more rigorous comparison of HRU- and GRID-based methodologies by isolating both the effects of the landscape routing and, in turn, the gridded spatial representation of the landscape. Therefore, for each resolution considered, three separate model simulations were assessed: (1) HRU-based spatial representation with standard routing (HRU); (2) GRID-based spatial representation with landscape routing (GRID-LAND or LAND); and (3) GRID-based spatial representation with standard routing (GRID-STD or STD).
2.1. Study Location
The Cedar Creek Watershed (CCW) is located in the southwestern region of the St. Joseph River Basin in northeastern Indiana. Covering approximately 707 km
2, the CCW drains two HUC (Hydrologic Unit Code) 10 subwatersheds (0410000306 and 0410000307). A monitoring network was established in the northern CCW subbasin in 2002 as part of the United States Department of Agriculture’s Agricultural Research Service’s (USDA-ARS) Source Water Protection Project (SWPP) and Conservation Effect Assessment Project (CEAP), and it is shown in
Figure 1. The topography is characterized by slopes of 2% to 10% with an overall gradient towards the southeastern outlet. Soils were formed from compacted glacial till and fluvial materials with silt loam, silty clay loam and clay loam textures. The CCW is primarily agricultural, with corn and soybean rotational planting, although winter wheat, oats, alfalfa and pasture are also present [
44,
45,
46]. Due to oversaturation of soils and the prevalence of agricultural land, tile drainage is commonly employed throughout the watershed. Climatic conditions for the CCW include historical average temperatures from −1 to 28 °C with 940 mm of annual precipitation. For this study, two weather stations were used, one in the lower center of the watershed and a second outside the watershed boundary to the north. Station data were obtained from the USDA SWAT climate dataset <
ars.usda.gov/Research/docs.htm?docid=19388>.
2.2. The SWAT Model
The SWAT model is a continuous-time, semi-distributed, process-based river basin model [
22]. Inclusion of hydrology, weather, crop growth, soil, nutrient, pesticide and agricultural management components allows SWAT to predict the impact of various agricultural management scenarios on hydrology and water quality. Fundamentally, the model is hydrologically driven, based upon the water balance of an individual landscape unit:
where SW is the soil water content, R is precipitation, Q is surface runoff, ET is evapotranspiration, P is percolation and QR is return flow in units of mm for each day, i. Each landscape unit consists of a possible snow layer, a soil profile divided into multiple layers and a shallow and deep aquifer.
Standard SWAT Routing
Hydrologic processes are first simulated at the landscape (HRU) level and subsequently routed at the subbasin level where each subbasin has an associated reach. Landscape discretization into HRUs provides a simple framework for simulation, but also represents a significant limitation of the model. Spatial connectivity and interactions among HRUs are ignored, instead aggregating the cumulative output of each spatially discontinuous HRU group at the subwatershed outlet, which is then directly routed to the channel. Water balance equations are applied at the HRU level. For complete documentation, see Neitsch et al. [
47].
2.3. SWATgrid
SWATgrid differs from HRU-based SWAT models in two key ways: (1) use of gridded, discretized landscape units and (2) linking these landscape unit grids and routing surface runoff, lateral subsurface and shallow groundwater flow. Gridded landscape units effectively function as both an HRU and subbasin and take advantage of standard SWAT landscape simulation algorithms.
As part of SWATgrid, Rathjens and Oppelt [
30] developed an approach that utilizes the Topographic Parameterization (TOPAZ [
48]) tool to generate gridded landscape units. The grids allow for spatial connectivity during the routing phase as they are at a defined location, as opposed to HRUs, which are not. Landscape routing is then implemented through these grids adapted from the approach used by Volk et al. [
31] and Arnold et al. [
20]. Spatial position for each grid is calculated from a modified topographic index adapted from the TOPMODEL (Topography Based Hydrologic Model) framework:
where
Ai is the upslope contributing area per unit contour length (m),
βi is the topographic slope angle,
Ki is mean saturated hydraulic conductivity (m/day) and
Zi is the soil depth (m) for each grid cell
i. As the index approaches one, it approximates a valley or floodplain and near zero, a divide. Effectively, the index partitions flow on a given grid cell between channelized and landscape flow. In SWATgrid, the flow separation index is normalized for either abrupt or gradual channel heads and is adjusted by the drainage density of the watershed to capture channel head locations more accurately. Thus, when the adjusted flow separation index is zero, there is only landscape flow, and conversely when it is one, the flow is completely channelized. For an in-depth description of the SWATgrid, see Rathjens et al. [
29].
2.4. Baseline Model
An initial 30-m spatial resolution model was developed and calibrated to serve as a baseline for comparison of model simulations across resolutions. The 30-m resolution was nominally selected as the finest scale considered, given that it is common for most spatial datasets used in SWAT modeling (e.g., National Elevation Data (NED), National Land Cover Database (NLCD)). The watershed was delineated using a 30-m digital elevation model (DEM) from the NED. Land use and soil properties were obtained from the Crop Data Layer (CDL, 2010) from the National Agricultural Statistics Service (NASS) at 30 m and State Soil Geographic (STATSGO) Dataset at 250 m, respectively. STATSGO data were selected given their native grid structure provided by the SWAT soils database and relatively fewer soil classes than other soil databases (e.g., SSURGO—Soil Survey Geographic Database) , which facilitate ease of writing the manually intensive SWATgrid input files. STATSGO data were resampled by the nearest neighbor approach to 30 m prior to rescaling at coarser resolutions. Various CDL agricultural land use categories constituting less than 1% of the watershed area (e.g., pumpkins) were reclassified to general agriculture.
The model was forced with precipitation and temperature data from 1990 to 2010 from the United States Department of Agriculture (USDA) Agricultural Research Service (ARS) SWAT U.S. Climate Dataset stations C120200 and C123207. Discharge data for calibration and validation at the watershed outlet were obtained from the United States Geological Survey (USGS) Water Data Server for Station Number 04180000, Cedar Creek near Cedarville, IN, over the same period. Geographic input datasets are shown in
Figure 1.
Delineation and initial model parameterization were achieved using the ArcSWAT interface with the input datasets previously described. Based on previous experience in modeling this watershed, 1000 ha were used to define the minimum area necessary for channel formation. As SWATgrid evaluates all unique possible combinations of land use and soil types in the model, the HRU threshold for soils and land use was set to 0%. Watershed delineation produced 43 subbasins consisting of 1082 HRUs.
Several initial parameters values were specified in accordance with previous research in similar agricultural watersheds, as well as using expert knowledge, including previously published SWAT modeling studies from this area. Management parameters were specified to fit dominant expected scenarios for a corn and soybean rotation, as shown in
Table 1 (see [
49] for details). Additionally, given the prominence of tile drained agricultural landscapes in the Midwest [
50], tile drainage was implemented in each HRU comprised of a corn-soybean rotation located on poorly- or very poorly-drained soils (soil hydrologic Groups C and D, respectively). The modifications required for tile drains are summarized in
Table 2 (see [
51] for details). Additionally, to allow water to drain into the tiles, the curve number parameters were adjusted by assuming that each C or D class soil would behave as one class above its level. Finally, the SWAT baseflow recession constant (ALPHA_BF) was estimated using the SWAT Baseflow Filter Program [
52,
53] and set to 0.0321.
Based on previous SWAT research in Northern Indiana [
54], six basin level parameters, summarized in
Table 3, were specified for model calibration. By selecting basin level, rather than distributed parameters, calibration was applied in a consistent way to all models, regardless of varying HRU relative area between different resolutions and model types. A multi-algorithm, genetically adaptive multi-objective method (AMALGAM) [
55] was adapted for use with SWAT input and output data structure. Calibration was performed with monthly simulations between 1990 and 2010. The first four years were used for model warm up, while 1994 to 2003 and 2004 to 2010 were used for calibration and validation. Five thousand simulations were sampled and evaluated over 50 generations. Parameter optimization was constrained by maximizing Nash–Sutcliffe efficiency (NSE [
56]).
2.5. Resampling Approach
The SWATgrid model was run at spatial resolutions of 30, 60, 90, 150, 250, 500 and 1000 m. Resolution of the input data has numerous implications in watershed simulations. However, the primary focus here is to select a scale that captures a realistic system response while achieving viable computational time, rather than extensively evaluate the implications of the spatial detail of input data. Land use, soil and topography data were resampled from 30-m base data for each set of conditions. Nearest neighbor resampling was used to preserve the “hard” classes of land use and soils. Similar to other research [
57,
58], resampling was achieved by bicubic convolution to provide a smoothly interpolated DEM that retains peaks and depressions in the dataset, which is critical in watershed delineation and channel definition. Resampled soil and land use data in the Cedar Creek Watershed had a <2% difference in proportional watershed area between the finest and coarsest resolutions for all classes, even though actual area changed across resolutions.
2.6. Simulation Effects
Input data at 30-m resolution were used to create a baseline SWAT model and calibrated as previously described. Parameters from this baseline model were then transferred to lower resolution HRU and GRID models. By maintaining a constant parameter set applied at the basin level, the effects of resolution and routing type were better isolated. Each HRU and GRID model at 30, 60, 90, 150, 250, 500 and 1000 m (14 total) was run during the same period used for the 30-m baseline model calibration. The same SWAT version (Ver. 574) was used for both HRU and GRID models. HRU models were only run with standard routing, while GRID models were simulated twice, once with landscape routing and once with standard routing used in HRU-based delineation.
Predictive uncertainty derived from varying the input data resolution and model type was based on relative error (
RE):
where
Simbase is the simulated value from the calibrated, baseline model and
Simi,j is the simulated value from the coarsened dataset of resolution
i (30, 60, 90 m, …) and model
j (HRU, GRID-STD, GRID-LAND). Simulated values refer to flow at the watershed outlet. A baseline model was used rather than USGS-measured streamflow to focus the analysis on differences between the modeling approaches.
2.7. Calibration Effects
Given the significant computation time required by the GRID models, the calibration procedures involved calibrating a standard SWAT HRU-based model and transferring the parameters [
33]. As SWATgrid resolution is user defined, it is important to understand the effect of modifying input resolution on calibrated parameter values. Therefore, for the six basin-level parameters used for calibration, HRU models were calibrated individually at all tested resolutions, and resulting parameter values were compared over the final 1000 simulations using the same calibration approach described for the baseline model.
4. Conclusions
Current SWAT applications utilizing the standard HRU-based approach are limited by the absence of a spatially-explicit landscape position for HRUs. SWATgrid, a gridded approach that utilizes a modified routing algorithm to allow interaction between grid cells, was developed to overcome this limitation. However, SWATgrid remains largely untested with limited understanding of the effects of user-defined input data resolution on simulation results. Further, computational overhead introduced by SWATgrid restricts direct calibration, instead relying on parameter transfer from HRU-based models. To better capture and understand SWATgrid simulations, three objectives were considered: (1) identifying an appropriate resolution for SWATgrid simulations; (2) comparing SWATgrid to commonly-used HRU-based SWAT models; and (3) evaluating the effect of input resolution on model calibration.
Several SWAT models were tested over a range of resolutions with HRU-based, as well as GRID-based approaches, with and without a modified routing technique. Corresponding to each objective, it was found that: (1) flow and water quality simulated outputs tended to decrease with reduced resolution; (2) the SWATgrid approach underpredicted flows and some nutrients; underprediction was greater when using modified routing; and (3) calibrated parameter ranges showed some stability below 90 m. Another key finding was a simulation limit for the GRID approach as a function of both watershed area and grid cell size, which for the watershed tested was 90 m. Differences between HRU-based and GRID models was due to differences in subsurface flows and increased opportunity for loss of water to hydrologic landscape processes in the GRID-LAND approach. Further, the different routing techniques employed by the models likely necessitate separate calibrations, rather than parameter transfer from HRU-based to GRID models, although the transfer may serve as an initial estimate.
Overall, these findings underscore the need for future testing and possible improvements to the SWATgrid approach. However, it is important to recognize that based on the design of this study, the output of the SWATgrid-based models is not necessarily incorrect, but rather different from HRU-based simulations, which were directly calibrated at the 30-m level. In its current state, SWATgrid is most useful in simulating relatively small watersheds (<500 km2) where it is beneficial to capture spatial detail and interaction. As SWATgrid models at 30 and 60 m were not simulated, it is difficult to explicitly define an appropriate scale, although 90 m may be appropriate for watersheds similar to Cedar Creek. Future efforts should focus on a more in-depth analysis of the drivers and causes of simulation differences between SWAT modeling approaches at other locations, improving the calculation and use of the flow separation index, testing the model at higher resolution, using higher resolution soils data and working to improve the computational efficiency of the SWATgrid model. A hybrid approach that utilizes HRU-based SWAT to identify critical subbasins and subsequently modeling these subbasins with the GRID approach may offer an opportunity to synergistically capitalize on the strengths while minimizing the weaknesses of both approaches. Although the choice of model structure and resolution inherently remains a function of the intended application, results from this study can provide guidelines to identify an appropriate approach for future SWATgrid applications in similar regions.