Review of Watershed-Scale Water Quality and Nonpoint Source Pollution Models

Watershed-scale nonpoint source (NPS) pollution models have become important tools to understand, evaluate, and predict the negative impacts of NPS pollution on water quality. Today, there are many NPS models available for users. However, different types of models possess different form and structure as well as complexity of computation. It is difficult for users to select an appropriate model for a specific application without a clear understanding of the limitations or strengths for each model or tool. This review evaluates 14 more commonly used watershed-scale NPS pollution models to explain how and when the application of these different models are appropriate for a given effort. The models that are assessed have a wide range of capacities that include simple models used as rapid screening tools (e.g., Long-Term Hydrologic Impact Assessment (L-THIA) and Nonpoint Source Pollution and Erosion Comparison Tool (N-SPECT/OpenNSPECT)), medium-complexity models that require detail data input and limited calibration (e.g., Generalized Watershed Loading Function (GWLF), Loading Simulation Program C (LSPC), Source Loading and Management Model (SLAMM), and Watershed Analysis Risk Management Frame (WARMF)), complex models that provide sophisticated simulation for NPS pollution processes with intensive data and rigorous calibration (e.g., Agricultural Nonpoint Source pollution model (AGNPS/AnnAGNPS), Soil and Water Assessment Tool (SWAT), Stormwater Management Model (SWMM), and Hydrologic Simulation Program Fortran (HSPF)), and modeling systems that integrate various sub-models and tools, and contain the highest complexity to solve all phases of hydrologic, hydraulic, and chemical dynamic processes (e.g., Automated Geospatial Watershed Assessment Tool (AGWA), Better Assessment Science Integrating Point and Nonpoint Sources (BASINS) and Watershed Modeling System (WMS)). This assessment includes model intended use, components or capabilities, suitable land-use type, input parameter type, spatial and temporal scale, simulated pollutants, strengths and limitations, and software availability. Understanding the strengths and weaknesses of each watershed-scale NPS model will lead to better model selection for suitability and help to avoid misinterpretation or misapplication in practice. The article further explains the crucial criteria for model selection, including spatial and temporal considerations, calibration and validation, uncertainty analysis, and future research direction of NPS pollution models. The goal of this work is to provide accurate and concise insight for watershed managers and planners to select the best-suited model to reduce the harm of NPS pollution to watershed ecosystems.


Introduction
Nonpoint source (NPS) pollution is a pervasive source of water pollution in the world. In practice, a watershed is a fundamental unit of monitoring and management of NPS pollution. To better understand the complex mechanisms and processes of NPS pollution, various watershed-scale NPS models and tools were developed to understand NPS pollution, and to evaluate water quality. These NPS models and tools are widely used to identify critical source areas of pollutants [1][2][3], evaluate the effects of NPS pollution on the water environment [4,5], future scenarios evaluation of hydrology and water quality [6,7], assist in the planning and implementation of best management practices (BMPs) [8][9][10], support development of water quality criteria or standards [11,12], and provide informed decision support for policy-makers [13][14][15]. Models have become essential tools in the effort to reduce the adverse effects of NPS pollution. However, the current diversity of the NPS pollution models makes it difficult to select the most suitable model for a given NPS pollution issue.
The diversity of current models stems from a history of watershed management and planning demands. With the development of model theory, computer technology, and environmental legislation, numerous water quality and NPS pollution models were developed or integrated into watershed models (i.e., hydrology models or rainfall-runoff models) [16]. These models use the watershed as the fundamental spatial unit to simulate various hydrologic, hydraulic, soil erosion, sediment transport, and nutrient dispersion processes that accounts for surface water, groundwater, and their interaction as a whole system [17]. The Stanford Watershed Model, developed in 1959Model, developed in -1966, was the first computer model to conduct watershed hydrology analysis and modeling that subsequently evolved to the well-known Hydrologic Simulation Program Fortran (HSPF) in the 1970s [18]. Most of the current NPS pollution models such as HSPF [19], Soil and Water Assessment Tool (SWAT) [20], Agricultural Nonpoint Source pollution model (AGNPS) [21], Long-Term Hydrologic Impact Assessment (L-THIA) [22], Generalized Watershed Loading Function (GWLF) [23], Source Loading and Management Model (SLAMM) [24], Stormwater Management Model (SWMM) [25], and Better Assessment Science Integrating Point and Nonpoint Sources uncertainty tools, as well as that the future challenges in the field of water quality modeling particularly related to large data and accurate predictions.
In spite of the aforementioned review articles, there is limited guidance on how to select an appropriate model for the purpose of watershed management and planning. In this article, 14 commonly used watershed-scale NPS pollution models that predict flow, sediment and nutrient concentrations or loads, and estimate watershed water quality were evaluated. These models include: AGNPS, BASINS, GWLF, HSPF, L-THIA, SLAMM, SWAT, SWMM, Automated Geospatial Watershed Assessment Tool (AGWA) [39], Loading Simulation Program C (LSPC) [40], Nonpoint Source Pollution and Erosion Comparison Tool (N-SPECT) [41], Watershed Analysis Risk Management Frame (WARMF) [42], Watershed Assessment Model (WAM) [43], and Watershed Modeling System (WMS) [44]. These models were evaluated and compiled with their key attributes including the primary intended use, model components, suitable land-use type, input parameter types, spatial and temporal scales, pollutants, strength and weakness, software developers or access, and availability. The above models were categorized into four groups: simple models, medium complexity models, complex nodels, and modeling system ( Figure 1). Details about each group are explained in Section 2. Most of these discussed NPS pollution models are no charge for public use. As this review was limited to a subset of well-known watershed-scale NPS pollution models that are in use, many other useful field-scale NPS pollution models and receiving water bodies or in-stream water quality models such as QUAL series [45] and Water Quality Analysis Simulation Program (WASP) [46], and other metamodels such as Smart prediction of the Impact of Management strategies (nutrient and water) on Phosphorus losses by LEaching (SIMPLE) [47] and others are not included in this review. In this work we develop a simple classification structure of model complexity and provided common examples of watershed-based water quality models in these complexity categories to provide insight for users and practitioners.
We attempted to not only reflect the main features of each watershed model, but also discussed current challenges in model selection such as spatial and temporal scale, calibration and validation, and uncertainty analysis. We also discussed the future potential direction of NPS pollution model research. This review will help modelers understand how these tools should be applied in practice and be useful for practitioners of watershed management and planning to make an informed decision while choosing an appropriate model for their application related to water quality and NPS pollution.

Watershed-Scale Nonpoint Source Pollution Model Evaluation
Watershed-scale NPS models can be classified according to a variety of criteria that include the methods used to quantify hydrologic processes (e.g., empirical, conceptual, or physically-based), spatial variability of input parameters (e.g., lumped, semi-distributed, or distributed), spatial and temporal scales (e.g., small, field, or watershed; event-driven or long-term simulation), or the uncertainty of model output (e.g., deterministic or stochastic) [48]. To simplify model selection, we categorized models into four groups: simple models, medium complexity models, complex models, and modeling systems. This classification considers the intrinsic structure and form of the tools and models as well as the data and Yuan et al. Page 4 calibration requirements necessary to apply the model. This provides a framework that also helps the end user understand the appropriate application and requirements for implementation of the model tool. Table 1 summarized the main characteristics of different NPS pollution model types.

Simple Models
Simple models have minimal data requirements that includes land-use, soil, precipitation, and digital elevation model (DEM) (optional for some specific areas). This type of model is built based on an empirical or statistical relationship between pollution loads or concentrations, land use, rainfall, and runoff volume. These models typically use the Soil Conservation Service-Curve Number (SCS-CN) method [49] to calculate overland runoff and employed the export coefficient or the Event Mean Concentration (EMC) method [50] to calculate pollution loading. Simple models are often used as quick screening tools to obtain the gross pollution loads at the outlet of a watershed or to evaluate long-term areal average pollution loads for a moderate or large size watershed [51]. Although the simple models described here can provide good generalized information with little calibration, some validation of models is preferred. Simple models do not consider the spatial route of flow, sediment and pollution transport, nor are pollutant fate and transport mechanism in water bodies included. Thus, these models possess limited capacity to simulate detailed hydrological and physicochemical processes. These tools alone may not be sufficient to assist in water pollution control decisions that require insight on opportunities or locations to implement prevention or regulatory measures. These tools are also limited to a small number of simulated pollutants. Simple models can be powerful tools for generalized understanding of the pollutant loads or concentrations but neglect details in physical processes of NPS pollution.
2.1.1. L-THIA-L-THIA (Long-Term Hydrologic Impact Assessment), developed by Purdue University, is designed to estimate the long-term effects of land-use change in urban or suburban areas on surface runoff, groundwater recharge, and NPS pollution [22]. L-THIA includes three versions: L-THIA basic model, a spreadsheet version of L-THIA runs on the internet [52]; ArcL-THIA, an ArcGIS-based extension tool [53]; L-THIA low impact development (LID) estimates the impacts of land-use changes and LID practices on urban runoff and water quality [10].
L-THIA was designed to quickly evaluate average annual runoff volumes and quantify NPS pollutant loads drained into water bodies in urban and suburban areas. L-THIA does not require intensive data input and can provide 'what-if' alternative future scenario analyses.
The L-THIA model consists of hydrology and water quality components. L-THIA uses the runoff curve number (CN) method to determine direct runoff and employs the EMC method to calculate annual pollutant loads. Lim,et al. [22] improved the model by using a single rainfall event to estimate average yearly NPS pollution for 15 pollutants. Contaminants incorporated by L-THIA include nitrogen, phosphorus, suspended sediment loads, bacteria, and metals [54].
L-THIA was successfully applied to evaluate NPS pollution in different places such as United States [55], South Korea [22], and China [56]. Zhang,et al. [56] used L-THIA to evaluate the spatial distribution of NPS pollution in Qingdao (10,654 km 2 ) China, and verified that L-THIA is a reliable model to provide informed decision for land use management and planning. You, et al. [57] simulated and validated nitrogen (N) and phosphorus (P) loads across different land-use using L-THIA model in the Sihu basin (11,547.5km 2 ) in China and found that the model performed well for estimating the average loads of N and P. Liu,et al. [58] used L-THIA-LID as an optimization tool to estimate the impact of land use and climate change on hydrology and water quality with future scenario analysis in Trail Creek Watershed (153.2 km 2 ) in northwest Indiana, United States.
Overall, L-THIA is a quick screening tool for NPS pollution and water quality evaluation in urbanized areas. Users with limited knowledge of hydrology and the environment can also utilize the L-THIA model. The data needed for model running are easily attainable, especially in the United States. The application of L-THIA needs no or only limited calibration when it is applied to watersheds across the Midwest U.S. [59]. However, L-THIA can only reflect the average annual runoff volume and pollutant loads, and ignores the spatial route of flow, sediment and pollutants.
The assumption of this model does not include pollutant fate and transport in receiving water bodies. Since the EMC method itself considers pollution concentration as a constant over time, L-THIA cannot reflect a change of management practice (e.g., fertilizer application) on land, nor express the dynamic relationship between seasonal or varied flow and concentration [36]. Validation of the L-THIA model remains a significant challenge without site specific data where the model is applied, especially in an application of a large watersheds that may have heterogeneity in land use, land cover, precipitation, soil, or locations that are not preloaded with topographic data.

N-SPECT/OpenNSPECT-N-SPECT (Nonpoint Source Pollution and Erosion
Comparison Tool), developed by National Oceanic and Atmospheric Administration (NOAA), allows managers to quickly examine relationships between potential water quality of water bodies, NPS pollution, and soil erosion in nearshore areas [41]. N-SPECT works as a plug-in extension of ArcGIS, and the latest version of N-SPECT is compatible with ArcGIS 9.3. However, users must have a license of ArcGIS and its spatial analysis tool before running this model. To expand access to users without an ArcGIS commercial license, the NOAA Coastal Services Center developed OpenNSPECT in 2011, which is a free and open-source version of N-SPECT. OpenNSPECT is a plug-in tool of open-source MapWindow GIS and has the same theoretical foundation as N-SPECT [60]. The software, user's manual and documentation of OpenNSPECT are available at [61].
N-SPECT/OpenNSPECT estimates runoff volume, sediment yield, and pollutant loads/ concentrations, identifies soil erosion susceptible areas, and estimates the impact of land-use changes on water quality [41]. The model can operate at an annual or event time step in a medium-to-large near-shore watershed, and support 'what-if' scenario analysis under different land use management practices. In N-SPECT/OpenNSPECT, the runoff CN method is used to estimate surface runoff; the EMC method is applied to calculate pollutant concentration; the Revised Universal Soil Loss Equation (RUSLE) and Modified Universal Soil Loss Equation (MUSLE) are employed to estimate erosion rate and sediment loads; the Sediment Delivery Ratio (SDR) method is used to evaluate sediment yield. Finally, the model generates an overall water quality rating and allocates this rating to river networks by comparing the simulated total pollutant concentrations to local water quality standards [41]. Evaluated pollutants include total nitrogen (TP), total phosphorus (TP), total suspended solids (TSS), zinc, and lead. Input data sets include DEM, land use, soil, slope, precipitation, rainfall and runoff erosivity (R-factor), local pollutant coefficients, and water quality standards. The model outputs have runoff volume and depth, accumulated pollutant loads and concentrations, soil erosion, and total sediment yield [41,62].
N-SPECT/OpenNSPECT was applied to understand and evaluate the effects of land use management practices on water quality, especially for nearshore ecosystem health. Maina,et al. [63] applied the N-SPECT to estimate annual sediment load per unit area in two catchments of Madagascar island in Australia to examine the relationship between coral reefs and environmental change in coastal watersheds. Butler,et al. [64] used N-SPECT model to calculate runoff volume and changes in the annual delivery of total N for each scenario analysis in Tully Murray catchment in Australia. Tulloch,et al. [65] employed OpenNSPECT model to simulate runoff and sediment discharge from New Ireland watershed (7404 km 2 ) while evaluating the impacts of the oil palm industry in the nearshore ecosystem health in Papua New Guinea.
N-SPECT/OpenNSPECT is a light-level screening tool for estimating water quality, soil erosion, sediment yield, and NPS pollution at a nearshore watershed. The model does not require intensive data input and is a simple model structure based on GIS raster calculation. Although the model manual claims that the model can be applied to any watershed, recent applications demonstrated that the model is often applied to ecosystem health assessment in coastal areas [63,66]. N-SPECT/OpenNSPECT can simulate the distinct contributions of sediment and pollutants from upstream areas. However, the model is not a sophisticated tool for watershed water quality and NPS pollution assessment. The model does not account for stormwater drainage, stream diversions, snowmelt, sediment redeposition and the dynamic processes of runoff, sediment, and pollutant transport [67].

Medium Complexity Models
Unlike the simple models discussed above, medium complexity models generally require relatively detailed data inputs such as topography, land use, soil, weather, and water quality. The type of model is typically used as a compromise between simple models and complex models. Medium complexity models not only account for the fundamental water and sediment movement processes, but also combine the empirical relationships to express nutrient loads. Meanwhile, it avoids applying a complicated physically-based watershed model that requires intensive data. The calculation of surface runoff is based on a water balance principle [68]. The temporal step of some models estimates a daily value of runoff, sediment, and nutrient loads. Although medium complexity models do not require rigorous calibration, model validation is also preferable. Most of these models can assist in the development of water quality criteria and are easy to operate compared to complex models.
However, these models have inherently restricted simulation capability for in-stream fate and transport of pollutants, and the number of simulated pollutants is limited as well.

GWLF-GWLF (Generalized Watershed Loading Function), initially developed by
Haith and Shoemaker [23], is used to estimate monthly flow, sediment, and nutrient loading, and assists TMDLs development from a medium-size urban or agricultural watershed with various land-use combinations [69]. The latest version of GWLF is integrated into the Mapshed model that currently maintained by Pennsylvania State University with online access at [70].
GWLF combines distributed/lumped parameters and estimates long-term continuous flow, sediment, and nutrient loads from the land surface into streams based on daily weather data input [1]. GWLF can simulate dissolved and solid-phase nutrient loads in streamflow and considers different pollution sources such as surface runoff, groundwater, and septic systems. In agricultural land, GWLF uses the CN method to calculate runoff and obtain dissolved nutrient loads by multiplying runoff volume by dissolved nutrient concentration from each land use type. GWLF computes solid-phase nutrient loads by multiplying monthly sediment yields by average sediment nutrient concentrations. Soil erosion can be obtained by applying a modification of the Universal Soil and Erosion Equation (USLE). Sediment yield is generated from soil erosion, runoff transport capacity, and sediment delivery ratio (SDR) [69]. In urban land, GWLF calculates runoff by SWMM [71] and STORM model [72]. Pollutant loads are regarded as entirely solid-phase, and calculated by exponential accumulation and washoff functions [69]. In groundwater, pollutant loads are spatially-lumped for a watershed and calculated as the product of subsurface flow and a watershed average nutrient concentration [73]. The inputs of GWLF include precipitation, temperature, runoff sources and transport, and chemical parameters on a daily time step. The model outputs include monthly flow, monthly soil erosion and sediment yield, monthly TN and TP loads in flow, monthly groundwater discharge to flow, annual erosion from each land-use type, and yearly TN and TP loads from each land use type [73].
GWLF was applied to assess the environmental impact of soil erosion and NPS pollution [27,74], estimate streamflow, sediment and nutrient loads [75], as well as simulate daily flow [68]. Niraula, et al. [1] applied GWLF and SWAT to identify the critical source areas of NPS pollution in the Saugahatchee Creek watershed (570 km 2 ) in east-central Alabama. Their findings showed that both models can accurately evaluate streamflow, but SWAT had a better performance in terms of predicting sediment yield, TN, and TP. Qi,et al. [76] also compared the performances of GWLF and SWAT in simulating monthly streamflow, sediment yield, and total N loads in the Tunxin catchment (2674 km 2 ) and the Hanjiaying basin (6736 km 2 ) of China. Their results support the fact that GWLF is easy to use since it requires fewer data to conduct simulations compared to the data needed for SWAT. Similarly, Gene [11] also demonstrated that GWLF is easy to use and is less complicated than SWAT and HSPF.
GWLF can be applied to an ungauged watershed, and it needs no calibration or only minimal calibration. If calibration is not conducted, all calibration parameters should be estimated with a similar method in the application of GWLF to a watershed [11]. However, GWLF only first identifies nutrient loads from different land-use types, then applies these results to the entire basin. The model needs distributed parameters as input for surface pollutant simulation but does not account for a spatial structure or channel routing of the flow component. In groundwater modeling, GWLF uses lumped parameters with a linear reservoir model that ignores the spatial variability of physical and chemical transport processes [76].

LSPC-LSPC
(Loading Simulation Program C ++ ), developed by Tetra Tech Inc. with funding from U.S. EPA, is a C ++ version of the HSPF model that can simulate hydrology and water quality [77]. The current version of LSPC is 5.0 that was released in 2015. Users can download LSPC 5.0 installable files and its user manual from the link [78].
Since LSPC rewrites the code of selected HSPF components in C ++ , the model has the same core algorithm of HSPF [40]. The main new features of LSPC include a modular structure, an input file organization form, model segmentation, meteorological linkage, data input/ output, routines and other capabilities that are not included in HSPC. LSPC can simulate flow, soil erosion and sediment transport, general water quality from both upland contributing areas and receiving water bodies, as well as modules for stream transport, total maximum daily loads (TMDL) calculation, and source allocations in an urban or agricultural watershed [79]. The model is driven by hourly precipitation and can model hourly surface runoff and subsurface flows. The spatial scale of LSPC application is applicable to a small or large size watershed composed of over 1000 sub-watersheds. A Microsoft Access database is used to manage data and weather files in ASCII format. The primary components of LSPC include hydrology, snow, temperature, irrigation, sediment, as well as water quality sub models like the general quality (GQUAL) and the reach quality (RQUAL) [77]. The simulated pollutants include sediment, nutrients, metals, dissolved oxygen (DO), biochemical oxygen demand (BOD), plankton, and other contaminants from pervious and impervious lands. The input of LSPC includes DEM, soil, land use, and meteorological data. The model produces a time series of nutrient loading, hydrographs, and impacts of predetermined Best Management Practices (BMPs) [79].
LSPC is an efficient and flexible watershed hydrology and water quality model [80]. Shen and Zhao [81] employed LSPC to model surface runoff and subsurface flow while estimating bacteria nonpoint source loading in Sandy Bottom Branch watershed (6.9 km 2 ) in Virginia. They found that LSPC can reasonably simulate streamflow over a 10-year period. Huang and Xiang [82] applied LSPC to investigate point source and NPS pollution of the Panjiakou Reservoir basin (42,443 km 2 ) in north China. They indicated the developed model based on LSPC had sufficient accuracy in representing the hydrological characteristics of the watershed.
LSPC was designed to facilitate data management, organization, and modeling for a large complex watershed such as a municipal boundary (e.g., state or province) scale. LSPC has no inherent limits on the size and spatial-temporal resolution of input data, and the model uses a Microsoft Access database to manage model configuration and parameter files [79]. LSPC incorporates point source and NPS pollution, and also can be applied to water quality criteria like TMDL development. The output from LSPC can also be easily linked with other in-stream water quality models such as Environmental Fluid Dynamics Code (EFDC), WASP, CE-QUAL-W2, and System for Urban Stormwater Treatment and Analysis Integration (SUSTAIN). However, LSPC does not allow multiple sub-basins to connect to a single reach, nor deal with complex groundwater routing to simulate interactive process between surface and subsurface water [79].

SLAMM-SLAMM (Source Loading and Management Model), initially developed by Robert Pitt and John
Voorhees, is a pollutant source area identification and management planning model for an urban stormwater runoff [24]. SLAMM/WinSLAMM is an event-based continuous urban stormwater quality model and planning tool that can predict runoff discharges and pollutant loads for each source area within each land use type [84]. The model exams the relationships between source areas of urban pollutants and stormwater runoff quality, including evaluating contributions of source areas, estimating stormwater outfall discharge, particulate washoff, stormwater controls practices, and predicting pollutant concentration and loads [85]. The types of urban land-use in SLAMM involve residential, institutional, commercial, industrial, open space, and freeways and corresponding 14 source areas for each land use. The SLAMM calculates pollutant concentration and loads by using the discharge volume and suspended solids concentrations at the outfall. The simulated pollutants include particulate solids, P, TKN, COD, metals, and fecal coliform bacteria. The input data include storm, pollutant probability distribution, runoff coefficient, particulate solids concentration, street delivery on different land-use types, particle size distribution on each source area and flow peak of average flow ratio. The output includes runoff and flow summary, outfall and source area totals, source areas by land use, and outfall for each rain [86]. The improvement of WinSLAMM includes tracking pollution particle size distribution through stream network and report pollutant reduction from each land-use type and control device.
SLAMM/WinSLAMM has shown its reliability in predicting the impacts of stormwater control practices on flows and pollutant loads. Hurley and Forman [14] used WinSLAMM to evaluate the potential reductions of phosphorus loading to the Charles River in the Allston Campus Institutional Site (0.75 km 2 ) and Zakim Industrial Area (0.81 km 2 ) in Boston. Selbig, et al. [87] employed WinSLAMM to analyze the impact and the spatial distribution of particles in stormwater on the required size of stormwater control measures intended to meet pollutant reduction target.
SLAMM/WinSLAMM primarily depends on field observations rather than pure theoretical estimates that have not been widely verified in practice. The model was built based on the theory of small storm hydrology [88] with the concept that stormwater contamination is caused by frequent, small or moderate rain events [85]. The model considers different stormwater controls for a long series of rains. The model highlights water quality simulation, rather than only treating it as a part of hydrology or physcial hydraulics [14]. However, the current versions of WinSLAMM do not consider snowmelt, baseflow conditions, or account for in-stream processes that can remove or transform pollutants. Also, the model cannot simulate mass erosion from urban construction sites. Additionally, WinSLAMM is not a public domain software and its help documentation is simple.

WARMF-WARMF (Watershed Analysis Risk Management Framework), mainly developed by System
Water Resources, Inc., is a comprehensive decision support system for watershed management and analysis, and can support water quality criteria development [89]. However, WARMF is not public domain software. Users can obtain the software from [89].
WARMF can predict short or long-term physical, biochemical processes in any river basin, and includes vertical and lateral flow, groundwater, sediment loads, the fate and transport of nitrogen, phosphorus, metals and pesticides, and supports watershed management criteria development such as TMDL calculations, point/nonpoint source pollutants analysis, and watershed water quality planning [90,91]. WARMF can continuously simulate hydrologic processes with a daily time step. WARMF is composed of five interlinked modules, that include the engineering module, consensus module, data module, manager module, knowledge module, and TMDL module. The engineering module is used to estimate hydrology and water quality. The consensus module is employed to evaluate management practices and conduct scenario analysis. The data module is used to edit database files, and represents outcome with graphs and spreadsheets. The manager module is designed for realtime watershed management. The knowledge module collects various laws/regulations and case studies regarding the watershed. The TMDL module instructs the user to calculate a single pollutant from a point source or nonpoint source for a watershed to meet designated criteria [42]. The simulated pollutants include inputs of acid mine drainage, inputs from septic systems, bacteria, DO, mercury loading and transport, sediment, periphyton in rivers, and algae in stratified reservoirs. The input of WARMF includes DEM, land use, soil, meteorology, air quality, point source discharge, monitored streamflow, and water quality. The final products have a TMDL implementation plan or watershed management plan to support control of point and NPS pollution.
Geza, et al. [92] used WARMF to a mountain watershed (126 km 2 ) and evaluates predictive uncertainty by using a calibration and uncertainty analysis algorithm. Dayyani,et al. [91] developed DRAIN-WARMF to improve the deficiency of WARMF, which was used to simulate subsurface flow and the nitrogen fate and transport of a small agricultural watershed (24.3 km 2 ) in eastern Canada.
WARMF is a based on the physical movement of water, sediments, and nutrients in a watershed. The tool is suitable for locations like acid mine drainage, mercury pollution, and on-site wastewater systems. WARMF is stand-alone software that possesses a GIS interface, so users do not need to procure additional GIS licenses to drive the model. However, WARMF does not account for a tile drainage system, deep groundwater aquifers, and groundwater quality [93].

Complex Models
Complex models can simulate NPS pollution processes based on intrinsic physical processes. These models generally integrate hydrology, erosion and sediment processes, and chemical pollutant fate and transport. The data required for these models are commonly large. Complex models consider the estimation of runoff, soil erosion, sediment and pollutant loading based on theoretical consideration of mass, momentum, and energy [94]. Complex models can not only address a wide range of watershed-scale hydrology and water quality issues, that include short-and long-term simulations of runoff, sediment, and pollutant loads, but also supply different algorithm options for various physical processes to more accurately describe the processes mathematically. These models must be calibrated and validated carefully before applying them in management decisions because poor data input or inappropriate application of algorithms can lead to large errors [95]. These models have detailed documentation that often require extensive training and/or experience to apply correctly. These complex models require intensive parameters inputs to drive the model calculations, which are often unavailable in some areas. Careful parameter sensitivity analyses are needed prior to calibration, and the uncertainties analysis of results need to be evaluated after validation. Calibration and validation of complex models are a labor intensive and time-consuming process [96]. For application of complex models end users need expertise, sufficient training, and experience to apply these models correctly.

AGNPS/AnnAGNPS-AGNPS (Agricultural Non-Point Source Pollution Model), initially developed by USDA-ARS cooperated with the Minnesota Pollution Control
Agency and the National Resources Conservation Service (NRCS), is a single event-driven NPS pollution model [21,97]. The latest AGNPS version 5.0 was released in 2018 [98]. AnnAGNPS (Annualized Agricultural Non-Point Source Pollution Model) is the upgraded product of AGNPS and did not focus on a single rainfall event but evolved into a modeling system that can conduct watershed-scale, continuous pollutant load simulations at a daily time step. AnnAGNPS 5.0 appends a pesticide component to the model. The latest version of AnnAGNPS is v5.5 [99], which has a GIS-based DEM analysis program and a Windowsbased input editor for writing and revising of input data.
AnnAGNPS can estimate the current or long-term effects of alternative management decisions on surface runoff, sediment, and nutrients loading on a daily time step from predominantly agricultural watersheds ranging from a few hectares to 300,000 hectares [21,100]. AnnAGNPS uses homogenous land areas as square cells or any hydrologicallybased shape that represents similar soil types, land use and management, and climate to discretize a watershed. Water, sediment, nutrients, and pesticides are generated from those homogenous land areas and then routed through stream networks finally to the watershed outlet [101]. The model uses the SCS-CN equation to estimate surface runoff [102] and soil moisture content to calculate potential evapotranspiration, applies the RUSLE method to estimate sheet and rill erosion [103], and adopts the hydro-geomorphic USLE to predict sediment delivery of the stream [104]. Core components of the model include hydrology, soil erosion, sediment, and nutrient transport. The simulated nutrients include nitrogen, phosphorus, organic carbon, fertilizer, pesticides, and chemical oxygen demand (COD), and point source loads [97]. The primary inputs of the model include precipitation, soils, land use, agricultural management, upland and channel networks, point source pollution (e.g., animal feedlots, streambanks, construction sites). The model outputs runoff volume, peak flow rate, erosivity, upland and channel erosion, sediment delivery ratio, sediment enrichment ratio, mean sediment concentration, total sediment yield, and pollutant concentration on an event, monthly, or yearly basis.
AnnAGNPS has been widely used to estimate runoff water quality and NPS pollution around the world. Li, et al. [105] applied AnnAGNPS to simulate yearly streamflow and monthly nutrient loading in the Taihu Lake watershed, China. Their results showed that the AnnAGNPS model can acquire a satisfactory accuracy for annual streamflow simulation, and the accuracy of the nutrient loading simulation is moderate, and monthly nitrogen loading evaluation has higher accuracy than monthly phosphorus loading. Karki, et al. [106] to evaluate runoff, nutrients, and sediment for an agricultural watershed of 30.3 ha in East-Central Mississippi by applying the AnnAGNPS model. They indicated that AnnAGNPS can perform better for runoff evaluation than sediment and phosphorus load assessments on a longer time scale. The accuracy of the model prediction dramatically depends on high quality available data for calibration.
AnnAGNPS is flexible and reliable tool to evaluate the amount of runoff, sediment, and nutrient generated and transported throughout a watershed. It can help to identify critical source areas and delivery routes of NPS pollution, support the determination of BMP locations, and cost/benefit analyses. The predictions of sediment and nutrient loads perform better at larger monthly, seasonal or annual time scales, than shorter daily simulations [107]. However, AnnAGNPS assumes a constant nutrient loading rate for point source pollution simulation for the entire simulation period and does not account for spatial variability of rainfall [101]. AnnAGNPS does not track nutrients and pesticides attached to sediments in stream reaches from one event to the next event [100]. AnnAGNPS may underestimate daily streamflow during a dry period, as it does not account for the baseflow. Additionally, AnnAGNPS need intensive parameter inputs to support model simulations, which may lead to increased computing demand for parameter optimization, calibration, and validation.

SWAT-SWAT (Soil and Water Assessment Tool)
, developed by Dr. Jeff Arnold for the USDA-ARS, is a continuous, semi-distributed, physically-based watershed model [20]. The SWAT model is regularly updated, SWAT2012 rev.670, at the time of this publication, was released in October 2018 [108].
SWAT can evaluate surface water, crop development, sediment, nutrient yield, pesticide transport, and the impact of climate change and land management based on hydrologic inputs in a complex, ungauged watershed under various soil, vegetation, and land use management conditions [109]. SWAT operates at daily or hourly time steps and can perform long-term, continuous-time simulation. Spatial scale of SWAT applications ranges from small size watershed to an entire continent [110,111]. Within this tool a watershed is divided into multiple sub-basins, then a sub-basin is further subdivided into hydrological response units (HRUs) where all land areas have homogeneous land use, soil characteristics, and slope combinations. In each HRU, hydrological components are calculated for surface water, soil, and groundwater [109]. SWAT uses an SCS-CN method or Green & Ampt infiltration method to determine runoff volume and applies a rational formula or TR-55 method to calculate the peak flow rate. Soil erosion generated by rainfall and runoff is computed with the Modified USLE equation [112]. Land use, soil, weather and topography are the required primary input to SWAT. The output of SWAT includes water volume that sums step and accumulated surface runoff and subsurface flow, sediment yield, soil water storage, evapotranspiration, and nutrients [113].
SWAT has been widely applied to many watersheds in the world, including hydrologic modeling, non-point source pollution control, surface or subsurface water quality evaluation, groundwater modeling, soil erosion prevention and control, the impact of land use management practices, and the prediction of hydropeaking [4,31,114]. According to the official SWAT literature database, the number of articles relevant to SWAT has exceeded 1300 from 2016-2019 [115].
SWAT can be applied to various spatiotemporal scales ranging from sub-daily to yearly, and from small watershed (e.g., 10 km 2 ) to a river basin and even a continent and enables users to study the long-term environmental impacts. SWAT can obtain higher prediction accuracy while predicting on a yearly or monthly scale than on the daily scale [76]. Data needed to support SWAT simulations are readily available from various sources including governmental agencies. SWAT can also be applied to a watershed with scarce or no monitoring data [109]. SWAT has an active and influential online user community [116]. Alternative calibration and validation approaches have been developed to simplify the often time consuming and difficult calibration process [117] for example SWAT-CUP (SWAT-Calibration and Uncertainty Programs), an auto-calibration program for SWAT [95]. SWAT-CUP is available at [118]. SWAT assumes vegetation growth is insensitive to season change, which often causes a low accuracy of SWAT prediction in the dry season [119]. One solution is to divide wet and dry seasons during calibration and validation, which can efficiently improve SWAT simulation accuracy [120,121].

SWMM-SWMM (Storm Water Management Model)
, developed by U.S. EPA, is a distributed, physically-based, dynamic stormwater runoff quantity and quality model [71]. The latest version of the model is 5.1, which is a Windows-based software released in 2018. User's manual and software installable files of the model are available at [122].
SWMM is developed for the evaluation of a single rainfall event or long-term continuous rainfall-runoff processes from primarily urban areas. Users can use SWMM to: (1) design drainage system component; (2) (5) estimate the effect of land-use change on urban hydrology [123]. SWMM can track runoff quantity and quality for hourly or sub-hourly time steps. The spatial scale of SWMM applications varies from separate lots up to hundreds of acres. SWMM divides a watershed into a collection of homogeneous subcatchments as a basic hydrological unit [71]. The model considers multiple physical processes such as surface runoff, infiltration, groundwater, flow routing, water quality routing, snowmelt, and surface ponding. SWMM applies hourly or more frequency rainfall data as input. SWMM requires inputs of buildup and wash off parameters to model stormwater quality, which produces pollutographs at any point in the watershed. Simulated pollutants include TN, TP, TSS, settleable solids, oil/grease, BOD, COD, total coliform, and other user-specified pollutants [124].
The SWMM model has been commonly employed in urban drainage flooding analysis, water quality and transport of contaminants [3], TMDL implementation plans [125], urbanization and climate change effects [7], and LID effectiveness [126]. Niazi,et al. [127] presented a synthetic overview of SWMM applications and gap analyses.
SWMM can efficiently simulate hydrological and contaminant transport in complex urban watersheds. SWMM can account for time-varying rainfall during the process of simulation [71]. The model has been continuously upgraded since 1971 through present [71]. One primary limitation of SWMM is that as an analytical tool (not a design tool), it does not simulate small outlets (e.g., manhole or inlets) loss directly, but rather can be aggregated [128].

HSPF-HSPF (Hydrological Simulation
Program-Fortran), jointly developed by the U.S. EPA and the USGS, is a comprehensive, continuous, physically-based watershed hydrology and water quality model [19]. The latest version 12.2 of HSPF can be downloaded via the BASINS model from [129]. This version is WinHSPF 3.0 that integrated into BASINS as a core module. The users can also access a standalone HSPF version 11.0 from USGS website [130].
HSPF is used to simulate water quantity and quality processes, conventional and toxic organic pollutants loads within a natural and developed watershed, and predicts nutrients fate and transport routing in-stream and well-mixed lakes and impoundments [19]. The simulated time scale of HSPF is from a few minutes to several hundred years by using time steps ranging from sub-hourly to daily. The spatial extent of HSPF application varies from a few acres to a large watershed (the Chesapeake Bay with roughly 160,000 km 2 ) [131]. In HSPF, a basin is divided into land units that can reflect a homogeneous hydrologic and water quality response [132]. HSPF has three primary modules: PERLND, IMPLND, and RCHRES. The PERLND controls runoff and water quality from pervious areas; the IMPLND module simulates water quantity and quality on impervious land segments; the RCHRES module reflects the route of flow and water quality constituents from the PERLND and IMPLND modules [133]. HSPF needs continuous time-series records as input including precipitation, potential evapotranspiration, air temperature, wind, solar radiation, humidity, tillage practices, and point sources. Water quality processes simulation also needs chemical pollutant (e.g., pesticide or fertilizer) application data. HSPF outputs include flow rate, sediment yield, nutrients, pesticides, toxic chemicals, and other water quality variables [134].
An early HSPF application literature summary can be found at [135]. Kim,et al. [5] integrated HSPF with a maximum likelihood filter to improve water quality forecasts in the Kumho catchment (23,384 km 2 ) in South Korea. Huo, et al. [136] used HSPF to evaluate nonpoint source water quality in the Feitsui reservoir watershed (303 km 2 ) in Taiwan.
HSPF is a comprehensive watershed model for agricultural or urban areas. HSPF adopted the flexible module design, thus users can use different simulation algorithms (empirical or physical) to analyze NPS pollution processes based on how much data available [132]. HSPF is suitable for predicting yearly and monthly streamflow and sediment yield, except for the months with extreme weather conditions. Daily streamflow simulations are reasonable except for flood events. HSPF is able to adapt to a wide range of watershed conditions involving various surface water and groundwater quantity and quality problems at multiple spatiotemporal scales [137]. The limitation of HSPF is that model calibration process requires expertise and experience to determine appropriate parameters, currently available documentation provides no uniform parameter estimation guide [131]. The data requirement of the model is extensive [138]. Due to use of a nonlinear flow or storage-based equation in routings, HSPF cannot simulate an intense single-event storm or flood waves, especially for a large sub-basin and long channel [33].

WAM-WAM (Watershed Assessment Model), developed by Soil and Water
Engineering Technology (SWET), Inc., is a GIS-based, deterministic, physically-based watershed-scale hydrology and water quality model [43]. WAM works as an extension in the ArcGIS environment. The latest WAM supports ArcGIS 10.4.1. Users can download WAM and its associated documents from the link [139].
WAM can estimate the hydrology and water quality response of land-use changes and associated management practices within a complex watershed including agriculture/urban land uses, natural or man-made stream networks, looped flow networks, multiple hydraulic structures, springshed groundwater systems, and tidally influenced boundary conditions [43]. It can run on both daily and hourly time steps to estimate surface water flow, groundwater flow, and pollutant loads. In WAM, a watershed is discretized into many single cells that represent the basic unit of hydrology and water quality simulation. WAM integrates GLEAMS (Groundwater Loading Effects of Agricultural Management Systems) and EAAMOD (Everglades Agricultural Area Model) to simulate soil and plant processes that affect water quality [140]. The calculations of daily surface and subsurface flow, nitrogen and phosphorus concentrations take place in each cell. The model then routes the surface and subsurface flow from cells to estimate the flow and phosphorus levels throughout the watershed by using the Basin Land Area to Stream Routing model (BLASROUTE). Water quality variables include TN, TP, TSS, and BOD. Primary inputs of the model include land use, soils, topography, climate data, and point source data. The outputs include surface and groundwater flow, pollutant loads, a ranking of land use by load source, daily time series of flows and pollutants, and displays of different BMPs in table, graph, and map [43].
The applications of WAM include the simulation of constituents that are important to predict eutrophication processes in water bodies, analysis of the hydrological and water quality effects on rivers and lakes for management scenarios, view and estimation of the simulated flow and concentrations for every source cell and stream reach, and prioritization of BMPs zones [141]. Bottcher [142] built a Suwannee River watershed assessment (25,770 km 2 ) model based on the early ArcView version WAM, displayed spatial loads of soluble N, soluble P, TSS, and BOD, and simulated flows, total P, and total N on the daily and seasonal time scale. Bottcher, et al. [140] applied WAM into a predominately agricultural watershed to hydrologic and hydraulic processes and NPS pollution loading simulation with 2025 and 2055 future scenario analysis. Corrales, et al. [13] used WAM in the northern Lake Okeechobee watershed (4150 km 2 ) to evaluate total P load at the source area, streams, and outfall levels.
WAM can simulate complicated surface and subsurface hydrologic processes and nutrient loading. It can describe spatial and hydraulic details and is flexible to accommodate varied hydrologic, water quality, land and water management processes, and support scenarios analysis. WAM accounts for the pollution contribution of upland landscape with deep water tables, lands with shallow water tables with and without artificial drainage, and wetlands, urban areas, and mining sites [141]. However, the drawback of WAM is limited numbers of applications based in Florida and it requires intensive physical characterization data, which might be hard to obtain for some places. WAM cannot simulate small-scale and short-term storm event impact.

Modeling Systems
A modeling system uses the concept of multiple modules to independently maintain and apply separate model structures or information to carry out complex decision analyses or synthesis. It integrates databases, tools, techniques, and models into a GIS platform. These data, tools, methods, and models have a close linkage and work together to perform various environmental simulations and analyses at multiple time and space scales, and different modules comply with the same data exchange protocols. Thus, the modeling system comprehensively uses compatible data types, various tools, and different models to assist in systematically solving complicated watershed water quality and NPS pollution issues. It can serve as a multipurpose decision support platform. The modeling system can simulate watershed-scale, hydraulic, hydrologic, water quality, and NPS pollution issues. However, the data requirement and model computation is enormous. These models have a steep learning curve for users as many individual models and various watershed analysis tools are integrated. This can involve extensive pre-processing and postprocessing of data and output. The user commonly needs extensive training prior to developing or running these models.

AGWA-AGWA (Automated Geospatial
Watershed Assessment Tool), codeveloped by the U.S. EPA, USDA-ARS, the University of Arizona, and the University of Wyoming, is a GIS-based distributed, light level hydrology modeling system [143]. AGWA underwent a number of upgrades from AGWA 1.5 for ArcView 3.x, AGWA 2.x for ArcGIS 9.x to AGWA 3.x for ArcGIS 10.
x. The user's manual, theoretical documentation, training, and software of AGWA are available at [144].
Using readily attainable GIS data, AGWA can evaluate the impacts of land-use change on soil erosion, water quantity and quality, and watershed-scale NPS pollution at different spatial and temporal scales, ranging from a small drainage area to a large size watershed [39]. AGWA facilitates the processes of parameterization, model execution and outcome visualization and packaged inside multiple sub-models such as RHEM (Rangeland Hydrology and Erosion Model), KINEROS2 (Kinematic Runoff and Erosion Model), KINEROS-OPUS, SWAT2000, and SWAT2005. RHEM is a rangeland soil erosion model by water for a single rainfall at hillslopes. RHEM is integrated into the overland flow component of KINEROS2. KINEROS2 is a physically-based, event-driven model that can simulate vegetation interception, soil infiltration, surface runoff, and soil erosion in small watersheds. KINEROS-OPUS is a combination of KINEROS and OPUS2 (not an acronym) with more sophisticated functionalities including simulations of runoff, sediment transport, nitrogen and phosphorus cycling, and chemical transport processes. Comparatively, SWAT is a hydrology and water quality model for long-term watershed modeling, and details on SWAT have been presented in Section 2.3. The modeling inputs of AGWA include DEM, land use, soil, and precipitation data. The outputs of the KINEROS sub-model include channel infiltration, plane infiltration, runoff volume, sediment yield, peak flow, channel erosion, and sediment discharge. The outputs of SWAT sub-model have runoff volumes, evaporation, percolation, transmission losses, water yield, sediment yield, and nitrate and phosphorous in surface runoff [145].
The applications of AGWA involve watershed management, water resource assessment, land-use change detection, and ecology evaluation. Yang and Li [146] employed AGWA-SWAT and KINEROS to evaluate the hydrological response of urban development in a Panther Creek watershed (94.2 km 2 ) with future land-use scenario analysis. Jones, et al. [15] used AGWA-KINEROS2 to estimate sediment delivery on the post-fire land cover to Strontia Springs Reservoir during a rainstorm with alternative scenarios.
AGWA provides a repeatable method to facilitate the setup and execution of multiple submodels. AGWA supports future scenario analyses, decision assist, and hydrology and soil erosion simulations at different spatiotemporal scales [143]. Further, AGWA can also supply a quick post-fire watershed evaluation by replacing the part of an existing land use map with a burn severity map [147]. However, AGWA does not integrate within the latest version (SWAT2012 or SWAT+) of SWAT. The tool does not consider the data acquirement component, so users cannot collect online data via AGWA.

BASINS-BASINS (Better Assessment Science Integrating Point and Nonpoint
Sources), developed by the U.S. EPA, serves as a multipurpose environmental analysis system [148]. The latest version 4.1 of BASINS was released in 2013. It is based on a userfriendly open-source MapWindow GIS and is available at [129].
BASINS is suitable for watershed management, water quality analysis and TMDLs development, and integrates environmental data, analysis tools, and various watershed and water quality models. BASINS facilitates examination of environmental issues, support analysis of watershed systems, and provides a platform for inspecting management alternatives [134,149]. As a comprehensive watershed model framework, BASINS integrates several watershed models such as HSPF, SWAT, SWMM, GWLF-E, PLOAD (Pollutant Loading Estimator), and instream and water quality models such as AQUATOX and WASP (Water Quality Analysis Simulation Program) as plug-ins and some analysis tools such as manual/automatic watershed delineation, land use reclassification, lookup tables, shapefile editor, PEST (the Parameter Optimization Program), time series functions, CAT (Climate Assessment Tool), GenScn (GENeration and analysis of model simulation SCeNarios), and WDMUtil (Watershed Data Management Utility) into a unified GIS interface [134,148]. The overarching objectives of BASINS are to: (1) identify and prioritize water bodies; (2) evaluate magnitude and potential significance of point and nonpoint source pollution; (3) simulate point source and nonpoint source nutrient loadings and fate and transport processes; (4) evaluate the relative value of potential control strategy; and (5) visualize environmental conditions to the public through tables, graphs, and maps [148]. Core data in BASINS includes four main types: base cartographic data, environmental background data, monitoring data, and point source data. The input data include DEM, land use, soil, weather data, and point source data, and output files include maps, graphs, and tables summarizing point and non-point source pollution in a watershed [148].
BASINS has been used to develop solutions for solving real-world problems in the past two decades. Saleh and Du [134] compared the application results of BASINS-SWAT and HSPF in the Upper North Bosque River watershed. HSPF performed better than SWAT in terms of evaluating daily flow and sediment; however, SWAT was a much better predictor for simulating daily and monthly nutrient loadings. Crossette, et al. [149] applied BASINS-HSPF to the data-scarce Shebelle watershed (2400 km 2 ) in central Ethiopia and presented the detailed steps of the BASINS application.
BASINS is a multipurpose environmental model system used to conduct watershed hydrology and water quality studies and develop TMDLs for water quality impaired water bodies [150]. The model system facilitates watershed and water quality studies through decreasing data collecting and processing time, reducing execution steps, and minimizing error caused by incompatible data format [149]. Due to integrating different models, BASINS can simulate water quality and NPS pollution issues at various spatiotemporal scales. Moreover, the model can analyze and develop TMDL guidelines that meet the need of the Clean Water Act. However, BASINS has a steep learning curve because of the involvement of environmental theory and technical knowledge.

WMS-WMS (Watershed Modeling System)
, developed by AQUAVEO Inc., is a comprehensive GIS-based modeling system for watershed hydrology and hydraulics [151]. The latest version of WMS is 11.0 and is available at [44]. WMS is not a public domain software.
WMS is designed to simulate hydrologic, hydraulic, storm drain, sanitary sewer, water distribution, and NPS pollution processes, including almost all phases of hydrology and hydraulics. WMS integrates eight modules under one GIS interface, including terrain, drainage, map, hydrologic, hydraulic (river), GIS, 2D grid, and 2D scatter point. Each module corresponds to one of the primary data types or modeling environments supported by WMS. WMS supplies a GIS-based framework to operate different models such as HEC-1, HEC-HMS, NSS, TR-20, TR-55, Rational Method, OC (Orange County, California) Rational, OC hydrographic, HSPF, SWMM, XP-SWMM, SMPDBK (Simplified Dam-Break Model), GSSHA (Gridded Surface Subsurface Hydraulic Analysis), CE-QUAL-W2, HY-8, HY-12, Hydraulic Toolbox, and EPANET [17]. The model integrates national streamflow statistics, supports conversion and comparison of the results from different sub-models, displays a comparison between observed and simulated hydrographs, provides terrain surface viewing, and exports images for reports and presentations [44].
WMS has been applied to watershed management [152], flood hazard analysis [153], water quality evaluation, and groundwater simulation [29]. Erturk,et al. [152] used WMS 7.1 for watershed management and NPS pollution modeling in Koycegiz Lake-Dalyan Lagoon watershed (1200 km 2 ) located at the southwest of Turkey. WMS supports a flexible watershed delineation method. For example, users can also control watershed boundaries created by WMS. The model system has detailed documentation along with an active user community [44]. However, it is not a public domain software, and users need to acquire a WMS license for use.

Model Selection
As described above, the NPS pollution model selection is rather challenging for model users based on the fact that numerous different watershed models are available currently. The users have to consider multiple factors while selecting an appropriate model for a specific project, including the nature of watershed issues that need to be solved, the processes that are interested to simulate, the complexity of the model, desired spatial and temporal scales, data requirement, expected output and simulation accuracy, user's knowledge and skills, and the budget of the project.
The desired output or management information needed should dictate the selection of a particular model. Some models are appropriate for estimating numerical values with minimal input while others are designed to inform complex decisions with intensive data input from multiple sources. For example, some tools may provide a quantitative estimate of the areal average annual pollution load from a watershed with existing data, while others are designed to inform a management practice or scenarios such as optimization of BMPs or determination of TMDLs for a watershed. A simple model or a medium complexity model may provide sufficient results for a specific unknown value in the former situation, but the latter problem might need to adopt a complex model or a modeling system combined with expert practitioners to apply multiple data inputs and management insight to generate meaningful or useful model output. The classification of NPS pollution models described here is a straightforward way to begin the model selection process and should help users choose the appropriate model or tool based on an understanding of how the models function and what the models require. Moreover, a detailed categorization and identification of individual model structure and function is helpful for those tasked in model selection in practice, but also provide insight to the manager charged with providing useful decisionmaking information. Table 2 lists the main features of each model discussed in the article to identify model selection criteria and compare different models.
Each NPS pollution model has unique function and characteristics for specific purposes. The more particular processes of simulation users expect, the less the range of model selection is. If a user needs to know the environmental contamination of mercury in a watershed, he/she might not face many difficulties from model selection since few models integrate mercury transport/transformations algorithm and codes. WARMF is a good option for that specific simulated objective. However, users need to be aware that each model is a simplification of the real world. Any simplification implies some physical processes are omitted in the model design. If the ignored processes are of significant importance for modelers, the user should investigate alternative approaches or at least be aware of the limitations of the selected model or tool.
Generally, the NPS model performs better on smaller spatial scales, and longer time steps (i.e., monthly or yearly) [154]. Some NPS pollution models can conduct an event-driven pollution simulation such as AGNPS, SWMM, and SLAMM; while others are more appropriate for the long-term continuous simulation including AnnAGNPS, GWLF, LSPC, L-THIA, N-SPECT/OpenNSPECT, WARMF, and SWAT. Few models can conduct both short-and long-term NPS pollution simulation such as HSPF, BASINS, WMS, and WAM. A lumped parameter, empirical model can quantitively calculate long-term average gross NPS pollution loads, but it does not consider the fate and transport, and spatial route of nutrients loads. A distributed, physically-based model presents the spatial distribution of pollutants, and can account for fate and transport of nutrient, but needs a lot of input parameters. Therefore, users should select a suitable watershed model in light of the watershed size and time duration of a simulated event.

Spatial and Temporal Considerations
Spatial and temporal scales affect many aspects of NPS pollution model applications such as model selection, watershed discretization, data pre-processing, as well as simulation accuracy [155]. Firstly, model selection will be determined mostly by the space and time scale of the model application. From the spatial perspective, the main transport route of runoff, sediment and nutrients is dominant by river networks system at a river scale basin. A large size watershed is less sensitive to short-duration, high-intensity rainfall, and the model structure will determine how well the model simulates these patterns with given lags and routing. In contrast, a small watershed has a sensitive response to high intensity, short duration rainfalls because overland runoff and streams control transport route of sediment and nutrients transportation and may be less sensitive to model structure or parameterization [147]. From the temporal point of view, NPS pollution processes may occur at any time scale (i.e., event, daily, monthly, or yearly), therefore different watershed-scale NPS models meet the requirement of the simulation of NPS pollution scenarios at various time scales. An event-driven NPS pollution model is typically used for the simulation of a short-duration, intensive rainfall process. A physically-based, long-term continuous-simulation NPS pollution model can estimate various hydrologic processes under multiple time scales, typically execute on a time span range from a few minutes to hundreds of years and its outcome can present a time-series of runoff, sediment, and pollutants loadings. Spatial and temporal considerations of NPS pollution model are discussed in detail by Baffaut,et al. [156].
Secondly, watershed space can be discretized into a grid cell, any shape subarea, HRU, or sub-watershed that represents the smallest spatially computed unit of the model. AnnAGNPS users are recommended to use a cell size of 40 acres to operate the model when the size of a watershed exceeds 2000 acres [157]. Much of the measured data related to NPS pollution are obtained from observations and artificial experiments at a laboratory or field, which is useful to verify understanding of NPS pollution underlying processes. However, when these experimental data are applied to watershed-scale NPS pollution simulation, it results in additional model uncertainty because they cannot correctly reflect the pattern of NPS pollution at a larger scale. If the selected scales do not match between observed data and model output variables, then upscaling or downscaling method might be applied to solve this problem. Upscaling refers to aggregate data from a smaller scale to represent a larger scale; and downscaling is used to disaggregate data at a larger scale to suit a smaller scale [158].

Calibration and Validation
The primary purpose of calibration is to build a mapping relationship between an NPS pollution model and the physical real-world. The target of calibration is to minimize the model error between simulated results and observed data by adjusting selected input parameter values. The calibration process will eventually find the optimized parameter combinations that make the model obtain higher accuracy and less uncertainty. Validation has the same process with calibration except for using an independent dataset from a different period and keeping calibrated-well parameters unchanged [159]. The purpose of validation is to assure that the calibrated model can produce properly evaluated results under similar hydrological conditions with calibration.
Simple models do not require calibration and validation processes (i.e., L-THIA, N-SPECT) since these models have been sufficiently verified at a specific region during the processes of model development. Simple models are typically used as a quick evaluation tool to present a relative rough estimation. Validation of these models is recommended when they are applied to other areas where they have not been verified. Lim,et al. [160] developed an automatic calibration system to search the optimized CN combinations while estimating runoff and pollutant in the Litter Eagle Creek watershed in Indiana. The authors found that the simulated results of hydrology and water quality can be significantly improved after calibration. Simple models are limited to calibrate SCS-CN and/or EMC parameters due to the simplicity of the model structure. On the other hand, the execution of complex NPS models or modeling system (i.e., AnnAGNPS, SWAT, SWMM, SWAT, WAM, BASINS, WMS) involves many parameter inputs, but most of the parameter values are not exactly known due to spatial variability and direct measurement is unavailable [161]. Thus, the model needs to obtain these parameter values through calibration before it may be practically applicable as a decision-making tool.
The modeler conducts parameter calibration processes by adjusting parameter value and running the model repeatedly. The calibration methods of a model have typically two: manual or automatic. The manual calibration depends on the experience of the modeler because it is difficult to search and try all parameter combinations. In most cases, the manual calibration process will cease once statistics from the objective function satisfied the target set by the modeler before running the model. The operation of manual calibration is complicated for non-expert users. In contrast, automatic calibration is an iterative procedure that can try more parameter combinations by using a grid search or other algorithms and setting searching value range of different parameters. Some standalone and professional automatic calibration procedures are developed for the specific model to facilitate the application of NPS pollution models and shorten the time of calibration. SWAT Calibration and Uncertainty Programs (SWAT-CUP) is a standalone, robust, and public software program to automatically calibrate parameters from SWAT model, which integrates Sequential Uncertainty Fitting 2 (SUFI2), Generalized Likelihood Uncertainty Estimation (GLUE), Markov Chain Monto Carlo (MCMC), Particle Swarm Optimization (PSO) and Parameter Solution (ParaSol) procedures [162]. The number of parameters that can be calibrated by SWAT-CUP is over 700. A detailed introduction of various calibration procedures with SWAT-CUP was provided by Yang, et al. [163]. SWAT-CUP facilitates sensitivity analysis, calibration, validation, and uncertainty analysis that almost involves the entire post-processing of SWAT simulation. A more detailed discussion on the use, calibration and validation of the SWAT model is provided by Arnold, et al. [109]. Additionally, HSPF Parameter (HSPFParm) database support software was designed for the HSPF model to identify initial value and expected ranges of parameters before running HSPF calibration progress. Expert System for Calibration of HSPF (HSPEXP) was developed to estimate statistical error for HSPF simulation of each time input parameter, which also provides specialist advice about optimized parameter combinations to improve the calibration [132]. According to Duda,et al. [131], HSPEXP procedure only works for manual hydrology calibration, not for water quality calibration, and does not consider snow simulation.
Notwithstanding, most NPS pollution models need to be calibrated but there are few professional and standalone procedure like SWAT-CUP and HSPEXP to assist the calibration process of other models. Moreover, even with the presence of an automatic calibration method, the model calibration is still a time-consuming process because of multiple involved procedures such as the preparation (formatting) of calibration file, parameter sensitivity analysis and iterative algorithm running. It should be noted that the calibration is often a hierarchical process beginning with flow volume, followed by sediment loads, and finally followed by NPS pollutant loading rates [131]. In practice, the calibration of nutrients concentration or loads is typically challenging to be accurate because of the existence of model uncertainty as well as a lack of observed pollution data in practice. According to Moriasi,et al. [158], model performance can be evaluated as 'satisfactory' if PBIAS is in range of ±70% for N and P while other statistical indicators can also satisfy the criterion. Detailed model calibration and validation strategies and steps are beyond the scope of this review article, but more calibration and validation information can be obtained from published literature [90,109,131,161,164].

Uncertainty Analysis
Uncertainty is a natural and inherent deficiency of models. Uncertainty of NPS pollution models comes from three aspects: structural uncertainty, data uncertainty, and parameter uncertainty [164]. The conceptual design of the model is simplified so it brings structural uncertainty. These simplifications or assumptions neglect some physical processes occurring in a watershed. For example, wind erosion and landslides processes are typically ignored at NPS pollution models. Additionally, it causes uncertainty of simulated results because the mechanism and process of pollution fate and transport in the stream might not have been entirely and correctly described in currently available NPS pollution models. Structural uncertainty also involves some unknown occurrences to a modeler although processes affecting water quality have been included in the model such as upstream reservoirs, dam operation, and farm fertilization [162]. In addition to structural uncertainty, some uncertainties come from errors of input data such as DEM, soils, rainfall, temperature, land use and point measurements, these data are used in complex models, but their values are rarely exactly known. For example, spatial uncertainty of rainfall could be large different especially in mountain areas, as well as the resolution of DEM and measuring error of point data have a significant impact on model uncertainty [165]. Finally, parameter uncertainty reflects the parameter non-uniqueness issue. Clearly, the optimized solution for model calibration is never unique and different parameter combination may output very similar simulation results [95]. Therefore, it is crucial for the modelers to clearly understand the intrinsically physical meanings and value ranges of major parameters at specific research areas. Furthermore, uncertainties from the interaction between model structure, input data, and parameters cause the more complex uncertainty of evaluated results from NPS pollution simulation.
Uncertainty analysis is a process to estimate the effect of structural uncertainty, data input and parameters on model results [164]. The evaluation of model structural uncertainty can be obtained by analyzing statistic characteristics from the error series of the objective function under the assumption without input and parameter uncertainty existing. The typical method is to compare the simulated error series of the same study area with different models [137]. Input uncertainty involves the spatial distribution of rainfall, the resolution of DEM and measurement errors of point data. The impact of rainfall input uncertainty on model output focuses on the adjustment of the density of rainfall stations and their spatial distribution by comparing different spatial interpolate methods. The influence of DEM resolution uncertainty on model output is studied by changing DEM grid size and observing the change of extracting terrain parameter and model output. Not much attention has been paid on how point data measurement error affects the model uncertainty so far [166]. The research method is to treat point data input as a random input, to get a prior probability distribution of measured values, which is then used to obtain the posterior probability distribution through statistic likelihood function by comparing this distribution to get model uncertainty due to input data errors. Another method is to add interference random number that follows with normal distribution into the original input series to observe the difference of model output. Currently, more researches focus on the impact of parameter uncertainty on model results. The research methods of parameter uncertainty involve parameter sensitivity analysis and probability analysis [163]. Two conventional methods of sensitivity analysis are: local (or One-at-a-time, OAT) and global (or All-at-a-time, AAT). OAT creates variability in one input parameter to determine the output result, while fixing remaining other parameters. In contrast, AAT can show the variance of all input parameters, with output results uncertainty averaged over the input parameters. The probability analysis method is used to describe the parameter uncertainty of a model. Since uncertainty is a structural and inevitable characteristic of NPS pollution models, uncertainty analysis is an important step that helps appropriately select and apply a model. It should be mentioned that in model evaluations, there exists an uncertainty range that shows to which extent the model can capture the observed signal.

Summary and Future Research Direction
This review article summarizes 14 watershed-scale NPS pollution models associated with essential model characteristics such as intended use, components/capabilities, model structure, applicable land use, spatial and temporal scale, availability as well as the strengths and limitations for each model. The categorizations, simple model (L-THIA, N-SPECT/ OpenNSPECT), medium complexity model (GWLF, LSPC, SLAMM, and WARMF), complex model (AGNPS/AnnAGNPS, SWAT, SWMM, HSPF, and WAM), and modeling system (AGWA, BASINS, and WMS), are used to guide users to identify the most appropriate model for his/her project. A simple model relies on minimal data, does not need rigorous calibration, and typically serves as a quick screening tool for NPS pollution evaluation. A medium complexity model is generally applicable topoint and nonpoint source pollution problems in diverse watersheds. These medium complexity models require detailed data input but minimal calibration requirements. They are commonly used as a compromise between simple models and complex models. A complex model can represent the complex processes of NPS pollution at different spatial and temporal scales and requires intensive data input. Parameters sensitivity analysis, calibration, validation, and uncertainty analysis need to be conducted while using complex models. A modeling system works as a multipurpose decision support platform that can provide an entire solution of watershed NPS pollution issue, but users need to possess profound professional knowledge and technique for appropriate use.
Omission or absence of any of the many useful and valuable models from our review should not be interpreted as a judgement and we encourage the readers to investigate beyond the limited number of models discussed here, including a recently published government report that describes popular models based on their publication frequency in scholarly literature (see [167]).
Model selection is a challenging task and is limited to different factors such as the nature of the watershed issue, the processes that need to be simulated, desirable spatial and temporal scale, data requirement, project cost, and so on. Table 2 is helpful to model selection by demonstrating a detailed classification and comparison of multiple models with key characteristics. Among these models, AGNPS, SWMM, and SLAMM are applicable to an event-driven simulation; AnnAGNPS, GWLF, LSPC, L-THIA, N-SPECT/OpenNSPECT, WARMF, and SWAT suited for a long-term continuous simulation; HSPF, WAM, BASINS, and WMS can conduct both short-and long-term simulation.
It is essential for most models to be well-calibrated and validated before they can be practically applied for the planning and management of watershed water resources. The manual calibration method is not recommended to non-expert users. Automatic calibration can significantly reduce model uncertainty, shrink model calibration time, and quickly find optimization parameter combinations. However, most NPS pollution models do not have a standalone and professional automatic calibration procedures for model calibration, validation, and uncertainty analysis. Therefore, it will be a promising direction to develop a universal, cross-model, and automatic calibration procedure that integrates advanced algorithms of parameter search.
Since uncertainty is an inherent and inevitable characteristic of all models, uncertainty analysis is still a frontier area in NPS pollution model future research. Current studies mainly focus on parameter uncertainty, so it should be strengthed on studies of reducing model uncertainty that stem from the model structure and data input, including to develop new research measures and innovate uncertainty evaluation techniques. Additionally, the future direction of NPS pollution models will involve strengthening the mechanism study of NPS pollution processes, especially for some specific pollutants, making the mathematic model more correctly describe the pollution processes. Furture work will also invovle developing new methods to improve the application accuracy of NPS pollution models, particularly in areas with limited data. Classification of water quality and NPS pollution models. Predict runoff and erosion rates on rangelands [39]; Conduct rapid, post-fire watershed assessment [147] Does