Continuous Water-Quality Data Collection

A Yellow Springs Instruments (YSI) EXO turbidity sensor, installed on a YSI EXO multiparameter water-quality sonde, was used to monitor turbidity in formazin nephelometric units at both study sites. The sondes were deployed within a protective polyvinyl chloride pipe with holes cut in it to allow for water flow around the sensors. The sonde at the St. Joseph study site is deployed on the left bank at the downstream abutment of the St. Joseph and Grand Island Railroad Bridge in St. Joseph, Missouri (not shown). The associated gage house that contains the signal transmission equipment is mounted on the bridge railing. The sonde at the Hermann study site is deployed along the bridge pier located on the State Highway 19 bridge at Hermann, Mo. (not shown), and the associated gage house containing the signal transmission equipment is located on the downstream side of the bridge. During some years, sonde deployments were seasonal owing to low water and ice conditions. Data were recorded every 15 minutes and transmitted hourly by way of satellite telemetry. Maintenance and calibration of water-quality sensors followed methods described in Wagner and others (2006). The turbidity sensor was calibrated using a three-point calibration and YSI’s calibration standards. The sonde was serviced approximately every 6 to 8 weeks, which included cleaning of the sonde and calibration checks. Missing periods of record are due to equipment failures, biofouling, and (or) data corrections that exceeded correction criteria (Wagner and others, 2006). All continuous turbidity data are available in the USGS National Water Information System (NWIS) database (U.S. Geological Survey, 2023).

Water Sample Collection and Analysis

Discrete water-quality samples were collected as part of a long-term routine monitoring at both sites funded by the USGS and Missouri Department of Natural Resources. The St. Joseph study site is part of the Missouri Department of Natural Resources and USGS cooperative Ambient Water-Quality Monitoring Network, which is a network of river and stream sites across Missouri that are monitored to assess water resources in Missouri (Markland, 2023). The Hermann study site is part of the USGS National Water Quality Network (NWQN), which is a network of 110 sampling sites established to monitor current water-quality conditions and trends in rivers and streams across the Nation and over time (Riskin and Lee, 2021). Since October 2007, an average of 12 samples per year were collected at the St. Joseph study site and 13 samples per year were collected at the Hermann study site. Sampling was completed according to procedures described by the “National Field Manual for the Collection of Water-Quality Data” (U.S. Geological Survey, variously dated). Samples were analyzed for total phosphorus and other constituents by the USGS National Water Quality Laboratory. All discrete sample data are available in NWIS (U.S. Geological Survey, 2023).

Surrogate Model Development

Surrogate models used simple ordinary least squares regression to relate total phosphorus concentrations in discrete samples to continuous turbidity data. Model development was completed using the statistical software R (R Core Team, 2024) in RStudio software (R version 4.1.3; RStudio Team, 2024). The calibration datasets included discrete total phosphorus concentration matched to the nearest continuous turbidity value. For each study site, a surrogate model for the logarithm (log)-transformed total phosphorous concentration was fitted to the log-transformed turbidity because log-transformation resulted in a better fit of the data. The model form was selected by comparing the Akaike information criterion (AIC) of several potential models (Helsel and others, 2020).

Model residuals were investigated for potential outliers exhibiting undue influence on the models. Plots were evaluated for low residual variance, normality and homoscedasticity of residuals, and a ratio of mean observed (sampled) concentration to computed (simulated) concentration near 1.0. The log-transformed total phosphorus values were retransformed back to the original units to allow for direct calculation of total phosphorus concentrations from continuous turbidity values. Because of the potential bias introduced to the calculated total phosphorus concentration by retransformation, a nonparametric smearing bias-correction factor was applied (Duan, 1983). To compute total phosphorus loads in U.S. short tons per day, the continuous computed total phosphorus concentrations (in milligrams per liter) are multiplied by continuous streamflow (in cubic feet per second) and a unit conversion factor. The daily loads are summarized as the daily mean of the continuous total phosphorus loads and are available in NWIS (U.S. Geological Survey, 2023).

With continued discrete sample collection and continuous turbidity monitor operation, the surrogate models may be used beyond the study period that was used to develop them. Collection of samples throughout the year and during a variety of streamflow conditions along with an annual validation of the regression models to ensure the relation between turbidity and total phosphorus has not significantly changed is warranted. The annual model validation may include plotting new observations with the original model calibration dataset, recomputing the regression coefficients, and determining whether the recomputed coefficients are significantly different than the original coefficients (using an analysis of covariance and probability value significance level of 0.05). The regression models would benefit from continued refinement by the addition of the new data collected each year (Rasmussen and others, 2009).

When a continuous turbidity record was unavailable or when turbidity values were outside of the specified range of the surrogate models, streamflow data were used with the LOADEST model to estimate daily loads of total phosphorus (Runkel and others, 2004). Initially, an alternate surrogate model for each study site was explored for use when a turbidity record was unavailable. Log-transformed streamflow data were fitted to log-transformed total phosphorus concentrations (fig. 1.1). Because of the poor fit between streamflow and total phosphorus concentrations (adjusted coefficient of determination [R²_a] of 0.243 for the St. Joseph study site and 0.453 for the Hermann study site), the alternate surrogate models for the calculation of continuous total phosphorus concentrations were not further pursued. Rather, the USGS rloadest package (Lorenz and others, 2024), which is the R version of the LOADEST program of Runkel and others (2004), was used to select a streamflow-based model to fill in the gaps of the computed daily total phosphorus load record. This process provided daily loads for days when loads from the turbidity-surrogate model were missing. The rloadest package assists the user in developing an appropriate regression model using explanatory variables that include various functions of streamflow and decimal time for the estimation of total phosphorus concentration and load (Runkel and others, 2004). The input variables for the LOADEST models included total phosphorus concentrations, streamflow, time, and seasonality. Coefficients with probability values less than 0.05 were considered significant. The best predefined model generated using LOADEST was selected based on the AIC value and the comparison of the observed and estimated loads across the range of values.

Turbidity-Surrogate Models

Because of curvature in the relation between turbidity and total phosphorus, a quadratic model was chosen for the St. Joseph study site (fig. 3, table 4) based on its low AIC compared to other models. The turbidity-surrogate model for the St. Joseph study site had an R²_a of 0.885, a root mean square error (RMSE) of 0.137, and a retransformation bias correction factor (BCF) of 1.005 (table 4; figs. 3, 4). The model is applicable to turbidity values ≤1,750, the highest turbidity value in the continuous dataset through water year 2022 (tables 2 and 4).

For the Hermann study site, the calibration dataset was split into two datasets and used to develop two surrogate models based on turbidity magnitudes. A low turbidity model was developed with calibration data ≤125 FNUs, and a high turbidity model was developed with calibration data greater than 125 and ≤890 FNUs (the highest turbidity value in the calibration dataset). Curvature between turbidity and total phosphorus also resulted in the selection of a quadratic model for the low and high turbidity datasets (fig. 3, table 4). The low turbidity model for the Hermann study site had an R²_a of 0.690, an RMSE of 0.057, and a BCF of 1.003. The high turbidity model had an R²_a of 0.740, an RMSE of 0.222, and a BCF of 1.004 (table 4; figs. 3, 4).

Quadratic models were chosen for both sites owing to larger residuals on the low and high ends that were not explained by other available parameters and that are likely controlled by different processes. When compared to linear models and models including additional explanatory variables (for example, temperature and streamflow), the quadratic models for both study sites had the lowest AIC. It is hypothesized that the curvature in the turbidity and total phosphorus relation at lower values is due to baseline total phosphorous contributions even during low turbidity and low flow conditions. Both study sites exhibited higher ratios of orthophosphate to total phosphorus during low turbidity conditions, with less of the total phosphorus attributed to orthophosphate during higher turbidity conditions (fig. 5). At the higher values of the turbidity and total phosphorus relation, it is hypothesized that the curvature is due to changes in the chemical or physical composition of the suspended sediment during high flow and high turbidity conditions. Both sites exhibited more consistent increased percentage of fine particles (less than 0.063 millimeter) at higher turbidity values (fig. 6), which are associated with higher total phosphorus concentrations, whereas larger particle sizes (described as percentage of coarse particles) are generally associated with lower total phosphorus concentrations (Capasso and others, 2020).

No outliers were excluded from the turbidity-surrogate model for the St. Joseph study site, whereas two outliers were excluded from the model for the Hermann study site. The outliers occurred on May 29, 2013 (turbidity=570 FNUs and total phosphorus less than 0.01 mg/L), and May 6, 2019 (turbidity=165 FNUs and total phosphorus=0.105 mg/L), and were excluded owing to undue influence on the model. The outliers did not represent extremes in turbidity or streamflow; thus, their exclusion did not constrain the models. These same outliers were also omitted from the input for the LOADEST calculations for the Hermann study site.

Supplemental LOAD ESTimator Model

The turbidity-surrogate models were used to estimate total phosphorus concentrations at the study sites (table 5). A regression model (LOADEST) that uses streamflow as the independent variable can supplement the surrogate model to obtain a complete total phosphorus record. When the turbidity record was unavailable or when values were outside of the specified ranges for the models (413 of the 5,479 days for the St. Joseph study site and 678 of the 5,479 days for the Hermann study site), rloadest regression equations were used to calculate daily total phosphorus loads for a more complete total phosphorus load record (table 5). Seasonal variables were significant and were included in the equation for the St. Joseph study site because total phosphorus followed a seasonal pattern at that site. The seasonal terms did not have the best fit for the Hermann study site and thus were not included in the LOADEST regression model (fig. 7). The rloadest regression equations had R²_a of 0.73 and 0.83 for the St. Joseph and Hermann study sites, respectively. Both sites had overall negative load biases, with a larger bias at the St. Joseph study site (−12.43 percent) than the Hermann study site (−2.195 percent) (table 5).

Computed Total Phosphorus Concentrations, Loads, and Yields

During the 15-year study period for the St. Joseph study site, daily mean total phosphorus concentrations calculated using the surrogate model ranged from 0.104 to 4.56 mg/L (median of 0.272 mg/L and mean of 0.417 mg/L; table 6). Computed total phosphorus daily loads from the turbidity-surrogate model (with gaps in record filled with LOADEST data) ranged from 5.19 to 1,760 tons/d (median of 36.5 tons/d and mean of 79.8 tons/d), and daily yields ranged from 0.0243 to 8.24 pounds per square mile per day (lb/mi²/d; median of 0.171 lb/mi²/d and mean of 0.374 lb/mi²/d). Depending on the year, 67.1 to 100 percent of the daily loads were computed using the surrogate model, with 92.4 percent of the total load for the 15-year period calculated using the surrogate model (table 6).

From the surrogate model with gaps filled with LOADEST models, annual total phosphorus loads at the St. Joseph study site ranged from 9,570 tons in water year 2022 to 50,500 tons in water year 2019 (total load for the study period of 437,000 tons), and annual yields ranged from 44.9 pounds per square mile per year (lb/mi²/yr) in water year 2022 to 237 lb/mi²/yr in water year 2019 (total yield for entire study period of 137 lb/mi²/yr; table 6, fig. 8).

For comparison purposes, total phosphorus loads computed only using LOADEST models were compared to the annual load calculated using the surrogate model for water year 2018 because it had a complete daily load record and did not require gaps to be filled with LOADEST data. The annual load calculated using the surrogate model for water year 2018 was 24,600 tons, and the annual load calculated using only LOADEST data was 27,600 tons (11.5-percent difference; Markland, 2024). Despite different inputs for these two models, the annual loads are within less than 15 percent of one another, which lends credence to the models. Total phosphorus concentrations, loads, and yields estimated using turbidity as the independent variable in the surrogate models are preferentially used instead of those using flow in the LOADEST models because total phosphorus concentrations are more highly correlated with turbidity than with flow. Estimates from the LOADEST models were used when turbidity data were not available.

For the Hermann study site, daily mean total phosphorus concentrations calculated using only surrogate regression models ranged from 0.183 to 1.97 mg/L (median of 0.319 mg/L and mean of 0.421 mg/L; table 6). Computed total phosphorus daily loads from the turbidity-surrogate model (with gaps in record filled with LOADEST data) ranged from 12.7 to 1,970 tons/d (median of 76.8 tons/d and mean of 152 tons/d), and daily yields ranged from 0.0486 to 7.55 lb/mi²/d (median of 0.294 lb/mi²/d and mean of 0.582 lb/mi²/d). Depending on the year, 54.1 to 100 percent of the daily loads were computed using the surrogate model, with 87.6 percent of the total load for the 15-year period calculated using the surrogate model (table 6).

Annual loads at the Hermann study site were approximately twice as much as those computed for the St. Joseph study site, ranging from 22,600 tons in water year 2022 to 101,000 tons in water year 2019 (total load for the study period of 833,000 tons; table 6, fig. 8). Annual yields ranged from 86.7 lb/mi²/yr in water year 2022 to 385 lb/mi²/yr in water year 2019 (total yield for entire study period of 213 lb/mi²/yr; table 6). For comparison purposes, total phosphorus loads computed using LOADEST were compared to the annual load calculated using the surrogate model for water year 2018 because it had a complete daily load record and did not require gaps to be filled with LOADEST data. The annual load calculated using the surrogate model for water year 2018 was 34,000 tons, and the annual load calculated using LOADEST was 30,700 tons (10.2-percent difference; Markland, 2024).

Because the Hermann study site is a USGS NWQN monitoring location, additional concentration and load estimates, as well as long-term trend information, are available using other methods for comparison to the concentrations and loads calculated in this report. Weighted regressions on time, discharge, and season (WRTDS) software has been used to calculate flow-normalized (removes effects of flow-related variability to aid in long-term trend analysis) total phosphorus loads at the Hermann study site since 1970, with results available through NWQN (U.S. Geological Survey, 2022). Discharge, as used in WRTDS, is the same as streamflow. For annual flow-normalized total phosphorus loads, a 3-percent decrease (attributed to changes in the watershed that were not streamflow-related) was observed for the period from 1982 to 2022. As for more recent trend periods, a 38-percent decrease and a 24-percent decrease were observed from 2002 to 2022 and from 2012 to 2022, respectively. As for annual loads of total phosphorus, approximately 716,000 tons of total phosphorus were calculated using the WRTDS model in water years 2008 through 2022 (U.S. Geological Survey, 2022), compared to the 833,000 tons calculated in this report (a 15.1-percent difference). The differences between various methods used to estimate total phosphorus loads have been explored, with total phosphorus being one of the more difficult constituents for which to accurately estimate annual loads (Lee and others, 2016, 2019). Because most of total phosphorus can often be transported during a small number of larger hydrologic events, it can be difficult to accurately estimate loads with discrete sampling methods when those events are not sampled (Lee and others, 2019). Conversely, the use of continuous turbidity monitor data and surrogate models can often better capture the range of hydrologic conditions and these larger events, which can result in the differences between methods in annual load estimations.

The sources of nutrient loads and yields from major subbasins within the Missouri River drainage basin were investigated using SPAtially Referenced Regressions On Watershed attributes (SPARROW models by Brown and others (2011). The total phosphorus load delivered to the Mississippi River from the Missouri River was estimated at approximately 43,000 tons per year (loads detrended to the base year 2002; Brown and others, 2011). The largest proportion of total phosphorus was contributed by the Middle Missouri and Lower Missouri subbasins (18.6 and 38.7 percent, respectively). The smallest proportions were delivered by the Upper Missouri and Yellowstone subbasins (0.1 percent each; Brown and others, 2011).

Although these studies of total phosphorus loading in the Missouri River provide important information, the use of turbidity time-series data in surrogate models improves the ability to estimate total phosphorus loads (Schilling and others, 2017). The use of continuous data allows for the capture of variability in nutrient transport during the full range of hydrologic conditions, which is not easily done by discrete sampling (Stutter and others, 2017; Villa and others, 2019). Using turbidity as a surrogate allows monitoring of total phosphorus, which is not practical to measure in a field setting, in a way that captures the concentrations and loads of total phosphorus in the Missouri River. The ability to calculate total phosphorus loads is critical in assessing progress towards nutrient reduction goals.