Generalized Likelihood Uncertainty Estimation and Markov Chain Monte Carlo Simulation to Prioritize TMDL Pollutant Allocations

The basic premise of GLUE is that there is not a single optimal set of parameters for any given model (i.e., equifinality), so multiple sets of parameters can be used to satisfactorily represent a watershed response (Beven and Binley 1992). In GLUE, Monte Carlo simulation is used by generating multiple sets of model parameters from parameter-specific PDFs that are determined a priori. In the majority of previous GLUE applications, a priori or prior parameter-specific PDFs were uniform. When using GLUE, the model is run using a set of parameters or a parameter vector $\theta$ sampled from prior parameter-specific PDFs. A likelihood weight is calculated for each parameter vector. Likelihoods can be computed using any goodness-of-fit metric that compares observed data and model-simulated response variables (Stow et al. 2007). For this application, likelihoods were a function of the variance of the residuals [Eq. (1)]. Hydrologic calibration compared daily flow volume to assess model performance, while water quality calibration compared log-transformed simulated daily average and observed instantaneous FC concentrations. In-stream FC concentration can generally vary by orders of magnitude and, therefore, log transformation reduced the influence of extreme FC concentrations on estimated uncertaintywhere $L_e$ is the likelihood; ${\sigma }_e^2$ = variance of the residuals or mean square error; $n$ = number of data points; $Y_i$ = observed data; $Q_i$ = simulated data; and $N$ is the shaping parameter chosen by the user. Eq. (1) has been used frequently with other GLUE applications, such as Beven and Binley (1992). As the value of $N$ increases, the magnitude of difference between the likelihood values of parameter sets with similar variance increases. Different values of $N$ can lead to different confidence intervals (Ratto et al. 2001). This application used an $N$ of 2 because the resulting model outputs bracket the observations (Mishra 2011).

All the parameter vectors that yield acceptable likelihood values, as determined subjectively by the analyst, are considered to be behavioral and are retained. Parameter vectors that produce unacceptable simulations are deemed nonbehavioral and are rejected. Previous applications of GLUE have reported a wide variety of parameter vector rejection criteria. For example, Beven and Binley (1992) considered parameter vectors with very low likelihood values to be nonbehavioral. In this study, the likelihood threshold used to differentiate behavioral and nonbehavioral parameter vectors was the inflection point on a cumulative distribution function (CDF) plot of likelihoods. The likelihood values of the retained parameter vectors are normalized to unity. The normalized likelihood values can be treated as a probabilistic weighting function for the simulated variables and can be used to assess the uncertainty associated with the simulations. The parameter-specific plots of the normalized behavioral likelihoods versus parameter values, also called dotty plots, defined the cumulative probability and posterior parameter distribution (Beven and Binley 1992). The normalized likelihoods for different parameter values are multiplied to the prior probability to calculate the posterior distribution of the model input parameters. As demonstrated in this application, the posterior distribution of parameters can be used to propagate uncertainty in the model outputs. The posterior PDFs provide a description of parameter uncertainty informed by the observed data, thereby leading to a better estimation of model output uncertainty.

MCMC is a broad class of algorithms which are used for sampling from complex distributions. To quantify the uncertainty in these parameter estimates, this study sampled from the Bayesian posterior [$\pi (\theta |Y)$] distribution of the formwhere $L(Y|\theta )$ is the likelihood function that expresses the dependence between observed data and unknown parameters $\theta$; and $\pi (\theta )$ denotes the prior distribution. Similar to GLUE, the MCMC likelihood is based on the squared deviation between the observed data ($Y_i$) and the simulation model ($Q_i$). Because the unknown, estimable parameters ($\theta$) control the physical specifications of the simulations, $Q_i=Q(Y_i|\theta )$ is used in Eq. (3) for clarity. Under this notation, likelihood can be expressed asEq. (3) assumes that the residuals between observed and simulated values are normally distributed. The resulting posterior distribution has this same form, up to proportionality, because $\pi (\theta )\propto 1$. Given the probability models sampling from the posterior requires an iterative algorithm, the Metropolis–Hastings algorithm (Hastings 1970), one of the most general MCMC variants, samples from $\pi (\theta |Y)$ by initializing ${\theta }^{(0)}$ and iterating over ($i=1,\dots ,k$):

The variable $\alpha$ is the acceptance probability for transitioning from our previous iteration to our proposed state and can take on any value between 0 and 1. Given that the proposal density function $g(\theta |{\theta }^{[i-1]})$ can be used to generate samples using $\pi (\theta |Y)$, this algorithm can be used to produce a sequence of random samples (${\theta }^{[1]},{\theta }^{[2]},\dots ,{\theta }^{[k]}$). Eventually, after a sufficient number of iterations (denoted $B$ for the burn-in period), these samples (${\theta }^{[B+1]},{\theta }^{[B+2]},\dots ,{\theta }^{[k]}$) will be distributed according to the stationary target $\pi (\theta |Y)$. Note that in this application, $\theta$ is a correlated, multidimensional parameter vector. Due to the correlations of parameters within the vector $\theta$, attempting to update every element of this vector for each iteration of the Metropolis–Hastings algorithm requires a large number of iterations for convergence. Therefore, in order to improve efficiency, this study applied the Metropolis within Gibbs (Gelfand and Smith 1990; Geman and Geman 1984) algorithm, which sequentially updates only single elements within $\theta$ at each iteration. Hence, for each scalar quantity, parameter update, the proposal density function ($\mathit{g}[·|·]$) is Gaussian, centered at the previous sample with the tunable scaling factor $s$, a Gaussian random walk. Choosing too large or too small of $s$ will result in poor convergence. This study used the Gelman–Rubin statistic (Gelman et al. 2014) to tune and assess Markov chain convergence. The statistic is based on generating multiple Markov chains and calculating the mixture of chain variance, within chain variance and Gelman–Rubin statistic. The Gelman–Rubin statistic is the ratio of the mixture of chain variance and within chain variance multiplied by a correction factor. The Markov chain is considered to be converged when these variances stabilize and the Gelman–Rubin statistic approaches unity. The Gelman–Rubin statistic is useful in determining the convergence of a MCMC algorithm by monitoring the variances, both between and within $m$ independent Markov chains. Letting $B$ and $W$ denote between and within chain variances for $m$ independent Markov chains, each of length $k$, these variances are computed aswhere elements $\overline{{\theta }_j}$, $\overline{\theta }$, and ${\theta }_{ij}$ represent the mean of $j$th chain, the total mean, and the $i$th sample draw from chain $j$, respectively. The convergence is determined by analyzing the number of iterations until the following ratio, $R$, approaches unity

To implement the MCMC approach using the HSPF model, a Visual Basic software utility was developed to populate the HSPF, fixed-format, user control input (UCI) file with a new parameter vector for each simulation; the utility also analyzed the model output prior to the next iteration in the Markov chain. The utility also checked whether each new sampled parameter value was within the bounds of the prior PDF. If not, another parameter value was selected before going through the acceptance or rejection process. Each parameter vector and the simulation results were stored in a SQL Server database.

A probabilistic approach was used to compare a suite of feasible TMDL allocations. Selected uncertain watershed model parameters were characterized with PDFs (instead of deterministic values) and the uncertainty was propagated to the estimated model outputs and the TMDL using Monte Carlo simulation. An ensemble of simulated daily average in-stream FC concentration time series was used to determine quantiles to estimate the confidence intervals. A time series was generated for a lower and upper quantile to estimate the corresponding confidence interval. Two quantile sets were chosen (2.5%–97.5% and 10%–90%) to estimate the corresponding confidence interval (i.e., 95% and 80%). The water quality criterion, instantaneous violation rate (IVR) or exceedance rate for each quantile time series, was calculated by dividing the number of violations (days when the criterion was exceeded) by the number of days in the simulation period. The confidence interval (difference between the lower and upper quantiles) can also be considered to be a measure of reliability of or confidence in achieving the water quality criterion. In the proposed framework, an optimal TMDL allocation is one that (1) minimizes water quality criterion violations (lowest IVR); and (2) demonstrates the greatest reliability (narrowest confidence interval). The latter criterion means that the TMDL allocation performs better in a variety of plausible modeling scenarios (i.e., multiple parameter sets). Typically, in TMDL development, no formal calculation is performed to quantify uncertainty. Instead, a fraction of the total point and nonpoint sources is included as the MOS. In the uncertainty analysis framework presented here, selected parameters in HSPF are characterized with PDFs (instead of deterministic values) and the uncertainty associated with those parameters is propagated to model outputs (e.g., FC concentration) and the TMDL. The probabilistic analysis framework presented here avoids arbitrarily assigning a MOS, by providing a systematic approach to include uncertainty in TMDL development. Comparing alternative TMDL allocations probabilistically allows stakeholders to make more informed decisions between alternative TMDL allocations by clearly understanding the tradeoff between the specified pollutant load reductions, which directly influences implementation costs and the risk of violating the water quality criterion.

Benham et al. (2004) developed the Mossy Creek bacterial TMDL using HSPF. The HSPF model is a comprehensive, lumped parameter model that continuously simulates the movement of water, sediment, and other water quality constituents on pervious and impervious surfaces, in soil profiles, and within streams and well-mixed reservoirs using three main modules PERLND, IMPLND, and RCHRES (Bicknell et al. 2005). The PERLND module represents pervious land, the IMPLND module represents impervious surfaces in which no infiltration occurs, and the RCHRES module represents the stream reaches and reservoirs in a watershed. When using HSPF to develop bacterial TMDLs, bacteria are typically simulated as a planktonic constituent. When modeling indicator bacterial with HSPF, the user typically specifies different build-up and wash-off relationships for the PERLND and IMPLND land uses. The daily FC loading rate to the land surface is specified by the user for each PERLND and IMPLND segment as a constant load (ACQOP) or a monthly variable load (MON-ACCUM). The concentration of FC in groundwater and interflow is specified for HSPF as a constant or monthly amount. Bacteria die-off on the land surface is represented indirectly by providing an asymptotic limit of bacteria build-up (SQOLIM-PERLND). The user specifies the asymptotic bacterial build-up limit for each PERLND and IMPLND land use. This limit is typically calculated using a first-order die-off rate. This asymptotic limit can vary dynamically by PERLND and IMPLND. Release of bacteria in overland runoff or FC wash-off potential is controlled by the WSQOP parameter that specifies the amount of runoff needed to wash off 90% of FC built up on the land surface. In-stream die-off is simulated using a temperature-dependent first-order relationship (Chick’s law).

For this study, Mossy Creek was delineated into eight subwatersheds using the digital elevation model with spatial resolution of 30 m from the USGS (2002). The watershed was monitored monthly near the watershed outlet by the Virginia Department of Environmental Quality (VADEQ) between July 1992 and March 2003 for FC concentration and other water quality constituents. The Department of Biological Systems Engineering at Virginia Tech monitored Mossy Creek semi-monthly between February 1998 and December 2001 near the VADEQ station for selected water quality constituents, including in-stream FC concentration. Daily flow data were also collected from May 1998 to December 2002 at the same site. Hydrologic calibration and validation, and water quality calibration was conducted using data from both monitoring sites. The other data sets required to simulate in-stream FC using HSPF include rainfall, land use, FC loading rates by cattle and wildlife, inflows from springs, solar radiation, and temperature. Model and data development are described in Benham et al. (2004).

Using HSPF to simulate FC bacteria requires information about several hydrologic and water quality parameters. GLUE and MCMC procedures require the modeler to define Bayesian prior PDFs of these model parameters. Many previous studies, such as Beven and Freer (2001) and Makowski et al. (2002), assigned a uniform PDF to most input parameters to avoid bias. For this investigation, assigned model parameter prior PDFs were based on the literature (Al-Abed and Whiteley 2002; Lawson 2003; USEPA 2001), expert opinion, and GIS data available for the Mossy Creek watershed. The uncertain hydrologic parameters were categorized into three groups on the basis of how each parameter potentially varied with land use and time. Table 1 lists the parameters that were not considered to be a function of land use or time. A uniform distribution was the most obvious choice for most parameters, because these are commonly used calibrated parameters, with no specific guidance, except the typical and possible limits in BASINS Technical Note 6 (USEPA 2000). The lower and upper limits of the uniform PDFs correspond to the typical minimum and maximum limits for these parameters from USEPA (2000). The interflow recession coefficient (IRC) was assigned a triangular PDF with lower limit, mode, and upper limit of 0.5, 0.7 (recommended value to start the model calibration process), and 0.8, respectively.

The second group was those parameters that were varied with land use and/or time (Table 2). The INFILT parameter is a function of land use only. An analysis of the histograms from gridded GIS soil and land use data suggested that a triangular PDF was appropriate for all the land uses except the loafing lot (an area where livestock linger as part of the farming operation) land use (Table 2). For those parameters that vary with time, i.e., including nominal upper zone soil moisture storage (UZSN), interception storage capacity (CEPSC), and index to lower zone evapotranspiration (LZETP), parameter PDF values (minimum and maximum) for the month of January were derived from USEPA (2000). The PDFs for the other months of the year were calculated by multiplying the January PDF by a monthly adjustment factor (Table S1) based on the TMDL (Benham et al. 2004).

Simulation of in-stream FC concentrations with HSPF requires the estimation of multiple water quality parameters including land-based constituent loading parameters (ACQOP and SQOLIM), the rate of surface runoff that results in a 90% wash-off in one hour (WSQOP), the first-order die-off rate for FC bacteria in the waterbody (FSTDEC), and the FSTDEC temperature correction coefficient (THFST). The FC loading rates depend on several factors, including species-specific feces production rates and fecal densities, die-off rates, animal density, and the fraction of time livestock are confined (Zeckoski et al. 2005). As cited in TMDL reports and in the literature (ASAE 2003; Geldreich 1978; Yagow 2001), the FC production rates, typically reported in colony forming units per day ($\mathrm{cfu}{\mathrm{day}}^{-1}$), for dairy cattle, beef cattle, and poultry can vary by several orders of magnitude. According to the Mossy Creek bacteria TMDL report, dairy cattle, beef cattle, and poultry are responsible for more than 94% of FC production in the watershed (Benham et al. 2004). Therefore, this study hypothesized that the uncertainty in production rates of dairy cattle, beef cattle, and poultry (present on the cropland and pasture) were likely to have the greatest impact on FC quantification uncertainty. To account for this uncertainty, the loading rates for pervious land areas were assigned a log-triangular PDF. The mode of the ACQOP PDF was determined from values generated using the Bacteria Source Load Calculator (BSLC), a software tool used in FC source characterization for TMDL development (Zeckoski et al. 2005), and the lower and upper limits of the ACQOP PDF were determined by multiplying the PDF mode by 0.1 and 10, respectively (Table 3).

Because the application of manure to cropland varies by month, the FC accumulation rate (MON-ACCUM) was varied monthly. To account for the temporal distribution of ACQOP, the FC cropland loading PDF for January was multiplied by an adjustment factor for each month (Table S1), with the factor developed using the trend of FC accumulation values generated by the BSLC for each month. Deterministic values for ACQOP calculated using the BSLC were used for the other land uses [forest; low-density residential (LDR); high-density residential (HDR); farmstead; and impervious areas].

The SQOLIM parameter is typically calculated by multiplying ACQOP by an adjustment factor of nine, based on the assumption that the die-off coefficient for FC bacteria on pervious land surface is $0.051{\mathrm{day}}^{-1}$ (Zeckoski et al. 2005). Crane and Moore (1986) reviewed several studies and reported a bacteria die-off rate ranging from 0.04 to $0.20{\mathrm{day}}^{-1}$, which translates to the SQOLIM adjustment factor of 2.5–11.5. The SQOLIM adjustment factor was assigned a uniform PDF between 2.5 and 11.5 and was used to calculate the SQOLIM values for each land use (Table 3).

There is no guidance available on estimating WSQOP, FSTDEC, and THFST when simulating FC bacteria with HSPF. The values of these parameters are generally adapted from previous studies and further calibrated. A review of values used in previous FC TMDLs shows a range of 0.5–2.4 for WSQOP (Lawson 2003). Based on these reported values, a uniform PDF bounded by these two values was used for WSQOP for all land uses. In a review by Bowie et al. (1985), FSTDEC values ranging from 0.12 to $2.52{\mathrm{day}}^{-1}$ were reported for various streams. The average of the reported values was $1.1{\mathrm{day}}^{-1}$. A similar value of FSTDEC has been used in TMDLs developed in Virginia, so a triangular PDF was assigned to FSTDEC with a mode of 1.1 and limits of 0.12 and $2.52{\mathrm{day}}^{-1}$. In this study, it was assumed that any uncertainty in THFST would be masked by the uncertainty in FSTDEC. Hence, a deterministic value of 1.07 was used for THFST (Lawson 2003).

Fecal coliform bacteria directly deposited in a waterbody is input as an hourly time series to HSPF; this study used the BSLC-generated time series for the Mossy Creek TMDL. Because the FC discharge from the sole point source was negligible and the FC production by wildlife and straight pipes (unlawful sewers emptying directly into the stream) was estimated to be less than 1% and 2%, respectively, of the total FC direct deposit, the uncertainty in direct deposit FC is primarily due to cattle. To be consistent with the other PDFs of bacteria loads, the cattle direct deposit PDF was assigned a log-triangular distribution. To define this distribution, the cattle direct deposit time series was multiplied by a factor from a log-triangular PDF with a mode of 1 and limits of 0.1 and 10 (Table 3).

The Mossy Creek model was calibrated and validated using the GLUE and MCMC techniques independently. For the GLUE application, this is a two-stage process in which the hydrologic calibration and validation is performed first, followed by water quality calibration and validation. Hydrologic calibration used the period from September 1998 to December 1999. The posterior parameter PDFs obtained using GLUE were used to validate the hydrologic model for the period from January 2000 to September 2002. Posterior PDFs of water quality parameters were estimated using GLUE for the period from October 1998 to December 2001. However, insufficient observed data prevented water quality validation (Benham et al. 2004).

For the MCMC application, each sampled parameter set is a function of the previous set, making separate hydrologic and water quality calibrations impractical. The Mossy Creek model for the MCMC application was calibrated using four years of concurrent hydrologic and water quality data (September 1998 to September 2002). In each iteration, two likelihood functions were calculated: one for hydrologic response based on the daily flow volume at the watershed outlet and one for water quality that compared observed instantaneous in-stream and simulated daily average FC concentrations. A log transformation was performed on the observed and simulated flow volume and in-stream FC bacteria before deriving the likelihood function to normalize the scale of residuals.

The Mossy Creek TMDL (Benham et al. 2004) lists six pollutant load allocations, with only two meeting the FC water quality criterion of zero violations (S5 and S6, Table 4). These allocations were run on a 3.5-year period that includes a range of representative hydrological events in Mossy Creek. The S5 and S6 allocations both require 100% reduction in the FC loading from the loafing lot and straight pipes, 98% reduction from pasture, and 95% from residential pervious land segments. The differences between the two allocations are the reductions in cattle direct deposit (99% versus 94%), loading from cropland (90% versus 95%), and wildlife direct deposit (30% versus 0%).

Benham et al. (2004) recommended the S6 allocation because it calls for a lesser reduction in wildlife direct deposit than the S5 allocation. The TMDL allocations that focus on the reduction of sources from human activities are generally favored. No effort was made in the Mossy Creek TMDL study to consider uncertainty when prioritizing allocation alternatives. The comparison of the S5 and S6 allocations presented in Table 4 considered only selected sources of uncertainty. Uncertainty in wildlife direct deposit was not considered owing to the limited production of FC bacteria by wildlife ($<1\%$ of total FC load).

Using GLUE, the posterior parameter-specific distributions were derived and were used to model Mossy Creek watershed. The prior and posterior CDFs of two hydrologic parameters, lower zone nominal soil moisture storage for pasture (LZSN-pasture) and fraction of infiltration into deep aquifers (DEEPFR) are presented in Fig. 2. Fig. 2(a) shows that the acceptable likelihoods for the LZSN-pasture occur at the lower end of the distribution, which results in the posterior distribution diverging from the prior. Fig. 2(b) illustrates that the acceptable likelihoods of DEEPFR occur for the entire range of the parameter, which results in a posterior distribution similar to the uniform prior distribution. The difference in prior and posterior distributions for the two parameters implies that the observed data provided greater information about LZSN-pasture than DEEPFR. The model output is more sensitive to LZSN-pasture than DEEPFR. A similar analysis was performed for all the uncertain hydrologic parameters in this study (Tables 1 and 2).

The posterior distributions of the hydrologic parameters were used to conduct Monte Carlo simulations for the validation period (January 2000 to September 2002). The daily flow volume resulting from each model parameterization was used to calculate the statistics typically used to assess HSPF calibration sufficiency (Lumb et al. 1994). The 2.5% and 97.5% quantiles for these measures were used to validate the posterior PDFs (Table 5). These dimensionless measures test the model performance in terms of mass balance, high- and low-flow distribution and seasonal variability. All the hydrologic calibration objectives, except 50% lowest flows and 10% highest flows, were met 95% of the time when using the posterior distributions. A comparison of the expert statistics produced by prior and posterior PDFs indicate that including expert statistics in addition to the residuals between observed and simulated daily flow volume in the likelihood evaluation can improve model calibration. The posterior distributions obtained from GLUE were used for uncertainty analysis.

The model performance in terms of bacteria simulation was also evaluated using various goodness-of-fit measures, including absolute error of the geometric mean, average, and median of the simulated FC concentrations. Kim et al. (2007) suggested that the absolute error between observed and simulated values for these measures should be less than 100% for bacteria calibration. Here, the absolute errors were 3.1%, 5.9%, and 10.8%, respectively. The IVR (i.e., FC concentration $>400\mathrm{cfu}/100\mathrm{mL}$) of the observed data is 60.0%, while the simulations exhibit a rate of 50.6%. Since the difference is smaller than 10.0%, this criterion is also satisfactorily met (Kim et al. 2007). In addition, 91.1% of the observed data fall between the minimum and maximum simulated hourly values on the day the observation was made. The observed FC concentration ranged over many orders of magnitude ($20–160,000\mathrm{cfu}/100\mathrm{mL}$), which could be due to various reasons such as the highly varying presence of livestock in the creek. The observed bacteria data are also subject to uncertainty (Hantush and Chaudhary 2014). Based on this assessment, we concluded that the model is sufficiently calibrated for bacteria modeling.

Following the hydrologic calibration and validation as well as water quality calibration, GLUE was used to determine the posterior distributions of water quality parameters. The likelihood function for each water quality parameter was calculated using log-transformed observed FC bacteria in-stream concentrations and simulated daily average FC bacteria concentration for the days for which observed data were available.

Fig. 3 shows the prior CDF and the GLUE-generated posterior CDFs for two water quality parameters: FSTDEC and ACQOP-pasture. The posterior CDFs differ from their prior distributions, implying the influence of observed data on the prior distributions. However, the difference is not as large as that for the hydrologic parameters, which is likely because water quality data were not measured as frequently as streamflow data. The posterior CDF of all water quality parameters were developed similarly.

The posterior distributions obtained via GLUE were used to conduct Monte Carlo simulation for the two Mossy Creek TMDL allocations S5 and S6 (Table 4). The simulations were conducted for a period of approximately 3.5 years (1,218 days). The daily average simulated FC concentration time series from each Monte Carlo simulation was used to determine the 2.5%, 10%, 90%, and 97.5% quantiles (80% and 95% confidence intervals) for the two valid allocations S5 and S6 (Fig. 4). Scenario S6 exhibited greater uncertainty than S5 (i.e., larger difference between the lower and upper quantiles). The S6 scenario reduced cattle direct deposit less than the S5 scenario, but allowed more FC loading from cropland, illustrating that cattle direct deposit was a greater source of uncertainty than cropland FC loading in the Mossy Creek simulation, given that both sources had similar input uncertainty (section “Uncertain Parameters”).

The ensemble average IVR (1.2% versus 1.3%) is comparable for the two TMDL allocations S5 and S6, suggesting that both allocations have a similar chance of improving water quality (Table 6). While the 80% confidence interval for S5 and S6 are similar, the 95% confidence interval for scenario S6 is wider than for S5. As a result, the probabilistic analysis based on the GLUE application suggests that S5 should be the preferred Mossy Creek TMDL allocation because S5 provides similar water quality improvement with greater reliability when compared with S6. This is not in agreement with the deterministic analysis by Benham et al. (2004). The inconsistency between deterministic and probabilistic frameworks when prioritizing TMDL allocations has also been documented in previous research (Borsuk et al. 2002; Langseth and Brown 2011). However, the focus of decision making in the probabilistic application is on water quality improvement and reliability; technical feasibility and cost of implementation was not considered as a decision criterion. With respect to the two alternative allocation scenarios, S5 will likely be costlier than S6 and less feasible to implement. More expensive because S5 requires a greater reduction in cattle direct deposits (i.e., more streambank livestock exclusion fencing to achieve the required reduction), while S6 requires a larger cropland bacteria load reduction that can be achieved in a variety of less costly ways, and less feasible because S5 requires a greater reduction from cattle direct deposits than S6 in addition to reductions in wildlife bacteria loading. The implementation of either of these scenarios will be challenging given the large load reductions that are needed to achieve the applicable water quality criteria.

Markov chain Monte Carlo was used to generate posterior distributions of hydrologic and water quality parameters for the Mossy Creek TMDL model to estimate the uncertainty associated with in-stream FC concentrations. After burn-in ($\sim 70\text{,}000$ iterations), 30,000 MCMC iterations were used to infer both hydrologic and water quality parameters. The distribution of hydrologic parameters exhibits less variation within the parameter space compared with the water quality parameters. The observed data had a greater effect on the hydrologic parameters LZSN-cropland and DEEPFR than the water quality parameters ACCUM-pasture and FSTDEC because the observed data for flow were available for the entire simulation period (1,218 days), whereas observed FC concentration data were available for only 90 days of the period.

When using the Gelman–Rubin statistic for assessing convergence, at least three independent Markov chains are required. To obtain these Markov chains, three independently initiated Markov chains were executed in parallel. Markov chain convergence is determined when the Gelman–Rubin statistic is approximately unity and estimates of the variance stabilize. This analysis was performed for all hydrologic and water quality parameters. All parameters converged around 70,000 iterations or earlier. Posterior distributions for the LZSN-cropland, DEEPFR, ACCUM-pasture, and FSTDEC parameters are shown in Fig. 5.

The posterior CDFs (solid lines in Fig. 5) differed from the prior CDFs (dashed lines in Fig. 5) showing the effect of observed data on the parameters. In general, the posterior CDFs of the hydrologic parameters differed more from the prior CDFs than did those of the water quality parameters. As expected, the GLUE analysis yielded similar results, but there are inconsistencies between the two techniques. The posterior CDFs for all hydrologic and water quality parameters were derived similarly. Using the posterior CDFs, Monte Carlo simulations were performed to yield average daily in-stream FC bacteria concentrations. As with the GLUE analysis, model output was used to generate 2.5%, 10%, 90%, and 97.5% quantiles (80% and 95% confidence intervals) for the Mossy Creek TMDL allocations S5 and S6 (Fig. 6).

The IVR for different quantiles for the S5 and S6 TMDL allocations (0.9% versus 1.1%) are similar suggesting that both allocations have a similar chance of improving water quality (Table 7). However, both the 80% and 95% confidence intervals for scenario S6 are wider than the respective bands for S5. As was the case with the GLUE analysis, this analysis suggests that S5 should be the preferred Mossy Creek TMDL allocation alternative since it provides similar water quality improvement with greater reliability when compared with S6. Similar to what was found in GLUE application, the optimal alternative is not in agreement with the deterministic analysis by Benham et al. (2004).