Evaluation of AnnAGNPS Model for Runo ﬀ Simulation on Watersheds from Glaciated Landscape of USA Midwest and Northeast

: Runo ﬀ modeling of glaciated watersheds is required to predict runo ﬀ for water supply, aquatic ecosystem management and ﬂood prediction, and to deal with questions concerning the impact of climate and land use change on the hydrological system and watershed export of contaminants of glaciated watersheds. A widely used pollutant loading model, Annualized Agricultural Non-Point Source Pollution (AnnAGNPS) was applied to simulate runo ﬀ from three watersheds in glaciated geomorphic settings. The objective of this study was to evaluate the suitability of the AnnAGNPS model in glaciated landscapes for the prediction of runo ﬀ volume. The study area included Sugar Creek watershed, Indiana; Fall Creek watershed, New York; and Pawcatuck River watershed, Rhode Island, USA. The AnnAGNPS model was developed, calibrated and validated for runo ﬀ estimation for these watersheds. The daily and monthly calibration and validation statistics ( NSE > 0.50 and RSR < 0.70, and PBIAS ± 25%) of the developed model were satisfactory for runo ﬀ simulation for all the studied watersheds. Once AnnAGNPS successfully simulated runo ﬀ , a parameter sensitivity analysis was carried out for runo ﬀ simulation in all three watersheds. The output from our hydrological models applied to glaciated areas will provide the capacity to couple edge-of-ﬁeld hydrologic modeling with the examination of riparian or riverine functions and behaviors.


Introduction
Excess nutrient (primarily nitrogen and phosphorus) losses from agricultural watersheds in glaciated settings of the Midwest and Northeast of USA are one of the greatest water quality problems tied to modern agriculture [1][2][3][4][5][6]. These water quality problems include eutrophication, harmful algae blooms, and fish kills in the Gulf of Mexico, the Chesapeake Bay, the Hudson River Estuary, and other coastal areas [7][8][9][10]. A substantial body of research on riparian zone hydrology and biogeochemistry has shown that riparian zones can serve as efficient best management practices (BMPs) for nutrient removal [11,12].
The functional efficiency of nutrient removal in a riparian zone can vary widely depending on the characteristics within the riparian zone (e.g., vegetation, soil texture, depth to water table) and on its location in the landscape, timing, characteristics and extent of contaminant and hydrologic loading and its hydrogeomorphic setting [11][12][13][14][15][16][17][18][19]. Thus, it is significant to explore the relationships area included Sugar Creek watershed in Indiana, Fall Creek watershed in New York and Pawcatuck River watershed in Rhode Island. All the watersheds are located in a glacial geomorphic setting.
Specifically, this paper aims to: (i) evaluate the performance of the AnnAGNPS model in simulating the runoff volume at three separate watersheds with glacial setting; (ii) improve the model's runoff prediction capacity through calibration; (iii) validate the model's runoff prediction with the improved calibrated parameters; (iv) conduct a parameter sensitivity analysis for runoff simulation; (v) conduct an analysis of the spatial distribution of runoff depth for three watersheds; (vi) provide a discussion of the model's performance in order to estimate event peak discharge.

Study Area
Three watersheds from three different states were modeled for this study. Sugar Creek watershed (39 • 43 21" N, 85 • 53 23" W), a part of the White River watershed in central Indiana, is about 69 km 2 ( Figure 1). The elevation of the watershed ranges from 241 m to 280 m, and the topography is nearly flat. The watershed consists largely of tile-drained agricultural lands (88% of the total watershed area, representative of agro-ecosystems of the glacial till plains from US Midwest [46]. This watershed is dominated by poorly drained soils where artificial drainage is usually used to lower the water table [47]. For the past 20 years, agricultural practices have been dominated by a corn/soybean rotation with either conventional or conservation tillage systems [48]. The temperature in the watershed is moderate, ranging from a 30-year (1982-2011)   Fall Creek watershed has an area of about 328 km 2 is located within the Finger Lakes region of New York State (42 • 28 N, 76 • 27 W) ( Figure 1). The most extensive source of parent material is glacial till, with additional parent materials that consist of glacio-lacustrine sediments and glacio-fluvial (outwash) deposits. The watershed is a mixed land use landscape located at the southern terminus of the Wisconsin glaciation [49]. The watershed is 4.8% urban/developed land use (residential, commercial and service, industrial, etc.), 45.3% forest (evergreen forestland, mixed forestland), and 49.4% agriculture (cropland and pasture, other agricultural land, shrub and brush rangeland) [50]. Soils in the watershed are dominated by Gravelly silt loam and Channery silt loam. These are typically very deep, well-drained soils. Elevations range from 270 m above mean sea level to 600 m [51]. The temperature in the watershed ranges from a 30-year (1982-2011) (Figure 1). The area of this watershed is about 258 km 2 . It consists mainly of forests (above 65% of the total watershed area) and agricultural fields (about 32% of the entire watershed area). The soil parent materials in the watershed are comprised mostly of glacial till, glacial outwash, and organic and alluvial deposit [52]. Agricultural lands (mostly turf farms) are predominately located on loess soils over glacial outwash. Forested settings are usually on till. The elevation of the watershed ranges from 16 m (shoreline) to 144 m (to gently rolling hills inland). It has a humid continental climate, with warm summers and cold winters.

Description of The AnnAGNPS Model
The Annualized Agricultural Non-Point Source Pollution Model (AnnAGNPS) [24] refers to a watershed scale, batch process, continuous and distributed simulation, daily time step, surface runoff, and pollutant loading computer model. The model has been designed to quantify and identify the source of pollutant loadings anywhere in the watershed for optimization and risk analysis. Hydrology, sediment, nutrient, and pesticide transportation are essential modeling components. This continuous version of the model is an improvement to the previously developed single-event Agricultural NonPoint Source model (AGNPS) watershed model [53]. The model uses and combines many modules of other commonly used models, such as Revised Universal Soil Loss Equation (RUSLE) [54], Chemicals, Runoff, and Erosion from Agricultural Management Systems (CREAMS) [55], Erosion Productivity Impact Calculator (EPIC) [56], and Groundwater Loading Effects on Agricultural Management Systems (GLEAMS) [57]. In this article, AnnAGNPS version 5.45 (United States Department of Agriculture (USDA)-Agricultural Research Service (ARS), National Sedimentation Laboratory, Oxford, MS, USA) (Official Release-21 December 2016) was used for all simulations. A full description of this model and its associated components are available in [58].

Hydrological Modeling Component in AnnAGNPS
The main components within AnnAGNPS are the combination of the Soil Conservation Service (SCS) curve number (CN) technique [59] used to generate daily runoff and RUSLE 1.05 tool (USDA-ARS, Washington, DC, USA) [54] to produce daily sheet and rill erosion from fields [60]. AnnAGNPS divides the watershed into drainage areas called 'cells' that can have any shape, and each cell is assumed to have homogenous management and soil [28]. These cells portray the spatial variability of land use, soil, and topography within the watershed. These simulated cells are then integrated by simulated streams and rivers, which route the runoff and pollutants from every single homogeneous area downstream.

AnnAGNPS Data Input
For the execution of an AnnAGNPS model, the major input data are climate, land characteristics (e.g., topography, soils), field operations, chemical characteristics, and feedlot operations. Topography information about the three studied watersheds was acquired from the United States Geological Survey (USGS)-The National Map Viewer (TNM Viewer version 2.0, USGS, Washington, DC, USA) 7.5-min digital elevation models (DEMs)-with a 10-m horizontal, 7-m vertical resolution. It was used to obtain the necessary input data for running the TOPAGNPS (Topographic AGNPS) program version 5.45.a.011 (United States Department of Agriculture (USDA)-Agricultural Research Service (ARS), National Sedimentation Laboratory, Oxford, MS, USA), a Geographic Information System (GIS)-based landscape analysis component of AnnAGNPS that is used to generate the input parameters of the model. TOPAGNPS requires a user-selected watershed outlet location to produce the prerequisite model input files from the DEM dataset. The DEM was used to identify and measure the topographic features, to define surface drainage channels, to subdivide watersheds into cells along drainage divides and also to calculate representative cell parameters (cell area, slope, and length). The size of the cells depends on the values of the Critical Source Area (CSA) and Minimum Source Channel Length (MSCL) [34]. The CSA is defined as the minimum upstream drainage area required for a channel to form, while the MSCL is the minimum acceptable length of concentrated flow in a cell before a stream channel can be defined [61]. The CSA and MSCL values are critical to determining the extent of the stream network and resulting AnnAGNPS cells. Various combinations of CSA and MSCL values were applied until an accurate representation of the stream network and of the land use of the studied watersheds was acquired. For the three sites, CSA ranged from 5-170 ha and MSCL from 30-130 m ( Figure 1). The number of cells per watershed ranged from 185 to 1800. The soil data are directly populated from United States Department of Agriculture (USDA)-Natural Resources Conservation Service (NRCS) Soil Survey Center's National Soil Information System (NASIS) data. NASIS data are associated with The Soil Survey Geographic (SSURGO) soil map. This soil map was overlaid onto the delineated watershed using the AGNPS GIS tool, and the dominant soil type for each subwatershed cell was determined. Then land use map obtained from National Land Cover Database (NLCD 2011)-United States Geological Survey (USGS) was also overlaid onto the delineated watershed using the AGNPS GIS tool. The six daily climate parameters needed for AnnAGNPS are (1) minimum air temperature; (2) maximum air temperature; (3) precipitation; (4) dew point; (5) solar radiation; and (6) wind speed. The data for three daily climate parameters-minimum air temperature, maximum air temperature, and precipitation-were acquired from the PRISM website at 4km spatial resolution (Parameter-elevation Regressions on Independent Slopes Model (PRISM) Climate Group, Oregon State University, created 4 February 2004). The remaining three daily climate parameters-dew point, solar radiation, and wind speed-were acquired from Texas A&M University's global weather data site [60]. The traditional manual baseflow filtering approach was applied to the streamflow record to obtain runoff by removing baseflow from streamflow before comparison with AnnAGNPS output, as baseflow is not considered in the model [29].

Model Assessment
The performance of model was evaluated by comparing observed and AnnAGNPS modeled data at the watershed outlet. The assessment of the model was accomplished for runoff on both daily and monthly time scales. Assessment of model performance for runoff included both qualitative and quantitative methods. Qualitative methods included comparing graphs of observed and modeled data. We followed the recommendation of [62] and used three quantitative statistics: Nash-Sutcliffe efficiency (NSE), percent bias (PBIAS), and ratio of the root mean square error to the standard deviation of measured data (RSR), along with the graphical techniques, to model performance evaluation. Generally, model simulation can be judged as satisfactory when NSE > 0.50 and RSR < 0.70, and also when PBIAS ± 25% for streamflow [62]. We also used the coefficient of determination (R 2 ) for quantitative evaluations; R 2 represents the variation in measured data explained by the model [62]. Values can range from 0 to 1, with 1 indicating that all variations in the measured data are explained by the model. Values greater than 0.5 are normally considered acceptable [62]. According to Nash and Sutcliffe [63], NSE is a normalized statistic that defines the relative magnitude of the residual variance when compared to the variance in the measured data. The statistic denotes how well the observed data fit the modeled data in the 1:1 line. The NSE value ranges from −∞ to 1 with 1 representing a perfect fit. Values between 0 and 1 are considered an acceptable performance level for the model [62].
NSE is computed as shown in Equation (1): where Y obs i is the ith observation for the constituent being evaluated, Y sim i is the ith simulated value for the constituent being evaluated, Y mean is the mean of observed data for the constituent being evaluated, and n is the total number of observations.
PBIAS is computed as shown in Equation (2): RSR is computed as shown in Equation (3):

Model Calibration, Validation and Sensitivity Analysis for Runoff Simulation
The SCS-CN, the most important parameter in the AnnAGNPS model for simulating runoff, is utilized in many studies to calibrate runoff [34][35][36][37]. For that reason, the SCS curve number was also used to calibrate runoff in this study. For Sugar Creek watershed, the AnnAGNPS model was In this study, we also evaluated the sensitivity value of the most sensitive parameter (SCS-CN) used for runoff estimation. The sensitivity analysis for CN was performed by the modification or adjustment of the curve number within the recommended range (30-100). The lower numbers indicate low runoff potential, whereas higher numbers signify increasing runoff potential. We utilized the integration of a local method into a global sensitivity method (the random one-factor-at-a-time) design proposed by [64]. This method consists of repetitions of a local method whereby the derivatives are calculated for each parameter by adding a small change to the parameter. The change in model outcome can then be measured by some lumped measure such as total mass export, sum of squares error between modeled and observed values or sum of absolute errors. The sensitivity analysis and the calibration of streamflow for the AnnAGNPS model were manually calibrated as in other studies [36].
Based on [65], the selection of initial SCS CNs for the different land use types was completed. Sugar Creek watershed consists of various land uses like cropland (only corn), cropland (corn-soybean rotation), fallow land, forested and urban area. Initially, the CN for a straight row crop with poor hydrological conditions was used for both corn and corn-soybean rotation during the growing season, while the CN for a fallow field with crop residue and good hydrological conditions was used after harvest during the non-growing season. The CN for woods with good hydrological conditions was used for forested areas. The CN for urban areas with 85% impervious cover was used for urban areas ( Table 1). The sensitivity analysis was done only for croplands. The CN for Row Crop (SR-Poor) was adjusted by running the model a number of times and by relatively changing the value of CN from its initial value by ± 9.8% to ± 2.8% during the calibration phase. Fall Creek watershed consists of various land uses like urban (developed), forested, some cropland (tall grass, squash, potato), and hay/pasture. At the start of the calibration, the CN for a straight row crop with poor hydrological conditions was used for potato, while the CN for a contoured row crop with poor hydrological conditions was used for tall grass and squash fields. The CN for woods with poor hydrological conditions was used for forested areas. The CN for urban land use with the newly graded condition was used for urban areas (Table 2). For this watershed, we focused our sensitivity analysis on the croplands (tall grass and squash fields). The CN for Row Crop (C-Poor) was adjusted by running the model for a number of times by relatively changing the value of CN from its initial value by ± 13.6% to ± 2.6% during the calibration phase.
The Pawcatuck watershed consists mostly of forested areas, some urban areas and agricultural fields (turf). At the beginning of the calibration the CN for woods with good hydrological conditions was used for forested areas. The CN for a straight row crop with good hydrological conditions was used for turf during the growing season, while the CN for crop residue cover with good hydrological conditions was used after harvest during the non-growing season. The CN for urban area with 85% impervious cover was used for urban areas (Table 3). We performed our sensitivity analysis for cropland and forested area for this watershed. The CN for Row Crop (SR-Good), Row Crop (C-Poor) and Woods (Good) were adjusted by running the model a number of times by relatively changing the value of CN from its initial value by ± 10% to ± 30% during the calibration phase.  After the initial run of the model without calibration, the model was calibrated to support a better estimation of runoff. The model performance improved for both daily and monthly runoff calculations after calibration. The results were evaluated using both graphical and statistical methods [64] until the best simulation results were obtained. For runoff validation, all model parameters after calibration were kept the same, and the simulated data were compared with the observed runoff data. Following calibration in the Pawcatuck watershed, the CN for the turf crop was substantially lower than for a straight row crop with good hydrological conditions, reflecting the higher infiltration rates of turf grass (Table 3).

Runoff Calibration and Validation
According to the classification tabulated in [31] for model correlations and efficiencies modified from [62], our calibrated model for the Sugar Creek watershed predicted the daily runoff volume of the watershed with good correlation and good agreement (R 2 = 0.57, NSE = 0.57 for daily and R 2 = 0.67, NSE = 0.63 for monthly calibration) between daily observed and daily modeled runoff volume (Table 4, Figure 2a). The calibrated model, when applied to the same watershed for the validation phase, predicted a daily runoff volume with good correlation and good agreement for both daily and monthly scales (R 2 = 0.58, NSE = 0.57 for daily and R 2 = 0.72, NSE = 0.68 for monthly) (Table 4, Figure 2b). Total runoff estimation by the model during the calibration phase differed from the observed runoff by only about 6.44%, whereas it differed by about 20.5% during validation. The calculated PBIAS value for calibration was less than 10, which indicated an excellent calibration performance rate. The model was biased to overestimate runoff volume during both calibration and validation phases. The observed runoff volumes from January 2000 to December 2013 at the watershed outlet were used for model calibration and validation at daily and monthly scales. The model over-predicted some runoff volumes during the drier months (December to February), whereas it under-predicted some during the wetter months (May to August) ( Figure 2).  In the case of Fall Creek watershed, the statistical evaluation of model performance for calibration and validation is presented in Table 5. The value of NSE for both daily and monthly time scales is greater than 0.5, so the model calibration performance can be rated as good. The positive PBIAS indicated the overall underestimation of runoff by the model compared to the observed runoff volume (calibration phase) and the negative PBIAS indicated the overall overestimation of runoff by the model compared to the observed runoff volume (validation phase). Total runoff estimation by the model during the calibration phase differed from the observed runoff by about 16.5%, whereas it differed by about 5.5% during validation. The calculated PBIAS value for validation was less than 10, which pointed to an excellent model performance. Figure 3a,b graphically illustrates observed and modeled daily runoff volume at the USGS 04234000 for calibration and validation phase, respectively, for Fall Creek watershed. The model over-predicted some runoff volumes during the months in which less precipitation occurred (December to March), whereas it under-predicted some during the months in which more precipitation occurred (April to August) (Figure 3). This tendency of the model could be improved by adjusting the evaporation rate associated with the interception of precipitation events.  In the case of Pawcatuck River watershed, graphical comparisons of daily observed and modeled runoff volumes at the USGS 01117500 were presented in Figure 4a,b for the calibration and validation phase, respectively. The statistical evaluation of model performance for calibration and validation is presented in Table 6. The range of NSE (0.51 to 0.54 for daily and 0.68 to 0.83 for monthly) showed a good agreement between daily observed and modeled runoff volume. The positive PBIAS indicates the overall underestimation of runoff by the model compared to the observed runoff volume for both calibration and validation. Total runoff estimation by the model during the calibration phase differed from the observed runoff by about 21.5%, whereas it differed by about 5.8% during validation. The calculated PBIAS value for validation was less than 10, which pointed to an excellent model performance. The results show a general tendency for AnnAGNPS to overestimate spring (March-May) and summer (June-August) runoff volumes compared to observed data for both calibration and validation periods (Figure 4).  The estimated RSR values varied from 0.65 to 0.70 for daily (fair) and 0.41 to 0.63 for monthly (very good to fair) for all three watersheds during calibration and validation periods.

Sensitivity Analysis
After the sensitivity analysis, it was quite clear how sensitive the CN was for runoff simulation in AnnAGNPS model. The sensitivity analysis demonstrated differences between sites in the response of modeled runoff to changes in CN. In Sugar Creek, the percent change in runoff volume (2000-2007 or 7-year average) from its initial condition was found to be up to 5.5 due to 9.8% changes in the CN ( Figure 5). In the case of Fall Creek watershed, during the entire phase of modification of CN, the percent change in runoff volume (2000-2007 or 7-year average) from its initial condition ranged from 11 to 149 due to 2.6% to 13.6% changes in the CN (Figure 6). In the case of Pawcatuck River watershed, during the entire phase of modification of CN for Row Crop (SR-Good), the percent change in runoff volume (2000-2004 or 5-year average) from its initial condition ranged from 1.4 to 180 due to −30% to 30% changes in the CN (Figure 7a). The percent change in runoff volume (2000-2004 or 5-year average) from its initial condition ranged from 0.41 to 9.5 due to −10% to 30% changes in the CN for Row Crop (C + CR-Good) (Figure 7b). The percent change in runoff volume (2000-2004 or 5-year average) from its initial condition ranged from 9 to 82 due to −10% to 30% changes in the CN for Woods (Good) (Figure 7c).

Spatial Distribution of Runoff Depth
Since the AnnAGNPS model is able to provide landscape spatial variability by representing a watershed with a number of land areas (cells), we also evaluated the average annual runoff depth for all the watersheds (Figure 8)

Spatial Distribution of Runoff Depth
Since the AnnAGNPS model is able to provide landscape spatial variability by representing a watershed with a number of land areas (cells), we also evaluated the average annual runoff depth for all the watersheds (Figure 8

Model Performance to Estimate Event Peak Discharge
After the successful calibration and validation of runoff volumes, the model performance for the estimation of event peak discharge was done. For this purpose, we first looked at the gage height data for the selected USGS gauging stations. Then, we picked only the events when the gage height exceeded the existing flood stage (determined by USGS) for the particular station. The existing flood stage is 2.4 meters, 1.8 meters and 1.5 meters at USGS 03361650, USGS 04234000 and USGS 01117500, respectively. Figure 9 presents a graphical comparison between observed and model simulated event peak discharge (cubic meter per second) for only those selected events for all three watersheds. The model generally underestimated the peak discharge compared to few overestimations (calibration period) for Sugar Creek Watershed. On the other hand, in the case of Fall Creek Watershed, the model performed well to capture the highest peak discharge from tropical storm Lee in September 2011 based on the entire simulation period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013). Pawcatuck River Watershed showed a similar model performance. The record peak discharge in Rhode Island from the historic flood in March 2010 was captured by the model.

Discussion
AnnAGNPS performed satisfactorily for runoff prediction and simulation. Both daily and monthly calibration and validation statistics (NSE > 0.50 and RSR < 0.70, and PBIAS ± 25% based on [62]) of the developed AnnAGNPS model was satisfactory for runoff simulation for all the watersheds in the study. The model is categorized as satisfactorily performing when the range of the NSE value falls between 0.36 and 0.75 [66].
This study tested the applicability of the AnnAGNPS model on a glaciated landscape, which is why the three chosen watersheds are located in a glacial geomorphic setting. The area of the studied watersheds varied from 69 to 328 km 2 . After the evaluation of developed model performance, it was learned that the range of NSE varied from 0.51 to 0.57 and 0.60 to 0.83 for daily and monthly scale, respectively. This range is quite similar or even better than the 0.69 to 0.75 (on a monthly scale) found by [34] for a Mediterranean agricultural watershed (2.07 km 2 ) in Spain, 0.73 (on a monthly scale) by AnnAGNPS found by [26] for a 289.3 km 2 watershed in Illinois, USA, 0.65 (on a monthly scale) by SWAT and 0.48 to 0.58 (on a monthly scale) by AnnAGNPS found by [36] for a Mediterranean watershed (506 km 2 ) in Southern Italy, 0.53 (on a daily scale) by SWAT found by [67] for a large 1110 km 2 agricultural watershed in southwest France, 0.67 to 0.84 (on a monthly scale; calibration phase) found by SWAT for a large 4000 km 2 watershed in the North Carolina coastal plain [68], and the 0.53 to 0.62 (on a daily scale) found by SWAT for two watersheds in a semiarid region of Iraq [69]. Even, if we look at the value of R 2 from our study, it ranged from 0.54 to 0.66 (daily) and 0.64 to 0.86 (monthly). This range is also better than the one (0.50 to 0.80 by AnnAGNPS and 0.62 to 0.81 by SWAT on a monthly scale) found by [31] for an agricultural watershed in south-central Kansas. We also compared our developed AnnAGNPS model evaluation results with the results found from another water quality model, GLEAMS, in a study by [70] for agricultural watersheds in Indiana. The runoff calibration results were reported by [70] as NSE = 0.62 and R 2 = 0.70 for a monthly scale, which is quite similar to our results (NSE = 0.63 and R 2 = 0.67 for Sugar Creek; NSE = 0.60 and R 2 = 0.64 for Fall Creek; NSE = 0.68 and R 2 = 0.75 for Pawcatuck) for our studied watersheds.
Simulated runoff followed a similar trend (seasonal fluctuation) to observed runoff. In general, the model performed better in capturing event peak discharge for Fall Creek and Pawcatuck River watersheds rather than Sugar Creek Watershed. This poor performance of the model at the Sugar Creek Watershed could be improved by testing the effect of different storm types for rainfall distribution. As the regression coefficients for calculating the unit peak discharge are determined by storm type, the storm type within the AnnAGNPS model significantly influences peak discharge [71]. Simulated peak runoff was underestimated during some flood periods such as a major January 2005 flood event, record December 2013 flooding in Sugar Creek watershed, and a major April 2005 flood in Fall Creek watershed. In most cases, the model could not accurately simulate the flood runoff when the river overflowed. This characteristic of a hydrological model such as, AnnAGNPS is similar to that found in studies carried out by other studies (AnnAGNPS and SWAT performed poorly due to several extraordinary floods recorded during the study period in [36]; the SWAT model underestimated the largest flood in the study of [67]). However, the model was able to capture the historic September 2011 flood caused by tropical storm Lee in the case of Fall Creek watershed. In addition to that, the model even successfully captured the peak runoff volume (3815724.36 m 3 on 30 March 2010) generated from the historic 2010 flood in Rhode Island during the validation phase. Daily simulated runoff was also overestimated for some periods. Larger errors occurred when simulated peak runoff and average runoff differed significantly from the observed runoff volume. The spatial variability of the runoff depth ( Figure 8) could be attributed to the differences in land use, topography, soil type, and soil physical characteristics [58]. The output (spatial variability, i.e., cell wise runoff depth within the watershed) of the AnnAGNPS model contributes the riparian model field input data required for further research.

Conclusions
Model performance during calibration and validation phases shows that AnnAGNPS can be successfully used to predict runoff from watersheds in the glaciated settings of the Northeast and Midwest United States. This provides the capacity to couple edge-of-field hydrologic modeling with models that examine riparian or riverine functions and behaviors. The AnnAGNPS model effectively estimated runoff volume and portrayed the seasonal pattern of runoff in all the studied watersheds. The developed AnnAGNPS model could not capture some peak runoff events during wet periods and formed some unnecessary over-prediction of runoff during dry periods of the year. This characteristic of the model could be improved by the readjustment of the evaporation rate in association with the interception of precipitation events. The sensitivity analysis was limited to one specific contributing land use only, which could be extended by considering all possible combinations of CN based on the mixed land use of the watershed.
Author Contributions: M.T. was responsible for the conceptualization, methodology, data processing, software, calibration, validation, formal analysis, writing-original draft preparation, and editing; S.M.P. was responsible for the conceptualization, funding acquisition, methodology, model/software supervision, and writing-review and editing; A.J.G., K.A. and P.G.V. were responsible for the conceptualization, funding acquisition, and writing-review and editing; R.L.B. was responsible for model/software supervision, as well as writing-review and editing. All authors have read and agreed to the published version of the manuscript.