Currents Status, Challenges, and Future Directions in Identifying Critical Source Areas for Non-Point Source Pollution in Canadian Conditions

: Non-point source (NPS) pollution is an important problem that has been threatening freshwater resources throughout the world. Best Management Practices (BMPs) can reduce NPS pollution delivery to receiving waters. For economic reasons, BMPs should be placed at critical source areas (CSAs), which are the areas contributing most of the NPS pollution. The CSAs are the areas in a watershed where source coincides with transport factors, such as runo ﬀ , erosion, subsurface ﬂow, and channel processes. Methods ranging from simple index-based to detailed hydrologic and water quality (HWQ) models are being used to identify CSAs. However, application of these methods for Canadian watersheds remains challenging due to the diversiﬁed hydrological conditions, which are not fully incorporated into most existing methods. The aim of this work is to review potential methods and challenges in identifying CSAs under Canadian conditions. As such, this study: (a) reviews di ﬀ erent methods for identifying CSAs; (b) discusses challenges and the current state of CSA identiﬁcation; and (c) highlights future research directions to address limitations of currently available methods. It appears that applications of both simple index-based methods and detailed HWQ models to determine CSAs are limited in Canadian conditions. As no single method / model is perfect, it is recommended to develop a ‘Toolbox’ that can host a variety of methods to identify CSAs so as to allow ﬂexibility to the end users on the choice of the methods.


Introduction
Deterioration of water quality has become a leading problem as human populations continue to grow, and industrial and agricultural activities are expanding [1][2][3]. Furthermore, climate change threatens alterations to the hydrological cycle that in turn will have serious implications on water quality [4][5][6]. Impacts are expected to be greater in cold and temperate climate regions, which are characterized by extreme temperature, frozen ground, freeze-thaw cycles, snowpack accumulation and melting, and mixed precipitation patterns (e.g., [7,8]). A leading water quality problem on a global scale is eutrophication due to high nutrient loads (i.e., phosphorus, P, and nitrogen, N) from agricultural land use activities [3,9,10]. Waters are classified as eutrophic when P levels exceed 16 mg/m 3 and N concentration is higher than 390 mg/m 3 [11]. The eutrophication 'crisis' of the North American Great The PI method was first developed in the early 1990s by the U.S. Department of Agriculture (USDA) [56]. The PI method estimates the risk of P loss from farmland through a simple combination of several factors related to source and transport conditions [28]. Different versions of PIs may vary in terms of the type of factors included and the way they are combined. In the original PI (Lemunyon and Gilbert [56]); however, soil P and rate, method, timing, and type of P applied are used as source factors. Runoff, erosion, and proximity to streams are used as transport factors. The P loss risk is calculated based on a simple arithmetic computation of all factors, and then usually is divided into four categories (low, medium, high, and very high). Since its initial development in the early 1990s [56], the PI method has been evolving to suit site specific conditions [57][58][59]. The PI method is the most widely applied method throughout the world [28]. It has been applied over the U.S. (in 47 States, e.g., [42,58,60]), Canada (e.g., [59,61,62]), China (e.g., [63]), Europe (e.g., [64,65]), Denmark (e.g., [66]), Sweden (e.g., [29]), Australia (e.g., [67,68]), and New Zealand (e.g., [69,70]). The broad acceptance of the PI method over wider regions is associated with the numerous advantages, which have been summarized by Maguire et al. [71]. It is relatively low demand for input data and manpower, compared to process based models, is the foremost factor [28]. Detailed reviews of the PI applications can be found from past review papers, such as Buczko and Kuchenbuch [55] for applications in the U.S. and European watersheds, and Heckrath et al. [44] on Nordic watershed applications. Both source and transport factors are considered in most of the studies. Importantly, some studies considered connectivity of the landscape to the stream system in applying the PI method (e.g., [59,61,66,68,[72][73][74][75]).
In Canada, development of a countrywide PI was first initiated by Agriculture and Agri-Food Canada (AAFC) in 1993 under the Agri-Environmental Indicator Project [76]. In 2000, the first adapted version, Indicator of Risk of Water Contamination by Phosphorus (IROWC-P), was developed and applied only for Quebec [77]. The IROWC-P is based on a P indexing system initially developed by Lemunyon and Gilbert [56] and adapted to Canadian conditions [77]. The IROWC-P was designed at the Soil Landscapes of Canada (SLC) polygon scale and considers P saturation, P soil test, excess, manure, and fertilizer as source factors. Soil erosion and overland flow are considered as transport factors. The IROWC-P indicator assigns a risk class (very low, low, moderate, high, and very high) to each Soil Landscapes of Canada (SLC) V3.0 agricultural polygon [78,79]. Van Bochove et al. [23] have suggested improving the transport component of the IROWC-P by incorporating additional transport factors to represent Canadian conditions. They proposed to include more transport factors for particulate P (erosion rate), dissolved P (surface and subsurface runoff), and connectivity (topographic index, tile drainage, preferential flow, and surface drainage density). Since then, several researchers have applied the modified IROWC to assess the vulnerability of Canadian watersheds Agriculture 2020, 10, 468 4 of 25 (e.g., van Bochove et al., [79][80][81]). In addition, Reid [59] has developed a PI for Ontario, while Schendel et al. [61] used the original PI for British Columbia, Canada.

Topographic Index (TI)
The TI method is primarily based on the idea that topography has strong control on spatial variations in hydrological conditions [82]. The concept of the TI was initially suggested as an indicator for surface runoff generating areas by Kirkby [83]. It has been used as the basis for the rainfall-runoff model called TOPMODEL [84]. The use of the TI as an indicator for P loss has been then suggested by Endreny and Wood [85]. The most commonly used form of the TI is defined as ln(α/tanβ), where α is the upslope contributing area to a given point in the catchment, and β is the local surface slope angle [86].
Even if the TI is a good indicator of the occurrence of surface runoff, high surface runoff does not necessarily indicate high P source availability. Therefore, the TI alone cannot accurately estimate CSAs for P loss, which is a function of both the occurrence of runoff and the potential for mobilization of P at a location [28]. The TI can only quantify 'potential' CSAs, but not necessarily actual CSAs. Thus, the method should be used with caution. In response, researchers suggested to use the TI in conjunction with other methods that can indicate source availability (e.g., [85]).
Though it is not as popular as the PI, the TI also is being applied in several places. The TI is mainly used in combination with other approaches that includes source factors. Endreny and Wood [85], for example, combined a topographic index with landscape based export coefficients to identify locations of pollutant loading risk areas. Similarly, Heathwaite et al. [87] combined the TI with field information to assess pollutants loss risk areas. Van Bochove et al. [23] incorporated the TI into the modified version of the PI for Canada called IROWC-P [77] to improve IROWC-P prediction capability in areas dominated by saturation excess runoff generation processes. Others combined the Generalized Watershed Loading Function (GWLF) model [88] with the TI to modify the USDA Soil Conservation Service Curve Number (SCS CN) equation for variable source area hydrology (e.g., [89][90][91]). Meals et al. [92] coupled TI with a spatially explicit mass-balance model to estimate CSAs for P loss over time for a Vermont watershed. Moreover, Djodjic and Villa [93] used the modified Unit Stream Power Erosion Deposition (USPED) model in combination with a high-resolution DEM to identify erosion-and overland flow-prone areas on arable land. In the USPED model, slope length (L) and steepness (S) factors of the Universal Soil Loss Equation (USLE) [94] and the Revised Universal Soil Loss Equation [95] are replaced by upslope contributing area.

Other Miscellaneous Methods
There are also several miscellaneous methods, which do not fall within either the PI or the TI methods, but are available for use to determine CSAs. Some of the simple inventory methods that fall within this category are the export coefficient (e.g., [96]), hydrologic unit codes (HUCs) (e.g., [97]), curve number (e.g., [98]), and others. They are mostly source or availability-based models, but they often do not consider transport factors. Such methods, therefore, identify 'potential' CSAs, but not actual CSAs. In reality, however, actual CSAs are affected by both source and transport factors. These methods, therefore, should be used cautiously while identifying CSAs. In the past, these methods have been used in several countries to identify CSAs (e.g., [49,[96][97][98][99]). More importantly, a method based on export coefficient modeling has been developed and applied in several watersheds in Ontario, Canada, to assess P loadings [100]. In this method, P export from cropland is based on a relation between soil loss estimated using the USLE and P export. They modified export coefficients to reflect site-specific soil and runoff characteristics for those land use classes that are most often changed with urbanization. Similarly, a regression model called SPAtially Referenced Regressions on Watershed attributes (SPARROW) model has been applied in the Red and Assiniboine river basins to identify hot spot areas of nutrient loadings [101].

Detailed HWQ Models for CSAs
The other approach that has been used by many researchers is the use of HWQ models to identify CSAs (e.g., [50]). Currently, numerous HWQ models are available for use in water pollution studies. Shoemaker et al. [102] provided a list of HWQ models that are currently available for use. In this section, however, only a few of them that have better potential for Canadian conditions were selected and reviewed. Booty and Benoy [103] assessed several HWQ models based on their potential to simulate Canadian watersheds. They implemented multi-criteria analysis for quantitative evaluation of HWQ models. The multi-criteria analysis used by Booty and Benoy [103] considers 21 model characteristics. Based on the multi-criteria quantitative evaluation, Booty and Benoy [103] identified five potential models for successful application for Canadian watersheds, and further discussion in this review will concentrate only on those five models. They included the Soil and Water Assessment Tool (SWAT) model [104], Annualized Agricultural Non-Point Source (AnnAGNPS) model [105], Better Assessment Science Integrating Point and Nonpoint Sources (BASINS) (U.S. Environmental Protection Agency [106]), Gestion Intégrée des Bassins versants à l'aide d'un Système Informatisé (GIBSI) [107], and Agricultural Non-Point Source pollution (AGNPS) model [108].
Applications of the selected five HWQ models to determine CSAs are listed in Table 1. It appears that the SWAT model is the most widely applied model followed by the AnnAGNPS and AGNPS models to identify CSAs as compared to other detailed HWQ models. As noted in the previous section, the wider applicability of SWAT for CSA identification is attributed to advantages of the SWAT model over other HWQ models. Despite the relatively wider applicability of the SWAT model as compared to other HWQ models, its application is still less popular as compared to the PI method.
The review of literature also shows that most detailed HWQ models were applied within the U.S. In contrast, application of detailed HWQ models to determine CSAs of P loss is very limited for Canadian watersheds [109]. Table 1. Overview of the top five hydrologic and water quality (HWQ) models applied to identify critical source areas (CSAs).

References
Model Used and Location of Study Area

Scale of Application and Modeling Approach
Arabi et al., [110] SWAT

Current State and Challenges of Determining CSAs
Substantial research efforts have been made to understand and develop techniques that can be used to identify CSAs of diffuse P loss since the mid-1990s [28]. However, identification of CSAs still remains challenging. Selection of appropriate model complexity level, scale issues, site-specific conditions representation (such as runoff generation mechanism, landscape depressions, tile drains, etc.), and modeling practices implemented prior to using models to identify CSAs are among the challenges that need to be addressed.

Level of Model Complexity and Scaling Issues
The issue of model complexity has been at the center of a debate in the hydrologic community, perhaps throughout all science [127]. These concerns become even more apparent when the issue moved to models of NPS pollution. In determining CSAs, the debate appears to be between the uses of simple index-based or detailed, complex HWQ models. However, there is no clear agreement among all the users regarding which method to use. One group argues to use of detailed HWQ models for identifying CSAs as they can provide the magnitude of P loss, and even may also provide the amount of different forms of P (such as soluble and particulate). More importantly, they are able to answer future management as well as climate change scenario analysis. In addition, development and application of a detailed HWQ model is becoming increasingly common with the advancement of computing resources over time [128]. On the other hand, the opposing group argues that detailed HWQ models are complex and require extensive input datasets, computational time, and manpower, unlike simple index-based methods. For a similar reason, detailed HWQ models are rarely used by conservation agents and farmers. They often implement simple index methods, such as the PI, as they are cost effective, user-friendly, and easy to apply [28,46].
Recently, there has been an increasing interest to use both methods in combination. Indeed, some researchers are trying to develop a 'Toolbox' that contains different methods (e.g., [28]) as it is difficult to answer all questions using a single method. As described by Sharpley et al. [28], such an approach has numerous merits as compared to the use of a single method. For instance, both simple indicators and detailed HWQ models could be included in the 'Toolbox'. This allows users to have an alternative for their problems. In such an approach, simple methods can be used in screening large area watersheds and detailed models can be used to answer complex questions, such as future management scenarios, like how long it takes after BMPs are implemented before meaningful improvements in water quality can be expected.
Scaling issues also are a challenge in CSA identification, as both source and transport factors vary with time so that P loss is variable over time. Variability of P loss with time leads to dynamic CSAs. CSAs are highly variable especially when they are considered over shorter time scales. However, static site assessment tools such as the PI are not able to account for the temporal variability of CSAs. They tend to make predictions on the long-term average annual temporal scale. In contrast, detailed HWQ models have the capability to simulate temporal variability of CSAs. However, these models are data intensive. In addition, spatial discretization is another issue in identifying CSAs. CSAs targeted BMPs are often implemented at the field level. However, semi-distributed models have been widely used to identify CSAs for a computational reason. The SWAT model, for example, is among the most widely applied semi-distributed models. In such semi-distributed models, a watershed is subdivided into smaller units such as the sub-basins and the heterogeneity of hydrological conditions within a sub-basin is handled using Hydrological Response Units (HRUs). The problem with such methods is that spatial location of HRUs within a sub-basin is not known. Thus, identification of CSAs at the HRU level is difficult. It is not easy to identify CSAs at the field scale level using semi-distributed HWQ models. To solve the problem, researchers suggest the development of an ArcGIS extension that is capable of tracking the location of HRUs. In this regard, researchers (e.g., [30,129]) have developed an ArcGIS tool to map SWAT outputs at the field level for conservation practice targeting.

Representation of Phosphorus Transformation and Transport Processes
As noted in the 'Introduction' section, the complexity of P in soils, sediment, water, and the phenomena involved in the transport of varying forms of P from the landscape to surface water bodies is complex and difficult to describe mathematically [52]. For instance in-stream (channel) processes, which include a combination of physical, chemical, and biological processes, regulate the downstream delivery of P [130]. Indeed, many studies (e.g., [131,132]) have shown that P entering rivers do not match with those measured at the watershed outlet. In-stream processes are influenced by several factors that include physicochemical controls, such as sorption/desorption reactions, mineral precipitation and dissolution, advection and diffusion and biological controls, such as periphyton and phytoplankton, microorganisms, and macrophytes [130]. However, most HWQ models often do not fully include the role of these factors in the algorithms used for in-stream retention and cycling processes. For example, most conventional models SWAT, Water Quality Analysis Simulation Program-5 (WASP5), and MIKE 11 focus on biological processes in which nutrient decomposition is considered as a first-order decay process without explicitly accounting for the role of the hyporheic zone processes [133]. In addition, the SWAT model does not consider mineral precipitation and dissolution of P in the river system [134].
The other issue in identifying diffuse P loss is the legacy P. Numerous studies have reported that legacy P significantly affects P dynamics in a watershed system. As described by Sharpley et al. [135], a portion of P can be accumulated at various locations along transport pathways from the landscape to the watershed outlet. Soils, downslope areas, ditches, streams, rivers, lakes, wetlands, and riparian areas are some of the important areas of P accumulation. Legacy P refers to a continuing supply of P to downstream water bodies from remobilization of the accumulated P [51]. Sharpley et al. [135] categorized legacy P into three groups: terrestrial, river, and standing water P legacies. The effects of land and nutrient management practices on the buildup of soil P to levels in excess of the crop needs are considered under the terrestrial P legacies. Though the legacies of P stored in land, river, and lake systems are clearly important issues, little has been done to include these processes into currently available methods of identifying CSAs.
Substantial research efforts have been made to understand and model P transport and chemical transformation processes [51]. There are three important stages of development that water quality models passed through since Phelps and Streeter [136] built the first water quality model. These include: an early stage , an improvement stage , and the broadening stage (after 1995) [137]. In these stages, HWQ models have made progress from point source models to the integrated models of point and NPS model, and from single variable water quality to multi-variable water quality, steady-state to the dynamic models, and one-dimensional to two and three-dimensional models [137]. The development of modeling theory and computing resources has contributed significantly to this progress in water quality modeling. Currently available HWQ models may represent the P cycle in soils and water in any of several ways, differing in the number of components (stores) and P transformation processes included [138]. The SWAT model, for instance, represents the P cycle in the soil using six different pools [134].
Some processes of P, however, are not yet fully integrated or are formulated in a simplified way in most HWQ models. Dissolved P transport in overland flow, for instance, is simulated in many HWQ models such as AnnAGNPS, SWAT, among others by assuming a single value of an extractable coefficient for the whole watershed to estimate surface runoff dissolved P [134]. In reality, however, the soluble P-extractable coefficient varies with soil type and management (e.g., [139]). Soil texture, aggregate diffusion, the degree of interaction between soil and water, organic matter content, vegetative soil cover, and P sorption capacities are some of the factors that influence P release [28]. Furthermore, the SWAT does not model dissolved P transport in subsurface flow. However, a concentration of soluble P in the shallow aquifer and groundwater flow can be specified to account for loadings of phosphorous from groundwater. The subsurface dissolved P transport algorithm also keeps the concentration constant throughout the simulation period, unlike in reality.

Landscape Connectivity to Stream System Accounted for
The link or connectivity between the landscape and the stream is an important factor in nutrient export and/or CSAs at the watershed scale [140]. In reality, all areas that generate nutrients are not necessarily contributing to the stream system. Furthermore, hydrologic connectivity is heterogeneous in space and often temporally transient [140]. However, connectivity between the source areas and the stream is not properly simulated by most existing HWQ models. The SWAT model, which is the most widely applied HWQ model [128], for instance, does not simulate landscape connectivity between nutrient generating areas to the stream. The SWAT model does not have a routine to route generated nutrients from landscape areas to the stream, but it assumes that all the nutrients generating areas would contribute to the stream or the watershed outlet. In SWAT, nutrients generated at any point within the sub-basin are treated identically regardless of position [128]. Such an assumption may not be valid for areas that have no pathway to the streams. Though the HRU concept used by the SWAT model reduces excessive computational requirements, the assumption about HRU connectivity may not be valid over some watersheds. Similarly, the use of the grid cell concept (such as that in the AnnAGNPS model) does not account for the position of the grid cell relative to the nearest water body in terms of connectivity between nutrient generation areas to the stream [140].

Modeling Practices
Accurate identification of CSAs depends on the reliability of the HWQ model used to identify CSAs. Therefore, HWQ models need to pass through good modeling practices prior to their application to identify CSAs. Nowadays, it is increasingly becoming common to observe misuse of HWQ models. Watershed modeling textbooks (e.g., [141]) present a scheme of good modeling practices. The suggested scheme contains sequences of activities that include sensitivity analysis, calibration, validation, and uncertainty analysis. According to Moriasi et al. [142,143], sensitivity analysis is the process of estimating the rate of change in output with respect to changes in model input parameters. The process is important to identify key parameters required for model calibration. It reduces the number of parameters that should be considered for the model calibration (parameter estimation stage). Otherwise, calibrating all the parameters of a detailed HWQ models could be time consuming. In the model calibration stage, the values of the selected model parameters, within acceptable ranges, Agriculture 2020, 10, 468 9 of 25 are estimated. This process involves identifying the 'optimal' or 'near optimal' acceptable values of the selected parameters. There are two types of calibration approaches: manual (heuristic) and automatic calibration. Automatic calibration, which involves the use of a search algorithm, is fast and less subjective as compared to manual calibration. It is particularly important for distributed models. Otherwise, it would have been time consuming to use manual calibration.
As model calibration does not guarantee the reliability of model predictions for other periods [144]; therefore, the calibrated model needs to be validated. During model validation stage, the model is run using the optimized parameters obtained during the calibration process. The other important component of a good modeling practice scheme is uncertainty analysis. Modeling uncertainty arises from multiple sources, which include input data used, parameter estimates, model structure, and even observed data used for calibration and validation [145]. Traditional modeling practices often involve fitting a single best simulation. However, there is no single best model as there is always an inevitable 'equifinality' problem [146]. It is, therefore, important to include uncertainty analysis as part of good modeling practices.
With better access to computational resources, a lot of research effort has been devoted to developing sensitivity analysis, optimization, and uncertainty analysis techniques in the watershed modeling literature. However, modelers often are not consistently following good modeling practices. A review by Wellen et al. [109] has shown that only 25% of water quality modeling studies reported results of sensitivity analysis. The number is even much lower for quantifying modeling uncertainty. Only 10% of water quality modeling studies tried to account for the uncertainty of the model predictions in a quantitative way. Similarly, only 17% of studies reported the use of any kind of optimization technique.
Researchers have recommended to use both graphical and statistical evaluation procedures while evaluating model performance during calibration and validation periods [143]. Indeed, model evaluation guidelines that consider the recommended model evaluation statistics with corresponding performance ratings and appropriate graphical analyses were established (e.g., [142,143]). According to Wellen et al. [109], 77% of water quality modeling studies quantified the goodness of fit of flow and water quality variables. However, such practices are rarely seen in estimating CSAs. Due to field data limitation, studies of water quality model application to identify CSAs often do not perform field verifications at all or not in a rigorous way.

Canadian Hydrological Conditions Conceptualization in Models
Numerous HWQ models have been developed, but most of them are not primarily designed for Canadian hydrological conditions. Canada exhibits highly diversified hydrology due to its landscape diversity that includes the Arctic tundra of the north, great prairies of the central area, the Rocky Mountains in the west to the Great Lakes, the St Lawrence River, and Niagara Falls in the southeast of Canada (Figure 1). The occurrence of infiltration and saturation excess runoff generation mechanisms, cold-climate hydrology, landscape depressions, and tile drainage system are among the diversified hydrologic conditions present in Canadian watersheds. The influence of each process is variable in time and space across Canada [23].
landscape diversity that includes the Arctic tundra of the north, great prairies of the central area, the Rocky Mountains in the west to the Great Lakes, the St Lawrence River, and Niagara Falls in the southeast of Canada (Figure 1). The occurrence of infiltration and saturation excess runoff generation mechanisms, cold-climate hydrology, landscape depressions, and tile drainage system are among the diversified hydrologic conditions present in Canadian watersheds. The influence of each process is variable in time and space across Canada [23].

Runoff Generation Mechanism
Saturation excess and infiltration excess runoff generation are the two types of generally accepted runoff generation mechanisms. The former, saturation excess, occurs when soil becomes saturated, and any additional precipitation causes runoff [148]. The concept of Variable Source Area (VSA) hydrology often is used to explain this type of runoff generation mechanism [149,150]. This type of runoff generation generally occurs in humid and thickly vegetated regions with permeable shallow soils underlain by an impervious layer [37,151]. The later, infiltration excess also called Hortonian overland flow, occurs when the net precipitation exceeds the infiltration capacity of the soil [40,152]. The infiltration rate depends on soil properties, land use, and landscape conditions [153][154][155] and on the rainfall intensity.
However, many HWQ models that have potential for use to identify CSAs, often are designed for infiltration the excess runoff generation mechanism [156] such as the AGNPS [108], Chemicals, Runoff, and Erosion from Agricultural Management Systems (CREAMS) [157], SWAT [104], and the GWLF [88]. These models use the USDA SCS CN equation [158] to predict runoff. In reality, both

Runoff Generation Mechanism
Saturation excess and infiltration excess runoff generation are the two types of generally accepted runoff generation mechanisms. The former, saturation excess, occurs when soil becomes saturated, and any additional precipitation causes runoff [148]. The concept of Variable Source Area (VSA) hydrology often is used to explain this type of runoff generation mechanism [149,150]. This type of runoff generation generally occurs in humid and thickly vegetated regions with permeable shallow soils underlain by an impervious layer [37,151]. The later, infiltration excess also called Hortonian overland flow, occurs when the net precipitation exceeds the infiltration capacity of the soil [40,152]. The infiltration rate depends on soil properties, land use, and landscape conditions [153][154][155] and on the rainfall intensity.
However, many HWQ models that have potential for use to identify CSAs, often are designed for infiltration the excess runoff generation mechanism [156] such as the AGNPS [108], Chemicals, Runoff, and Erosion from Agricultural Management Systems (CREAMS) [157], SWAT [104], and the GWLF [88]. These models use the USDA SCS CN equation [158] to predict runoff. In reality, both infiltration excess and saturation excess runoff generation mechanisms contribute to overland flow. However, often only one of the processes dominates in a specific time and location. The SCS CN should only be used 'if it is known a priori that the runoff is produced by the infiltration excess mechanism [156]. However, Steenhuis et al. [159] suggested that the theoretical basis of the CN method is valid for both infiltration and saturation excess runoff generating mechanisms [160].
In Canada, both infiltration and saturation excess runoff generation mechanisms occur; though the dominant process varies between seasons. In summer when the soils are dry, the infiltration excess runoff generation mechanism is dominant. As such, the original CN method can be used to simulate runoff generation. However, using the curve number method will have some serious implications [156] while trying to simulate runoff generation in areas where the saturation excess mechanism prevails, such as the case in winter and early spring months in Canada as soils are nearly saturated. Therefore, most of the currently available HWQ models, which use the CN method, need to be modified to accommodate the saturation excess runoff generation mechanism [90].
In general, appropriate representation of the runoff generation mechanism is important for accurate estimation of NPS pollution export to the water system. In this regard, there have been some research efforts to upgrade the existing HWQ models to include the saturation excess runoff generation process. Winchell et al. [50], for instance, upgraded the SWAT model to identify CSAs in Lake Champlain's Missisquoi Bay watershed, which is dominated by saturation excess runoff generation processes. In order to properly represent the runoff generation process in the region, Winchell et al. [50] modified the SCS curve number method of SWAT. In their work, they developed and applied an approach for adjusting curve numbers based on the local compound topographic index. Similarly, Mishra and Singh [161] recommended the use of SCS CN method in cases when the potential maximum retention is less than or equal to twice the amount of precipitation. Withstanding the need to include the saturation excess mechanism in rural, humid regions, Easton et al. [90] re-conceptualized the SWAT model. As such, they modified the way the CN and available soil water content are defined in the SWAT model.

Landscape Depressions
Landscape depressions are important features that have an influence on both runoff and pollutant generation and transport across a watershed. They are more important in watersheds that fall within the Canadian Prairie Pothole Region (PPR), where millions of landscape depressions exist. In this region, the runoff and pollutant generation and transport occur under the existence of numerous landscape depressions varying in storage capacity and having dynamic connectivity to one another. This leads to a dynamic area contributing to streamflow and pollutant export, which invalidates the application of conventional HWQ models requiring a basin's contributing area to be a fixed value [162].
For computational reasons, semi-distributed HWQ models have wider applicability as compared to fully distributed models. However, semi-distributed HWQ models such as SWAT allow only a single depression storage per sub-watershed. In reality, numerous landscape depressions exist within a sub-basin, which invalidates the applicability of semi-distributed models in the Canadian PPR. To overcome this problem, several researchers aggregated all the depressions within a sub-basin and represented them by a single synthetic storage. Kiesel at al. [163] and Almendinger et al. [164], for instance, applied the SWAT model in depression-dominated watersheds by aggregating depressions and representing them by a single storage. Such an approach, however, does not represent the real processes as landscape depression response depends on depression's storage capacity and connectivity. Alternatively, Mekonnen et al. [165] have proposed and applied new modeling techniques to handle numerous landscape depressions within a watershed (e.g., [166]). In the upgraded SWAT model, the numerous landscape depressions within a sub-basin are incorporated using a probability distribution function. The modified SWAT model, also called SWAT with Probability Distributed Landscape Depressions (SWAT-PDLD), was tested over two Canadian prairie watersheds and an improved streamflow simulation was reported [166].
In addition, several researchers have also confirmed the existence of significant chemical transformations within landscape depressions. In a Canadian prairie watershed, for instance, the transformation of pollutants is expected to be significant in potholes [167]. HWQ models, however, often neglected such transformation processes. The SWAT model, for instance, assumes, as there is a transformation within the storage modules called 'Pothole', 'Pond', and 'Wetland'. For a better simulation of P processes in depression-dominated watersheds, existing HWQ models need to be improved to include transformation processes in the landscape depressions.

Cold-Climate Conditions
Most Canadian watersheds exhibit cold-climate hydrology that includes frozen ground, snowpack accumulation, melting, and freeze-thaw cycles. The cold-climate condition has a significant influence on streamflow and water quality [31]. Surface runoff from frozen or partially frozen soils increases due to reduced infiltration [168,169]. Spring runoff is an important component of streamflow in Canadian watersheds. The influence is not limited to streamflow, but also to water quality variables. In cold-climate watersheds, sediment export is expected to be high during the spring [7,170,171]. In Canada, more than 75% of the stream flow and sediment loads occur during late winter and early spring [31,171]. This is because of three main reasons: (a) increased surface runoff enhanced by reduced infiltration in frozen or partially frozen soils [168,169]; (b) an increased soil erodibility during freeze-thaw cycles [171][172][173]; and (c) the longer duration of the snowmelt-runoff period as compared to individual rainfall-runoff events [174]. Indeed, several researchers have revealed reduced erosion resistance of frozen soils both in laboratory (e.g., [171,175,176]), and field conditions (e.g., [7,172]). The erosion resistance of the soil is thought to be reduced by the freeze-thaw cycles as the ice that forms in soil voids during freezing pushes the soil grains apart [175,176]. The influence of cold-climate conditions is not only limited to the soil erodibility factor, but also to other variables such as hydraulic conductivity of soil as well as to the source, transformation, and P loss as the latter processes are disturbed by soil frost, freeze-thaw cycles, and snowpack accumulation and melting [5,173].
Cold-climate conditions are missing or not yet well incorporated in most existing watershed models [5]. The original SWAT model, for instance, uses a constant annual value of soil erodibility throughout the simulation period. To take into account seasonality of soil erodibility, a recent upgrading of the sediment module of the SWAT model was made by Mekonnen et al. [165] to incorporate cold-climate conditions into the sediment module of SWAT. The improved model was tested for two Canadian watersheds for sediment and nutrient export, and improved simulations are reported by Mekonnen et al. [53,165]. The modification considers a variation of the weighting factor of seasonal soil erodibility between soil types. An experimental study, however, shows that the weighting factor varies between seasons [173]. Therefore, the improved upgraded sediment module needs to be further reformulated to incorporate variations of the weighting factor with soil type.

Tile Drains
In some parts of Canada, indeed in some parts of North America, tile drainage is a common practice in agricultural areas with high water tables [177]. The USDA [178] claims that settlement in North America is land drainage facilitated. Currently, about 62 million ha of agricultural areas are tile drained over the North American landscape [179]. In Ontario, about 43% of the total 3.5 million ha of land classified as cropland (OMAFRA, [180]), have been tile drained. Such practice changes the hydrology of the watershed in several ways as compared to natural, undrained lands. For instance, tile drains can change the way P is exported from the landscape into the streams [181]. Studies have shown that tile drains in Canada are significant contributors of P export from agricultural land [181][182][183][184][185][186][187].
For accurate simulation of areas dominated by tile drains, it is important to understand all processes occurring in tile drains system. However, understanding tile drains process is complex as tile outflow consists of water that has drained quickly through macropores (also called preferential flow) and water that has percolated through the soil matrix to the drain [188]. Furthermore, often tile drain systems were built so long ago that the actual details of these systems are unknown making it difficult to model these systems in detail. The differences in transformation and transport processes of P between the preferential flow and soil matrix pathways are significant. In a preferential flow pathway, which constitutes channels of earthworm burrows, root channels, and others, the interaction between the channel walls and P losses being transported through them is relatively minimal. Due to the minimal interaction, researchers often assume similar proportions of particulate and dissolved P in tile flow as in surface runoff (e.g., [188]). Field data support similar findings (e.g., [181]). In contrast, P has plenty of time to interact when it moves through the soil matrix often is precipitated out of solution [188].
In the glacial soils of Canada and the northern U.S. states, the movement of P through the soil matrix to tile drains often is assumed negligible due to high levels of calcium, iron, and aluminum available to react with and precipitate P out of solution. The P portion that reaches the tile drains continue its movement through the tile drains, which represent a direct conduit from the field to the outlet. Thus, P that reaches the tile outlet can be carried from a larger part of the landscape than would otherwise be possible. Tile drainage system design factors, such as tile spacing and depth, also influence tile drain outflows. Furthermore, for Canadian conditions, the movement of P losses through the different pathways need to be adjusted seasonally to account for the effect of many factors including antecedent moisture, snow accumulation and melting, and freeze-thaw cycles. Importantly, macropore flow also is highly affected by these factors. Indeed, they can range from an increase in macropores from freeze-thaw actions to a complete blockage of macropores with ice that blocks this pathway until after spring thaw.
Generally, for the application of different methods to determine CSAs, the proportion of P delivery through tile drain system is either neglected or incorporated in a simplified way without explicit representation of the process. For example, a simplified concept of the P movement mechanism through the tile drains, as a component of the PI method, was introduced by Reid [59]. Similar attempts also are made to incorporate tile drains into detailed HWQ models. For example, Arnold et al. [104] enhanced the capability of the SWAT model by introducing a subsurface drainage component. The modification, however, does not allow change in tile spacing and tile size. To address this problem, further improvements were made by Moriasi et al. [142] by incorporating the steady-state Hooghoudt's [189] and Kirkham [190] tile drain equations into the SWAT model. The improved tile flow simulation approach considers tile spacing and drain tube size as well as the depth of tile drain. Golmohammadi et al. [191] explicitly incorporated tile drainage into the SWAT model by replacing the subsurface hydrology of SWAT with the subsurface hydrology of DRAINMOD. Several other researchers have also reported an improved simulation of tile drain dominated watersheds using an improved tile drainage module in SWAT (e.g., [192]), including a module that can explicitly simulate dissolved reactive phosphorus through tile drains [193]. These improved methods still need a thorough evaluation to test their ability in simulating P loss and identifying CSAs in Canadian conditions.
In conjunction with tile drainage, Water and Sediment Control Basins (WASCoBs) are becoming increasingly popular in Canadian agricultural watersheds [194,195]. A WASCoB generally consists of a berm, surface inlets, and a drainage system which diverts surface runoff ponded behind the berm in underground tiles and traps NPS pollution within the berm [196,197]. A WASCoB is generally placed in depressions in order to impede gradual formation of ephemeral or permanent gullies [198].
Withstanding the drainage mechanism of a WASCoB, the upstream areas draining to the WASCoB, previously non-contributing, are, thus, becoming contributing area for runoff. Hence, such widespread use of WASCoBs in Canadian agricultural watersheds is profoundly altering hydrology of these watersheds. Most of them are not designed to represent complex hydraulic interactions, and feasible integration with hydraulic models is, therefore, sought [199]. Proper representation of such complex interactions occurring between the berm, surface inlets, and drainage system of a WASCoB is a challenging task. As such, modeling of a WASCoB system in popular HWQ models, such as the SWAT, has rarely been practiced. Her et al. [196] simulated WASCoBs in a similar way to the terraces on the contour by adjusting related SWAT parameters such as CN, P-factor of the USLE, average slope length of the HRU, etc. Similarly, Yang et al. [200] developed a new SWAT module to explicitly represent hydraulic behavior of a WASCoB. The module is similar to a reservoir implementation in SWAT, and, thus, is a simplification of reality. A more robust representation of hydraulic as well as the NPS pollution entrapment behavior of WASCoBs is warranted for more accurate simulation of P loss to identify CSAs in Canadian conditions.

Future Research Directions
Encouraging progress has been made to identify CSAs of diffuse P where runoff and high P sources coincide. This was also the topic of a keynote paper at the Fifth International Phosphorus Workshop (IPW5) [201]. However, numerous challenges remain to identify CSAs under Canadian conditions. The following future research directions are recommended to further improve currently available methods to identify CSAs for Canadian conditions: • A better database is needed to improve the model prediction capability. More importantly, the TI method will benefit from better resolution data, such as Light Detection and Ranging (LiDAR) data, which could also be helpful to incorporate the saturation excess runoff generation into HWQ models.

•
Most BMPs are implemented at the field scale. Therefore, it is important to improve field or sub-field scale modeling capability of the existing methods. This is particularly important while using semi-distributed HWQ models such as SWAT. This can be achieved by developing an ArcGIS extension that can identify an HRU location within a sub-basin. This will also facilitate evaluation of field and sub-field level output.

•
There also is a need for improvements in the predictions of dissolved P in surface runoff from soil, which is now typically based on simple extraction coefficients relating a measure of soil P (usually routine soil tests) to dissolved P concentrations in runoff. Most existing methods use a constant value for all soil types. It is, therefore, necessary to conduct research to identify the variations in extraction coefficients with soil properties such as soil texture.

•
With significant variations at field and watershed scale estimation of P loss, it is important to further investigate if different algorithms for sediment and P loss should be used.

•
Canadian hydrology exhibits a distinct seasonal pattern. It is, therefore, necessary to quantify and include the impact of seasonal hydrology of P generation and transport processes in HWQ modeling approaches.

•
With better understanding of physical, chemical, and biological processes affecting water quality, future research should also focus on improving the representation of processes, such as P stratification and legacy sources in HWQ models. Both vertical stratifications of P in no-till soils and legacy sources of P need to be incorporated as they affect the short-term benefits of remedial practices applied to CSAs.

•
Another research area that needs attention is nutrient processes in the wetland modules. For example, the 'Wetland' and 'Pond' modules of the SWAT model do not consider nutrient transformation. However, several field studies have indicated that transformation of nutrients indeed occurs within the ponds, wetlands, or potholes. • Development of a 'Toolbox' to identify CSAs also needs serious attention and could be a step forward. In fact, Sharpley et al. [28] envisioned the need and value of such a 'Toolbox'. Especially given the fact that it is very difficult to answer all pertaining questions related to identifying CSAs of NPS pollution using a single model/method. The envisioned 'Toolbox' can host a range of different models/methods, ranging from a simple method such as topographic index to a complex HWQ model. This will provide the end users with a variety of options to try to use different methods/models to identify CSAs. Users may opt to use the simple methods for screening larger watersheds while detailed HWQ models may be used to provide absolute values of P loss and long-term effects of future management scenarios on P loss. The 'Toolbox' can be hosted in a dedicated platform or can be made available as a web-based modeling framework such as the Hydrology and Water Quality System (HAWQS) using SWAT framework [202]. The 'Toolbox' essentially offers features such as an automated workflow of input data preparing and processing for an area of user's choice and output data repository to store model and scenario runs. We believe that provision of such a 'Toolbox' renders the redundancies associated with pre-processing of large volume of input data (e.g., spatial and meteorological) which is often time-intensive and error-prone. Recent advances in computing facilities and web technologies also support the idea of making the envisioned 'Toolbox' available to the wider public.

•
With advances in computational resources and techniques, future research on CSAs should focus on the sources of uncertainty, uncertainty quantification methods, and incorporation of all possible uncertainties in HWQ modeling approaches.

•
With advances in instrumentation technology, future research should also focus on field research to identify CSAs. Special attention is needed on the quantification on the seasonal variability of CSAs. Data collected from such experiments will be useful to evaluate the CSA prediction capability of HWQ models.

•
It has been shown that the majority of NPS pollution would be exported in the late winter and early spring months [31,171,203] and, therefore, constitute hot-moments of NPS pollution export. Sampling campaigns should be effective enough to capture the variability in such important time/seasons. Identification of the proper locations for water quality monitoring is equally important too. Thus, a cost-effective and efficient water quality monitoring framework, as suggested by Alilou et al. [204] is desired. Moreover, detailed HWQ model ideally needs continuous water quality data for proper calibration (and/or validation). However, such data are rarely available because of logistic issues. Often, sporadic water quality measurements are available which hinders proper model calibration (and/or validation) and model results may deem to be unreliable. This is a general problem in NPS pollution modeling. One possible solution is to carryout dedicated sampling campaigns to cover wide range of streamflow/field conditions (e.g., dry, wet, 24 h, etc.) which helps to verify a model's ability in respective conditions [205]. However, this requires a team of dedicated crew waiting for such an event to occur, which may not always be available [206]. Another possibility is to install automatic samplers to measure some explanatory variables (e.g., specific conductance, turbidity) and using an established relationship (e.g., regression model), other water quality variables (e.g., total suspended solids, total phosphorus) are predicted [207]. However, uncertainties associate with such estimations render their use in calibrating (and/or validating) detailed HWQ models. The use of remote sensing technologies in validating HWQ model based CSAs can also be an alternative as shown by Shrestha et al. (in press). In the study, they used oblique aerial images to qualitatively validate CSA of phosphorus in a watershed. Such a technique may be useful for a large watershed, which needs water quality monitoring at multiple locations.

Conclusions
Currently, different types of methods, ranging from simple indicators to detailed HWQ models, are available to identify CSAs of NPS pollution. Selecting the right method has been a challenge for many years. In general, it appears that the PI method is the most widely used method amongst the simple methods. On the other hand, the SWAT model is ranked first as a continuous simulation model for long-term studies, whereas AGNPS is used as an event-based model for more detailed studies of specific events.
As no single method/model is perfect, there has been a growing tendency towards the development of a 'Toolbox' that could host a variety of methods/models to identify CSAs. Such an approach also facilitates the use of simple index-based methods in conjunction with detailed HWQ models. Results from detailed methods will also help to evaluate and improve the accuracy of simple methods. Otherwise, it is difficult to verify the results, as field monitoring of CSAs is almost non-existent. Furthermore, such a 'Toolbox' allows users to compare the performance of different methods. However, a direct comparison between different methods often is challenging due to differences in the scale of applications.
For Canadian watersheds, applications of both simple index-based and detailed HWQ models are limited. There are, however, some notable recent efforts to determine CSAs in Canadian watersheds. Most of these studies focus on modifying and applying PI and TI methods. Despite their limited applications in the past, HWQ models such as SWAT, AnnAGNPS, BASINS, GIBSI, and AGNPS have potential to use in Canadian watersheds. Direct application of these HWQ models for Canadian watersheds has proven to be a limited success, as Canada's unique hydrological conditions are not fully incorporated in these models. As a result, currently available models should be adapted to represent specific conditions of the region prior to their application. Some of the processes that should be incorporated include the cold-climate conditions, landscape depressions, tile drains, saturation, and infiltration excess runoff generation mechanisms, and incorporation of more detailed processes in phosphorus cycle, such as source legacy, among others. Furthermore, a good modeling practice, including quantification of different sources of uncertainties, should be followed.