Quantifying the Role of Large Floods in Riverine Nutrient Loadings Using Linear Regression and Analysis of Covariance

This study analyzes the role of large river flow events in annual loads, for three constituents and for up to 32 years of daily data at multiple watersheds with different land-uses. Prior studies were mainly based on simple descriptive statistics, such as the percentage of nutrient loadings transported during several of the largest river flows, while this study uses log-regression and analysis of covariance (ANCOVA) to describe and quantify the relationships between large flow events and nutrient loadings. Regression relationships were developed to predict total annual loads based on loads exported by the largest events in a year for nitrate plus nitrite nitrogen (NO3-N + NO2-N, indicated as total oxidized nitrogen; TON), total phosphorus (TP), and suspended solids (SS) for eight watersheds in the Lake Erie and Ohio River basins. The median prediction errors for annual TON, TP, and SS loads from the top five load events for spatially aggregated watersheds were 13.2%, 18.6%, and 13.4%, respectively, which improve further on refining the spatial scales. ANCOVA suggests that the relationships between annual loads and large load events are regionally consistent. The findings outline the dominant role of large hydroclimatic events, and can help to improve the design of pollutant monitoring and agricultural conservation programs.


Introduction
Globally, increased nutrient and sediment loading (flux) from agricultural and urban activities is a primary factor in declining water quality in receiving streams, lakes, and estuaries [1,2].Excess pollutant leaching, specifically nitrogen (N) and phosphorus (P) from intensively cropped and artificially drained watersheds and other nonpoint sources, has caused the severe impairment of surface water quality worldwide.Increased enrichment of these pollutants leads to eutrophication followed by hypoxia, thereby severely disturbing ecological balances and reducing biodiversity.This renders surface water resources unsuitable for designated uses such as fishing, drinking, irrigation, and recreation, and causes human health issues such as methemoglobinemia [3].Within the United States, over the last few decades, this phenomenon has adversely affected parts of the Gulf of Mexico [4,5] and Lake Erie [6] due to pollutant leaching from the primarily agricultural Midwestern watersheds.
The detrimental impacts of hypoxic conditions in lakes and estuaries, attributed to pollutant delivery from agricultural watersheds, make the export of pollutants from such watersheds a matter of great interest.It is essential to monitor streams and rivers draining into these water bodies, in order to estimate annual pollutant delivery rates.Such frequent long-term monitoring helps to: (I) assess trends in pollutant loading over decadal timescales; (II) evaluate the efficacy of best management practices and pollutant abatement programs in these watersheds; (III) assess the impacts of climate change, urbanization, and land-use change on pollutant delivery rates; (IV) gauge the relative loading contributions of the streams to a receiving water body; and (V) provide a better understanding of pollutant export dynamics, thereby informing future strategies for reducing pollutant transport to downstream water bodies [7].
In estimating annual pollutant loads, water managers and monitoring agencies are usually concerned with improving estimation accuracy and minimizing uncertainty.However, resource constraints such as limited funding and workforce can inhibit the collection of large numbers of samples [8], leading to a dearth of long-term and complete water quality datasets.Over the past several decades, water quality sampling strategies have also changed, moving from traditional fixed-frequency water sampling (bi-weekly, monthly, quarterly) to mixed strategies, such as fixed-frequency water sampling supplemented with some storm-based sampling, to capture high flows [9,10].These approaches are time-and resource-intensive, as it is essential to characterize all flow regimes and capture most of the high load periods throughout the year in order to avoid the underestimation of annual pollutant loading values.The application of several empirical or statistically-based approaches is less data-and resource-intensive, such as the Universal Soil Loss Equation (USLE) [11], SPAtially Referenced Regression on Watershed attributes (SPARROW) [12], Geographic Information System (GIS) Pollutant Load (PLOAD) [13], or the Spreadsheet Tool for Estimating Pollutant Loads (STEPL) [14], but these approaches lack the high temporal resolution of load estimation and only estimate long-term average loadings.Therefore, there is a need to develop efficient alternative methodologies in order to estimate annual pollutant loadings that can maximize monitoring efficiencies.
Numerous studies [15][16][17][18][19][20] have reported that in various Midwestern U.S. watersheds, a few high-flow events carry significant amounts of the total annual pollutants exported over a year.Applying similar findings, Markus and Demissie [21] developed an innovative methodology to estimate annual sediment loads for 27 watersheds in Illinois by using regression relationships between high-flow event loads and total annual loads.
The present study utilizes multiyear long-term water quality datasets from 10 Midwestern watersheds to understand the role of large events in total annual pollutant export.Unlike previous studies, which focused on short-term datasets from similar watersheds and typically analyzed sediment export, this study is more comprehensive as it assesses total oxidized nitrogen (TON), total phosphorus (TP), and suspended solids (SS) load exports from watersheds with varied land-uses and sizes over periods ranging from seven to 29 years.A modified pollutant-specific baseflow separation technique was developed to identify high-flow events in a water year by comparing the local-minimum and recursive digital base-flow filters obtained through the Web-based Hydrograph Analysis Tool (WHAT) [22].Pollutant loads exported during these events were computed and correlated to total annual loads exported using regression relationships at different scales of spatial aggregation (separate watersheds, all watersheds together).The analysis of covariance (ANCOVA) was used to compare regression lines from different watersheds to evaluate whether they can be regarded as statistically equivalent or different.To evaluate the spatial transferability and the practical applicability of this method, the developed regressions were applied to watersheds in the same geographical region with similar land-use patterns.

Site and Data Descriptions
The data used in this study were collected by the Water Quality Lab (WQL), located at the National Center for Water Quality Research, Heidelberg College, Ohio.The sampling program was designed to collect daily samples during low flows with additional samples collected during storm events.Data were available for the eight watersheds in the Lake Erie and Ohio River basins (Ohio watersheds), for monitoring periods ranging from 8 to 32 years.Data preprocessing was conducted to eliminate water years with considerable missing concentration and flow data.In-depth details about the sampling program and chemical analysis are described in [23][24][25].The study analyzed three pollutants: nitrate plus nitrite nitrogen (NO 3 -N + NO 2 -N; hereafter total oxidized nitrogen, TON), suspended solids (SS), and total phosphorus as P (TP).Daily monitored flow and TON concentration data were also available from the Illinois Environmental Protection Agency (IEPA), spanning six years for the Upper Sangamon watershed and 11 years for the Vermilion watershed in Illinois.Six watersheds, namely Cuyahoga, Grand, Maumee, Raisin, Sandusky, and Vermilion, drain into Lake Erie.Two watersheds, Great Miami and Muskingum, drain southwards into the Ohio River basin, and two watersheds in Illinois, Vermilion, and Upper Sangamon, drain into the Upper Mississippi River basin as we can see from Figure 1.
Sustainability 2018, 10, 2876 3 of 19 practical applicability of this method, the developed regressions were applied to watersheds in the same geographical region with similar land-use patterns.

Site and Data Descriptions
The data used in this study were collected by the Water Quality Lab (WQL), located at the National Center for Water Quality Research, Heidelberg College, Ohio.The sampling program was designed to collect daily samples during low flows with additional samples collected during storm events.Data were available for the eight watersheds in the Lake Erie and Ohio River basins (Ohio watersheds), for monitoring periods ranging from 8 to 32 years.Data preprocessing was conducted to eliminate water years with considerable missing concentration and flow data.In-depth details about the sampling program and chemical analysis are described in [23][24][25].The study analyzed three pollutants: nitrate plus nitrite nitrogen (NO3-N + NO2-N; hereafter total oxidized nitrogen, TON), suspended solids (SS), and total phosphorus as P (TP).Daily monitored flow and TON concentration data were also available from the Illinois Environmental Protection Agency (IEPA), spanning six years for the Upper Sangamon watershed and 11 years for the Vermilion watershed in Illinois.Six watersheds, namely Cuyahoga, Grand, Maumee, Raisin, Sandusky, and Vermilion, drain into Lake Erie.Two watersheds, Great Miami and Muskingum, drain southwards into the Ohio River basin, and two watersheds in Illinois, Vermilion, and Upper Sangamon, drain into the Upper Mississippi River basin as we can see from Figure 1.Most of these watersheds have primarily agricultural land-use, except for Cuyahoga, which is approximately 47% urbanized and Grand, which is approximately 52% forested in Table 1 [26].Most of this agricultural land is given over to row-cropping with rotational schedules mainly composed of corn, soybean, and some wheat crops.Large parts of these watersheds are characterized by poorly draining but fertile soils and therefore have to be artificially drained using subsurface tile drain systems to make them conducive to agriculture.Muskingum is the largest watershed analyzed in this study, and drains an area of approximately 19,200 km 2 upstream of the sampling station.Among the watersheds draining into the Lake Erie basin, Maumee is the largest, with approximately 16,400 km 2 of watershed area upstream of the sampling station.Vermilion is the smallest watershed analyzed here, with a drainage area of approximately 700 km 2 [27].The study watersheds are strongly affected by subsurface tile drainage systems.The tile drainage system is general in northwest Ohio, eastern Indiana, and southern Michigan.The Cuyahoga and Grand watersheds contain higher precipitation and runoff than the other watersheds because of the lake-effect [26].Most of these watersheds have primarily agricultural land-use, except for Cuyahoga, which is approximately 47% urbanized and Grand, which is approximately 52% forested in Table 1 [26].Most of this agricultural land is given over to row-cropping with rotational schedules mainly composed of corn, soybean, and some wheat crops.Large parts of these watersheds are characterized by poorly draining but fertile soils and therefore have to be artificially drained using subsurface tile drain systems to make them conducive to agriculture.Muskingum is the largest watershed analyzed in this study, and drains an area of approximately 19,200 km 2 upstream of the sampling station.Among the watersheds draining into the Lake Erie basin, Maumee is the largest, with approximately 16,400 km 2 of watershed area upstream of the sampling station.Vermilion is the smallest watershed analyzed here, with a drainage area of approximately 700 km 2 [27].The study watersheds are strongly affected by subsurface tile drainage systems.The tile drainage system is general in northwest Ohio, eastern Indiana, and southern Michigan.The Cuyahoga and Grand watersheds contain higher precipitation and runoff than the other watersheds because of the lake-effect [26].

Separation of Large Events
To assess the role of large events in predicting annual pollutant loads, it is necessary to separate them from an annual hydrograph.This is even more vital in the Midwestern watersheds, as most of the annual loads are exported during a few large events.Given the large drainage areas of the watersheds studied, it is easy to identify the beginning of a large event, as a precipitation event leads to a steep rise in the hydrograph.Identifying the exact end of a large event is more complicated in these watersheds, as the receding limb of the hydrograph is generally more gradual.Additionally, after a large event, the flow stabilizes at a higher value than that at the beginning of the event.There is some overlap of closely occurring events, with two or more hydrograph peaks separated by a flow value between peak-and base-flow levels.Similarly to Richards et al. [25], snowmelt-driven events were treated as precipitation events because the flow volumes for both kinds of events are similar.
A combination of two base-flow separation filters from the WHAT were used to identify the beginning and end of large events.The points of intersection of a local-minimum filter curve and a recursive-digital filter curve defined the beginning and end of a high-flow event.The local-minimum filter linearly connects local minima to separate the base-flow from a hydrograph, while the recursive-digital filter separates base-flow and direct runoff similarly to separating high-and low-frequency signals in signal processing.For other intermediate local minima points, the local minimum value was compared to the recursive-digital filter value.This ratio was defined as the baseflow filter ratio (BFR) and was used to further separate an identified large event into smaller events based on the differences in pollutant export characteristics.If the BFR exceeded a set threshold at any local minima point in an identified large event, the event was further split at that point into two smaller events.Setting a high threshold for the BFR thereby split an identified large event containing multiple peaks into separate events; conversely, setting a low BFR did not split the hydrograph with multiple peaks into separate events [22].
In the watersheds analyzed, SS and TP were primarily exported through overland flow [25,28], and thereby their concentration peaks coincided with or just preceded the hydrograph peaks.On the other hand, TON was exported through subsurface flow [25,[28][29][30][31]. Thus, TON concentration peaks trailed the hydrograph peaks by an order of days, depending on the size and land-use characteristics of the watershed.The delay in TON concentration peaks is likely caused by dilution, as most of the increase in flow post-precipitation is contributed by surface flow.In addition, subsurface drainage system outlets are submerged due to high water levels in drainage ditches, further slowing TON export.For fast-response pollutants such as SS and TP, a high BFR of 1.75 was used, which very closely identified separate large events, as each peak in the hydrograph was separated and interpreted as a different event.For slow-response pollutants such as TON, a lower BFR of 1.25 closely identified separate large events, as events occurring in close succession (two or more close peaks) were considered as one to better characterize large event TON loads.

Predictability of Annual Loads
To test the predictability of annual pollutant loads using large events, actual pollutant and watershed specific annual loads were first computed using a simple numeric integration approach.In cases where nearly continuous data are available, a numeric integration approach is the most desirable as it is free from statistical assumptions of normality and is simple to implement [25].Similarly to previous studies [16,25,[32][33][34][35][36], annual loads were computed as: where C is the pollutant concentration, Q is the flow, and dt is the duration represented by the sample.The duration associated with each sample comprised half the time between the preceding and present sample plus half the time between the present and succeeding sample.
Similarly, pollutant loads were also calculated for each large event separated from an annual hydrograph.These large events for each water year (1 October-30 September) and watershed were then ranked based on decreasing loads exported during these events (load events), for all three pollutants respectively.Further analysis was performed to check whether the top load events were also the events with the highest peak flows.Markus and Demissie [21] reported that annual sediment loads and loads exported during the largest events based on peak flows have a high correlation.They also established that annual loads could be predicted based on large load events using regression relationships.In this study, regression relationships were developed for annual TON, SS, and TP loads and loads exported during the top few load events.The following equation describes the relationship between annual loads and top load events: where L a and L le are pollutant loads exported annually during a water year and cumulatively during the top load events in a year, respectively (le = 1, 2, 3, 4, 5, . . ., n); and a and b are regression parameters.The load values were normalized per unit watershed area for easier comparison.They were log-transformed in order to normalize residual variances and reduce skewness resulting from extreme values [37][38][39].Distinct regression relationships were developed for each combination of watershed and pollutant.In addition, regression relationships were developed by spatially aggregating data for a given pollutant from all eight Ohio watersheds and similar watersheds, for annual loads and the top few load events in a water year.

Similarity of Predictions
To test the similarity between regression lines developed for different watersheds and pollutants, ANCOVA was used [21,40].The purpose was to examine the similarity between linear regressions of annual pollutant loads and top load events among watersheds.These regressions can differ in slope, intercept, and residual variances.ANCOVA was used to compare the following model: where X is an independent variable, high flow loads in this case; Y is a dependent variable, annual load; i is watershed 1, watershed 2, and all watersheds together; j represents the observations, with pollutants being evaluated; α is the intercept; β is the slope; and ε represents the residuals.This approach compares the slope, intercept, and residual variances for two regression equations based on a two-tail F-test.If heterogeneous slopes, intercepts, or residual variances are found, then the regression lines are regarded as statistically different.The slope represents a rate of increase in annual loads based on the top load events.For similar slopes, the intercept generally represents the overall magnitude of annual loads compared to the top high-flow event loads.Regression relationships from watersheds having statistically similar regression slopes and intercepts can possibly be used interchangeably to predict annual loads based on top load events, or combined to produce statistically stronger relationships.If the regression lines were similar, then regional water quality monitoring would likely be beneficial for new monitoring sites for which there is limited data.Secondly, it also helps identify watersheds that account for similar contributions of top load events to the annual loads for a given pollutant.This would aid in designing uniform sampling strategies for such watersheds and help to identify representative watersheds for long-term continuous monitoring [41].

Applicability of Regression Lines in Other Regions
The proposed methodology for estimating annual pollutant loads by only using the loads exported in large events would be a useful tool if developed regression relationships can be successfully applied to physiologically and climatically similar watersheds in the current year, rather than using data collected on past events.To evaluate whether the developed regression relationships were applicable to watersheds with limited water quality data, regression relationships from all Ohio watersheds were used to estimate loads for the Upper Sangamon and Vermilion watersheds in central Illinois.Daily flow and TON concentration data for a number of years were available for both these Illinois watersheds, which were used as a benchmark to evaluate the accuracy of load prediction by just using loads exported in the largest events and regression relationships developed using the data from similar Ohio watersheds.To separate all major hydrographic events in a year, a flow threshold was identified based on the historical flow records for both Illinois watersheds.Consequently, the separated large events were ranked based on loads exported, and the cumulative loads exported in the top few load events were used as an input for regression relationships developed from Ohio watersheds to predict annual TON loads.

Large Event Characteristics
An analysis of the contributions of top load events to the annual loads for all combinations of pollutants, water years, and watersheds indicated that, typically, the five biggest load events cumulatively carried a large percentage of the annual loads.These percentages declined rapidly after the top six to seven events. Figure 2 shows the contribution of the top 25 load events to annual loads for a characteristic agricultural watershed (Maumee) and an urbanized watershed (Cuyahoga).This was substantial for all pollutants in Maumee and other agricultural watersheds.In contrast, for the urbanized watershed Cuyahoga, the annual loads were more evenly distributed amongst the top 10 load events.Overall, the top five load events were chosen for predicting annual loads, as in most cases, they carried a large percentage of annual loads.For some scenarios, the contributions from the 6th-10th largest events were also substantial, but such cases were infrequent.Therefore, the total cumulative loads exported by the top one to five load events each year were computed for all watersheds and pollutants.Data from all Ohio watersheds, excluding Cuyahoga, indicated that, on average, 56% of total annual TON loads were exported during the top five load events, which lasted for approximately 17% of the days in a year.Similarly, data from all Ohio watersheds indicated that, on average, for TP and SS, the top five load events carried 57% and 67% of annual loads during just 15% and 14% of the days in a year, respectively.For TP, the average duration of the top load event ranged from just eight days (2%) for Grand to approximately 18 days (5%) for Muskingum, carrying a minimum of 15% of annual TP load for Cuyahoga and a maximum of 28% for Sandusky.The top five load events carried a minimum of approximately 41% and a maximum of approximately 71% of the annual TP load for Cuyahoga and Sandusky, respectively, jointly lasting for a minimum of 36 days (10%) for Cuyahoga and a maximum of 76 days (21%) for Sandusky.As expected, the results were very similar for SS loads, with the top load event lasting from eight days (2%) to 19 days (5%) for Cuyahoga and Muskingum, respectively, carrying a minimum of 22% of annual SS load for Cuyahoga and a maximum of 31% for Sandusky.The top five load events carried a minimum of approximately 57% and a maximum of approximately 75% of annual SS load for Cuyahoga and Sandusky, respectively, jointly lasting for a minimum of 34 days (9%) for Cuyahoga and a maximum of 72 days (20%) for Sandusky.
Among the three pollutants analyzed, the least variance among the loads contributed by the top one to five load events was observed for TON export, followed by TP and SS.Additionally, the top load events lasted longer while carrying a smaller percentage of annual TON loads than for TP and SS loads.Among the watersheds, for TON, Cuyahoga had the smallest variance, while Sandusky and Vermilion had the largest variance, in top load event contributions.For TP and SS, Great Miami had the least variance, while Maumee and Vermilion showed the greatest variance, in the contributions of the top one to five load events towards annual loads.
Figure 3 shows the average percentage of common top events based on total loads exported and on peak flows for all the Ohio watersheds over their complete monitoring durations.The analysis was performed for both BFR values used in this study, which were 1.25 for TON and 1.75 for SS and TP.The size of the marker indicates the watershed size (not to scale), and shows that, for most of the agricultural watersheds, more than 75% of the top five load events matched the top five events based on peak flows over the whole year.This value increased to more than 80% when the analysis was limited to the winter and spring months, when, typically, the highest pollutant export occurs in the Midwestern watersheds [27].Specifically for TON, the average duration of the top load event ranged from approximately 14 days (4% of the number of days annually) to 26 days (7%) for Cuyahoga and Muskingum, respectively.On average, the minimum contribution of the top load event was 7% of the annual TON load for Cuyahoga, and the maximum contribution was approximately 23% for Raisin.The top five load events combined lasted for a minimum of 60 days (16%) for Maumee and a maximum of approximately 90 days (25%) for Sandusky, carrying a minimum of 24% of annual TON load for Cuyahoga and a maximum of 67% for Raisin.
For TP, the average duration of the top load event ranged from just eight days (2%) for Grand to approximately 18 days (5%) for Muskingum, carrying a minimum of 15% of annual TP load for Cuyahoga and a maximum of 28% for Sandusky.The top five load events carried a minimum of approximately 41% and a maximum of approximately 71% of the annual TP load for Cuyahoga and Sandusky, respectively, jointly lasting for a minimum of 36 days (10%) for Cuyahoga and a maximum of 76 days (21%) for Sandusky.As expected, the results were very similar for SS loads, with the top load event lasting from eight days (2%) to 19 days (5%) for Cuyahoga and Muskingum, respectively, carrying a minimum of 22% of annual SS load for Cuyahoga and a maximum of 31% for Sandusky.The top five load events carried a minimum of approximately 57% and a maximum of approximately 75% of annual SS load for Cuyahoga and Sandusky, respectively, jointly lasting for a minimum of 34 days (9%) for Cuyahoga and a maximum of 72 days (20%) for Sandusky.
Among the three pollutants analyzed, the least variance among the loads contributed by the top one to five load events was observed for TON export, followed by TP and SS.Additionally, the top load events lasted longer while carrying a smaller percentage of annual TON loads than for TP and SS loads.Among the watersheds, for TON, Cuyahoga had the smallest variance, while Sandusky and Vermilion had the largest variance, in top load event contributions.For TP and SS, Great Miami had the least variance, while Maumee and Vermilion showed the greatest variance, in the contributions of the top one to five load events towards annual loads.
Figure 3 shows the average percentage of common top events based on total loads exported and on peak flows for all the Ohio watersheds over their complete monitoring durations.The analysis was performed for both BFR values used in this study, which were 1.25 for TON and 1.75 for SS and TP.The size of the marker indicates the watershed size (not to scale), and shows that, for most of the agricultural watersheds, more than 75% of the top five load events matched the top five events based on peak flows over the whole year.This value increased to more than 80% when the analysis was limited to the winter and spring months, when, typically, the highest pollutant export occurs in the Midwestern watersheds [27].

Annual Load Predictions
For all of the Ohio watersheds, Figure 4 shows the total annual TON loads plotted against the loads contributed by the top five load events.The first panel of Figure 4 shows the individual regression lines for each watershed, while the second panel shows spatially aggregated regression lines for the contributions of the top one and top five load events to annual loads.All loading values were normalized by watershed area (yields).It is evident from the figure that, for TON, the Cuyahoga watershed displays a distinctly different relationship to that seen in the other watersheds.Regression lines from the other watersheds align closely to one another, having visibly similar slopes and residuals despite their large differences in total annual TON loads and top five event loads.The regression residuals can be caused by seasonal variations in antecedent hydroclimatic conditions, as well as agricultural and other management practices.

Annual Load Predictions
For all of the Ohio watersheds, Figure 4 shows the total annual TON loads plotted against the loads contributed by the top five load events.The first panel of Figure 4 shows the individual regression lines for each watershed, while the second panel shows spatially aggregated regression lines for the contributions of the top one and top five load events to annual loads.All loading values were normalized by watershed area (yields).It is evident from the figure that, for TON, the Cuyahoga watershed displays a distinctly different relationship to that seen in the other watersheds.Regression lines from the other watersheds align closely to one another, having visibly similar slopes and residuals despite their large differences in total annual TON loads and top five event loads.The regression residuals can be caused by seasonal variations in antecedent hydroclimatic conditions, as well as agricultural and other management practices.
Among the Ohio watersheds, the urbanized Cuyahoga had the highest point-source inputs to total riverine loads, at approximately 33% [26], and the lowest contributions to annual TON loads from the top load events.Climate-and land-management-driven variations in annual loads and their corresponding event contributions were seen in most watersheds, with the latter being particularly important for Maumee and Sandusky, resulting in the highest variation in annual TON loads.Regression lines were also developed to describe the contributions of the top two, three, and four load events toward annual loads, similarly to those for the top one and top five events, as shown in the second panel of the figure, but were omitted for clearer representation.Two regression lines were created for the contributions of the top five load events: one with, and the other without, data from the urban Cuyahoga watershed.A correlation with an adjusted R 2 of 0.77 was observed for spatially aggregated data from all watersheds including Cuyahoga, which further increased to 0.91 when data from Cuyahoga were excluded.
watershed displays a distinctly different relationship to that seen in the other watersheds.Regression lines from the other watersheds align closely to one another, having visibly similar slopes and residuals despite their large differences in total annual TON loads and top five event loads.The regression residuals can be caused by seasonal variations in antecedent hydroclimatic conditions, as well as agricultural and other management practices.Among the Ohio watersheds, the urbanized Cuyahoga had the highest point-source inputs to total riverine loads, at approximately 33% [26], and the lowest contributions to annual TON loads from the top load events.Climate-and land-management-driven variations in annual loads and their Similarly, Figure 5 shows the regression lines developed for TP and SS exports from these watersheds.In the case of TP export, individual regression lines from all watersheds lie within the same region, with potential differences in slopes and intercepts.Most watersheds showed similar contributions from the top load events towards annual loads, except for the Raisin watershed, which had the lowest contributions from the top load events.However, it had the strongest watershed-scale regression relationship for TP export.Over the monitoring period, Raisin showed relatively consistent total annual and top event TP loads, possibly due to the presence of impoundments in the upper watershed.These impoundments tend to slow the river flow, acting as sediment sinks, thereby altering the amount and timing of TP load exports from Raisin [42].Impoundments also possibly minimize the impacts of any land management and nutrient abatement programs on the relative contributions of the top load event loads to annual loads, leading to robust regression relationships.corresponding event contributions were seen in most watersheds, with the latter being particularly important for Maumee and Sandusky, resulting in the highest variation in annual TON loads.Regression lines were also developed to describe the contributions of the top two, three, and four load events toward annual loads, similarly to those for the top one and top five events, as shown in the second panel of the figure, but were omitted for clearer representation.Two regression lines were created for the contributions of the top five load events: one with, and the other without, data from the urban Cuyahoga watershed.A correlation with an adjusted R 2 of 0.77 was observed for spatially aggregated data from all watersheds including Cuyahoga, which further increased to 0.91 when data from Cuyahoga were excluded.Similarly, Figure 5 shows the regression lines developed for TP and SS exports from these watersheds.In the case of TP export, individual regression lines from all watersheds lie within the same region, with potential differences in slopes and intercepts.Most watersheds showed similar contributions from the top load events towards annual loads, except for the Raisin watershed, which had the lowest contributions from the top load events.However, it had the strongest watershed-scale regression relationship for TP export.Over the monitoring period, Raisin showed relatively consistent total annual and top event TP loads, possibly due to the presence of impoundments in the upper watershed.These impoundments tend to slow the river flow, acting as sediment sinks, thereby altering the amount and timing of TP load exports from Raisin [42].Impoundments also possibly minimize the impacts of any land management and nutrient abatement programs on the relative contributions of the top load event loads to annual loads, leading to robust regression relationships.Unlike for TON, TP export from the Cuyahoga watershed was more closely aligned with the exports from other watersheds, and was therefore included to create the spatially aggregated regression lines.Climatic variations, along with watershed-scale changes in phosphorus management policies coupled with phosphorus abatement programs implemented in these watersheds during the 1970s and 1980s, potentially accounted for most of the variation in TP export during the monitoring periods.For spatially aggregated regression lines, the adjusted R 2 values for Unlike for TON, TP export from the Cuyahoga watershed was more closely aligned with the exports from other watersheds, and was therefore included to create the spatially aggregated regression lines.Climatic variations, along with watershed-scale changes in phosphorus management policies coupled with phosphorus abatement programs implemented in these watersheds during the 1970s and 1980s, potentially accounted for most of the variation in TP export during the monitoring periods.For spatially aggregated regression lines, the adjusted R 2 values for loads exported during top load events and total annual TP loads were high, increasing from 0.62 to 0.82 for the top one to the top five load events, respectively.The strength of the regression relationship for TP loads for predominantly forested Grand River was comparable to those of other agricultural watersheds.
For SS export, individual regression lines are shown in Figure 6.Visually, these were the most similar among the three monitored pollutants.Consistent patterns were observed for SS export during top load events.It can be speculated that the absence of long-term trends in SS concentration in some of the Ohio watersheds between the years 1975 and 1990 [24], possibly led to fewer variations in both annual loads and top load event loads.This resulted in the better predictability of SS.For spatially aggregated regression lines, a high adjusted R 2 , increasing from 0.79 to 0.94, was observed for annual SS load contributions from the top one and five load events, respectively.These adjusted R 2 values for the spatially aggregated regressions were the highest among the three monitored pollutants.
Sustainability 2018, 10, 2876 10 of in some of the Ohio watersheds between the years 1975 and 1990 [24], possibly led to fewer variations in both annual loads and top load event loads.This resulted in the better predictability of SS.For spatially aggregated regression lines, a high adjusted R 2 , increasing from 0.79 to 0.94, was observed for annual SS load contributions from the top one and five load events, respectively.These adjusted R 2 values for the spatially aggregated regressions were the highest among the three monitored pollutants.
(a) (b) Most of the predominantly agricultural watersheds showed strong regression relationships for predicting annual SS loads.However, in a few of the watersheds, changes in land management practices possibly affected the performance of the regression relationships.For example, increased conservation tillage and changes in other agricultural practices were implemented in the Maumee and Sandusky watersheds in the 1990s.These changes have been linked to a decline in SS loads exported from these watersheds [25], and possibly caused greater variation in predicted annual loads as compared to other agricultural watersheds, over the span of the few decades assessed in this study.Cuyahoga, which is the least agricultural and most urbanized of the watersheds, had the weakest watershed-scale regression relationship for predicting annual SS loads from top load events.Despite having the lowest percentage of agricultural land, Cuyahoga had one of highest total annual SS load contributions to Lake Erie.This hydrological behavior, coupled with the least agricultural and most urbanized land cover, potentially led to a weak regression relationship with the highest maximum relative error.

ANCOVA
ANCOVA [40] was used to compare the regression lines developed for total annual load estimates using loads exported during the top five load events for all eight Ohio watersheds.Furthermore, for each pollutant, these individual watershed-specific regression lines were compared with a spatially aggregated regional regression line, for similarity.ANCOVA analysis was performed to gauge the statistical similarity between linear regression slopes, intercepts, and residual variances Most of the predominantly agricultural watersheds showed strong regression relationships for predicting annual SS loads.However, in a few of the watersheds, changes in land management practices possibly affected the performance of the regression relationships.For example, increased conservation tillage and changes in other agricultural practices were implemented in the Maumee and Sandusky watersheds in the 1990s.These changes have been linked to a decline in SS loads exported from these watersheds [25], and possibly caused greater variation in predicted annual loads as compared to other agricultural watersheds, over the span of the few decades assessed in this study.Cuyahoga, which is the least agricultural and most urbanized of the watersheds, had the weakest watershed-scale regression relationship for predicting annual SS loads from top load events.Despite having the lowest percentage of agricultural land, Cuyahoga had one of highest total annual SS load contributions to Lake Erie.This hydrological behavior, coupled with the least agricultural and most urbanized land cover, potentially led to a weak regression relationship with the highest maximum relative error.

ANCOVA
ANCOVA [40] was used to compare the regression lines developed for total annual load estimates using loads exported during the top five load events for all eight Ohio watersheds.Furthermore, for each pollutant, these individual watershed-specific regression lines were compared with a spatially aggregated regional regression line, for similarity.ANCOVA analysis was performed to gauge the statistical similarity between linear regression slopes, intercepts, and residual variances at α = 0.05 and α = 0.01 significance levels, and the results are shown in Table 2. Watershed names indicate watershed-specific regression lines, while "ALL" in Table 2 indicates spatially aggregated data from all the watersheds.A comparison of watershed A and watershed B would be similar to a comparison of watershed B and watershed A for slope and intercept, while it would be different for residual variances.This difference arises because the slope and intercept comparisons between two regression lines employ pooled data (which would be similar in both cases), whereas comparisons of residual variance employ individual data and their associated degrees of freedom.
Table 2. ANCOVA-based statistical comparison for annual and top five TON, TP, and SS load event regression relationships, for individual watershed and spatially aggregated data.R, S, and I indicate residual variances, slope, and intercept, respectively.Significant differences at both α = 0.05 and α = 0.01 are indicated by a value of 2, while differences significant at α = 0.05 and insignificant at α = 0.01 are indicated by a value of 1.A value of 0 indicates no significant difference between two lines at both α = 0.05 and α = 0.01.

ANCOVA
Cu Gr Gm Ma Mu Ra Sa Ve ALL
For TP, the Grand, Maumee, Muskingum, and Vermilion watersheds had regression lines that were not significantly different from the spatially aggregated regional regression line.Among the individual watersheds, a few groups were found to have statistically similar regression lines, such as Great Miami and Grand; and Maumee and Sandusky.Among the three pollutants studied, the regression lines for SS showed the fewest statistical differences.Five out of the eight watersheds had similar regression lines to that of the spatially aggregated regional trend.A few groups of watersheds were found to have statistically similar regression lines, for example, Maumee, Raisin, and Sandusky.

Error Analysis
The various regression lines developed in the study were used to predict total annual loads, and the distributions of relative errors in predicting annual loads are presented using box and whisker plots.Regression analysis was carried out for different levels of watershed aggregation, ranging from individual watersheds to aggregation by similarity, for example: by similar land use; or by proximity, such as all watersheds in the Lake Erie and Ohio River basins (Ohio watersheds).Firstly, for all three pollutants, the spatially aggregated regional regression line using all eight Ohio watersheds was used to predict annual loads from the top one to five load events.The resulting TON predictions are shown in Figure 7. Secondly, based on the ANCOVA, watersheds with statistically similar regressions (predominantly agricultural watersheds) were grouped together to create aggregated regressions for just the top five load events, and these grouped regressions were used to predict annual loads.Thirdly, individual watershed-specific regressions for the top five load events to predict annual loads were used to calculate prediction errors (Figure 8). Figure 7 shows the relative errors in annual loads predicted using the spatially aggregated regression lines.Relative errors were defined as [35]: Sustainability 2018, 10, 2876 12 of 19

Error Analysis
The various regression lines developed in the study were used to predict total annual loads, and the distributions of relative errors in predicting annual loads are presented using box and whisker plots.Regression analysis was carried out for different levels of watershed aggregation, ranging from individual watersheds to aggregation by similarity, for example: by similar land use; or by proximity, such as all watersheds in the Lake Erie and Ohio River basins (Ohio watersheds).Firstly, for all three pollutants, the spatially aggregated regional regression line using all eight Ohio watersheds was used to predict annual loads from the top one to five load events.The resulting TON predictions are shown in Figure 7. Secondly, based on the ANCOVA, watersheds with statistically similar regressions (predominantly agricultural watersheds) were grouped together to create aggregated regressions for just the top five load events, and these grouped regressions were used to predict annual loads.Thirdly, individual watershed-specific regressions for the top five load events to predict annual loads were used to calculate prediction errors (Figure 8). Figure 7 shows the relative errors in annual loads predicted using the spatially aggregated regression lines.Relative errors were defined as [35]: The median relative errors for annual load predictions from the top one, two, three, and four load events for TON were 24.3%, 23.6%, 23.0%, and 23.1%, respectively, for all Ohio watersheds.For predictions using the top five load events, two regression relationships were developed, consisting of one including and one excluding the results from Cuyahoga.The median relative errors for these regressions were 21.9% and 13.1%, respectively.The regression relationship for the top five load events, excluding data from Cuyahoga, considerably improved annual load predictions for the remaining watersheds, as was evident by the reduction in median relative error.Its maximum relative error and standard deviations were also the smallest, at 52.3% and 11.5%, respectively.
The median relative errors for annual load predictions were lower for TP and SS compared to TON.For TP, the median relative errors for annual load predictions for the top one to five load events were 27.5%, 23.6%, 21.9%, 20.1%, and 18.6%, respectively.The maximum relative errors and standard deviations also declined for the top one and top five load events: from 117% to 56.9%, and 24.3% to 13.2%, respectively.SS annual load predictions had the lowest relative errors among the three pollutants, with median values of 20.8%, 19.0%, 17.5%, 15.4%, and 13.5% for the top one to five load events.The maximum relative errors and standard deviations for SS annual load predictions declined from 118% to 46%, and from 23.4% to 10.6%, for the top one and top five load events, respectively.Overall, the results showed that the top load events for TON tend to last longer, carry a smaller percentage of annual loads, and have higher prediction errors than TP and SS load events.Based on the ANCOVA, for every pollutant, a few watersheds were found to have statistically similar regression relationships for annual load predictions using loads exported during the top five load events.Data from these similar watersheds were aggregated, and new regressions were developed for annual load predictions from the total loads exported during the top five load events.The distributions of the relative errors in the predicted loads using these new regressions are presented using box and whisker plots in Figure 8.

Figure 8.
Relative error distributions for observed annual loads and annual loads estimated using the top five nitrate plus nitrite nitrogen (TON), total phosphorus (TP), and suspended solids (SS) load events for spatially aggregated data from a few watersheds with statistically similar regression relationships between annual loads and top load events, followed by data from individual watersheds.A cross (+) symbol beyond the length of whiskers indicates outliers, and a small square symbol within the boxes indicates the mean of the distributions.
Data from the four agricultural watersheds (Maumee, Muskingum, Raisin, and Sandusky) were aggregated as their individual watershed-specific regressions for TON were statistically similar.The median and maximum relative errors for this new regression were 12.0% and 45.3%, with a standard deviation of approximately 10.9%.The median and maximum relative errors were 13.2% and 52.2%, respectively, for spatially aggregated regression from all watersheds, excluding Cuyahoga.For TP, a Figure 8. Relative error distributions for observed annual loads and annual loads estimated using the top five nitrate plus nitrite nitrogen (TON), total phosphorus (TP), and suspended solids (SS) load events for spatially aggregated data from a few watersheds with statistically similar regression relationships between annual loads and top load events, followed by data from individual watersheds.A cross (+) symbol beyond the length of whiskers indicates outliers, and a small square symbol within the boxes indicates the mean of the distributions.
The median relative errors for annual load predictions from the top one, two, three, and four load events for TON were 24.3%, 23.6%, 23.0%, and 23.1%, respectively, for all Ohio watersheds.For predictions using the top five load events, two regression relationships were developed, consisting of one including and one excluding the results from Cuyahoga.The median relative errors for these regressions were 21.9% and 13.1%, respectively.The regression relationship for the top five load events, excluding data from Cuyahoga, considerably improved annual load predictions for the remaining watersheds, as was evident by the reduction in median relative error.Its maximum relative error and standard deviations were also the smallest, at 52.3% and 11.5%, respectively.
The median relative errors for annual load predictions were lower for TP and SS compared to TON.For TP, the median relative errors for annual load predictions for the top one to five load events were 27.5%, 23.6%, 21.9%, 20.1%, and 18.6%, respectively.The maximum relative errors and standard deviations also declined for the top one and top five load events: from 117% to 56.9%, and 24.3% to 13.2%, respectively.SS annual load predictions had the lowest relative errors among the three pollutants, with median values of 20.8%, 19.0%, 17.5%, 15.4%, and 13.5% for the top one to five load events.The maximum relative errors and standard deviations for SS annual load predictions declined from 118% to 46%, and from 23.4% to 10.6%, for the top one and top five load events, respectively.Overall, the results showed that the top load events for TON tend to last longer, carry a smaller percentage of annual loads, and have higher prediction errors than TP and SS load events.
Based on the ANCOVA, for every pollutant, a few watersheds were found to have statistically similar regression relationships for annual load predictions using loads exported during the top five load events.Data from these similar watersheds were aggregated, and new regressions were developed for annual load predictions from the total loads exported during the top five load events.The distributions of the relative errors in the predicted loads using these new regressions are presented using box and whisker plots in Figure 8.
Data from the four agricultural watersheds (Maumee, Muskingum, Raisin, and Sandusky) were aggregated as their individual watershed-specific regressions for TON were statistically similar.The median and maximum relative errors for this new regression were 12.0% and 45.3%, with a standard deviation of approximately 10.9%.The median and maximum relative errors were 13.2% and 52.2%, respectively, for spatially aggregated regression from all watersheds, excluding Cuyahoga.For TP, a statistical similarity was found in a few watersheds, such as Grand and Great Miami.The median and maximum relative errors for this regression reduced from 18.6% to 12.0%, and from 56.9% to 39.6%, respectively, compared to spatially aggregated regression from all watersheds.For SS, a new regression was developed for data from Maumee, Raisin, and Sandusky.The median relative error of this regression was less than that from the spatially aggregated analysis, at 11.3%, but the maximum relative error and standard deviations declined from 46.0% to 31.5% and from 10.6% to 8.10%, respectively.
Figure 8 also shows the distributions of relative errors for annual load predictions for watershed and pollutant-specific regression equations from the top five load events.These watershed-specific regressions performed better compared to spatially aggregated regressions.For TON, the median relative errors ranged from 7.3% for Muskingum to 18.7% for Vermilion, while the average median relative error for all watersheds was 11.6%.The maximum relative error ranged from 36.9% for Sandusky to 22.5% for Muskingum, while the average maximum relative error was 29.4%.For TP, the median relative errors ranged from 4.9% for Great Miami to 16.2% for Maumee, while the average median relative error for all watersheds was 11.9%.The maximum relative error ranged from 49.0% for Cuyahoga to 11.3% for Great Miami, while the average maximum relative error was 29.8%.For SS, the median relative errors ranged from 6.3% for Great Miami to 16.4% for Maumee, while the average median relative error for all watersheds was 11.3%.The maximum relative error ranged from 44.1% for Cuyahoga to 17.6% for Great Miami, while the average maximum relative error was 25.8%.The average standard deviation observed for these individual watershed specific regressions was 7.19%.

Spatial Transferability
The regression relationship for predicting annual TON loads, developed using spatially aggregated data from all Ohio watersheds (except Cuyahoga, the only watershed with high urbanization), was used to predict annual TON loads for two non-urbanized agricultural watersheds in Illinois, namely Upper Sangamon and Vermilion in Figure 1.This study clearly demonstrates that loads exported during top load events during a year can be used to predict annual loads for watersheds with complete water quality and flow records, but it is challenge to predict, a priori, which large events should be sampled for accurate annual load predictions in the future.As suggested by Markus and Demissie [21], a few monitoring "triggers" can be established for each watershed, based on past flow records.A majority of these top load events also have the highest peak flows in a water year, and as there is better forecasting of flows, triggers could be based on flow thresholds, which would separate the biggest hydrographic peaks in a year.The exact percentile selected can be optimized for a better accuracy, and data mining can be used to determine which auxiliary parameters (such as the antecedent moisture index, time of year, total predicted precipitation, number of dry days, fertilizer application time, etc.) are most critical and should be included in determining the events to be sampled.However, this study focuses solely on demonstrating the spatial transferability of the proposed methodology, and therefore a simple flow threshold value was used.It was found that the 75th percentile threshold separated the largest few events (10-15 events per year) by peak flow for both the Illinois watersheds.TON loads exported through these separated events for each year of record were calculated using only the water samples above the 75th percentile flow threshold, and were then ranked by load export.The top five events each year by load export were then used as inputs for the regression relationship established using the spatially aggregated data from the Ohio watersheds.For comparison, a similar analysis was performed using the complete flows (no thresholds).
The results using the 75th percentile flow threshold indicated that the median absolute prediction error for annual TON loads was approximately 6.42% and 7.52% for the Upper Sangamon and Vermilion watersheds, respectively in Figure 9.The maximum errors were 33.6% for Upper Sangamon and 33.9% for Vermilion.When using the complete flow data from Upper Sangamon and Vermilion, the median absolute prediction errors changed to 7.99% and 11.4%, respectively, and maximum error changed to 24.9% and 13.6%, respectively.

Discussion
The predictability of annual pollutant loads from top load events was evaluated at three different levels of spatial aggregation.Firstly, a spatially aggregated and pollutant-specific regional regression relationship was developed using data from all six Ohio watersheds.Secondly, spatially aggregated regression relationships were developed using data from a few similar watersheds (i.e., all agricultural watersheds).ANCOVA was used to assess the similarity of regression relationships developed for different watersheds individually and for spatially aggregated data from multiple watersheds.The regression relationships from watersheds with similar land-uses were found to be statistically similar to each other, indicating that they can be transferred spatially without the loss of

Discussion
The predictability of annual pollutant loads from top load events was evaluated at three different levels of spatial aggregation.Firstly, a spatially aggregated and pollutant-specific regional regression relationship was developed using data from all six Ohio watersheds.Secondly, spatially aggregated regression relationships were developed using data from a few similar watersheds (i.e., all agricultural watersheds).ANCOVA was used to assess the similarity of regression relationships developed for different watersheds individually and for spatially aggregated data from multiple watersheds.The regression relationships from watersheds with similar land-uses were found to be statistically similar to each other, indicating that they can be transferred spatially without the loss of predictive accuracy.Thirdly, individual watershed-specific regression relationships were used to predict annual loads from top load events.As might be expected, the results indicated that the predictability of annual pollutant loads from top load events increases as the level of spatial aggregation decreases.Regression relationships were weakest for an aggregation of all watersheds by proximity at the regional level, and were strongest at the individual watershed level.Watersheds differ in hydrological characteristics, geology, soils, climate, point sources, land use, and cropping patterns, and thereby have different pollutant loading patterns.Thus, a larger spatial aggregation of watersheds is expected to have a greater variation in pollutant loading patterns, leading to weaker regression relationships between annual loads and top load events.These variations are smaller in similar watersheds, and are smallest at the individual watershed-scale, as shown by the relatively stronger regression relationships with smaller prediction errors.
Among the three pollutants analyzed, annual loads were most predictable for SS, with the prediction errors for all three pollutants declining with a smaller scale of spatial aggregation.Additionally, for all regression relationships evaluated, the variation and mean errors for annual load predictions declined as more top load events were used in calculating the regression relationships, indicating a trade-off in overall monitoring costs and annual load prediction accuracies.However, despite the expected weaker performance of the regional-level spatially aggregated regression relationships, they still had median relative errors below 25% for annual load predictions.These errors are comparable to or lower than those reported in previous studies, which have utilized the commonly used six-to eight-week monitoring frequency for estimating pollutant loads for short monitoring periods [9,33,36,43].
To test the applicability of the developed approach for predicting annual loads for watersheds with limited water quality data, annual TON loads for two watersheds in Illinois were estimated using regression relationships from Ohio.The results provided reasonable estimates, with absolute median errors of less than 12% for both watersheds obtained using the regression relationship from all Ohio watersheds.These results confirmed the merits of the approach, demonstrating that a regional regression relationship may be utilized for predicting annual TON loads for similar watersheds in which the top load events are monitored.
Lastly, there are a few limitations of the proposed approach.In cases with short monitoring records, variations in annual load contributions from top load events during extremely wet or dry years might affect the regression relationships due to extreme values.Similarly, changes in nutrient application patterns and land management practices within a watershed might skew a regression relationship developed during a short monitoring period.Nonetheless, our analysis clearly outlines the role of a few large events in predicting annual pollutant loads.By reducing regular monitoring during low flow periods in a given water year and increasing the focus on big load events, which are generally linked to peak flows, annual pollutant loads may be predicted with a greater efficiency.

Conclusions
The research results provide an additional understanding of riverine pollutant dynamics in relation to hydroclimatic variations and land-use.In particular, this study evaluated the role of large load events in annual pollutant load export by assessing the predictability of annual TON, TP, and SS loads based on top load events in a year.It also provided a comprehensive analysis of long time-series of multiple constituents at eight watersheds with different land-uses, located in the Midwestern United States.A modified pollutant-specific baseflow separation technique was designed to identify high-flow events, and consequently, loads exported during these high-flow events.The results confirm

Figure 1 .
Figure 1.Locations of watersheds and sampling stations.The Cuyahoga, Grand, Maumee, Raisin, Sandusky, and Vermilion watersheds drain into western and central Lake Erie, while Great Miami and Muskingum watersheds drain into the Ohio River basin.The Vermilion, IL and Upper Sangamon watersheds drain west into the Upper Mississippi River basin [7].

Figure 1 .
Figure 1.Locations of watersheds and sampling stations.The Cuyahoga, Grand, Maumee, Raisin, Sandusky, and Vermilion watersheds drain into western and central Lake Erie, while Great Miami and Muskingum watersheds drain into the Ohio River basin.The Vermilion, IL and Upper Sangamon watersheds drain west into the Upper Mississippi River basin [7].

Figure 2 .
Figure 2. Contribution of top high-flow events to annual SS, TP, and TON loads for Maumee and Cuyahoga watersheds.

Figure 2 .
Figure 2. Contribution of top high-flow events to annual SS, TP, and TON loads for Maumee and Cuyahoga watersheds.

Figure 3 .
Figure 3. Percentage of common top events based on loads exported and peak flows for BFR = 1.25 (TON) and BFR = 1.75 (SS, TP).For example, the open circle for the Muskingum BFR = 1.25 case shows that 80% of the top five events, based on TON loading and peak flows, were common over the 15 years of monitoring.Point size indicates relative watershed size (not to scale).

Figure 4 .
Figure 4. Regression relationships for total annual and top event nitrate plus nitrite nitrogen (TON)

Figure 3 .
Figure 3. Percentage of common top events based on loads exported and peak flows for BFR = 1.25 (TON) and BFR = 1.75 (SS, TP).For example, the open circle for the Muskingum BFR = 1.25 case shows that 80% of the top five events, based on TON loading and peak flows, were common over the 15 years of monitoring.Point size indicates relative watershed size (not to scale).

Figure 4 .
Figure 4. Regression relationships for total annual and top event nitrate plus nitrite nitrogen (TON) loads on (a) an individual watershed scale and (b) a spatially aggregated regional scale.Most watersheds in the study are predominantly agricultural, unlike the mostly urbanized Cuyahoga River watershed.

Figure 4 .
Figure 4. Regression relationships for total annual and top event nitrate plus nitrite nitrogen (TON) loads on (a) an individual watershed scale and (b) a spatially aggregated regional scale.Most watersheds in the study are predominantly agricultural, unlike the mostly urbanized Cuyahoga River watershed.

Figure 5 .
Figure 5. Regression relationships for total annual and top event Total Phosphorus (TP) loads on (a) an individual watershed scale and (b) a spatially aggregated regional scale.

Figure 5 .
Figure 5. Regression relationships for total annual and top event Total Phosphorus (TP) loads on (a) an individual watershed scale and (b) a spatially aggregated regional scale.

Figure 6 .
Figure 6.Regression relationships for total annual and top event Suspended Solids (SS) loads on (a) an individual watershed scale and (b) a spatially aggregated regional scale.

Figure 6 .
Figure 6.Regression relationships for total annual and top event Suspended Solids (SS) loads on (a) an individual watershed scale and (b) a spatially aggregated regional scale.

Figure 7 .
Figure 7. Relative error distributions for observed annual loads and annual loads estimated using top one to five nitrate plus nitrite nitrogen (TON) load events for spatially aggregated data from all watersheds.A cross (+) symbol beyond the length of whiskers indicates outliers, and a small square symbol within the boxes indicates the mean of the distributions.

Figure 7 .
Figure 7. Relative error distributions for observed annual loads and annual loads estimated using top one to five nitrate plus nitrite nitrogen (TON) load events for spatially aggregated data from all watersheds.A cross (+) symbol beyond the length of whiskers indicates outliers, and a small square symbol within the boxes indicates the mean of the distributions.

Figure 9 .
Figure 9. Relative error distributions for observed annual loads and annual loads estimated using the top five nitrate plus nitrite nitrogen (TON) load events for watersheds in Illinois using regression equations from Ohio watersheds.Small square symbol within the boxes indicates the mean of the distributions.A cross (+) symbol beyond the length of whiskers indicates outliers, and a small square symbol within the boxes indicates the mean of the distributions.

Figure 9 .
Figure 9. Relative error distributions for observed annual loads and annual loads estimated using the top five nitrate plus nitrite nitrogen (TON) load events for watersheds in Illinois using regression equations from Ohio watersheds.Small square symbol within the boxes indicates the mean of the distributions.A cross (+) symbol beyond the length of whiskers indicates outliers, and a small square symbol within the boxes indicates the mean of the distributions.

Table 1 .
Characteristics of studied watersheds for modeled areas upstream of the United States