Modeling and Prioritizing Interventions Using Pollution Hotspots for Reducing Nutrients, Atrazine and E. coli Concentrations in a Watershed

: Excess nutrients and herbicides remain two major causes of waterbody impairment globally. In an attempt to better understand pollutant sources in the Big Sandy Creek Watershed (BSCW) and the prospects for successful remediation, a program was initiated to assist agricultural producers with the implementation of best management practices (BMPs). The objectives were to (1) simulate BMPs within hotspots to determine reductions in pollutant loads and (2) to determine if water-quality standards are met at the watershed outlet. Regression-based load estimator (LOADEST) was used for determining sediment, nutrient and atrazine loads, while artiﬁcial neural networks (ANN) were used for determining E. coli concentrations. With respect to reducing sediment, total nitrogen and total phosphorus loads at hotspots with individual BMPs, implementing grassed waterways resulted in average reductions of 97%, 53% and 65% respectively if implemented all over the hotspots. Although reducing atrazine application rate by 50% in all hotspots was the most effective BMP for reducing atrazine concentrations (21%) at the gauging station 06883940, this reduction was still six times higher than the target concentration. Similarly, with grassed waterways established in all hotspots, the 64% reduction in E. coli concentration was not enough to meet the target at the gauging station. With scaled-down acreage based on the proposed implementation plan, ﬁlter strip led to more pollutant reductions at the targeted hotspots. Overall, a combination of ﬁlter strip, grassed waterway and atrazine rate reduction will most likely yield measureable improvement both in the hotspots (>20% reduction in sediment, total nitrogen and total phosphorus pollution) and at the gauging station. Despite the model’s uncertainties, the results showed a possibility of using Soil and Water Assessment Tool (SWAT) to assess the effectiveness of various BMPs in agricultural watersheds.


Introduction
There is a growing concern about sediments and nutrients entering both surface water and groundwater bodies from nonpoint sources, and causing deterioration in water quality. Agricultural nonpoint source pollution (NPS) is the primary source of water quality degradation in different countries and regions [1][2][3][4]. It is also a major source of impairments to wetlands, and a key contributor to contamination of estuaries and groundwater according to the United States Environmental Protection Agency [5]. Activities such as overgrazing, improper fertilizer and pesticide application, poorly located or managed animal feeding operations and manure, as well as plowing too frequently and at the wrong time, may adversely impact water bodies. The activities which cause NPS in agriculture mostly occur due to the nonexistence or lack of adoption of conservation plans [6].
The movement of sediments, nutrients and pesticides from agricultural NPS into receiving waters can be successfully controlled by implementing best management practices (BMPs) [7,8]. Usually, BMPs are implemented at a field scale even if the desired water-quality goals are required at a sub-watershed or watershed scale. The goal of any BMP implementation plan is to reduce the pollutant loads to allowable levels in a watershed while minimizing the infrastructural and maintenance cost [9]. However, due to the spatio-temporal variability of sources and transport pathways, it is difficult to effectively mitigate agricultural pollution [10,11]. In addition, important factors such as the hydrology, meteorology, and geomorphology at field, sub-watershed and watershed scales also complicate the response of BMPs to NPS pollution [12,13]. Thus, the BMPs to be implemented in a watershed need to be of the right type and magnitude (acreage), and also be in the right place, in the right combinations, and at the right adoption level within different spatial scales [14][15][16].
Streams in the Big Sandy Creek Watershed (BSCW), which is located within the Little Blue River watershed, Nebraska, are impaired by high concentrations of the bacterium Escherichia coli (E. coli), the herbicide Atrazine and high concentrations of nitrogen and phosphorus. In order to better understand the pollutant sources and the potential to successfully implement BMPs, the Natural Resources Conservation Service (NRCS), Nebraska Department of Environment and Energy (NDEE) and Little Blue Natural Resources District (LBNRD) conducted a detailed assessment of the BSCW which resulted in the identification of target sub-watersheds (or hotspots). According to the Big Sandy Creek National Water Quality Initiative Watershed Implementation Plan (BSCW-NWQI) by the United States Department of Agriculture-Natural Resources Conservation Service (USDA-NRCS) [17], if the pollutant loads from these target sub-watersheds are reduced, this could have the most significant impact on water quality in Big Sandy Creek in the near future. Although Nebraska has no criteria for sediment and nutrients in flowing waters, the NWQI project set a goal of 20% reduction of these pollutants while total maximum daily loads (TMDLs) of 126 colony-forming unit (cfu)/100 mL and 12 µg/L were developed for E. coli and Atrazine respectively [17].
Though the SPAtially Referenced Regression on Watershed attributes (SPARROW) model was used to identify the location and origin of these pollutants for the Little Blue River watershed, there is a need for an alternative modeling tool that considers both hydrological and agricultural activity factors in a watershed to quantify the reduction of pollutant loads due to BMPs. The Soil and Water Assessment Tool (SWAT) is a watershedscale hydrological and water quality modeling tool used for predicting streamflow and pollutant losses from watersheds made up of mixed land covers, soils and slopes [18]. The model was developed to help water resource managers in their assessment of water quantity and/or quality in watersheds and as a tool for evaluating the impact on BMPs' implementation.
The SWAT model has been extensively used worldwide for over three decades [19,20]. The model is process-based and can simulate the hydrological cycle, crop yield, soil erosion, and nutrient transport. It has been used globally to successfully simulate sediment [21], nitrogen [22,23], phosphorus [24,25], Atrazine [26] and E. coli [27] as well as BMPs [15,28,29]. Mittelstet et al. [28] used the SWAT model to identify a combination of potential land management practices in Oklahoma to meet the water-quality standard for total phosphorus (0.037 mg/L) in three rivers in Oklahoma. The impact of various BMPs in reducing the sediment and phosphorus rates in Grand River's discharge to Lake Erie in southern Ontario was investigated by Hanief and Laursen [25] using SWAT. Ikenberry et al. [23] used SWAT to simulate daily flow, nitrate concentrations in tile-drains, and soil-nitrogen dynamics.
While most of the previous SWAT studies are hypothetical and involve blanket BMPs that are implemented across entire watersheds, there are very few studies that considered implementing BMPs at hotspot sub-watersheds [15,30,31]. Lamba et al. [15] studied the effect of BMP implementation at sub-watershed and watershed scales on sediment and phosphorus load reductions at hotspots in south-central Wisconsin. They observed 66-99% reductions in sediment and total phosphorus loads at hotspot sub-watersheds after converting croplands to grasslands, and 14-25% reductions after implementing no-till practices at these hotspots. However, at the watershed outlet, sediment and phosphorus reduction was ≤15% after implementing both BMPs because only about 8% of the watershed area was targeted for BMPs. In order to reduce nitrate and total suspended sediment yields from hotspot areas in Raccoon River watershed in west-central Iowa, USA, Teshager et al. [30] analyzed fourteen BMPs based on systematic combinations of five strategies: fertilizer/manure management, changing row-crop land to perennial grass, vegetative filter strips, cover crops and shallower tile drainage systems. Their results indicated that sufficient reduction of nitrate load may require either implementing multiple BMPs or converting extensive areas into perennial grass in order to meet and maintain the drinking water standard. Himanshu, et al. [31] identified erosion hotspots in Marol watershed, India and evaluated the implementation of BMPs in those critical erosion-prone areas. They observed sediment and nutrient load reductions with alternate cropping of groundnut and soybean when compared to paddy and maize cultivation. The very few studies that consider hotspots usually evaluate different BMPs and give best-scenario recommendations based on the assumption that funding for implementation would not be a limiting factor.
Although few recent studies have shown that the cost-effectiveness of BMPs in mitigating pollution is greatly boosted by using BMP selection and placement tool at watershedscale [32,33], studies which incorporate hotspots and fund limitations, not only for sediment and nutrient load reduction, but also for Atrazine and E. coli, are rare. Thus, based on the financial constraint of $500,000 annually for three years of the NWQI, and the detailed preliminary assessment conducted by the NDEE, the objectives of this study are to: (i) simulate BMPs proposed by the NDEE within target sub-watersheds to determine reductions in pollutant loads and (ii) determine if water-quality standards are met at the watershed outlet.

Study Area
The Big Sandy Creek Watershed (BSCW) is located along the northern edge of the Little Blue River watershed, Nebraska ( Figure 1). The watershed is located between longitude 97 • 17 29 and 98 • 18 24 W and latitude 40 • 12 36 and 40 • 34 40 N. The area of BSCW is 1610 km 2 and the mean annual precipitation is 640 mm. Most of the precipitation occurs in the months of May, June and July. The climate is dry and continental, characterized by hot summers and cold winters. High winds, periodic droughts and flooding are climatic factors in this region. The land surface slopes from west-northwest to east-southeast with a maximum elevation of 1909 m to 1344 m above mean sea level. The topography is a relatively flat upland and slightly rolling hills, with narrow valley regions of low relief along rivers and streams. About 44% of the watershed has less than 3% slope, 29% has between 3-6% slope and the remaining 27% has more than 6% slope. The watershed lies within the Central Loess Plain Major Land Resource Area. Water-deposited sand and gravel occur as a heterogeneous mixture. Soils of the Jansen and Meadin series formed in loess overlaying the sand and gravel. Below is semi-consolidated bedrock or consolidated bedrock. The unconsolidated materials are primarily windblown clayey silts and loess overlying sands and gravels of alluvial origin. There are also alluvium soils that consist of valley sediments recently washed from the uplands and deposited on flood plains and stream terraces. There are 19.2 km 2 of pasture in the target sub-watersheds. About 25,000 cattle (out of 36,000 in the whole BSCW) were in target sub-watershed 1a. Most pastures are used for summer grazing, but some are used for wintering cows prior to calving and staging pastures where cow-calf pairs are held after calving and prior to the summer pasture. Pastures often are overgrazed and void of vegetation along the streambanks. Moreover, cattle usually have direct access to the streams and they deposit fecal material while crossing, drinking from or loafing in the stream.  BSCW is predominantly an agricultural watershed with approximately 76% row crops (mainly corn and soybeans) and 23% pasture/hay. The remaining 1% includes urban, forest, water and other non-agricultural uses. Crop rotation primarily is corn and soybeans which accounts for 49% of the watershed while 27% is in rotation with at least two years of continuous corn. On average, approximately 15 km 2 of perennials crops such as alfalfa are grown annually.
The total area of the targeted sub-watersheds ( Figure 2) is 240 km 2 and the mean annual precipitation is about 737 mm per year with approximately 80% coming as rain in high-intensity, convective thunderstorms. May and June are considered the critical erosion period following spring tillage and planting. Based on estimates from land use and crop rotation data for the target sub-watersheds, 49% of cropland is in a corn-soybean rotation and 36% is in a rotation with at least two years of continuous corn.
There are 19.2 km 2 of pasture in the target sub-watersheds. About 25,000 cattle (out of 36,000 in the whole BSCW) were in target sub-watershed 1a. Most pastures are used for summer grazing, but some are used for wintering cows prior to calving and staging pastures where cow-calf pairs are held after calving and prior to the summer pasture. Pastures often are overgrazed and void of vegetation along the streambanks. Moreover, cattle usually have direct access to the streams and they deposit fecal material while crossing, drinking from or loafing in the stream.

SWAT Model Set-Up
ArcSWAT (v2012) was used to divide the BSCW into 53 sub-watersheds, which were further divided into 1860 hydrologic response units (HRUs) based on the following thresholds: 5% land use, 0% soil and 0% slope. Each HRU was defined as a sub-watershed with similar land use, soil type and slope. Data used in SWAT model set-up/calibration/validation are presented in Table 1. The outlet of the BSCW was specified by a point downstream of the NDEE water quality monitoring station and the Nebraska Department of Natural Resources (DNR) streamflow gauging station 06883940 (latitude 40.235 and longitude −97.388). A lake at the US Meat Animal Research Center (Lake MARC) close to the headwater sub-watershed was included as a reservoir in the model though there is usually no outflow from this reservoir year round ( Figure 1).   The NDEE provided monthly and/or bimonthly water-quality data on sediment (total suspended solids, n = 146), nitrogen (n = 146), phosphorus (n = 146), Atrazine (n = 105), and E. coli. (n = 39). Since daily data was not available for the water quality parameters, Load Estimator (LOADEST) was used to estimate daily loads for sediment, Sustainability 2021, 13, 103 6 of 22 nutrients and Atrazine, while artificial neural network (ANN) was used to estimate daily E. coli concentrations. LOADEST estimates loads for a user-defined time period from discrete water quality samples and measured daily flow using a formulated regression model. Runkel et al. [34] provided detailed information on the LOADEST program and its techniques and methods. An ANN is a mathematical model comprising of a network of computation nodes with established connections between them [35]. Hsu et al. [36] and Sarle [37] provided detailed descriptions of ANN and its structure. An advantage of ANN is that it requires no a priori assumptions about the relationships between the independent and dependent variables as well as the mathematical functions to be used [38].
The prediction accuracy of an ANN model is strongly dependent of its generalization ability. A well-trained ANN model should have a very low root mean squared error at the end of the training and testing phases in order to avoid over-fitting and to generalize the model. An overfitted ANN model will give good predictions only for the training set (most of the time) but will perform poorly when predicting the test set (out-of-sample E. coli values). To develop the ANN model for predicting E. coli concentrations, the available dataset was randomly divided into a training dataset (70% of the total data), a validation dataset (15% of the total data), and a test dataset (15% of the total data). However, seeing that there were only thirty nine E. coli concentration observations for our study, the splitting of the dataset was done several times in order to get reasonable distributions for the training, validation, and test datasets which at least covered more than 80% of the range of E. coli values. This was an important step because any random split of our small E. coli dataset could give training, validation and test datasets with totally different data distributions based on the range of E. coli values covered. A model trained on a vastly different data distribution than the test set will perform inferiorly.
The fitting ANN was created using five hidden layers and Levenberg-Marquardt backpropagation algorithm in MATLAB (MathWorks Inc. Product R2019b). Training was done according to training parameters of the algorithm with their default values. Before ANN model development, the E. coli concentrations were log-transformed due to the presence of high concentration values. The initial predictor variables considered for predicting E. coli concentrations were cumulative precipitation on the day of sampling, cumulative precipitation a day before sampling, 7-day and 10-day moving averages of daily precipitation prior to sampling day, average streamflow on the day of sampling, and average streamflow a day before sampling. Later, the cumulative precipitation a day before sampling was removed due to low correlation with E. coli concentration. The coefficient of determination (R 2 ) and root mean squared log error (RMSLE) statistics were used for checking model performance (Equations (1) and (2), respectively). RMSLE was used instead of root mean squared error (RMSE) because outliers can greatly increase the prediction error. While RMSLE considers only the relative error between predicted and actual values, and the scale of the prediction error is nullified by the log-transformation, in contrast, RMSE value increases in magnitude if the scale of error increases.
where x i is the observed E. coli concentration (colony-forming unit (cfu)/100mL),x i is the predicted E. coli concentration (cfu/100 mL), x is the mean of the observed E. coli concentration (cfu/100 mL), and n is the total number of samples considered. To improve modeling results and reduce the uncertainty in our SWAT model findings, management practices within the watershed were characterized. The management practices included the following: tillage operations, yearly planting and harvesting dates, nutrient and pesticide application rates and timing, irrigation and cattle density. Data on location of terraces and irrigation were also used in developing a detailed watershed model. Atrazine was simulated as loads from cropland only while E. coli was simulated as loads from both corn fields and pasture. Survey showed that producers commonly apply 1.7 kg/ha Atrazine to corn fields at a rate of 1.1 kg/ha at or before planting and 0.6 kg/ha post emergence. Management input files were created to simulate corn/soybeans rotation on all lands classified as soybeans in the land use layer file (Table S1 in Supplementary Materials). A grazing rate of 3.23 ha/cow was used for all pastures. SWAT was run at monthly time step from 1989 to 2016 with a warm-up period of five years (from January 1989 to December 1993).

Model Calibration and Validation
In order to parameterize the SWAT model to the given set of local conditions and to reduce uncertainty in model predictions, the auto-calibration program SWAT-CUP SUFI-2 was used to calibrate the model for streamflow, sediment, total nitrogen (TN), total phosphorus (TP), Atrazine and E. coli at the NDNR gauging station 06883940 with several repetitions of 1000 simulations. The SUFI-2 method takes into account all sources of uncertainties such as those from the input data, conceptual model, model parameters and observed data [39]. Manual calibration was also applied to improve the results of auto-calibration [40,41]. Before calibration, the ranges of targeted parameters for calibrating flow and pollutants were first defined by physically meaningful values. The values of these parameters were then allowed to change by replacement, addition, and relative change to the initial values. The model was calibrated (01/2004 to 12/2016) and validated (01/1994 to 12/2003) on a monthly time step using observed data at the NDRN gauging station 06883940. Streamflow was first calibrated, then sediment, nutrients, Atrazine and E. coli in that order.
In order to ensure a good calibration procedure, we used multiple statistics, each covering a different aspect of the hydrograph, so that the whole hydrograph is covered [42]. This reduced the unnecessary emphasis that using a single statistic can place on matching one aspect of the hydrograph at the expense of other aspects while recognizing potential errors in the observed data [43]. The model performance was assessed using R 2 , Nash-Sutcliffe efficiency (NS), root mean square error (RMSE)-observations standard deviation ratio (RSR), and percent bias (PBIAS) as recommended by Moriasi et al. [42]. These statistics were calculated using Equations (1) and (3)-(5) respectively. R 2 values range from 0 to 1 and higher values indicate a better fit of the simulated data to the observed data. According to the general performance ratings for statistics for a monthly time step recommended by Moriasi et al. [42], values of R 2 > 0.5 are considered satisfactory. NSE ranges from −∞ to 1, where a value of 1 indicates an optimal model fit to the observed data. An NSE value between 0.75 and 1 shows a very good model fit, 0.65-0.75 is good, 0.50-0.65 is satisfactory, and <0.5 is considered unsatisfactory. RSR standardizes the RMSE with the standard deviation of the observed data. This ranges from 0 to large positive values with lower values indicating better model fit. RSR values between 0 and 0.5 indicate a very good model fit, 0.5-0.6 is good, 0.60-0.7 is satisfactory, and >0.7 is unsatisfactory. PBIAS is used to assess the average tendency of simulation to over-or under-predict the observation data. The optimal value of PBIAS is 0, with low values indicating better model fit. For streamflow, sediment, and nutrients, absolute PBIAS values <10%, <15%, and <25% respectively indicate a very good model fit, 10-15%, 15-30%, and 25-40% indicate good, 15-25%, 30-55%, and 40-70% indicate satisfactory, and >25%, >55%, and >70% indicate unsatisfactory. Positive PBIAS values show the model is underestimating bias while negative values indicate bias overestimation [44]. Visual comparisons between observed and simulated graphs were also carried out in addition to the statistics to see how peaks aligned in both observation and simulation.
where X obs i is the observed variable (streamflow or pollutant load), X sim i is the simulated variable, X mean is the mean of n variables, STDEV obs is the standard deviation of observations, and n is the number of observations.

Implementing BMPs in Hotspots
Hotspots for prioritizing installation of BMPs to reduce sediment, nutrients, Atrazine and E. coli runoff were identified from data collected during BSCW windshield surveys as well as data obtained from Agricultural Conservation Planning Framework (ACPF) and Agronomy technical notes. Although several potential BMPs were highlighted in the NWQI project, this study focused on six BMPs. These BMPs include no-till farming practice (NT), filter strip (FS), crop rotation (CR), grassed waterway (GW), terrace (TR), and reduced Atrazine application rate (RAR). The calibrated and validated SWAT model was used to evaluate the BMPs in terms of their effectiveness in reducing pollutant loads. Slope, soil type, and crop type were considered major factors in implementing the BMPs since these factors could influence pollutant runoff under normal farming and grazing operations.
The expected load reductions from BMPs as well as the total areas considered by the NWQI project for implementation of these BMPs in three years are shown in Tables 2 and 3 respectively. The load reductions for sediment, nitrogen and phosphorus are for the targeted subbasins ( Figure 2) while the target concentrations of E. coli and atrazine are for the Nebraska Department gauging station. Pollutant loads were simulated before and after the implementation of BMPs. The locations and type of BMPs were based on the data provided by the NWQI project. Although the configurations of the BMPs may be variable in the real world, we assumed that they were generally representative and were meant to reflect the general operations of the BMPs in the target sub-watersheds in BSCW (Table 4).   Table 4. BMPs in target sub-watersheds and representation in SWAT.

RAR
Reduced Atrazine application rate Application rates in corn reduced by 50% NT increases the amount of residue on soil surface after crop harvest and before planting the next crop. To represent no-till farming in the SWAT model, conventional tillage operations with field cultivators and tandem disk plows were replaced with no tillage by modifying the till.dat and management files (Table 4). Crop rotation switched the conventional continuous corn to the two-year corn-soybean rotation (Table S1 in Supplementary Materials). According to the SWAT conservation practice modeling guide [45], representing filter strip in SWAT included adjusting the fraction of total runoff from a field entering the most concentrated 10% of a filter strip (VFSCON), field area to filter strip area ratio (VFSRATIO), and the fraction of flow within the most concentrated 10% of the filter strip that is fully channelized (VFSCH). Grassed waterway was represented by adjusting channel Manning's roughness (CH_N2), channel cover factor (CH_COV2) and the channel erodibility factor (CH_COV1) [46][47][48]. Terrace was simulated by modifying the runoff curve numbers of cropland according to Strauch et al. [49] and modifying the Universal Soil Loss Equation practice factor (USLE_P) by using the values 0.12 (for slopes <3% and >6%) and 0.10 (for slopes between 3% and 6%) according to Haan et al. [50]. RAR was simulated by reducing Atrazine application rate by 50% in both rainfed and irrigated continuous corn and corn/soybean rotation.
It is important to note that whenever an empirical approach is used for simulating BMPs, it usually involves applying a linear approach to non-linear problems. In reality, due to the reason that one cannot expect any ratio or factor to be time-invariant irrespective of the hydro-meteorological, soil or land cover conditions, the results of BMP simulations and the direct comparisons amongst BMPs (in terms of pollutant reduction criteria) should be interpreted with caution.
Changes in the amount of sediment, TN and TP were evaluated at the outlets of the target sub-watersheds while changes in Atrazine and E. coli loads were evaluated at the NDNR gauging station 06883940 according to the NWQI project plan. The changes in pollutant loads due to BMPs' implementation were compared with the baseline scenario which is the current practice carried out in the watershed. The baseline scenario was established with the calibrated SWAT model. The percent change was calculated using Equation (6): where Load i is the load before implementing BMP and Load f is the load after implementing BMP (negative values imply reduction).

Daily E. coli Load Prediction
This study used ANN to predict daily E. coli concentrations at the gauging station. Figures S1 and S2 in the Supplementary Materials show the performance of the ANN model for training, validation and testing. The blue, green and red bars in Figure S1 (left) represent the training data, the validation data, and the testing data, respectively. The error from the neural network ranged from −0.52 logcfu/100 mL (leftmost bin) to 0.85 logcfu/100 mL (rightmost bin). The test set error and the validation set error had similar characteristics, and no significant overfitting occurred after iteration 3 where the best validation performance occurred (Figure S1 right). The training stopped when the validation error increased for six validation checks, which occurred at epoch 9. The RMSLE values for training, validation and testing at iteration 3 were all less than 0.45 (MSLE < 0.2) with an overall R 2 value of 0.83 ( Figure S2). Table S2 in the Supplementary Materials lists the parameters for SWAT model calibration of streamflow and pollutants in the BSCW. The default values of these parameters were extracted from SWAT documents and used to set the initial parameter ranges. The SWAT-CUP program with the SUFI-2 algorithm was used for calibrating the model on monthly time step with 1000 runs. Simulated sediment loss was also compared with the observed values at the gauging station. The calibration and validation results also showed good agreement between simulations and observations with acceptable statistical criteria (NSE = 0.53, RSR = 0.68, PBIAS = 17.0). Similarly for TN and TP, the calibration and validation results indicated good performance of the model. The calibration and result for Atrazine indicated good performance although the validation objective functions are slightly below the satisfaction levels. The poor E. coli calibration and validation results were probably due to the uncertainty involved in E. coli source estimation from wildlife and wastewater treatment facilities, coupled with exclusion of failing septic systems in our SWAT model. A satisfactory calibration and validation result may be achieved if detailed wildlife records, which include daily information on forage type, number and type of swine and sheep, number and type of cattle grazing, the number of days each group stayed in each pasture or feedlot for the entire study period, unrestricted cattle access to streams, the number of migrating birds, etc., were considered as model inputs.

Streamflow and Water Quality Simulation
In general, the underestimation of LOADEST peaks by our SWAT model (Figure 3) may be attributed to poor LOADEST estimates due to the unavailability of water quality data from major storm events. In a study focused on simulating sediment and nutrient loads in several gauging stations, Park and Engel [51] observed that, although regression models based on LOADEST simulated pollutant loads reasonably with reliable model coefficients, some of their results also indicated that LOADEST overestimated pollutant loads and had a bias. With nine sub-sampled water quality data sets, they observed that some regression models in LOADEST overestimated sediment load by more than ten times and overestimated nitrogen and phosphorus loads by several hundred times. However, estimates of all pollutants significantly improved with the inclusion of water quality data for maximum streamflow. Therefore they concluded that avoiding water quality data extrapolation, coupled with the availability of water quality data from storm events, were critical in annual pollutant load estimation using LOADEST.  Moreover, the underestimation of some peaks for the five pollutant loads (Figure 3) could be dependent on whether the SWAT inputs such as weather (especially precipitation) and management operations (fertilizer type, application timing, application amount, and application location) which vary spatially and temporally were true representations of reality in the watershed. For instance, regardless of how complex or detailed a model is in simulating hydrological processes, the accurate representation of precipitation over the watershed (and sub-watersheds) is important for accurate runoff and pollutant load simulation. Since SWAT uses the minimum distance from a sub-watershed centroid to the available weather stations for selecting its representative station, the PRISM weather data used in our study may not be totally representative of the reality in the sub-watersheds of BSCW. This could also have contributed to the underestimation of simulated peaks. Furthermore, in order to reduce the errors involved in predicting E. coli concentration using data-driven methods, a recent study by Abimbola et al. [52] showed that more detailed information about animal management (e.g., animal density, locations or distances from streams, and grazing patterns) needs to be coupled with hydro-meteorological variables.
Another potential reason for underestimating the pollutant loads for some storm events is the overflowing of the reservoir. The reservoir, located in the headwaters of Big Sandy Creek, has high E. coli and nutrient concentrations [53]. From a separate study, we observed that the reservoir overflowed in 2019, but there is no documentation of previous reservoir water depths and overflowing.

Effect of Individual BMP Scenarios
This study evaluated the performance of five individual BMPs using the 'best-case' scenarios where the BMPs were implemented individually over the entire target subwatersheds. The plan of the NWQI project is to implement all or combinations of these BMPs on portions of the target sub-watersheds for three years (Table 3). However, since there were many uncertainties in the future adoption of theses BMPs by agricultural producers, this made it difficult to select the exact future locations of the BMPs in the target sub-watersheds. Thus, the 'best-case' scenarios were based on the total areas of the target sub-watersheds, and the reductions of these 'best-case' scenarios were later scaled down to estimate reductions based on the acreage proposed by NWQI ( Table 3). The effectiveness of the BMPs in mitigating sediment loss, total N loss, total P loss, Atrazine loss, and E. coli loss was determined under baseline climate conditions (2004-2016, calibrated). This baseline scenario formed the basis for evaluating the individual BMP scenarios. The differences in average annual losses of pollutants associated with BMP implementation were computed for the calibration period and aggregated by BMP. Summaries of the longterm average annual percentage reductions at the target sub-watersheds and gauging station are presented in Tables 5-9. The distributions of annual reductions (box plots) and intra-annual reductions (bar charts) at the target sub-watersheds are shown in the Supplementary Materials (Figures S5-S16). Table 5. Average annual reduction (standard deviation) in loadings using no-till at the target subwatersheds. Negative values indicate an increase.

No-Tillage Scenario
When compared to the existing baseline conditions, results showed that change of tillage practices (e.g., field cultivator, tandem disk) to no-till resulted in sediment load reductions in all the hotspots (target sub-watersheds) and total phosphorus load reductions in most of the hotspots. Reductions in sediment and total phosphorus loads at subwatershed level ranged from 17% to 39% and 1% to 16% respectively, with little or no increase in total phosphorus loads in two sub-watersheds ( Table 5). Distributions of annual and intra-annual reductions are shown in the Supplementary Materials. In general, the reduction in sediment and total phosphorus loads with no-till could be due to the lower mixing efficiency (0.05) and lower mixing depth (~25 mm) of the no-till when compared with field cultivator (mixing efficiency = 0.3; depth of mixing 100 mm) and tandem-disk regular (mixing efficiency = 0.6; depth of mixing = 75 mm) in SWAT modeling. No-till results in no disturbance of the soil through tillage, and leaves previous crop residue on the soil surface. This reduces soil erosion due to surface runoff, especially during the early part of the growing season when fields are bare and more likely to be eroded [54]. Similar to our study, several SWAT studies have shown percentage reductions in sediment and total phosphorus due to the application of no-till at watershed scale [55][56][57]. Table 6. Average annual reduction (standard deviation) in loadings using crop rotation at the target sub-watersheds. Negative values indicate an increase.

Sediment
Total N Total P   77 (17) 15 (8) 29 (11) Similar to some SWAT studies showing mixed results on the effects of no-till on total nitrogen and nitrate in literature [55,58,59], our no-till BMP predicted slight load reductions in only two sub-watersheds while other sub-watersheds had slight increase in loads. Mittelstet et al. [59] also found increase in simulated nitrate loads for three hydrologic soil groups in some Nebraska watersheds. In our study, the relatively high reductions in total nitrogen and total phosphorus loads in sub=watershed 1a was because this sub-watershed had feedlots with about 25,000 cattle. Table 6 shows sediment and nutrient losses from crop rotation in the target subwatersheds. In this scenario analysis, results showed that sediment was not sensitive to crop rotation management with values ranging between −2% and 5%. In a previous study focusing on crop rotation scenarios which included continuous corn, continuous soybean, and corn-soybean rotation in Mississippi, Ni and Parajuli [60] also reported similar insensitivity of sediment to crop rotation.

Crop Rotation Scenario
The total nitrogen loss would be expected to be lower in the corn-soybean rotation BMP than that in the continuous corn baseline due to the lower amount of nitrogen fertilizer application rate to soybeans. However, since the baseline land use in all of the target subwatersheds contained about two times more corn-soybean rotation than continuous corn on the average, the conversion of continuous corn to corn-soybean rotation showed little effect in total nitrogen loads with mixed results. In addition, because total nitrogen uptake by corn and soybean may be similar [61,62], the only difference between the total nitrogen load reductions for sub-watershed 1a and other hotspots was probably due to the large number of cattle in feedlots present in sub-watershed 1a.
For total phosphorus loads, there were different magnitudes of reductions in all hotspots in the corn-soybean rotation scenario (Table 6). Compared to sub-watershed 1a (22% reduction), relatively lower total phosphorus reductions (11%) occurred in subwatersheds 1b and 1c, and small reductions (between 1 and 5%) in the remaining hotspots. The 1-5% reductions in these hotspots could be attributed to the differences in cropland and pasture areas, as well as areas with slopes greater than 6%. In the baseline scenario, hotspots with small pasture areas (<22%) had more corn-soybean rotation (>53%) than continuous corn (<23%), which resulted in small reductions when continuous corn was converted to corn-soybean rotation. On the other hand, since hotspots with larger areas with >6% slope had pasture on these slopes close to the streams (see land use map for the inset in Figure  S3 in Supplementary Materials), we could only speculate that the total phosphorus slightly reduced, probably because of phosphorus trapping and/or lower fertilizer demand in corn-soybean rotation. Distributions of annual and intra-annual reductions are shown in the Supplementary Materials.

Filter Strip Scenario
High reductions in sediment and nutrient loads were achieved by implementing filter strips in the agricultural areas on slopes greater than 6% (Table 7). Under this BMP, the highest reduction in pollutants was observed for total phosphorus with an average of 23.5%, followed by total nitrogen with an average of 20%. The higher reduction for total phosphorus may be caused by differences among nutrient sources considered in this study. Moreover, reductions were slightly higher in hotspots with less than 20% pasture (sub-watersheds 2a and 2b) than for those with pasture areas greater than 20% since most of the pastures were on slopes greater than 3% in most of the hotspots in the baseline scenario. Distributions of annual and intra-annual reductions are shown in the Supplementary Materials.
A number of studies have reported high removal rates of sediment and nutrients in surface flows through the use of filter strips [63,64]. In their study conducted in Iran, Babaei et al. [63] found filter strip to be a very effective BMP in reducing total nitrogen load (42.61%), and total phosphorus load (39.57%). In a field-scale modeling conducted by Merriman et al. [64], their results indicated that constructing filter strips was an effective BMP for reducing sediment (97% reduction), total nitrogen (13%) and total phosphorus (63%) loads.

Grassed Waterway Scenario
The highest reductions in pollution loads were achieved by implementing grassed waterways in the hotspots (Table 8). Under this BMP, the highest reduction in pollutants was observed for sediment (>91%). These high reductions were due to fact that the target hotspots already had pastures on slopes greater than 6% and close to stream reaches in the baseline scenario ( Figures S3 and S4 in the Supplementary Materials). These pastures acted as filter strips. Coupling the existence of pastures with the implementation of grassed waterways in all tributaries in the hotspots, such high reductions would be expected. Furthermore, the results indicated that the total phosphorus loads were reduced slightly more than the total nitrogen (Table 8). Distributions of annual and intra-annual reductions are shown in the Supplementary Materials. Merriman et al. [64] observed reductions of 60%, 9% and 39% for sediment, total nitrogen and total phosphorus loads, respectively, with grassed waterways. When they combined grassed waterways with filter strips, they observed load reductions of 99%, 9% and 64% for sediment, total nitrogen and total phosphorus, respectively.

Terrace Scenario
Constructing terraces across field slopes would be an effective BMP for reducing erosion in any watershed. SWAT predicted that terraces, if constructed across slopes greater than 6% in hotspots, would reduce annual sediment, total nitrogen, and total phosphorus loads by 46-77%, 1-51% and 8-58%, respectively (Table 9). Although subwatersheds 2d and 2e both had more than 60% land areas with slopes >6%, when compared to the other seven hotspots with an average of 41% area and slopes >6% (slope area in all nine hotspots ranging from 25.8% to 53.7%), the relatively lower reductions in sediment and nutrients in sub-watersheds 2d and 2e could be attributed to the bigger pasture areas that already existed on these slopes in the baseline scenario. With >6% slope area in all nine hotspots ranging from 25.8% to 53.7%, sub-watersheds 2a and 2b had an average of 33%, slightly above sub-watershed 1a (25.8%). Moreover, both sub-watersheds also had more land areas (average of 37.8%) with slopes <3% compared to other hotspots (average of 23.5%). For these two hotspots, terrace was not so effective probably because there was not enough slope gradient. The dominant flow process in both hotspots was most likely saturation excess overland flow because the field areas were relatively flat. In a field study in southwestern Iowa, Gassman et al. [65] found ten-fold reductions in sediment and nutrients loads with terrace BMP.

Atrazine Concentration at Gauging Station
Results indicated that reducing the application rates of Atrazine to 50% on all rainfed and irrigated corn in the 'best-case' scenario only reduced Atrazine concentration by 21% (from 92 µg/L to 72 µg/L) at the NDNR gauging station 06883940. This reduced concentration was still six times higher than the target concentration of 12 µg/L. As shown in Figure 4, baseline Atrazine concentrations in the reaches of sub-watersheds 2a and 2b were three times higher compared to other hotspot reaches. This was because subwatersheds 2a and 2b were headwater sub-watersheds (Figure 2), and they had the highest combined percent of continuous corn and corn-soybean, 80.2% and 79.6%, respectively. Thus, the effect of dilution was minimal in both reaches. Since the predominant land use types in BSCW are corn-soybeans and continuous corn ( Figure S3), there is a need to reduce Atrazine application in several upland sub-watersheds. With respect to the nine hotspots, there is need for more application rate reductions in sub-watersheds 2a and 2b, while in the BSCW in general, reducing or eliminating the use of Atrazine may be the way forward in order to achieve the target concentration.

E. coli Concentration at Gauging Station
Although grassed waterway had the highest E. coli reduction of 64%, followed by crop rotation (29% reduction), none of the BMPs reduced the E. coli concentrations at the gauging station to the target concentration of 112 cfu/100 mL ( Figure 5). Other BMPs had less than 10% reductions. Even though the reach of sub-watershed 1a had very high average E. coli concentration because of the feedlots with about 25,000 cattle, this high concentration was diluted by streamflow from upland reaches flowing through reach 1b and continuing downstream to reach 1c ( Figure 2). However, this dilution was not sufficient to reduce the concentration to target value. Moreover, the relatively high E. coli concentration in reach 2c did not originate from reach 3a. Rather, it was due to high concentrated streamflow from upland reaches with feedlots combining with streamflow from reach 3a before going through reach 2c (Figure 2).
Due to the uncertainty involved in E. coli source estimation from wildlife, coupled with exclusion of failing septic systems in our SWAT model, an effective manure management plan needs to be put in place based on the application amount, timing, and strategy, most especially in sub-watersheds with the greatest E. coli source contribution. This would go a long way to further reduce E. coli concentration at the gauging station to the target value. The target value can also be met by repairing failing septic tanks. Weed management alternatives to Atrazine should also be promoted for adoption in the BSCW. In field trials in major sweet corn production areas which included three US states, Arslan et al. [66] demonstrated that Atrazine-free weed management systems were comparable in performance to standard Atrazine-containing weed management systems based on input from sweet corn farmers. These systems included treatments with tembotrione either with or without interrow cultivation, and treatments with topramezone in locations where small-seeded weed species were prevalent. A published meta-analysis of existing literature by Rohr and McCoy [67] and other studies [68][69][70][71][72] showed the effects of different concentrations of Atrazine on amphibians and freshwater fish, which led to the conclusion that the adverse effects of Atrazine must be weighed against any of its benefits. It is important to note that whereas the European Union (EU) announced a ban on using Atrazine in 2003 because of unpreventable water contamination, the US Environmental Protection Agency (EPA) approved the continued use of Atrazine that same year [73].

E. coli Concentration at Gauging Station
Although grassed waterway had the highest E. coli reduction of 64%, followed by crop rotation (29% reduction), none of the BMPs reduced the E. coli concentrations at the gauging station to the target concentration of 112 cfu/100 mL ( Figure 5). Other BMPs had less than 10% reductions. Even though the reach of sub-watershed 1a had very high average E. coli concentration because of the feedlots with about 25,000 cattle, this high concentration was diluted by streamflow from upland reaches flowing through reach 1b and continuing downstream to reach 1c ( Figure 2). However, this dilution was not sufficient to reduce the concentration to target value. Moreover, the relatively high E. coli concentration in reach 2c did not originate from reach 3a. Rather, it was due to high concentrated streamflow from upland reaches with feedlots combining with streamflow from reach 3a before going through reach 2c (Figure 2). less than 10% reductions. Even though the reach of sub-watershed 1a had very high average E. coli concentration because of the feedlots with about 25,000 cattle, this high concentration was diluted by streamflow from upland reaches flowing through reach 1b and continuing downstream to reach 1c ( Figure 2). However, this dilution was not sufficient to reduce the concentration to target value. Moreover, the relatively high E. coli concentration in reach 2c did not originate from reach 3a. Rather, it was due to high concentrated streamflow from upland reaches with feedlots combining with streamflow from reach 3a before going through reach 2c (Figure 2).
Due to the uncertainty involved in E. coli source estimation from wildlife, coupled with exclusion of failing septic systems in our SWAT model, an effective manure management plan needs to be put in place based on the application amount, timing, and strategy, most especially in sub-watersheds with the greatest E. coli source contribution. This would go a long way to further reduce E. coli concentration at the gauging station to the target value. The target value can also be met by repairing failing septic tanks.

Scaling down the Best-Case Scenario BMPs
In order to examine how the BMPs would perform based on the control areas considered by the NWQI for their implementation in three years, we aggregated the nine hotspots to reflect the three major hotspots in the NWQI plan. Hotspots with numbers 1-, 2-and 3-in Figure 2 correspond to the three target sub-watersheds delineated by US Geological Survey (USGS) 12-digit hydrologic code (HUC-12). These target sub-watersheds 1, 2 and 3 are Outlet Dry Sandy Creek, Big Sandy Creek, and City of Belvidere, respectively. The aggregated reductions in these three HUC-12 sub-watersheds were calculated using the weighted average of the contributing hotspot areas. The load reductions were scaled down based on an area ratio method. Although there may not be linear relationships between increasing BMP implementation scale and Atrazine or E. coli concentration reductions at the gauging station due to the contributions from other 44 sub-watersheds, a linear relationship in the form of area ratio may be useful for sediment and nutrient loads transported with water into a reach at sub-watershed scale.
The area ratio method used in this study is similar to the drainage-area ratio method commonly used for estimating streamflow for ungauged sites. It was based on the assumption that the reduction in load for a hotspot of interest could be estimated by multiplying the ratio of the BMP implementation area for the hotspot of interest and the BMP Due to the uncertainty involved in E. coli source estimation from wildlife, coupled with exclusion of failing septic systems in our SWAT model, an effective manure management plan needs to be put in place based on the application amount, timing, and strategy, most especially in sub-watersheds with the greatest E. coli source contribution. This would go a long way to further reduce E. coli concentration at the gauging station to the target value. The target value can also be met by repairing failing septic tanks.

Scaling down the Best-Case Scenario BMPs
In order to examine how the BMPs would perform based on the control areas considered by the NWQI for their implementation in three years, we aggregated the nine hotspots to reflect the three major hotspots in the NWQI plan. Hotspots with numbers 1-, 2-and 3-in Figure 2 correspond to the three target sub-watersheds delineated by US Geological Survey (USGS) 12-digit hydrologic code (HUC-12). These target sub-watersheds 1, 2 and 3 are Outlet Dry Sandy Creek, Big Sandy Creek, and City of Belvidere, respectively. The aggregated reductions in these three HUC-12 sub-watersheds were calculated using the weighted average of the contributing hotspot areas. The load reductions were scaled down based on an area ratio method. Although there may not be linear relationships between increasing BMP implementation scale and Atrazine or E. coli concentration reductions at the gauging station due to the contributions from other 44 sub-watersheds, a linear relationship in the form of area ratio may be useful for sediment and nutrient loads transported with water into a reach at sub-watershed scale.
The area ratio method used in this study is similar to the drainage-area ratio method commonly used for estimating streamflow for ungauged sites. It was based on the assumption that the reduction in load for a hotspot of interest could be estimated by multiplying the ratio of the BMP implementation area for the hotspot of interest and the BMP area for the 'best-case' scenario in that hotspot by the load reduction for the 'best-case' scenario. To simplify the evaluation and to ensure that the estimated load reduction was an unbiased estimate of the order of magnitude of the actual load reduction, the average load reduction for the calibrated years was used rather than load reductions for individual months. It is important to note that while this area ratio method may overestimate or underestimate loads for some hotspots, it is a simple method that may be sufficient to capture the order of magnitude of load reductions. To remove all uncertainties and correct the bias in load estimations due to area ratio, a BMP simulation in SWAT would require having all the information and details about the exactness of all input parameters, including specific locations of each BMP, timing, proportion, magnitude, and adoption of the BMP by farmers.
However, with the exception of filter strip, no individual BMP was effective in achieving the targeted load reductions after down-scaling the implementation acreage (Table 10). Therefore, to achieve a 20% sediment and nutrient load reduction in the hotspots, we suggested that multiple BMPs be implemented. Overall, a combination of filter strip, grassed waterway and Atrazine rate reduction is most likely to yield measureable improvement both in the hotspots and at the gauging station. Table 10. Scaled-down load reductions from BMPs in the three HUC-12 sub-watersheds aggregated from the nine hotspots.

Conclusions
With limited funding and acreage for BMP implementation, the NWQI has set a target of 20% reduction in sediment and nutrients losses from target sub-watersheds (hotspots) in three years. With the focus of the NWQI in these hotspots, our results indicated that individual conservation practices were effective in reducing a particular pollutant load when implemented all over the hotspots (best-case scenario). Among the BMPs modeled, grassed waterway was the most effective individual BMP in reducing sediment and nutrient loads in the hotspots, and reducing E. coli concentrations at the gauging station. Reducing application rate by 50% was most effective for Atrazine reduction at the gauging station. It is noteworthy that, even in the best-case scenario, Atrazine and E. coli target concentrations at the station were not met. Besides the hotspots, a reduction or elimination of Atrazine use in upland sub-watersheds, and treatment of manure in sub-watersheds with high number of animals and feedlots, would benefit the whole BSCW.
When implementation areas were scaled down, a combination of filter strip and grassed waterway was sufficient for meeting the 20% reduction target for sediment and nutrients. In general, a combination of BMPs could yield measureable improvement in the hotspots and at the gauging station. Implementing a BMP such as filter strip may require additional management effort from farmers, and may reduce yields and profits due to high operating costs, especially in the first few years of adoption while the farmers are still learning how to successfully implement it. In addition, BMPs such as filter strips which displace cropland may create high opportunity costs for agricultural producers since there will not be profits on the plots of land with filter strips.
Moreover, future work needs to consider how to achieve these reduction objectives by exploring barriers to, and incentives for adopting and implementing the BMPs so that BMP solutions meet the target needs going forward.
Supplementary Materials: The following are available online at https://www.mdpi.com/2071-105 0/13/1/103/s1, Figure S1: Error histogram with 20 bins and best validation performance, Figure S2: Observed and ANN-simulated E. coli concentrations at the gauging station, Figure S3: Land use map of BSCW, Figure S4: Slope map of BSCW, Figure S5: Box plots showing the distributions of annual reductions (n = 13) using no-till, Figure S6: Box plots showing the distributions of annual reductions (n = 13) using crop rotation, Figure S7: Box plots showing the distributions of annual reductions (n = 13) using filter strips, Figure S8: Box plots showing the distributions of annual reductions (n = 13) using grassed waterway, Figure S9: Box plots showing the distributions of annual reductions (n = 13) using terrace, Figure S10: Box plots showing the distributions of annual reductions (n = 13) using reduced atrazine rate, Figure S11: Average monthly (growing season) reduction in loadings using no-till at the target subwatersheds, Figure S12: Average monthly (growing season) reduction in loadings using crop rotation at the target subwatersheds, Figure S13: Average monthly (growing season) reduction in loadings using filter strips at the target subwatersheds, Figure S14: Average monthly (growing season) reduction in loadings using grassed waterway at the target subwatersheds, Figure S15: Average monthly (growing season) reduction in loadings using terrace at the target subwatersheds, Figure S16: Average monthly (growing season) reduction in loadings using reduced atrazine rate at the target subwatersheds, Table S1: Management operation for land in corn/soybeans rotation and pasture management, Table S2: SWAT parameters and their final values used in calibration for the Big Sandy Creek Watershed, Nebraska SWAT model.