City of Durham

Water-quality data for the DSS monitoring stations on Ellerbe Creek and Sandy Creek were obtained through the City of Durham’s Water Quality Data Web Portal (City of Durham, 2024). Only results for samples collected by DSS during January 2009–December 2020 as part of their ambient monitoring program and analyzed by the South Durham Water Reclamation Facility (SDWRF) laboratory were considered for use in this study (City of Durham Public Works Department, 2019). The SDWRF is listed as laboratory number 176 under the North Carolina Department of Environmental Quality-certified municipal and industrial laboratories (North Carolina Department of Environmental Quality, 2023) certified to run analytical methods to meet NPDES permit requirements (U.S. Environmental Protection Agency [EPA], 2021). DSS manually collects samples on a fixed monthly basis and occasionally rotates monitoring sites. Ellerbe Creek was not sampled during 2014; therefore, no water-quality data were available for that year. Water-quality data were available for Sandy Creek for the entire period.

DSS water-quality data included remark codes (RCs) associated with analytical results, such as “less than” (<) for values censored at the laboratory reporting level (LRL) and “estimated” (E) for values that fall between the LRL and method detection level. Constituents with estimated concentrations were used as reported for the load analyses in this study. Laboratory analyses were performed on unfiltered samples and included concentrations for unfiltered phosphorus by EPA method 365.1, unfiltered ammonia + organic N by EPA method 351.2, and unfiltered nitrate + nitrite by EPA method 353.2 (table 2; EPA, 1993a, b, c, respectively). Hereafter, filtered constituents will be referred to as “dissolved,” and unfiltered constituents will be referred to as “total.” Total ammonia + organic N also will be referred to as “TKN.”

USGS

Water-quality data for the USGS monitoring stations on Ellerbe Creek and Sandy Creek included historical data and new data collected during this study. At Sandy Creek, nutrient concentration data for dissolved nitrate + nitrite, total N, and total P were available for 10 surface-water samples collected weekly from April 8 to June 12, 2014, as part of a regional USGS study (Van Metre and Journey, 2014; Journey and others, 2015). Results for only 3 of these 10 USGS samples collected at higher streamflows (April 8 and 30, 2014; and May 28, 2014) than the DSS monthly samples during that time period were included in the calibration files for the load estimation analyses. Water-quality datasets also included the nutrient results for surface-water samples collected during 2018–20 by the USGS for this study to supplement the fixed-interval data collected through the DSS ambient monitoring program. The EPA-updated approved laboratory methods for NPDES monitoring include the USGS and DSS methods for nutrients used in this study (table 2), indicating a degree of compatibility in the results from these methods (EPA, 2021).

The objective of the USGS sampling effort was to collect nutrient data across a wider range of streamflow conditions than previously attempted and at times of the diel cycle that previously had been underrepresented by the DSS’s monitoring program. At the beginning of the study, the USGS installed ISCO automated samplers to collect surface-water samples at the Ellerbe Creek and Sandy Creek study sites. The USGS has existing continuous streamgages at the selected sites that provide streamflow data from August 1, 2008, to the present day. Streamflow at these small, urban stream sites was described as being highly flashy, having a rapid, short-lived response to stormwater runoff events. Peak streamflows from heavy rainfall, especially late-day convective thunderstorms, commonly occurred at times outside of a normal work week, such as during nights and weekends. The automated samplers were equipped with two-way cellular telemetry, which allowed them to be queried and controlled remotely and provided flexibility in water-quality sample collection during periods of high streamflows, regardless of the time of day or day of the week.

A total of 24 high-flow samples were collected at each site by the USGS, mainly during stormwater runoff events by using automated samplers, from October 2016 to January 2018 for laboratory analysis of nutrients (table 2). During stormflow, the USGS high-flow sample collections targeted different parts of the stormflow hydrograph to evaluate water chemistry changes during rising limb, peak, and falling limb. The samplers were programmed following manufacturer’s instructions to automatically collect an individual sample once a predetermined stream stage condition was met during a stormwater runoff event. The stream sample was collected for about 1 minute and deposited into a 1-liter collection bottle stored within the automated sampler. Samples were retrieved as soon as possible following collection for subsequent processing at the USGS office laboratory in Raleigh, North Carolina, typically within 12–16 hours of collection, which met the 28-hour hold time for unprocessed nutrient samples (EPA, 2009). The samples were collected and processed using techniques described in the USGS “National Field Manual for the Collection of Water-Quality Data” (U.S. Geological Survey, variously dated). Specific conductance and turbidity were measured in the office laboratory. Processed samples were shipped on ice overnight within 1–2 days of collection and were received by USGS National Water Quality Laboratory (NWQL) in Denver, Colorado, for TKN, dissolved nitrate + nitrite, and total P analyses. The NWQL analytical methods for nutrients remained consistent for all USGS samples analyzed during the study (table 2). All USGS data compiled in the water-quality datasets were retrieved from the USGS National Water Information System (NWIS) database (U.S. Geological Survey, 2021a).

Quality Assurance and Quality Control

The City of Durham data were quality-assured by DSS prior to being downloaded from their Water Quality Data Web Portal. The USGS also worked in consultation with DSS staff to further review the City of Durham nutrient data to identify any potential issues or questionable results prior to finalizing the dataset for use in computing nutrient loads. Results for quality assurance samples were retrieved and used to evaluate the environmental sample results. During the review, it was noted that the City of Durham nutrient data from the SDWRF had fluctuating LRLs, producing multiple censoring levels (table 2). Beginning around October 2012, the SDWRF started reporting TKN and nitrate + nitrite concentration values with a higher precision than before, to two decimal places rather than one. The LRLs for the City of Durham nutrient data ranged from 0.12 to 0.9 mg/L for TKN, from 0.008 to 0.1 mg/L for nitrate + nitrite, and from 0.02 to 0.1 mg/L for total P concentrations. Sixteen of 17 samples collected during a limited span from January 2009 to May 2010 at Sandy Creek had total P concentrations reported below the LRL of 0.03 mg/L. The total P detection rate of 5.9 percent during that period was considerably lower than the 88-percent detection rate observed over the entire study period at this site. Conversely, only 3 of 17 samples collected during the same period at Ellerbe Creek had total P concentrations <0.03 mg/L. Because there was no basis to determine whether the lack of detectable total P concentrations at Sandy Creek during January 2009–May 2010 was reflective of true site conditions or some type of laboratory issue that may have occurred during sample analysis, these data were retained in the final dataset used for the January 2009–December 2020 period.

The USGS collected quality-control samples, including field blanks and replicate samples, in association with the surface-water samples to document potential bias and variability in data that may result during the collection, processing, shipping, and handling of environmental samples (U.S. Geological Survey, variously dated). Field blanks were collected and processed using inorganic-free blank water with the same equipment used for the environmental samples, including the automated sampler, bottles, filters, and preservatives. Field blanks help identify potential contamination resulting from improperly cleaned equipment, field sampling activities and exposure, and laboratory practices. Replicate samples help document the variability in data results associated with sample collection, processing, and laboratory analysis. Analytical results of the quality-assurance samples collected by USGS during this study are summarized in table 3.

Four field equipment blanks were collected and analyzed for specific conductance, TKN, dissolved nitrate + nitrite, and total P (table 3). None of the nutrient constituents were detected above LRLs in the field blanks; thus, there was no indication of systematic positive bias that would affect the environmental data.

Eight paired environmental and replicate sample sets were collected during the study (table 3). The two environmental replicate sample sets collected on June 26, 2017, were used to test for potential temperature effects on surface-water samples following collection and storage in the automated samplers. The “temperature hold test” helped to evaluate whether nutrient constituents in samples experienced any degradation effects during the time they were held in the automated samplers, under ambient conditions, and subsequently retrieved for processing and submittal to the NWQL. For this test, streamwater was collected through the automated sampler into a churn. The composited water was split into two sample collection bottles, one for the regular environmental sample and one for the replicate sample. The environmental sample was processed soon after collection, whereas the replicate sample was placed back in the automated sampler for overnight storage; air temperatures were in the 70- to 80-degree Fahrenheit range during the test. The replicate sample was retrieved for processing the next day, about 18 hours following collection, which was above the average retrieval and processing time for the environmental samples. Overall, the temperature hold test results indicated very good agreement, with minimal to no differences between the nutrient constituent concentrations for the paired replicate samples (table 3). These results indicate nutrient concentrations in collected surface-water samples were not negatively affected when held for periods of up to 18 hours before retrieval from the automated samplers and subsequent processing. Each of the other six paired environmental-replicate sample sets were split from a single collection bottle and processed as soon as possible following collection and retrieval from the automated samplers.

The overall variability of sampling and analysis was assessed by examining the relative percent differences (RPDs)—calculated as the absolute difference multiplied by 100 and divided by the average—between the constituent concentrations in all eight replicate sample sets. A statistical summary of the RPDs determined for each analyte for all paired environmental and replicate samples is provided in table 4. The RPDs in constituent concentrations exceeded 10 percent in only one instance. The mean and median RPDs were less than 5 percent between the environmental and replicate samples, which indicated very good reproducibility for the measured water-quality constituents. No quality-assurance problems were identified for the USGS environmental dataset on the basis of the replicate samples.

Water-Quality Datasets for Load Estimation

The following discussion describes the differences in the nutrient concentration data obtained from both data sources (City of Durham and USGS) and the steps taken to prepare the final water-quality datasets for subsequent use in developing the nutrient data calibration files for the R-LOADEST program. Both agencies reported concentrations of the various nitrogen fractions in milligrams per liter as N and concentrations of total P in milligrams per liter as P.

The nutrient data retrieved from the City of Durham included concentrations for TKN, total nitrate + nitrite, and total P. Because of nitrate’s and nitrite’s high solubility in water, concentrations reported by City of Durham for total nitrate + nitrite were assumed to be equivalent to concentrations of dissolved nitrate + nitrite reported by USGS. Values of total N for the City of Durham data were computed by summing the reported concentrations for TKN and total nitrate + nitrite. The USGS data for this study included laboratory results for TKN, dissolved nitrate + nitrite, and total P; results for 2014 USGS regional study data included laboratory results of total N but not TKN (Journey and others, 2015). Values of total N for USGS samples collected during this study were computed using the measured concentrations for TKN and dissolved nitrate + nitrite. The approach used in this study for calculating total N concentrations when a constituent concentration was below the LRL is discussed in Harden and others (2020).

The handling of estimated concentrations reported for TKN and (or) nitrate + nitrite when computing values of total N concentrations is described as follows. If both TKN and nitrate + nitrite had estimated concentrations, then an “E” RC also was assigned to the total N concentration. If only one of the underlying constituents was estimated and that concentration represented 25 percent or more of the computed total N value, then an “E” RC was assigned to the total N concentration; an “E” was not assigned if the concentration was less than 25 percent.

For each study site, the sampling frequency of compiled data was reviewed. On five occasions, analytical results for samples collected by both the City of Durham and USGS were available on the same date. In these cases, only the sample associated with the highest instantaneous streamflow at the time of collection was retained for use in this study. The final water-quality datasets compiled for the Ellerbe Creek and Sandy Creek study sites are presented in Harden and others (2024).

R-LOADEST Models for Estimating Constituent Loads

The R-LOADEST program (Runkel and others, 2004; Runkel, 2013; Lorenz and others, 2017) was used for estimating annual stream loads of TKN, nitrate + nitrite, total N, and total P at Ellerbe Creek and Sandy Creek. The selected R version of the LOADEST program (Lorenz and others, 2017) used for this study also was consistent with the Microsoft Windows based version of the LOADEST software package (Runkel and others, 2004; Runkel, 2013) used by DSS for computing nutrient loads as part of their stream ambient water-quality monitoring program. Instantaneous streamflow values, in cubic feet per second, were used with the sample nutrient concentrations in the calibration input files. The nitrate concentrations in the calibration files actually represented the reported concentrations of nitrate + nitrite, but because nitrite typically represented a small fraction of the total concentration, the concentrations and subsequent loads computed with the R-LOADEST program will be presented and discussed as nitrate. Daily mean streamflow values, in cubic feet per second, were used in the streamflow prediction files that were combined with the calibrated models in R-LOADEST for estimating nutrient loads.

R-LOADEST runs nine predefined regression models that use up to seven explanatory variables to estimate stream loads of each constituent (table 5; Runkel and others, 2004). Within R-LOADEST, streamflow and constituent concentrations were transformed to natural logarithms prior to regression analysis to improve the fit of the regression model. The simplest predefined models include only streamflow variables (models 1 and 2). These simple models often work well with estimation of particulate-dominant constituents like suspended sediment (Crawford, 1991). Models 3 and 5 account for streamflow variability and long-term temporal changes in streamflow or concentrations (Helsel and others, 2020). The models that account for streamflow and seasonality only (4 and 6) or streamflow, time trends, and seasonality (7, 8, and 9) often perform well for nutrient load estimation (Helsel and others, 2020). The R-LOADEST package allows the user to manually select any of the nine regression models on the basis of multiple diagnostic measures of bias or to use an automated option that selects the best model on the basis of one diagnostic measure.

Initially, the automated option was used to select the best regression model for each constituent at each site (Runkel, 2013). The calibration and estimation procedures used with the selected regression models within R-LOADEST were based on the adjusted maximum likelihood estimation (AMLE) statistical estimation method (Cohn and others, 1989; Cohn, 2005). The AMLE method is appropriate when the model calibration errors are normally distributed and is the method of choice when the calibration datasets contain censored data; most calibration data compiled for the study sites typically contain censored data at multiple levels (between 8.5 and 15.7 percent of values censored). The ESTIMATOR routine within LOADEST assigns a “censoring threshold” for a constituent, even where there are multiple censoring levels prior to the AMLE routine. Runkel (2013) reported that the AMLE method is relatively insensitive to the precise value of the censoring threshold, even when the associated observation is uncensored. When LOADEST is executed, the first censored observation is used to assign an initial value to compute the censoring threshold until another censored observation is encountered.

Model evaluation criteria presented with the automated option include statistical values for both the Akaike information criterion (AIC) and the Schwarz posterior probability criterion (SPPC) that are computed and ranked for each of the calibrated models (Judge and others, 1988; Runkel and others, 2004). The best-fit model is selected on the basis of its lowest ranking of the AIC statistic, indicating relative goodness of fit between observations and predictions. The SPPC, which is also known as the Bayesian information criterion, is not used directly as part of the automated model selection process; however, it is included with the model output for comparative purposes. The SPPC has a higher penalty than does AIC for additional parameters in the model and will generally favor a simpler model with fewer terms (Runkel, 2013).

Subsequent to using the automated option to select the best-fit model, additional predefined models that are standard in the LOADEST package were considered to develop 1–3 candidate models on the basis of a review of their rankings for the AIC and SPPC statistics. As candidate models were developed for each scenario and time step, other diagnostics, including regression statistics, load bias statistics, and residuals analysis plots published in the R-LOADEST update (Runkel, 2013; Lorenz and others, 2017), were utilized to select the final model for that scenario and time step to be included in the analysis of nutrient loads at the study sites.

The probability plot correlation coefficient (PPCC) values included with the regression statistics and the residual analysis plots provided a check on the assumption that the model residuals were normally distributed. A PPCC value of 1.00 represents a perfect linear relationship. If the reported p-value for the PPCC test was <0.05, then the assumption that the data were normally distributed can be rejected; thus, p-values greater than (>) 0.05 were preferred.

The load bias statistics, including the load bias in percent (Bp) and the Nash-Sutcliffe Efficiency Index (NSE), were based on a comparison of observed and estimated loads for the dates and times in the model calibration data (Runkel, 2013). Values of Bp are equal to the sum of estimated loads minus the sum of observed loads, divided by the sum of observed loads, and expressed as a percentage (Runkel, 2013); positive values indicate overestimation of loads, and negative values indicate underestimation of loads. Values of NSE can theoretically range from negative infinity to 1.0, with a value of 1.0 indicating a perfect fit between observed and estimated loads, and a value of 0.0 indicating that the load estimates are only as accurate as the mean (Runkel, 2013). Runkel (2013) indicates that calibrated models having an absolute value of Bp exceeding 25 percent or a negative NSE value should not be used for load estimation.

The R-LOADEST program output generated for each constituent and year were compiled in Harden and others (2024). Results included the regression model used, the number of days in the year, the annual mean daily load (in kilograms per day), the standard error of prediction (SEP) of the load, and the lower and upper 95-percent confidence intervals (L95s and U95s) of the load (Runkel and others, 2004). Although not discussed in this report, the cumulative constituent load for a given year, in kilograms per year, can be determined by multiplying the mean daily load by the number of days in the year. The SEP represented the uncertainty associated with each estimate of daily load (Gilroy and others, 1990; Cohn and others, 1992) and incorporated both variability attributed to the model calibration (parameter uncertainty) and unexplained variability about the model (random error). The SEP indicated how closely estimated loads correspond to actual loads and was used within R-LOADEST to develop the 95-percent confidence intervals for each daily load estimate.

Load Testing Scenarios

Lee and others (2016) evaluated the relative performance of different methods for estimating constituent loads across a broad range of water-quality constituents, stream types, and sampling regimes. The tested sampling strategies included the “scenario B” monitoring approach in this study; that is, the collection of fixed monthly samples. One of the main findings in Lee and others (2016) and Robertson and Roerish (1999) was that sampling strategies that incorporate high-flow samples produced more accurate estimates of constituent loads with lower error levels than those that did not include high flow samples, regardless of the load estimation method.

For this study, subsets of the water-quality data compiled for Ellerbe and Sandy Creeks were used to test different load-estimation scenarios with R-LOADEST. The test scenarios included different combinations of the DSS and USGS nutrient data for the study period of January 2009–December 2020 for calibrating models and predicting annual nutrient loads for five select time periods during 2009–20 (table 6). Tests 1, 2, 3, 4, and 5 represented decreasing time steps of 12, 9, 7, 6, and 5 years, respectively. Additionally, two scenarios for each time-step test were run, where “scenario A” represented load estimations based on both DSS (fixed-monthly sampling) and USGS (targeted high-flow sampling) calibration data, and “scenario B” represented load estimations based on the DSS data only (table 6). The tests were designed to examine how uncertainties associated with the load predictions may be influenced by the length of record analyzed and the inclusion or exclusion of high-flow samples. The individual nutrient calibration files and streamflow prediction files used for the different test scenarios are provided in Harden and others (2024).

Influence of High-Streamflow Samples

Comparison of the test scenario results illustrated the influence of high-flow samples on the load estimates and uncertainties. For all tests having paired scenario A and scenario B results, A tests produced lower mean annual SEPs relative to their paired B tests at Ellerbe and Sandy Creeks (table 10). In other words, the inclusion of high-flow samples with fixed frequency samples decreased the load uncertainties for all nutrient fractions at both Ellerbe Creek (fig. 11) and Sandy Creek (fig. 12). Overall, SEPs increased with higher nutrient loads among the time steps for scenario A at Ellerbe Creek and but were slightly more uniform at Sandy Creek. Additionally, ranges in SEPs for scenario A at Ellerbe Creek and Sandy Creek were more consistent among time steps compared to the SEPs for scenario B. Even when the mean annual load for a constituent at a given time step in scenario A was similar to or higher than that for scenario B (with exception of estimated nitrate loads at Sandy Creek), the corresponding SEPs for scenario B were higher, with some constituents up to 2.0 to 2.8 times higher than scenario A at Ellerbe Creek (total N; table 10 and fig. 11C) and up to 2.5 to 3.2 times higher than scenario A at Sandy Creek (nitrate; table 10 and fig. 12B).

At Ellerbe Creek and Sandy Creek, the mean annual constituent loads were not significantly different among time steps between scenarios A and B, except for nitrate at Sandy Creek, which was statistically lower for the A tests relative to the B tests (table 11). The corresponding SEPs for TKN, nitrate, and total N at both sites were all statistically lower for scenario A compared to scenario B (figs. 11 and 12).

A comparison of scenario A results between Ellerbe and Sandy Creeks indicated that mean annual loads and SEPs were not significantly different for TKN, nitrate, or total N (table 11). However, comparison of scenario B results between the two study sites indicated that the mean annual load for nitrate and both the mean annual load and SEP for total P were statistically higher for Sandy Creek compared to Ellerbe Creek.

The influence of the USGS high-flow samples on the model calibration datasets was reflected by the mean instantaneous streamflows for samples included for each test scenario (table 12). The mean sample streamflow for Ellerbe Creek was about 5.6 times higher for test 1A than for test 1B and increased to about 16.6 times higher for test 5A relative to test 5B. The mean sample streamflow for Sandy Creek was about 4.4 times higher for test 1A than for test 1B and increased to about 10.9 times higher for test 5A relative to test 5B.

As previously indicated, 77.2 percent and 71.4 percent of the study period streamflows for Ellerbe Creek and Sandy Creek, respectively, were contributed from the uppermost 10 percent (90–100th percentiles) of daily mean streamflows (table 7). The percentage of days sampled by both DSS and USGS during the highest 10 percent of flows was 17.6 percent for Ellerbe Creek and 19.2 percent for Sandy Creek. These percentages were reduced by more than half when the USGS high-flow samples were not included.

Exclusion of the USGS high-flow samples from the model calibration datasets typically resulted in significantly higher uncertainties in the load estimates. Although individual nutrient loads were relatively similar between calibration files (with or without higher flow samples), the calibration data with higher flow samples were more representative of the observed flows in prediction files, which resulted in lower uncertainties in predicted loads. The calibrated regression models based primarily on lower flow sample data had a higher degree of uncertainty in the estimated loads when extrapolated to the higher flow days in the streamflow prediction files.

Because most instream nutrient transport occurs during periods of higher streamflows, it is important that water-quality samples collected for computing stream nutrient loads are representative of all streamflow conditions at a site, especially for the upper flow regime (Lee and others, 2016; Aulenbach and others, 2017; Lee and others, 2019). The load testing scenario evaluations in this study are in agreement with Lee and others (2016) in that sampling programs that incorporated high-flow samples yielded better estimates of constituent loads having lower uncertainties than sampling programs that did not include high-flow samples.

Influence of Load Prediction Time Step

The test scenario results also were used to examine the influence of prediction time step on the load uncertainties by evaluating changes in mean annual loads and SEPs or 95-percent confidence intervals (developed from SEPs; Runkel, 2013) across time steps. The five prediction time steps reflected different subsets of the overall dataset and included test 1 (all 12 years), test 2 (9 years), test 3 (7 years), test 4 (6 years), and test 5 (5 years).

Thus far, the test evaluations have been based on mean annual loads and SEPs, which represented all the years in the prediction period for an individual test scenario or time step (table 6). Each test scenario represented a different snapshot of the study where, for example, test 1 covered the entire study period of 2009–20 and test 5 covered the last 5 years of the study, 2016–20. The mean annual load and SEP results between the five tests were not expected to be the same because scenarios having smaller time steps contained less record with fewer calibration data points than those with longer time steps and did not adequately represent the annual variability in streamflows throughout the study period. Therefore, smaller time steps can produce increased uncertainty in the load estimation compared to larger, more data-rich time steps.

The following discussion focuses on the year-to-year results for each test scenario (based on different time steps) that were computed for the same 5-year period (2016–20) when all five test scenarios were available (table 6). To better understand the influence of time step used for model calibration and prediction, load estimates and uncertainties for the different time steps were compared for the same year(s). The annual mean loads and their lower and upper 95-percent confidence intervals (L95s and U95s) for each nutrient from 2016 through 2020 are shown by test scenario and time step in figure 13A–D for Ellerbe Creek and figure 14A–D for Sandy Creek. There were no test 3 constituent results for Ellerbe Creek and no scenario A results for total P at Sandy Creek.

Examination of the test results for each site indicates that the confidence intervals of the annual mean nutrient loads overlapped across scenarios A and B for a given year. There was notably higher variability among the loads and wider confidence intervals within a given year for the B tests relative to the A tests at both sites (figs. 13 and 14). For example, at Ellerbe Creek, TKN test 2A for 2018 (fig. 13A) had a mean load of 40.84 kg/d, with an L95 of 31.18 kg/d and U95 of 52.55 kg/d (spread = 21.37). The TKN result for test 2B for 2018 had a mean load of 40.98 kg/d, with an L95 of 27.20 kg/d and U95 of 59.34 kg/d (spread = 32.14). In other cases, large differences were noted for both the mean daily load and confidence intervals between the paired test scenarios. At Sandy Creek for example, nitrate test 3A for 2020 (fig. 14B) had a mean load of 13.02 kg/d, with an L95 of 9.90 kg/d and U95 of 16.83 kg/d (spread = 6.93). The nitrate result for test 3B for 2020 had a mean load of 18.90 kg/d with an L95 of 11.01 kg/d and U95 of 30.35 kg/d (spread = 19.34).

Annual nutrient loads also reflected the variability in mean annual daily mean streamflow (fig. 3, table 13). For the 5-year period of 2016–20, annual daily mean streamflows at Ellerbe and Sandy Creeks were lowest in 2016 and 2017, intermediate in 2019, and highest in 2018 and 2020. Correspondingly, loads of all nutrient fractions were lowest during 2016 and 2017, intermediate in 2019, and highest during 2018 and 2020. Similarly, the confidence intervals were smaller during lower flow years and larger during higher flow years.

Given the higher variability and uncertainties for the scenario B results relative to the scenario A results, only the scenario A results were used to explore potential relations among year-to-year loads for individual nutrients and the length of record, or time step, analyzed. The TKN loads and confidence intervals for both Ellerbe Creek (fig. 13A) and Sandy Creek (fig. 14A) are mostly consistent among the time steps for each year with the exception of 2020 for Sandy Creek.

At Ellerbe Creek, nitrate load and confidence intervals generally increased with decreasing time steps during 2016 and 2020 and decreased with decreasing time steps during 2018 and 2019 (fig. 13B). Although there are no nitrate load results for test 1A at Sandy Creek (fig. 14B), there is a consistent pattern of greater nitrate loads and confidence intervals at the 9-year time step (test 2A) relative to the 7-, 6-, and 5-year time steps (tests 3A to 5A), which could be influenced by the changed LRLs in 2012 (fig. 7D).

The year-to-year total N loads for Ellerbe Creek (fig. 13C) had a similar distribution as described above for nitrate, where for each year the total N loads and confidence intervals tended to be more similar between the 12- and 9-year time steps relative to the 6- and 5-year time steps. For Sandy Creek, the yearly total N loads (fig. 14C) generally followed the consistent pattern noted for TKN (fig. 14A), being relatively uniform among the time steps for years 2016–19 and more variable in 2020. Scenario A data for total P are only available for Ellerbe Creek (fig. 13D), which shows that the annual loads and confidence intervals are mostly consistent across time steps for each year. The total P load confidence intervals are substantially largest during 2018, which had the highest mean annual streamflow (table 13).

The evaluations of the year-to-year test scenario results during 2016–20 indicated how the nutrient load estimates and confidence intervals at Ellerbe and Sandy Creeks varied depending on the time step used for the load model calibrations and predictions. In general, prediction errors for most nutrient fractions tended to be higher for shorter time steps, especially when targeted high-flow samples were not included. Although not a definitive pattern, the longer 12-year and 9-year time steps tended to provide more consistent annual load estimates and lower confidence intervals than the shorter 6- and 5-year time steps. This finding emphasizes the importance of analyzing a study period that adequately represents longer-term variations in streamflow to better understand the temporal variability among annual nutrient loads and uncertainties at the study sites.