Nitrate and biochemical oxygen demand change in a typical Midwest stream in the past two decades

Nitrate and organic contamination fromMidwest rivers, including theWhite River at Muncie, IN, has been an on-going concern and contributes to the hypoxic zone in the Gulf. Despite rich data, recent water quality changes have rarely been investigated. This study employed 16 years of continuous monitoring data, including biochemical oxygen demand (BOD), dissolved oxygen (DO), and nitrate–nitrite as nitrogen (NN) from five sites near Muncie, and analyzed thewater quality trend and pollution sources. A novel approach,Weighted Regression onTime, Discharge and Seasons (WRTDS) that allows for the representation of long-term water quality patterns by considering seasonal variance and discharge-related effects over time, is adopted. Flow-normalized BOD and NN concentration and flux both increased, and DO concentration and flux decreased. However, the changes vary among sites. Muncie wastewater treatment plant and combined sewage outflows (CSOs) contribute remarkably to NN pollution during low-flow seasons. Urban and agricultural runoff, and CSOs impact BOD levels. Agricultural runoff contribution to BOD is increasing in recent years. Seasonal patterns of nitrate and BOD in the river are also analyzed. The results are helpful for watershed managers to re-think conservation practices and have indications to water quality management beyond the study area.


INTRODUCTION
Water quality degradation has been a significant concern in the US Midwest. As a consequence of industrial, domestic, and urban/agricultural runoff, the White River was heavily contaminated prior to the 1970s (Tiffany 2018). Despite improvements, much of the White River is still impaired by organic compounds and nutrients (USEPA 2017). Excess nutrients and the induced eutrophication are among the major on-going concerns for the White River (Lowe et al. 2008) and farther downstream (Bargu et al. 2019), as the White River flows into the Mississippi River and drains into the Gulf of Mexico. As such, downstream aquatic environment degradation including the so-called 'dead zone' in the Gulf of Mexico where fish and marine life is affected by the low level of dissolved oxygen (DO) level could potentially be traced back to the nutrient and biochemical oxygen demand (BOD) load in the upstream river basins including the White River watershed (Bargu et al. 2019).
An appropriate DO level in water is critical for maintaining good water quality and a healthy aquatic environment (Banerjee et al. 2019). Oxygen in water bodies could be produced by aquatic plant photosynthesis, absorbed from the atmosphere by dissolving oxygen through the air-water interface, or mixed in from turbulence. However, DO can be consumed or even depleted due to plant and animal respiration and from the decomposition of organic matter as measured by BOD (Wen et al. 2017;Suplee et al. 2019). BOD originates from various point and non-point sources, including the organic-rich discharge from domestic wastewater treatment plants (WWTPs); stormwater runoff from agricultural land, parking lots, and livestock operations; combined sewage outflows (CSOs), and septic and wastewater pipeline leakages. Nitrates are essential nutrients for plants which can cause plants and algae to grow rapidly, and also contribute to high BOD levels and depletion of oxygen (Poulsen et al. 2018). The presence of high concentrations of nitrate-nitrite nitrogen (NN for the remainder of the text) indicates poor water quality. Potential sources include surface drainage, discharge from WWTPs, CSOs, agricultural runoff, decomposition of plant and animal tissue, and the complex nitrogen cycle (Ghane et al. 2016).
Following passage of the Clean Water Act (CWA) of 1972 and the establishment of the U.S. Environmental Protection Agency (EPA), pollutant discharge from point sources was strictly maintained with the promulgation of the National Pollutant Discharge Elimination System (NPDES) which established the effluent standards and monitoring and reporting requirements of individual dischargers. For example, 5-day BOD (i.e., BOD 5 ) is one of the secondary treatment effluent standards for publicly owned treatment works (POTWs) to reduce organic contents discharged into water bodies under the requirements of NPDES permitting. Pollutant discharge from non-point sources are mainly managed by states as assigned by the CWA through various programs and practices, e.g., via promoting best management practices in agricultural production and establishment of educational programs (Shortle et al. 2020). However, BOD and nutrient levels remain high in many freshwater bodies, leading to eutrophication, oxygen depletion, and degraded water quality (Bargu et al. 2019).
Currently, water quality data are collaboratively monitored by local, state, and national agencies such as the U.S. Geological Survey (USGS), U.S. EPA, Indiana Department of Environmental Management (IDEM), and Muncie Bureau of Water Quality (BWQ). Water quality data are readily available from various sources such as the National Stream Quality Accounting Network (NASQAN) of USGS, the Water Quality Surveillance and Response system of the EPA, and organizations at state and local levels. Nationwide coverage in the EPA STORET/WQX and USGS National Water Information System (NWIS) allows for about 1.8 million observation sites in the United States (Beran & Piasecki 2009). However, these data are not analyzed systematically nor in a timely manner, especially at the local level. With the abundance of data, it is important to understand long-term water quality changes to inform decision-makers about the effects of conservation practices and climate change. Nonetheless, the analysis of long-term water quality data is complicated due to the inconsistencies in data collection frequency, accuracy, completeness, detection level change over time, and dilution effects at different discharge levels. Using a simple statistical approach may not be suitable for addressing natural variations and inconsistent data variations for predicting trends of water quality parameters. Weighted Regression on Time, Discharge and Seasons (WRTDS) is a novel approach recently developed by USGS to analyze long-term water quality data by considering long-term trends, seasonal variation, and discharge effects (Hirsch et al. 2010Hirsch 2014). WRTDS is designed for large and long-term datasets to overcome challenges by estimating concentration and flux with and without the influence of discharge variation (Sprague et al. 2011). By looking into concentration and flux trends, WRTDS can also help identify pollutant sources (Murphy & Sprague 2019).
With long-term monitoring data and the employment of WRTDS, the current study aims to investigate the long-term trends and seasonal variations of White River water quality and analyze the sources of organic and nitrogen pollutants upstream of the White River at Muncie.

Study area and data sources
The White River is an upstream headwater of the Mississippi River, flowing through south and central Indiana, which is dominated by agricultural land. The river also passes major Indiana cities, including Muncie, Anderson, and Indianapolis (Crawford 2001). As a result, the river may be adversely affected by various contaminant sources, including point and non-point sources.
This study analyzes the BOD, DO, and NN data of the White River between 2002 and 2018 at five sites collected by the Muncie BWQ. The BWQ of Muncie was established in 1972 and is among the oldest local water pollution testing and enforcement agencies in the United States. Since its establishment, BWQ has monitored water quality variables of the White River around Muncie and contributed greatly to local water quality improvement. However, monitoring is less frequent and occurs at fewer sites in earlier years, and the detection limits have improved over time. This research analyzed the water quality data from 2002 to 2018 because the data are relatively complete and consistent at the sites of interest, and also fulfill the data requirement for the WRTDS method. The five sites are located immediately before and after the city boundaries of Muncie and Yorktown ( Figure 1). Specifically, the Memorial site is the most upstream and the 575 W site the most downstream of the study. The Walnut and Tillotson sites occur near the downtown area of Muncie. The Muncie WWTP is located about 200 m upstream of the Nebo site which is one of the main points pollution sources of the river, and a smaller WWTP, in Yorktown, is located downstream of 575 W. There is also a series of CSO outlets along the White River which release combined sewage directly into the river during storm events ( Figure 1).
The discharge data of the White River at Muncie were downloaded from the USGS NWIS website (USGS station 03347000) and are applied to all five sites.
Data processing and modeling with WRTDS WRTDS combines time expressed in years (the Trend component), time of year (the Season component), and discharge into a single regression equation (Hirsch et al. 2010Hirsch 2014). For any location in time-discharge space (t and Q), WRTDS assumes that the log-transformed concentration of a variable is a function of time, flow, and season, and is expressed in the following form: where c is the concentration of the variable of interest (i.e., BOD, DO, or nitrate); β 0 , β 1 , β 3 , and β 4 are coefficients which change while moving through the time-discharge space; Q is the stream discharge; t is the time of year; and ε is the unexplained variation or random component. In WRTDS, the coefficients are estimated for every combination of Q and t within the period of record. For every combination of Q and t, the coefficients in Equation (1) are estimated using weighted regression, where observations within a half window for each regression are weighted relative to annual, seasonal, and flow distances from the observation at the center of the window. Observations with distances farther from the center (i.e., greater time and different flow values from the center) have less weight during parameter estimation for each regression (Beck et al. 2018). This approach differs from previous trend-analysis methods where fitted coefficients are applied to the entire dataset and avoids the assumption that the flow versus concentration relationship is constant with time. Daily flux is then calculated by multiplying the estimated concentration with the respective daily mean discharge. This three-dimensional matrix is used to estimate the final concentration values for each day over the entire study period. The flux of a parameter estimated in WRTDS is calculated as where f is the expected value of flux in kg/day, c is the expected value of concentration in mg/L, Q is the daily discharge in m 3 /s, and 86.40 is the unit conversion factor. These time series of estimates can then be summarized into time series of monthly averages and then into annual averages. The default model fitting systems are used in this study utilizing the half-window widths approach as 6 months for seasonal weights, 10 years for annual weights, and half the range of flow in the input data for discharge weights (Beck et al. 2018).
Another important feature of WRTDS is the calculation of flow-normalized concentration and flux which reflect long-term trends by eliminating the influences of random discharge variation. Streamflow discharge has a great influence on the concentration and flux of constituents in water. For example, sediment concentration and flux will generally increase with increasing discharge during  stormwater events. Changes of extent and frequency of stormwater events will lead to changes of sediment concentration and flux. From the perspective of environmental conservation, it is of interest to understand how conservation practices at the watershed-scale have affected sediment loading without changes in climate and hydrologic conditions. As such, it is desirable to eliminate the effects of random variations in discharge while maintaining the effects associated with seasonal and longterm trends in streamflow . As a result, WRTDS introduces the concept of flow-normalized concentration and flow-normalized flux on any given day (Equations (3) and (4) where FN is the flow-normalized concept, t is the date of any point in the record T, n is the number of years in the record, C is the concentration, F is the flux, T i is the set of days in the record for a given calendar day, e.g., Jan 1st, and k is a unit conversion factor. Annual flow-normalized concentration and flux are then calculated by taking the mean value of estimated daily values over a water year. Long-term water quality analysis usually suffers from frequent non-detectable values, variable methods used over the analysis period, and inconsistency of data sampling. In this study, for example, the total number of BOD data ranges from about 4,000 samples at the Memorial site to about 200 samples at the 575 W site. There are also non-detectable values which were considered as left-censored data following through survival analysis (Tobin 1958;Moyer et al. 2012;. This study used the R package Exploration and Graphics for RivEr Trends (EGRET) developed by Hirsch & De Cicco (2019) to organize the dataset and run WRTDS. EGRET has built-in functions to extract discharge data directly from the USGS NWIS. This research followed the instructions of EGRET to format water quality data for the package to read and interpret.
The application of WRTDS to describe trends could reveal new insights given the disproportionate effects of physical drivers on water quality. The effects of biological drivers may also be more apparent because hydrological effects can be removed by WRTDS. As such, the application of WRTDS models for trend analysis can facilitate a broader discussion on the need to focus beyond nutrients to develop integrated plans for water quality management. More detailed information about WRTDS could be referred to Hirsch et al. (2010). BOD, DO, and NN were analyzed and discussed over the entire period from 2002 to 2018, and then for each 5-year interval during of 2003-2008, 2008-2013, and 2013-2018, to understand details of changing patterns.

Discharge data processing and modeling
River discharge is a critical component to understand the water cycle and climate system, as the amount of freshwater movement portraits the potential relationship of the water quality. As the study sites are within the proximity of the one USGS station (03347000), the same discharge rate was applied to each site. The analyte of interest is subsequently paired with flow data for each observational date . Streamflow variation over the data collection period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) and the long-term  was analyzed to understand the long-term discharge trend and variation in the White River. This streamflow variation broadens the understanding of the physical environment of water quality variation in this watershed.

Streamflow conditions
Discharge plays a major role in concentration and flux variance. Between 2002 and 2018, the mean annual discharge fluctuated but generally presents a decreasing trend. The mean winter discharge has a similar decreasing trend, while the mean summer discharge increases. However, this trend is different from the long-term discharge change of the White River at Muncie from the early 1930s, which shows a clear increasing trend annually, including during both winter and summer seasons (Figure 2).

Model diagnostics
The observed concentration and flux are plotted against estimated concentration and flux for each variable with a 1:1 line (line of identity) to demonstrate goodness of fit (Figure 3). The plots indicate that simulated results fit well with observed results (Figure 3). While the sites with fewer data are affected by outliers (e.g., Nebo), the observation versus concentration points are evenly distributed along the line of identity, indicating no major bias in the population component of the model. Overall, the flux simulation predicts the observations better than concentration simulation. Figures 3-5 show the diagnostic plots of BOD, DO, and NN. The coefficient of determination (R 2 ) between observations and model estimates for concentration is lower, ranging from an average of 0.38 for BOD, 0.53 for NN, to 0.79 for DO. The R 2 for flux is high, ranging from an average of 0.83 for BOD, 0.86 for NN, to 0.99 for DO. Overall, R 2 reflects a good match between observations and model estimates, especially for flux. Keep in mind that the model is not supposed to perfectly predict observational concentration at any specific date as the model eliminates the effects of discharge and season.

Changes in concentration and flux over time
Flow-normalized BOD concentration increased for all the five monitoring sites between 2002 and 2018 ( Figure 6). Specifically, the flow-normalized BOD concentration increased by 55.0% at Memorial, 19.0% at Walnut, 12.0% at Tillotson, 11.0% at Nebo, and 6.4% at 575 W (Table 1; Figure 6).  Flow-normalized BOD flux also shows an increasing trend for all the monitoring sites over the analysis period. Specifically, the flow-normalized BOD flux increased by 19.0% at Memorial, 57.0% at Walnut, 6.7% at Tillotson, 51.0% at Nebo, and 2.6% at 575 W (Table 1; Figure 6).
The flow-normalized BOD concentration and flux changes for all the five monitoring sites over the 2003-2008, 2008-2013, and 2013-2018 analysis periods (Table 2). In particular, the flow-normalized    (Table 2). The flow-normalized DO concentration decreased slightly in the monitoring sites over the period between 2002 and 2018. Specifically, the flow-normalized DO concentration changed by À0.5% at Memorial, À1.9% at Walnut, À4.0% at Tillotson, and À2.0% at 575 W. However, the DO concentration at Nebo increased by 6.9% during this period (Table 1; Figure 7).    2018 at four of the five sites (Table 3) (Table 3).  Flow-normalized NN concentration increased overall during the period of 2002-2018. Specifically, the flow-normalized NN concentration increased by 13.2% at Memorial,9.3% at Tillotson,32.2% at Nebo,and 27.1% at 575 W, and decreased 11.0% at Walnut over the entire analysis period (Table 1; Figure 8).  The flow-normalized NN flux increased for four out of the five monitoring sites over the analysis period. Specifically, NN flux concentration increased by 5.5% at Memorial, 6.1% at Tillotson,15.2% at Nebo,and 20.9% at 575 W, while it decreased by 14.2% at Walnut (Table 1; Figure 8).   (Table 4).

Overall and seasonal water quality change
The White River at Muncie experiences marked daily, seasonal, and annual variations in discharge, distinct water quality patterns at different locations and sampling dates; therefore, a simple time-series analysis of observational data is not suitable for characterizing water quality and identifying pollution sources. WRTDS provides flow-normalized long-term results and trend analyses and provides insights into the spatial and temporal variation of constituents, and is applied in BOD, DO, and NN analysis for five sites covering the urban area of this agricultural dominated watershed.
The WRTDS analysis of BOD, DO, and NN indicates that water quality regarding nutrient and organic contaminants of the White River has deteriorated over the past two decades. The BOD and NN concentrations and flux experience increasing trends, and DO concentration and flux have been decreasing. This trend has been especially pronounced during the last 5-year segment, i.e., from 2013 to 2018. The increase of NN is consistent with studies for the larger Mississippi River basin-scale which concluded that little consistent progress has been made in controlling nitrate concentration and flux since the 1980s (Sprague et al. 2011).
The BOD and NN increase at the Memorial site (24% and 15%, respectively) during 2013-2018 is fairly significant. Since Memorial is the most upstream of this study area and receives drainage from agricultural land, conservation practices may have underperformed in reducing nutrients and organic materials in streams in recent years. This inference differs from the findings of Koltun (2019), i.e., that changes in flow-normalized NN concentrations and flux were mostly decreasing in the Upper White River. However, the analysis period of Koltun's work was different from this study and only the Walnut site at Muncie was analyzed. Further examination of this study's results shows that before 2013, the NN increase was most notable at downstream sites and was marginal at the Memorial site and decreased at the Walnut site ( Figure 8; Table 4). As a result, it is important to interpret the conclusions from different studies under different research time frame and sites.
The monthly BOD analysis shows that BOD concentration is generally higher in spring and summer and lower in fall and winter (Figure 9). This trend is somewhat inconsistent with the general knowledge that BOD levels tend to be lower in summer and fall due to the increased metabolism of aerobic bacteria in warm temperatures. A possible reason is the frequent heavy rainfall events and high runoff during spring and early summer in the watershed, which carry over huge quantities of organic materials and CSOs into streams. The monthly DO analysis shows that DO concentration is consistently higher in spring and winter, and lower in summer and fall ( Figure 10). The oxygen level is mainly affected by temperature and does not show notable differences between sites. While high BOD could potentially reduce DO, the BOD influence on DO in small streams and rivers like the White River is generally not a determining factor for DO depletion due to the constant flow of water and turbulence.
The monthly NN concentration shows distinct patterns between upstream and downstream sites: values are higher in spring and summer at the upstream Memorial, Walnut, and Tillotson sites; in contrast, concentrations are higher in fall at the downstream Nebo and 575 W sites ( Figure 11). NN at the upstream sites originates mainly from non-point sources, and the high values in spring and early summer are consistent with fertilizer applications in those seasons. Fertilizer application, CSOs, and non-point source runoff are lower in late summer and fall seasons, but effluents from the WWTPs remain high, which taken together, lead to higher NN concentrations at downstream sites  in late summer and fall. The seasonal pattern of NN levels is consistent with a recent study conducted at the Upper White River site (Walnut site of this study) which also revealed the same seasonal patterns of NN at the Walnut site (Koltun 2019).

Effects of stream discharge on water quality
Streamflow variation has a dramatic influence on water quality by affecting the concentration of dissolved and suspended constituents, the degree of erosion, the deposition of solids, and carryover of constituents from runoff. The White River occurs at the headwater region of the Mississippi River Basin, and the dominant source of discharge during dry seasons is base flow. During wet seasons and storm periods, agricultural and urban runoff, rainfall, and CSOs are major contributors to streamflow of the river. August and September are usually the months of lowest river flow.
Contrary to the long-term  increasing trend of annual, summer, and winter discharge, the annual and winter discharge in the White River has decreased over the past two decades (Figure 2). A slight increase in summer discharge occurred from 2002 to 2018, but the increase  was marginal and may be overshadowed by extreme rainfall events. The BOD and NN concentration would have increased with a decreasing discharge assuming no change in BOD and NN loading. DO levels would also be affected by the discharge. Higher volumes of faster-moving water with more turbulence will incorporate more atmospheric oxygen, possibly over-saturating the water with oxygen. In contrast, smaller volumes of slow-moving water may be heated in summer and retain less DO than would cooler water. Overall, the DO concentration and flux have decreased along with a decreased annual discharge from 2002 to 2018 (Figure 7).

Spatial variation and potential sources of contamination
Although overall trends suggest degrading water quality, considerable differences of pollutant concentrations and fluxes occur among the sampling locations. Flow-normalized concentration and flux rates may be used to infer potential pollution sources. Flux increase is more related to the storm events which carry large volumes of non-point source contaminants, while concentration increase is more  related to low flows when point source pollution and groundwater recharge might be the major factors (Hirsch et al. 2010). The major CSOs downstream may be the main point source contributors to BOD and nitrogen contaminants, which is evidenced by the increase of BOD and NN flow-normalized concentration and flux rates at Nebo (Figures 6 and 8). After Nebo, the BOD level decreased to background levels of other sites, while the NN concentration remained at a high level until 575 W. This result could be explained that BOD may be quickly degraded in freshwater, while NN uptake may require longer distances. The Yorktown WWTP is located nearby 575 W, and the treatment capacity of that WWTP is small and its effluent does not markedly affect BOD and NN levels. The increased rate of BOD concentration is highest at Memorial, which implies that the contribution from agricultural runoff increases during low-flow seasons. The BOD flux increasing rate is highest at Walnut and Nebo, which indicates that the contribution of urban runoff and CSOs increases during high-flow seasons. BOD concentration increases as the river flows downstream, and peaks at Nebo after the Muncie WWTP and a major CSO. Data from the Muncie WWTP have shown that BOD levels have been kept under control and are not notably higher than river background levels. As a result, the CSOs along the river near urban regions may be the main contributor of BOD, especially in recent years. The NN concentration increasing rate is highest at Nebo and 575 W, which represents the contribution of the WWTP increases during low-flow seasons. The NN flux increase rate is not as clear as in the concentration increase. The NN from non-point source runoff during high-flow seasons increased slightly. The NN concentration increase does not differ much in the first three sites but is much higher at the Nebo and 575 W sites, and the results indicate that WWTPs effluent contribution is significant during low-flow seasons.

Implications for improved water management
Significant efforts have been made in recent years to reduce nitrogen loading from regional WWTPs, given the disproportionate contribution of nutrients relative to other sources (Beck et al. 2018). WWTPs in Delaware County need to take sufficient action in nutrient reduction, considering evidence that it is a major contributor of nitrogen to the White River especially during the low-flow seasons. CSOs are the major contributors of BOD, and the rapid BOD increase at the Memorial site indicates that agricultural runoff contribution is increasing. In addition, the high NN concentration in spring and early summer at upstream sites also indicate the contribution from agricultural runoff. Together, findings of the study show that water and soil conservation in the upstream watershed must be improved. This study underscores the rapid and multifaceted changes that are facing White River. To improve the overall water quality of the White River, a broader, more holistic assessment of current conditions is needed. Likewise, an understanding and emphasis on preventing potential future degradation are important to the regional planning approach.

SUMMARY AND CONCLUSIONS
Two decades of BOD, DO, and NN data from five sites along the White River near Muncie were analyzed using the WRTDS method. Results show that flow-normalized BOD concentration increased 7-55% over the past 16 years, while flow-normalized DO concentration decreased 1-4% and NN concentration increased 10-33%. The analysis of the sources indicates that agricultural runoff contributes to NN during high-flow seasons, and local WWTPs and major CSOs have significantly contributed to NN inputs in low-flow seasons. Likewise, upstream agricultural runoff and CSOs along the urban corridor have made a big impact on river BOD levels. Seasonal differences, which are naturally occurring, reflect fertilizer contribution to NN levels. The findings from this study should be helpful for watershed managers to re-think conservation practices in the agriculture-dominated watershed; this work has useful implications for water resources and water quality management across a broad geographic domain.

AUTHOR CONTRIBUTIONS
Md S.A.: data processing, manuscript writing, and original draft preparation; B. H.: conceptualization, methodology, data visualization, and manuscript writing and editing; A. G.: conceptualization, reviewing, and proofreading; and J. P.: reviewing and editing.