Continuous Simulation of Highly Urbanized Watershed to Quantify Nutrients’ Loadings

: The United States has witnessed various extreme land use changes over the years. These changes led to alterations in watersheds’ characteristics, impacting their water quality and quantity. To quantify this impact in highly urbanized watersheds such as the Chicago Metropolitan Area, it is crucial to examine the characteristics and imperviousness distribution of urban land uses and available point and non-point sources. In this paper, the effect of urban runoff and nutrient loadings to water bodies in the Chicago River Watershed resulting from level (III) detailed urban land uses is investigated. A watershed scale hydrologic and water quality simulation using BASINS/HSPF model was developed for the highly urbanized watershed. Appropriate considerations were given to the effective impervious area (EIA). The results from the ﬁve-year calibrated water quality simulation were reasonably reﬂected with observed data in the study area and nutrient loadings of both point and non-point sources for 44 different land uses were found. The export coefﬁcients (EC) values obtained are site-speciﬁc depending on conditions and variables at the watershed level such as physical characteristics, land use management practices, hydro-meteorological and topographical data, while using a continuous simulation approach and watershed perspective analysis. This is the ﬁrst attempt to measure and model nutrients’ loadings using detailed land use types in the Chicago River Watershed. The proposed continuous calibrated and validated model can be used in the investigation and analysis of different scenarios and possible future conditions and land utilization.


Introduction
Land use changes in highly urbanized watersheds can adversely affect runoff quantity and transport patterns. Transforming land uses into impermeable surface-such as buildings, roads, parking lots, and pavements-result in the change of water quantity and quality available for direct runoff, stream flow, and groundwater flow [1][2][3][4], leading to alterations in watershed characteristics and increase of pollutants' quantities from both point and non-point sources [3,[5][6][7]. Several studies have shown how water quality of natural water bodies in urban watersheds have drastically changed over the years [1,4,[8][9][10]. To analyze a water body, it is essential to estimate nutrient loads and concentrations transported by a stream during a given period of time to identify source areas. These are particularly important when considering developing total maximum daily loads (TMDLs) as mandated by the Clean Water Act (CWA), land use management and best practices, and mitigation strategies. Using models to simulate a simplified form of reality has been known for a long time. Physical models are known integral modeling tools to investigate contaminants and contaminant concentrations. Along with monitoring data, physical models estimate concentrations and exposures providing information and evidence that can support risk analysis and management decisions in any specific site especially if drinking water resources are at risk [11,12]. Unsafe levels of contaminants in drinking water can cause health effects and chronic illnesses. The type of contaminant and its concentration in the water are factors that can greatly influence whether a contaminant will lead to adverse health outcomes and ecological effects. The impacts of small range urban land uses on both water quantity and quality were the focus of many watershed scale studies on aspects such as hydrology, climate, and ecology [13][14][15][16][17][18][19]. Moreover, several monitoring and modeling studies have consistently shown that urban pollutant loads increase with increased imperviousness degree of urban watersheds [20,21]. The significance of the effect of imperviousness on water quality degradation of aquatic resources and surface waters were reported as well [22]. Similarly, impacts on watersheds where spatial variation of urbanization was considered, showed high impact on runoff and nitrogen levels directly proportional to the level of urbanization investigated [23]. Understanding those impacts of land uses and how they affect water quality is essential for planners, regulatory authorities, and decision makers to develop comprehensive plans to tackle, address, and mitigate environmental issues [24,25]. However, quantifying those impacts in a watershed based on detailed urban land use is still considered developmental [1,3].
Hydrological models that are coupled with geographic information systems (GIS) and remote sensing proved to be powerful tools in conducting watershed studies [26][27][28][29]. Integrated approaches that involve the use of statistical and spatial analyses plus hydrologic modeling to examine the effects of land use on water quality were successfully utilized [3,30]. However, most studies depend on local field scale studies and focus on a small range of land use patterns to view the problem [31][32][33][34][35]. Watershed management practices are analyzed using event-based rainfall-runoff models [36][37][38][39]. These are temporal scale models that simulate single storm events not taking into account the hydrologic cycle process. Continuous models are used to investigate long-term processes, such as fate and transport of pollutants [40][41][42]. The continuous hydrological models consider the longterm effects of the hydrological cycle, hydrological changes, and watershed management practices.
In order to accurately estimate pollutant loads using hydrological models, accurate estimates of runoff volumes are essential. Since impervious area is a rough indication of the total watershed utilized, the effective impervious area (EIA) as a portion of the total impervious area (TIA) should be considered [43][44][45]. The EIA is the portion of the TIA within a watershed that is partially or totally connected to the drainage collection system. It is the most important and hardest parameter to determine for a watershed [43]; some ways to determine EIA are field measurements, empirical equations, and calibrated computer models [45][46][47][48]. Rooftops, parking lots, street surfaces, and paved driveways that are directly connected to the storm sewer system, are all included in the EIA [46]. The EIA for a given basin is usually less than the TIA in hydrologic analysis. However, in highly urbanized basins, EIA values can approach and equal TIA ones [44]. In this study BASINS and HSPF models were used to assess the water quality in the Chicago River Watershed. BASINS is a multi-purpose environmental analysis system that integrates geographical information system (GIS), national watershed data, and different modeling tools (such as HSPF, SWAT, SWMM, etc.) into one convenient package. The system is a widely accepted watershed-based water quality assessment tool that promotes better assessment, integration, and management of point and non-point sources [1,3,30,49,50]. HSPF is a comprehensive watershed scale model that performs continuous simulation of hydrology and water quality. It is extensively used to model urbanized watersheds [35,[51][52][53][54][55]. It performs continuous simulation of nonpoint source hydrology and water quality, combines point source contributions, and performs water quantity and quality routing in the watershed reaches [40]. HSPF lacks the capability to simulate storm sewer networks (Mohamoud et al., 2010). However, studies show that, among reviewed models that simulate storm water quantity and quality in urban environments, HSPF is the most comprehensive and flexible hydrology and water quality model available [29,56,57]. The HRUs resulted from the delineation process requires input data such as temperature, precipitation, potential evapotranspiration, and parameters related to land use, soil characteristics, and agricul-tural practices to simulate hydrology, sediments, and nutrients. For accurate hydrologic analysis, EIA parameters are very important to determine for the urban watershed [46]. Considering non-point sources for nutrients' load estimations can give unrealistic results because the land covers in urban areas are mostly impervious and drainages are routed to water reclamation plants WRPs (which may or may not be in the same basin), then are discharged into streams as point sources [16]. Calculating the EIA for all impervious areas in each HRU is a novel method to consider these limitations using HSPF in order to simulate and predict the impact of urban land use on nutrient loadings into watershed water bodies.
Export coefficients (EC) that result from the continuous modeling process, represent the concentration of a specific pollutant in stormwater runoff discharging from a specific land use type within a watershed [58]. They are required in several water quality models to calculate runoff pollutant loads into watersheds and are reported as mass of pollutant per unit area per year (e.g., lb/ac/yr) [58]. EC values are the combination of several site-specific variables and conditions at the watershed level including physical characteristics, land use management practices, hydro-meteorological, and topographical data [58,59]. Total nitrogen (TN) and total phosphorus (TP) are the most common pollutants for which export coefficients are usually generated [58]. Calculating site-specific export coefficients using locally collected data can be cost-prohibitive; thus, researchers or regulators often use values that are already available in the literature [58]. ECs are estimated in many studies but only for a limited range of land use types [58][59][60].
For the Chicago River Watershed, many studies were conducted, but generally as part of studies to investigate the flow and water quality for the Upper Illinois River Basin system, and not the individual highly urbanized watershed [61][62][63][64][65][66]. The limited land use categorization used could not explain the more detailed behavior of the highly urbanized watershed. In this paper, detailed urban land use effect on nutrients runoff to water bodies in the Chicago River Watershed was simulated. Five years simulation of water quality using the BASINS/ HSPF resulted in TN and TP export coefficients for level (III) land use. Pollutants simulated for both pervious and impervious land segments are total ammonium (NH 3 + NH 4 ) as N, total nitrate (NO 3 + NO 2 ) as N and ortho phosphorus (PO 4 ). The EC values presented in this paper along with the suggested method to calculate impervious areas in HSPF represent a novel way to quantify the pollutant loads. It is the first attempt at measuring and modeling nutrient using detailed land use, watershed perspective analysis, and a continuous simulation approach in the Chicago River Watershed. The proposed continuous calibrated and validated model provides a planning tool for regulatory environmental agencies and can be used in the investigation and analysis of different scenarios and possible future conditions in the watershed. The approach can be applied unmodified into any other watershed analysis studies.

Study Area
The Chicago River Watershed is the smallest part (6%) of the Upper Illinois River Basin (UIRB), which is the second-largest drainage basin in the world. The watershed area is located in northern Illinois and drains approximately 645 square miles. The watershed is confined within latitudes 41 • 11 and 42 • 20 N and longitudes 87 • 32 and 88 • 46 W. The North Branch Chicago River originates as three tributary streams-West Fork, Middle Fork, and the Skokie River-that flow and joins the South Branch of the river in downtown Chicago. The river flows westwards into the Chicago Sanitary and Ship Canal joins the Des Plaines River as a tributary of the Illinois River which flows southwest across the state and join the Mississippi River system. The uppermost bedrock of the basin is comprised of mainly undifferentiated Silurian Devonian dolomite and limestone, and Ordovician shale [67]. The Chicago River and the Des Plaines Basins are naturally divided by a subcontinental fault. The origin of the fault is either because volcanic activity or from meteoric impact [67]. The average elevation in the Watershed is 443 ft above sea level and the average basin slope is 0.001. Because of the nature of the cool, dry winters and warm, humid summers in the basin, its climate is classified as humid continental. Most of the large daily fluctuations of temperature and precipitation in the basin are results of the combinations of cool, dry and warm, moist air. The average annual temperature ranged from 46 • F to 51 • F in winter and 77 • F to 82 • F in summer. The average annual precipitation is 16 to 18 in and the average snowfall is 50 in. Evapotranspiration returns 70% of the average annual precipitation to the atmosphere [67]. The watershed has a moderately slow infiltration rate with very poorly drained areas along the western border of the watershed and the rest of the watershed is considered highly altered and mainly impervious.
The Chicago River Watershed is significant for its navigable systems, and particularly for the Chicago Sanitary and Ship Canal that connects Lake Michigan and the Mississippi River. The basin witnessed a steadily growth in population over the years due to the urban and industrial growth in the area. Major changes took place in the region such as the construction of navigable waterways, diversion of Lake Michigan water, and construction of wastewater-treatment plants. Land use, urbanization, and population change resulted in numerous inputs of contaminants and nutrients from manmade sources including municipal and industrial releases, urban runoff, and atmospheric deposition. The urban wastewater disposal, with storm runoffs had significantly affected the quantity and quality of surface waters in the watershed. Currently, the Chicago River Watershed is approximately 82% urban land use and considered highly urbanized area. Per the Census Bureau's urban-rural classification, urban areas represent densely developed territory, and encompass residential, commercial, and other non-residential urban land uses. The Chicago Metropolitan Agency for Planning (CMAP) land use Inventory for the years 2005, 2010, and 2015 divided Chicago River the urban land use in the study area as approximately 56% residential, 10% commercial, 10% industrial, 10% institutional, 15% Transportation/utilities, and 21% open space, agriculture, vegetation, wetland, and water. Detailed land use percentages compiled from CMAP will be shown later in Figure 6a. The land use inventory data was created using digital aerial photography and supplemented with data from numerous governments and private-sector sources [68].
Much of the pollutant load in the runoff originates from impervious surfaces, particularly roadways and parking lots. Some of the more common water quality impacts of stormwater runoff are sediment contamination, nutrient enrichment and toxicity to aquatic life, bacterial contamination, salt contamination, impaired aesthetic conditions, and elevated water temperatures [69]. Development also alters runoff patterns by changing the lay of the land, and thus, drainage patterns, resulting in a dramatic increase in the rate and volume of stormwater runoff and a reduction in groundwater recharging. In general, nutrient loads for nitrogen and phosphorus were greatest from the urban center of the Chicago Metropolitan Area, reflecting the effect of wastewater return flows to the Chicago River and the Chicago Sanitary and Ship Canal.

Data Sources and Types
For this study, as part of building a framework for modeling the Chicago River Watershed, a data warehouse was created to aggregate available data from various agencies in the Chicago River Watershed. This helps to integrate the data for watershed assessment and physical modeling. Water quantity, water quality and land use data were compiled from different sources. The sources include U.S. Geologic Survey (USGS), Metropolitan Water Reclamation District of Greater Chicago (MWRD), Chicago Metropolitan Agency for Planning (CMAP), US Army Corps of Engineers-Chicago District (USACE), BASINS data store, and data from several facilities, treatment plants, and combined sewer overflows (CSOs) sources that are permitted through National Pollutant Discharge Elimination System (NPDES) permits. Figure 1 shows some of the data sources used in the model and their locations. There are 18 active USGS stations where water quantity data were obtained. They measure daily flow and gage heights in the watershed. Water quantity data were obtained from the MWRDGC available 41 stations within the watershed that measures up to 65 different water quality parameters monthly or bi-monthly. Land use inventory for 2001 and 2005 acquired from CMAP were utilized. Meteorological data comprised precipitation, air temperature, potential evapotranspiration, wind speed, solar radiation, dew point temperatures, and cloud cover were compiled from Chicago O'Hare Airport metrological station. Description of data sources are shown in Table 1.

Watershed Simulation
In the following subsections, steps to develop the water quality model to quantify the effect of detailed (level III) land use on nutrient loading in the Chicago River Watershed using BASINS/HSPF are outlined. It explains types and sources of data used and how hydrologic and water quality models were constructed and used in the BASINS/HSPF model environment.

Watershed Simulation
In the following subsections, steps to develop the water quality model to quantify the effect of detailed (level III) land use on nutrient loading in the Chicago River Watershed using BASINS/HSPF are outlined. It explains types and sources of data used and how hydrologic and water quality models were constructed and used in the BASINS/HSPF model environment.
BASINS provides different data layers to HSPF including: National Land Cover Data (NLCD or GIRAS) land use data to calculate land use distribution within watershed, meteorological data to provide meteorological data requirements, digital elevation model (DEM) grid data to determine boundaries of watershed, reach files to determine stream networks, permit compliance system (PCS) data to provide loading information in the watershed, and monitoring data such as STORET and USGS data to provide water quality and quantity data. Detailed land use data for the watershed, including 44 types of urban land use, were acquired from CMAP and added to BASINS. BASINS land use available are either Level (I) or Level (II) land use type. To fulfill the objective of the study level (III) land use type, CMAP's 2005 land use inventory, in shapefile format, was added. The CMAP land shapefile was clipped into smaller shapefile to fit the watershed study area and because of limitations in processing detailed land use classifications in BASINS for larger area. To run HSPF, data was formatted to watershed data management (WDM) files that contain time series data required by HSPF. A user control input (UCI) file contains all the needed parameter values and control specifications to run the HSPF model. Calibration and validation analysis to evaluate the model are done using the GenScenario tool in the BASINS package. Delineation is a requirement process by HSPF and is either performed automatically using DEM grids or manually where existing streams and basins are manually selected and used to determine the watershed. The automatic delineation process used in this study ended up determining the streams, sub-basins, and outlets layers used in the simulation. The Upper Chicago River subbasin is divided into homogeneous land areas known as hydrologic response units (HRUs) that were used to define six reaches and seven subwatersheds. HRUs can be impervious or pervious areas, which, once determined, would be modeled independently. HSPF main simulation modules are PERLND, IMPLND, and RCHRES and they simulate pervious land segments, impervious land segments and free flow respectively [70]. The hydraulic characteristic of each reach was defined by parameters that represent volume-discharge relationships for each reach. Relationship between water level, surface area, volume and discharge are fixed. In order to simulate hydrology, sediments, and nutrients for each individual HRU, input data such as meteorological data and parameters related to land use, soil characteristics for each HRU are needed [70].

Impervious Area Assumptions
To determine the percentage of impervious area assigned for each HRU that represent the land use in the area, some estimation and consideration had to be made. The EIA as a percentage of TIA was determined for basins that are directly connected to the drainage systems [46]. TIA is determined using land use maps and aerial photography [71]. The majority of impervious surface studies correlate percentage of impervious surface to land use by using estimates of the proportion of imperviousness within each class. The TIA and EIA percentages and equations, based on the literature [45], used in this study are shown in Tables 2 and 3, respectively.

Flow Simulation
The first component to be simulated in HSPF is the flow. IWATER and PWATER modules are used for the flow simulation. IWATER module simulates retention, routing, and evaporation of water from impervious segments. PWATER module computes the water budget and predicts the total runoff from pervious segments. The hydraulic behavior in the streams is simulated by the HYDR module. A fixed relationship is assumed among water level, surface area, volume, and discharge in each reach. In-stream simulation is assumed to be mixed-system with unidirectional and longitudinal flow simulation. The hydraulic characteristics of reaches are defined by parameters such as nominal upper zone storage, nominal lower zone storage, soil moisture infiltration rate, percent vegetation cover of each land use type, and groundwater recession rate. The parameters' values were populated with BASINS default values and literature values and later adjusted during hydrologic calibration.

Water Quality Simulation
Water quality simulation in HSPF is done using IQUAL and PQUAL modules for impervious and pervious land segments respectively. Pollutants are simulated in the IQUAL and PQUAL modules using either direct washoff by overland flow where pollutant is simulated based on basic depletion and accumulation rate; or by washoff of detached sediments that are simulated as a function of sediment removal. The direct washoff by overland flow approach was adopted for all simulated constituents since the study area is largely impervious. Several biological, physical, and chemical processes within stream reaches are simulated in HSPF using the RCHRES module. The reaches are assumed to be completely mixed and the flow is unidirectional. Point sources added to the simulation are the two known NPDES in the watershed, which are North Side WRP and Calumet Water WRP. NPDES were included as direct inputs to the main reaches in the watershed. Pollutant species considered are the total nitrates as nitrogen (NO 2 + NO 3 ), total ammonia (NH 3 + NH 4 + ), and TP as phosphorus. For this study only the North shore WRP was considered as Calumet WRP was out of the scope of the study area.

Model Calibration and Validation
Models need to be evaluated to provide a quantitative estimate of performance and predictive ability of the model [72]. For HSPF evaluation, a weight of evidence approach that considers multiple graphical and statistical model comparisons, is commonly accepted and used to assess model performance [73]. The calibration-validation process is a hierarchal process where the parameters are first developed, then hydrology calibrationvalidation is done, and finally, water quality calibration-validation is performed. In this study Pearson's correlation (r) and determination (r 2 ) coefficients were used for evaluation. For model performance, r values range from −1 to 1 and for r 2 values range from 0 to 1, a value >0.5 are considered acceptable [73]. A major drawback for r 2 is its dispersion property, since systematically over-or under-predicted values still result in good r 2 values. Another model evaluation criterion is the Nash-Sutcliffe efficiency coefficient (NSE). NSE lies between 1 (perfect fit) and −∞. The main disadvantage of NSE is that the differences between the observed and simulated values are calculated as squared values resulting in larger values being overestimated and lower values being neglected [72]. Other statistical indices used to evaluate the model performance in this study are root mean square error (RMSE), normalized root mean square error (NRMSE), and mean absolute error (MAE). RMSE and MAE measure the aggregated differences between simulated values and observed values where values close to zero indicate better model performance. The last evaluation index used in the study is the percent mean error (PME). PME is a general guidance measure that has been specifically used to evaluate the HSPF model performance [73]. It provides percent mean errors or differences between simulated and observed values, so that HSPF users can determine the level of agreement or accuracy (very good, good, and fair) that might be expected using the model [73]. PME target values for different modeling process [73] are shown in Table 4. Using all the above measures fulfilled the weight of evidence approach adopted in this study to evaluate the model.

Flow Calibration and Validation
Meteorological data from Chicago O'Hare airport station was used to simulate the stream flow. To calibrate and measure flow sensitivity to impervious land segments, the EIA equations proposed by Alley, Laenen, and Sutherland [46][47][48] shown Table 3 were adopted. The percentages of pervious and impervious land segments were found accordingly. The calculated land uses were used in the flow simulation iterative procedure. Due to limited availability of data and input parameters for pervious and impervious land segments, BASINS's default input parameters were initially applied. Parameters used were lower storage nominal (LZSN), upper zone storage nominal (UZSN), mean soil infiltration rate (INFILT), lower zone evapotranspiration (LZETP), ground water detention storage (INTFW), and interflow recession coefficient (IRC). These parameters are used for calibration and can be estimated and adjusted until calibration process is complete. LZSN, UZSN and INFILT are parameters that affect the total annual flow volume. LZETP, INTFW, and IRC are parameters that affect the base flow conditions of the river and hydrograph shape and peak flow conditions. The various module parameters were repeatedly adjusted and simulated and observed values were compared until reasonable r and r 2 were obtained. Sensitivity analysis results for EIA equations for the study area are shown in Table 5. Lanen's EIA equation (equation 4) showed acceptable performance and it was the one adopted for calculating EIA and determining percentages of pervious and impervious land segments for the watershed.  Figure 2. The statistical indicators are fair based on graphical representation and according to the guidelines given by Donigian [73]. The model showed better performance in the validation period relative to the calibration period except for the poor r and fair to acceptable r 2 . According to Table 6, the model's performance is considered acceptable based on statistical indicators and acceptable ranges published in the literature for hydrologic simulation given the collective criteria measure.

Water Quality Calibration and Validation
The calibration and validation process in HSPF model is a hierarchical process that starts with the hydrology and ends with water quality constituents [73]. Constituents simulated were total ammonium (NH3 + NH4) as N, total nitrate (NO3 + NO2) as N, and ortho phosphorus (PO4). Total nitrogen and phosphorous loads were later calculated utilizing scripts provided in the model. Several nutrient parameters were added in the model for both pervious and impervious land segments including initial storage for each constituent, constituent washoff factor, and monthly constituent accumulation factor. These parameters were adjusted and calibrated until reasonable model behavior was seen.
The simulation results were compared to observed values at the North Branch of the Chicago River monitoring station at Grand Ave, Chicago. The location was chosen to represent the outlet for the sub-basin. Two factors limited the time period for the calibration and validation of the model used for the water quality simulation. The observed flow at the monitoring station was limited to the period 2002-2010 (with some missing data in the period 2003-2004), while the meteorological data was only available for the period 2002-2006. The initial simulation values were consistently over predicted, mostly during the wet season. To address this, calibration parameters-including the monthly accumulation factors and monthly values for limiting storage for both pervious and impervious land segments-were repeatedly adjusted. The instream simulation process parameters, nitrification and denitrification parameters (KNO320), along with oxidation rate (KTAM20) and algal growth rate parameters were all adjusted until reasonable modeling results were obtained. The calibration and validation results for total nitrates, total ammonium and ortho phosphate simulation are shown in Figures 3-5. Table 7 summarizes the water quality calibration and validation statistical for the simulated nutrients. Calibration results show that there is an acceptable agreement between the observed and simulated data.

Water Quality Calibration and Validation
The calibration and validation process in HSPF model is a hierarchical process that starts with the hydrology and ends with water quality constituents [73]. Constituents simulated were total ammonium (NH 3 + NH 4 ) as N, total nitrate (NO 3 + NO 2 ) as N, and ortho phosphorus (PO 4 ). Total nitrogen and phosphorous loads were later calculated utilizing scripts provided in the model. Several nutrient parameters were added in the model for both pervious and impervious land segments including initial storage for each constituent, constituent washoff factor, and monthly constituent accumulation factor. These parameters were adjusted and calibrated until reasonable model behavior was seen.
The simulation results were compared to observed values at the North Branch of the Chicago River monitoring station at Grand Ave, Chicago. The location was chosen to represent the outlet for the sub-basin. Two factors limited the time period for the calibration and validation of the model used for the water quality simulation. The observed flow at the monitoring station was limited to the period 2002-2010 (with some missing data in the period 2003-2004), while the meteorological data was only available for the period 2002-2006. The initial simulation values were consistently over predicted, mostly during the wet season. To address this, calibration parameters-including the monthly accumulation factors and monthly values for limiting storage for both pervious and impervious land segments-were repeatedly adjusted. The instream simulation process parameters, nitrification and denitrification parameters (KNO320), along with oxidation rate (KTAM20) and algal growth rate parameters were all adjusted until reasonable modeling results were obtained. The calibration and validation results for total nitrates, total ammonium and ortho phosphate simulation are shown in   Table 7 summarizes the water quality calibration and validation statistical for the simulated nutrients. Calibration results show that there is an acceptable agreement between the observed and simulated data. Statistical results and model performance are considered good to very good for all the constituents as it fell within the accepted tolerances suggested by Donigian (Table 4). Statistical results and model performance are considered good to very good for all the constituents as it fell within the accepted tolerances suggested by Donigian (Table 4).     Statistical results and model performance are considered good to very good for all the constituents as it fell within the accepted tolerances suggested by Donigian (Table 4).     Statistical results and model performance are considered good to very good for all the constituents as it fell within the accepted tolerances suggested by Donigian (Table 4).

Total Annual Loads of Nutrients
Different studies have done significant progress, using long term modeling and machine learning, towards estimating inputs of total nitrogen and total phosphorus loads to receiving waters in metropolitan areas [74,75]. However, there have been fewer studies that explore the long-term hydrological conditions and detailed land use in the context of high urbanized watershed. For this study, HSPF modules PQUAL and IQUAL were used to estimate annual loadings of total nitrogen and total phosphorus from forty-four different land use types in the Upper Chicago River Basin. Based on the results from the calibrated and validated water quality model, the total annual loads from the Upper North Chicago River sub-basin were computed using scripts provided by HSPF that transfer simulated total nitrates and total ammonium into total nitrogen loads and simulated ortho phosphate into total phosphorus loads. To compare different land use areas, total nitrogen and total phosphorous associated with each land use type percentage values are shown in Figure 6. The simulation results show that the highest total nitrogen and total phosphorus loads produced during the period of 2000 to 2005 in the Upper Chicago River sub-basin was for residential single-family land use. This is an expected result since residential single-family land use is the dominant land use type in the basin, as shown in Figure 6a. It was difficult to determine the accuracy of the loads simulated by this since no information available that can relate the contribution of a detailed land use, level (III), to the total nitrogen and total phosphorus loads to the Chicago River Watershed or any similar highly urbanized watersheds. However, based on the results of water quality model calibration and validation, it can be assumed that this a unique work in estimating total nitrogen and total phosphorus loads from a detailed land use segments in the watershed.

Export Coefficients of Detailed Land Uses
It is essential to adopt approaches to estimate different region-specific nutrient export coefficients other than the ones available in literature [74]. Although the literature is showing limited ranges of EC values for general land uses, the ranges of EC values presented in this section were found comparable with ranges established in literature [45,58,76]. The export coefficients presented in this study are believed to be the first attempt to measure and model TN and TP using detailed land use types in the Chicago River Watershed or any similar highly urbanized watersheds using a continuous simulation approach and watershed perspective analysis. ECs are very useful indicators that allow the prediction of the possible yield of nutrients reaching receiving water bodies. The average obtained export coefficients for total nitrogen and total phosphorous, respectively, for each land use type investigated in this study are shown in Figure 7.

Export Coefficients of Detailed Land Uses
It is essential to adopt approaches to estimate different region-specific nutrient export coefficients other than the ones available in literature [74]. Although the literature is showing limited ranges of EC values for general land uses, the ranges of EC values presented in this section were found comparable with ranges established in literature [45,58,76]. The export coefficients presented in this study are believed to be the first attempt to measure and model TN and TP using detailed land use types in the Chicago River Watershed or any similar highly urbanized watersheds using a continuous simulation approach and watershed perspective analysis. ECs are very useful indicators that allow the prediction of the possible yield of nutrients reaching receiving water bodies. The average obtained export coefficients for total nitrogen and total phosphorous, respectively, for each land use type investigated in this study are shown in Figure 7.

Conclusions
Simulation models at the watershed scale offer an effective watershed management tools to estimate nutrients yields for a wide spectrum of problems dealing with surface waters. Since the Chicago River Watershed is 82% urban land use (i.e., highly urbanized area), examining effect of land use on water quality requires a detailed level of land use.

Conclusions
Simulation models at the watershed scale offer an effective watershed management tools to estimate nutrients yields for a wide spectrum of problems dealing with surface waters. Since the Chicago River Watershed is 82% urban land use (i.e., highly urbanized area), examining effect of land use on water quality requires a detailed level of land use. The five years of water quality simulation done using BASINS, coupled with the continuous simulation watershed scale model HSPF, resulted in unique export coefficients for level (III) detailed land use for the Chicago River Watershed. This can be considered as a successful contribution to measure and model nutrients in the watershed and allow the prediction of the possible yield of nutrients reaching receiving water bodies. This can provide a planning tool for regulatory environmental agencies in the Chicago River Watershed to develop better management programs especially if total maximum daily loads (TMDL) and water quality standards (WQS) are developed for the Chicago River Watershed. The model can be utilized in the investigation and analysis of future scenarios in the watershed. The range of land use export coefficients obtained from long term continuous simulation reflects the different conditions of watershed and different meteorological and physical variables included in the simulation and hence provide a perfect input for a multi-objective optimization approach to evaluate multiple scenarios that seek to find optimal land use change and distribution in highly urbanized developed watershed. Based on different detailed land use types, scenarios that take into account different combinations of pervious and impervious land use segments and tradeoffs between them (e.g., changing an impervious parking lot land use into pervious, etc.), along with factors such as environmental, social, and economical factors can be investigated as part of planning and decision-making tool. The multi-objective optimization approach will allow the optimizing of independent objectives to find the best land use combination while the high priority goal is to meet certain water quality standards regarding nutrient loadings of total nitrogen (TN) and total phosphorus (TP). The model can be successfully applied to any highly urbanized watershed with appropriate consideration given to EIA.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.