Effect of Watershed Delineation and Climate Datasets Density on Runoff Predictions for the Upper Mississippi River Basin Using SWAT within HAWQS

: The quality of input data and the process of watershed delineation can affect the accuracy of runoff predictions in watershed modeling. The Upper Mississippi River Basin was selected to evaluate the effects of subbasin and/or hydrologic response unit (HRU) delineations and the density of climate dataset on the simulated streamﬂow and water balance components using the Hydrologic and Water Quality System (HAWQS) platform. Five scenarios were examined with the same parameter set, including 8- and 12-digit hydrologic unit codes, two levels of HRU thresholds and two climate data densities. Results showed that statistic evaluations of monthly streamﬂow from 1983 to 2005 were satisfactory at some gauge sites but were relatively worse at others when shifting from 8-digit to 12-digit subbasins, revealing that the hydrologic response to delineation schemes can vary across a large basin. Average channel slope and drainage density increased signiﬁcantly from 8-digit to 12-digit subbasins. This resulted in higher lateral ﬂow and groundwater ﬂow estimates, especially for the lateral ﬂow. Moreover, a ﬁner HRU delineation tends to generate more runoff because it captures a reﬁned level of watershed spatial variability. The analysis of climate datasets revealed that denser climate data produced higher predicted runoff, especially for summer months. effects of different hydrologic response unit (HRU) and/or subbasin delineation schemes with SWAT. Over 60% of the studies were performed for watersheds in the United States, and the majority of that subset were reported for locations in the Corn Belt region. Other studies were conducted for watersheds in China, Ethiopia, Germany, and phosphorus represent different constituent forms reported in some studies; e.g., nitrate and mineral phosphorus [17]. e NR = not reported; NA = not applicable. f These studies included Soil and Water Assessment Tool (SWAT) conﬁgurations based on a grid approach or other non-HRU/subbasin methods. g These three studies report that sediment losses at the landscape (HRU) level were sensitive to shifts in HRU/subbasin delineations. h Three of these four river systems are portions of larger river basins (the 4776 km 2 system is the Maquoketa River Watershed); the sediment loss predicted at the outlet of the largest river system (17,491 km 2 ) was not sensitive to the delineation scenarios.


Introduction
Watershed models haves become an important technology to explore the effects of climate change and human activities on water resources and hydrological cycles [1][2][3]. The ability of a model to accurately represent the spatial variability of a catchment is the key element for the hydrological process. The accuracy required to achieve reliable simulation results relies on both the quality of the input data and also the watershed delineation in modeling [4][5][6].
The Soil and Water Assessment Tool (SWAT) ecohydrological model [7,8] is a popular watershed model and has proven to be an effective tool for evaluating agricultural management simulations for complex landscapes and varying climate regimes worldwide [9,10]. Studies have been conducted to evaluate the effects of watershed delineation on SWAT simulated streamflow and other outputs . Table 1 lists the majority of studies that have investigated the effects of different hydrologic response unit (HRU) and/or subbasin delineation schemes with SWAT. Over 60% of the studies were performed for watersheds in the United States, and the majority of that subset were reported for locations in the Corn Belt region. Other studies were conducted for watersheds in China, Ethiopia, Germany, India, South Korea, Tunisia and Turkey. The majority of studies reported that the delineation scenarios resulted in only minor or negligible effects on SWAT-predicted streamflow (Table 1). However, some studies described streamflow effects estimated with SWAT that were influenced by the various simulated delineation scenarios [11,15,20,22,[30][31][32]. Several of these studies noted that the effects on streamflow were relatively stable for more refined delineation schemes, but that the accuracy of the streamflow estimates declined for coarser HRU and/or subbasin subdivisions (e.g., [11,15,20]).
SWAT-predicted effects on sediment, nitrogen or phosphorus losses were reported to be sensitive to HRU and/or subbasin delineation effects in virtually all of the studies that described impacts on pollutant losses (Table 1). Sediment loss was found to be very sensitive at the landscape (HRU) level in some studies [13,14,16], indicating the transportlimited nature of the watershed and model structure issues that need to be considered in applications of SWAT. Chiang and Yuan [27] further assessed SWAT pollutant estimates as a function of both different HRU/subbasin subdivisions and stream order. They found that the impact of HRU delineation greatly increased as higher stream orders were used, and that the use of stream orders > 3 should be avoided for best results.
The size of systems analyzed in these studies ranged from 6.2 km 2 /6.3 km 2 to in Indiana/Germany to 17,941 km 2 in Iowa ( Table 1). The level of detail incorporated in these different assessments also varied considerably, with large ranges of HRUs and/or subbasins used in some studies (e.g., [13,14,16,27]) versus limited detail in other studies [12,31]. For small watersheds, the landscapes may not be complex enough to detect the runoff variation caused by different delineation schemes. It should be also noted that flow outputs were regarded as insensitive to different delineation schemes in some studies because sediment or nutrient outputs are much more sensitive in those studies [17,27]. Since a large basin has great variations of land-uses, soils or topography, it is necessary and important to evaluate the effects on the hydrologic cycle and to balance the level of model representation and computational limitations.
Watershed delineations were conducted in the Upper Mississippi River Basin (UMRB) in this study. Chen et al. [35] reviewed previous SWAT applications reported for the UMRB. Most of the SWAT applications for the UMRB relied on delineation approaches in which the subwatershed boundaries were aligned with what are commonly referred to as 8-digit watershed boundaries [36]; e.g., see previous SWAT-based studies reported by Jha et al. [37], Demissie et al. [38], Srinivasan et al. [39] and Qi et al. [40]. A refined delineation scheme has been incorporated into some reported UMRB SWAT models [41][42][43][44], which use subwatershed boundaries that are coincident with the much smaller USGS 12-digit watersheds [36]. Panagopoulos et al. [43] stated that the use of 12-digit watersheds should allow the model to more accurately capture key climate and topography data, relative to a coarser 8-digit delineation. They suggested that 12-digit watersheds are preferable for a large-scale system such as the UMRB compared to 8-digit watersheds and provided some other recommendations for large scale SWAT modeling. However, the hydrologic response of 8-digit and 12-digit watersheds have not been compared for the UMRB, and it is not clear how similar delineation scheme comparisons affect other river systems of comparable size. several of the studies did not report specific numbers of HRUs and a few also did not report numbers of subbasins. c Some studies reported separate calibration and validation periods; the simulation lengths reported here for those studies are based on the calibration period. d Nitrogen and phosphorus represent different constituent forms reported in some studies; e.g., nitrate and mineral phosphorus [17]. e NR = not reported; NA = not applicable. f These studies included Soil and Water Assessment Tool (SWAT) configurations based on a grid approach or other non-HRU/subbasin methods. g These three studies report that sediment losses at the landscape (HRU) level were sensitive to shifts in HRU/subbasin delineations. h Three of these four river systems are portions of larger river basins (the 4776 km 2 system is the Maquoketa River Watershed); the sediment loss predicted at the outlet of the largest river system (17,491 km 2 ) was not sensitive to the delineation scenarios.
In addition to watershed delineation, climate data is of ultimate importance to achieve a better match between observed and simulated streamflow. Since the climate spatial variability is the major driving force of hydrological process, water balance and subsequently water quality estimations across the basin, the density of climate observations can affect the simulation results. A comprehensive review of the literature addressing the relation between the density of rainfall observations and the accuracy with predicted streamflow can be found in Chaubey et al. [45] and Bárdossy and Das [46]. Chaubey et al. [45] assessed the variability induced in calibrated model parameters solely due to rainfall spatial variability. Bárdossy and Das [46] used different spatial resolutions of rainfall input in the model calibration and found that the model performance radically worsens with an excessive reduction of rain gauges. Given the fact that only one set of climate data can be input per subwatershed in SWAT, it is important to assess the model performance with denser or sparser climate data input.
In this study, SWAT models constructed by Hydrologic and Water Quality System (HAWQS) for UMRB were used to evaluate the impact of watershed delineation and climate data density on the runoff predictions. The objectives of this research include: (1) exploring the variations of land-use, soil and characterization of geometric properties due to watershed delineation; (2) evaluating the effect of watershed delineation on model predicted streamflow and on different water components of the hydrological process; and (3) evaluating the effect of climate data input with different densities on runoff predictions at spatial and temporal scales.

Model Description and Configuration
The development of the UMRB SWAT model was set up in Hydrologic and Water Quality System (HAWQS) platform [47]. HAWQS is a web-based interactive water quantity and quality modeling system that employs SWAT as its core modeling engine. Details of the SWAT modeling methods are described in the Theoretical Documentation [48]. HAWQS provides users with interactive web interfaces and pre-loaded input data and maps, including stream network, land use and land management, soil, historical climate, point sources, atmospheric deposition and reservoir data. The sources of these input data are listed in HAWQS [49]. Moreover, users can specify preferred values of parameters or climate data input instead of using the pre-defined sets, although some default values of parameters are preliminarily calibrated in HAWQS. It is quick and convenient for users to set up a SWAT project in HAWQS by only locating the ending Hydrologic Unit Code (HUC) of the watershed or river basin. Then, the watershed is subdivided into subbains/subwatersheds, and smaller drainage units referred to as hydrologic response units (HRUs) within each subbasin. For each HRU, spatially and temporally varying physical parameters are assumed to be homogeneous. At the subbasin level, the SWAT project provided by HAWQS can be conducted at 8-digit, 10-digit and 12-digit HUC [50]. At the HRU level, users can apply the HRU threshold to further discretize each subbasin considering land use and soil landscape heterogeneity.
The SWAT-based HAWQS simulations were performed from 1981 to 2005; the first 2 years served as an initialization period. Scenarios in this study were executed with the HAWQS default input parameters. It should be noted that the Hargreaves method was applied in the SWAT model to calculate ET, which was recommended for the UMRB [35]. The files created for UMRB SWAT model were also downloaded from HAWQS after the initial model construction, which allows additional parameter modification using the SWAT editor program or other external software.

Study Area
The UMRB is a headwater basin of the Mississippi River, and originates from Lake Itasca in northern Minnesota and outlets at the confluence of the Ohio and Mississippi Rivers near the town of Cairo in southern Illinois. It drains approximately 491,700 km 2 , Water 2021, 13, 422 6 of 23 which includes large portions of five states (Illinois, Iowa, Minnesota, Missouri and Wisconsin) and small portions of three other states (Indiana, Michigan and South Dakota). The basin outlet is often assumed to be a gauge site located near Grafton, Illinois, in SWAT modeling studies and is located just upstream of the confluence of the Mississippi and Missouri Rivers (Figure 1). Three USGS gauge sites were selected for model evaluation in this study: St. Paul (USGS 05331000), Clinton (USGS 05420500) and Grafton (USGS 05587450) with reported drainage areas of 95,312 km 2 , 221,703 km 2 and 443,665 km 2 , respectively. Observed flow data measured at the St. Paul, Clinton and Grafton stations were obtained from the United States Geological Survey (USGS) Water Data Server over the same period from 1980 to 2005 [51]. The major land use in the UMRB is cropland (44.7%), which is dominated by rotations of corn (27.5%) and soybean (17.2%). There are smaller areas of wheat, oats and other crops in the UMRB but those are excluded in this HAWQS modeling framework. Other important land-use categories include forest (20.2%), grassland (16.2%), water and wetlands (9.8%), and urban/developed areas (9.1%). The soil types range from heavy, poorly drained clay soil to light, well-drained sands. Among those soil types, silty loam and loam soils cover about 66% of the total UMRB area [52]. The topography is characterized by flat to gently rolling terrain, with 55% of the area having less than a 2% slope and an average elevation of 280 m. Major spatial maps applied in the model are presented in Figure 2, including a (a) Digital elevation model (DEM) map, (b) land-use map and (c) soil map. Topographic information is provided in the form of 30 m digital elevation model (DEM) data from the National Elevation Dataset at a resolution of 1 arc second [53]. Agricultural and non-agricultural land use data were derived from the Cropland Data Layer [54] and National Land Cover Database [55]. The soil data comes from the State Soil Geographic (STATSGO) database, which contains soil maps at a 1:250,000 scale [56].  Major spatial maps applied in the model are presented in Figure 2, including a (a) Digital elevation model (DEM) map, (b) land-use map and (c) soil map. Topographic information is provided in the form of 30 m digital elevation model (DEM) data from the National Elevation Dataset at a resolution of 1 arc second [53]. Agricultural and nonagricultural land use data were derived from the Cropland Data Layer [54] and National Land Cover Database [55]. The soil data comes from the State Soil Geographic (STATSGO) database, which contains soil maps at a 1:250,000 scale [56].

Watershed Delineation
The number of subbasins configured for a SWAT application is normally decided by the minimum drainage area (MDA) of the river network. The smaller the MDA is, the denser the river network becomes and the more subbasins are generated. Instead of MDA, the delineation of the UMRB into subbasins in this study was based on 8-digit Hydrologic Unit Codes (HUCs) and 12-digit HUCs from the Watershed Boundary Dataset (WBD). The WBD maps the full areal extent of surface water drainage for the U.S. using a hierarchical system of nesting hydrologic units at various scales [36]. The WBD contains eight levels of progressive hydrologic units identified by unique 2-to 16-digit codes. The number of digits for a respective HUC indicates the level of classification in the hydrologic unit system from the largest geographic area to the smallest geographic area. HUC 8 is analogous to medium-sized river basins, while HUC 12 is a more local sub-watershed level that captures tributary systems ( Figure 3). The National Hydrography Dataset Plus (NHDPlus) is a national geospatial surface water framework, which is also used to provide the vector stream network connecting subbasins hydrologically in HAWQS [57]. Major spatial maps applied in the model are presented in Figure 2, including a (a) Digital elevation model (DEM) map, (b) land-use map and (c) soil map. Topographic information is provided in the form of 30 m digital elevation model (DEM) data from the National Elevation Dataset at a resolution of 1 arc second [53]. Agricultural and non-agricultural land use data were derived from the Cropland Data Layer [54] and National Land Cover Database [55]. The soil data comes from the State Soil Geographic (STATSGO) database, which contains soil maps at a 1:250,000 scale [56].

Watershed Delineation
The number of subbasins configured for a SWAT application is normally decided by the minimum drainage area (MDA) of the river network. The smaller the MDA is, the denser the river network becomes and the more subbasins are generated. Instead of MDA, the delineation of the UMRB into subbasins in this study was based on 8-digit Hydrologic Unit Codes (HUCs) and 12-digit HUCs from the Watershed Boundary Dataset (WBD). The WBD maps the full areal extent of surface water drainage for the U.S. using a hierarchical system of nesting hydrologic units at various scales [36]. The WBD contains eight levels of progressive hydrologic units identified by unique 2-to 16-digit codes. The number of digits for a respective HUC indicates the level of classification in the hydrologic unit system from the largest geographic area to the smallest geographic area. HUC 8 is analogous to medium-sized river basins, while HUC 12 is a more local sub-watershed level that captures tributary systems ( Figure 3). The National Hydrography Dataset Plus (NHDPlus) is a national geospatial surface water framework, which is also used to provide the vector stream network connecting subbasins hydrologically in HAWQS [57].
After the subbasins are defined in SWAT, those subbasins are apportioned into HRUs by setting thresholds of land-use and soil type to capture relatively small areas. Two levels of HRU thresholds (1 km 2 /1 km 2 and 70%/70%) were tested in this study. For a large basin like UMRB, the threshold of 1 km 2 is expected to be a refined representation of the watershed. Conversely, considering the total area of URMB, the 70% threshold is definitely much coarser than 1 km 2 . Thus, the number of HRUs for 70% threshold is expected to shrink significantly, which means the comparison between thresholds of 1 km 2 and 70% is more like the comparison between multiple HRUs and dominant HRUs.

Spatial Climate Datasets
The PRISM (Parameter-elevation Relationships on Independent Slopes Model) dataset was developed by Oregon State University's PRISM Climate Group [58]. The highspatial-resolution climate map products of PRISM, created at 30-arcsec (~800 m) grid resolution, serve as the official spatial climate datasets of the U.S. Department of Agriculture. PRISM is an analytical tool that uses point data, a digital elevation model, other spatial datasets, and an encoded spatial climate knowledge base to generate gridded estimates of climate elements, such as precipitation and temperature [59].
In HAWQS, original PRISM data were processed by averaging their gridded cell values within 8-digit watersheds or 12-digit watersheds. The mean values of the processed climate data were assigned to the virtual gauge sites, which are located in the geometrical center of each sub-watershed. Therefore, the density of climate data is directly affected by the number of subbasins. To compare the effect of spatial density of climate data on runoff, two climate datasets were used in this study ( Table 2). The "PRISM" dataset represents After the subbasins are defined in SWAT, those subbasins are apportioned into HRUs by setting thresholds of land-use and soil type to capture relatively small areas. Two levels of HRU thresholds (1 km 2 /1 km 2 and 70%/70%) were tested in this study. For a large basin like UMRB, the threshold of 1 km 2 is expected to be a refined representation of the watershed. Conversely, considering the total area of URMB, the 70% threshold is definitely much coarser than 1 km 2 . Thus, the number of HRUs for 70% threshold is expected to shrink significantly, which means the comparison between thresholds of 1 km 2 and 70% is more like the comparison between multiple HRUs and dominant HRUs.

Spatial Climate Datasets
The PRISM (Parameter-elevation Relationships on Independent Slopes Model) dataset was developed by Oregon State University's PRISM Climate Group [58]. The high-spatialresolution climate map products of PRISM, created at 30-arcsec (~800 m) grid resolution, serve as the official spatial climate datasets of the U.S. Department of Agriculture. PRISM is an analytical tool that uses point data, a digital elevation model, other spatial datasets, and an encoded spatial climate knowledge base to generate gridded estimates of climate elements, such as precipitation and temperature [59].
In HAWQS, original PRISM data were processed by averaging their gridded cell values within 8-digit watersheds or 12-digit watersheds. The mean values of the processed climate data were assigned to the virtual gauge sites, which are located in the geometrical center of each sub-watershed. Therefore, the density of climate data is directly affected by the number of subbasins. To compare the effect of spatial density of climate data on runoff, two climate datasets were used in this study ( Table 2). The "PRISM" dataset represents the original PRISM data processed for 8-digit subbasins, so there are 119 virtual gauge sites in the study area. The alternative "PRISM*" dataset represents the original PRISM data processed for 12-digit subbasins, resulting in 5163 virtual gauge sites in the study area. Daily precipitation and temperature data from 1981 to 2005 were used from both climate datasets to simulate streamflow in this study.

Simulation Scenario
A total of five scenarios of various combinations of subbasins, HRU definitions and climate data were simulated in this study (Table 2), including (1) two levels of subbasins (8-digit and 12-digit HUCs), (2) two levels of HRU thresholds (1 km 2 for both land-use and soil, 70% for both land use and soil) and (3) two levels of climate data densities (PRISM and PRISM*). In the HUC8 scenario, the SWAT model was configured for the UMRB at the 8digit watershed level within the HAWQS Version 1.1 platform (https://hawqs.tamu.edu/#/ (accessed on 10 February 2020)), resulting in 119 8-digit subbasins that encompass the previously described 447,802 km 2 area that drains to Grafton, IL (the outlet is the 8-digit watershed identified as HUC07110009). HRU thresholds of 1 km 2 were then applied to the land use and soil type in each subbasin. The application of the thresholds resulted in a total of 30,812 HRUs. Climate data processed for 8-digit subbasins, namely "PRISM", was selected as input.
When the SWAT model was configured for the UMRB at the 12-digit watershed level within the HAWQS Version 1.1 platform, there are 5163 12-digit subbasins with an area of 445,358 km 2 draining to Grafton, IL (the outlet is located in the 12-digit watershed at HUC071100090401). It should be noted that the total area of 12-digit subbasins is not identical with the 8-digit subbasins total area. Some closed-basins are excluded in the drainage areas to Grafton at the 12-digit subbasins level. This bias is caused by the hydrology dataset provided by the WBD and is acceptable for such a large basin like the UMRB. The number of subbasins delineated by 12-digit HUCs for UMRB is about 43 times greater than the number of subbasins delineated by 8-digit HUCs. With the threshold of 1 km 2 , the 12-digit watershed model generated 120,454 HRUs for the UMRB, which is four times as many as the threshold of 70%. Generally, as the HRU threshold increases, fewer HRUs are subdivided in a subbasin as greater portions of minor land use and soil were regrouped into the major HRUs. This is the reason why the HUC12 and HUC12* scenarios with 70% threshold only have 28,823 HRUs. Even though the HUC8 scenario and the New-HUC12 scenario have the identical HRU threshold during the process of HRU definition, the number of HRUs in the 12-digit watershed is about four times as many as the 8-digit watershed. This indicates that the number of HRUs can be increased significantly with the increase in subbasins for a large basin. When simulating, all of the delineation scenarios adopted the same set of parameters, which were preliminarily calibrated by HAWQS.

Land-Use, Soil and Landscape Discrepancies from Different Delineations
Delineating a particular modeled area has an impact on both terrain complexity and topographic attributes. These characteristics are related to the input elevation, soils and land use data. Land-use distributions from various delineations are presented in Table  3 to show the discrepancies among them. The original scenario refers to the true GIS land-use distribution with a 0% threshold. As shown in Table 3, more minor land uses were regrouped into major land uses as the HRU thresholds become greater. The landuse distribution of the HUC8 scenario is the closest to the original distribution. For the New-HUC12 scenario, changes in all land uses are as small as no more than 1% at Grafton and St. Paul. Land-use distribution of the HUC12 scenario deviated significantly from the original situation at three gauges due to the 70% HRU threshold. Cropland increased from 40.67% to 47.50% at St. Paul, from 29.42% to 30.80% at Clinton and from 44.69% to 50.05% at Grafton. The percentage of the forest also increased at all three gauges in the HUC12 scenario compared to the original situation. Forest increased by 4.80%, 8.68% and 5.07% at St. Paul, Clinton and Grafton, respectively. Urban area decreased significantly on average by 6% for three gauge sites. It should be noticed that even though the variation of land use is small, the redistribution areas can be great since the total drainage areas are large in this study. St. Paul, Clinton and Grafton are reported to control drainage areas of 95,312 km 2 , 221,703 km 2 and 443,665 km 2 , respectively. There are more than 160 soil types identified for the whole UMRB with the STATSGO database. The predominant soil types in the watershed cover 3.64%, 2.48% and 2.06%, respectively. In addition, 16 other types of soil cover >1% of the UMRB land area. Hydrologic soil groups (HSGs) serve an important role in the determination of surface runoff [60]. Soils are classified into HSGs by national agencies based on the minimum rate of infiltration obtained for bare soil after prolonged wetting [61]. Soils in group A are characterized by the low runoff potential and high infiltration rates when thoroughly wet, while soils in group D are characterized by high runoff potential and very low infiltration rates when thoroughly wet. Dual hydrologic soil groups (i.e., A/D, B/D and C/D) are placed in group D based solely on the presence of a high water table. Since there are too many soil types in the watershed, only hydrologic groups among different delineation scenarios are listed (Table 4). When applying a 0% threshold for the soils, 69.6% of the area belongs to group B, which represents moderately low runoff potential and moderate infiltration rates. Group C soils form the second largest portion of soil types in the watershed, which accounts for 16.7% of the area. Group A covers 8.7% and group D covers 5.0% in the original situation. The compositions of hydrologic soil groups for HUC-8 scenario and New-HUC12/New-HUC12* scenario, which both used 1 km 2 HRU threshold, are almost the same as the original situation ( Table 4). The biggest deviation from the original situation occurred in the HUC-12/HUC-12* scenario when applying 70% threshold, with group A 7.7%, group B 69.7%, group C 17.1% and group D 5.5%. The composition of the hydrologic soil groups varied little among different delineation scenarios. Differences between the four hydrologic groups were <1% for the whole UMRB versus the original land use distribution, for the threshold set to 1 km 2 for either the HUC8 scenario or the New-HUC12 scenario. Soils in group A decreased from 8.67% to 7.70% compared to the original land use distribution in response to the 70% threshold applied for the HUC12 scenario, which was the largest discrepancy (Table 4). The individual subbasin areas varied in size within each subdivided watershed, as shown in Figure 3. Descriptive statistics of subbasin sizes are presented in Table 5. Areas of 8-digit subbasins varied from 1606.1 km 2 to 8415.9 km 2 , while areas of 12-digit subbasins varied from 18.0 km 2 to 940.5 km 2 . The number and size of subbasins have a direct impact on the river or main channel simulated for each subbasin. Table 5 summarizes the physical characteristics of channels that affect water flow, and transport of sediment and nutrients. The average length of the channel decreased from 189.3 km to 14.9 km, while the total length of channels increased two times in 12-digit subbasins versus 8-digit subbasins. The average width of a channel decreased sharply from 177.1 m to 18.4 m. There was an obvious increase in the slope of the channel from 0.2% to 0.5% when the number of subbasins increased. The maximum slope of the channel in each subbasin also increased from 2.8% to 4.6%. This increase in slope could result from a better accounting of spatial variation in elevation when smaller subbasins were used, which is a fundamental property of natural terrain [62]. Studies demonstrated that the drainage density can affect the accuracy of runoff predictions since it characterizes the scale of landscape forms [63][64][65][66]. Table 5 shows that a wider range in drainage densities of 0.002 to 0.497 was observed for the 12-digit subbasins, as compared to a drainage density range of 0.005 to 0.145 for the 8-digit subbasins. For the whole study area, drainage density increased from 0.05 to 0.17. Table 5 shows that drainage density increased as the number of subbasins increased. A similar trend was also reported by Jha et al. [17]. Higher drainage density in the simulation process implies that scenarios based on the 12-digit subbasin may be characterized by a stronger tendency to generate runoff.

Evaluation of Monthly Streamflow at Gauge Sites
The model performance for monthly streamflow at St. Paul, Clinton and Grafton was presented for all scenarios in terms of the Nash-Sutcliffe Efficiency (NSE), percent bias (PBIAS), R 2 and Kling-Gupta efficiency (KGE) statistics in Table 6. NSE is a normalized statistic that determines the relative magnitude of the residual variance compared to the measured data variance, which is widely used throughout the modeling literature [67,68]. PBIAS measures the average tendency of the simulated data to be larger or smaller than their observed counterparts. R 2 describes the degree of collinearity between simulated and measured data. As for KGE, it is designed to be an improved criterion by incorporating error compensation for the bias and variability components [69]. Generally, NSE values ≥ 0.5, R 2 values ≥ 0.6 and PBIAS values ≤ ±25% [67] or ≤±15% [68] can be regarded as satisfactory. Patil and Stieglitz [70] implied that simulated values could be judged to be satisfactory with a KGE value > 0.6.
Since HAWQS was preliminarily calibrated based on several 8-digit subbasins, monthly streamflow performance at the St. Paul gauge was satisfactory for the HUC8 scenario with NSE 0.66, PBIAS 12.23%, R 2 0.74 and KGE 0.79 based on the suggested criteria by Moriasi et al. [67] or Patil and Stieglitz [70]. NSE at Clinton was 0.32 for the HUC8 scenario, which was the weakest among three gauge sites. PBIAS at Clinton was 29.82% (>25%) for the HUC8 scenario, while R 2 was 0.76 (>0.6) and KGE was 0.57 (<0.6). The uncalibrated flow performance at Grafton can be considered acceptable for the HUC8 scenario shown in Table 6, except that PBIAS was 30.32% (>25%). Overall, the uncalibrated flow performance at St. Paul and Grafton can be considered acceptable for the HUC8 scenario shown in Table 6.
The statistics of monthly streamflow at the St. Paul gauge site were satisfactory for all 12-digit subbasin scenarios ( Table 6) in terms of NSE (0.59-0.64), PBIAS (2.85-16.77%), R 2 (0.68 to 0.71) and KGE (0.77 to 0.80). PBIAS decreased slightly in most of the 12-digit subbasin scenarios versus the HUC8 scenario at Clinton and Grafton, which indicated an improvement of bias between observed and simulated streamflow. KGE values were considerably strong at the URMB outlet of Grafton, as evidenced by values that were close to 0.6 for all five scenarios. However, the values of NSE and R 2 declined considerably at the Clinton and Grafton gauges when using 12-digit subbasin scenarios (Table 6). NSE values decreased to <0 in all 12-digit subbasin scenarios at Clinton and decreased to <0.5 at Grafton. R 2 values decreased to <0.6 at Clinton and Grafton in the 12-digit subbasin scenarios. The shift from 8-digit to 12-digit subbasins for the drainage area to Clinton resulted in the worst simulation of monthly flow among the three gauge sites in this study (Table 6). This is partly due to the weak representation of forest from the SWAT parameterization, which has been reported in previous SWAT applications [71,72]. Since the area draining to Clinton has the largest percentage of forest (29.8%) among three gauge sites, this weakness could be magnified with finer delineation.
The last column in Table 6 provides the average monthly flow on an annual basis for the different delineation scenarios at the three gauge sites. The average monthly flow is underestimated for all of the scenarios compared to the observed flow data, which is consistent with the positive PBIAS values. The average monthly flow increased for most of the 12-digit subbasin scenarios as compared to the HUC8 scenario, except for the HUC12 scenario at St. Paul and Clinton. The increase in flow is attributable to variation in the number of HRUs and corresponding variation in the characteristics of HRUs. Since the 12-digit model has considerably more HRUs than the 8-digit model, more detailed information on the spatial variability of rainfall is captured. Although HUC8 scenario and New-HUC12 scenario both used 1 km 2 as the HRU threshold, the New-HUC12 scenario generated 120,454 HRUs, 4 times greater than the HUC8 scenario of 30,812 HRUs. The average monthly flow increased from the HUC8 scenario to the New-HUC12 scenario by 2.22%, 1.79% and 7.17% at St. Paul, Clinton, and Grafton, respectively. The New-HUC12 scenario of 12,454 HRUs had higher average monthly flows compared to the HUC12 scenario of 28,823 HRUs. Average monthly flows increased from HUC12 scenario to New-HUC12 scenario by 7.79%, 6.47% and 3.80% at St. Paul, Clinton and Grafton, respectively. Similarly, streamflow increased from the HUC12* scenario to New-HUC12* scenario. Increases in streamflow from coarse to fine watershed delineation were also reported in previous studies [13,22]. HRU modifications that affect the distribution of simulated land use, soils, and other landscape characteristics will have the greatest impact on the predicted runoff [17]. The monthly streamflow decreased from the HUC8 scenario to the HUC12 scenario at St. Paul and Clinton, which might be due to the smaller number of HRUs for the HUC12 scenario, resulting from the 70% threshold (Table 2). In general, there will be less information loss as the HRU definition becomes more refined. It is perceived that this will result in more accurate model prediction [28]. However, this may not always be the better outcome of statistics, especially regarding SWAT-predicted hydrologic results.
Time series plots of observed versus simulated streamflow on the monthly time scale for the three gauge sites (Figure 1) in response to different delineation scenarios are presented in Figure 4. The solid blue dots represent the measured monthly streamflow obtained from USGS [50]. The simulated SWAT model streamflow generally tracked the seasonal variance pattern including the peaks and recessions, although there is an obvious underestimation of the observed streamflows for all three gauge sites (St. Paul, Clinton and Grafton). The St. Paul station (Figure 4a) demonstrates strong agreement between simulated and observed monthly flows, which is consistent with the statistics criteria presented in Table 6. Winter low flows were generally underpredicted during November to February, especially for Clinton. Based on the enlarged picture of winter flow shown in Figure 4, the HUC8 scenario (represented by the solid black line with hollow black circles) usually had the lowest flow among the different scenarios, again confirming that the 12-digit watersheds generate higher flow than 8-digit watersheds.   Table 7 presents average annual values of hydrologic components predicted with SWAT across the UMRB for the five different delineation scenarios, including precipitation, evapotranspiration (ET), potential evapotranspiration (PET), surface runoff, lateral runoff, groundwater flow and water yield. Annual values are averaged for the whole UMRB during the period from 1983 to 2005. The results show that higher lateral flow and groundwater flow estimates occurred for the 12-digit subbasin scenarios relative to the 8digit subbasin simulations. This was especially true for lateral flow, as reflected by an increase of nearly 100%. Groundwater flow for the HUC8 scenario was 18.3 mm, while it ranged from 18.6 to 19.2 mm in all 12-digit subbasin scenarios. The lateral flow for the HUC8 scenario was 24.0 mm, smaller than the 12-digit subbasin scenarios (45.0-47.9 mm), reflecting the greatest difference of any of the estimated hydrologic components. Increases in these components revealed that transmission gains from shallow groundwater increased as the subbasins decreased in size, which could be due to increases in drainage density ( Table 5). The spatial variation in elevation can be better accounted for as drainage density increases [17]. Surface runoff and water yield decreased in the HUC12 and New-HUC12 scenarios as compared to the HUC8 scenario but increased in the HUC12* and New-HUC12* scenarios. Predicted sediment yields also increased, with a greater number  Table 7 presents average annual values of hydrologic components predicted with SWAT across the UMRB for the five different delineation scenarios, including precipitation, evapotranspiration (ET), potential evapotranspiration (PET), surface runoff, lateral runoff, groundwater flow and water yield. Annual values are averaged for the whole UMRB during the period from 1983 to 2005. The results show that higher lateral flow and groundwater flow estimates occurred for the 12-digit subbasin scenarios relative to the 8-digit subbasin simulations. This was especially true for lateral flow, as reflected by an increase of nearly 100%. Groundwater flow for the HUC8 scenario was 18.3 mm, while it ranged from 18.6 to 19.2 mm in all 12-digit subbasin scenarios. The lateral flow for the HUC8 scenario was 24.0 mm, smaller than the 12-digit subbasin scenarios (45.0-47.9 mm), reflecting the greatest difference of any of the estimated hydrologic components. Increases in these components revealed that transmission gains from shallow groundwater increased as the subbasins decreased in size, which could be due to increases in drainage density ( Table 5). The spatial variation in elevation can be better accounted for as drainage density increases [17]. Surface runoff and water yield decreased in the HUC12 and New-HUC12 scenarios as compared to the HUC8 scenario but increased in the HUC12* and New-HUC12* scenarios. Predicted sediment yields also increased, with a greater number of subbasins represented by the 12-digit scenarios. This is due to changes in channel length and slope, which affect the deposition and degradation of sediments. To investigate the spatial variation of hydrologic components, Figure 5 depicts outputs for the HUC8 scenario versus the New-HUC12 scenario, respectively. The results are presented at the respective subbasin scale for both scenarios. Considerably more spatial details are observable in the 12-digit subbasins compared to 8-digit subbasins. A distinct area with a long ribbon shape came out as compared to the HUC8 scenario, showing lower surface runoff (Figure 5a). Surface runoff is related to soil type, where soils with higher infiltration rates (A and B hydrologic soil groups) tend to generate lower runoff. This indicated that more soil characteristics were captured by 12-digit subbasins than 8-digit subbasins when delineating the watershed. The steeper areas in the northern and northeast parts of the basin showed higher lateral flow than the flatter areas in the middle subbasins (Figure 5b). Unlike surface runoff, groundwater flow was higher in regions of higher infiltration (Figure 5c). This is why the ribbon-shaped area showed lower surface runoff but higher groundwater flow within the 12-digit subbasins (Figure 5a,c). With respect to water yield, it is a composite value related to surface runoff, lateral flow, groundwater flow, etc. Water yield did not show obvious spatial characteristics in the comparison between 8-digit and 12-digit subbasins (Figure 5d). Overall, the results demonstrate that not only does the landscape position matter in simulation, but that the spatial distribution of land use and soils is also important.

Response of Hydrologic Components
Many SWAT studies report analyses of the effects of land-use change [73], and dozens of these studies describe how simulated streamflow varies as a function of land-use change [74][75][76]. The shifts in average annual values of key hydrologic components in response to different land-use types simulated during 1983 to 2005 are presented in Figure 6. The highest PET and ET levels were predicted for the water bodies; i.e., lakes or ponds. ET estimated for forest was greater than the other non-water land uses. In the SWAT hydrologic process, water areas do not generate runoff to the streamflow. Cropland and urban areas were predicted to produce more surface runoff than other land-uses. The average annual amount of surface runoff for the forest was the lowest. Differences in the surface runoff between land-uses were primarily caused by the runoff curve number (RCN), a factor used to estimate surface runoff in the SCS curve number procedure [77]. The forest curve number was lower than cropland, which reflected the effects of the forest canopy and ground cover, and resulted in a decrease in surface runoff. Wetland was predicted to generate the largest lateral flow in all of the scenarios, while grassland was predicted to produce the highest groundwater flow. Lateral flow was obviously greater for the HUC12 scenarios relative to the HUC8 scenario for all land uses. This is in accordance with the variation of lateral flow for the whole UMRB in Table 7. Among all land use types, the predicted water yield was lowest for forest. In summary, comparisons among the different land uses indicated that forest will consistently result in the lowest levels of generated runoff, which was also reported in Owuor et al. [78] and Shang et al. [79].

Effect of the Spatial Density of Climate Dataset on Runoff
In this section, the effects of spatial density of the climate datasets (including precipitation and temperature) on hydrologic processes are presented. Two climate datasets were used in this study ( Table 2). The PRISM dataset was referred to as a processed dataset, which was averaged from original gridded-PRISM data to the 8-digit watersheds. Similarly, the PRISM* dataset provided average daily climate data at the 12-digit watershed level. More descriptions of the climate datasets are shown in Section 2.3. The PRISM* dataset resulted in increases of SWAT-predicted streamflow relative to the PRISM climate data on an average annual basis at the three gauge sites (Table 6). From the HUC12* scenario to the HUC12 scenario, monthly streamflow increased by 8.36%, 6.98% and 6.26% at St. Paul, Clinton and Grafton, respectively. Similarly, average monthly streamflow increased by 8.29%, 6.86% and 6.34% at the three gauge sites. On average, monthly streamflow increased by 8.32%, 6.92% and 6.30% at St. Paul, Clinton and Grafton, respectively, from denser climate data to sparser climate data.
From the view of the entire basin, differences in average annual precipitation and daily temperature are quite small (Table 7). Among the five scenarios, average annual precipitation ranged from 830.7 mm to 831.2 mm, and the daily mean temperature ranged The investigation of land use showed that cropland resulted in a relatively overall stronger water production capacity (Figure 6c,e,f). Although the change of land-use percentage is small, the total changed land area can be considerable for a large basin like the UMRB. Thus, the streamflow performance for different stations can be very different in response to the impact of the division of the basin into subbasins and accompanying HRUs. To conclude, the changes in land use and soil characteristics at the HRU level, due to the implementation of HRU thresholds, can have a great influence on simulated surface runoff. Average annual flow decreased at the St. Paul and Clinton gauges but increased at Grafton station for the HUC12 scenario compared to the HUC8 scenario. This may be partly due to the fact that land-use compositions are different for the areas draining to three gauges (Table 3). To be specific, the percentage of cropland at Grafton is the highest among three gauges.

Effect of the Spatial Density of Climate Dataset on Runoff
In this section, the effects of spatial density of the climate datasets (including precipitation and temperature) on hydrologic processes are presented. Two climate datasets were used in this study ( Table 2). The PRISM dataset was referred to as a processed dataset, which was averaged from original gridded-PRISM data to the 8-digit watersheds. Similarly, the PRISM* dataset provided average daily climate data at the 12-digit watershed level. More descriptions of the climate datasets are shown in Section 2.3. The PRISM* dataset resulted in increases of SWAT-predicted streamflow relative to the PRISM climate data on an average annual basis at the three gauge sites (Table 6). From the HUC12* scenario to the HUC12 scenario, monthly streamflow increased by 8.36%, 6.98% and 6.26% at St. Paul, Clinton and Grafton, respectively. Similarly, average monthly streamflow increased by 8.29%, 6.86% and 6.34% at the three gauge sites. On average, monthly streamflow increased by 8.32%, 6.92% and 6.30% at St. Paul, Clinton and Grafton, respectively, from denser climate data to sparser climate data.
From the view of the entire basin, differences in average annual precipitation and daily temperature are quite small (Table 7). Among the five scenarios, average annual precipitation ranged from 830.7 mm to 831.2 mm, and the daily mean temperature ranged from 8.01 • C to 8.03 • C. In the comparisons of HUC12 versus HUC12* scenario and New-HUC12 versus New-HUC12* scenario, ET decreased by 9.4 mm and 9.5 mm from the sparser dataset (PRISM) to the denser dataset (PRISM*). Differences in PET between the two climate datasets are increments of 2.6 mm and 2.9 mm (Table 7). Surface runoff, lateral flow, groundwater flow and water yield also increased to varying degrees between the results obtained with PRISM versus PRISM*. Among these hydrologic components, larger increases were predicted for surface runoff and water yield as compared to lateral runoff and groundwater flow. Surface runoff increased by 9.9 mm for the HUC12 scenario relative to the HUC12* scenario and 10.7 mm from the New-HUC12 scenario to the New-HUC12* scenario. With respect to water yield, increases between the two sets of scenarios were up to 11.7 mm and 12.1 mm, respectively. In contrast, lateral runoff only increased by 1.1 mm and 0.8 mm, and groundwater increased by just 0.2 mm (in both cases), for the scenarios executed with PRISM*-versus the PRISM -based scenarios. Figures 7a and 8a show the spatial variations in the precipitation and daily mean temperature for the 1983 to 2003 period for the HUC12 scenario with the PRISM dataset, which had 119 virtual gauges. Average annual precipitation was found to vary widely across the study area from <600 mm in the northwest part of the basin to >1000 mm in the southeast area of the basin. The average mean temperature varied from <4 • C in the north to >12 • C in the south. Figures 7b and 8b depict the absolute changes that were predicted for the PRISM dataset relative to PRISM* dataset in precipitation (mm) and mean temperature ( • C). There was no distinct geographic pattern regarding the absolute changes of annual precipitation between the denser dataset and sparser dataset. Highs and lows of absolute changes are distributed in the study area randomly. In some 12digit subbasins, annual precipitation with PRISM* was <40 mm versus PRISM. In other subbasins, annual precipitation with PRISM* was 40 mm > than PRISM. This can be explained by the calculation method used in SWAT for the precipitation and temperature data input, which was weighted on an areal basis using Thiessen polygons. Thiessen polygons are irrelevant to the geographical characteristics but are related to the number and size of subbasins. Similarly, differences in daily mean temperature between the two datasets did not show an obvious relationship with land-use, soil or landscape.
In addition to comparing the spatial distribution of precipitation and temperature, mean monthly values for selected hydrologic components across the UMRB in response to the different scenarios are presented in Figure 9. As shown in Figure 9a, precipitation amounts across the different scenarios were very close, which demonstrated that precipitation changed little both on an annual basis and monthly basis as a function of different densities of climate datasets. The PRISM* dataset generated lower ET as compared to the PRISM dataset across the year from January to December (Figure 9b). With respect to surface runoff and water yield, the higher levels generated by the PRISM* data are very obvious, especially during the May to September growing season (Figure 9c,d). The cumulative difference between PRISM and PRISM* during the growing season accounted for 75% of the total annual difference on average. To conclude, there are distinct seasonal variations between the monthly differences. The PRISM* dataset, which has denser gauges, tended to generate more runoff during the summer months, as compared to the PRISM dataset. Figure 9 further shows that the New-HUC12* scenario resulted in the highest runoff among the different scenarios in most months.
annual precipitation between the denser dataset and sparser dataset. Highs and lows of absolute changes are distributed in the study area randomly. In some 12-digit subbasins, annual precipitation with PRISM* was <40 mm versus PRISM. In other subbasins, annual precipitation with PRISM* was 40 mm > than PRISM. This can be explained by the calculation method used in SWAT for the precipitation and temperature data input, which was weighted on an areal basis using Thiessen polygons. Thiessen polygons are irrelevant to the geographical characteristics but are related to the number and size of subbasins. Similarly, differences in daily mean temperature between the two datasets did not show an obvious relationship with land-use, soil or landscape.   In addition to comparing the spatial distribution of precipitation and temperature, mean monthly values for selected hydrologic components across the UMRB in response to the different scenarios are presented in Figure 9. As shown in Figure 9a, precipitation amounts across the different scenarios were very close, which demonstrated that precipitation changed little both on an annual basis and monthly basis as a function of different densities of climate datasets. The PRISM* dataset generated lower ET as compared to the PRISM dataset across the year from January to December (Figure 9b). With respect to surface runoff and water yield, the higher levels generated by the PRISM* data are very obvious, especially during the May to September growing season (Figure 9c,d). The cumu-

Conclusions
The SWAT model was developed for the UMRB using the data provided by the HAWQS platform to analyze the hydrologic response to watershed delineations and the density of climate data. Five delineation scenarios combined with two levels of subbasins (8-digit HUC and 12-digit HUC), two levels of HRU thresholds (1 km 2 /1 km 2 and 70%/70%) and two levels of climate data (PRISM and PRISM*) were simulated. The same set of default parameters generated within HAWQS was used for all five delineation scenario simulations, which included some "partial built-in calibration" for some UMRB 8digit subbasins in HAWQS. The results of the analyses lead to the following conclusions.
1. The hydrologic response can be varied in different parts of the area for a large basin when changing delineation schemes, since the watershed delineation affects the landscape position, land-use distribution and soil distribution spatially. Results showed that groundwater flow and surface runoff are more sensitive to the areas with soils of higher infiltration rates, and the steeper areas generated a greater amount of lateral runoff than flat areas when applying denser delineation schemes. The statistics of monthly streamflow at St. Paul were satisfactory for all scenarios. However, several evaluation criteria of monthly streamflow declined considerably at the Clinton and Grafton gauges when using 12-digit subbasin scenarios.
2. The hydrologic cycle changed in the shift from 8-digit to 12-digit subbasin. The average length of channels and width of channels decreased in the 12-digit subbasins versus 8-digit subbasins, while the average slope of channels and drainage density increased.

Conclusions
The SWAT model was developed for the UMRB using the data provided by the HAWQS platform to analyze the hydrologic response to watershed delineations and the density of climate data. Five delineation scenarios combined with two levels of subbasins (8digit HUC and 12-digit HUC), two levels of HRU thresholds (1 km 2 /1 km 2 and 70%/70%) and two levels of climate data (PRISM and PRISM*) were simulated. The same set of default parameters generated within HAWQS was used for all five delineation scenario simulations, which included some "partial built-in calibration" for some UMRB 8-digit subbasins in HAWQS. The results of the analyses lead to the following conclusions.
1. The hydrologic response can be varied in different parts of the area for a large basin when changing delineation schemes, since the watershed delineation affects the landscape position, land-use distribution and soil distribution spatially. Results showed that groundwater flow and surface runoff are more sensitive to the areas with soils of higher infiltration rates, and the steeper areas generated a greater amount of lateral runoff than flat areas when applying denser delineation schemes. The statistics of monthly streamflow at St. Paul were satisfactory for all scenarios. However, several evaluation criteria of monthly streamflow declined considerably at the Clinton and Grafton gauges when using 12-digit subbasin scenarios.
2. The hydrologic cycle changed in the shift from 8-digit to 12-digit subbasin. The average length of channels and width of channels decreased in the 12-digit subbasins versus Water 2021, 13, 422 20 of 23 8-digit subbasins, while the average slope of channels and drainage density increased. For the whole study area, average slope of channels increased from 0.2% to 0.5%, and drainage density increased from 0.05 to 0.17. This resulted in higher lateral flow and groundwater flow estimates for the 12-digit subbasin scenarios relative to the 8-digit subbasin simulations. This was especially true for lateral flow, as reflected by an increase of nearly 100% from 8-digit to 12-digit subbasins. Predicted sediment yields also increased with a greater number of subbasins represented by the 12-digit scenarios. This is due to changes in channel length and slope, which affect the deposition and degradation of sediments.
3. A finer HRU delineation can result in increased streamflow because it captures a refined level of rainfall spatial variability. The New-HUC12 scenario generated higher averages of monthly streamflow than the HUC12 scenario at three gauge sites during the period from 1983 to 2005, because the New-HUC12 scenario had a finer delineation with 120,454 HRUs. Similarly, streamflow increased from the HUC12* scenario to the New-HUC12* scenario. Although the evaluation statistics of monthly flow are not always improved with a finer delineation, the magnitude of simulated streamflow increased and was closer to the observed flow in the finer watershed delineation scenarios for this UMRB study.
4. Denser climate data resulted in a SWAT-predicted increase in streamflow. Monthly streamflow increased on average by 8.32%, 6.92% and 6.30% at the St. Paul, Clinton and Grafton, respectively, from the PRISM* to the PRISM dataset. Spatial deviations in rainfall of no more than 60 mm and in temperature of <0.8 • C were observed between these two climate datasets. Spatial deviations in climate data input led to a 9 mm decrease in ET, 10 mm increase in surface runoff and 12 mm increase in water yield for the entire UMRB. There are distinct seasonal variations between monthly differences of surface runoff and water yield. The cumulative difference between PRISM and PRISM* during the growing season accounted for 75% of the total annual difference on average. The PRISM* dataset, which has denser gauges, tended to generate more runoff during the summer months as compared to PRISM dataset.
The New-HUC12* scenario had the finest watershed delineation and denser climate data, so it produced the highest streamflow at all three gauge sites. This indicated that the New-HUC12* scenario captured the most spatial information. In this study, it is not aimed at the better choice between 8-digit and 12-digit subbasins, or at the best HRU threshold, but at the hydrologic response to different delineations for a large basin like UMRB. The extension of this work involves the investigation of more local hydrology stations at a tributary to further assess the hydrologic response of the delineation schemes. Those local stations can stand for various smaller drainage areas contained within the overall large basin. It is also of interest to calibrate and validate 8-digit and 12-digit models separately and explore the portability for sensitive parameters. Above all, this study provides modelers insights to know how watershed delineation affects hydrologic cycle and to balance the level of model representation and the computational limitations.