{'fraction_non_alphanumeric': 0.03626346517556122, 'fraction_numerical': 0.01907737850505715, 'mean_word_length': 4.793830395426394, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
River Water‐Quality Concentration and Flux Estimation Can be Improved by Accounting for Serial Correlation Through an Autoregressive Model

Accurate quantification of riverine water‐quality concentration and flux is challenging because monitoring programs typically collect concentration data at lower frequencies than discharge data. Statistical methods are often used to estimate concentration and flux on days without observations. One recently developed approach is the Weighted Regressions on Time, Discharge, and Season (WRTDS), which has been shown to provide among the most accurate estimates compared to other common methods. The main objective of this work was to improve WRTDS estimation by accounting for the autocorrelation structure of model residuals using the first‐order autoregressive model (AR1). This modified approach, called WRTDS‐Kalman Filter (WRTDS‐K), was compared with WRTDS for six constituents including nitrate‐plus‐nitrite (NOx), total phosphorus, total Kjeldahl nitrogen, soluble reactive phosphorus, suspended sediment, and chloride. Near‐daily concentration records at nine sites were used to generate subsets through Monte Carlo sampling for five different sampling scenarios. Results show that WRTDS‐K provided generally better daily estimates of concentration and flux than WRTDS under these sampling scenarios for all constituents, especially NOx. The degree of improvement is strongly affected by the underlying sampling scenario, with WRTDS‐K gaining more advantage when more samples are available, and hence more residuals can be exploited. The performance of WRTDS‐K depends on the AR1 coefficient (ρ) and that relationship varies with constituents and sampling scenarios. These results provided recommendations on the optimal ρ for each constituent and sampling scenario. Overall, WRTDS‐K has the potential for broad applications to monitoring records elsewhere, as demonstrated by a pilot application to Chesapeake Bay tributaries.


Introduction
Surface water-quality degradation has negative effects on the health of ecosystems and abundance of biological living resources. A fundamental challenge to the effective management of surface water quality is the accurate quantification of the flux of materials (e.g., sediment and nutrients) that are transported through rivers and streams. Constituent flux, representing the mass of materials passing a particular location in a given period of time, is critical for the establishment of ecosystem restoration targets (e.g., total maximum daily loads; U.S. Environmental Protection Agency, 2010), calibration of watershed models (e.g., SPARROW, HSPF; Ator et al., 2011;Ator et al., 2019;Linker et al., 2013;, and evaluation of water-quality trends and drivers (Bowes et al., 2008;Green et al., 2014;Hirsch, 2012;Zhang et al., 2015).
Accurate quantification of constituent concentration and flux is challenging because cost considerations usually dictated that water-quality observations are typically only made at frequencies such as weekly or monthly even though discharge data are available at a daily or more frequent time step. For example, sites in the Chesapeake Bay watershed (USA) are sampled at least 20 times per year, including at least 12 regular monthly samples and at least eight stormflow samples (Chanat et al., 2016). In comparison, concentration data at many other locations are even more sparse and often do not contain targeted stormflow samples.
To estimate complete records of daily concentration and daily flux (the product of daily concentration and daily discharge), a number of statistical approaches have been developed, including but not limited to Cohn et al. (1989), Cohn et al. (1992), Johnes (2007), and Park and Engel (2015). These approaches are based on statistical relations between measured concentrations (the response variable) and a set of explanatory variables that typically include time, season, and discharge associated with the date of sample collection. In virtually all these models, a single regression model developed with all available data in the period of record is used to estimate concentrations and fluxes for unsampled days.
These statistical approaches rely on assumptions of homoscedasticity of model errors and fixed relations between concentration and each covariate, which can be frequently violated in real water-quality datasee examples in Hirsch et al. (2010) and Zhang et al. (2016). To relax such model assumptions, Hirsch et al. (2010) developed Weighted Regressions on Time, Discharge, and Season (WRTDS). Like many of its predecessors, WRTDS uses time, discharge, and season as explanatory variables. However, WRTDS differs in that it functionally develops a separate model for each day of the observed record by re-evaluating the dependencies of concentration on time, discharge, and season using samples most relevant to the day of estimation (Hirsch et al., 2010;. The method also removes the constraint of homoscedastic residuals that is typical of regression-based models. Consequently, WRTDS can better represent the temporally-varying seasonal and discharge-related water-quality patterns and thus produce more accurate estimates (Chanat et al., 2016;Hirsch, 2014;Lee et al., 2016Lee et al., , 2019Moyer et al., 2012). For example, Chanat et al. (2016) compared WRTDS and the seven-parameter LOADEST model (Cohn, 2005;Cohn et al., 1989;Cohn et al., 1992) for 80 sites in the Chesapeake Bay watershed and observed generally greater explanatory power by WRTDS. Since its publication, WRTDS has been adopted in a range of surface waterquality studies in the United States and elsewhere, including nitrogen (Ator et al., 2019;Chanat & Yang, 2018;Sprague et al., 2011;Strickling & Obenour, 2018;Van Meter et al., 2017), phosphorus (Fanelli et al., 2019;Rankinen et al., 2016), sediment (Hirsch, 2012;Zhang et al., 2016), chloride (Corsi et al., 2015;Stets et al., 2018), organic carbon (Zhang & Blomquist, 2018), organic matter (Stackpoole et al., 2017), and major ions (Bird et al., 2018).
Although WRTDS is more likely to compute accurate load estimates than other commonly used regression methods, there are potential opportunities to further improve its model performance. These improvements can be realized by investigating model errors that are not accounted for by the three model covariates (i.e., time, discharge, and season). In this regard, one potential approach is to consider new explanatory variables to further account for variability in sampled concentrations. For example, it is often reported that riverine transport can be affected by watershed antecedent conditions due to the buildup of momentum and activation of new pathways (Davis et al., 2014;Gray et al., 2015;Macrae et al., 2010;Warner et al., 2009). To quantify antecedent conditions, Vecchia (2003) has developed the "flow anomaly" concept, which can be computed based on daily discharge records that are widely available. In the context of WRTDS, this concept was first considered by Murphy et al. (2014) to focus on nitrate transport. Expanding upon the effort of Murphy et al. (2014), Zhang and Ball (2017) developed and tested a set of modified WRTDS models with flow-anomaly variables and reported generally improved estimates of constituent concentrations and fluxes by the modified models.
While it can be useful to consider new explanatory variables within WRTDS, another potential opportunity for improvement is to investigate the autocorrelation structure of the WRTDS model residuals and exploit that information to adjust estimates. In other words, WRTDS estimates could be adjusted through an empirical residual correction approach to better approximate observations. Such a concept has been previously tested in the context of LOADEST (Appling et al., 2015;Aulenbach, 2013;Aulenbach & Hooper, 2006;Verma et al., 2012). These methods are collectively referred to as the "composite method" and have been found to provide more accurate estimates through the empirical residual correction. However, these studies were limited to LOADEST and limited to simple autocorrelation structures for concentration residuals. Specifically, Aulenbach and Hooper (2006) used a piece-wise linear interpolation of residuals, whereas Verma et al. (2012) used rectangular-and triangular-shaped patterns to capture the autocorrelation structure. Similarly, Schwarz et al. (2006) developed the FLUXMASTER method, which also has the capacity of adjusting estimates based on model residuals. However, no time series model has been invoked in this context. Thus, the overall objective of this work is to test the hypothesis that WRTDS estimation can be improved by accounting for the autocorrelation structure of model residuals. Previous studies of time-series residuals often exhibit strong serial correlation at time scales of days and up to many months (Gray, 2018;Hirsch & Slack, 1984;Kirchner & Neal, 2013;. Serial correlation is generally related to specific characteristics of individual hydrologic events (high flows or drought periods) that are not captured by the statistical model being used. These characteristics may include differences in antecedent conditions (see Murphy et al., 2014;Vecchia, 2003) or increased (or decreased) hydrologic contributions from subwatersheds that differ in terms of sources or dominant transport processes. An extreme example of this is the Missouri River flood of 2011, as documented by Kalkhoff (2013).

Water Resources Research
To limit the scope of the work, we focused on the first-order autoregressive (AR1) model for capturing residual autocorrelation. We call this modified approach "WRTDS-Kalman Filter (WRTDS-K)" because it uses the residual information on the sampled days to produce adjusted estimates for the unsampled days. We compared the performance between WRTDS-K and WRTDS for nitrate-plus-nitrite (NO x ), total phosphorus (TP), total Kjeldahl nitrogen (TKN), soluble reactive phosphorus (SRP), suspended sediment (SS), and chloride (Cl). Method performance was compared using near-daily data at nine monitoring sites. For each site and constituent, concentration subsets were generated through Monte Carlo subsampling for five sampling scenarios and used for model implementation and evaluation. Specific goals of this work included the following: 1. To test the hypothesis that WRTDS-K can improve the estimation of concentration and flux of common water-quality constituents, particularly NO x for which discharge is often a rather poor predictor. 2. To search for the optimal value(s) of the lag-1 coefficient (ρ) for different water-quality constituents and sampling scenarios. 3. To apply WRTDS-K (with appropriate ρ) to major tributaries to Chesapeake Bay and compare the WRTDS-K estimates with the regular WRTDS estimates. The results from this research will inform the design and management of riverine water-quality monitoring networks as well as the analysis of water-quality data for concentration and flux estimation. Thus, this work is expected to have broad implications toward traditional water-quality monitoring sites that are distributed across the world.
It should be noted that WRTDS-K is only intended for use when the goal is to provide the best possible estimates on the actual daily values of concentration (or flux) or their mean values aggregated over periods such as months, seasons, or years. The original WRTDS method also produces estimates of the flow-normalized concentration and flux values, which are used in evaluating long-term trends in water quality (Hirsch et al., 2010;. Such trends will continue to be computed using the original model.

WRTDS
The implementation of the typical WRTDS model is illustrated in Figure 1a. The core equation of the model can be mathematically described as follows (Hirsch et al., 2010): where C i is the constituent concentration on day i, in mg/L; Q i is the daily mean discharge on day i, in m 3 /s; T i is the time expressed as decimal year; β 0,i, β 1,i, β 2,i, β 3,i, and β 4,i, are the fitted coefficients that vary smoothly over the model domain; σ i is the fitted value of the conditional standard deviation of the error of the model; and z i is the standardized model residual on day i (standard deviation = 1).
The ln(C)~ln(Q) relationship is adopted for three reasons: (1) it provides a generally better fit than a linear relationship, (2) it eliminates the possibility of obtaining negative estimates for concentration, and (3) the residuals from this model are more nearly normal .
Functionally, WRTDS develops one regression model for each day of the record to estimate concentration and flux. With the known time, discharge, and season of the estimation day, WRTDS prescreens the entire concentration record and selects samples that are sufficiently "close" to the estimation day with respect to time, discharge, and season. The selected samples are used to build a weighted regression model using equation (1) and a unique set of model coefficients is obtained. These coefficients are used to estimate ln(C) on the estimation day by substituting known values of time and discharge, which is then transformed to concentration in real space using a bias correction factor (BCF; equation (3)) The BCF arises because we seek to estimate the conditional mean of concentration, but the model is designed to estimate the conditional mean of the natural log of concentration. This BCF is approximately correct when the errors in log space are normal and sigma is relatively small, and the sample size is large (e.g., >50).
To expedite the estimation process, the above steps are conducted on a grid network formed by t and ln(Q). For the t axis (x axis), grid values are spaced 1/16th of a year apart from the beginning year to the end year of the record. For the ln(Q) axis (y axis), 14 grid values are spaced with equal distance for the discharge range from 5% below the minimum discharge to 5% above the maximum discharge in the record. For each grid point, WRTDS develops a separate weighted regression model, which results in an estimated concentration "surface" as a function of t and ln(Q). Daily concentration is then estimated using a bilinear interpolation of the surface, which is multiplied by daily discharge to compute daily flux. For more details, see Hirsch and De Cicco (2015).

WRTDS-K
WRTDS-K is a variation of the regular WRTDS method for estimating concentration as a function of time, discharge, and season ( Figure 1b). First, the method develops the regular WRTDS model using the available concentration record. Second, the method calculates the residuals in the logarithm space for days with observed concentrations, that is, These residuals are the error in the model, expressed as the observed ln(C) minus predicted ln(C). Third, the method standardizes these model residuals by dividing them by the standard deviation appropriate to that day (see equation (1)), resulting in the standardized residuals, z i : In this work, WRTDS-K assumes that the standardized residuals have a serial correlation structure of AR1 with a lag-1 correlation coefficient (ρ) that falls in the range of (0, 1). It is recognized that AR1 may not be the correlation structure that best represents the behavior of the standardized residuals, but simulation of more complex correlation structures is generally not feasible given the irregular and sparse sampling frequencies of most water-quality data sets. Previous study of selected daily water-quality records have shown that an AR(1) model may understate the correlation at lag times such as a week or a month; however, the authors believe that the AR(1) model is the most feasible to incorporate the residual information and still produce useful results.
WRTDS-K breaks the whole data record into nonoverlapping gap intervals based on the days with observations. For a gap interval, assume that the last observation before the gap is on day 1 and the first observation after the gap is on day n, then there are n-2 days that need to be estimated by WRTDS-K for this gap interval.
Because of the complexity of the process (e.g., log transformations, time-varying model coefficients, and time-varying variances), the expected value for each of the n-2 missing values are estimated using a Monte Carlo simulation with 50 replicates. For each replicate, the method generates the time series of missing values, that is, c 2 ,c 3 ,c 4 ,…. ,c n − 3 ,c n − 2 ,c n − 1 , based on the known values of c 1 and c n , the parameters of equation (1), and the estimate of σ for each day in the gap period. The generation of a single replicate is as follows: 1. Generate n-2 values of e k , which are independent standard normal random variables (mean 0, variance 1). 2. Conditioned on the estimated values of z 1 and z n (determined from the data and equations (2) and (3)), the remaining n-2 z k values are generated based on the recursive relationship of an AR1 process: Note that unlike the usual way that an AR1 process is generated, the generating process used here is conditioned on two known (nonrandom) values, that is, z 1 and z n .When ρ is 0, the generated values are solely based on e k , that is, the independent standard normal random variables. As ρ approaches 1 in the limit, the generated values are solely based on z k , that is, the standardized residuals. This represents the special case where the algorithm conducts a linear interpolation between the observed residuals (i.e., z 1 and z n ), which is equivalent to the composite method (Aulenbach & Hooper, 2006).
3. This process is repeated for every gap in the record. This process is designed so that the entire time series of values of z (including both the values calculated from the data and all the generated values in between them) have a lag-1 correlation coefficient ρ. 4. This time series of z values is transformed to a set of concentration values: The above process is repeated M times and the expected value of concentration for each day i is the mean of the M replicates of C i , which is the expected value of concentration for day i, C i . Note that the BCF, defined in Equation (3), is not used in WRTDS-K. It is used in the original formulation of WRTDS because it estimates a mean value based on the conditional moments of the logarithms. In WRTDS-K, in equation (7), the estimates are for a single value rather than a mean value. The estimate of the mean is computed only after retransforming each of the individual replicate values for the day. To identify an appropriate value for M, we explored the sensitivity of C i to M and found that C i values based on 50 replicates were very close to values based on 1,000 replicates-see a case study in Figure S1 in the supporting information. Thus, 50 replicates were used for all our WRTDS-K model runs in this study. There are data sets elsewhere for which using 50 replicates in the WRTDS-K calculations may not be sufficient to achieve stable results. It can be useful to run the analysis a second time with a different random number seed and determine if the new set of results (when aggregated to monthly or annual totals) differs substantially from the first set of results. If that is the case, then one could increase the number of replicates to several hundred. The computational burden of using a larger number of replicates is relatively small.
For sampled days, the estimate of concentration is equal to the measured value. By contrast, for the regular WRTDS model the estimate of concentration obtained using equation (2) is the unbiased estimate of concentration. Thus, WRTDS-K will always be more accurate than WRTDS for the sampled days. For unsampled days, the estimate of concentration by WRTDS-K is also expected to be different from that by WRTDS due to the step of empirical residual correction. In general, when the sampling is sparse, such as a typical gap of more than 60 days, WRTDS-K should produce estimates that are similar to those determined in the regular WRTDS model. When the typical gaps are small, such as seven days or less, WRTDS-K estimates can be quite different from WRTDS estimates because there are many measured data values that WRTDS-K can use and the serial dependence of the data at short lags (e.g., one to seven days) have increased influence on the estimates.
Another way to conceptualize the operation of WRTDS-K is to consider that on any given day there are three sources of information available to make our daily estimate. One source is the most recent residual from before the day being estimated. The second is the subsequent residual after the day being estimated. The third is zero, which represents the regression model (with no error), which provides the best available representation of water-quality transport processes. The task then is to come up with a single estimate that is the weighted sum of these three bits of information. The weights of the adjacent samples are rather high when they are very close to the sampled days and low when there is a long period before or after them. The approach presented is viewed as a reasonable and objective way to set those weights.

Testing Sites and Data
The development and testing of the WRTDS-K approach was done using near-daily water-quality data for several tributaries to Lake Erie and the Ohio River, which are available through the National Center for Water Quality Research at Heidelberg University, OH, USA (National Center for Water Quality Research, 2015). The data sets cover nine monitoring sites that had at least 15 years of concentration samples, located within the Maumee River, Scioto River, Great Miami River, Sandusky River, River Raisin, Cuyahoga River, Grand River, Honey Creek, and Rock Creek (Figure 2). These sites cover a range of drainage area (34.6-6,330 mi 2 ), agricultural land (30.4%-89.9%), wooded land (7.3%-50.1%), and urban land (0.6%-9.6%). Consequently, the results from this research are expected to be not limited to a single watershed or watershed type within this geographic setting (Table S1 in the supporting information). The Lake Erie and Ohio River tributaries that have been monitored by the National Center for water quality research at Heidelberg University, OH, USA. Details of the selected sites and available water-quality data are summarized in tables S1 and S2, respectively. This region has the following general characteristics: (1) it is located in the north temperate climatic zone, (2) its physiography is dominated by the Central Lowlands of the Interior Plains Division, and (3) it is underlain with thousands of feet of sedimentary rocks formed during the Paleozoic Era (Schiefer, 2002). For characteristics of each tributary, see Schiefer (2002).
For each site and constituent, the original concentration data have between 256 and 349 sampled days per year on average (Table S2). The median of sampled concentrations was used for any day in which multiple concentration values were recorded. Daily river discharge data at each site are available for these sites through the U.S. Geological Survey (USGS) National Water Information System database (U.S. Geological Survey, 2018). Concentration and discharge previously downloaded and processed by Zhang and Ball (2017) were used here.

Data Subsampling
To test the performance of WRTDS-K, we selected five sampling scenarios and performed Monte Carlo subsampling of the original concentration data to produce concentration subsets (30 subsets per sampling scenario) to serve as the "training data" for subsequent model run. These scenarios are as follows: S1: One concentration sample randomly selected for each calendar month for which concentration samples were available (i.e., 12 concentration samples per year), S2: Flow-stratified sampling with 12 regular samples and 8 storm samples per year (i.e., 20 samples per year). Within each specific year of the record, samples with daily discharge exceeding the 90th percentile of daily discharge distribution of the year were marked as "storm sample"; all other samples were considered as "regular samples." Because a similar sampling scenario has been adopted for many sites in the Chesapeake Bay watershed-see details in Chanat et al. (2016), evaluation of model performance for this strategy can provide important insights toward the Chesapeake Bay watershed monitoring network and similar programs elsewhere, S3: Twelve random monthly samples (i.e., S1 above) plus additional eight random samples per year (i.e., 20 samples per year), S4: Two samples randomly selected for each calendar month for which concentration samples were available (i.e., 24 samples per year), and S5: One sample randomly selected for each week for which concentration samples were available (i.e., 52 samples per year).
For illustration, an example of NO x in the Maumee River at Waterville, OH is provided in Figure 3, which shows the original concentration data and the five types of subsamples.
To quantify model performance, concentration and flux values in the original data set that were not selected in the subsamples (section 3.2) were used for model validation. Based on these data, the modified Nash-Sutcliffe efficiency (MNSE) was calculated: where Y obs i and Y est i are the observed and modeled values for the ith day, respectively. Y mean i is the average of observed data in the validation data, and n is the total number of observations in the validation data. MNSE quantifies the absolute difference between the estimated and observed values normalized by the sum of the absolute deviations of the observations. MNSE is less sensitive to large errors than the standard NSE, because it uses absolute difference whereas the latter uses squared difference (Krause et al., 2005;Moriasi et al., 2007). In general, a larger MNSE means a better model and a model with MNSE = 1 implying that the estimates match perfectly with observations.

Water Resources Research
We followed the approach of Zhang and Ball (2017) to repeat the concentration subsampling experiment 30 times for each site-constituent combination. For each subset, we quantified the MNSE for WRTDS and WRTDS-K (with varying ρ) for each of the 30 replicate runs and the median MNSE was recorded and used for subsequent model comparison. This procedure helps reduce the required computational time, while providing a reasonable number of replicates to describe overall uncertainties in the model performance.
The above analysis resulted in MNSE values for WRTDS and WRTDS-K (with varying ρ) for each sampling scenario for each of the 54 site-constituent combinations (i.e., 9 sites × 6 constituents). To more directly compare the results, the ratio between the MNSE values of the two methods (R) was calculated for each case of site, constituent, and sampling scenario: Note that R values greater than 1 indicate that the WRTDS estimation is improved by WRTDS-K. These R values were used to determine the optimal ρ value(s) for the six constituents under different sampling scenarios.
While our comparison focused on MNSE, we used percent bias (PB) as an auxiliary metric for evaluating model performance for the entire period of record: PB could be negative, positive, or zero; a value of zero means that the model is unbiased.

Ancillary Analysis of Continuous-Monitoring Data
To further characterize optimal ρ values, we conducted an independent analysis on continuous-monitoring data at three additional monitoring stations. These stations are North Raccoon River near Sac City, IA (USGS ID: 05482300); Potomac River near Washington, D.C. (Little Falls; 01646500); and Mississippi River at Baton Rouge, LA (07374000), and they vary in 3 orders of magnitude in terms of drainage area. This preliminary analysis was limited to NO x due to data availability. For each site, 15-min measurements of NO x were downloaded from the USGS National Water Information System database (U.S. Geological Survey, 2018). These data were processed as follows: 1. Subsample the 15-min record to a daily record by selecting the values recorded at the time 12:00 of each day, 2. Run the regular WRTDS model on the daily record to obtain daily concentration estimates, 3. Calculate the standardized residuals (z i ) using equations (4) and (5), 4. Extract the longest contiguous sequence (i.e., gap free) in z i , 5. Estimate ρ and the autocorrelation function (ACF) for the extracted sequence using the R function ar (order.max = 1, method = "yule-walker") (R Core Team, 2019).
These estimates of ρ were then compared with the optimal ρ identified using Monte Carlo analysis of the Heidelberg data sets (section 3.3).

Model Application to Chesapeake Tributaries
A pilot application of the WRTDS-K approach was implemented with optimal ρ values on water-quality data collected from major tributaries draining to Chesapeake Bay, the largest estuary in the United States. We focused on the USGS River Input Monitoring Program stations that collectively account for a majority of nutrient and sediment inputs from the nontidal Chesapeake watershed (Table S3 and Figure S2). For these sites, sampling scenario S2 is considered to best approximate actual sampling procedures (Chanat et al., 2016). Water-quality and discharge data at these sites were available through Moyer et al. (2017). Cl was excluded because it had not been a focus of Chesapeake watershed monitoring. For each river and constituent, flux estimates from WRTDS and WRTDS-K were compared on daily, monthly, and annual scales. Mean absolute difference (MAD) was used to quantify the level of difference between the two sets of estimates: where W est i and K est i represent the WRTDS and WRTDS-K flux estimate for time index i, respectively. This metric was computed for annual (where n = number of years × 365), monthly (n = number of years × 12), and annual scales (n = number of years), respectively.

WRTDS-K Versus WRTDS: The NO x Example
To guide readers through our results, we start with a preliminary comparison of the two models by focusing on NO x and a preselected ρ (i.e., 0.95) for WRTDS-K (Figure 4). For NO x concentration (Figure 4a), two key patterns emerged. First, WRTDS-K performed better than WRTDS for almost all sites and sampling scenarios. Second, the degree of improvement achieved by WRTDS-K, represented by the vertical distance of the data points from the 1:1 line, varied considerably among the different sampling scenarios. Among the five sampling scenarios, S1 generally had the least amount of improvement, whereas S5 showed the greatest improvement. Since both scenarios are based on random sampling, this difference is attributable to the different frequencies of concentration samples-that is, 12 versus 52 per year for S1 and S5, respectively. As expected, this result shows that WRTDS-K becomes more accurate relative to WRTDS when more concentration samples are available, and thus, more residuals can be utilized in the AR1 interpolation procedure (equation (6)). For NO x flux (Figure 4b), it is also evident that WRTDS-K performed generally better than WRTDS, but the degree of improvement was lower than in the case of NO x concentration. In addition, the degree of improvement by WRTDS-K was less variable among the sampling scenarios. Overall, this preliminary comparison provides three key messages. First, WRTDS-K can provide improved estimates of NO x 10.1029/2019WR025338 Water Resources Research compared to WRTDS for almost all sites under the five different sampling scenarios. Second, the degree of improvement in concentration estimation is strongly affected by the underlying sampling scenario, with most improvement achieved by S5. Third, the degree of improvement by WRTDS-K is more evident for concentration than flux.

Sensitivity of WRTDS-K Performance to ρ
Here we consider results from all sites and constituents to examine the sensitivity of the WRTDS-K performance to the choice of ρ. For better communication of the results, the median of the computed R values (equation (9)) of all sites was quantified and plotted as a function of ρ. Recall that R values greater than 1 indicate that WRTDS-K estimates have a higher MNSE than WRTDS estimates (equation (9)). Figure 5 shows the site-integrated results for NO x , which reinforce the messages presented above. Generally, R values are greater than 1.0, meaning that WRTDS-K produces more accurate estimates than WRTDS for almost the entire range of ρ. As expected, among all sampling scenarios, WRTDS-K does not show improvement over WRTDS when ρ is zero, since the estimation is not informed by residuals on sampled days. Second, as illustrated in Figure 4, the greatest improvement by WRTDS-K is achieved with sampling scenario S5 and the least improvement with S1. Also, as illustrated previously, the degree of improvement by WRTDS-K is more evident for concentration than flux. The relative improvement by S2 with respect to NO x flux estimation is even greater than that by S4, stressing the value of stormflow sampling in the context of residual utilization by WRTDS-K. Finally, Figure 5 shows that R gradually increases with ρ and the optimal ρ is roughly in the range of 0.9-0.99 for concentration estimates and 0.95-0.99 for flux estimates. This pattern holds true for all five sampling scenarios and shows that the WRTDS-K performance relative to WRTDS is sensitive to the choice of ρ for NO x .  Figure S7) to examine the sensitivity of WRTDS-K performance to ρ. These results demonstrate similar results as presented for NO x and show that the performance of WRTDS-K relative to WRTDS varies with ρ. Importantly, the optimal value (or range) of ρ varies among sampling scenarios and among constituents, which should be considered when developing recommendations on the optimal ρ for WRTDS-K. The patterns of R versus ρ for these constituents are generally different from the pattern for NO x (Figure 5). For example, TP results showed that R gradually increases with ρ between 0 and 0.875 ( Figure S3). Through visual inspection of these plots, the optimal range of ρ was identified to be roughly 0.8-0.875, 0.8-0.9, 0.8-0.925, 0.8-0.925, and 0.8-0.95 for TP, TKN, SRP, SS, and Cl, respectively.
Overall, these results demonstrate that (1) WRTDS-K can lead to improved estimation for a large range of ρ through the utilization of model residuals on sampled days, (2) the most improvement with WRTDS-K is generally achieved with ρ values greater than~0.8, and (3) the optimal ρ varies among constituents and among sampling scenarios.

Recommended ρ for Different Constituents and Sampling Scenarios
While the above site-integrated results provide insights on the optimal ρ, they do not speak to the differences among the nine study sites. In this regard, we developed a simple and objective algorithm to search for the optimal ρ for concentration estimation. For each combination of site, constituent, and sampling scenario, we ranked the R values for the nine sites and selected the ρ values corresponding to the three largest R values. We then counted the occurrence of the selected ρ values for each constituent and each sampling scenario and plotted the counts as functions of ρ and sampling scenario for each constituent (Figure 6).
The recommended values of ρ (Table 1) are further inspected against the R versus ρ plot of each individual site-constituent case (plots not shown). In most cases, the optimal range of ρ is generally consistent with the range listed in Table 1. Also, in most cases, R values are greater than 1.0, indicating improved performance by WRTDS-K than WRTDS. Moreover, the shape of the R versus ρ plots are also comparable among the nine different sites in most cases-that is, R starts at a value of~1.0 at ρ = 0, increases steadily with ρ, and reaches

10.1029/2019WR025338
Water Resources Research a plateau at ρ = 0.7 or 0.8, which may continue through ρ = 0.99 or decrease slightly in the range of ρ = 0.9-0.99, depending on the constituents and sampling scenarios.
Broadly speaking, the recommended values of ρ (Table 1) should be most applicable to monitoring data collected in geological and hydrological settings that are similar to the Lake Erie and Ohio River region. Caution should be exercised if these recommended values are applied to monitoring records obtained in settings that are substantially different from this region. Unfortunately, long records of near-daily water quality Figure 6. Count of ρ values that resulted in the greatest improvement in daily concentration estimation by WRTDS-K, shown as functions of ρ and sampling scenarios. Note that the counts cannot exceed 9 (i.e., nine sites).  To further support the validity of these results, we estimated the AR1 coefficient for continuously monitored NO x records at three selected sites. Figure 7 shows the longest contiguous sequence of standardized residuals for these sites, as obtained using the procedure described in section 3.4. These results show that standardized residuals behave differently among the sites for the selected time periods. For the North Raccoon River, which has the smallest drainage area among the three sites, the standardized residuals shift from generally positive values to generally negative values at about the 200th day of the 439-day period (Figure 7a) and the ACF exhibits a long memory (i.e., high correlation even at a lag of 30; Figure 7b). For the Mississippi River, which has the largest drainage area among the three sites, the standardized residuals remain generally positive in this 266-day period (Figure 7e) and the ACF switches from positive correlations at small lags to

10.1029/2019WR025338
Water Resources Research negative correlations at large lags ( Figure 7f). For these two sites, ρ is estimated to be 0.93 and 0.956, respectively, which are consistent with the recommended range for NO x that were derived from the Monte Carlo analysis (Table 1 and Figure 6a). Standardized residuals from the Potomac River near Washington D.C. site have positive values on most days in the sequence except for several events (Figure 7c). The associated ACF exhibits a long memory (>30 lags) and the estimated ρ is 0.868 (Figure 7d). On the one hand, this analysis demonstrates that the recommended range of ρ for NO x (Table 1) is generally supported with these independent, continuous-monitoring data sets. On the other hand, this analysis suggests that the AR1 model may not always be the most appropriate assumption for capturing the ACF of the standardized residuals, since the correlation can persist at lags much longer than one day. As noted above, simulation of more complex correlation structures is generally not feasible given the irregular and sparse sampling frequencies of most water-quality data sets. Analysis of the correlation structure of model residuals of water-quality sensor data is an important topic for future research, which is only beginning to be possible as long-term sensor records become available.

Summary of Improvement by WRTDS-K (With Optimal ρ)
With the optimal ρ identified for the different sampling scenarios, the improvement by WRTDS-K relative to WRTDS can be quantified accordingly. Figure 8 compares the MNSE of WRTDS and WRTDS-K (with optimal ρ) for concentration for the six constituents. The optimal ρ value used here corresponds to the largest count in Figure 6. As expected, WRTDS-K generally shows a higher MNSE than WRTDS. (For comparison on flux estimates, see Figure S8.) Table 2 summarizes the percent improvement in MNSE by WRTDS-K relative to WRTDS. For simplicity, only the median values of the nine study sites are presented, which are 56-119%, 13-36%, 23-61%, 16-45%, 10-20%, and 21-59% for NO x , TP, TKN, SRP, SS, and Cl, respectively ( Table 2). The lower and upper ends of these values correspond to S1 and S5, respectively. Comparing the six constituents, the smallest improvement by WRTDS-K is associated with SS, whereas the greatest improvement is associated with NO x . Again, these results demonstrate that WRTDS can be further improved by WRTDS-K through the exploitation of the sample-associated residuals, especially for monitoring records with high sampling frequencies.
The two models were briefly compared using the percent bias (PB) metric (equation (10)) using optimal ρ values identified above. The results from this comparison are summarized in Figures S9 and S10. For daily concentration estimates, WRTDS-K is less biased than WRTDS in almost all cases (i.e., five sampling scenarios × nine study sites) for TP, TKN, and especially SRP ( Figure S9b-S9d) and a majority of cases for NO x ( Figure S9a). For SS, WRTDS-K appears to be less biased than WRTDS in about one half of the cases but more biased in the other half ( Figure S9e). For Cl, both methods are nearly unbiased in all cases ( Figure S9f). For daily flux estimates, the two approaches appear to have similar yet small biases in all cases for NOx, TKN, and Cl (Figures S10a, S10c, and S10f). In addition, WRTDS-K appears to be less biased than WRTDS in most cases for TP, SRP, and especially SS flux (Figures S10b, S10d, and S10e).

Model Application to Chesapeake Tributaries
In a pilot application, we implemented WRTDS-K (with optimal ρ) on water-quality data in the rivers draining to Chesapeake Bay. The goal of this effort was to evaluate the difference of the estimated fluxes from WRTDS-K and WRTDS, which is already used for flux estimation at these sites. These differences are evaluated at daily (Figures S11-S14), monthly (Figures S15-S18), and annual (Figures S19-S22) time scales.
We illustrate the differences using a case study of NO x in the Potomac River (Figure 9). On the daily scale, concentration estimates from WRTDS appear to follow the pattern of regression toward the mean and fail to capture the extreme values ( Figure 9a). Notably, the WRTDS estimates were consistently lower than the observations for several months after the peak concentration observed in early 1996, which is associated with one of the biggest flood events during the period of water-quality data collection. Concentration estimates of the regular WRTDS model were persistently below the true values for a couple of months after the flood event. This suggests that the flood event (caused by a large regional rainfall and warm weather that happened at a time of a large regional snowpack) may have mobilized nitrate sources that are not typical of winter conditions. This anomaly, of several months' duration, is captured well by WRTDS-K but not by the regular WRTDS model.

Water Resources Research
The differences between the WRTDS and WRTDS-K estimates depend on the temporal scales at which estimates are made (Figures 9b-9d). On the daily scale, the two sets of estimates have considerable differences, but generally scatter around the 1:1 line (Figure 9b). Differences in daily estimates produced a MAD (equation (11)) of 5,470 kg/day for the entire period. At a monthly scale, the two methods have fewer large differences (Figure 9c), and thus a smaller MAD (4,613 kg/day) than daily estimates. At an annual Figure 8. MNSE of WRTDS-K (with optimal ρ) versus MNSE of WRTDS for daily concentration estimates for the six constituents for the nine Lake Erie/Ohio River sites. The gray-shaded area indicates MNSE WRTDS-K ≥ MNSE WRTDS . All panels share the same scale for both axes.

10.1029/2019WR025338
Water Resources Research scale, the two sets of estimates are almost always on the 1:1 line (Figure 9d), corresponding to a MAD (2,738 kg/day) that is even smaller than that of the monthly estimates. This scale dependence was also observed with all the other site-constituent combinations (Table S4)-that is, in all cases, the MAD value gets smaller from the daily scale to the monthly scale (i.e., MAD d /MAD m > 1) and gets smaller from the monthly scale to the annual scale (i.e., MAD m /MAD y > 1).
Overall, these results show that WRTDS-K and WRTDS estimates tend to be more comparable at coarse scales (e.g., annual) and that they become more distinguishable on short-term scales (e.g., daily), which is related to the residual adjustment in WRTDS-K. Such short-term based differences may have important implications to the modeling, research, and management efforts that rely on riverine flux estimates, including the calibration of Chesapeake Bay watershed and estuarine models, and exploration of terrestrial drivers of riverine flux patterns.  Constituent  S1  S2  S3  S4  S5   NO x  56  60  69  72  119  TP  13  20  17  19  36  TKN  23  34  29  37  61  SRP  16  20  21  28  45  SS  10  14  11  11  20  Cl  21  30  31  31  59 Note. For simplicity, the median improvement of the nine Lake Erie/Ohio River sites is presented.

Conclusions and Recommendations
This work provided statistical evidence to support the hypothesis that WRTDS daily estimates of concentration and flux can be improved by accounting for the autocorrelation structure of model residuals. Our results show that the WRTDS-K method generally produces more accurate estimates of daily concentration and flux values than WRTDS for all constituents, especially NO x . The degree of improvement by WRTDS-K is strongly affected by the underlying sampling scenario; WRTDS-K estimates tend to be more accurate with increased water-quality sampling because more residuals can be exploited in the AR1 interpolation procedure. The performance of WRTDS-K depends on the selection of an appropriate ρ value, which can vary depending upon the constituent being estimated and the sampling scenario being deployed. Model estimates were used to recommend optimal ρ values for different water-quality constituents and sampling scenarios. In a pilot application, the optimal ρ values were used to quantify nutrient and sediment fluxes from major tributaries to Chesapeake Bay, where it was shown that estimates from the two approaches were comparable at relatively coarse scales (e.g., annual) but became more distinguishable on short-term scales (e.g., daily).
Overall, the WRTDS-K method has the potential for broad applications for estimating concentration and flux from sparsely collected water-quality monitoring records, particularly those from geological and hydrological settings similar to the Lake Erie and Ohio River region. From a monitoring perspective, our results also demonstrate that flux estimation can benefit from both increased sampling frequency and targeted stormflow sampling. Software to compute WRTDS-K daily estimates of concentration and flux are available at http://usgs-r.github.io/EGRET/articles/WRTDS-K.html. This software is designed as an add-on to the existing EGRET software and it also includes some graphics capabilities to help visualize the differences between the original WRTDS estimates and the WRTDS-K estimates.
While this work demonstrates the potential of WRTDS-K, we recognize several areas that are appropriate for future research. First, it would greatly expand the applicability of the approach if the AR1 coefficient can be directly estimated from water-quality data, which are typically irregularly sampled. Second, an alternative approach is to add a cross-validation procedure in WRTDS-K to identify the optimal ρ for any specific data set and use that ρ value to run WRTDS-K. However, this may be difficult to implement because of the low sampling frequencies of traditional water-quality monitoring data. Third, we assumed in this work that WRTDS model residuals follow an AR1 model, but that may not always be the most appropriate model. In this regard, future efforts may consider more complex time series models such as the autoregressive moving average model or the autoregressive fractionally integrated moving average model. However, scientific advancement is needed to directly estimate the model coefficients for these models for the irregularly sampled water-quality data and to develop the analytical approach to fill the gaps between adjacent samples. Finally, as high-frequency monitoring sites and records become more available (Fovet et al., 2018;Halliday et al., 2015;Outram et al., 2014;Pellerin et al., 2012;Pellerin et al., 2016), WRTDS-K has the potential to estimate missing days when the sensors become inoperable, but this would require further research to evaluate the usefulness of WRTDS-K and its optimal settings for varying gap lengths within continuous records. Moving forward, scientists and researchers need to advance the understanding on these aspects and to translate the understanding to software tools that can be broadly applied.