Assessing the Impact of BMPs on Water Quality and Quantity in a Flat Agricultural Watershed in Southern Ontario

Abstract

Non-point source pollution poses a continuous threat to the quality of Great Lakes waters. To abate this problem, the Great Lakes Agricultural Stewardship Initiative (GLASI) was initiated in Ontario, Canada, with the primary aim of reducing phosphorus pollution. Therefore, a case-study analysis of the Wigle Creek watershed, one of the six priority watersheds under the GLASI program, was undertaken to evaluate the effectiveness of various existing and potential future Best Management Practices (BMPs) and to identify BMPs that might aid in mitigating the watershed’s contribution to phosphorus loads reaching Lake Erie. Given the watershed’s very flat topography, hydrological/nutrient modeling proved an extremely challenging exercise. The Soil and Water Assessment Tool (SWAT) model was used in this evaluation. Several digital elevation model (DEM) options were considered to accurately describe the watershed and represent flow conditions. A 30 m resolution DEM, implementing a modified burning in of streams based on ground truthing, was finally employed to develop the SWAT model’s drainage framework. The model was first calibrated for flow, sediment, and phosphorus loads. The calibrated model was used to evaluate the ability of potential BMPs (minimum tillage, no-till, retiring croplands into pasture, retiring croplands into forest, winter wheat cover crop, and vegetative filter strips) to reduce phosphorus loads compared to implemented practice. Converting all croplands into pasture or forest significantly decreased P loads reaching Lake Erie. Comparatively, a winter wheat cover crop had minimal effect on reducing phosphorus loading.

1. Introduction

Ontario’s population principally resides within catchments draining into the Great Lakes Basin (GLB), thereby making it an extremely valuable resource for Ontario’s residents [1]. Likewise, a sizable proportion of the United States’ north-central population also resides in the GLB. However, in the last few decades, the water quality of Lake Erie (one of the Great Lakes) has faced continuous threats due to increased levels of harmful pollutants, particularly phosphorus (P) [1,2,3]. Non-point source (NPS) pollution engendered by agricultural activities in the Western Lake Erie Basin (WLEB) has been identified as the dominant limiting nutrient source for Lake Erie [4].

Generally accepted methods to address hydrology and water quality problems in agricultural sites [5,6,7], agricultural BMPs are mainly focused on addressing nutrient management along with erosion and runoff control. The prevention of soil erosion and attendant reduction in the transport of nutrients can be achieved through the installation of barriers and buffers to intercept sediments and nutrients leaving the field (Agriculture and Agri-food Canada, Ottawa, ON, Canada, 2000). BMPs can be implemented at the field-scale or watershed-scale. The impact of BMPs in North America has been the subject of numerous studies [8,9]. One study [10] used an ACT (avoid, control, and trap) approach for BMPs implemented under an Agricultural Conservation Planning Framework (ACPF), thereby facilitating the acquisition of consistent input data and enabling consistent planning analyses in different regions. A few studies have assessed the impact of the implementation of BMPs in Ontario watersheds, particularly those draining into Lake Erie [4,9,11]. Although field monitoring can provide adequate data for evaluating these BMPs, it is a time- and labor-intensive procedure, further complicated by the substantial additional costs and logistical difficulties involved in monitoring during Ontario’s sub-zero winters and frequent freeze/thaw cycles. These limitations could only be overcome through the provision of extensive financial resources. Hydrological modeling constitutes a practical alternative, as models can simulate several land management scenarios. Provided the model adequately represents the agricultural management setting, this approach saves time, labor, and other expenses associated with field investigations [9,12,13,14]. In the present study, the Soil and Water Assessment Tool SWAT 2012 [15] developed by the U.S. Department of Agriculture-Agricultural Research Service (USDA-ARS), Texas A & M university, Texas, USA was selected to evaluate flow, sediment, and phosphorus transport from the Wigle Creek watershed under the potential implementation of various BMPs.

The SWAT model has been widely used to simulate flow/water quality within American watersheds contributing to Lake Erie’s waters. Bosch et al. [16] simulated six watersheds contributing to Lake Erie using the SWAT model and evaluated streamflow and sediment load along with N and P concentrations contributed by the different watersheds. Moreover, several studies [17,18,19,20,21] focused on large-scale watersheds (e.g., Western Lake Erie basin; Maumee basin) that drain into Lake Erie have successfully employed SWAT to investigate the impacts of changes in land use, agricultural practices, and climate on the effectiveness of BMPs in reducing NPS pollution [1]. These studies have offered valuable insights for policymakers and decision-makers. On the Canadian side, application of the SWAT has been limited to small to medium-sized watersheds [11,12,22,23]. For instance, the use of the SWAT model by Liu et al. [11] simulated the impact of several BMPs in Southern Ontario’s Grand River watershed (6800 km²). Nutrient management and wetland restoration were shown to offer significant impacts on watershed-scale nutrient reduction. However, modeling exercises for watersheds having their outlets on the Canadian shores of Lake Erie are limited [24]. Hence, Daggupati et al. [1] conducted an assessment of the various input data types to analyze their impact on simulating hydrological processes and streamflow. The SWAT model was employed for this investigation, covering the entirety of the Lake Erie Basin.

The implementation of the SWAT model for water quality simulation remains a challenge when modeling watersheds with a flat topography. Donmez et al. [25] improved the SWAT model’s ability to accurately simulate daily flow and NO₃ concentrations in a flat, data-scarce agricultural watershed situated on Turkey’s Lower Seyhan Plain (LSP). Another study by Schmalz et al. [26] evaluated the capacity of the SWAT model to model water quality parameters for flat, low hydraulic gradient watersheds in the north German lowlands. These studies concluded that topography is one of the most sensitive inputs for parameterization, directly affecting model behavior. In a further study, Habeck et al. [27], using the SWIM model to simulate nitrogen flows in northeastern Germany’s Nuthe lowlands, drew similar conclusions. Employing SWAT to analyze the combined effect of topography and land cover on runoff generation in a low-slope/low variance watershed, Al-Khafaji et al. [28] concluded that the quality of the digital elevation model (DEM) greatly altered the accuracy of watershed and stream network delineation.

The Wigle Creek watershed has extreme flat topography, thereby making phosphorus modeling extremely challenging. To the best of our information, evaluations of the impact of existing/proposed BMPs on water quality in flat watersheds in Ontario are sparse [29]. Therefore, a modeling approach employing SWAT was developed to evaluate the effectiveness of BMPs in improving water quality in an Ontario low-slope watershed draining into Lake Erie. The specific issues associated with this novel modeling approach, addressed in this paper, are to (a) evaluate the suitability of the SWAT model’s application to the extremely flat Wigle Creek watershed in southern Ontario, for which limited datasets are available, and (b) use the model to simulate the potential impact of a watershed-scale implementation of BMPs on sediment and phosphorus loading.

2. Materials and Methods

2.1. Description of the Study Area

The Wigle Creek watershed—one of six GLASI project priority watersheds—was selected for this study. To be precise, a sub-basin of the Wigle Creek watershed, lying on the western branch of Wigle Creek and draining an area of approximately 1949 ha., stayed focused on water quality and quantity modeling and BMP evaluations. The study watershed lies between 82°47′30″ W to 82°43′30″ W longitude and 42°07′30″ N to 42°03′30″ W latitude. The watershed is relatively flat, with the elevation varying from 185 to 203 m (Figure 1). Agriculture being the predominant land use within the watershed, the territory is extensively tile-drained. The main crops are corn (Zea mays L.) and soybeans (Glycine max (L.) Merr.). The predominant soil type in the study watershed is Brookston clay. This watershed was selected because it is representative of watersheds in this area, having a flat topography, clay soils, and row crop agriculture. In addition, the GLASI program included a cost-sharing program for farmers to implement BMPs as well as an outreach component. This sub-watershed is home to environmentally conscious farmers who help promote conservation practices.

Data Collection

The spatial datasets like DEM, land use maps, and soil map were collected from OMAFRA, Windshield survey report, and the Soil Landscapes of Canada, respectively (

Table 1. Details of data input for the Wigle Creek watershed.

Data	Source	Description
DEM	OMAFRA	Digital elevation model (30 m × 30 m)
Land use	Windshield survey report	Land use for years 2004 and 2006
Soil	Soil Landscapes of Canada (SLC)	Daily data (July 2016–July 2017)
Precipitation	‘Jack Miner’ Environment Canada weather station, lacunae filled from the ‘Harrow CDA’ station	Daily data (July 2016–July 2017)
Crop/crop management	Essex Region Conservation Authority	List of crops grown: 2012–2016: from a 5-year survey 2016–2017 from a 2-year windshield survey

). The climate datasets were collected from Environment Canada weather station ‘Jack Miner’, with gaps filled from the nearby environment station ‘Harrow CDA’. The Water Quality and Quantity dataset samples were collected by the Essex Region Conservation Authority (ERCA) using an ISCO autosampler twice a week year-round, plus during rain/snow melt events. During the events, discrete samples were taken on the rising limb, at the peak, and on the falling limb of the hydrograph, and these averaged samples were used for modeling purposes. Level loggers and bubblers were used to measure water depth at a 15 min interval. Most of the discrete flow measurements were taken during the base flow to midflow conditions due to safety constraints, which limit the accuracy of a rating curve. Later, these samples were analyzed for total suspended sediment (TSS), total phosphorus (TP) concentration, and soluble reactive phosphorus (SRP) concentration, as well as various nutrient forms, by Caducean Environmental Laboratories (http://www.caduceonlabs.com/resources.php (accessed on 26 January 2018)).

2.2. Data Preparation for Model Inputs

Geospatial data used for setting up the SWAT model included a DEM, weather datasets (mainly precipitation and minimum and maximum temperatures), soil type, land use, and land management data. Most of this data were obtained from the Essex Region Conservation Authority (ERCA), which assembled such information as part of the GLASI priority sub-watershed initiative. A detailed description of the input data (Table 1), its sources, and resolution follows.

2.2.1. Digital Elevation Model (DEM)

Since Wigle Creek is an extremely flat watershed, two DEMs were drawn upon to set up the SWAT model. Initially, a Lidar Hydro-Enforced DEM (0.5 m × 0.5 m) processed by the Essex Region Conservation Authority was used to set up the model (ERCA, 2016); however, some issues with this high-resolution DEM for SWAT setup were observed (see Section 3). Therefore, a 30 m × 30 m DEM was obtained from OMAFRA [30] and used in setting up the model (Figure 2a).

2.2.2. Soil Data

A provincial-level soil database (Ontario), included within Soil Landscapes of Canada (SLC) version 3.2 (Soil Landscapes of Canada Working Group 2007, Ontario, Canada. https://sis.agr.gc.ca/cansis/nsdb/slc/index.html accessed on 20 May 2018), was employed. The SLC data contain a soil map of Canada, which provides major characteristics of the soils for the whole country. The SLC data were compiled at a scale of 1:1 million, with each polygon on the map describing a distinct soil type and its associated characteristics (Figure 2b). Almost the entire watershed contains one type of soil, Ontario Brookston clay (ONBKN), with a small portion of the watershed containing Ontario Caistor Soil (ONCTRA).

2.2.3. Land Use and Land Management Datasets

Land use and cropping pattern data for the watershed were collected from a windshield survey report provided by the ERCA and used to prepare the land-use map for the study (Figure 2c).

Land management data were obtained from the 2-year windshield survey report of farmers conducted by ERCA. However, the survey’s report contains the details of crops grown, tillage practices, amount of fertilizer used, application method, and harvesting efficiency for these operations. About 29% of farmers in this watershed participated in the survey, and the response was partial. But there was no information about the timings for various operations, and not much data was available from the local survey report given by ERCA. Therefore, data related to timing of operations such as planting, tillage operation, fertilizer application, harvesting, etc. were taken from Liu et al. [11], as shown in

Table 2. Crop management operations based on observations and time for the Wigle Creek watershed.

Crop	Corn	Soybean	Winter Wheat
Year	1	1	1	2
Tilling	25 October	12 May	24 October
Planting	2 May	15 May	25 October
Fertilizer-I	2 May	15 May	25 October
Fertilizer-II	29 May	—	—	25 April
Harvest and kill	20 October	12 October	—	20 July

. As for the fields/plots, no crop data record was found in the watershed. Thus, in the model, a standard corn–soybean rotation was assigned according to the information given by ERCA and the data used by Liu et al. [11] for southern Ontario.

2.2.4. Weather Datasets

The Jack Miner weather station, located within the Wigle Creek watershed, was used to procure the required climatic datasets. However, only precipitation data for a brief period (18 July 2016–12 July 2017) were available from the Jack Miner weather station. Therefore, other weather datasets (including missing daily precipitation, maximum and minimum daily temperatures, and other climatic parameters) were obtained from a nearby Environment Canada weather station (‘Harrow CDA’), which was used to fill in the missing data between the simulation periods (2012 and 2017), respectively.

2.3. SWAT Model Setup

In this study, a SWAT model setup was developed with combinations of land use, soil, slope, and weather datasets for the Wigle Creek watershed, accordingly. Wigle Creek watershed was delineated into 78 sub-watersheds based on 30 × 30 m DEM by burning the stream (streams that were prepared based on ground truth from field visits) with a threshold value of 10 hectares, as defined by SWAT recommended threshold values in a study by Neitsch et al. [31]. As a result, the total delineated area of the watershed is found to be 1949.13 ha. During the model development, single slopes (0–2%) were used for land use (78 polygons) and soil (2 groups) to delineate the HRUs (Hydrological Response Unit). Altogether, the model has generated 396 HRUs for the Wigle Creek watershed.

Here, before flow simulation, crop yields were compared to determine whether the simulated yields were acceptable or not based on the calibrated parameters, as endorsed by Moriasi et al. [32]. Crop yields simulated for 2012 exceeded observed yields, with winter wheat yields also being lower but closer to observed values (

Table 3. Average annual crop yield in the Wigle Creek watershed from 2012–2016.

Year	Yield (kg ha−1)
	Winter Wheat		Soybeans		Corn
	Obs.	Sim.	Obs.	Sim.	Obs.	Sim.
2012	4202	4038	2297	2863	6170	9985
2013	—	—	2745	2248	8812	6533
2014	—	—	2690	2434	—	—
2015	—	—	3194	2807	9492	9801
2016	—	—	2465	2742	—	—

Obs. = observed, Sim. = SWAT-simulated.

). However, from 2013 onward, simulated yields were closer to the previous year’s observed yields and significantly lower than the current year’s observed yields. Only soybean yields were available in 2014, with the simulated value being approximately 12% less than the observed yield (Table 3). The simulated yields for soybean in 2015 followed the same trend, with the simulated results displaying a roughly 10% reduction compared to the observed value. Meanwhile, corn yield displayed an increase of around 3% compared to the measured yield. In 2016, simulated soybean yields displayed an increase of 11% over the observed values (Table 3). However, the SWAT model has been actively used by a wide range of researchers worldwide [14,20,32]. Further details can be found in the SWAT documentation (http://swat.tamu.edu/documentation/ (accessed on12 March 2017)).

2.4. Evaluation of BMP Scenarios

The evaluation of BMP scenarios for the Wigle watershed was conducted using a calibrated SWAT model. During the calibration phase, two current BMPs (minimum-till (Curr_MT) and no-till (Curr_NT)), representing management practices applied in the watershed for the 2016–2017 season, were modeled. The coefficient of determination (R2), the percentage of bias (PBIAS), and Nash-Sutcliffe efficiency (NSE) suggested by Nash et al. [33] served as goodness-of-fit statistics to evaluate model prediction accuracy against observed values. The ranges of these statistics indicating excellent, good, poor, and unacceptable model performance used to assess calibration phase model performance were drawn from Moriasi et al. [32].

In the post-calibration phase, a total of six proposed BMP scenarios were developed by implementing BMPs onto a worst-case scenario no-BMP ‘control’ scenario (i.e., all existing BMPs currently implemented in the watershed are eliminated (Curr_-BMP)) as follows: (i) minimum tillage (PropMT), (ii) no-till (PropNT), (iii) retire croplands to pasture (PropRTP), (iv) retire croplands to forest (PropRTF), (v) cover crop after winter wheat (PropCCWW), and (vi) vegetative filter strips (PropVFS). While the PropRTP and PropRTF scenarios do not represent practical practices, they serve here as a basis for comparison with other possible scenarios. The Curr_-BMP scenario removes all the existing BMPs applied in the Wigle Creek watershed during the calibration phase of the modeling exercise (Figure 2c). However, for some BMPs (e.g., phosphorus management, variable rate fertilizer application, and in-field erosion control structures), only fragmentary data were available regarding their implementation and the location thereof. The subsequent section describes the various BMPs hypothetically implemented in the study (

Table 4. Existing and proposed BMPs implemented in the modeling phase.

BMP			SWAT Simulation
Description	Name	Type	Phase	BMPs Applied
Minimum tillage	CurrMT	Current	Calibration	Current practice
No-till	CurrNT	Current	Calibration	Current practice
No BMP	Curr-BMP	Control	Post-calibration	None
Minimum tillage	PropMT	Proposed	Post-calibration	Minimum tillage, and all croplands in watershed.
No-tillage	PropNT	Proposed	Post-calibration	No till, and all croplands in watershed.
Retire croplands to pasture	PropRTP	Proposed	Post-calibration	No agricultural management all crop land converted to pasture
Retire croplands to forest	PropRTF	Proposed	Post-calibration	No agricultural management, and all crop land converted to forest
Cover crop after winter wheat	PropCCWW	Proposed	Post-calibration	No agricultural management, oats planted late summer, and winter killed in January.
Vegetative Filter Strips	PropVFS	Proposed	Post-calibration	No agricultural management, and vegetative filter strips are applied at the edge of all agricultural fields

).

3. Challenges to Model Built-Up for the Flat Land Watershed (Wigle Creek Watershed)

The drainage network was prepared from the DEM, then corrected by ground truthing of overland flow pathways. This involved field work by visiting all the fields/plots in the watershed on several occasions to identify flow paths. In an initial step, an attempt was made to employ a 0.5 m × 0.5 m DEM to delineate the stream network in the SWAT model using the model’s default stream burning application. However, as automatic delineation of streams/flow paths from the DEM did not work and the resulting stream networks were mostly disconnected, in Figure 3a, these were manually prepared and modified based on ground truthing and then overlaid on a 30 m × 30 m DEM to detect the flow pattern of the drainage network. This limitation could be attributed to the fact that although an extremely fine-resolution DEM was used, certain man-made features like bridges, culverts, and field-side ditches might not have been properly captured, or errors occurred when the Lidar image was processed. To address this issue, in the second step, an externally derived drainage network (obtained from field visits) overlayed upon the stream network created with the 0.5 m × 0.5 m DEM was used for the modeling exercise. Henceforth, changes were detected in the flow path pattern, thereby demonstrating that the generation of drainage patterns after burning the streams was much more accurate for analysis with high resolution (0.5 m Hydro-Enforced DEM).

However, there were some flow issues noted in the northeastern portion of the watershed. These were attributed to an improper representation of the drainage network near the highway, even after delineation of the stream network with the 0.5 m Hydro-Enforced DEM (Figure 3b). In the third step, a 30 m × 30 m resolution DEM was used for stream delineation. The low resolution of the DEM resulted in another issue: the drainage network in the middle and upper parts of the watershed disappeared, draining rather towards the west (Cedar Creek watershed) (Figure 3c). Therefore, the automatically generated drainage network was found to be unsuitable for modeling purposes. Accordingly, in an ultimate step, the model was set up with an externally derived and corrected drainage network (obtained from field visits) overlayed upon the stream network created with the 30 m × 30 m DEM. The drainage network pattern thus obtained was found to mimic the observed drainage network realistically and accurately in the middle part of the watershed (Figure 3d). Therefore, the stream network thus obtained was finally selected for the SWAT model setup.

4. Results and Discussion

4.1. Calibration of Flow, Sediment, and Phosphorus

To accurately simulate watershed processes, the SWAT model was built with state-of-the-art components, as suggested by Santhi et al. [34]. Nonetheless, being a physically based watershed model, SWAT has several empirical components. Certain variables, such as the curve number and the cover and management factor, are not physically fixed, according to a study by Santhi et al. [34]. Therefore, the SWAT was first calibrated with observed data from Kalin and Hantush [35].

Only 47 days of observed data were available for 2016 and 2017. Accordingly, the SWAT model was simulated for a period between 2012 and 2017, with 4 years (2012–2015) as a warm-up period and 2 years (2016–2017) as a calibration period. In addition to daily climate data, continuous flow and a significant amount of water quality data were required to properly calibrate the watershed and water quality models. As far as we know, there are no recommendations available related to the adequacy of these data. In general, it is recommended to have five or more years of flow and water quality data to properly capture climate variability and calibrate the model. However, as mentioned above, only data obtained through spot measurements of flow and water quality were available. To address this challenge, model calibration and validation efforts were focused on the days for which data were available. Specific attention was given to model performance for flow, sediment, and phosphorus on the days with observed data. Hence, the list of various parameters in the model that were changed during the SWAT model calibration with existing BMPs is shown in Table S1.

4.1.1. Calibration of the Flow

Daily observed and simulated flows at the Wigle Creek outlet were used for SWAT model calibration. To procure statistically accepted values for calibration, a range of sensitive parameters were provided from SWAT-CUP (CUP: Calibration and Uncertainty Programs), where 1000 simulations were run with a uniform prior distribution of the parameters based on the upper and lower limits (Table S1). Also, if required, more runs were simulated as per the ‘new parameter sets’, as recommended by the SUFI-2 (sequential uncertainty fitting) algorithm.

Figure 4 provides a graphical comparison of streamflow for the calibration period during the period of 2016–2017. However, to assess the agreement between observed and simulated datasets, a scatter plot is drawn, as shown in Figure 5. The values of NSE and PBIAS of 0.52 and 6.71 obtained for the total stream flow were within their respective acceptable ranges, according to Moriasi et al. [32], indicating that the model calibration was adequate (

Table 5. Model performance for flow simulation at the watershed outlet.

Calibrated Parameter		Sampling			Model Accuracy Statistics
Type	Units	Start	End	Number	PBIAS (%)	R2	Daily NSE
Flow	m3 s−1	January 2016	May 2017	47	6.71	0.56	0.52
Sediment concentration	mg L−1	January 2016	December 2017	122	−1.19	0.13	0.13
Sediment load (calculated)	ton	January 2016	December 2017	47	−15.94	0.31	0.30
Phosphorus (calculated)	kg	January 2016	May 2017	47	−82.57	0.17	0.08

).

4.1.2. Sediment Load Calibration

Sediment load calibration was performed after flow calibration. Like flow calibration, SWAT-CUP was used to calibrate sediment loads.

High flow periods were given greater consideration during the calibration process since they produced large sediment loads. The time series graphical comparison of sediment for the calibration period (Figure 6) suggests that, based on the PBIAS, model performance is satisfactory for sediment loads and concentrations, according to Moriasi et al. [32]. However, NSE and R2 values indicate the model’s performance is inadequate. Notably, the model performed better in predicting sediment loads than sediment concentration (Table 5). Figure 7 shows the scatter plot between observed and simulated sediment loads in the Wigle Creek watershed.

4.1.3. Phosphorus Load Calibration

After streamflow and sediment calibration, the SWAT model was calibrated for phosphorus using SWAT-CUP (Figure 8). Modeled phosphorus loads suggest that SWAT prediction during the calibration period was not satisfactory (Table 5). The value for NSE is below 0.35, and the PBIAS is well above/below ±30%, both of which are not satisfactory based on the criteria suggested by Moriasi et al. [32]. However, to evaluate the discrepancy between the observed and simulated sets of (event) data, a scatter plot is studied, as represented in Figure 9.

Thus, the model limits indicate that most of the data used in the calibration of the model was instantaneous data, with some days having an average of many samples. Though some days have average daily data, the times at which observations were taken were not spread uniformly throughout the day. This ultimately skews the observed data, which affects calibration due to low data quality. A comparison was made between the continuous simulated results and the instantaneous observation. As a result, less than suitable results are produced from the simulation after calibration, which does not accurately represent the watershed behavior. However, Wigle Creek watershed average annual flow, sediment, and phosphorus loads at the outlet of the watershed are given in

Table 6. Average annual flow, sediment, and phosphorus of the watershed model.

Year	Average Flow (m3 s−1)	Sediment (ton)	Phosphorus (kg)
2016	0.157	527.5	239.2
2017	0.222	746.5	413.9

.

4.2. Effectiveness of Existing BMPs

The calibrated SWAT model was used to assess the effectiveness of BMPs presently implemented on the watershed, namely Curr_MT and Curr_NT. To estimate their effectiveness, existing BMPs were removed from the calibrated model (Curr_-BMP), and phosphorus loads at the watershed outlet were simulated (

Table 7. Effectiveness of existing BMPs on reducing flow and P loading at the watershed outlet for GLASI land management BMPs during 2016–2017 (units are percentage).

Existing BMPs Reductions (%)
Year	Season	Flow	TSS	P Organic	P Mineral	P Total
2016	Non-Growing	−0.47	−1.01	−4.69	0.70	−2.57
	Growing	2.13	3.74	11.55	100.22	43.56
	Year	0.25	−0.11	−4.07	3.66	−1.05
2017	Non-Growing	0.16	0.36	−11.69	9.49	−4.29
	Growing	0.99	1.23	0.00	−0.04	−0.03
	Year	0.41	0.58	−11.12	8.45	−4.01

Negative numbers represent a reduction whereas positive numbers are an increase.

). Results were analyzed on a yearly and seasonal (non-growing and growing season) basis. The non-growing season spans from November to April (when most of the crops are not growing), and the growing season begins in May and ends in October. In this study, all scenarios were compared to the Curr_-BMP scenario (Figure 10). Table 7 and Figure 10 demonstrate the change in phosphorus loads for the existing BMPs. The current BMPs implemented in the Wigle Creek watershed show an annual phosphorus reduction of 1.05% and 4.01% in 2016 and 2017, respectively. The phosphorus reductions are higher in the non-growing season. The reduction in total phosphorus is primarily due to the reduction of organic phosphorus, with a smaller increase in mineral phosphorus. Also, current BMP implementation demonstrates no significant change in either flow or sediment yield at the watershed outlet.

4.3. Effectiveness of Proposed BMPs

The calibrated model was used to evaluate the effectiveness of potential future BMPs (PropMT, PropNT, PropRTP, PropRTF, PropCCWW, and PropVFS) (Table 4). These BMPs were individually applied in all fields, and the total phosphorus loading in each scenario was compared with the Curr_-BMP condition. The BMP effectiveness was also computed at various temporal scales (e.g., annual, conventional seasons, and growing/non-growing seasons).

4.3.1. The Minimum Tillage Scenarios

The PropMT BMP includes the modification of tillage operations in all agricultural fields in the watershed area to “min-till.” The results of this input setup were compared with the Curr_-BMP (fall moldboard plow, spring cultivate) model input scenario. Phosphorus reduction for the PropMT scenario occurred as expected, with annual load reductions exceeding 40% and 30% for 2016 and 2017, respectively (

Table 8. Effectiveness of proposed BMPs compared to Curr_-BMP (units are in percentage).

BMP Imposed	Year	Season	$\mathbf{\%}\mathsf{\Delta } \mathbf{C}\mathbf{o}\mathbf{m}\mathbf{p}\mathbf{a}\mathbf{r}\mathbf{e}\mathbf{d} \mathbf{t}\mathbf{o} {\mathbf{C}\mathbf{u}\mathbf{r}\mathbf{r}}_{\mathbf{-BMP}}$
			Flow	TSS	Phosphorus
					Organic	Mineral	Total
PropMT	2016	Non-growing	1.99	2.71	−50.17	−25.9	−40.59
		Growing	0.37	−0.17	−40.61	−28.32	−36.17
		Year	1.58	2.16	−49.82	−25.97	−40.43
	2017	Non-growing	0.54	0.53	−46.18	−12.25	−34.35
		Growing	0.5	0.57	−3.55	−3.87	−3.74
		Year	0.54	0.55	−44.08	−11.38	−32.21
PropNT	2016	Non-growing	2.24	4.00	−65.87	−42.29	−56.55
		Growing	0.34	−0.02	−64.01	−48.41	−58.38
		Year	1.77	3.22	−65.79	−42.47	−56.61
	2017	Non-growing	0.46	1.30	−64.24	−32.54	−53.19
		Growing	0.41	0.48	−3.66	−4.11	−3.94
		Year	0.45	1.11	−61.27	−29.54	−49.75
PropRTP	2016	Non-growing	−23.73	−82.45	−80.59	−48.99	−68.1
		Growing	19.81	−65.79	−63.9	−17.32	−47.08
		Year	−12.23	−79.26	−80.01	−48.03	−67.41
	2017	Non-growing	−31.77	−84.03	−80.24	−35.01	−64.48
		Growing	36.41	−67.42	−9.97	10.22	1.24
		Year	−12.48	−80.15	−76.81	−30.14	−59.87
PropRTF	2016	Non-growing	−27.53	−83.61	−83.43	−63.89	−75.71
		Growing	25.39	−63.11	−68.53	−46.45	−60.55
		Year	−13.62	−79.7	−82.91	−63.38	−75.21
	2017	Non-growing	−37.03	−85.77	−87.85	−62.72	−79.09
		Growing	40.94	−66.58	−8.59	10.12	1.8
		Year	−14.86	−81.29	−83.98	−54.9	−73.42
PropCCWW	2016	Non-growing	−0.36	−0.2	−2.41	−2.42	−2.41
		Growing	−3.22	−4.01	−0.49	0.08	−0.28
		Year	−1.14	−0.93	−2.28	−2.33	−2.34
	2017	Non-growing	−1.31	−2.24	−3.56	−3.14	−3.41
		Growing	−0.8	−0.71	−0.01	0.01	0
		Year	−1.17	−1.86	−3.38	−2.79	−3.17
PropVFS	2016	Non-growing	0	−73.54	−45.56	−33.01	−40.6
		Growing	0	−73.57	−35.48	−23.43	−31.14
		Year	0	−73.55	−45.2	−32.72	−40.3
	2017	Non-growing	0	−73.73	−43.26	−30.77	−38.9
		Growing	0	−73.46	−1.91	−1.91	−1.91
		Year	0	−73.66	−41.24	−27.68	−36.31

). This was to be expected given the reduced surface runoff, which would decrease the sediment transport rate and consequently the organic phosphorus transport to the stream. This would reduce the phosphorus loads during the growing season. Reductions in total phosphorus loads can be largely attributed to reductions in organic phosphorus loads. However, the reductions in mineral phosphorus under this BMP were also significant and contributed to the overall reduction of phosphorus loads. Furthermore, while sediment yield fluctuates with flow, its increase could be considered insignificant (Table 8). Results on a seasonal basis indicate phosphorus load reductions of 40, 23, 7, and 33 percent during the winter, spring, summer, and fall seasons, respectively (Figure 11). The PropMT scenario resulted in 11% and 37% reductions in total phosphorus loads during the growing and non-growing seasons, respectively.

4.3.2. The No Tillage Scenarios

Like the PropMT BMP, for the PropNT BMP, tillage operations in all agricultural fields in the watershed were changed to no-till. No-till practices over the entire watershed could provide a 50% annual reduction in phosphorus loads (Table 8). However, the 2017 growing season, when phosphorus loads were relatively small compared to the non-growing season, did not show any major reduction in phosphorus. Sediment yields under this scenario were greater than the previous ones (e.g., PropMT) but are unlikely to have a significant impact on the overall phosphorus yield. Similar to the trend when PropMT was applied throughout the entire watershed, PropNT provided relatively high reductions in phosphorus yield at the watershed outlet. The reductions in annual total phosphorus load were 56% and 50% in 2016 and 2017, respectively (Figure 12). Undisturbed soils under no-till show enhanced biological activity, resulting in an increase in infiltration and a concomitant decrease in surface runoff. Reduced surface runoff and relatively undisturbed soil would obviously decrease soil erosion rates and organic phosphorus transport within the watershed. The reduction in total phosphorus loss was more evident in the simulated non-growing season (54%) than the growing season (15%). However, results split up according to the seasonal total percentage indicate that a maximum reduction of 50% was observed during winter, followed by spring (23%), fall (58%), and summer (8%) (Figure 12).

4.3.3. Retiring Land Scenarios

The retiring of land scenarios is based on the reintroduction of native flora (trees, shrubs, grasses, etc.) to the area on agricultural land. This would be useful for land that does not produce or yield much of a crop, either due to erosion or soil degradation. As the data for the Wigle Creek watershed was lacking to determine the fields that are not profitable, the operations were given to the entire watershed as a best-case scenario, even though the operation of such a BMP is unlikely.

Retire Pasture Scenario

Phosphorus reduction under the PropRTP BMP (conversion of all the agriculture land in the watershed into pasture) was significantly greater than that achieved by the other BMPs simulated, suggesting the outcome of a near-best-case scenario, although the implementation of such a BMP is unrealistic. The average annual reduction in phosphorus losses was above 60%, with higher reductions occurring during the non-growing season (Table 8). Retiring land where no crops are grown, no land is tilled, and no fertilizer is applied, although impractical, was found to be a best-case scenario. As expected, phosphorus loads were significantly reduced at both seasonal and annual scales compared to PropMT or PropNT operations. Flow also shows a significant reduction during the simulation period, with a greater reduction occurring during the non-growing season. Though overall yearly average flow was reduced, there was an increase in flow during the growing period for each year. This scenario displays tremendous reductions in sediment at the outlet for all years. For flow and phosphorus, the reductions occur primarily during the non-growing season, when flow and phosphorus are at their peak. Figure 13 represents the results of the conventional four seasons along with the growing and non-growing seasons. The effectiveness of this BMP was much less pronounced during the summer than during the winter and fall seasons. Further, it was more effective during the non-growing season months (66%) than during the growing season months (9%).

Retire Forest Scenario

The PropRTF BMP (conversion of all the agriculture land in the watershed into forest), an utterly utopian scenario, shows a further reduction in phosphorus loads compared to the PropRTP BMP. The average annual reduction in phosphorus loads compared to Curr_-BMP was estimated to exceed 70%. Similar results were found for the non-growing seasons (Table 8 and Figure 14). The flow reduction was significant but showed only a minimal increase compared to the PropRTP BMP. However, total phosphorus reduction results indicate that mineral phosphorus reductions were greater than those for organic phosphorus. This was true for sediment yield as well. On a seasonal basis, average annual phosphorus reduction was 78% during the non-growing season and 11% for the growing season (Figure 14).

4.3.4. Cover Crop after Winter Wheat Scenarios

The PropCCWW BMP includes an evaluation of establishing a cover crop after winter wheat in years where winter wheat is grown in rotation. For the PropCCWW BMP, oat was chosen as the cover crop based on information provided by ERCA. In this scenario for the Cover Crop BMP, the oats were planted after winter wheat was harvested and applied to all the fields in the Wigle Creek watershed. Compared to the conservation tillage scenarios (PropNT and PropMT), cover crops showed a minimal decrease in flow and sediment yield. An average annual reduction of 3% was observed in annual total phosphorus loads. This BMP demonstrates some phosphorus reduction because of cover crop growth prior to winter, which killed the crop, and the residue left over from winter or snowfall months. A 3% decrease in total phosphorus loads during the non-growing seasons was estimated, with no change during the growing season. Most of the reductions in phosphorus load were achieved during the winter and fall seasons, at rates of 4% and 1%, respectively (Figure 15).

4.3.5. Vegetative Filter Strips Scenarios

Under the PropVFS BMP, vegetative filter strips (VFS) were applied along the edge of a field. Based upon the parameters used, for every unit area of filter strip, there were 40 units of crop-seeded area. As such, the area and width of the filter strips varied based on the size of the HRU to which they were applied. Table 8 shows that VFS resulted in a 38% reduction in total annual phosphorus loads at the watershed outlet. Reductions in both mineral and organic phosphorus loads were similar to those under other BMPs. During the winter, spring, summer, and fall seasons, reductions in total phosphorus of 45%, 12%, 4%, and 40%, respectively, were observed (Figure 16). Further, decreases of 8% and 40% were observed in total phosphorus load during the growing season and non-growing season, respectively.

5. Conclusions

In this study, the SWAT model simulated potential future Best Management Practices (BMPs) for the phosphorus reduction in a flat agricultural watershed in Ontario. The calibration of the model focused on streamflow, sediment, and total phosphorus loads at the watershed outlet. Prior to flow calibration, simulated crop yields were compared with observed yields to ensure the average annual simulated yields fell within an acceptable range. This study also addresses the difficulties in modeling a watershed with flat terrain and limited observed data for flow and water quality. Various combinations of DEM resolutions were tested, revealing that a 30 m resolution DEM, incorporating modified streams from ground truthing, was effective for SWAT model drainage framework development. The study emphasizes that a finer DEM resolution may not always enhance drainage network delineation in flat terrains. To overcome limited calibration data, the study focused on model performance during days with available observed data, underscoring the potential of grab sampling for calibration. However, the study highlights the importance of collecting sufficient flow and water quality data for improved modeling confidence. Challenges also arose from the scarcity of land management information, including fertilizer use in the watershed.

Projected phosphorus load reductions were examined by converting agricultural fields to pasture or forest land use, but this is impractical due to agriculture’s vital role in the regional economy. Applying PropMT or PropNT BMPs across the watershed holds potential for substantial phosphorus reduction, but it may also introduce challenges such as herbicide dependence for weed control and the need for ongoing land maintenance investment. Cover crops, such as post-winter wheat, have a limited impact on phosphorus reduction. Vegetative filter strips along field edges were found to offer significant phosphorus reduction, but implementing this BMP may reduce crop-growing land, causing economic losses to producers unless incentives or compensation payments are established.

The results of this study provide preliminary information about the effectiveness of agricultural BMPs to address high phosphorus losses within a flat watershed. This study reveals the importance and challenges associated with proper watershed delineation, accurate identification of flow paths, the importance of the collection of crop management data, and the efforts required to collect continuous flow and water quality data to support watershed modeling. However, farmers frequently hesitate to adopt BMPs due to potential conflicts with their crop needs or local environmental conditions. Economic constraints further impede BMP adoption, as some practices involve initial investments that farmers find challenging to afford. The effectiveness of BMPs relies on widespread and consistent adoption, highlighting the crucial role of education and outreach. Addressing these challenges demands a nuanced approach, considering the intricacies of individual farming systems and striving for a balance between environmental sustainability and economic feasibility.