Evaluation of BMPs in Flatland Watershed with Pumped Outlet

Abstract

The effectiveness of existing and potential best management practices (BMPs) to cropped lands in the Jeannette Creek watershed (Thames River basin, Ontario, Canada) in reducing P loads at its pumped outlets was assessed using the Soil and Water Assessment Tool (SWAT). Existing BMPs consisted of banded, incorporated, and variable phosphorus (P)-rate application, conservation tillage, cover crops, and vegetative buffer strips. Potential BMPs consisted of banded P application, no-till, and a cover crop following winter wheat. Two separately delineated sub-watersheds, J1 and J2, characterized by a flat topography and distinct pumped outlets, were selected for analysis. Despite challenges in delineation, the SWAT model was successfully set up to assess the impact of BMPs in reducing P loads in these sub-watersheds. Each BMP was systematically removed, and the resulting simulated P loads were compared with the baseline scenario. Compared to cover crops or vegetative buffer strips, the implementation of conservation tillage and no-till, along with altering the mode of P application, offered superior effectiveness in reducing the P load. On average, the annual reduction in total P (P_tot) loads under existing BMPs was 9.2% in J1 and 11.3% in J2, whereas, under potential BMPs, this reduction exceeded 60% in both watersheds.

1. Introduction

An unmatched treasure of Ontario, Canada, the Great Lakes bear one-fifth of the world’s fresh water [1]. The Great Lakes’ health, particularly that of Lake Erie, is under a serious threat due to increased levels of harmful contaminants and escalating levels of phosphorus [2,3]. Presenting a major environmental issue, elevated levels of phosphorus (P) in their waters has been linked to the rapid growth (blooming) of blue-green algae (cyanobacteria) [4,5]. Phosphorus loading has been reported to be associated with runoff, and to originate primarily from agricultural areas [6,7,8].

Inadequate agronomic management practices result in serious environmental threats, especially in terms of promoting blue-green algal blooms in downstream lakes (e.g., Lake Erie) [5,9,10]. Various agronomic best management practices (BMPs) have been reported to be effective in limiting P loading [11,12,13,14,15,16]. While the implementation of appropriate BMPs can reduce P_tot losses [17], the optimal manner to evaluate the effectiveness and viability of BMPs within a given watershed remain poorly documented [18]. Moreover, the application of non-integrated practices within a watershed can generate the majority of the P load the watershed contributes to downstream bodies of water [19,20]. Accordingly, collaborative efforts between Canada and the United States have targeted the reduction of P loads reaching Lake Erie. A recent initiative, the Great Lakes Agricultural Stewardship Initiative (GLASI), has placed a major emphasis on P load reduction [21].

Designed to simulate complex hydrological processes and evaluate BMPs, the Soil and Water Assessment Tool (SWAT) model [22,23] has seldom been applied to simulate watersheds with a flat topography [24,25] and never those combining a flat topography and pumped outlets. Donmez et al. [26] used the SWAT model to simulate a flat watershed. To address the associated conditions, they subdivided their watershed into sub-basins based on topography, then further split the sub-basins into hydrologic response units (HRUs). Accordingly, the SWAT model was selected for use in this study given its ability to simulate different watershed topologies and scales [27].

In this study, the Jeannette Creek watershed, acknowledged as a priority watershed under GLASI, displays distinctive features, including a flat topography, Brookston clay soils, cash cropping, and a tile drainage management system. The flat topography hinders natural outflow of water, compelling the use of a sump pump structure to pump water out of the watershed and into the stream. Recognizing the potential of using the SWAT model to effect a time-saving analysis of the impact of BMPs on P load reduction in a flatland watershed with a pumped outlet, the following objectives were targeted: (i) build a SWAT model using hydro-meteorological datasets for a watershed with a pumped outlet, (ii) calibrate the SWAT model(s) with available but limited monitoring data, (iii) investigate the effectiveness of existing and potential BMPs in reducing P loads.

2. Materials and Methods

2.1. Description of Study Area

The Jeannette Creek watershed, located near Chatham, ON, has priority status in the GLASI program. Covering 37,600 ha, with 1866.92 ha studied, it comprises two sub-watersheds, J1 (914.34 ha) and J2 (952.58 ha), draining into Lake Erie. Predominantly flat with Clyde soils, 97% is dedicated to agriculture (corn and soybeans). The flat topography requires pumped outlets to discharge water into the Thames River. Because of flat topography, the watercourse outlets are below the creek level, and the watershed outlet is not a natural one. Pumped outlet systems are employed to discharge water from the Jeannette Creek watersheds into the Thames River. This summary provides a concise overview of the J1 and J2 watersheds, detailing their location, extent, topography, soil, land use, and water management practices. (Figure 1).

Data required for simulating the Jeannette Creek Watershed were obtained from the Lower Thames Valley Conservation Authority (LTVCA), Ontario Ministry of Agriculture, Food and Rural Affairs (OMAFRA, https://omafra.gov.on.ca/English) (accessed on 5 May 2017), and Environment Canada Climate Change (ECCC). To set up the SWAT model, various datasets were required, including climate data, spatial data (viz., elevation data, soil data, land use, and crop and management practices data) (

Table 1. Data available for Jeannette Creek watershed.

Data	Type	Source	Description
Digital Elevation Model	Raster file	LTVCA, 2017	0.5 m × 0.5 m LIDAR image
Precipitation, Relative Humidity, Solar Radiation, Maximum and Minimum Temperature	Excel	LTVCA, 2017; ECCC, 2017	Obtained from light house Cove, Merlinb (LTVCA stations) and Chatham (ECCC station)
Soil	Shape file	OMAFRA soils	Soil Landscapes of Canada (SLC) version 3.2
Land Use	Shape file	LTVCA, 2017	Plot-wise crop date spatial map
Stream Network	Shape file	University of Guelph (UOG), 2017	Prepared based on ground truth survey
Land Management	Shape file and Excel	LTVCA, 2018	5-year farmer survey report
Stream Flow	Excel	LTVCA, 2017	Instantaneous data
Sediment	Excel	LTVCA, 2017	Instantaneous concentration data
Phosphorus	Excel	LTVCA, 2017	Instantaneous concentration data

). To validate the model for surface water hydrology required observed streamflow, sediment, and P load datasets for nutrient modelling (Table 1).

2.2. Model Input Datasets

The precision of hydrological model predictions depends on the input of variables that accurately characterize watershed features. Accordingly, to successfully run the SWAT model, a wide set of data is vital. The fact that the J1 and J2 watersheds were flat, and each equipped with a pumping outlet, presented a challenge in developing a representative SWAT model. The data needed to implement the SWAT model includes a Digital Elevation Model (DEM), climate data (precipitation and temperature), soil, and land use, along with flow, sediment, and P load data at the watershed pump outlet [28]. Accordingly, the input data were analyzed and processed prior to their use as a model input, then used to create various thematic maps and their associated databases.

2.2.1. DEM

In lieu of the Jeannette Creek watershed (J1 and J2), a 0.50 m × 0.50 m (resolution) DEM was developed from contour data. Elevation of the full (J1 + J2) watershed varies from 174 m to 186 m above mean sea level (AMSL), with the J1 sub-watershed ranging from 173 m to 179 m AMSL south-east to north-west, and the J2 watershed ranging from 174 m to 186 m AMSL in the opposite direction (Figure 2). To validate the projection of horizontal and vertical units of the procedures, DEM properties were set up for the J1 and J2 watersheds.

2.2.2. Climate Data

Meteorological data are essential as inputs in setting up the SWAT model. Climate data for seven years (2011 to 2017) were obtained from two sources: (i) five years (2011–2015) from nearby Environment Canada Climate Change stations and (ii) the remaining two years (2016–2017) from the Lower Thames Valley Conservation Authority (LTVCA) for the Jeannette Creek watershed. Two-year data (2016 to 2017) of precipitation and temperature for J1 and J2 were collected from the Lighthouse Cove and Merlinb stations, respectively.

Preliminary analysis showed roughly 55% of observed precipitation occurred during the crop growing season (April to October) and 45% fell as rain/snow over the non-growing seasons. Over the 2011–2017 period, some precipitation data values were missing; these were generated using SWAT’s weather simulator. Over the period of 2011–2017, the average (±standard deviation) annual observed precipitation on the J1 and J2 watersheds was 792 ± 183 mm and 797.6 ± 173 mm, respectively, indicating that these watersheds showed substantial variation (Coeff. of Variation > 20%) in total annual precipitation during the period of study. Over the same period, the maximum annual precipitation was 827 mm on the J1 watershed and 823 mm on the J2 watershed, both in 2013. In contrast, 2012 was a relatively dry year with annual precipitation of 636 mm and 652 mm for J1 and J2, respectively.

The average annual temperature (T_mean) for J1 was 10.1 ± 1.2 °C [range from 8.2 °C (2014) to 11.9 °C (2017)] and, for J2, 9.9 ± 1.1 °C [range from 8.1 °C (2014) to 11.4 °C (2017)]. A summary was established of monthly average precipitation and T_mean, for the Jeannette Creek watershed from January 2016 to December 2017

2.2.3. Soils

To prepare the Jeannette Creek (J1 and J2) watersheds’ soil input data for the SWAT model, soil data was drawn from Soil Landscapes of Canada (SLC) version 3.2 (Soil Landscapes of Canada Working Group, 2007, https://sis.agr.gc.ca/cansis/nsdb/slc/index.html, accessed on 20 May 2018) and used as a field database. The SLC data encompasses a soil map of Canada as well as key features of soils for the whole country. The SLC exists at a scale of 1:1 million (~1 km spatial resolution), with each polygon on the map describing a distinct type of soil and its associated characteristics. Figure 3 displays the soil coverage map of the J1 and J2 portions of the Jeannette Creek watershed.

Soil properties are important factors in controlling infiltration and soil water movement and play a key role in surface runoff, groundwater recharge, evapotranspiration, soil erosion, and the transport of chemicals. For the Jeannette Creek watersheds (J1 and J2), soil map and data were downloaded from Soil Landscapes of Canada (SLC) and, using ArcGIS 10.3.1 software, subsequently integrated to create the soil attribute datasets necessary to build up the SWAT model. J1 watershed major type soil is Ontario Clyde (ONCYDA) and J2 watershed is Ontario Brookston clay (ONBKN).

2.2.4. Land Use and Land Cover (LULC)

To create a land use spatial map for both watersheds, a plot-wise crop classification was undertaken based on ground-truth cropping pattern data drawn from a 5-year crop rotation survey report provided by the LTVCA. On this basis, the J1 watershed was classified into 49 land use polygons/specific classes, whereas the J2 watershed was categorized into 54 land use classes.

When modelling agricultural watersheds with the SWAT model, accurate land management data inputs play a significant role in ensuring the accurate simulation of runoff, sediment, and nutrient loads. Land management information includes date/month of tillage operations, sowing/planting dates, type and rates of fertilizer application, harvest/kill date, etc. Field-by-field information from the farmers of the cropped areas in the J1 and J2 (Jeannette Creek) watersheds (Table S1) was similarly drawn from a 5-year farmer land management survey report provided by the LTVCA. To represent the land management data on a field-wise basis, land use maps were again classified, such that each field was assigned a particular land use code (Figure 4).

For the J1 watershed, data collection between 2013 and 2016 was only carried out for 5 of the 49 specified fields (viz., J161, J162, J163, J165, J175). For the J2 watershed, over the course of five years, the LTVCA collected detailed information for 13 fields (specifically, J204, J2E4, J205, J207, J212, J213, J224, J232, J233, J238, J239, J241, J244) out of a total of 54. Based on the data collected, a seven-year rotation database was created to serve as input for the SWAT model. Land management data for the remaining fields in both the watersheds was taken from the land-use data description files (windshield survey) of 2016–2017, a product of the LTVCA. Information collected from the land management survey was extrapolated to fields with missing data using the adjacency principle. This approach enabled the creation and synthesis of wide-ranging land management datasets. Preparation of land management datasets involved 4 steps: (i) planting details for 3 growing seasons, (ii) fertilizer application details, (iii) harvest/kill details, and (iv) tillage operation details (Table S1).

3. SWAT Model Built-up and BMP Implementation

A process-based watershed model, SWAT is widely recognized for its capacity to accurately assess land management impacts on water, sediment, nutrient (phosphorus), and associated agrichemicals arising in a watershed under conditions of changing soil, land use, and management. The model performs semi-distributed, continuous-time hydrological simulations at a daily time step [29,30,31].

3.1. Delineation of Flow Paths

The initial challenge in establishing the model centered around delineating flow paths. Traditionally, watershed models require a natural free-flow outlet, but the Jeannette Creek J1 and J2 watersheds presented an obstacle with a dam/wall and a pumped outlet system. This posed challenges, such as defining streams and sub-watersheds, accommodating flat topography, incorporating pump stations, and assuming a pumped outlet. Similar challenges were encountered by Donmez et al. [26]. To address these issues, the study independently set up the two sub-watersheds (J1 and J2) within the larger Jeannette Creek watershed. The delineation of natural streams was performed at different resolutions, but due to limitations, pre-defined streams and sub-watersheds were applied to ensure comprehensive coverage. The methodology aligns with Donmez et al.’s approach, where pre-defined streams were crucial for accurately delineating the stream network in the model. This process was repeated for both the J1 and J2 watersheds to overcome the complexities arising from flat topography and multiple outlets.

For instance, Figure 5a,b show that the model was clearly unable to delineate the pattern of streams with the 30 m × 30 m DEM. In the case of Figure 5c, the natural streams were burned based on a 0.5 m × 0.5 m Hydro-Enforced DEM, but the model remained unable to delineate some portions of the watershed, as highlighted by the red circles. Figure 5d displays an appropriate delineation of the stream network, with pre-defined streams within sub-watersheds, allowing the routing of streams node to node. Thus, the pre-defined approach was applied to delineate the watershed.

3.2. Incorporation a Pumping Outlet in SWAT Model

The J1 and J2 watersheds function as pump catchments, with Dauphin in J1 and Dreary and Boudreau in J2 serving as pumped outlets. Each outlet is equipped with a sump and pump system. To ensure efficient operational control of the pumps, these pumps were designed to operate through an automated float mechanism. In the case of J2 Dreary, the pump outlet was functional throughout, whereas the Boudreau pump station was not functional due a misfunctioning automated float sensor. Therefore, among the three stations monitored under the GLASI program, two stations provided the majority of data for this study. The system could not be simulated as a continuous flow system, as, for these watersheds, sediment and P load data were not available on the channel. Pumped discharge, flow, sediment, and P load data were used during the model calibration.

Within the SWAT environment, pump outlet chambers were represented as reservoirs positioned at ground level. The implementation of a reservoir in SWAT requires the reservoir surface area and reservoir volume to be specified for different stages (i.e., emergency and principal spillway) [32,33]. To address this difficulty, the maximum surface area and corresponding volume of the reservoir were evaluated using the ArcGIS spatial analyst tool. To compute streamflow, sediment, and P loads at the reservoir (i.e., the pump) outlet, we assumed that the pumped discharge corresponded to the depth of the reservoir. Therefore, the amount of flow, sediment, and P discharged by the pump was taken as equal to the flow, sediment, and nutrients entering the channel from the watershed. In technical terms, we calibrated the J1 and J2 watersheds by estimating the reservoir area through a combination of Google Earth analysis and ground truthing, employing ArcGIS 10.3.1 software. Later, applying sensor depth (m) obtained from level logger data, we estimated the volume of the reservoir (m³) as the product of depth (m) and reservoir area (m²). This information enabled us to estimate the flow rate by dividing the volume by the time duration (m³ sec⁻¹) for a specific day.

3.3. Delineation of Watersheds and Sub-Watersheds, and Generation of HRUs

The J1 watershed was delineated into 24 sub-watersheds based on stream networks. To achieve this delineation, the ground truth flow paths (streams) were overlaid on a 0.5 m × 0.5 m DEM through the SWAT model’s burning process. The total delineated area of the watershed for the Dauphin pump outlet was 900.2 ha (Figure 6, J1). Single slope, land use (containing 186 polygon data), and soil (3 soil groups, detailed data) input maps were then overlain to generate HRUs. With a combination of all these inputs, the model generated 200 HRUs for the J1 watershed. In contrast, the J2 watershed was delineated into 32 sub-watersheds based on a 0.5 m × 0.5 m DEM by burning the ground-truthed flow paths (streams). This included three outlets, two of which were pumping outlets and a third that operated in a particular sub-watershed (Field ID 30) (Figure 6, J2).

3.4. Model Calibration

For the J1 and J2 portions of the Jeannette Creek watershed, the SWAT model was calibrated with data from the period of 1 January 2011 to 31 December 2017, followed by a 5-year (2011–2015) warmup period and 2 years (2016–2017) for the simulation of streamflow, sediment, and ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ loads at the watersheds’ pumped outlets. Streamflow calibration focused on improving model performance for the flow measured by monitoring stations at pump outlets. Various parameters related to flow hydrology (i.e., snow and snowmelt-related parameters, groundwater parameters, and relative parameters like CN2, SOL_K, SOL_ALB, SOL_AWC, etc.) were considered. Special emphasis was given to the parameters that may have spatial patterns that differ from HRU to HRU at the field level. To analyse the model’s sensitivity, parameter ranges were assigned to the SWAT CUP function to run for 500 simulations based on the upper and lower limits. Then, based on the Nash–Sutcliffe model efficiency coefficient (NSE), along with p- and r-stats pertaining to the uncertainty band, a final range of parameter values was obtained for both watersheds (J1 and J2).

To simulate the sediment at the pump outlet, the SWAT model implemented for the J1 and J2 (Jeannette Creek) watersheds was calibrated for the period of 1 January 2016 to 31 December 2017. The availability of sediment data was a challenge as only 71 and 85 sediment concentration sampling data (on daily basis) were available for the J1 and J2 pump outlets, respectively. The available data were used to accomplish this daily basis calibration through manual and auto SWAT Calibration and Uncertainty Procedures (SWAT-CUP) by comparing the simulated sediment loads to the measured loads. Special attention was given to the calibration for high flow periods during which large sediment loads were generated. Prior to manual calibration, auto calibration and sensitivity analysis were also performed though the SWAT-CUP model by giving the parameters ranges as described in the SWAT-CUP manual. The SWAT-CUP was then set to run for 500 simulations, with a uniform distribution of the parameters based on the upper and lower limits. Simulation was also conducted, if required, as per the new parameter sets recommended by the SUFI-2 algorithm [34]. This exercise narrowed down the ranges of new parameter sets and confirmed fast conversion.

3.5. Best Management Practices (BMPs)

A land management survey was conducted by LTVCA to collect field-by-field information from the farmers of the cropped area in the J1 and J2 watersheds. Information available in the land management survey was extrapolated to the fields with missing data. Land management datasets were completed by synthesizing data based on the adjacency principle, to create synthesized and complete datasets. The land management datasets’ preparation involved four steps: (i) planting details for three growing seasons, (ii) fertilizer application details, (iii) harvest/kill details, and (iv) tillage operations details. SWAT land use codes were allocated to represent field-wise land management data. The SWAT modelling gauged the main BMPs promoted and applied in the watershed. The main BMPs included conservation tillage, fertilizer incorporation, precision nutrient management, cover cropping, vegetative buffer strip establishment, and extensive tile drainage. Some possible BMPs (banding P, no-till, cover crop) were applied in the J1 and J2 watersheds to test their effectiveness. The emphasis was to test their impact on reducing P loads at the watershed outlets.

3.5.1. Implementation of BMPs

The baseline model was set up with BMPs implemented in the J1 and J2 watersheds for the period of 2016–2017 (i.e., banding and incorporation of P, variable rate P application, below canopy application of P, variable rate of P + incorporation of P, and conservation tillage and cover crops). Historical land management practices and crop rotation data for seven years (2011–2017) provided to the model for such fields as J165, 61, 62, 63, and 75, and J204, J2E4, J205, J207, J212, J213, J224, J232, J233, J238, J239, J241, and J244, while, for the remaining fields, the two-year land use datasets (2016–2017) were used. Details of the BMPs that were practiced in the Jeannette Creek watershed in 2016 and 2017 are given in Table S2. This table also provides a field-wise description, along with their modified parameters used during the model built-up for existing BMPs.

3.5.2. Potential BMPs for J1 and J2 Watersheds

An additional scenario was evaluated to show the possible impacts on P load reduction of expanding the existing BMP effort. Hence, some possible BMPs (i.e., banding P, no-till, cover crops) were also applied on the J1 and J2 watersheds to test their effectiveness. Firstly, fertilizer banding was applied to all the fields in both watersheds, to test its potential impact on P load reduction. For this, the parameter FERT_SURFACE (a fraction of the fertilizer applied to the first 10 mm of the soil layer) was set to 0.01, indicating that 99% of the fertilizer was applied beneath the top 10 mm of soil. It was hypothesized that a lesser availability of fertilizer in the topsoil would lead to less P being transported from flat lowland fields/watersheds to the stream.

Secondly, tillage practice was such that soybean [Glycine max (L.) Merr.] and winter wheat (Triticum aestivum L.) were planted under no-till practices, whereas corn (Zea mays L.) would be planted under vertical tillage (a.k.a. zero-till). A new tillage ID (#108) was assigned to vertical tillage operations, with the depth of tillage and mixing efficiency both set to zero. Lastly, cover crops were applied to all the fields after the harvest of winter wheat to assess how much amount of reduction in P exports was achieved just before reaching (i.e., just upstream of) the watershed’s pumping outlet.

4. Results and Discussion

4.1. Flow Calibration

Streamflow calibration for the Jeannette Creek J1 and J2 sub-watersheds focused on improving model performance at the flow monitoring stations at pump outlets. Table S3 lists the various parameters related to flow hydrology (i.e., snow and snowmelt-related parameters, groundwater parameters, and relative parameters like CN2, SOL_K, SOL_ALB, SOL_AWC, etc.) which may have spatial patterns that vary from HRU to HRUs at the field level. Hence, these parameters were selected for the stream flow model calibration, and the final specified parameter best values, ranges, and sensitivity are listed in Table S3.

In this case, the SWAT model default value was 0.5, and the range was −2 to 2 for both J1 and J2. The best values were −1.38 and −0.6 for J1 and J2, respectively. For the case of the soil and evaporation compensation factor, the most sensitive SWAT parameter, its default value and range were ESCO, 0.95, and 0.9–1.0, respectively, while the best values were 0.9795 for J1 and 0.955 for J2. These results fell into the high-sensitivity category. In the calibration exercise, the baseline for the relative changes ranged from −0.1% to +0.1% or −10% to +10% of the parameter value. A similar training was conducted for multiple case parameters until the calibration of the SWAT model was completed.

Figure 7 shows a time series plot for the J1 watershed, showing stream flow and precipitation with a 95% parameter predictive uncertainty (PPU) band, generated by SWAT-CUP and based on the range of the parameter value for the J1 watershed. The low and high peaks of streamflow in the watershed are commensurate with precipitation amounts. Moriasi et al. [35] suggests that the overall qualitative rating can be considered unsatisfactory if PBIAS > 15%, NSE < 0.5, and R² < 0.6. However, while the model predictions for J1 tended to follow the trend of observations, the goodness-of-fit statistics reflected the central tendency—the PBIAS was quite high (48.20%), while values of NSE = 0.23 and R² = 0.33 were found to be reasonable.

The unreasonable values (i.e., substantial (48%) over estimations) generated by the model could be attributed to the non-continuity in observed streamflow at the watershed’s pump station outlet. On the other hand, the 95% PPU and a p-stat equal to 0.33 indicates that 33% of the observations were compressed within a comparatively thin band (r-stat 0.22). This further indicates that although model performance statistics were in the unacceptable range, the stochastic-based calibration indicated a reasonable prediction of streamflow.

Figure 8 displays a time series plot of stream flow and precipitation, with the 95% PPU band, generated from the SWAT-CUP based on the range of parameter value for the J2 watershed. In this plot, the final range of streamflow parameter values were used in an uncertainty analysis. The low and high peaks of streamflow from the watershed corresponded with precipitation amounts. Showing model performance parameter statistics, Table S4 reflects an overall qualitative rating that can be considered as ‘good,’ based on the values of PBIAS (9%), NSE (0.6), and R² (0.71) [35]. In addition, the 95% PPU indicates that 30% of the observations (p-stat = 0.3) were compressed within a comparatively thin band (r-stat 0.12), which further confirms that the model performance statistics are within the acceptable range.

4.2. Sediment Calibration

The results displayed in Figure 9 and Figure 10 represent time series plots featuring simulated sediment concentration and precipitation, with a 95% PPU band for both watersheds. The low and high peaks of streamflow in the watershed corresponded with a similarly high precipitation amount. While a thorough model validation was hindered by the unavailability of continuous observations, the plot indeed shows that the model simulations captured the trend of observations during the calibration periods.

The simulation results show a very low (0.02 Mg ha⁻¹) annual sediment load, which means the simulated sediment concentration was near zero. A low annual sediment load can be linked with low moisture conditions during the growing season, such that summer rainfall events generate localized and insignificant runoff (Figure 9). The high observed sediment concentrations were indeed unexpected, but might be attributed to a very localized event covering a small area around the rain gauge or potential issues in the sampling techniques employed. Further detailed investigation is necessary to address and understand this issue.

The model performance statistics displayed a ‘very good’ PBIAS value of −5.5% (<10%), coupled with ‘unsatisfactory’ NSE (0.14) and R² (0.1) values (Table S5). Overall, based on the PBIAS sediment simulation, the results for the pump outlet were deemed ‘reasonable’. Lower values of NSE and R² were indeed expected, as the observations were instantaneous, while the SWAT simulations were generated on a daily time frame. As for the uncertainty band, 44% of the observations (p-stat = 0.44) were encapsulated in a relatively thick (r-stat = 1.05) uncertainty band. The calibrated J1 model indicated the average annual simulated sediment load to be 0.76 Mg ha⁻¹ over the simulation period (2016–2017).

The SWAT model was also set up and run for the J2 watershed for the period of 1 January 2016 to 31 December 2017. The simulated average annual sediment yield was 0.39 Mg ha⁻¹ at the pump outlet. Over the simulation period, the high and low sediment loads corresponded with the simulated streamflow pattern. The simulated sediment loads for J2 were found to be very low during the summer season (Figure 10). Model performance statistics showed a ‘very good’ PBIAS value of 8.3% (<10%, Moriasi et al. [35]) and ‘unsatisfactory’ NSE (0.13) and R² (0.14) values at the pump outlet. Based on the PBIAS and instantaneous observed sediment concentration, the performance of SWAT to simulate sediment loads appeared to be reasonable.

4.3. Phosphorus Calibration

Supplementary Table S2 presents the SWAT model parameters for P loads for the J1 and J2 watersheds. Beyond listing the parameters, the table provides comprehensive information such as descriptions, model ranges, optimal values, and sensitivities.

Following the calibration for streamflow and sediment, the SWAT model was further calibrated for P in the J1 and J2 watersheds using SWAT-CUP. The calibrated J1 model results indicated that the 2016 annual organic and mineral P (${\mathrm{P}}_{\mathrm{o}\mathrm{r}\mathrm{g}}$ and ${\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$) were about 0.14 kg ha⁻¹ and 0.09 kg ha⁻¹, respectively, whereas, in 2017, ${\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ was 0.18 kg ha⁻¹ and ${\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ was 0.11 kg ha⁻¹. Low and high ${\mathrm{P}}_{\mathrm{o}\mathrm{r}\mathrm{g}}$ and ${\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ corresponded with increases and decreases in streamflow during the simulation period (Figure 11 and Figure 12). The simulated average annual ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ load from this watershed was 0.26 kg ha⁻¹ (Table S6).

Model performance statistics for ${\mathrm{P}}_{\mathrm{o}\mathrm{r}\mathrm{g}}$ and ${\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ calibration at the pump outlet of the J1 watershed monitoring station showed a PBIAS of 15.7% for ${\mathrm{P}}_{\mathrm{o}\mathrm{r}\mathrm{g}}$ and of −17.4% for ${\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$. The PBIAS value for ${\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ fell within the range of ‘good’ qualitative rating criteria, whereas NSE and R² were deemed unsatisfactory [35,36,37]. For ${\mathrm{P}}_{\mathrm{o}\mathrm{r}\mathrm{g}}$, the 95% PPU band summarized 42% of the observations in a relatively thick band (r-stat = 1.04). For ${\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$, the uncertainty band encapsulated only 48% of the observations (p-stat = 0.48) in a relatively thick band (r-stat = 1.18), which again highlights the difficulty in comparing observed instantaneous P concentrations with daily simulated SWAT results (Figure 13, Figure 14 and Table S6).

For the J2 watershed, the time series plot was generated by comparing simulated and observed ${\mathrm{P}}_{\mathrm{o}\mathrm{r}\mathrm{g}}$ and ${\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ and plotting a 95% PPU band. The simulated values of ${\mathrm{P}}_{\mathrm{o}\mathrm{r}\mathrm{g}}$ were 0.10 kg ha⁻¹ for 2016 and 0.11 kg ha⁻¹ for 2017. The simulated values of ${\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ were 0.08 kg ha⁻¹ for 2016 and 0.13 kg ha⁻¹ for 2017.

The overall average annual ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ yield over the period of two years was 0.20 kg ha⁻¹. The low and high ${\mathrm{P}}_{\mathrm{o}\mathrm{r}\mathrm{g}}$ and ${\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ corresponded well with high and low streamflow patterns at the pump outlet. The model performance statistics showed a very good PBIAS (6%) for ${\mathrm{P}}_{\mathrm{o}\mathrm{r}\mathrm{g}}$ and a satisfactory one (22%) for ${\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$. The NSE for ${\mathrm{P}}_{\mathrm{o}\mathrm{r}\mathrm{g}}$ (0.15) and ${\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ (−0.03), along with the R² for ${\mathrm{P}}_{\mathrm{o}\mathrm{r}\mathrm{g}}$ (0.16) and ${\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ (0.03), were found to be ‘unsatisfactory.’ Thus, for ${\mathrm{P}}_{\mathrm{o}\mathrm{r}\mathrm{g}}$, the uncertainty band summarized only 62% of the observations (p-stat = 0.62) and was compressed in a relatively thick band (r-stat = 1.82), whereas, for ${\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}},$ the 95% PPU band summarized 77% of the observations in a relatively thick band (r-stat = 1.94). Based on the PBIAS, the model’s performance in predicting temporal trends was satisfactory, but it was unsatisfactory for concentrations. Again, this is due to the comparison of instantaneous observed concentration with the daily simulated concentrations. Overall, based on the trends, the performance of the model to simulated P was poor but approaching a reasonable range [35,36,37].

Considering existing land management conditions, average yearly sediment and nutrient yield were simulated. Sediment loads at the J1 and J2 outlets were 82.77 Mg yr⁻¹ and 84.9 tons, respectively, representing sediment yields of 0.54 and 0.39 Mg yr⁻¹ ha⁻¹. The ${\mathrm{P}}_{\mathrm{o}\mathrm{r}\mathrm{g}}$ for J1 and J2 were 261.5 kg yr⁻¹ and 93.17 kg yr⁻¹, respectively, corresponding to yields of ${\mathrm{P}}_{\mathrm{o}\mathrm{r}\mathrm{g}}$ of 0.165 kg yr⁻¹ ha⁻¹ and 0.1 kg yr⁻¹ ha⁻¹. Additionally, the ${\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ for J1 and J2 were 99.7 kg yr⁻¹ and 269.5 kg yr⁻¹, respectively, with associated yields of 95 g yr⁻¹ ha⁻¹ and 90 g yr⁻¹ ha⁻¹. The ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ loads for J1 and J2 were 361.2 kg yr⁻¹ and 362.6 kg yr⁻¹, respectively, with corresponding yields of 0.26 kg yr⁻¹ ha⁻¹ and 0.19 kg yr⁻¹ ha⁻¹.

4.4. BMP Scenarios

4.4.1. Effectiveness of Current BMPs

Assessing the effectiveness of current BMPs using the SWAT model was undertaken by considering the baseline (existing) BMP scenarios employed in building up the model for the J1 and J2 (Jeannette creek) watersheds. The effectiveness of each individual BMP scenario was assessed by comparing it against the baseline scenario with respect to flow, sediment, and P yields (Figure 15).

Table 2. Modelling results showing effectiveness of existing BMPs just upstream of the Dauphin (J1) and Deary pump (J2) outlets (2016–2017).

Scenario	Year	Watershed	Flow	Sediment	${\mathbf{P}}_{\mathbf{o}\mathbf{r}\mathbf{g}}$	${\mathbf{P}}_{\mathbf{m}\mathbf{i}\mathbf{n}}$	${\mathbf{P}}_{\mathbf{t}\mathbf{o}\mathbf{t}}$
			(m3 s−1)	(Mg yr−1)	(kg yr−1)	(kg yr−1)	(kg yr−1)
All BMPs Applied	2016	J1	0.08151	93.02	243.5	89.9	333.4
		J2	0.052	63.21	76.25	199.6	275.85
	2017	J1	0.06243	72.54	279.5	109.5	389
		J2	0.080	106.6	110.1	339.4	449.5
All BMPs Removed	2016	J1	0.082	93.18	249.8	122.1	371.9
		J2	0.053	64.66	105	204.6	309.6
	2017	J1	0.064	73.97	287.2	135.8	423
		J2	0.082	110.3	169.2	340.4	509.6
Reduction
Effective Reduction	2016	J1	0.0002	0.16	6.3	32.2	38.5
		J2	0.001	1.45	28.75	5	33.75
	2017	J1	0.0015	1.43	7.7	26.3	34
		J2	0.003	3.7	59.1	1	60.1
Relative Reduction (%)
Effective Relative Reduction	2016	J1	0.22	0.17	2.52	26.37	10.35
		J2	1.50	2.24	27.38	2.44	10.90
	2017	J1	2.38	1.93	2.68	19.37	8.04
		J2	3.15	3.35	34.93	0.29	11.79

shows the simulation results for the change in P loads under the existing BMPs applied in the J1 and J2 watersheds and when all present BMPs are removed. The 2016 and 2017 mean simulated ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ discharge (i.e., upstream, just before pump outlet) from the J1 and J2 watersheds was 361.2 kg yr⁻¹ and 362.7 kg yr⁻¹, respectively, representing ${\mathrm{P}}_{\mathrm{o}\mathrm{r}\mathrm{g}}$ loads of 261.5 kg yr⁻¹ and 93.2 kg yr⁻¹ (J2), and ${\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ loads of 99.7 kg yr⁻¹ and 269.5 kg yr⁻¹. For the J1 watershed, the average non-growing season (November to April) ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ load was 196.8 kg yr⁻¹, significantly greater than that over the growing season (May to October; 164.3 kg yr⁻¹). Simulation results suggest that the existing BMPs were more effective during the non-growing season months (17% reduction in the transportation load of ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$) compared to growing season months (2%).

The J2 watershed saw average ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ loads of 214.2 kg yr⁻¹ transported during the non-growing season (November to April) compared to 148.5 kg yr⁻¹ during the growing season (May to October). The simulation results suggest that the existing BMPs were more effective during the non-growing season months (21% reduction in the transportation load of TP) compared to growing season months (7% reduction). The simulation results also suggest that the effectiveness of the existing BMPs was greater during spring and summer compared to winter and fall.

4.4.2. Scenarios for Existing BMPs

Retiring the P Application BMP from the Existing BMPs for J1 and J2 Watersheds

In the J1 watershed, P application as a BMP was implemented on 14 of 49 fields. Figure 16 (J1) shows the result of retiring this BMP on P loads in the J1 watershed. On an annual scale, the effectiveness of this BMP was found to be 0.4%. Simulation results also suggest that under a seasonal retirement of P application, P_tot would increase by 1% during the growing season, whereas no change would occur during the non-growing season.

On the J2 watershed, the banding of P was applied to all 54 fields, but P incorporation was implemented in only 6 fields, while variable application with incorporation was applied to only 3 fields. One remaining field employed below-canopy application as a BMP. Figure 16 (J2) shows that retiring existing P application BMPs from the baseline model resulted in a moderate effect on the transportation load of ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$. Retiring these BMPs led to no change in P loads during the non-growing season and a 10% increase during the growing season. However, the greatest increase in P loads occurred during the winter season.

Retiring Conservation Tillage from the Existing BMPs

In the J1 watershed, 7 out of 49 fields were subject to conservation tillage. Retiring existing conservational tillage from the baseline model increased the annual ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ losses to 375.1 kg from 361.2 kg. This limited increase (3.84%) can be linked with the fact that only 14% (by area) of the fields in the watershed were under conservation tillage. Thus, retiring this BMP increased the ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ loads by 2% during the growing season and 5% during the non-growing season. The retirement of conservational tillage in the J2 watershed increased the annual P loads to 366.2 kg from 362.6 kg. Roughly 1% increases in the P loads were observed during the growing and non-growing seasons. This small increase can be linked to the fact that only 4% of the fields in this watershed were under conservation tillage.

Retiring the Cover Crop from the Existing BMPs

In the J1 watershed, 2 of 49 fields during 2016 and 5 of 49 fields during 2017 were under cover crops. Retiring the existing cover crops in the J1 watershed did not show any significant impact (change was within the range of ± 0.11%) on the total annual and seasonal P loads. The change in annual and seasonal ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ loads achieved by retiring the existing cover crop BMPs in the J2 watershed was very similar to that observed for the J1 watershed.

Retiring Vegetative Buffer Strips from the Existing BMPs

In 2016, the J1 watershed had 2 of 49 fields under Vegetative Buffer Strips (VBS), while watershed J2 had none. Retiring the existing VBS resulted in a less than 2% increase in total annual P loads and no more than a 3% increase in seasonal ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ loads. Retiring VFS in the J1 watershed led to an increase in the export of annual ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ loads from 361.2 kg to 367.7 kg. This slight increase in P load can be linked with the small fraction of fields (around 4%) under the VFS BMP.

4.4.3. Scenarios for Possible BMPs

The calibrated model was used to evaluate the potential effectiveness of possible BMPs (e.g., banding P, no-till, and cover crops) on reducing P loads just upstream from the pump outlets. These BMPs were applied in all the fields in the J1 and J2 watersheds, and the export of ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ was documented for both watersheds. Each potential scenario was compared with the baseline scenario, and its effectiveness was computed at various temporal scales (e.g., annual, seasonal). The application of banded fertilizer, undertaken for all fields in both watersheds, was compared with broadcast fertilizer application. Under the banded application, the fertilizer was placed 10 mm below the soil layer at a rate inferior to that of a surface application by spreading [21]. In the J1 watershed, export of ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ was reduced from 361 kgyr⁻¹ (base case value) to 343.1 kg yr⁻¹ (banded application of fertilizer case), reflecting an estimated reduction of up to 5% on an annual scale, with reductions of 8% and 3% in the growing and non-growing seasons, respectively. The reductions in ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ loads were of 4%, 8%, 3%, and 1% for the winter, spring, summer, and fall seasons, respectively. For the J2 watershed, the banded fertilizer resulted in almost the same reduction (8%) in the export of ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ on an annual scale, growing season, and non-growing season. The reduction in ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ loads during the four seasons were 6% in winter, 7% in spring, and 6% in summer, with the maximum reduction (11%) occurring during the fall season.

The simulation results presented in Figure 17 (J1 and J2) for no-tillage show a more than 60% reduction in annual ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ loads in both watersheds, making it the most efficient BMP among all the BMPs under consideration. The annual ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ loads were reduced from 361 kg to 103 kg in watershed J1 and from 362 kg to 114 kg in watershed J2. No-till reduced soil compactibility (Proctor bulk density) and increased the stability of wet aggregates and water infiltration [38]. No-till also enhances biological activities [39,40]. These physical properties of undisturbed soil might have contributed to reducing surface runoff, decreased the rate of soil erosion, and thereby reduced the transport of ${\mathrm{P}}_{\mathrm{o}\mathrm{r}\mathrm{g}}$ to the streams.

The ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ load reduction pattern in the growing and non-growing seasons was similar for both watersheds. For the non-growing season, the reduction in P transport was 76% in J1 and 66% in J2, whereas, in the growing season, the rate of reduction was 67% for J1, and 52% for J2. In the J1 watershed, the maximum reduction (82%) occurred during winter, followed by the spring (76%), fall (75%), and summer (37%) seasons. For the J2 watershed, the maximum reduction (69%) occurred during the fall, followed by the spring (64%), winter (55%), and summer (24%) seasons (Figure 17).

The cover crop BMP—cereal rye/red clover planted after winter wheat—was applied to all fields in both the J1 and J2 watersheds. For J1, the cover crop offered a small reduction (2%) in annual ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ loads (Figure 18). The seasonal pattern in ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ load reduction was different than for other BMPs: a 15% decrease in ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ loads during the non-growing season, but a 15% increase in ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ loads during the growing season. Most of the increase (93%) occurred during the summer season. For other seasons, the maximum absolute change in ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ loads under cover cropping occurred during the winter (−39%), followed by the spring (−12%) and fall (−2%) seasons. Moreover, the annual ${\mathrm{P}}_{\mathrm{o}\mathrm{r}\mathrm{g}}$ and ${\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ loads just upstream of the pump outlet decreased on average by 6% and 7%, respectively.

5. Conclusions

The SWAT model can evaluate the performance of BMPs in reducing P loads in flatland watersheds with pumped outlets. The model was applied to two (J1 and J2) small flatland watersheds in southern Ontario, in the Canadian Great Lakes basin. The BMP evaluation included banded fertilizer and variable rate application, conservation tillage, no-till, and a post-winter-wheat cover crop. The calibration performance for stream flow was in the acceptable range for both the watersheds, with watershed J2 performing slightly better than J1. However, the performance for the simulation of sediment concentration remained below the acceptable range due to the limited availability of measured sediment concentration data. The SWAT model successfully predicted temporal P trends but had difficulty estimating ${\mathrm{P}}_{\mathrm{o}\mathrm{r}\mathrm{g}}$ and ${\mathrm{P}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ concentrations. Again, this was due to the limited availability of water quality data. The existing BMPs resulted in an average annual total P load reduction of 9% (36.25 kg) for J1 and 11% (46.9 kg) for J2. The P application rate and conservation tillage proved to be more effective in reducing ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ loads than cover crops after winter wheat. Altering the P application rate resulted in 0.5% (J1) to 4% (J2) reductions in ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ loads. Conservation tillage decreased ${\mathrm{P}}_{\mathrm{t}\mathrm{o}\mathrm{t}}$ loads by 3.8% in J1 and 1% in J2. Among the BMPs that could potentially be applied in the watersheds under study, banded P application resulted in a 5% reduction in P loads for J1 and 8% for J2. However, the implementation of no-till was found to be the most effective in reducing P loads (>60%) in both watersheds. This evaluation is critical for successful implementation and assessment of effectiveness of BMPs in flatland agricultural landscapes where water management through a series of pumped outlets plays a pivotal role. The findings of this research have far-reaching implications for future agricultural practices and water resource management, such as (a) enhancing the water resource management in similar watersheds by identifying and recommending suitable BMPs for improving water quality; and (b) developing appropriate land management strategies to target resources in areas contributing relatively higher sediment and nutrient losses.