Threshold storm approach for locating phosphorus problem areas: An application in three agricultural watersheds in the Canadian Lake Erie basin

Journal of Great Lakes Research xxx (xxxx) xxx

Contents lists available at ScienceDirect

Journal of Great Lakes Research journal homepage: www.elsevier.com/locate/ijglr

Threshold storm approach for locating phosphorus problem areas: An application in three agricultural watersheds in the Canadian Lake Erie basin B. Zhang a, N.K. Shrestha a,⇑, R. Rudra a, R. Shukla a, P. Daggupati a, P.K. Goel b, W.T. Dickinson a, N. Allataifeh a a b

School of Engineering, University of Guelph, Guelph, Ontario, Canada Ministry of the Environment, Conservation and Parks, Etobicoke, Ontario, Canada

a r t i c l e

i n f o

Article history: Received 25 March 2019 Accepted 13 November 2019 Available online xxxx Communicated by Craig Stow

Keywords: Threshold storm AGNPS Critical source area Phosphorus

a b s t r a c t Hot-spots and hot-moments of phosphorus loads in an agricultural watershed depend not only on the watershed characteristics but also on the type and intensity of storms. Not all storms will generate phosphorus that can be considered problematic. A threshold storm is thus proposed and defined as the maximum storm intensity in which the phosphorus generated in a watershed is below seasonal phosphorus tolerance limit. To evaluate the threshold storm approach, separate Agricultural Non-point Source (AGNPS) models for three diverse small agricultural watersheds in southern Ontario, Canada were calibrated for runoff volume, sediment yield, and total phosphorus and run for representative storms with increasing return periods (2-year through 100-year). Results showed that in an upland watershed (Holtby), a 4.8-year early spring storm tend to generate phosphorus load above the threshold limit for the season. The same for low-land watersheds (Wigle and Jeannette) were, respectively, 14.9-year and 12.4-year. In all three watersheds, summer storms up to 100-year will fail to reach the seasonal tolerance limit for phosphorus. The critical source areas, identified based on the threshold storms, were distributed uniformly across the watersheds. As a phosphorus problem is essentially a source problem, such a simple yet robust approach to identify critical source areas of phosphorus can be useful in designing costeffective best management practices. Crown Copyright Ó 2019 Published by Elsevier B.V. on behalf of International Association for Great Lakes Research. All rights reserved.

Introduction Increased phosphorus loss from agricultural areas has led to water quality impairments of receiving water bodies (Paerl et al., 2011). In Ontario, Canada, an increase in phosphorus loss from agricultural areas is attributed to extensive fertilizer use and has posed serious environmental threats such as the blue-green algal blooms in Lake Erie (McQueen and Lean, 1987; Michalak et al., 2013; Scavia et al., 2014). In order to lessen further degradation, various agronomic best management practices (BMPs) have been shown to be effective (Liu et al., 2016; Merriman et al., 2019; Merriman et al., 2018a,b; Rodriguez et al., 2011; Shrestha and Wang, 2019; Stang et al., 2016; Yang et al., 2013). Due to lack of a targeted approach, BMPs are often implemented over the entire watershed (Daggupati et al., 2011) while only certain areas in a watershed generate the majority of phosphorus load. These areas ⇑ Corresponding author. E-mail address: [email protected] (N.K. Shrestha).

are often termed as critical source areas (CSAs) (Sharpley et al., 2011). Proper identification of CSAs for phosphorus is thus crucial for economic efficacy of BMPs (Meals et al., 2012; White et al., 2009). Watershed models are increasingly being used for identification of CSAs (Winchell et al., 2015). Choice of a suitable model is also crucial as it can affect the locations of CSAs (Niraula et al., 2013). In general, a watershed model is calibrated and validated against observations. Model outputs are then used to create spatial maps of phosphorus loss. Based on arbitrarily selected classes of phosphorus loss, CSAs are identified as those areas having the highest phosphorus loss. For example, in the Grand River Basin, Ontario, Liu et al. (2016) divided annual phosphorus yield into five classes with breaks at 0.45, 0.63, 0.84, 1.1 and 1.8 kg P/ha. For the same river basin, Hanief and Laursen (2019) also used five classes, however, phosphorus yield ranged from 0.1 to 0.5 with an increment of 0.1 kg P/ha/yr. Merriman et al. (2019) and Merriman et al. (2018a) opted for seven classes (0.5 to >3.0 with an increment of 0.5 kg P/ha/yr) in Upper East River, Wisconsin and Eagle Creek

https://doi.org/10.1016/j.jglr.2019.12.003 0380-1330/Crown Copyright Ó 2019 Published by Elsevier B.V. on behalf of International Association for Great Lakes Research. All rights reserved.

Please cite this article as: B. Zhang, N. K. Shrestha, R. Rudra et al., Threshold storm approach for locating phosphorus problem areas: An application in three agricultural watersheds in the Canadian Lake Erie basin, Journal of Great Lakes Research, https://doi.org/10.1016/j.jglr.2019.12.003

2

B. Zhang et al. / Journal of Great Lakes Research xxx (xxxx) xxx

Watershed, Ohio, respectively. In another study, Merriman et al. (2018b) divided phosphorus yield into five classes (0.3 to 1.5 with an increment of 0.3 kg P/ha/yr). Similarly, Winchell et al. (2015) partitioned annual phosphorus yield of Missisquio Bay Watershed, Vermont, into four classes (0.5 to 2.0 with an increment of 0.5 k gP/ ha). These example cases clearly illustrate non-uniformity or randomness in assigning severity classes of phosphorus yield. Devoid of a well-defined phosphorus tolerance limit, such an approach is merely useful in identifying CSAs of phosphorus. Therefore, there is a need to standardize the phosphorus severity class based CSA identification. A tolerance limit based approach will thus be useful. A study by Schertz (1983) provided a basis as to how such a tolerance limit for sediment can be defined. The study underlined the need to consider important factors such as soil formation rate, erosion’s effect on crop productivity, etc. while determining a soil loss tolerance limit. The study also stressed the need to make assumptions based on the knowledge gained through observation when ‘research was lacking or nonexistent’. Following this, Dickinson et al. (1986) suggested an annual sediment tolerance limit of 1.25 t/ha for agricultural watersheds in Southern Ontario. Owing to seasonal variability of sediment load (Wall et al., 1988), the tolerance limit should also be defined for different seasons. As such, identification of CSAs should be carried out in a seasonal basis (Shrestha et al., 2019). To our knowledge, seasonal tolerance limits of phosphorus have not been reported to date. However, due to high association between sediment and phosphorus load in agricultural watersheds (Wall et al., 1996), seasonal tolerance limits for phosphorus can be estimated. While accounting for seasonality in identifying CSAs, it should be noted that a storm of the same intensity but occurring in different seasons will have different potentials of pollutant delivery. Based on a seasonal phosphorus tolerance limit, a storm for which phosphorus loss from a watershed will exceed the limit can be determined. Such a storm can be termed as a threshold storm and the approach can be termed as the threshold storm approach. Testing the applicability of such an approach requires a model which can simulate the hydrology and water quality dynamics of selected events. Hence, an event based model, the Agriculture Non-Point Source (AGNPS) (Young et al., 1989) is selected as the model which has been widely tested in agricultural watersheds across the globe (Grunwald and Norton, 1999; Haregeweyn and Yohannes, 2003; León et al., 2004; Liu et al., 2008; Ma and Bartholic, 2003; Mohammed et al., 2004; Otieno and Onyando, 2012; Shrestha et al., 2019; Sun et al., 2009). In this study, the threshold storm approach is tested in three diverse small agricultural watersheds (Wigle, Jeannette and Holtby) of the Lake Erie basin located in southern Ontario, Canada. Three separate AGNPS models of the watersheds are built-up and calibrated against streamflow, sediment and phosphorus for selected storm events. Then, CSAs of phosphorus are identified based on representative storms for increasing return periods (2-year through 100-year). We believe that such an approach will help in standardizing the procedure in identifying of CSAs for phosphorus loss in small agricultural watersheds.

Materials and methods The threshold storm approach As highlighted above, annual or seasonal tolerance limits of phosphorus have not previously been suggested. However, based on the initial works of Schertz (1983), an annual sediment tolerance limit (1.25 t/ha) was suggested by Dickinson et al. (1986). The annual tolerance limit was used to derive seasonal sediment tolerance limits (Shrestha et al., 2019). It has been reported that

watershed scale sediment and phosphorus yield are highly correlated (Wall et al., 1996; McDowell et al., 2000; Reid et al., 2018). Hence, similar annual and seasonal tolerance limits for phosphorus can be estimated. Based on experimental works at various agricultural fields of southern Ontario watersheds, Wall et al. (1996) estimated a ratio between annual sediment (t) and phosphorus (kg) to be 1.23 ± 0.20 (mean ± 1.0 standard deviation). A non-negligible value of standard deviation indicates uncertainty/variability in the ratio estimation. This was indeed expected as the agricultural fields had different: (a) phosphorus levels; (b) forms of the phosphorus; (c) land-use and land-management practices (Wall et al., 1996). Using all the field measurements, including at both small and large watersheds, the study estimated the coefficient determination (R2) between annual sediment (t) and phosphorus (kg) to be just 0.62. Such a low value of R2 can affect reliability of the ratio estimation. Furthermore, Wall et al. (1996) conducted the experiments in mid-1990s. Moreover, the annual sediment tolerance limit was suggested in mid-1980’s (Dickinson et al., 1986). However, devoid of a more recent estimate, we had to use the sediment tolerance limit and the ratio. This is indeed one of the limitations of this study. These issues will need further investigations based on carefully designed field experiments. In order to account for the uncertainty in the ratio estimation, we calculated the upper limit and lower limit of the ratio using mean ± 1.0 standard deviation. These limits will encapsulate 68% confidence interval if the samples follow the normal distribution. We then calculated the annual phosphorus tolerance limit to be 1.55 ± 0.25 kg/ha (upper limit = 1.80 kg/ha and lower limit = 1.2 9 kg/ha). In order to calculate seasonal phosphorus tolerance limits, the annual tolerance limit was adjusted seasonally, weighted against the average seasonal storm erosion index (EI) to the annual EI for London, Ontario, Canada (Rudra et al., 1986). This resulted in seasonal phosphorus tolerance limits of 0.11 ± 0.02 kg/ha, 0.42 ± 0. 07 kg/ha, 0.74 ± 0.12 kg/ha and 0.28 ± 0.05 kg/ha, respectively, for winter, spring, summer and fall (autumn) seasons. These seasonal tolerance limits were then used to calculate their respective threshold storms. A threshold storm can be defined as the maximum storm intensity in which phosphorus generated in a watershed is below the seasonal phosphorus tolerance limit. Storms with different rainfall intensity and duration will occur in different seasons. In general, Southern Ontario experiences shorter duration summer, milder duration fall and long duration spring storms (Dickinson, 1976, 1977; Rudra et al., 2015). We have chosen representative storm durations for summer, fall and spring seasons to be 2-hour, 4hour and, 6-hour, respectively. Then, rainfall intensity for return periods of 2-year through 100-year are ascertained from monthly rainfall extreme statistics (mean and standard deviation) developed by Dickinson (1976), utilizing the double exponential Gumbel’s distribution. Two distinct periods in spring season; early and late spring are considered to account for changing rainfall pattern in spring season (Rudra et al., 2015). Fig. 1 shows representative seasonal storm intensities for different return periods (2-year through 100-year) used in this study. It should however be noted that the proposed threshold storm approach is one of a number of possible ways to derive seasonal thresholds. An alternative could be that thresholds were not the annual P or sediment loss divided into seasons but, instead, the ‘storm for each season’ for which the model predicts a certain amount of phosphorus/sediment loss. Moreover, the use of average seasonal storm EI to derive seasonal thresholds can be debated. The EI is only one metric of potential erosion, does not include a number of other factors (soil eroded, connectivity, proximity to the stream, vegetation along the flow path, runoff generation, etc.) that directly affect phosphorus/sediment loss. The proposed approach may be improved upon in the future.

Please cite this article as: B. Zhang, N. K. Shrestha, R. Rudra et al., Threshold storm approach for locating phosphorus problem areas: An application in three agricultural watersheds in the Canadian Lake Erie basin, Journal of Great Lakes Research, https://doi.org/10.1016/j.jglr.2019.12.003

B. Zhang et al. / Journal of Great Lakes Research xxx (xxxx) xxx

an area of about 374 ha. This watershed has extensive agricultural farms (86.5%). The rest of the land-use/land-cover include pasture, residential area, and woodland with coverage of 6.3%, 4.2%, and 3%, respectively. The majority of the watershed is drained by Brant silt loam and poorly-drained Muriel silt clay loam. The Holtby watershed is an upland watershed with the elevation range between 254 and 285 m (ESM Fig. S1) and slope up to 7% (ESM Fig. S2).

40 Early Spring (6-hour)

Rainfall intensity (mm/hr)

Late Spring (6-hour) Summer (2-hour)

30

3

Fall (4-hour)

20

AGNPS model build-up and calibration

10

0 1

10

1 00

Return Period (year) Fig. 1. Representative seasonal storms intensity for return periods 2-year through 100-year.

Agriculture Non-Point Source (AGNPS) model AGNPS is an event-based distributed hydrologic and water quality (H/WQ) model, primarily used for agriculture watersheds ranging from a few hectares to more than 20,000 ha (Young et al., 1989). The model requires spatial data sets such as digital elevation model, land-use and soil types, and precipitation data. The model discretizes a watershed into uniform user-defined grids/cells. The model then calculates various watershed properties such as slope, slope length, etc. as well as runoff, sediment and phosphorus related parameters. The model uses the soil conservation service curve number (SCS CN) method to calculate runoff generation. As such, the model utilizes different CN values as per soil hydrologic group and land cover. The peak flow is then computed using an empirical relationship developed by Smith and Williams (1980). As for the sediment, the model first computes soil erosion from each cell using the universal soil loss equation (USLE). The cell erosion rate is then compared with sediment transport capacity of the cell overland flow to determine the cell sediment yield. The sediment yield is routed to a downstream cell, using a steady-state continuity equation utilizing the Bagnold (1966) approach. Finally, for phosphorus, particular and soluble forms are estimated using empirical equations which are adapted from the Chemical, Runoff and Erosion from Agricultural Management Systems (CREAMS) model (US-SEA, 1980).

Study watersheds Three small agricultural watersheds; Wigle, Jeannette and Holtby (Fig. 2) are considered in this study. These watersheds drain to Lake Erie and represent three typical watersheds in southern Ontario, Canada. The first watershed, Wigle creek watershed has drainage area of about 1949 ha and has poorly drained Brookston clay (99.7%). Agricultural fields (84.3%) are the dominant landuse/land-cover type in the watershed (Fig. 2). Residential area, woodland, pasture, and water accounts for 8.2%, 4.6%, 2.7%, and 0.2% of the watershed area, respectively. The watershed has a low-lying flat topography (elevation range: 190 to 200 m, Electronic Supplementary Material (ESM) Fig. S1 and slope < 1.3%, ESM Fig. S2). Jeannette creek watershed drains a total area of about 911 ha (Fig. 2). Two poorly drained soils, Clyde silty clay loam and Brookston clay, cover the watershed. Land-use a\/land-cover in the watershed is mostly agricultural area (98.2%). This is a mildly sloped watershed with elevation range of 174 to 186 m (ESM Fig. S1) and slope reaching up to 3% (ESM Fig. S2). Finally, Holtby watershed is a part of Kettle Creek Paired Watershed and drains

An AGNPS model set up requires three spatial datasets: digital elevation model, soil and land-use/land-cover. For the Wigle creek watershed, a 0.5 m  0.5 m resolution hydro-enforced DEM, provided by Essex Region Conservation Authority (ERCA), a soil map obtained from Soil Landscapes of Canada (SLC, 2010), and a landuse/land-cover map provided by the ERCA, were used. Choosing a uniform cell size of 400 m has resulted in a total number of 124 AGNPS cells. For the Jeannette Creek watershed, a 0.5 m  0.5 m DEM was prepared from a contour map at 50 cm interval which was provided by Lower Thames Valley Conservation Authority (LTVCA). Similarly, a land-use/land-cover map prepared by the LTVCA based on field to field survey and a soil map obtained from Soil Landscapes of Canada (SLC, 2010) were also used. A cell size of 200 m was chosen which resulted in a total of 188 AGNPS cells. Similarly, for the Holtby watershed, a 1 m  m light detection and ranging (LIDAR) based DEM (OMAFRA, 2013), a soil map created based on a soil survey (AAFC, 2012), and a land-use/landcover map based on field surveys conducted by Upper Thames Conservation Authority (UTCA) were used. A uniform cell size of 100 m was used in this case, which resulted in a total of 371 AGNPS cells. Based on these spatial maps, a functional AGNPS model then calculates watershed properties (e.g. flow direction, land slope, slope length, channel slope) and derives parameters related to runoff (e.g. curve number, Manning’s roughness coefficient), sediment (e.g. soil erodibility factor, cover factor, conservation practice factor), and phosphorus (e.g. soil phosphorus concentration) for each cell. For all watersheds, the DEM-based stream network was manually verified after several field visits and was used to identify correct flow path. As with other Ontario agricultural watersheds, some fields are tile drained. The case of the Holtby watershed is shown in ESM Fig. S3. While the current version of the AGNPS is not capable of modelling any type of tile drains, it is evident from the literature that tile drains may influence the phosphorus transport profoundly. Tile drainage is shown to have influence on soil permeability and soil matrix, thus affecting macropore flow (Kleinman et al., 2015), and macropore flow is often the main pathway of sub-surface flow (Brady and Weil, 1999). The inability to comprehensively represent tile drains in the AGNPS, is one of the limitations of the work. As the AGNPS is an event-based model, it has to be calibrated against selected storm events. Only a small number of storm events were qualified as the number is constrained by the need to include all four seasons (early and late spring, summer and fall) and availability of precipitation, runoff, sediment and phosphorus data. A total of 7, 8, and 11 storm events were carefully chosen to calibrate AGNPS models for Wigle, Jeannette, and Holtby watershed, respectively (Table 1). Each storm is assigned with one of the four storm types (IA, I, II and III) and one of the three antecedent Soil Moisture Conditions (AMC) (I, II and III) using the Antecedent Precipitation Index (API), as suggested by Perrone and Madramootoo (1998). Each storm type has different intensity (type IA: 0–1 mm/hour, type I: 0.8–2 mm/hour, type II: 1.6–2.6 mm/ hour, and type III: >2.4 mm/hour). Similarly, the AMC class of ‘I’ represents a dry soil, the AMC class II represents a soil with average moisture content and the AMC class III represents a nearly saturated soil (Young et al., 1989). For this, precipitation data at 15-,

Please cite this article as: B. Zhang, N. K. Shrestha, R. Rudra et al., Threshold storm approach for locating phosphorus problem areas: An application in three agricultural watersheds in the Canadian Lake Erie basin, Journal of Great Lakes Research, https://doi.org/10.1016/j.jglr.2019.12.003

4

B. Zhang et al. / Journal of Great Lakes Research xxx (xxxx) xxx

Fig. 2. Location of three watersheds; (a) Wigle, (b) Jeannette, and (c) Holtby with respect to Lake Erie and its basin.

Table 1 Characteristics of selected storm events for model calibration of Wigle, Jeannette and Holtby watershed. Season

Storm

(a) Wigle Creek Watershed Early Spring 1 2 Summer 1 2 Fall 1 2 3 (b) Jeannette Creek Watershed Early Spring 1 2 3 Late Spring 1 2 3 Fall 1 2 (c) Holtby Watershed Early Spring 1 2 Late Spring 1 2 Summer 1 2 3 Fall 1 2 3 4

Amount (mm)

Duration (h)

Intensity (mm/h)

Storm Type

Antecedent Soil Moisture Condition

9.5 8.39 14.69 26.24 38.93 10.87 18.71

13 7 4 3 31 9 9

0.73 1.20 3.67 8.75 1.26 1.21 2.08

IA I III III I I II

Nearly Saturated Nearly Saturated Average Average Nearly Saturated Nearly Saturated Average

33.6 12.6 13.2 14.2 59.6 17 27.4 26.6

31 4 5 18 36 9 8 12

1.08 3.15 2.64 0.79 1.66 1.89 3.43 2.22

I III II IA II II III II

Nearly Saturated Nearly Saturated Average Average Average Average Nearly Saturated Average

14.5 9.2 10 23.9 19.3 6.2 11.9 34 8.1 8.5 12.9

5 6 17 21.5 1.5 3 9 4 7 11.5 16

2.90 1.53 0.59 1.11 12.87 2.07 1.32 8.50 1.16 0.74 0.81

III II IA I III II I III I IA IA

Nearly Saturated Nearly Saturated Nearly Saturated Average Dry Dry Dry Average Average Average Average

Please cite this article as: B. Zhang, N. K. Shrestha, R. Rudra et al., Threshold storm approach for locating phosphorus problem areas: An application in three agricultural watersheds in the Canadian Lake Erie basin, Journal of Great Lakes Research, https://doi.org/10.1016/j.jglr.2019.12.003

B. Zhang et al. / Journal of Great Lakes Research xxx (xxxx) xxx

15- and 30-minute interval were used, respectively for Wigle, Jeannette, and Holtby watershed, as provided by their respective conservation authorities. Each AGNPS model was calibrated against surface runoff, sediment and total phosphorus concentration for the selected storm events. The surface runoff data were measured at 5- and 30-minute time interval at the outlet of Jeannette Creek and Holtby watershed, respectively. However, for the Wigle Creek watershed, only instantaneous measurements of surface runoff were available. As for sediment and phosphorus concentrations, only sporadic samples were collected for all watersheds. While a continuous monitoring of sediment and phosphorus data would have been preferred as they would allow a more robust calibration, such a case is possible only when automated samplers are installed (van Griensven et al., 2000). However, water quality samples were collected at important instances such as at pre-event, early-rise, midrise, peak-flow, early-drop, and late-drop of a storm hydrograph thereby addressing some of the issues in such sporadic water quality sampling (Shrestha et al., 2013). While the AGNPS model allows calibration of 22 parameters, not all the parameters were subjected to calibration as it is widely reported that each of the input parameters is not equally important and sensitive (Liu et al., 2008; Mohammed et al., 2004). Therefore, parameters such as curve number (CN), Manning’s roughness coefficient (n), soil erodibility factor (K), cover factor (C), fertilizer application and availability were selected for the AGNPS model calibration. Seasonal variabilities in these factors were considered also. For example, CN number associated with a particular storm depends on the antecedent moisture condition (AMC), soil and land-use types. While the calculation the API using 14-day precipitation value prior to a storm event facilitated in choosing a proper AMC type, judgement was also used to ascertain the realistic soil moisture condition at the initiation of a storm event. In southern Ontario watersheds, a near-saturation soil moisture condition (AMC type III) will prevail during spring periods while summer season will be generally dry (AMC type I). However, a summer storm after several preceding storms will have a high API value and have potential to generate higher surface runoff than in a typical summer soil moisture condition. Hence, the antecedent soil moisture condition in such a case should be elevated to ‘average’ (AMC type II) from ‘dry’ (AMC type I). Similarly, literature values were used to assign seasonal dependent ‘n’ values at each cell respecting existing land cover. For example, agricultural/cultivated areas are left as a fallow during early spring season. As such, an ‘n’ value of 0.05 following NRCS (1986) was used for the season. For the same land-use, an ‘n’ value of 0.2 was used in summer season as vegetation growth during the season offers more resistant to surface runoff. As the K factor, default values as given in the AGNPS manual (Young et al., 1989) were used, respecting different soil types in the watersheds. However, during spring freeze-thaw period condition, measurements at similar southern Ontario agricultural watersheds showed higher K values (Coote et al., 1978; Wall et al., 1988). Hence, it was increased by a factor of 2 during spring season. Similarly, the C factor also varies seasonally due to plant growth and its canopy cover against erosive potential of rainfall. To incorporate this, we used a relationship between normalized difference vegetation index (NDVI) and the C factor as suggested by van der Knijff et al. (2000). As such, the NDVI values for different crops grown in these watersheds were ascertained from the study of Hatfield and Prueger (2010). For other land-use types, default values (forest: 0.008, meadow: 0.3 and urban: 0.01) as per the AGNPS manual (Young et al., 1989), were used. Three widely used goodness-of-fit statistics; coefficient of determination (R2), percentage of bias (PBIAS), and the NashSutcliffe Efficiency (NSE) (Nash and Sutcliffe, 1970) were used for model performance evaluation. Moreover, as per suggestion of

5

Moriasi et al. (2015), four qualitative ratings (very good, good, satisfactory and unsatisfactory) were assigned to model results based on the range of values of the selected goodness-of-fit statistics (ESM Table S1). Results and discussion Calibration results for streamflow, sediment and phosphorus Fig. 3 shows scatter plots of simulated and observed surface runoff, sediment and phosphorus at the outlet of the Wigle, Jeannette, and Holtby watersheds for selected storm events (Table 1). The goodness-of-fit statistics with qualitative rating is presented in Table 2. In general, there is a good agreement between simulated and observed variables, especially for surface runoff. As such, at least a ‘good’ quality of surface runoff results was observed for all watersheds. The PBIAS value calculated for surface runoff at the outlet of the Wigle Creek watershed is quite high (>95%). Such serious overestimation of simulated surface runoff can be related to the nonavailability of continuous streamflow measurements. Such instantaneous measurements of streamflow obviously hindered calculation of observed surface runoff for a storm event. Hence, a triangular outflow hydrograph was constructed utilizing the AGNPS model outputs such as the time to the peak for which the AGNPS assumes to occur at 37.5% of the total duration of a storm (termed as the shape factor; Young et al., 1989) and peak runoff rate. Then, the instantaneous observation was plotted in the constructed triangular hydrograph in order to compare the results and to compute the goodness-of-fit statistics. For the Jeannette and Holtby watersheds, the streamflow simulation results were better, albeit due to the availability of continuous streamflow data. All goodness-of-fit statistics (R2 > 0.95, NSE > 0.85 and PBIAS < 20%, Table 2) estimated ‘good’ and ‘very good’ quality of streamflow results at these watersheds, respectively. Furthermore, the data points are also scattering on both sides of the 1:1 line showing there is no systematic trend of either over- or under-estimation (Fig. 3). As expected, a lower quality of model results for sediment and phosphorus, as compared to surface runoff, has been obtained for all three watersheds. For the Wigle Creek watershed, model seriously overestimated both sediment and phosphorus concentrations, which is partly due to the lower quality of surface runoff simulation result, as described above. For the Jeannette and Holtby watersheds, contrary to streamflow observations, sediment and phosphorus observations, as depicted already, were taken at rising limb, peak, and falling limb of the storm hydrograph. Hence, devoid of continuous measurements, model simulated sediment and phosphorus concentrations were compared against storm-eventaveraged observed concentrations. As such, the sediment concentrations at the Jeannette and Holtby watersheds have been simulated with a ‘satisfactory’ and ‘good’ accuracy while the phosphorus concentrations have been simulated with a ‘good’ accuracy for both watersheds. Our results are indeed comparable to those reported in similar southwestern Ontario watersheds, e.g. Booty et al. (2005). Considering the apparent difficulties in water quality modelling of small agricultural watersheds with sporadic measurements, the modelling results for phosphorus can be considered reasonable to carry out further analyses. Identified threshold storms in different seasons Fig. 4 shows threshold storms corresponding to seasonal phosphorus tolerance limits for three watersheds. As stated, the cali-

Please cite this article as: B. Zhang, N. K. Shrestha, R. Rudra et al., Threshold storm approach for locating phosphorus problem areas: An application in three agricultural watersheds in the Canadian Lake Erie basin, Journal of Great Lakes Research, https://doi.org/10.1016/j.jglr.2019.12.003

6

B. Zhang et al. / Journal of Great Lakes Research xxx (xxxx) xxx

Fig. 3. Comparison of simulated and observed surface runoff (mm), sediment (mg/l), and total phosphorus (mg/l) at Wigle, Jeannette, and Holtby creek watershed, respectively.

Table 2 Goodness-of-fit statistics and qualitative rating of phosphorus simulation results at Wigle, Jeannette, and Holtby watersheds. Watershed

Variables

R2

PBIAS

NSE

Overall Rating

Wigal

Runoff (mm) Sediment (mg/L) Phosphorus (mg/L) Runoff (mm) Sediment (mg/L) Phosphorus (mg/L) Runoff (mm) Sediment (mg/L) Phosphorus (mg/L)

0.98(VG) 0.49(S) 0.1(NS) 0.99(VG) 0.84(VG) 0.71(G) 0.95(VG) 0.94(VG) 0.81(VG)

95.36(NS) 110.06(NS) 102.28(NS) 19.80(NS) 20.05(NS) 23.85(S) 0.31(VG) 5.6(VG) 49.98(NS)

0.76(G) 0.23(NS) 0.03(NS) 0.87(VG) 0.69(S) 0.52(G) 0.91(VG) 0.59(S) 0.65(G)

G NS NS G S G VG G G

Jeannette

Holtby

VG (very good); G (good); S (satisfactory); NS (not satisfactory). R2: Coefficient of Determination; PBIAS: Percentage of Bias; NSE: Nash-Sutcliffe Efficiency.

Please cite this article as: B. Zhang, N. K. Shrestha, R. Rudra et al., Threshold storm approach for locating phosphorus problem areas: An application in three agricultural watersheds in the Canadian Lake Erie basin, Journal of Great Lakes Research, https://doi.org/10.1016/j.jglr.2019.12.003

7

B. Zhang et al. / Journal of Great Lakes Research xxx (xxxx) xxx

Fig. 4. Identified threshold storms during different seasons at three watersheds (Wigle, Jeannette and Holtby).

brated models with representative storms were run for increasing return periods (2-year through 100-year), and the storm event that generated phosphorus load equal to the seasonal tolerance limit was identified as the threshold storm. As expected, phosphorus load increased with storms of higher intensities. The increments for the upland watershed (Holtby) were more severe than that for two flat watersheds (Wigle and Jeannette). Despite the lowest rainfall intensity during early spring season (Fig. 1), phosphorus load was the highest (Fig. 4). In early spring season, soils have high moisture level and less cover. It implies that soil moisture and cover exert greater control to phosphorus loss compared to the rainfall intensity in such agricultural watersheds. It was found that threshold storms of return 12.4 (8.1–18.2) [mean (lower–upper)], 14.9 (9.0–21.5), and 4.8 (3.7–6.1) years (Table 3) will generate phosphorus load equal to early spring seasonal tolerance limit of 0.42 ± 0.07 kg/ha for the Wigle, Jeannette and Holtby watershed, respectively. The severity of phosphorus loss problem in the upland watershed (Holtby) has been duly reflected in the threshold storm of very small return period (mean value of 4.8year, Table 3). Threshold storms at the two low relief watersheds are fairly comparable. Hence, it is evident that topographical variations play an important role in phosphorus load generation in such agricultural watersheds. During late spring period, threshold storms were becoming getting increasingly rare as their return period increased to 29.9 (20.6– 41.8), 39.9 (25.5–59.5) and 10.3 (8.2–13.6) a year for the Wigle,

Jeannette and Holtby watershed, respectively. Interestingly, summer storms up to 100-year return period will not generate phosphorus that will reach seasonal tolerance limit of 0.74 ± 0.12 kg/ ha. Although, representative summer storms have the highest rainfall intensity (Fig. 1), such low phosphorus load (<0.1 kg/ha, Fig. 4) can be attributed to soil moisture and vegetative cover conditions. During summer season, agricultural fields will have higher canopy cover and soils have depleted soil moisture level. Such a condition was indeed reported in a similar watershed of southern Ontario (Golmohammadi et al., 2017). Hence, summer storms will not pose any serious phosphorus problems. Similarly, during fall season, storms with a return period of 33.5 (20.9–48.0) 44.3 (26.9–67.4) and 13.2 (9.2–17.0) years, respectively for the Wigle, Jeannette and Holtby watershed, will generate phosphorus load equal to seasonal threshold limit of 0.28 ± 0.05 kg/ha. These storms are less frequent than that observed for both early and late spring periods. Hence, from the results, it was found that early spring season is the hot moment of phosphorus load, and that the threshold storm in the upland watershed is more frequent than that in flat watersheds. Furthermore, summer storms are found not to be problematic at all. Location of critical source areas Based on the derived seasonal threshold tolerance limits, critical source areas of phosphorus yield were located. As such, the

Table 3 Seasonal maximum acceptable total phosphorus and corresponding threshold storm events for three watersheds (Wigle, Jeannette and Holtby). Seasons

Maximum acceptable Total Phosphorus (kg/ha)

Threshold storm return period of Wigle creek watershed (yr)

Threshold storm return period of Jeannette creek watershed (yr)

Threshold storm return period of Holtby watershed (yr)

Early Spring Late Spring Summer Fall

0.42 0.42 0.74 0.28

14.9 (9.0–12.5) 39.9 (25.5–59.5) – 44.3 (26.9–67.4)

12.4 (8.1–18.2) 29.9 (20.6–41.8) – 33.5 (20.9–48.0)

4.8 (3.7–6.1) 10.3 (8.2–13.6) – 13.2 (9.2–17.0)

(0.35–0.49) (0.35–0.49) (0.62–0.86) (0.23–0.33)

Please cite this article as: B. Zhang, N. K. Shrestha, R. Rudra et al., Threshold storm approach for locating phosphorus problem areas: An application in three agricultural watersheds in the Canadian Lake Erie basin, Journal of Great Lakes Research, https://doi.org/10.1016/j.jglr.2019.12.003

8

B. Zhang et al. / Journal of Great Lakes Research xxx (xxxx) xxx

severity of cell phosphorus yield was partitioned into below and above threshold tolerance limit for a particular season. For instance, the mean phosphorus tolerance limit for spring season was found to be 0.42 kg/ha, and cell phosphorus yield less than or equal to the limit was grouped a class (‘Below Threshold’). Four more classes for cell phosphorus yield more than the tolerance limit were defined as ‘Moderate’ (0.43–0.84 kg/ha), ‘High’ (0.85–1.68 kg/ha), ‘Very High’ (1.69–3.36 kg/ha) and ‘Extremely High’ (>3.36 kg/ha). Figs. 5–7 show spatial distribution of cell phosphorus yield as generated by the representative storms of increasing return periods (2-year through 100-year) in different seasons. It is clear from the spatial plots that early spring period is the hot moment of phosphorus yield for all the watersheds whereas summer storms up to 100-years of return period, will generate phosphorus yield that will always be below the seasonal threshold limit. Even a 2-year early spring storm will have more than 60% (>70% in the Wigle, 77% in the Jeannette and 60% in the Holtby watershed, ESM Table S2) of the areas generating phosphorus yield more than the mean seasonal threshold limit of 0.42 kg/ha (Fig. 8). For such lower return period storms, the uncertainty in percentage of cells exceeding the seasonal phosphorus threshold limit was non-negligible. For example, in the upland watershed (Holtby), while considering the lower limit of the threshold value (0.35 kg/ ha), 83.6% of the cells have become contributing areas while they

reduced to 37.7% when considering the upper limit (0.49 kg/ha). As expected, for higher return periods (e.g. 25-year), the percentage of contributing area will increase steadily, and will reach to a pseudo-equilibrium level around 95%. As such, the uncertainty in percentage of cells exceeding the seasonal phosphorus threshold limit is reduced. For the most extreme storm in early spring season (100-year), the phosphorus yield at the outlet of two low lying watersheds (Wigle and Jeannette) is almost the same (0.82 kg/ha and 0.90 kg/ha, respectively) whereas the yield was almost double (1.79 kg/ha; Fig. 6, ESM Table S2) in the upland watershed (Holtby). It is thus evident that topographic features (e.g. slope) exert a greater control on phosphorus yield dynamics as the upland watershed (Holtby) has cells which have slope more than 6% while two low lying watersheds have slope less than 1.5% and 3.5%, respectively. This was further illustrated in Fig. 9 which shows that cell phosphorus yield and slope are positively correlated in the upland watershed. Contrary to early spring, there will be significant decrease in phosphorus yield at both cell and watershed outlets during late spring season in all three watersheds, especially for less extreme storms (up to 10-year). In all watersheds, the 2-year storm will not generate phosphorus yield that would exceed the seasonal phosphorus limit (mean value of 0.42 kg/ha). Therefore, none of the cells will constitute critical source areas for phosphorus. However, as the severity of the storm increases, the contribution area

Fig. 5. Cell phosphorus yield for representative storms of increasing return periods in different seasons in Wigle creek watershed.

Please cite this article as: B. Zhang, N. K. Shrestha, R. Rudra et al., Threshold storm approach for locating phosphorus problem areas: An application in three agricultural watersheds in the Canadian Lake Erie basin, Journal of Great Lakes Research, https://doi.org/10.1016/j.jglr.2019.12.003

B. Zhang et al. / Journal of Great Lakes Research xxx (xxxx) xxx

9

Fig. 6. Cell phosphorus yield for representative storms of increasing return periods in different seasons in Jeannette creek watershed.

fraction will increase steadily (Fig. 8, ESM Table S2). For more extreme storms (>50-year), there are virtually no differences in percentage of contributing areas. However, there exist significant differences on phosphorus yield at watershed outlet such as 1.33 kg/ha in late spring compared to 1.79 kg/ha in early spring season at the outlet of the Holtby watershed. During summer season, phosphorus yield always remained below the seasonal threshold limit (mean value 0.74 kg/ha) for all three watersheds. Phosphorus yield at the outlet of the Wigle, Jeannette and Holtby watersheds will only reach to a value of 0.04, 0.06 and 0.09 kg/ha, respectively (Fig. 8, ESM Table S2). Hence, any summer storm (up to 100-year) would not result in CSAs for phosphorus. During fall season, similar to the pattern observed during late spring, a lower percentage of cells will contribute phosphorus yield higher than the seasonal tolerance limit for less extreme storms (<10-year) which will increase steadily. It will reach to the same level as that observed during the early spring case for more extreme storms (>50-year). However, the total phosphorus yield at the watershed outlet will be significantly lower compared to that during early spring period (Fig. 8, ESM Table S2). For all seasons, cell phosphorus yield doesn’t show any spatial pattern of phosphorus hot-spots (Figs. 5–7). This is interesting as a previous study (Shrestha et al., 2019) reported that cells in vicin-

ity of the streams have had higher sediment yield on account of higher connectivity of those cells compared to faraway cells. This may be related to model’s inability to model tile drains explicitly. As stated earlier, provision of tile drains generally increase soil permeability which in turn would decrease surface runoff. As such, an enhanced macropore flow, the main connection pathway of subsurface flow, would lead to increased mineral phosphorus loss from tile drained portion of a watershed. Explicit consideration of tile drains in the AGNPS model might be able to show distinct patterns of phosphorus loading from the watersheds. In fact, such a uniform pattern of CSAs is only possible if contribution of the mineral phosphorus is higher than that of the organic phosphorus. In contrast, a significant and positive correlation between cell phosphorus and sediment yield (Fig. 9) indicates a non-negligible contribution of the organic phosphorus. Moreover, an area with higher topographic wetness index (TWI) is believed to generate high surface runoff and as such high phosphorus yield. However, we found that such areas may not necessarily generate high phosphorus load, as poor correlation coefficients (R2 < 0.02, between the TWI and cell phosphorus yield, Fig. 9) are observed for all three watersheds. In hindsight, the exact location of CSAs of phosphorus yield is difficult in such upland and flat watersheds, and that for a BMP to be effective, it should be placed virtually all over the watershed.

Please cite this article as: B. Zhang, N. K. Shrestha, R. Rudra et al., Threshold storm approach for locating phosphorus problem areas: An application in three agricultural watersheds in the Canadian Lake Erie basin, Journal of Great Lakes Research, https://doi.org/10.1016/j.jglr.2019.12.003

10

B. Zhang et al. / Journal of Great Lakes Research xxx (xxxx) xxx

Fig. 7. Cell phosphorus yield for representative storms of increasing return periods in different seasons in Holtby watershed.

Fig. 8. Cells (%) that exceeds phosphorus seasonal tolerance value and total phosphorus (kg/ha) at watershed outlet during early spring, late spring, summer and fall respectively.

Please cite this article as: B. Zhang, N. K. Shrestha, R. Rudra et al., Threshold storm approach for locating phosphorus problem areas: An application in three agricultural watersheds in the Canadian Lake Erie basin, Journal of Great Lakes Research, https://doi.org/10.1016/j.jglr.2019.12.003

B. Zhang et al. / Journal of Great Lakes Research xxx (xxxx) xxx

11

Fig. 9. Cell phosphorus yield as function of sediment yield, slope and topographic wetness index (TWI) for three watersheds.

Conclusions

Declaration of Competing Interest

Identification of hot-spots and hot moments of phosphorus load in small agricultural watersheds is crucial in designing and implementing best management practices. While phosphorus load from a watershed depends on the watershed’s topographic features and other watershed characteristics, the type and intensity of storms are important, too. As not all storms will generate phosphorus which can be considered problematic, a standardized approach is thus necessary for systematic identification of critical source areas of phosphorus yield. In this scope, a threshold storm is thus defined as the maximum storm intensity in which the phosphorus generated in a watershed is below seasonal phosphorus tolerance limit. In order to evaluate the threshold storm approach, three separate AGNPS models for three diverse small agricultural watersheds in southern Ontario, Canada were calibrated for runoff volume, sediment yield, and total phosphorus. The models were then run for representative storms with increasing return periods (2-year through 100-year). Results showed that in an upland watershed, more frequent storms (4.8-year) in an early spring period will generate phosphorus load which can be considered problematic. The seasonal threshold tolerance phosphorus limit is only exceeded by rarer storms (>12.4-year) in the lowland watersheds. In all three watersheds, up to 100-year summer storms will fail to reach the seasonal tolerance limit for phosphorus, implying that summer storms are not problematic at all. The critical source areas, identified based on the threshold storm approach, are found to be located all over the watersheds. As the phosphorus problem is essentially a source problem, such a simple yet robust approach in identifying critical source areas of phosphorus can be useful in designing costeffective best management practices.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement We would like to thank Essex Region Conservation Authority (ERCA), Lower Thames Valley Conservation Authority (LTVCA) and Upper Thames Conservation Authority (UTCA), respectively for providing spatial, field and crop management, hydrometeorological and water quality data of Wigle, Jeannette and Holtby watersheds. We would also like to appreciate the funding of Ontario Soil and Crop Improvement Association (OSCIA), and Ministry of the Environment, Conservation and Parks (MECP) for partial funding. Furthermore, we would like to thank personnel of the Ontario Ministry of Agriculture, Food and Rural Affairs (OMAFRA) and the MECP for technical support. Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.jglr.2019.12.003. References AAFC, 2012. Soil Survey Reports for Ontario. Agriculture and Agri-Food Canada (AAFC), Canada Bagnold, R.A., 1966. An approach to the sediment transport problem from general physics. Geological Survey Professional Paper 422-I, U. S. Govt. Print. Off, pp. I1– I37

Please cite this article as: B. Zhang, N. K. Shrestha, R. Rudra et al., Threshold storm approach for locating phosphorus problem areas: An application in three agricultural watersheds in the Canadian Lake Erie basin, Journal of Great Lakes Research, https://doi.org/10.1016/j.jglr.2019.12.003

12

B. Zhang et al. / Journal of Great Lakes Research xxx (xxxx) xxx

Booty, W., Lam, D., Bowen, G., Resler, O., Leon, L., 2005. Modelling changes in stream water quality due to climate change in a Southern Ontario watershed. Can. Water Resour. J. 30, 211–226. Brady, N.C., Weil, R.R., 1999. The Nature and Properties of Soils. Prentice Hall, New Jersey. Coote, D.R., MacDonald, E.M., Dickinson, W.T., 1978. Agricultural watershed studies in the Canadian Great Lakes drainage basin; Final summary report. PLUARG Technical Report Series. pp. 1–78. Daggupati, P., Douglas-Mankin, K.R., Sheshukov, A.Y., Barnes, P.L., Devlin, D.L., 2011. Field-level targeting using SWAT: mapping output from HRUs to fields and assessing limitations of GIS input data. T. ASABE. 54, 501–514. Dickinson, T., 1976. Seasonal variability of rainfall extremes. Atmosphere 14, 288– 296. Dickinson, T., 1977. Rainfall intensity-frequency relationships from monthly extremes. J. Hydrol. 35, 137–145. Dickinson, W.T., Rudra, R.P., Wall, G.J., 1986. Identification of soil erosion and fluvial sediment problems. Hydrol. Process. 1, 111–124. Golmohammadi, G., Rudra, R., Dickinson, T., Goel, P., Veliz, M., 2017. Predicting the temporal variation of flow contributing areas using SWAT. J. Hydrol. 547, 375– 386. Grunwald, S., Norton, L.D., 1999. An AGNPS-based runoff and sediment yield model for two small watersheds in Germany. Am. Soc. Agr. Eng. 42, 1723–1731. Hanief, A., Laursen, A.E., 2019. Meeting updated phosphorus reduction goals by applying best management practices in the Grand River watershed, southern Ontario. Ecol. Eng. 130, 169–175. Haregeweyn, N., Yohannes, F., 2003. Testing and evaluation of the agricultural nonpoint source pollution model (AGNPS) on Augucho catchment, western Hararghe. Ethiopia. Agr. Ecosyst. Environ. 99, 201–212. Hatfield, J.L., Prueger, J.H., 2010. Value of using different vegetative indices to quantify agricultural crop characteristics at different growth stages under varying management practices. Remote Sens. 2, 562–578. Kleinman, P.J.A., Smith, D.R., Bolster, C.H., Easton, Z.M., 2015. Phosphorus fate, management, and modeling in artificially drained systems. J. Environ. Qual. 44, 460–466. León, L.F., Booty, W.G., Bowen, G.S., Lam, D.C.L., 2004. Validation of an agricultural non-point source model in a watershed in southern Ontario. Agr. Water Manage. 65, 59–75. Liu, J., Zhang, L., Zhang, Y., Hong, H., Deng, H., 2008. Validation of an agricultural non-point source (AGNPS) pollution model for a catchment in the Jiulong River watershed, China. J. Environ. Sci. 20, 599–606. Liu, Y., Yang, W., Leon, L., Wong, I., McCrimmon, C., Dove, A., Fong, P., 2016. Hydrologic modeling and evaluation of Best Management Practice scenarios for the Grand River watershed in Southern Ontario. J. Great Lakes Res. 42, 1289– 1301. Ma, Y., Bartholic, J.F., 2003. GIS based AGNPS assessment model in a small watershed. Nat. Sci. 1 (1), 50–56. McDowell, R., An, S., Folmar, G., 2000. Phosphorus export from an agricultural watershed: Linking source and transport mechanisms. J. Environ. Qual. 30, 1587–1595. McQueen, D.J., Lean, D.R.S., 1987. Influence of water temperature and nitrogen to phosphorus ratios on the dominance of blue-green algae in Lake St. George, Ontario. Can. J. Fish. Aquat. Sci. 44, 598–604. Meals, D.W., Sharpley, A.N., Osmond, D.L., 2012. Lessons learned from the NIFACEAP: identifying critical source areas. NC State University, Raleigh, NC. 1–7, Raleigh. Merriman, K.R., Daggupati, P., Srinivasan, R., Hayhurst, B., 2019. Assessment of sitespecific agricultural Best Management Practices in the Upper East River watershed, Wisconsin, using a field-scale SWAT model. J. Great Lakes Res. 45, 619–641. Merriman, K.R., Daggupati, P., Srinivasan, R., Toussant, C., Russell, A.M., Hayhurst, B., 2018a. Assessing the impact of site-specific BMPs using a spatially explicit, field-scale SWAT model with edge-of-field and tile hydrology and water-quality data in the Eagle Creek Watershed. Ohio. Water 10 (10), 1299. Merriman, R.K., Russell, M.A., Rachol, M.C., Daggupati, P., Srinivasan, R., Hayhurst, A. B., Stuntebeck, T.D., 2018b. Calibration of a field-scale Soil and Water Assessment Tool (SWAT) model with field placement of best management practices in Alger Creek, Michigan. Sustainability 10 (3), 851. Michalak, A.M., Anderson, E.J., Beletsky, D., Boland, S., Bosch, N.S., Bridgeman, T.B., Chaffin, J.D., Cho, K., Confesor, R., Daloglu, I., DePinto, J.V., Evans, M.A., Fahnensteil, G.L., He, L., Ho, J.C., Jenkins, L., Johengen, T.H., Kuo, K.C., LaPorte, E., Liu, X., McWilliams, M.R., Moore, M.R., Posselt, D.J., Richards, R.P., Scavia, D., Steiner, A.L., Verhamme, E., Wright, D.M., Zagorski, M.A., 2013. Record-setting algal bloom in Lake Erie caused by agricultural and meteorological trends consistent with expected future conditions. P. Nat. Acad. of Sci. 110 (16), 6448– 6452. Mohammed, H., Yohannes, F., Zeleke, G., 2004. Validation of agricultural non-point source (AGNPS) pollution model in Kori watershed, South Wollo, Ethiopia. Int. J. Appl. Earth Obs. Geoinf. 6, 97–109. Moriasi, D.N., Gitau, M.W., Pai, N., Daggupati, P., 2015. Hydrologic and water quality models: performance measures and evaluation criteria. T. ASABE. 58, 1763– 1785. Natural Resources Conservation Services (NRCS), 1986. Technical release-55: Urban hydrology for small watershed. Nash, J.E., Sutcliffe, J.V., 1970. River flow forecasting through conceptual models part I — a discussion of principles. J. Hydrol. 10, 282–290.

Niraula, R., Kalin, L., Srivastava, P., Anderson, C.J., 2013. Identifying critical source areas of nonpoint source pollution with SWAT and GWLF. Ecol. Model. 268, 123–133. OMAFRA, 2013. LIDAR – Digital Terrain Model. Ontario Ministry of Agriculture, Food and Rural Affairs. Otieno, H., Onyando, J., 2012. Coupling agricultural non-point source (AgNPS) model and geographic information system (GIS) tools to predict peak runoff and sediment generation in the Upper River Njoro catchment in Kenya. Int. J. Water. Res. Environ. Eng. 4 (12), 397–403. Paerl, H.W., Hall, N.S., Calandrino, E.S., 2011. Controlling harmful cyanobacterial blooms in a world experiencing anthropogenic and climatic-induced change. Sci. Total Environ. 409, 1739–1745. Perrone, J., Madramootoo, C.A., 1998. Improved curve number selection for runoff prediction. Can. J. Civil Eng. 25, 728–734. Reid, K., Schneider, K., McConkey, B., 2018. Components of phosphorus loss from agricultural landscapes, and how to incorporate them into risk assessment tools. Front. Earth Sci. 6, 1–15. Rodriguez, H.G., Popp, J., Maringanti, C., Chaubey, I., 2011. Selection and placement of best management practices used to reduce water quality degradation in Lincoln Lake watershed. Water Resour. Res. 47, 1–13. Rudra, R.P., Dickinson, W.T., Ahmed, S.I., Patel, P., Zhou, J., Gharabaghi, B., Khan, A.A., 2015. Changes in rainfall extremes in Ontario. Int. J. Environ. Res. 9, 1117–1126. Rudra, R.P., Dickinson, W.T., Clark, D.J., Wall, G.J., 1986. GAMES—a screening model of soil erosion and fluvial sedimentation on agricultural watershed. Can. Water Resour. J. 11, 58–71. Scavia, D., David, A.J., Arend, K.K., Bartell, S., Beletsky, D., Bosch, N.S., Brandt, S.B., Briland, R.D., Daloglu, I., DePinto, J.V., Dolan, D.M., Evans, M.A., Farmer, T.M., Goto, D., Han, H., Hook, T.O., Knight, R., Ludsin, S.A., Mason, D., Michalak, A.M., Richards, R.P., Roberts, J.J., Rucinski, D.K., Rutherford, E., Schwab, D.J., Sesterhenn, T.M., Zhang, H., Zhou, Y., 2014. Assessing and addressing the reeutrophication of Lake Erie: central basin hypoxia. J. Great Lakes Res. 40, 226– 246. Schertz, D.L., 1983. The basis for soil loss tolerances. J. Soil Water Conserv. 38, 10– 14. Sharpley, A.N., Kleinman, P.J.A., Flaten, D.N., Buda, A.R., 2011. Critical source area management of agricultural phosphorus: experiences, challenges and opportunities. Water Sci. Technol. 64, 945–952. Shrestha, N.K., Allataifeh, N., Rudra, R., Daggupati, P., Goel, P.K., Dickinson, T., 2019. Identifying threshold storm events and quantifying potential impacts of climate change on sediment yield in a small upland agricultural watershed of Ontario. Hydrol. Process. 33, 920–931. Shrestha, N.K., Leta, O.T., De Fraine, B., van Griensven, A., Bauwens, W., 2013. OpenMI-based integrated sediment transport modelling of the river Zenne, Belgium. Environ. Modell. Softw. 47, 193–206. Shrestha, N.K., Wang, J., 2019. Water quality management of a cold climate region watershed in changing climate. J. Environ. Inform. 35 (1), 56–80. SLC, 2010. Soil Landscapes of Canada version 3.2: digital map and database at 1:1 million scale. Soil Landscapes of Canada Working Group (SLCWG), Agriculture and Agri-Food Canada (AAFC), Canada, Ontario. Smith, R.E., Williams, J.R., 1980. Simulation of the surface hydrology. In: Knisel, W. (Ed.), CREAMS, a field scale model for chemicals, runoff, and erosion from agricultural management systems. USDA ARS Conservation Research Report 26, Washington DC, USA. Stang, C., Gharabaghi, B., Rudra, R., Golmohammadi, G., Mahboubi, A.A., Ahmed, S.I., 2016. Conservation management practices: success story of the Hog Creek and Sturgeon River watersheds, Ontario, Canada. J. Soil Water Conserv. 71, 237–248. Sun, J.H., Zhu, Q.D., Yan, Z.J., Lu, H.M., Wang, H.R., 2009. A review of research and application of AGNPS model. Adv. Wat. Sci. 20, 876–884. US-SEA, 1980. CREAMS: A Field Scale Model for Chemicals, Runoff, and Erosion from Agricultural Management Systems, Washington, D.C. van der Knijff, M.J., Jones, R.J.A., Montanarella, L., 2000. Soil Erosion Risk Assessment in Europe. European Soil Bureau, European Commission. van Griensven, A., Vandenberghe, V., Bols, J., De Pauw, N., Goethals, P., Goethals, P., Meirlaen, J., Vanrolleghem, P.A., Vooren, L.V., Bauwens, W., 2000. Experience and organisation of automated measuring stations for river water quality monitoring. 1st World Congress of the International Water Association, July 3– 7, 2000, Paris, France. Wall, G.J., Bos, A.W., Marshall, A.H., 1996. The relationship between phosphorus and suspended sediment loads in Ontario watersheds. J. Soil Water Conserv. 51, 504–507. Wall, G.J., Dickinson, W.T., Rudra, R.P., Coote, D.R., 1988. Seasonal soil erodibility variation in Southwestern Ontario. Can. J. Soil Sci. 68, 417–424. White, M.J., Storm, D.E., Busteed, P.R., Stoodley, S.H., Phillips, S.J., 2009. Evaluating nonpoint source critical source area contributions at the watershed scale. J. Environ. Qual. 38, 1654–1663. Winchell, M.F., Folle, S., Meals, D., Moore, J., Srinivasan, R., Howe, E.A., 2015. Using SWAT for sub-field identification of phosphorus critical source areas in a saturation excess runoff region. Hydrol. Sci. J. 60, 844–862. Yang, W., Liu, Y., Simmons, J., Oginskyy, A., McKague, K., 2013. SWAT Modelling of Agricultural BMPs and Analysis of BMP Cost Effectiveness in the Gully Creek Watershed. OMAFRA, UoG. Young, R.A., Onstad, C.A., Bosch, D.D., Anderson, W.P., 1989. AGNPS: a nonpointsource pollution model for evaluating agricultural watersheds. J. Soil Water Conserv. 44, 168–173.

Please cite this article as: B. Zhang, N. K. Shrestha, R. Rudra et al., Threshold storm approach for locating phosphorus problem areas: An application in three agricultural watersheds in the Canadian Lake Erie basin, Journal of Great Lakes Research, https://doi.org/10.1016/j.jglr.2019.12.003