Modeling Agricultural Watersheds with the Soil and Water Assessment Tool (SWAT): Calibration and Validation with a Novel Procedure for Spatially Explicit HRUs

Applications of the Soil and Water Assessment Tool (SWAT) model typically involve delineation of a watershed into subwatersheds/subbasins that are then further subdivided into hydrologic response units (HRUs) which are homogeneous areas of aggregated soil, landuse, and slope and are the smallest modeling units used within the model. In a given standard SWAT application, multiple potential HRUs (farm fields) in a subbasin are usually aggregated into a single HRU feature. In other words, the standard version of the model combines multiple potential HRUs (farm fields) with the same landuse/landcover, soil, and slope, but located at different places of a subbasin (spatially non-unique), and considers them as one HRU. In this study, ArcGIS pre-processing procedures were developed to spatially define a one-to-one match between farm fields and HRUs (spatially unique HRUs) within a subbasin prior to SWAT simulations to facilitate input processing, input/output mapping, and further analysis at the individual farm field level. Model input data such as landuse/landcover (LULC), soil, crop rotation, and other management data were processed through these HRUs. The SWAT model was then calibrated/validated for Raccoon River watershed in Iowa for 2002–2010 and Big Creek River watershed in Illinois for 2000–2003. SWAT was able to replicate annual, monthly, and daily streamflow, as well as sediment, nitrate and mineral phosphorous within recommended accuracy in most cases. The one-to-one match between farm fields and HRUs created and used in this study is a first step in performing LULC change, climate change impact, and other analyses in a more spatially explicit manner.


Introduction
Watershed heterogeneity and limitations associated with monitoring make it impractical to measure every aspect of the hydrological system (Pechlivanidis et al. 2011). Hydrological models are frequently used to overcome these limitations and extrapolate information from available measurements in both time and space to the watershed scale. Ecohydrological models are important for a wide range of applications such as water resources planning, watershed development and management, flood prediction and design, water quality evaluation, hydro-ecology, and climate change analyses. The Soil and Water Assessment Tool (SWAT) watershed-scale ecohydrological model (Arnold et al. 1998;Arnold and Fohrer 2005;Gassman et al. 2007) is currently one of the most widely used ecohydrological models and it has been extensively tested for a wide variety of watershed scales and environmental conditions worldwide (Gassman et al. , 2014Douglas-Mankin et al. 2010;Tuppad et al. 2011;Krysanova and White 2015;Bressiani et al. 2015).
Applications of SWAT typically involve delineation of a watershed into subwatersheds/subbasins that are then further subdivided into hydrologic response units (HRUs). HRUs are homogeneous areas of aggregated landuse, soil, and slope and are the smallest modeling units used in the model. The incorporation of HRUs in SWAT has provided flexibility for simulating a broad spectrum of conditions and supports adaptation of the model for watershed scales ranging from small field plots to entire river basins .
ArcSWAT, the standard ArcGIS input interface for SWAT, can spatially identify potential HRUs within a subbasin based on the above-mentioned HRU definition principles. However, ArcSWT-generated HRUs do not necessarily correspond spatially to individual farm fields in a given standard SWAT application. In other words, multiple potential HRUs (farm fields) with the same landuse/landcover, soil, and slope, but located at different places of a subbasin, are considered as one HRU. This characteristic of HRUs is referred, herein, as spatial non-uniqueness. Moreover, if lumping thresholds for HRU definition are used, spatial identification of individual HRUs will often be impractical. This is a key weakness in finding a one-to-one match between farm fields and HRUs, and hence assigning and presenting various inputs (e.g., rotation, tile drainage, manure and fertilizer application rates, etc.) and outputs at the HRU level. In this study, however, we introduce a data pre-processing procedure, discussed below, to create HRUs that are spatially unique so that there is a one-to-one match between HRUs and farm fields in a subbasin.
There have been other recent efforts in mapping SWAT outputs to a field level for identifying priority pollutantcontributing areas (e.g., Ghebremichael et al. 2010;Daggupati et al. 2011;Pai et al. 2012) so that appropriate management actions could be taken for specific land parcels. However, establishing spatially explicit (unique) HRUs prior to a SWAT simulation rather than after the model has been run could be extremely useful in certain applications. Arabi et al. (2008) subdivided the 7.3 km 2 Smith Fry watershed in Indiana into 97 subbasins covering an average of 3 ha each, and used dominant landuse and soil type combinations such that each subbasin consisted of a single HRU, in order to more accurately simulate the impacts of different best management practices (BMPs) on sediment loss. Bekele and Nicklow (2005) used a similar delineation scheme for the 133 km 2 Big Creek watershed in southern Illinois to assess the role of agricultural landscapes in reducing non-point source pollution while maximizing annual gross margin earned through agricultural production. However, the average subbasin size they used (175 ha) was much larger than those used by Arabi et al. (2008). The use of dominant HRUs (i.e., subbasins equivalent to HRUs) is not practical for many SWAT applications that are conducted at larger scales because key landscape, soil, and other details cannot be accounted for at the level of accuracy needed for the type of analyses being performed.
The ability to defining spatially unique HRUs in delineated subbasins prior to a SWAT simulation could be useful in simplifying assignment of inputs (rotation, tile drainage, manure application, etc.) and mapping of inputs/outputs at HRU level. It is important to note though that such ''spatially defined HRUs'' in the standard SWAT interfaces are still not simulated in a spatial manner within a simulation. A similar study on spatial one-to-one correspondence between HRUs and fields was done recently by Kalcic et al. (2015) for a small watershed (56 km 2 ) in west-central Indiana. They used the Common Land Unit (CLU) field boundary data developed by the U.S. Department of Agriculture along with assigned unique soil names to create spatially unique HRUs in a subbasin. After clean-up of the CLU boundaries using ArcGIS and updates of the soil database for the new unique soil names, ArcSWAT was used to configure the model. Kalcic et al.'s work was successful in demonstrating how ArcSWAT could be utilized to create a one-to-one correspondence between HRUs and CLU fields. However, one must purchase the CLU data, distributed prior to the 2008 Farm Bill, to utilize their approach, and their methodology is of course not replicable for watersheds for which shape files of field boundaries are not available. Our procedure extends this previous work by relying only on raster data, which is generally freely available.
The SWAT model also requires many input parameters related to landuse, soil, weather, topography, water quantity and quality, which may need to be calibrated and validated prior to using the model for specific analyses. Calibration and validation of a SWAT model for a watershed is essential for reducing uncertainties and increasing the confidence of the user for effective and efficient analysis (White and Chaubey 2005;Jha 2011). SWAT can be calibrated and validated at the daily, monthly or annual time scales depending on the purpose of the specific modeling exercise. The most commonly calibrated SWAT output is streamflow, especially at annual and monthly time steps, although an increasing number of SWAT studies are reporting testing of daily streamflow results (Arnold et al. 2012;Gassman et al. 2007Gassman et al. , 2014Douglas-Mankin et al. 2010;Tuppad et al. 2011;Bressiani et al. 2015). Streamflow is calibrated more often than water quality in part because it is essential for the other water quality components of the model (Gikas et al. 2005) and also because observed flow data are relatively abundant. On the other hand, sediment and nutrient parameters are not calibrated and validated as often, especially at the daily time scale (Arnold et al. 2012;Gassman et al. 2007Gassman et al. , 2014Douglas-Mankin et al. 2010;Tuppad et al. 2011;Bressiani et al. 2015). Calibration and validation of water quality parameters (sediment and nutrients) of SWAT at coarser time scales is mainly attributed to scarcity of observed water quality data at finer time scales (Wu and Chen 2009).
For water quantity and quality analysis in SWAT, there are three groups of parameters: Flow, Sediment, and Nutrients, which could be calibrated either separately (e.g., Muleta and Nicklow 2005;Gikas et al. 2005;Jha et al. 2007;Chahinian et al. 2011, etc.) or simultaneously (Kaur et al. 2004;Tolson and Shoemaker 2008;Wu and Chen 2009, etc.). Though challenging, the latter procedure seems to be preferable for improved calibration and validation results (Chahinian et al. 2011) due to the fact that certain parameters, such as the curve number at moisture condition II (CN2), affect all flow, sediment, and nutrient concentrations. Proliferation of computers with high processing capacity is likely to lead to greater use of simultaneous calibration. In this study, we primarily perform simultaneous calibration/validation and also perform separate calibration/validation procedures for comparative purposes.
The objectives of this study are, therefore, to (1) present an ArcGIS-based NASS Cropland Data Layer (CDL) (NASS 2012) landuse/landcover data and soil data preprocessing method to create spatially unique HRUs in a subbasin, (2) show the significance of the pre-processing procedure in assigning and visualizing inputs and/or outputs, and providing a suitable platform for future landuse change impact analysis at the farm field level, and (3) utilize the results of the pre-processing procedure during ArcSWAT configuration and calibration/validation of the SWAT model for two Midwestern U.S. watersheds at daily, monthly, and annual time scales.

Study Area
Based on availability of data and previous research, two watersheds in the midwestern U.S. have been selected for this study (Fig. 1): (1) the Big Creek River watershed (BRW) in Illinois, which is a combination of two HUC12 or 12-digit watersheds (USGS (U.S. Geological Survey) 2012), and (2) the Raccoon River watershed (RRW) in Iowa, which is a combination of two HUC8 or ''8-digit'' watersheds (USGS (U.S. Geological Survey) 2012). Basic characteristics of each watershed are shown in Table 1.
The RRW is intensively agricultural (primarily corn and soybean), with agricultural land cover accounting for 76.2 % of the watershed. The BRW agricultural area, which is most frequently cropped in soybean, covers only 25.2 % of the watershed area. Both watersheds receive ample annual rainfall and therefore irrigation is not a common practice.

SWAT Ecohydrological Model
SWAT is a watershed-scale model developed to estimate the impacts of various landuse and management practices on water quantity and water quality over a continuous long period of time ). The model is proven to be efficient in using readily available data and in studying long-term impacts (Neitsch et al. 2011;Arnold et al. 2012). The model is usually executed at the daily temporal scale, although sub-daily time step applications can also be done. ArcSWAT (Olivera et al. 2006;SWAT 2015) and MWSWAT (MapWindow SWAT; George and Leon 2008) are the two GIS-based graphical input interfaces which could be used to configure a SWAT model in a GIS environment, which are commercial versus open source software, respectively.

SUFI-2 (Sequential Uncertainty Fitting Version 2)
The SUFI-2 calibration tool is part of the stand-alone, public domain SWAT Calibration and Uncertainty Program software (SWAT-CUP), which has its own interface and is used for sensitivity analysis, uncertainty analysis, and calibration and validation of SWAT model parameters . SUFI-2 is widely used (e.g., Abbaspour et al. 2007;Schuol et al. 2008;Faramarzi et al. 2009;Gong et al. 2011) mainly due to the relatively fewer required number of runs to reach an acceptable calibration results. SUFI-2 requires 2-30 times fewer runs than the other programs of the SWAT-CUP . In this study, therefore, SUFI-2 has been used for sensitivity analysis, and calibration and validation of the SWAT models.
Uncertainties from driving variables (e.g., rainfall), model, parameters, and observed data (e.g., flow) are all accounted for in parameter uncertainty during a SUFI-2 calibration/uncertainty analysis. The measure of the degree to which these uncertainties are accounted for is defined as the P-factor (the percentage of measure of data bracketed by the 95 % prediction uncertainty, 95PPU) and the Rfactor (the ratio of average thickness of the 95PPU band to the standard deviation of the measured data). These two measures work together to bracket most of the measured data with the smallest possible uncertainty band. The regression correlation coefficient (R 2 ) and the Nash-Sutcliffe coefficient of efficiency (NS e ) between observed data and the final best simulation output are further quantifications of goodness of fit during SUFI-2 calibration/uncertainty procedure. While R 2 (ranges from 0 to 1) measures how well the observed data are correlated to simulated values, NS e (ranges from -? to 1) measures how well the observed and simulated data match. For both statistics, values of close to one indicate a well-calibrated and validated model. See Krause et al. (2005) for further description of the R 2 and NS e , as well as statistical implications of using the two statistics for hydrologic and water quality model evaluations.
Land Use and Land Cover (LULC) and Soil Pre-processing Procedure As discussed in ''Introduction'' section, HRUs are not spatially unique in a given standard SWAT application. Even though ArcSWAT can identify potential HRUs spatially, these HRUs are not spatially unique within a subbasin for a one-to-one match between agricultural fields  and HRUs. Hence, an ArcGIS-based data processing methodology has been developed to overcome this problem before using ArcSWAT to create HRUs. Here, ArcGIS tools have been utilized to develop a landuse/landcover (LULC) and soil pre-processing procedure to define spatially unique HRUs with in a subbasin. These HRUs are still simulated in a nonspatial manner in SWAT but the soil, landuse, and topographic data used to define the HRUs are geo-referenced, which allows more realistic spatial identification of the model inputs and outputs, and easier linkages with other spatially explicit models such as economic or ecological ones. This procedure requires freely available landuse and soil raster data, river networks, road networks, and pre-defined subbasins for the watershed to be processed. The general procedure is as follows: 1. Convert LULC to polygons (based on LULC type) 2. Exclude roads and ''very small'' urban lands from LULC polygons and re-assign them with the neighboring LULC type. In this procedure, roads are considered to be dividing lines between LULC types. Hence, this step helps to eliminate ''unnecessary and small'' polygons due to roads, very small towns, and isolated buildings in the middle of farmland or forest.

Split new LULC polygons with road and river
networks. In addition to roads, river networks are also considered to be dividing lines between LULC types. Like roads, LULC types on opposite sides of river networks could potentially be different. However, in this particular procedure, forests and urban areas are considered intact throughout the analysis; hence, forest and urban area polygons are excluded from the split. 4. Eliminate ''small'' polygons. Here, small may be subjective but generally the size of the polygons eliminated depends on the size and LULC pattern of the watershed. Moreover, elimination has to be made so that the number of polygons created is less than the number of HRUs that could be managed by the SWAT simulation processing facilities. The result of this step can be exported to GoogleEarth and compared with the actual boundaries of fields for an accuracy check. 5. Intersect the final polygons, after elimination, with predefined subbasin polygons (from ArcSWAT). One must ensure that the subbasins used here are exactly the same as those created during the ArcSWAT watershed delineation. 6. Re-label LULC polygon codes by assigning different codes for the same LULC types within a subbasin. For example, if there are three CORN polygons within subbasin 1, then these polygons are re-assigned with names like COR1, COR2, and COR3. These same new crop types should also be added into the SWAT crop database but should all have the same crop parameters, i.e., crop parameters for CORN. The same is true for other LULC types. 7. Finally, convert the re-labeled LULC polygons into raster-based data which will be used for the actual analyses. 8. Use the updated LULC polygons to assign a dominant soil for each polygon for the watershed analyses.
The main objective of the above procedure is to define spatially unique HRUs within a subbasin using ArGIS tools before an ArcSWAT model is developed. This will assist in assigning management input data in a spatially explicit manner for each HRU (Figs. 2 and 3), presenting output results at the HRU level, and assisting in explicitly analyzing future LULC change possibilities at the HRU level in SWAT modeling (but keeping in mind the limitation that the HRUs are not simulated spatially internally in SWAT).
The raw LULC layers, 2000 LULC for BRW and 2010 LULC for RRW, were utilized to define the HRU boundaries from which the new gridded LULC data were constructed ( Fig. 4). Once the HRU boundaries (polygons) were set using a single year of data, the polygons were utilized in determining the dominant LULC type (i.e., cropland or non-cropland) for each polygon (HRU). Crop rotations (Fig. 2) for the cropland HRUs were then determined by overlaying multiple years of crop land use data on each polygon, using crop data obtained from the USDA Cropland Data Layer (CDL; USDA-NASS, 2012) for the Raccoon (2002) and Big Creek (2000-2003 watersheds. Using multiple crop years to construct the land use is crucial in watersheds where farmers rotate annual crops, because the various crops can have very different water quality impacts, and using only 1 year of data may mask long-term trends in land cover and affect the accuracy of the calibration process.

Load Estimator (LOADEST)
The Load Estimator (LOADEST) is a software package developed by USGS to generate water quality parameters through regression based on observed water quality data (grab samples) and their corresponding flow rate and time of observation (Runkel et al. 2004). LOADEST has eleven pre-defined regression models to choose from. Moreover, the user-defined option of LOADEST can be utilized to incorporate additional parameters to improve the prediction capability of the regression equations. The regression model can be developed for daily, monthly, or annual data generation. In this study, LOADEST has been utilized to interpolate water quality data to be used during calibration and validation procedures.

Data
Recent 30-m gridded Digital Elevation Model (DEM) data were downloaded from the United States Geological Survey (USGS) National Map Viewer and Download platform (USGS (U.S. Geological Survey) 2012). Likewise, 30-m gridded landuse (CDL; USDA-NASS, 2012) and countybased soil data from the Soil Survey Geographic Database (SSURGO) (USDA-NRCS, 2012) were downloaded from the USDA respective geospatial data gateways. Climate data (daily precipitation and maximum and minimum temperature) from the National Oceanic and Atmospheric Administration's (NOAA) National Climatic Data Center (NCDC) (NOAA-NCDC, 2012) were obtained for ten weather stations for the RRW and one weather station for the BRW (Fig. 1).
Data were collected from one RRW stream gauge and two BRW stream gauges (Fig. 1) to perform the calibration/validation processes. While about 95 % of the RRW drains to the gauging station, about 64 and 17 % of the BRW area drains to its downstream and upstream gauging stations, respectively. Multiple sources were utilized in collecting daily flow and water quality (Total Suspended Sediment, TSS, Nitrate, NO 3 , Mineral Phosphorous, MINP) data, most of which were obtained from the Environmental Protection Agency STOrage and RETrieval (EPA-STORET) online database. Sources such as USGS, the Des Moines Water Works (DMWW), the Des Moines River Water Quality Network (DMRWQN), the Illinois Water Survey website, and individual e-mail requests were also utilized for flow and water quality data collection. The simulation periods for each watershed were determined based on the availability of water quality data as shown below (Table 2).
For the watersheds considered in this study, urban, forest, and pasture landuse/landcover types were assumed to remain the same for the entire simulation periods during the calibration and validation processes. Muleta and Nicklow (2005) assumed that corn, soybean, pasture, and hay are grown without tillage in the BRW based on interviews with personnel from the Southern Illinois District office of the NRCS (Natural Resources Conservation Services). However, there is no evidence that manure applications occur or that tile drains are used in the BRW. In the RRW, on the other hand, information about tillage practices, manure applications, and tile drainage were available based on previous SWAT model studies reported by Jha et al. (2007Jha et al. ( , 2010 and Schilling et al. (2009). Jha et al. (2010) and Schilling et al. (2010) reported that approximately 77.5 % of the row cropland in Northern Raccoon watershed and 42 % in the South Raccoon watershed were tile drained (Fig. 3), and that tile drains were installed (1200 mm below the ground) on about 48 % of the cropland overall, based on data obtained from the Iowa Department of Natural Resources (IDNR 2013).
Estimated manure application distribution on cropland, at a rate of 160 lbs N/acre (*179 kg N/ha) from animal feeding operations in Iowa in 2006 (Fig. 3), was also  1997, 1998, 1999, 2000, 2001, 2005, and 2010. Data from the years within and closest to the simulation periods (2000, 2001, 2005, and 2010) were used for the RRW analysis. The average values of rates of nitrogen and phosphorous fertilizer applications, and percentages of three tillage practices [conventional (Cv), conservation (Cs), and no-till] for corn and soybeans were estimated for RRW using data from those years mentioned above ( Table 3).
The BRW is located south of the main southern Illinois crop production region, and thus, basing the management practice on state average values could not be adequately justified. Therefore, the estimated nutrient application rates and tillage practices used for the BRW SWAT simulations were based on expert opinion provided by local agency personnel (Table 3).

Model Baseline Setup
The SWAT2009/rev. 481 and ArcSWAT 2009.93.7b were utilized as a tool and interface, respectively, in this study. Pre-processed landuse and soil data along with the DEM data were used to construct subbasins and HRUs for each 0 20 40 60 80 10 Km   Accordingly, 2481 HRUs in 154 subbasins and 1644 HRUs in 79 subbasins were constructed for the RRW and BRW, respectively. In addition to the simulation periods indicated in Table 2, 2-year warm-up periods were introduced for each watershed. Crop rotations were constructed using multiple years of LULC data for both watersheds (Fig. 2)   Since all of the above-mentioned rotation percentage values were results of estimations based on the pre-processing procedure, the proportion of agricultural/pasture/forest/developed/open water areas could be slightly different from what was reported in Table 1 or the raw LULC type proportions of the base year used for pre-processing (see ''LULC and Soil Pre-processing Procedure'' section). Management operations such as fertilizer and/or manure applications, tillage practices, and tile drainage were implemented for each pre-defined HRU according to the availability of data for the watershed. For the RRW, tile drainage and manure applications were distributed throughout HRUs based on the acreage of the management operation in each HRU. If an operation (tile drainage or manure application) covered more than 50 % of HRU's area, then the operation was applied for the HRU. Otherwise, the operation was excluded. Figure 3 shows original tile drainage and manure application distribution data along with processed distribution results used in this study for the RRW. Fertilizer application rates estimated based on data from the USDA-ERS were applied for each HRU in accordance with crop rotations (Table 3).

SWAT Calibration and Validation
Both automatic calibration (using SUFI-2) and manual calibration procedures were used in this study. Initially, the curve number at moisture condition II (CN2), saturated hydraulic conductivity of the soil (SOL_K), and available water capacity of the soil (SOL_AWC) values for various LULC types were adjusted manually to calibrate, at least partially, the annual water/sediment/nutrient balances and the baseflow contributions to the average annual total water yields of both watersheds. In addition to the above parameters, the threshold depth of water in the shallow aquifer required for return flow to occur (GWQMN) for BRW, and fraction of porosity from which anions are excluded (ANION_EXCL) for the RRW were found to be sensitive to annual water yield and its baseflow fraction, and annual nitrate output, respectively, and hence were adjusted manually as well.
Previous studies on RRW show that baseflow is an important component of the water balance for the watershed (Jha et al. 2007(Jha et al. , 2010Schilling et al. 2009). Hence, baseflow needs to be calibrated before stream flow and other components to enhance modeling accuracy. The automated method of estimating baseflow, developed by Arnold and Allen (1999), is used to estimate the fraction of baseflow contribution for the total annual flow. Accordingly, for the RRW, baseflow is estimated to contribute about 61.5 % on an average annual basis for 2002-2010. This value is consistent with Jha et al. (2007) estimates of 58 % for 1981-2003 using the same method of estimation. Similar estimates for the BRW resulted in baseflow contributions of about 28 % on an average annual basis for 2000-2003. Since previous studies on the BRW did not consider baseflow separation during calibration/validation of the watershed, there was no information to compare the result to. The BRW baseflow component is expected to be lower than the RRW due to the lack of tile drainage in the watershed. However, the BRW LULC is mainly forest and pasture which would increase the baseflow contribution. Therefore, it is possible that the baseflow contribution for BRW may be underestimated due to the short duration (4 years) of flow record considered in this study.
The calibration/validation periods for RRW and BRW were 2002and 2000-2001/2002-2003  respectively, except that the validation period for the BRW nutrient simulation was only the year 2000 due to data limitations. The upstream gauging station of the BRW was used as further verification of the model performance for the watershed. Both simultaneous and separate calibration/validation procedures were performed by changing relevant SWAT parameters (Appendix Table 6). These calibration and validation processes were performed at daily and monthly time steps for each watershed, and also at an annual time step for the RRW. The model predictions, during both calibration and validation periods, were evaluated using graphical and statistical comparisons. The R 2 and NS e , statistical measures of SUFI-2 discussed in ''SUFI-2 (Sequential Uncertainty Fitting version 2)'' section, were used for evaluation. Evaluations were performed at daily, monthly, and annual time scales. Though there are no firm criteria for minimum required values of R 2 and NS e , an NS e value of at least 0.5 is recommended by Moriasi et al. (2007) for the SWAT model to be satisfactory for hydrologic and pollutant loss evaluations on a monthly time scale. Therefore, Moriasi's recommendation was considered as the minimum goal for both R 2 and NS e at monthly time scales, and the minimum values of R 2 and NS e for daily and annual time scales were expected to be less and greater than that of Moriasi's recommendation, respectively.

LULC and Soil Pre-processing Procedure
The percentage of each LULC type changed as a result of the pre-processing procedure (Fig. 4). Decreases from 11.5 to 10.0, 7.3 to 5.9, and 0.7 to 0.5 % for pasture, developed, and open water areas, respectively, and increases from 75.6 to 79.2 and 4.1 to 4.4 % of agricultural and mixed forest areas, respectively, were observed for the RRW. For the BRW, decreases from 43.2 to 37.7, 29.9 to 27.5, and 0.3 to 0.2 % for agricultural, pasture, and open water areas, respectively, and increases from 23.7 to 32.6 and 0.9 to 1.5 % for mixed forests and developed areas, respectively, were observed due to the pre-processing procedure.
Some of the LULC changes above may look objectionable at first glance; however, it should be noted these changes are similar to those that are inevitable if threshold levels are used in SWAT to eliminate minor LULCs during HRU definition. Detailed descriptions of the processes during the HRU definition in SWAT modeling can be found in the ArcSWAT user's guide by Winchell et al. (2007). The data pre-processing procedure introduced in this study performs a similar function in defining the HRUs except that the pre-processed data helps define spatially unique HRUs in a subbasin during ArcSWAT configuration. It is also possible that the quality of the LULC data may have an influence on the results of the pre-processing procedure, though exploring this issue goes beyond the scope of this work, and therefore, we cannot offer firm conclusions on this issue. According to Johnson and Mueller (2010), ''the overall quality and accuracy of the CDLs has steadily increased over time due to better ground truth and greater access to imagery.'' Therefore, using a more recent (2010) LULC data for the RRW may be the reason for lower changes in LULC between raw and processed data, compared to the 2000 LULC data used for BRW that produced higher changes in LULC. However, this could also be due to the fact that the RRW landscape is much more homogenous than that of the BRW.
One of the applications of results of this procedure is that multiple years of LULC data could be utilized to determine HRU-level crop rotations explicitly (Fig. 2). This will in turn help determine annual fertilizer application rates more accurately as shown in Table 3. In addition to this, HRU boundaries from the pre-processing procedure were used to determine the soil type for each HRU. The dominant soil type for each HRU polygon was assigned, and the resulting soil data were used for the SWAT modeling. Pre-defined HRUs were further utilized to determine which area was likely tiled and received manure. Accordingly, processing raw tile drainage and manure application data from the IDNR through the pre-defined HRUs resulted in a tiled area of 57 % (as compared to the original 48 %) of the entire RRW and manure application in 20 % of the watershed, respectively.
Comparisons between observed and simulated flow, sediment, nitrate, and mineral phosphorous (Table 4; Fig. 5 and Appendix Figs. 9 and 10) show that the calibration processes produced representative SWAT models for both watersheds, with more accurate results for the RRW. In the RRW, annual and monthly comparisons (Table 4) were comparable with previous results by Jha et al. (2007;, except the MINP values. However, at daily time scale, much better results of R 2 /NS e were obtained in this study compared to the daily calibration results of Jha et al. (2010). The R 2 /NS e values for daily calibration were 0.81/0.80 for flow, 0.54/0.53 for sediment concentration, 0.69/0.66 for NO 3 load, and 0.47/0.45 for MINP load.
Having a longer and more recent time series of data to calibrate and validate to, it was expected that the RRW would do much better during calibration and validation than the BRW. These expectations were proven accurate by the BRW calibration results (Table 4; Fig. 6 and Appendix Fig. 11). Generally, in the BRW, less accurate calibration results were obtained. Due to the very short series of data available for the BRW, annual calibration procedures using SUFI2 were not performed. Daily flow and sediment concentration calibration results were consistent with previous studies by Muleta and Nicklow (2005) and Bekele and Nicklow (2007) for the watershed. Moreover, the monthly calibrations (R 2 /NS e value of 0.66/0.54 for flow, 0.74/0.70 for sediment concentration, 0.56/0.30 for NO 3 load, and 0.85/0.68 for MINP load), which were performed only in this study, resulted in significantly improved calibration values, especially during water quality calibration.
Validations of the RRW were performed for years 2008-2010 for all the calibrated outputs at the same gauging station that the calibration was performed for. Even though the R 2 /NS e values during validation were generally less than the calibration period, most of the values were satisfactory (  Fig. 8 and Appendix Fig. 14). The BRW model was also validated for flow at the upstream gauging station for years 2000-2003. The R 2 /NS e values for the upstream gauging station were 0.45/0.42 at daily and 0.58/0.46 at monthly time steps.
Finally, comparisons were made between R 2 /NS e results of simultaneous and separate (stepwise) calibration procedures. As can be deduced from the results presented in Table 5, neither procedure provided a clear advantage on the basis of the R 2 /NS e values. Calibrated SWAT parameter values are presented in Appendix Table 7 in the appendices for convenience and future reference.

Conclusions
The LULC pre-processing procedure introduced in this study was the first step in defining spatially unique HRUs for the SWAT model before the simulations. Its applications in assigning and presenting certain management data inputs and presenting model outputs at the HRU level in a spatially explicit manner including accounting for multiyear crop rotations were demonstrated in this study. The procedure could be used to predict LULC changes in a watershed at the HRU level in studies such as multi-objective management of ecosystem services (Bekele and Nicklow 2005), coupling agent-based LULC change modeling with SWAT model (Ng et al. 2011), and strategic targeting and prioritization of areas that need certain kind Environmental Management (2016) 57:894-911 905 of management practice implementation (Daggupati et al. 2011). Moreover, the approach could be utilized anywhere in the world as long as landuse, soil, DEM, and road network data are available. In the absence of road network data, the procedure could still be used but the accuracy of the results will likely be lower.
The R 2 and NS e results (Table 4) for the RRW showed that the SWAT model was able to replicate annual, monthly, and daily streamflow, as well as sediment, nitrate, and mineral phosphorous, within recommended accuracy as suggested by Moriasi et al. (2007) except for a few cases. Due to limited and older observed data, less accurate results were obtained for the BRW simulations. However, except for nitrate, the BRW monthly calibration results were also satisfactory as per Moriasi's performance ratings. Despite the comparability of both simultaneous and separate calibration procedures, the importance of simultaneous calibration with respect to its accuracy presented in this study should also be noted. It is also important to remember that there are bias issues with estimating pollutant loads with LOADEST, especially for daily load estimates (Stenback et al. 2011). Thus, future research is needed to improve upon current available load estimation methods.
Finally, the HRUs simulated in this study still do not overcome the problem that HRUs are simulated nonspatially within the actual SWAT simulations. An alternative approach is currently being developed for the SWAT model that allows spatial representation of different landscapes within a subbasin as well as simulated routing of flow and pollutants between the different landscapes (e.g., see Bosch et al. 2010;Arnold et al. 2010;Bonumá et al. 2014). HRUs are delineated within the landscape units in this approach; these HRUs are still not spatially identified, although a single HRU could be used to represent a given landscape unit depending on the application [Arnold (2014) Personal communication. USDA-ARS Grassland, Soil and Water Research Lab., Temple, TX]. Interfacing the HRU approach described in this study with the landscape structure being developed for SWAT will provide an even more enhanced method for spatially representing landuse, management practices, and other landscapespecific characteristics within the model.