Two-Phase Monte Carlo Simulation for Partitioning the Effects of Epistemic and Aleatory Uncertainty in TMDL Modeling

A two-phase Monte Carlo simulation (TPMCS) uncertainty analysis framework is used to analyze epistemic and aleatory uncertainty associated with simulated exceedances of an in-stream fecal coliform (FC) water quality criterion when using the Hydrological Simulation Program–FORTRAN (HSPF). The TPMCS framework is compared with a single-phase or standard Monte Carlo simulation (SPMCS) analysis. Both techniques are used to assess two total maximum daily load (TMDL) pollutant allocation scenarios. The application of TPMCS illustrates that cattle directly depositing FC in the stream is a greater source of epistemic uncertainty than FC loading from cropland overland runoff, the two sources specifically targeted for reduction in the allocation scenario. This distinction is not possible using SPMCS. Although applying the TPMCS framework involves subjective decisions about how selected model parameters are considered within the framework, this uncertainty analysis approach is transparent and the results provide information that can be used by decision makers when considering pollution control measure implementation alternatives, including quantifying the level of confidence in achieving applicable water quality standards. DOI: 10.1061/(ASCE)HE.1943-5584.0001731. This work is made available under the terms of the Creative Commons Attribution 4.0 International license, http://creativecommons.org/licenses/by/4.0/. Author keywords:Water quality modeling; Uncertainty analysis; Two-phase Monte Carlo simulation; Total maximum daily load (TMDL); Non-point-source pollution management.


Introduction
According to the USEPA (2001), billions of dollars are required to develop and implement total maximum daily loads (TMDLs) in the United States.Watershed-scale water quality models are widely used to develop TMDLs.According to the ASCE-EWRI TMDL Analysis and Modeling Task Committee (2017) and Borah et al. (2018), commonly used models for TMDL development include the Soil and Water Assessment Tool (SWAT) (Neitsch et al. 2011) and the Hydrological Simulation Program-FORTRAN (HSPF) (Bicknell et al. 2005).The latter is often used to develop TMDLs (Singh and Woolhiser 2002).In spite of the potentially significant costs associated with developing TMDLs, limited science-based guidance is available for estimation of the uncertainty associated with a TMDL.Without some measure of uncertainty, the probability of achieving a given water quality criterion or the risk of violating it cannot be accurately assessed (Hession et al. 1996b).
A margin of safety (MOS) is included in the TMDL by making conservative assumptions or adding a fraction of the TMDL to account for the uncertainties (Franceschini and Tsai 2008;Reckhow 2003).This approach is not scientifically reliable because the MOS is often determined arbitrarily (Dilks and Freedman 2004;Shirmohammadi et al. 2006).The National Research Council (NRC 2001) suggested that the USEPA should end the process of accounting for the uncertainty by using an arbitrarily selected MOS, and instead establish a more rigorous, less subjective approach.Although suggested by many agencies and researchers, the current TMDL practice is still based on the MOS approach without appropriately and adequately considering uncertainty (Ames and Lall 2008;Chaudhary and Hantush 2017;Gronewold and Borsuk 2009).
In general, two types of uncertainty exist in the simulation of a natural system: (1) epistemic or knowledge uncertainty, which is due to the limited knowledge of the processes being modeled; and (2) aleatory or natural uncertainty, which refers to the inherent natural variability within the system (Beck 1987;Hall and Solomatine 2008;Suter et al. 1987).Whereas the latter cannot be reduced, the former can be reduced through additional information (Der Kiureghian and Ditlevsen 2009).The ability to quantitatively estimate the contribution of these two categories of uncertainty can add power to a modeling analysis because an estimation of the impact of aleatory uncertainty provides the modeler with a benchmark that defines the best model performance achievable by using the complete information available in a given data set (Gong et al. 2013).Superimposing both types of uncertainty can be misleading and can yield inaccurate inferences and mask important information (Cullen and Frey 1999).Therefore, differentiating these two types of uncertainty during a computer-based simulation can potentially reduce prediction errors and yield useful decision-support information.
In water quality modeling, uncertainty can stem from the input data sets (climatic, hydrologic, geologic, geomorphologic, and pedologic), spatiotemporal resolution, model parameters and structure, and the observed data used to calibrate/validate the model (Cullen and Frey 1999;Harmel et al. 2006;Harmel and Smith 2007;Nguyen and Willems 2016;Shirmohammadi et al. 2006;van Straten 1998;Willems 2008).Although uncertainty analysis has advanced with respect to hydrologic modeling, there has been much less effort to quantify the uncertainty of complex water quality models.Most of the past efforts to analyze the uncertainty of water quality models were performed for sediment and nutrients (Gardner et al. 2011;Han and Zheng 2016;Mitsova-Boneva and Wang 2007;Park and Roesner 2012;Raat et al. 2004;Van Griensven and Meixner 2007;Wu et al. 2006;Zhang and Yu 2004;Zheng and Han 2016;Zheng and Keller 2007b); the uncertainty associated with predicting in-stream bacteria has not received much attention.Paul et al. (2004) used first-order error analysis to estimate the contribution of sensitive parameters to the fraction of variance in simulated in-stream fecal coliform (FC) concentrations using HSPF.They recommended Monte Carlo simulation to evaluate the uncertainty in HSPF-based in-stream FC simulation because of the inherent limitations of the first-order variance methods (e.g., assumptions of linearity and values of input variances).Shirmohammadi et al. (2006) applied Monte Carlo simulation, Latin hypercube sampling (LHS) (McKay et al. 1979), and generalized likelihood uncertainty estimation (GLUE) (Beven and Binley 1992) for uncertainty analysis of SWAT-based simulation of several pollutants.Jia and Culver (2008) and Mishra et al. (2018) assessed the uncertainty of HSPF-based predicted in-stream FC via GLUE.To date, studies have focused on assessing overall uncertainty, examining the influences of both epistemic and aleatory uncertainty combined.Rode et al. (2010) suggested that novel uncertainty analysis techniques are crucial to rigorously evaluating the uncertainty of water quality models, but no attempt has been made to apply techniques capable of separating the effects of uncertainty types.Partitioning the effects of epistemic and aleatory uncertainty sheds insight into the model applications by providing a more coherent, informative, and transparent picture of the role of uncertainty (Helton 1994;Macintosh et al. 1994;Merz and Thieken 2005).Although there have been efforts to tackle the complex task of assessing uncertainty in other environmental areas such as flood frequency analysis (Merz and Thieken 2005) and hydrologic modeling (Gong et al. 2013), no effort has been made to separately assess uncertainties in the context of developing TMDLs, despite the recommendations of previous research (Langseth and Brown 2011).Adequately assessing uncertainty serves decision makers with a sound basis for more-informed water quality management decision-making, which could potentially lead to a more achievable TMDL.
This study uses two-phase Monte Carlo simulation (TPMCS) (Helton 1994;Macintosh et al. 1994) to independently evaluate the epistemic and aleatory uncertainty associated with simulated in-stream FC concentrations using HSPF.The results are compared with the standard or single-phase Monte Carlo simulation (SPMCS) to illustrate the benefits of TPMCS.The uncertainty analysis framework presented here has the potential to allow stakeholders to develop achievable TMDLs and make moreinformed pollution control-measure implementation decisions.This is the first application of a technique to separate the epistemic and aleatory uncertainties in a watershed-scale water quality model application.

Uncertainty Analysis Techniques
Two uncertainty analyses were performed: (1) single-phase Monte Carlo simulation (SPMCS), which combines the effects of the epistemic and aleatory uncertainty; and (2) two-phase Monte Carlo simulation, which separates these two uncertainty types.In the application of both techniques, the focus was on model parameter uncertainty.Other sources of uncertainty such as model structure and observed data were not considered.

Single-Phase Monte Carlo Simulation
Standard Monte Carlo simulation or the SPMCS is a method in which repeated simulations of the model are performed using randomly selected parameter sets.Each parameter involved in the uncertainty analysis is characterized via a probability density function (PDF).In each random scenario, the SPMCS randomly populates parameter values from predetermined PDFs for each included parameter and executes the model.This process is repeated for a number of simulations sufficient to converge on an estimate of the PDF of output variables (Gardner and O'Neill 1983).In other words, a complex probabilistic method is solved by dividing it into multiple deterministic scenarios (Madani and Lund 2011).The combined results can be aggregated to obtain mean values, confidence intervals, and probabilistic estimates.
In SPMCS applications, a covariance among the uncertain parameters can be used to improve the sampling efficiency depending upon the model and existing knowledge about the parameters.However, in this type of analysis, parameters are often assumed to be independent.This assumption has been made by many researchers, such as Wu et al. (2006).Although some (e.g., Stow et al. 2007) have suggested that considering the correlation between the parameters may improve the accuracy, others (e.g., Omlin et al. 2001) disagree.For this study, the parameters were assumed to be independent and covariance relationships were not considered.

Two-Phase Monte Carlo Simulation
Two-Phase Monte Carlo simulation propagates the epistemic and aleatory uncertainty separately based on the methodology proposed by Helton (1994) and Macintosh et al. (1994).In a TPMCS procedure, model parameters are classified as either epistemic or aleatory.Parameters in which knowledge is limited or in which insufficient field data are available to estimate their values are considered epistemic.Values for these parameters are typically obtained through a model-calibration process.On the other hand, aleatory parameters are generally derived from observed data.Some parameters can be grouped as both epistemic and aleatory, implying that the categorization depends on the circumstances and is decided by the modeler, and is therefore somewhat subjective (Der Kiureghian and Ditlevsen 2009).For the HSPF application presented here, the parameters that were estimated using GIS and field data were considered to be aleatory, whereas the parameters that could be only estimated through model calibration were considered to be epistemic.The rationale was that if a model parameter can be measured, the uncertainty associated with that parameter does not stem from lack of knowledge, and therefore that parameter is classified aleatory.In contrast, when a model parameter cannot be or is not measured, its value is typically determined through model calibration, and the associated uncertainty stems from a lack of knowledge.Therefore, such a parameter is classified epistemic.Classifying the type of uncertainty can depend on factors such as data availability and spatial scale.The potential uncertainty type classification for selected HSPF parameters used in FC bacteria modeling application presented in this study is shown in Table 1.
To clarify TPMCS, suppose that a model has parameters a, b, x, and y, of which the parameters a and b are epistemic and parameters x and y are aleatory.To perform a TPMCS analysis, m sets of epistemic parameters (a and b) are generated by randomly sampling from their predefined PDFs (Fig. 1).For each set of a and b, a set of n random values is generated for the aleatory parameters

Study Watershed
Mossy Creek, located in Rockingham and Augusta counties in Virginia (Fig. 2), was selected for this research.Mossy Creek is an upland tributary of the Potomac River contributing to the Chesapeake Bay via the Shenandoah River, the South Fork Shenandoah River, and the North River.Mossy Creek was listed as impaired in 1996 due to violations of the single-sample FC criterion, and a TMDL was developed for Mossy Creek by Benham et al. (2004).With an area of 40.8 km 2 , the watershed is characterized as a rolling valley with the Blue Ridge Mountains to the east and the Appalachian Mountains to the west.The predominant land uses in Mossy Creek are pasture, forest, and cropland, which cover more than 96% of the watershed.Other land uses include  farmstead, low-density residential (LDR), high-density residential (HDR), and loafing lot.
The sources responsible for FC direct deposit in Mossy Creek were livestock (poultry, sheep, goats, horses, and cattle), wildlife, straight pipes (illegal discharging of sewage waste directly from homes), and one permitted point source.The primary sources of FC identified in the Mossy Creek TMDL were direct deposition of feces in the stream by cattle (cattle loitering and defecating in the stream, i.e., direct deposit) and runoff from pastures where grazing animals defecate.The cattle were present on six dairy farms and one permitted concentrated animal feeding operation. Mossy

Watershed Model Implementation
HSPF was used to develop the Mossy Creek bacterial impairment TMDL (Benham et al. 2004).It is a comprehensive, continuous, lumped-parameter, watershed-scale model that simulates the movement of water, sediment, and a wide range of water quality constituents on pervious and impervious surfaces, in soil profiles, and within streams and well-mixed reservoirs (Bicknell et al. 2005).HSPF is composed of three main modules: PERLND, IMPLND and RCHRES.The PERLND module represents pervious land, the IMPLND module represents impervious surface area where no infiltration occurs, and the RCHRES module represents the stream reaches and reservoirs in a watershed.When using HSPF to develop TMDLs for bacteria impairment, FC is typically simulated as a planktonic (dissolved) constituent.When modeling FC with HSPF, the user typically specifies different build-up/washoff relationships for PERLND and IMPLND land uses.The daily FC loading rate to the land surface is input by the user for each PERLND and IMPLND segment as a constant input (ACQOP) or varying monthly (MON-ACCUM).The concentration of FC in groundwater and interflow can be input to HSPF as a constant or a monthly variable.FC 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 value for each PERLND and IMPLND land use.This limit is typically calculated assuming a first-order die-off rate relationship.This asymptotic limit can vary by PERLND, IMPLND, and time.Release of bacteria in overland runoff or FC wash-off potential is controlled by the WSQOP parameter, which specifies the amount of runoff needed to wash off 90% of FC.In-stream die-off is modeled using a temperature-dependent first-order relationship (Chick's law).
In this application, HSPF was run at an hourly time step and yielded time series for various hydrologic (e.g., flow) and water quality (e.g., FC concentration) outputs.The model has been successfully used to develop many FC impairment TMDLs in Virginia [e.g., Benham et al. (2004), Mostaghimi et al. (2004), andYagow (2001)].Mossy Creek was delineated by Benham et al. (2004) into eight subwatersheds using the digital elevation model from the USGS National Elevation Dataset.The other data sets required by the model include rainfall, land use, direct FC loading from cattle and wildlife, inflows from springs, solar radiation, and temperature.TMDL modeling data development and acquisitions are described in the Mossy Creek TMDL (Benham et al. 2004).The Bacteria Source Load Calculator (BSLC) (Zeckoski et al. 2005) was used to calculate the FC loading rates for PERLND and IMPLND land uses, and the FC direct deposit loads at an hourly time interval.

Uncertain Hydrologic Parameters
Simulation of FC by HSPF requires information about several hydrologic and water quality parameters.USEPA (2000) described the hydrologic parameters and provided typical values and plausible ranges for all the hydrologic parameters.Typically, the values of these parameters are refined through the model calibration.Al-Abed and Whiteley (2002) and Lawson (2003) listed several hydrologic parameters that are typically calibrated and thus considered highly sensitive.This subset of parameters, described in subsequent sections, was used for the SPMCS and TPMCS analyses.For the TPMCS, the uncertain parameters were further classified into aleatory and epistemic as described in subsection "Two-Phase Monte Carlo Simulation (TPMCS)."Aleatory Hydrologic Parameters.The index to mean infiltration rate (INFILT), which assigns the portion of available soil moisture from precipitation (after interception) into surface runoff, is the only model parameter considered to be aleatory because it can be estimated using GIS soil data.Again, the selection of a single aleatory parameter does not imply that other parameters could not be classified as aleatory.INFILT was classified as the only aleatory parameter to illustrate the TPMCS procedure in this example.Because INFILT can be estimated using GIS data, a histogram of INFILT values for 30 × 30-m cells for each Mossy Creek watershed land use was plotted to estimate land use-specific INFILT PDFs.The resulting INFILT histograms suggested a triangular PDF for all the land uses except loafing lot (Table 2).A uniform PDF, the best choice in the absence of additional information (Beven and Binley 1992), was used for the loafing lot land use.The parameters of the PDF were set as the range of the observed loafing lot INFILT.Epistemic Hydrologic Parameters.The hydrologic parameters that were considered epistemic are listed in Table 3.The variation of these parameters with land use and time were considered to be minor in the study watershed, and thus it was assumed that these parameters were not a function of land use or time.Aside from the interflow recession coefficient (IRC), all the parameters listed in Table 3 were assigned a uniform PDF.The lower and upper limits of the PDFs correspond to the typical minimum and maximum limits for these parameters from USEPA (2000).For the IRC, although the typical upper limit is 0.7, it was also chosen as the starting value for the calibration process and can be taken to be the mode value.Therefore, the IRC was assigned a triangular PDF with lower limit, mode, and upper limit of 0.5, 0.7, and 0.7, respectively.
Table 4 lists the PDFs of the epistemic parameters that vary with land use and time.Again, upper and lower parameter limits of the PDFs were assigned following USEPA (2000) guidance.The parameter PDF limits for the month of January are listed in Table 4. Parameter values for the other months were calculated by multiplying the January PDF by a monthly adjustment factor that was generated based on similar TMDLs developed by experts' opinions (Benham et al. 2004).This implies that the distribution of uncertainty over time was not explored.
Uncertain Water Quality Parameters Simulation of in-stream FC concentrations with HSPF requires the estimation of several water quality parameters.The parameters which have been reported as sensitive in the literature and are typically calibrated in water quality modeling, were taken as uncertain.All these uncertain parameters were considered epistemic, including the ACQOP (FC loading rate to the land surface), SQOLIM (asymptotic limit of FC accumulation on the land surface), FC loading rate directly deposited in the streams (direct deposit time series), WSQOP (FC wash-off potential), first-order decay rate for FC in the waterbody (FSTDEC), and the FSTDEC temperature-correction coefficient (THFST).
The FC loading rates depend upon several factors, including species-specific feces production rates, species-specific fecal densities, die-off rates, animal density, and the fraction of time livestock are confined (Zeckoski et al. 2005).As cited in various TMDL reports and in the literature (ASAE 2003;Geldreich 1978;Yagow 2001), the FC production rates in colony-forming units per day (cfu 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, it was assumed that the uncertainty in production rates of dairy cattle, beef cattle, and poultry (in the cropland and pasture land uses) have the greatest impact on FC concentration uncertainty in the Mossy Creek simulation.To incorporate the uncertainty in the FC application rates for the appropriate land uses, the loading rates for pervious land areas were assigned a log-triangle PDF.The mode of the ACQOP values were the values obtained via BSLC and the lower and upper limits of the Note: LDR = low-density residential; and HDR = high-density residential.a Numbers in parentheses show lower limit, mode, and upper limit of the triangular PDF, respectively.b Numbers in parentheses show lower and upper limits of the uniform PDF, respectively.b Numbers in parentheses show lower limit, mode, and upper limit of the triangular PDF, respectively.ACQOP PDF were determined by multiplying the PDF mode by 0.1 and 10, respectively.
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, FC cropland loading PDF for January (Table 5) was multiplied by an adjustment factor for each month.This factor was developed using the trend of FC accumulation values obtained by the BSLC for each month.Deterministic values of ACQOP calculated via the BSLC were used for the other land uses (forest, LDR, HDR, farmstead, and impervious areas).
The SQOLIM parameter is typically calculated by multiplying ACQOP values by a factor of 9.This SQOLIM adjustment factor is based on the assumption that the die-off coefficient for FC on pervious land surface is 0.051 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 day −1 , which translates to a SQOLIM adjustment factor of 2.5-11.5.The SQOLIM adjustment factor was assigned a uniform PDF between 2.5 and 11.5, and it was used to calculate the SQOLIM values for each land use.
There is no guidance available on estimating parameters WSQOP, FSTDEC, and THFST when simulating FC via 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 between these two values was used for WSQOP for all land uses.Bowie et al. (1985) reported FSTDEC ranging from 0.12 to 2.52 day −1 for various streams.The average of the reported values was 1.1 day −1 .Similar FSTDEC values were used in several 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 day −1 .For 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).
FC directly deposited in a waterbody was input into HSPF as an hourly time series (generated by the BSLC) for the Mossy Creek TMDL.Because the FC discharge from the sole point source was negligible and the FC production by wildlife and humans (e.g., straight pipes) was estimated to be less than 1% and 2%, respectively, of the total FC direct deposit, the uncertainty in direct deposit FC was assumed to be primarily due to cattle.To be consistent with the other bacteria load PDFs, the cattle direct deposit load PDF was assumed to be log-triangular.A log transformation was performed because cattle direct deposit load values were large.To obtain this PDF, the cattle direct deposit time series was multiplied by a factor that had a log-triangular PDF with a mode of 1 and limits of 0.1 and 10.

Hydrologic Calibration and Validation
The hydrologic calibration period was September 1998-December 1999 and the validation period was January 2000-September 2002 (Benham et al. 2004).The HSPF expert system (HSPEXP) tool (Lumb et al. 1994) was used to assist with model calibration.HSPEXP is a standard tool for HSPF calibration and has been successfully used in many TMDL studies such as Chin (2009) and Jia and Culver (2008).The calibration objectives were to meet all the expert statistics suggested by Lumb et al. (1994).This meant minimizing the error in predicted total volume, 50% lowest flows, 10% highest flows, storm peaks, seasonal volume, and summer storm volume (Table 6) based on the daily average flow volume and rate time series output.
To calibrate HSPF, the TPMCS was conducted with 300 epistemic and 40 aleatory iterations (i.e., 12,000 HSPF simulations).The guidance by Lumb et al. (1994) was used to adjust the limits of the parameters' PDFs during the calibration.Care was taken not to violate the maximum possible parameters limit values suggested by Bicknell (USEPA 2000).The process was repeated until the model performance was satisfactory.After calibration, the same number of iterations was conducted for both SPMCS and TPMCS to validate the model.The number of aleatory parameters, INFILT for each land use, was far less than the number of epistemic parameters, and therefore a greater number of simulations was required to sample the parameter space of epistemic parameters.The selection of 300 and 40, respectively, was considered to give an adequate convergence and acceptable balance to sample the epistemic and aleatory parameters effectively.Table 6 presents the values of the expert statistics for the calibration and validation periods.The PDFs of the parameters obtained following the calibration are given in Table 7.For parameters varying with month and land use (MON INTER-CEP, UZSN, and LZETP), the values of the remaining months and land uses were obtained by multiplying a predetermined factor by corresponding values of forest land use in January.These PDFs were used for both the SPMCS and TPMCS application.Log-triangular (1.12 × 10 11 , 1.12 × 10 12 , 1.12 × 10 13 ) a Cropland (for January) Log-triangular (2 × 10 6 , 2 × 10 7 , 2 × 10 8 ) a SQOLIM adjustment factor All Uniform (2.5, 11.5) b SQOLIM-PERLND All ACQOP-PERLND (for each land use) times SQOLIM adjustment factor WSQOP-PERLND All Uniform (0.5, 2.4) b FSTDEC (day −1 ) All Triangular (0.12, 1.1, 2.52) c a Numbers in parentheses show logarithm of lower limit, mode, and higher limit of the log-triangular PDF, respectively.b Numbers in parentheses show lower and upper limits of the uniform PDF, respectively.c Numbers in parentheses show lower limit, mode, and upper limit of the triangular PDF, respectively.

Water Quality Calibration
For water quality calibration, a TPMCS consisting of 300 epistemic and 40 aleatory iterations (12,000 HSPF simulations) was conducted for the period October 1998-December 2001 by using the PDFs of hydrologic parameters from calibration (Table 7).
The output from HSPF runs included ensembles of daily maximum, daily minimum, and daily average FC concentration time series.These ensemble values for each day were averaged for all 12,000 simulations and plotted against the observed data alongside the single-sample FC water quality criterion (400 cfu=100 mL) (Fig. 3).Because the data were observed by collecting a sample once per day, it cannot be expected that the simulated daily average FC concentration will exactly match the observed data.However, it is reasonable to assume that the observed data will fall between the minimum and maximum simulated values for a specific day.For the Mossy Creek watershed, 72.2% of the observed data fell between the average minimum and maximum FC concentrations.According to Kim et al. (2007), the percentage of observed values within a 5-day minimum-maximum range should be greater than 70%, and thus the model was found to be sufficiently calibrated for FC simulations.The percentage of observations above the simulated minimum-maximum range was slightly greater than that below the range, suggesting that overall, the model tends to slightly overestimate the FC concentration.As is often the case when modeling FC, there were insufficient observed data to permit model validation.

Alterative TMDL Allocation Scenarios
A given TMDL scenario allocates the pollutant load among different sources, and hence suggests the amount of reduction in pollutant loading from each source needed to meet the applicable water quality criteria.To simulate the Mossy Creek TMDL allocation scenarios in HSPF, a reduction factor was applied to the pollutant load from each source.The Mossy Creek TMDL (Benham et al. 2004) listed six pollutant allocation scenarios, with only two meeting the FC water quality criteria (S5 and S6; Table 8).Both the S5 and S6 scenarios require 100% reduction in the FC loading from straight pipes (i.e., unlawful sewer connections emptying directly into the stream).The major difference between the two scenarios is the reduction in cattle direct deposit, wildlife direct deposit, and loadings from cropland.
Of these two successful TMDL allocation scenarios, S6 was recommended through the deterministic analysis of Benham et al. (2004).The S6 scenario was recommended because it calls for a lower reduction in wildlife direct deposit than does S5.TMDL allocation scenarios that focus on the reduction of anthropogenic sources are generally favored.No effort was made to incorporate a consideration of uncertainty into the TMDL allocation scenario prioritization.A comparison of the S5 and S6 scenarios was repeated here to investigate whether the consideration of selected sources of uncertainty in allocation scenarios would lead to the selection of the S5 scenario instead of S6.Because of the low production of FC by wildlife, uncertainty in wildlife direct deposit was not considered.A 3.5-year period that included a range of representative hydrologic events in the Mossy Creek was selected to simulate the in-stream FC concentration under the two allocation scenarios.Two statistical analyses, a two-sample Kolmogorov-Smirnov (KS) test (Massey 1951) and an ANOVA, were also performed to further analyze the effect of epistemic and aleatory uncertainties in the two TMDL allocation scenarios.The KS test determines whether two CDFs are identical at some significance level (α).The null hypothesis is that two data sets are derived from identical continuous distribution.If the maximum difference between the two data sets is sufficiently large (i.e., small p-value), the null hypothesis can be rejected, and if the maximum difference is sufficiently small, the null hypothesis is accepted.This analysis compared a family of CDFs of the two TMDL allocation scenarios.In constructing such CDFs, the median values of the associated violation rate of the ensemble simulations were used (Merz and Thieken 2005).The goal of the KS test was to determine if the epistemic uncertainty of the simulated FC concentrations was greater than the aleatory uncertainty in the two alternative allocation scenarios.Moreover, a two-way ANOVA without replication tested whether the epistemic uncertainty of the simulated FC concentrations dominated the aleatory uncertainty in the two allocation scenarios.The test examines the null hypothesis that the means of two given samples are equal.The null hypothesis can be rejected if p-value > α or F > F Cr .

Results
As discussed previously, two uncertainty analyses were performed.The metric violation rate was the fraction of days when the predicted simulated daily average FC concentration was greater than the single-sample FC water quality criterion (400 cfu=100 mL) divided by the total number of days in the simulation period.

Single-Phase Monte Carlo Simulation
For the SPMCS analysis, a daily average in-stream FC concentration time series was calculated from the ensemble of daily average FC concentrations produced by the Mossy Creek watershed model.Quantile time series (10% and 90%, 80% probability interval; 2.5% and 97.5%, 95% probability interval) were also calculated using the same model output.The average and quantile time series were plotted for the two TMDL allocation scenarios.Representative plots show results for first six months of the simulation period (Fig. 4).Table 9 presents the violation rate for each time series in the prediction period.The single-sample FC water quality criterion was violated less than 1% (0.7%-0.8%) of the time during the prediction period for both TMDL allocation scenarios.However, the violation rate increased to 3.5% and 14.6% when the 80% and 95% probability intervals of the output FC concentration time series were used.
The violation rate (calculated for the ensemble average) was similar for the two TMDL allocation scenarios.However, the violation rate was different for S5 and S6 for the 80% and 95% probability intervals.The S6 scenario exhibited a greater uncertainty compared with S5, especially on the upper bound.The upper bounds of violation rate for the 95% probability interval increased as much as 10 times for S6 compared with S5.This suggests that the S5 scenario would more reliably achieve the water quality goal than would S6, which is not in agreement with the deterministic analysis by Benham et al. (2004).The inconsistency between deterministic and probabilistic frameworks when prioritizing TMDL allocations was also documented in previous research (Borsuk et al. 2002;Langseth and Brown 2011).However, the focus of decision-making in the probabilistic application was on water quality improvement and reliability; technical feasibility was not considered as a decision criterion.Despite this limitation, the implementation of either of these scenarios will be challenging because of the needed high load reductions.
Direct deposit from cattle was reduced 99% in S5, compared with 94% in S6 (Table 8), whereas loading from cropland was reduced in S5 by 90%, compared with a 95% reduction in S6.The S6 scenario exhibited greater uncertainty than S5 because it had greater input of FC direct deposit than S5 and lower input of loading from cropland.These results illustrate that the FC direct deposit is a greater source of uncertainty than cropland FC loading because the input uncertainty in the two sources was similar (log-triangular PDF, spread over two orders of magnitude).However, the SPMCS can only present the overall uncertainty and is unable to differentiate the impacts of the epistemic and aleatory uncertainty.

Two-Phase Monte Carlo Simulation
Similar to the SPMCS, the overall uncertainty of the FC concentration was estimated (Table 10).As with the previous analysis, the single-sample FC water quality criterion was violated less than 1% (0.7%-0.8%) of the time during the prediction period for both TMDL allocation scenarios.However, the violation rate increased to 2.1% and 7.8% when the 80% and 95% probability intervals were considered.These are smaller than the percentages obtained via the SPMCS, but are still similar.The TPMCS analysis also suggests that S5 is a more reliable alternative allocation scenario than S6.Again, this result contradicts the deterministic analysis by Benham et al. (2004), but the difference between the decision-making in the two studies should be kept in mind because the present study focused only on water quality improvement and reliability of the allocation a Single-sample violation rate refers to the number of days during which the water quality criterion (400 cfu day −1 ) is exceeded divided by the total number of days of simulation.
scenarios, and technical feasibility (lower reduction in wildlife direct deposit) was not considered.
In addition to the overall predictive uncertainty, which was also estimated via the SPMCS, the effects of the epistemic and aleatory uncertainty on the FC water quality criterion violations were assessed.Using the FC ensemble, the cumulative probability of number of violations for each epistemic iteration was calculated.For example, Table 11 shows the cumulative probability for the number of violations for a selected epistemic simulation (40 HSPF iterations) for the S5 allocation scenario.The cumulative probabilities of the violations were plotted as a CDF for each epistemic HSPF iteration, in which each CDF shows the probability of the number of violation incidences.The complete TPMCS analysis included 300 epistemic iterations yielding 300 CDFs.The maximum number of violation incidences for any HSPF run was 38 for the S5 scenario and 283 for the S6 scenario.Each individual CDF is the result of aleatory uncertainty.These CDFs illustrate the effect of epistemic and aleatory uncertainty on the number of violations.Visual analysis of these CDFs suggests that the water quality criterion violation CDF curves are nearly vertical for both TMDL allocation scenarios, implying a very little impact of the aleatory uncertainty (a vertical CDF implies no effect of aleatory uncertainty).It is evident from the spread of CDFs that the S6 scenario [Fig. 5(b)] exhibited a greater effect of epistemic uncertainty than did Scenario S5 [Fig.5(a)].
As mentioned in subsection "Comparative Analyses" the median of the combined epistemic simulations was used to develop the CDF for each allocation scenario (Fig. 6), and a KS test (α ¼ 0.05) was performed to compare the CDFs.The KS test  resulted in a small p-value (<0.001), suggesting that the two data sets (derived from the two allocation scenarios) are significantly different (i.e., the null hypothesis can be rejected).
An ANOVA (α ¼ 0.05) of the two scenarios showed that there is a significant effect of epistemic uncertainty on the number of violations for both allocation scenarios, but there is no such effect of aleatory uncertainty (Table 12).Therefore, overall uncertainty of the HSPF-based water quality modeling output could potentially be reduced through collection of more information about the epistemic parameters.In the present research, this result is likely an artifact of assigning only one parameter as aleatory because the focus was on framework illustration.A different categorization of parameters could affect these results, implying the impact of subjectivity in the TPMCS-based uncertainty analysis.This is, however, not a drawback of the TPMCS, but is the case for all the uncertainty analysis methods (Pappenberger and Beven 2006).
To illustrate the uncertainty in predicted in-stream FC concentrations at the watershed outlet, the average, lower (2.5% and 10%), and upper (90% and 97.5%) quantiles of the average daily FC concentration were plotted (Fig. 7).Percentage violations of the singlesample FC criterion were calculated as explained previously.The average time series of the violation rate were similar for the two TMDL allocation scenarios.However, the violation rates were different for S5 and S6 for the 80% and 95% probability intervals.Scenario S6 exhibited greater uncertainty than S5, especially with respect to the upper bound.The average time series of the violation rate was the same for the SPMCS and the TPMCS for both allocation scenarios.Because both techniques used the same parameter PDFs, agreement between the average time series was expected.The range of violation rate for the S5 scenario was similar for SPMCS and TPMCS, but for the TPCMS simulation, the range was slightly smaller for S6.

Discussion
Both SPMCS and TPMCS estimated a comparable overall uncertainty of the HSPF-predicted in-stream FC concentration for Mossy Creek (Tables 9 and 10).The same allocation scenario was also selected through both methods, suggesting that the separation of the epistemic and aleatory uncertainties does not affect the TMDL scenario prioritization in the study watershed.However, the TPMCS analysis provided a more informative picture of the uncertainties in the water quality modeling process by providing information about the effects of both epistemic and aleatory uncertainty on model predictions.In the present study, the TPMCS analysis indicated that there was a significantly greater effect of epistemic uncertainty (Fig. 5).This result is, in part, a function  of the subjective selection and classification of uncertain parameters as either epistemic or aleatory.On the other hand, the SPMCS failed to distinguish between epistemic and aleatory uncertainty.Information about the effect of epistemic uncertainty can serve as an indicator of the potential benefits of performing additional field measurements to gather more detailed information to potentially reduce epistemic uncertainty (Hession et al. 1996a).Decision makers can use this kind of information to evaluate the trade-off between gathering more information about epistemic parameters and the costs of field measurements.This raises an interesting research question about the extent to which having additional information is economical in both the TMDL development phase and the implementation phase when decisions are made about which pollutant controls will be put in place to address a given water quality impairment.The authors emphasize that efforts to answer this question are vital to making more-informed water quality-based watershed management decisions.The ultimate decision about the use of an uncertainty analysis technique depends on the project objectives.If the goal is to estimate the overall uncertainty, the SPMCS may be sufficient.
The ability to assess the uncertainty of pollutant allocation scenarios allows decision makers to make more-informed decisions about the selection of a given pollutant reduction scenario or the choice of the appropriate pollution control measures.In general, if a given TMDL allocation scenario meets water quality standards with a greater degree of confidence, it is likely to be more expensive to implement because it would require more extensive pollution control measures.Such a scenario might be preferable if regulations are strict or the ecosystem is fragile (i.e., the risk of a water quality violation has potentially greater consequences).A less expensive allocation scenario with a greater uncertainty might be preferred if the ecosystem is more robust, there are conflicts of interest, or the funding for watershed management program is limited.Additionally, an uncertainty analysis may afford decision makers the opportunity to prioritize their efforts to control the pollutant sources and choose to address those that are responsible for a greater uncertainty first.For instance, in the Mossy Creek case study presented here, although the input uncertainty in the FC loadings from cropland and direct deposit was similar, direct deposit was responsible for a greater uncertainty than runoff from croplands.
The present study was a first effort to separate the impacts of the epistemic and aleatory parameters in uncertainty analysis of watershed-scale water quality modeling.This is an early-stage demonstration paper, and the findings are not conclusive.Care should be taken in generalizing these findings, and the limitations of the study should be considered as other applications of the approach presented here are attempted.First, assignment of epistemic and aleatory parameters in the TPMCS is subjective, and some parameters can be grouped as either or both.Here, the parameters that could only be estimated via calibration were considered epistemic.This suggests that epistemic uncertainty might not only stem from the modeler's knowledge, but also arise from observed data sets and Monte Carlo random sampling (Merz and Thieken 2005).Thus, the conclusion about the impact of either of these two types of uncertainty is sensitive to how model parameters are categorized with respect to uncertainty type (Merz and Thieken 2005;Mishra 2011).This limitation is common among most uncertainty analysis techniques (Zheng and Keller 2007a).Second, this study focused only on the uncertainty that arose from the model parameters.Other uncertainties, such as those resulting from the model structure, were not taken into account.It is recommended that, where possible, other significant sources of uncertainty be incorporated into future uncertainty analyses to draw more general conclusions.Additionally, the study was conducted  under the assumption of stationarity, and the uncertainty triggered by nonstationarity, such as climate change and land use change, which might have a major impact on water quality (Fonseca et al. 2015;Whitehead et al. 2009), was not considered.Third, the TPMCS does not have a specific algorithm to assist with calibration of a probabilistic model.Hence, it provides no information about the potential for equifinality by classifying the generated parameters into behavioral and nonbehavioral sets.The utility of this technique could potentially be improved if it is used in conjunction with Bayesian methods, as proposed by Borsuk et al. (2003), Patil and Deng (2011), Hantush and Chaudhary (2014), Camacho et al. (2018), andMishra et al. (2018) for TMDLs, which can be used to estimate the posterior PDFs of the model parameters.Therefore, it is suggested that such techniques be coupled with the TPMCS to enhance its efficiency.Fourth, estimating the uncertainty of the model outputs is one step in uncertainty analysis of TMDL projects.The next step is to transform the uncertainties into the planning and decision-making phase.When prioritizing a pool of feasible pollutant allocations, uncertainty analysis should be incorporated into the analysis (Jia and Culver 2006).In addition, the uncertainties must be effectively communicated to stakeholders and decision makers.In practice, the decision makers are often poorly informed about the uncertainties of modeling (Stow et al. 2007).
Efforts should be made to advance developing efficient methods for communicating the impacts of uncertainty to decision makers (Ahmadisharaf et al. 2016;Ocampo-Duque et al. 2013;Xie and Huang 2014).

Summary
This study illustrated the application of the TPMCS uncertainty analysis framework to assess the effects of epistemic and aleatory uncertainty associated with violations of an in-stream fecal coliform water quality criteria when using the HSPF.The authors compared both SPMCS and TPMCS approaches using the Mossy Creek bacterial impairment TMDL (Benham et al. 2004) as a case study.Two FC pollutant allocation scenarios (S5 and S6) were considered.The scenarios differed in the specified FC reductions from three dominant sources in the watershed: cattle directly depositing FC in the stream, FC loadings from cropland, and wildlife contributions.Both the SPMCS and TPMCS techniques reported similar uncertainty estimates for the S5 pollutant reduction scenario, which reduced FC loading from cattle direct deposit by 99% and loading from cropland runoff and wildlife by 90% and 30%, respectively.However, the TPMCS reported a lower uncertainty than did SPMCS for the S6 scenario, which required 0% reduction in wildlife loads and reduced FC loading from cattle direct deposit by 94% and loading from cropland runoff by 95%.Both SPMCS and TPMCS suggested that the S5 scenario is a more reliable alternative, which contradicts the deterministic analysis by Benham et al. (2004).This implies that overlooking uncertainty can be misleading when developing TMDLs.This study involved an initial application of an alternative uncertainty analysis framework, the TPMCS, for TMDL development studies.Although the TPMCS framework allows separating the effects of epistemic and aleatory uncertainties in watershed-scale water quality modeling, applying the framework involves making subjective decisions about how selected model parameters are considered within the framework.The application presented here demonstrated that the TPMCS framework is transparent and the results provide information that can be used by decision makers when considering pollution control alternatives, including quantifying the level of confidence in achieving applicable water quality standards.The ultimate decision about the use of any uncertainty analysis technique depends on the project objectives, and if the goal is to estimate the overall uncertainty, the more straightforward SPMCS may be sufficient as it is less subjective and less computationally intensive.
(x and y) from the predefined PDFs for each parameter.Along with other deterministic inputs and model parameters, each generated set of epistemic and aleatory inputs represents a model iteration.The model is run for the n aleatory parameter values and the ensemble of output is plotted as a cumulative distribution function (CDF).Each CDF represents the output PDF due to the aleatory uncertainty.Next, similar CDFs are generated for the m sets of epistemic parameter random values.The resulting family of CDF curves describes both the epistemic and aleatory uncertainty.A vertical CDF illustrates no aleatory uncertainty effect, whereas overlapping CDFs indicate no epistemic uncertainty.

Fig. 1 .
Fig. 1.Two-phase Monte Carlo simulation schematic for separation of impacts of epistemic and aleatory uncertainty.(From Hession et al. 1996a.Used with permission.) Creek was monitored monthly by the Virginia Department of Environmental Quality (VADEQ) between July 1992 and March 2003 for FC concentration and other selected water quality constituents at Station 1BMSS001.35located near the outlet of the watershed.The Department of Biological Systems Engineering at Virginia Tech monitored Mossy Creek semimonthly between February 1998 and December 2001 for selected water quality constituents, including FC concentration near the VADEQ station.Daily flow data were also collected from May 1998 to December 2002 at the same site.HSPF hydrologic calibration and validation and water quality calibration were performed using measured data at this site.

Table 1 .
Potential uncertainty type classification of selected uncertain parameters in HSPF Parameter could be classified as epistemic and/or aleatory.For this study, INFLIT was the only parameter classified as aleatory.

Table 2 .
Probability density function of aleatory parameter INFILT (mm h −1 ) for different land uses

Table 3 .
Probability density function of epistemic hydrologic parameters for all land uses a Numbers in parentheses show lower and upper limits of the uniform PDF, respectively.

Table 4 .
Probability density function of epistemic hydrologic parameters that vary with land use and time for the month of January a Note: LDR = low-density residential; and HDR = high-density residential.a Numbers in parentheses show lower and upper limits of the uniform PDF, respectively.

Table 5 .
Probability density function of uncertain water quality parameters

Table 7 .
Calibrated probability density function for hydrologic parameters

Table 8 .
TMDL allocation scenarios for Mossy Creek watershed

Table 9 .
Violation rate (%) of single-sample fecal coliform criteria for the two TMDL allocation scenarios using single-phase Monte Carlo simulation based on daily time series Lower and upper quantiles of violation rate over the simulation period.

Table 11 .
Example of cumulative probability for numbers of singlesample fecal coliform criterion violations for a given epistemic simulation