Stochastic Residual-Error Analysis for Estimating Hydrologic Model Predictive Uncertainty Mohamed M. Hantush, A.M.ASCE1; and Latif Kalin, A.M.ASCE2 Abstract: A hybrid time series-nonparametric sampling approach, referred to herein as semiparametric, is presented for the estimation of model predictive uncertainty. The methodology is a two-step procedure whereby a distributed hydrologic model is first calibrated, then followed by brute force application of time series analysis with nonparametric random generation to synthesize serially correlated model residual errors. The methodology is applied to estimate uncertainties in simulated streamflows and related flow attributes upstream from the mouth of a rapidly urbanizing watershed. Two procedures for the estimation of model output uncertainty are compared: the Gaussianbased l-step forecast and the semiparametric ensemble forecast. Results show that although both methods yielded comparable uncertainty bands, the Gaussian l-step forecast underestimated the width of the uncertainty band when compared to the semiparametric method. An ensemble of streamflows generated through Latin-hypercube Monte Carlo simulations showed relatively larger values of the coefficient of variation for long-term average annual maximum daily flows than for long-term daily, monthly maximum daily, and monthly median of daily flows. Ensemble of flow duration curves is generated from the error-adjusted simulated flows. The computed low flows displayed greater values of the coefficient of variation than flows in the medium and high range. The ensemble flow durations allow for the estimation of daily flow range upstream from the outlet with 95% confidence for a specified design recurrence period. The computed uncertainties of the predicted watershed response and associated flow attributes provide the basis for communicating the risk to stakeholders and decision makers who are involved in the future development of the watershed. DOI: 10.1061/共ASCE兲1084-0699共2008兲13:7共585兲 CE Database subject headings: Watersheds; Hydrologic models; Uncertainty principles; Time series analysis; Monte Carlo method; Stochastic processes.

Introduction Extant practices in environmental modeling are primarily concerned with model calibration and validation, but often pay much less attention to the assessment of model predictive uncertainty. The widespread use of hydrologic and water quality models in water resources management and environmental decision making, however, has raised concerns about the reliability of model predictions for decision making. Models are imperfect representations of the real world and data/observations are often incomplete or insufficient and invariably subject to uncertainty. In real world problems, uncertainties are unavoidable and should be rigorously addressed in the development and application of models. The National Research Council 共NRC兲 共2001a兲 even challenged the U.S. EPA to end the practice of arbitrary selection of the margin of safety 共MOS兲 in the total maximum daily loading program 共TMDL兲 and promoted conducting uncertainty analysis as the basis for MOS determination. In a state-of-the-science review of 1

Research Hydrologist, U.S. EPA NRMRL, 26 W. Martin Luther King Dr., Cincinnati, OH 45268 共corresponding author兲. E-mail: [email protected] 2 Assistant Professor, School of Forestry and Wildlife Sciences, Auburn Univ., Auburn, AL 36849. E-mail: [email protected] Note. Discussion open until December 1, 2008. Separate discussions must be submitted for individual papers. To extend the closing date by one month, a written request must be filed with the ASCE Managing Editor. The manuscript for this paper was submitted for review and possible publication on April 11, 2007; approved on September 21, 2007. This paper is part of the Journal of Hydrologic Engineering, Vol. 13, No. 7, July 1, 2008. ©ASCE, ISSN 1084-0699/2008/7-585–596/$25.00.

the role of research in confronting the nations’ water problems, The NRC 共2004兲 called for the explicit recognition of uncertainty occurrence, measuring its importance, and incorporating it into decision making. Imperfect model structure due to incomplete knowledge of reality, insufficient data and measurement/ observation errors, spatiotemporal variability of model parameters and climate, and inadequate description of boundary conditions are all sources of uncertainty. Failure to account for these sources of uncertainty can lead to inadequate environmental decisions that can be far from satisfactory in outcome. Several research papers and reports have also emphasized uncertainty analysis as essential for hydrologic and water quality modeling 共Beck 1987; Reckhow 1994a,b; NRC 2001a,b; Singh and Woolhiser 2002; NRC 2004; Refsgaard et al. 2005; Shirmohammadi et al. 2006兲. While the scientific literature is replete with approaches and methods for estimating predictive uncertainty, careful examination of the hydrologic and environmental modeling literature reveals five dominant approaches including Bayesian uncertainty estimation, sampling-based methods, Pareto optimality, stochastic analysis of model residuals, and analysis of variance. Among these approaches, the sampling approach and Bayesian inferences are gaining increasing popularity in watershed model applications. The Bayesian approach recasts a deterministic model into a standard regression form and conducts model simulations based on Bayesian statistics to estimate uncertainties assuming zeromean and normally distributed residual errors 共e.g., Sorooshian and Dracup 1980; Bates and Cambell 2001; Vrugt et al. 2003b; Kavetski et al. 2006; Kuczera et al. 2006; Samanta et al. 2007; and Ajami et al. 2007兲. In the sampling approach, the Monte Carlo 共MC兲 method is used to obtain the distribution of model JOURNAL OF HYDROLOGIC ENGINEERING © ASCE / JULY 2008 / 585

output共s兲 by sampling model parameters from a priori probability distributions derived from literature or new knowledge gained from experience and model calibration 共Binley et al. 1991; Beven and Binley 1992; Dilks et al. 1992; van der Perk and Bierkens 1997; Schulz et al. 1999; Beven and Freer 2001; Carpenter and Georgakakos 2004; Muleta and Nicklow 2005; Hossain and Anagnostou 2005; Hantush and Kalin 2005; and Arabi et al. 2007兲. The generalized likelihood uncertainty estimation 共GLUE兲 methodology 共Binley et al. 1991; Beven and Binely 1992; Beven 1993; Freer and Beven 1996; Schulz et al. 1999; Beven and Freer 2001; Hossain et al. 2004; Hossain and Anagnostou 2005; and Arabi et al. 2007兲 is one of the most widely used Bayesian-based, sampling approaches. By pooling precision-weighted information from site-specific observed data, the application of MC simulation with Bayesian inference in GLUE could lead to reduced posterior output uncertainty band and improved model parameter distributions. Rather than seeking an “optimal parameter set,” the methodology embraces the concept of “equifinality” 共Beven 1993兲 whereby many different parameter combinations 共behavior sets兲 are considered as equally likely good simulators 共mimic the behavior兲 of the modeled system. The separation of parameter sets into behavior and nonbehavior sets in the GLUE methodology was originally introduced by Spear and Hornberger 共1980兲. The Pareto optimality approach is another approach that is inherently deterministic and multiobjective in nature 共e.g., Gupta et al. 1998; Madsen 2000; Wagener et al. 2001; Vrugt et al. 2003a兲; it shares a similar notion with the “equifinality” concept of GLUE methodology in the sense that the number of the 共good or behavior兲 parameter sets, namely, Pareto sets, defines model output uncertainty. A fourth approach involves direct stochastic analysis of model residual errors 共e.g., Kim and Valdés 2005; Kim et al. 2006; and Bruen and Yang 2006兲. In this approach, observed or residual-error time series are described by a stochastic model. The fifth strategy for quantifying uncertainty is the well known first-/ second-order analysis where the mean and variance of model output are expressed in terms of means and variances of the input random variables 共Zhang and Yu 2004; Melching and Bauwens 2001; Haan 2002; and Ang and Tang 2007兲. While the first-/ second-order approximation has a particular appeal to a wide range of engineering applications, it is limited to small parameters’ variance and requires significant coding in complex watershed models. Among these five strategies, stochastic modeling of residual errors has been given the least attention. Although this approach requires traditional manual/automated model calibration as a first step, it nonetheless applies to serially correlated nonGaussian residual errors and accounts for both parametric and model structure uncertainties. The main premise of the methodology is that residual errors of a calibrated model are serially correlated and can be described by a quasi-parametric stochastic model. This paper applies a hybrid time series-nonparametric probabilistic sampling approach to residual errors of a calibrated watershed model 共SWAT兲 to estimate its predictive uncertainty. The SWAT model is calibrated and validated for daily flows in the Pocono Creek watershed which is located in eastern Pennsylvania and is experiencing rapid population growth and urbanization 共Kalin and Hantush 2006a,b兲. The objectives of this paper are to: 共1兲 demonstrate a novel, yet simple approach for quantifying model predictive uncertainty; 共2兲 compare two procedures for the estimation of uncertainty band, namely, the Gaussian-based l-step forecast and the semiparametric method; 共3兲 evaluate model forecast performance; and 共4兲 compute ensemble of duration curves 586 / JOURNAL OF HYDROLOGIC ENGINEERING © ASCE / JULY 2008

and estimate long-term statistics of flow attributes that have management and ecological significance. The remainder of the paper is organized in sections and subsections. It devotes an entire section to the development of a time series model for model residual errors, construction of the l-step forecast, and description of the semiparametric approach. The subsequent section deals with the application of the methodology to the calibrated SWAT model for the Pocono Creek watershed, followed by discussion of the results. The paper ends with the conclusions.

Methodology Prediction Error We start by noting that while model uncertainty is actually the deviation of a calculated value from the true value, the analysis that follows rather focuses on the deviation of calculated value from the observation. Model error, ␧t, is thus defined by this equation ␧ t = o t − q t,

t = 1,2, . . . T

共1兲

where ot = measured streamflow; qt = model computed streamflow; t = time index 共days兲; and T = length of the time record. The ␧t accounts for model structural uncertainty, parametric uncertainty, measurement error, and errors in the precipitation and other model inputs 共e.g., surface elevation, climatic data, channel morphology兲. It should be noted that an attempt to predict ␧t will yield a forecast for streamflow excluding the effect of measurement errors. The smaller the measurement error, the closer the forecast is to actual flow values. Fig. 1共a兲 shows the time series of residual errors ␧t of a SWAT hydrologic model calibrated and validated by Kalin and Hantush 共2006b兲 for the Pocono Creek watershed for the period July 1, 2002–April 30, 2005. The 120 km2 Pocono Creek watershed is located between the latitudes 40° 59⬘N – 41° 06⬘N and longitudes 75° 14⬘W – 75° 26⬘W in Monroe County, Pa., situated within the Delaware River Basin. The watershed valley is 26 km long and drains from the Pocono Plateau in its headwaters eventually into the Brodhead Creek which is a tributary to the Delaware River. The model is constructed for the area above a United States Geological Survey 共USGS兲 stream gauge located about 6.4 km upstream from the mouth near the city of Stroudsburg, Pa., and drains an area of about 98.87 km2. The model was calibrated and validated using a split data set approach and the Nash–Sutcliffe efficiency for the simulated outlet daily streamflows used in this analysis was 0.74 for the calibration and 0.71 for the validation runs. The coefficient of determination 共R2兲 was 0.74 and 0.72 for the calibration and validation runs, respectively 共Kalin and Hantush 2006b兲. These performance evaluation values are well above the threshold values recommended by Santhi et al. 共2001兲 for SWAT calibration. However, both performance evaluation measures were relatively low for the simulated base flow, but with less than 5% mass balance error at the monthly time scale. A close look at the daily ␧t time series reveals systematic errors shown as large positive or negative spikes 关Fig. 1共a兲兴 that are symptomatic of SWAT performance at the daily time scale and appears to be not random in nature. In many instances a large 共⫺兲 error is immediately followed by a large 共⫹兲 error. One reason might be timing errors in recording either the precipitation or the streamflow. Another, probably more reasonable, explanation is the inability of the SWAT model to use subdaily rainfall data. When a

Fig. 1. Time series of model prediction errors 共observed—SWAT simulated兲 for calibration and validation time period 共July 1, 2002–April 30, 2005兲: 共a兲 daily simulations; 共b兲 3-day averages

large rainfall event happens at the very end of a day, the actual watershed response will be about 1 day delayed compared to what SWAT actually simulates with daily rainfall data. This will clearly result in overprediction 共negative error兲 and underprediction 共positive error兲 on 2 successive days, respectively. Hence, to minimize such types of errors to some extent, we decided to use the average over 3 consecutive days 共i.e., 3-day moving average兲 as a surrogate to daily flows. Two compelling reasons drove the need for consideration of simulations at the daily scale. First, monthly median daily flow is a required flow metric 共or attribute兲 for wild trout habitat condition, as one of the primary objectives of the calibrated watershed model is to assess the impact of projected future buildout in the watershed on the sustainability of the wild brown trout population in the Pocono Creek. Second, there is not sufficiently long monthly streamflow data to conduct a meaningful time series analysis at the monthly time scale. Therefore, in addition to smoothing the effect of successive large 共⫺兲 and 共⫹兲 errors during some significant events, 3-day averaged streamflows provide surrogates to daily flows for use in the wild trout habitat assessment analysis. Later in the analysis, it will be shown through regression relationships derived from model simulated daily flows that estimates of 3-day averaged flow attributes can be converted back to the daily ones with high efficiency 共R2 ⬎ 0.97兲. Time Series ARIMA Model Development Fig. 1共b兲 shows the 3-day average of model residual errors for the calibration and validation period 共July 1, 2002–April 30, 2005兲. The errors are not totally eliminated but their magnitudes are significantly reduced. The Minitab statistical computer package is used to identify a time series model and construct a reasonable forecast for future ␧t values. In the most general form, the sequence ␧t may be described by a multiplicative seasonal autoregressive integrated moving average model denoted by ARIMA

共p , d , q兲 ⫻ 共P , D , Q兲S. The parameters p, d, and q are, respectively, the orders of the autoregressive, difference, and moving average operators; P, D, and Q are, respectively, the orders of the seasonal autoregressive, seasonal difference, and seasonal moving average operators; and S is the seasonal period. Autocorrelation function 共ACF兲 and partial autocorrelation function 共PACF兲 are utilized and the procedure outlined by Shumway 共1988兲, which derives its basis from Box and Jenkins 共1970兲, is followed to identify the best model. The objective of the model identification process is to produce identically distributed, independent residuals wt, with a minimum variance ␴2w arising from fitting an ARIMA 共p , d , q兲 ⫻ 共P , D , Q兲S to the ␧t time series. Figure 2共a兲 shows the computed the computed ACF and PACF of ␧t. It is obvious that the series display nonstationarity accentuated by the slowly decaying ACF as a function of lag and the large positive value of PACF at Lag 1. The ACF of the first difference in Fig. 2共b兲 contains a fairly strong negative peak 共−0.44兲 at Lag 3 and zero thereafter, indicating that a seasonal moving average with S = 3 共days兲 and Q = 1 might be appropriate. The decreasing peaks of the PACF at multiples of three are due to the seasonal moving average component. Therefore, ARIMA 共0 , 1 , 0兲 ⫻ 共0 , 0 , 1兲3 appears to be, at the moment, a suitable selection among, perhaps, other competing models. We experimented by adding autoregressive and moving average components progressively. Table 1 lists some of the examined ARIMA models, with computed variance of residuals, ␴ˆ 2w, and various goodnessof-fit measures: final prediction error 共FPE兲, Akaike’s information criterion 共AIC兲, and Bayesian information criterion 共BIC兲. It is clear that the most appropriate model choice may be ARIMA 共1 , 1 , 1兲 ⫻ 共0 , 0 , 1兲3, which has the smallest values of ␴ˆ 2w, FPE, AIC, and BIC. The model was fitted to the error time series depicted in Fig. 1共b兲 共recall the 3-day averages of the actual model errors兲, and the equation describing ARIMA 共1 , 1 , 1兲 ⫻ 共0 , 0 , 1兲3 is JOURNAL OF HYDROLOGIC ENGINEERING © ASCE / JULY 2008 / 587

Fig. 2. Autocorrelation 共ACF兲 and partial autocorrelation 共PACF兲 functions of: 共a兲 prediction error ␧t; 共b兲 first difference of residual error ⵜ␧t = ␧t − ␧t−1; and 共c兲 residuals 共wt兲 for ARIMA 共1 , 1 , 1兲 ⫻ 共0 , 0 , 1兲3. Lag is in units of days.

共1 − ␾B兲 ⵜ ␧t = 共1 − ⌰B3兲共1 − ␪B兲wt

共2兲

where ⵜ␧t = ␧t − ␧t−1 = first difference operator; B␧t = ␧t−1⫽backward shift operator; and ␾, ␪, and ⌰ = fitted parameters. Equation 共2兲 can be expanded to yield the following representation for ␧t ␧t = 共1 + ␾兲␧t−1 − ␾␧t−2 + wt − ␪wt−1 − ⌰wt−3 + ␪⌰wt−4

共3兲

␾, ␪, and ⌰ are estimated by maximizing the log likelihood function for wt, t = 1 , . . . , T, with maximum likelihood estimates of ˆ = 0.98, respectively. The maximum likeˆ = 0.78, ␪ˆ = 0.62, and ⌰ ␾ lihood estimate for the variance ␴2w is ␴ˆ 2w = 0.629. The ACF and PACF of wt from this model is plotted in Fig. 2共c兲, and shows no prominent peaks at 5% significance limits. Therefore, it seems reasonable to regard the residuals as being white noise. The Box–Pierce Q for Lags 48, 100, and 250 measuring randomness are satisfied with 95% confidence for the selected ARIMA 共1 , 1 , 1兲 ⫻ 共0 , 0 , 1兲3 model. In the following sections, we present two methods for the estimation of uncertainty limits of the simulated streamflows. The 588 / JOURNAL OF HYDROLOGIC ENGINEERING © ASCE / JULY 2008

first is referred to as the l-step forecast and the second as the semiparametric method. l-Step Forecast The l-step forecast is given by the conditional mean 共Shumway 1988兲 t ␧t+l = E关␧t+l兩␧t,␧t−1, . . . 兴

共4兲

and the forecast variance l−1

t t 2 = E关共␧t+l − ␧t+l 兲 兩␧t,␧t−1, . . . 兴 = ␴2w Pt+l

⌿2k 兺 k=0

共5兲

where ⌿k is to be determined later by noting that ␧t can be expressed in terms of wt

Table 1. Values of Residual Variance 共␴w2 兲, Final Prediction Error 共FPE兲, Akaike’s Information Criterion 共AIC兲, and Bayesian Information Criterion 共BIC兲 for Various Models Applied to Model Prediction Error Data, ␧t ARIMA 共p , d , q兲 ⫻ 共P , D , Q兲s model 共0 , 1 , 0兲 ⫻ 共1 , 0 , 0兲3 共0 , 1 , 0兲 ⫻ 共2 , 0 , 0兲3 共0 , 1 , 0兲 ⫻ 共0 , 0 , 1兲3 共0 , 1 , 1兲 ⫻ 共0 , 0 , 1兲3 共0 , 1 , 2兲 ⫻ 共0 , 0 , 1兲3 共1 , 1 , 0兲 ⫻ 共0 , 0 , 1兲3 共1 , 1 , 1兲 ⫻ 共0 , 0 , 1兲3

␴ˆ w2

FPE

AIC

BIC

0.926 0.805 0.674 0.653 0.641 0.647 0.629

0.928 0.808 0.675 0.656 0.645 0.650 0.633

−0.075 −0.213 −0.393 −0.422 −0.439 −0.432 −0.458

−0.070 −0.203 −0.388 −0.413 −0.425 −0.422 −0.443

⌿kwt−k, 兺 k=0

⌿0 = 1

⬁

␧t = ⌿共B兲wt = 共1 + ⌿1B + ⌿2B2 + . . . 兲wt =

共6兲 where B again= backward shift operator. We may construct the forecast using the following expectation identities E共␧t兩␧s,␧s−1, . . . 兲 =

再 再

␧t , t 艋 s ␧st , t ⬎ s

wt , t 艋 s E共wt兩␧s,␧s−1, . . . 兲 = 0, t ⬎ s

冎 冎

共7兲

共8兲

Applying Eqs. 共7兲 and 共8兲 to Eq. 共3兲 yields the following l-step forecasts t = 共1 + ␾兲␧t − ␾␧t−1 − ␪wt − ⌰wt−2 + ␪⌰wt−3 ␧t+1

共9兲

t t = 共1 + ␾兲␧t+1 − ␾␧t − ⌰wt−1 + ␪⌰wt−2 ␧t+2

共10兲

t t t = 共1 + ␾兲␧t+2 − ␾␧t+1 − ⌰wt + ␪⌰wt−1 ␧t+3

共11兲

t t t ␧t+4 = 共1 + ␾兲␧t+3 − ␾␧t+2 + ␪⌰wt

共12兲

t t t = 共1 + ␾兲␧t+l−1 − ␾␧t+l−2 , ␧t+l

l艌5

共13兲

Note that the one-step forecast Eq. 共9兲 has the terms wt, wt−2, wt−3 which are not directly observed in the data sample ␧1 , ␧2 , . . . , ␧T. These values can be approximated by assuming that w2 = w1 = w0 = w−1 = 0 and rewriting Eq. 共3兲 wk = ␪wk−1 + ⌰wk−3 − ␪⌰wk−4 + ␧k − 共1 + ␾兲␧k−1 + ␾␧k−2 共14兲 for k = 3 , 4 , . . . , t. The coefficients ⌿1, ⌿2 , . . . in Eq. 共6兲 can be found explicitly by substituting Eq. 共6兲 into Eq. 共2兲 and equating coefficients of Bk, k = 1 , 2 , 3 , . . . in the relation 共1 − ␾B兲共1 − B兲共1 + ⌿1B + ⌿2B2 + . . . 兲 = 共1 − ⌰B3兲共1 − ␪B兲 共15兲 This leads to ⌿1 = 1 + ␾ − ␪

共16兲

⌿2 = 共1 + ␾兲⌿1 − ␾

共17兲

⌿3 = 共1 + ␾兲⌿2 − ␾⌿1 − ⌰

共18兲

⌿4 = 共1 + ␾兲⌿3 − ␾⌿2 + ␪⌰ ⌿k = 共1 + ␾兲⌿k−1 − ␾⌿k−2,

k = 5,6, . . .

共19兲 共20兲

The recursive relationships Eqs. 共16兲–共20兲 can be substituted into Eq. 共5兲 to compute the l-step forecast variance as the mean-square error of prediction. For normally distributed wt, the 共1-␣兲 probability interval 共confidence band兲 for the residual forecast value is t t ␧t+l ⫾ z␣/2冑Pt+l ,

l = 1,2, . . .

共21兲

where z␣/2 = upper ␣ / 2 tail on the standard normal distribution. Although not shown here, plotting positions reveal that the ˆ t are non-Gaussian. While this violates the computed residuals w Box–Pierce goodness of fit, which is based on normally distributed wt, the ACF and PACF in Fig. 2, ␴ˆ 2w, and comparison of FPE, AIC, and BIC criteria altogether point toward a suitable model fit. For non-Gaussian wt, the l-step forecast Eqs. 共9兲–共13兲 and its variance Eq. 共5兲 are no longer sufficient to infer the 共1-␣兲 confidence band; however, a nonparametric method can be employed to construct an empirical probability density function for wt. Eq. 共21兲, nonetheless, still can be used to provide a crude estimate for the confidence band. Semiparametric Method For wt with arbitrary PDF, we use the parametric relationship Eq. 共3兲 as a base for constructing the SWAT prediction uncertainty band. When commonly used parametric probability distributions 共e.g., normal, log-normal, exponential, etc.兲 poorly fit the frequency of the observed stochastic series, the probability density function of the observed random variable can be locally approximated by a nonparametric model 共Lall 1995兲. We therefore explore a nonparametric model fit to the observed probability density functional 共PDF兲 of wt. In parametric methods, the density function is estimated by assuming that data are drawn from a known parametric family of distributions. The methods of moments, maximum likelihood estimation, or any other methods are commonly used to estimate the parameters of the chosen PDF. In nonparametric approaches, a kernel function is often used to generalize the density function estimation. Given a set of n observations w1, w2 , . . . , wn, a mathematical expression of a univariate kernel probability density estimator is 共e.g., Kim and Valdés 2005兲 n

冉 冊

x − wk ˆf 共x兲 = 1 K X nh k=1 h

兺

共22兲

where x = given value of the random variable X denoting wt; K = kernel function; n = number of observations; and h⫽bandwidth that controls the variance of the kernel function. Kim and Valdés 共2005兲 provide a list of kernel functions typically used in hydrology; the most widely used one being the Gaussian K共x兲 =

1

冑2␲ e

−共x2/2兲

共23兲

An optimal estimate for the bandwidth for a Gaussian kernel is provided by Kim and Valdés 共2005兲 citing Silverman 共1986兲 hˆ =

冉 冊 4 3n

1/5

␴

共24兲

where ␴ = standard deviation of the observed record, which in this case is equal to the standard deviation of observed residuals, JOURNAL OF HYDROLOGIC ENGINEERING © ASCE / JULY 2008 / 589

Application Model Forecast Evaluation To examine model forecast performance, the period from July 1, 2002 to April 30, 2005 共1,035 days兲 is divided into two parts. The first 900 days 共July 1, 2002–December 16, 2004兲 were used to construct the residual error time series model and the remaining 135 days from 共December 16, 2004–April 30, 2005兲 were used for validation. The 共1-␣兲 confidence l-step forecast is obtained by combining the SWAT simulation and ␧t forecast t t ⫾ z␣/2冑Pt+l , qˆt+l = qt+l + ␧t+l

Fig. 3. Frequency histogram 共%兲, Gaussian PDF, N共x兲, and nonparametric model fit, f共x兲, to computed series of residuals wk. Both histogram and Gaussian distribution have zero mean and standard deviation, ␴ˆ 2w

␴ˆ w ⬇ 0.79. n = 1,028= approximately the number of streamflow measurements in the calibration and validation period July 1, 2002–April 30, 2005兲. Although the residuals wt are zero mean and statistically independent, they do not follow any of the known parametric distributions. Fig. 3 shows that Gaussian PDF⬃ N共0 , ␴ˆ 2w兲 inadequately fits the observed histogram of wt and the nonparametric fit 关Eqs. 共22兲–共24兲兴 provides a remarkably improved local fit with a longer tail. However, one should bear in mind that a nonparametric approach is limited by the hypothesis that future wt has a similar nonparametric functional form of the fitted PDF. Both distributions appear to be symmetric with almost zero mean and medians. Random sampling of wt is achieved in two steps. First, the cumulative distribution function 共CDF兲, FX共x兲 is constructed by integrating Eq. 共22兲 共with u replacing x in the integrand兲 from u = −⬁ to u = x to yield n

冉 冊

1 1 x − wk erf FX共x兲 = + 2 2n k=1 冑2hˆ

兺

共25兲

where erf共u兲 = 共2 / 冑␲兲兰u0e−y dy = well known error function. Second, a procedure for sampling x from ˆf X共x兲 starts with generating a random number R from the uniform distribution in the interval 关0, 1兴 共e.g., Haan 2002兲, then setting R = FX共x兲 and solving for x = F−1 X 共R兲 using any of the finding searching techniques. A large sequence of independent, identically distributed model residuals wt are generated randomly using the foregoing procedure and subsequently fed into Eq. 共3兲 to synthesize an ensemble of ␧t time series. Latin-hypercube sampling 共LHS兲 is employed as an efficient and effective alternative to conventional Monte Carlo sampling 共MCS兲. The LHS 共McKay et al. 1979兲 divides the CDF 共i.e., the FX共x兲 function兲 into segments of equal width from each of which a random variate R is generated. This way the whole CDF is covered, but with a smaller number of replications than the MCS. LHS is more efficient than MCS at both ends of the CDF. Finally, the streamflow forecast is obtained by adding each of the simulated daily streamflow rates. 2

590 / JOURNAL OF HYDROLOGIC ENGINEERING © ASCE / JULY 2008

l = 1,2, . . .

共26兲

For the semiparametric forecast, the procedure outlined in the previous section was implemented to synthesize 500 sequences of ␧t values which then were added to the corresponding SWAT model simulation conducted over the same time period to obtain the ensemble forecast 共500 time series of streamflows兲 qˆt = qt + ␧ˆ t

共27兲

qˆt = estimate of 3-day average of the streamflow; and ␧ˆ t = sampled sequence of ␧t. For each t 共i.e., day兲, the median and uncertainty band are obtained from the corresponding ensemble of 500 values of qˆt. As pointed out earlier, the uncertainty band computed with Eqs. 共26兲 and 共27兲 actually provides a degree of confidence in the computed value relative to the observation rather than the true value. Figs. 4 and 5 compare observed streamflow during the validation period with the 95% confidence band and the median of qˆt computed by the two methods. The l-step forecast with an assumed Gaussian white noise 共wt兲 produced a relatively narrower confidence band 共Fig. 4兲 than that constructed based on the semiparametric model, with nearly twice the number of observations falling outside the 95% confidence band for the former method. About 7% of the observed values fall outside the 95% confidence band for the semiparametric method. Although not shown at the depicted scale, two large values of the observed streamflows fell way outside the confidence band during the time period March 26–April 15, 2005 for both methods. However, the bulk of the observed streamflows fell within the uncertainty band for both methods. The median of the l-step forecast and the semiparametric method compared fairly well with the measurements. To further examine the model forecast during low-flow conditions, the forecast evaluation period is extended from May 1, 2005 to September 30, 2005. This period of data was not available during the initial model calibration/validation efforts by Kalin and Hantush 共2006b兲. During this period, no significant storm events were recorded. Although the simulated median daily flows generally compared well with measured counterparts further validating the calibrated model, the 95% confidence bands, however, were wide for both the l-step forecast and nonparametric ensemble forecast 共Fig. 6兲; the former still underestimates the uncertainty bandwidth. Ensemble Streamflow Attributes In this section and the following one, we estimate distributional characteristics of key streamflow attributes that have management and ecological significance, and construct ensemble streamflow duration curves. We use SWAT and the stochastic residual-error model as a base to generate an ensemble of 500 20-year long time series of daily streamflows 共actually 3-day averages兲 through MC simulations. Table 2 lists the mean, standard deviation 共SD兲, co-

Fig. 4. l-step forecast 共assuming normality of wt兲 of 3-day averaged streamflows and measured counterparts during validation period 共December 16, 2004–April 30, 2005兲. Total number of observed flows that are outside 95% confidence band is 18 共13.3%兲.

efficient of variation 共COV兲, median, and 2.5th and 97.5th percentiles 共95% confidence interval兲 for average daily flow, average monthly median of daily flows, average monthly maximum daily flow, and average annual maximum daily flow, each of which is computed from the ensemble of 500 time series of the synthesized daily streamflows. In obtaining the ensemble of 20-year streamflows, rather than using historical data and in order to account for uncertainty in future precipitation values, we used the SWAT built-in weather generator WXGEN to generate 500 sets of 20-year long records of daily precipitation at the two gauge stations that were used in the model calibration 共see Kalin and Hantush 2006a兲. Each of the synthesized precipitation sequences is fed into SWAT to obtain a 20-year long time series of daily streamflows. A 30 year warmup period is used to eliminate the effect of unknown initial conditions. For each MC simulation, 3-day averages are computed from the simulated daily streamflows. Similarly, an ensemble of 500 20-year long sequences of ␧t is generated as outlined above, and each sequence 共realization兲 is superimposed on each SWAT MC simulation to obtain an ensemble forecast of streamflows at the USGS stream gauge, 6.4 km upstream from the mouth of the watershed. In Table 2 summary statistics related to high-flow conditions as well as long term averages are given. This method did not perform well for low-flow index 共e.g., 7Q10兲, because generated

large negative errors superimposed on streamflows resulted in negative flow values 共about 12% of all the synthesized flows were negative兲. Therefore, low-flow indices were not computed. The term “average” in each column title denotes the arithmetic average over the 20-year simulation period. For example, for each realization 共i.e., time series兲 out of 500, daily flows are averaged over the 20-year simulation record to yield “average daily flow.” Thus, there are 500 such “average daily flow” values from which the mean, median, COV, and 95% confidence limits were computed. The average monthly median of simulated daily streamflow in column two represents the arithmetic average over the simulation period of the monthly median values. As indicated earlier, this is a required metric to a fish habitat model used by the Pennsylvania Fish and Boat Commission. In column three the monthly maximum daily streamflows averaged over the simulation period are given. Column four shows the annual maximum daily streamflows averaged over the simulation period. The average daily flow, average monthly median of daily flows, and average monthly maximum daily flow have small COV values and relatively narrower 95% confidence bands. The average annual maximum daily flow, on the other hand, shows a relatively greater uncertainty, as reflected by the larger COV and wider 95% confidence band. Note that the mean and median of the four streamflow characteristics are almost equal, indicating symmetric distributions. Compared to absolute predictions of daily

Fig. 5. Semiparametric ensemble forecast 共median and 95% confidence band兲 and measured streamflows during validation period 共December 16, 2004–April 30, 2005兲. Total of nine measurements lie outside 95% confidence band 共6.7%兲. Simulated and measured values are 3-day averages. JOURNAL OF HYDROLOGIC ENGINEERING © ASCE / JULY 2008 / 591

Fig. 6. SWAT model forecast 共median and 95% confidence band兲 and measured streamflows during postvalidation period 共May 1, 2005– September 30, 2005兲: 共a兲 l-step forecast; 共b兲 semiparametric ensemble forecast. Simulated and measured values are 3-day averages.

streamflows—which can be highly uncertain as shown in Fig. 6—these flow attributes appear to be fairly predictable and therefore more reliable and suitable for design purposes. Ensemble Flow-Duration Curves A flow duration curve 共FDC兲 is constructed for each time series qˆt in the ensemble forecast. Fig. 7共a兲 plots the median of the computed FDCs for the 3-day averaged streamflow rates and the 95% confidence interval. A FDC is a plot of the flow rate 共say, x兲 versus probability of exceedance, P共X 艌 x兲. The return period, Tr, defined as the average recurrence interval between events produc-

ing flow rates equal to or greater than a specified flow magnitude, x, is the reciprocal of P共X 艌 x兲: Tr = 1 / P共X 艌 x兲. High flows typically are associated with small probability of exceedance, whereas low flows are associated with large probability of exceedance or, equivalently, with small return period. The ensemble of FDCs can be interpreted as follows. As an example, from Fig. 7共a兲 it can be seen that for P共X 艌 x兲 = 0.002, daily streamflow at the gauge station is between 16 and 25 m3 / s with 95% confidence. Since a return period Tr = 1 / 0.002= 500 days corresponds to 0.002 probability of exceedance, then with 95% confidence, the Tr = 500 day is between 16 and 25 m3 / s. Similarly, from the fig-

Table 2. Summary of Ensemble Statistical Parameters of Selected Streamflow Characteristics Obtained from 500 Time Series of Synthesized 3-Day Average Streamflows; Values within Parentheses Denote Daily Statistics Converted from 3-Day Statistics Using Regression Relationships Shown in Fig. 9 Average daily flow 共m3 / s兲

Average monthly median daily flow 共m3 / s兲

Average monthly maximum daily flow 共m3 / s兲

Average annual maximum daily flow 共m3 / s兲

Mean

2.42

SD COV Median

0.203 0.084 2.42

2.10 共2.06兲 0.182 0.087 2.10 共2.06兲 1.75–2.45 共1.72–2.41兲

7.12 共8.06兲 0.411 0.058 7.13 共8.06兲 6.35–7.91 共7.12–8.99兲

17.39 共20.27兲 1.93 0.111 17.43 共20.32兲 13.91–21.82 共16.55–24.84兲

95% C.I.

2.02–2.80

592 / JOURNAL OF HYDROLOGIC ENGINEERING © ASCE / JULY 2008

Fig. 8. Coefficient of variation computed from ensemble of flow duration curves versus probability of exceedance P共X 艌 x兲

Fig. 7. Median and 95% confidence band for MC simulated duration curve: 共a兲 3-day averaged flows; 共b兲 daily flows. Flow duration curves are generated from 500 replicas of 20 years daily streamflows, oˆt.

ure, a Tr = 2 yr return flow 共corresponding to exceedance probability= 0.14%兲 can be expected to be between 17 and 28 m3 / s at 95% confidence. The increasingly wider 95% confidence band with decreasing probability of exceedance can be visually illusive as to the degree of uncertainty. Close inspection of Fig. 8, however, reveals greater COV values for low flows and, thus, shows much higher uncertainty than medium range and high flows. It is interesting to note that computed streamflow rates that vary from about 2 to 11 共m3 / s兲 in the median CDF—corresponding to return periods of 2.5 and 50 days, respectively—have the smallest COV values and, therefore, the least uncertainty. Higher flows with a recurrence period longer than 50 days have relatively higher uncertainty but much smaller than that associated with low flows, with return periods smaller than 2 days. Also note that the COV of 1 year return period flow compares well with the COV of the average annual maximum flow listed in Table 2. Of course, these results are specific to Pocono Creek watershed.

Discussion The presented analysis was preceded by an attempt to fit the ARIMA model to the daily residual errors. The shear number of spikes in the residual errors, which are inherent to the SWAT model during storm events, resulted in numerous synthesized

negative streamflow values, thus precluded any meaningful uncertainty analysis. We attempted to transform the model errors, ␧t, yet the coexistence of negative and positive values prevented the commonly used log 共i.e., ln ␧t兲 and Box–Cox transformations 兵关共␧t + 1兲␭ − 1兴 / ␭其. We also tried transforming observed and model simulated flows, and defined the model error as the difference between the transformed quantities. The log transformation yielded magnified synthesized flows and the Box–Cox transformation was limited by the condition of having a positive base in the power term. Thus, analyses are carried out without transformation. However, the smoothing produced by the 3-day moving average resulted in reasonably synthesized surrogates to the simulated daily streamflows, as bulk of the observed 3-day averaged streamflows fell well within the constructed 95% confidence bands for the validation and postvalidation periods. This is notwithstanding the relatively wide confidence band obtained for the postvalidation period 共Fig. 6兲 during which no significant rainfall events have occurred. The 3-day average flows used in the analyses in fact could be considered as transformation of daily flows. Fig. 9 provides linear regression relationships, which are obtained from model simulation results, to convert streamflow attributes of 3-day averages back to daily ones. Values shown in parentheses in Table 2 are daily flow statistics obtained through regression relationships shown in Fig. 9. Similar regression relationships could be obtained for other flow attributes and their statistical properties. For instance, the FDC for the daily flow could be obtained from the 3-day average FDC as follows: 1. Divide the FDC into segments by using the exceedance probabilities as the breakpoints; 2. From the model simulations compute the flows 共Q兲 corresponding to each exceedance probability 共p兲, Q p, for daily and 3-day average flows. Repeat this for each MC simulation result; 3. Regress Q p values obtained from daily and 3-day average flows for each break point 共p兲; and 4. Use the regression equations obtained in step 共3兲 as a base to convert the synthesized 3-day average flows to daily flows. By following the above procedure we converted the ensemble of 3-day average FDC shown in Fig. 7共a兲 to the ensemble of daily FDC 关Fig. 7共b兲兴. The exceedance probabilities of 0.001, 0.01, 0.05, 0.1, 0.2, . . . , 0.8, 0.9, 0.95, 0.99, 0.999 are used as the breakpoints, and regression equations are generated at these locations with R2 艌 0.93. Of course, the ideal procedure is to construct such regression relationships from measured streamflows, provided that adequately long historical records are available. Unfortunately, this is not the case here. Comparison of the two JOURNAL OF HYDROLOGIC ENGINEERING © ASCE / JULY 2008 / 593

Fig. 9. Conversion of 3-day averages to daily values of flow attributes

figures reveals that daily flows with probability of exceedance 0.01 and less 共i.e., return period 艌100 days兲 tend to be relatively higher than the corresponding 3-day averages, whereas no significant differences are observed among the two FDCs for probability of exceedance greater than 0.01. It is evident from Fig. 7 and Table 2 that the 3-day averaging affects high flows more so than the more persisting medium to low flows. Because of their relative simplicity, classical linear time series models 共e.g., Shumway 1988; and Box and Jenkins 1970兲 have been widely used in hydrology, economics, and many other fields. However, the often cited requirement of normally distributed white noise, wt, has precluded the use of times series models, especially in the class of hydrologic problems characterized by non-Gaussian random terms. This may be partly perpetrated by the inadequacy of forecasting/synthesis algorithms in readily 594 / JOURNAL OF HYDROLOGIC ENGINEERING © ASCE / JULY 2008

available statistical packages to handle non-Gaussian wt. This limitation is relaxed in this paper by fitting a nonparametric distribution to the random component of the residual error model 关wt in Eq. 共3兲兴 and inverting the resulting CDF 关Eq. 共25兲兴 to generate independent wt values sampled from generic distributions. Excluding synthesized negative streamflow values, the computed 95% confidence bands, nonetheless, were large enough to encompass most of the observed values. The pitfall in the blind application of the l-step forecast without challenging the normality assumption can lead to inadequate estimation of uncertainty. Although, not strongly accentuated by this application, underestimation of uncertainty band can lead to a false sense of model reliability, as depicted by the narrower uncertainty band for the l-step forecast 共Fig. 4兲. Figs. 5 and 6 show that for the semiparametric method most of the observed streamflows fell well within their respective 95% confidence intervals. The higher COV values for low flows 共Fig. 8兲 point toward significant uncertainties in the SWAT simulated low flows. This corroborates the higher COV values obtained by Kalin and Hantush 共2006a兲 for the low-flow attributes, base flow, and 7Q10 relative to the COVs of median and high-flow metrics. The relatively high uncertainty was also observed for small rainfall events in a different model study in which the physically based runoff and sediment erosion model KINEROS yielded higher COV values during small rainfall events 共Hantush and Kalin 2005兲. Since baseflow dominates streamflow during low-flow periods, the relatively high uncertainty in the simulated low flows might be attributed to inadequate calibration of the groundwater component of the SWAT model and/or errors arising from the baseflow separation method. Different methods of baseflow separation can give significantly different baseflow estimates; this clearly affects both calibrated model parameters and model performance. It should be noted that baseflow separation is the first step recommended in SWAT model calibration, and the less stringent ⫾15% of baseflow volume error requirement recommended by Santhi et al. 共2001兲 with no threshold values for the Nash–Sutcliffe efficiency and coefficient of determination highlight the dilemma in attempting to calibrate a quantity whose measurement itself is uncertain. Although not so evident, the proposed methodology may also be inadequate to handle low flows. In general, the relatively low COV values 共Table 2兲 indicate that the calibrated SWAT model is capable of predicting the streamflow attributes fairly reliably, especially when compared to predictions of real-time flow values. However, the relatively large COV of model forecasted annual maximum daily flow indicates a level of uncertainty to be communicated in terms of a greater risk factor for the performance flood control measures designed based on this metric, whereas the relatively low uncertainty associated with median monthly daily flow indicates that the calibrated model reliably estimates this ecologically significant indicator. The average daily flow, which also measures base flow during interstorm periods, showed an almost similar COV value to that for the average monthly median of daily flows.

Conclusions A two-step methodology was presented for the estimation of predictive uncertainty of a distributed hydrologic model calibrated for the rapidly urbanizing Pocono Creek watershed. The first step requires adequate model calibration and the second, which is the topic of this paper, deals with stochastic residual-error analysis. The aim of the stochastic error analysis was the estimation of

uncertainty associated with daily streamflows and related longterm attributes simulated by the distributed hydrologic model. The methodology combined classical linear time series analysis 共ARIMA兲 with a nonparametric probabilistic sampling methodology to synthesize an ensemble of serially correlated model residual errors of simulated streamflows by means of Monte Carlo-type Latin-hypercube simulations. The ensemble residualerror time series was used to construct the forecast and long-term daily flow attributes of management and ecological significance, and flow duration curves for stream gauge station upstream from the mouth of the watershed. However, systematic prediction errors generated by the watershed model during significant events precluded the analysis using daily streamflow data. Alternatively, 3-day averages 共as a surrogate to daily flow兲 were used to mitigate the effect of large errors during significant events and construct the time series model. Conversion of flow metrics from 3-day averages to daily was achieved by means of simple linear regression relationships. An ARIMA time series model was fit to 3-day averaged model residual errors and a split-sample approach was implemented to construct and validate the error model. Two methods were used to construct the uncertainty band for model forecast, the l-step forecast and the semiparametric model. The former relies on the assumption of Gaussian white noise 共random term in the linear time series model兲, while the latter is independent of the distribution of the noise term and, therefore, is more versatile. It was shown that the application of the l-step forecast could lead to an underestimation of the uncertainty band when compared to the wider band generated by the semiparametric model; however, both methods produced comparable results. Monte Carlo simulations showed that long-term annual maximum daily flow had the highest uncertainty 共measured by COV兲 among four flow attributes. The long-term monthly median of daily flows 共a wild trout habitat metric兲 and average daily flow 共base flow during interstorm periods兲 showed lower uncertainty. The higher coefficient of variation 共COV兲 values associated with simulated highflow attributes indicate a level of uncertainty to be communicated in terms of a greater risk factor for the design and performance flood control measures. The versatility of the synthesis methodology—by means of deterministic SWAT and stochastic residual error model simulations—was further highlighted by the construction of ensemble flow duration curves, which allowed for the estimation of the 共range兲 of flow for a specified design recurrence 共return兲 period with 95% confidence. The simulated low flows 共⬍2 m3 / s兲 at the site showed highest uncertainty, which may be attributed to a combination of both inadequate model account for groundwater discharge and errors in the baseflow separation method, or inadequacy of the methodology for low flows. The computed streamflow prediction uncertainties and the associated flow attributes provide the basis for communicating the risk to water resources managers and decision makers to examine alternative management plans for sustainable water resources in the watershed. The presented results, while specific to the study watershed, demonstrated a methodology that combines the relative simplicity of classical time series models with the versatility of nonparametric probabilistic analysis to account for all sources of model and data uncertainty. The approach, however, is suitable for continuous time simulations and contingent upon a successful model calibration.

Acknowledgments The U.S. Environmental Protection Agency through its Office of Research and Development partially funded and managed the research described here through in-house efforts and in part by an appointment to the Postgraduate Research Program at the National Risk Management Research Laboratory. This program is administered by the Oak Ridge Institute for Science and Education through an interagency agreement between the U.S. Department of Energy and the U.S. Environmental Protection Agency. It has not been subjected to Agency review and therefore does not necessarily reflect the views of the Agency, and no official endorsement should be inferred.

References Ajami, N. K., Duan, Q., and Sorooshian, S. 共2007兲. “An integrated hydrologic Bayesian multimodel combination framework: Confronting input, parameter, and model structural uncertainty in hydrologic prediction.” Water Resour. Res., 43, 1. Ang, A. H.-S., and Tang, W. H. 共2007兲. Probability concepts in engineering, emphasis on applications in civil and environmental engineering, Wiley, New York. Arabi, M., Govindaraju, R. S., and Hantush, M. M. 共2007兲. “A probabilistic approach for analysis of uncertainty in evaluation of watershed management practices.” J. Hydrol., 333, 459–471. Bates, B. C., and Campbell, E. P. 共2001兲. “A Markov chain Monte Carlo scheme for parameter estimation and inference in conceptual rainfallrunoff modeling.” Water Resour. Res., 37共4兲, 937–948. Beck, M. B. 共1987兲. “Water quality modeling: A review of the analysis of uncertainty.” Water Resour. Res., 23共8兲, 1393–1442. Beven, K. 共1993兲. “Prophecy, reality and uncertainty in distributed hydrologic modeling.” Adv. Water Resour., 16, 41–51. Beven, K., and Binley, A. 共1992兲. “The future of distributed models: Model calibration and uncertainty prediction.” Hydrolog. Process., 6, 279–298. Beven, K., and Freer, J. 共2001兲. “Equifinality, data assimilation, and uncertainty estimation in mechanistic modeling of complex environmental system using the GLUE methodology.” J. Hydrol., 249, 11–29. Binley, A. M., Beven, K. J., Calver, A., and Watts, L. G. 共1991兲. “Changing responses in hydrology: Assessing the uncertainty in physically based model predictions.” Water Resour. Res., 27共6兲, 1253–1261. Box, G. E. P., and Jenkins, G. M. 共1970兲. Time series analysis, forecasting, and control, Holden Day, San Francisco. Bruen, M., and Yang, J. 共2006兲. “Combined hydraulic and black-box models for flood forecasting in urban drainage systems.” J. Hydrol. Eng., 11共6兲, 589–596. Carpenter, T. M., and Georgakakos, K. P. 共2004兲. “Impacts of parametric and radar rainfall uncertainty on the ensemble streamflow simulations of a distributed hydrologic model.” J. Hydrol., 298, 202–221. Dilks, D. W., Canale, R. P., and Meier, P. G. 共1992兲. “Development of Bayesian Monte Carlo techniques for water quality model uncertainty.” Ecol. Modell., 62, 149–162. Freer, J., and Beven, K. 共1996兲. “Bayesian estimation of uncertainty in runoff prediction and the value of data: An application of the GLUE approach.” Water Resour. Res., 32共7兲, 2161–2173. Gupta, H. V., Sorooshian, S., and Yapo, P. O. 共1998兲. “Toward improved calibration of hydrologic models: Multiple and noncommensurable measures of information.” Water Resour. Res., 34共4兲, 751–763. Haan, C. T. 共2002兲. Statistical methods in hydrology, 2nd Ed., Iowa State Press, A Blackwell Publishing Company, Ames, Iowa. Hantush, M. M., and Kalin, L. 共2005兲. “Uncertainty and sensitivity analysis of runoff and sediment yield in a small agricultural watershed with KINEROS2.” Hydrol. Sci. J., 50共6兲, 1151–1171. Hossain, F., and Anagnostou, E. N. 共2005兲. “Assessment of a probabilisJOURNAL OF HYDROLOGIC ENGINEERING © ASCE / JULY 2008 / 595

tic scheme for flood prediction.” J. Hydrol. Eng., 10共2兲, 141–150. Hossain, F., Anagnostou, E. N., Dinku, T., and Borga, M. 共2004兲. “Hydrological model sensitivity to parameter and radar rainfall estimation uncertainty.” Hydrolog. Process., 18, 3277–3291. Kalin, L., and Hantush, M. M. 共2006a兲. “Effect of urbanization on sustainability of water resources in the Pocono Creek watershed.” V. P. Singh, and Y. J. Xu, eds., Coastal hydrology and processes, American Institute of Hydrology, Water Resources Publications, Highlands Ranch, Colo., 59–70. Kalin, L., and Hantush, M. M. 共2006b兲. “Hydrologic modeling of the Pocono Creek watershed with NEXRAD and rain gauge data.” J. Hydrol. Eng., 11共6兲, 555–569. Kavetski, D., Kuczera, G., and Franks, S. W. 共2006兲. “Bayesian analysis of input uncertainty in hydrologic modeling. I: Theory.” Water Resour. Res., 42, 1. Kim, T.-W., and Valdés, J. B. 共2005兲. “Synthetic generation of hydrologic time series based on nonparametric random generation.” J. Hydrol. Eng., 10共5兲, 395–404. Kim, T.-W., Valdés, J. B., and Yoo, C. 共2006兲. “Nonparametric approach for bivariate drought characterization using palmer drought index.” J. Hydrol. Eng., 11共2兲, 134–143. Kuczera, G., Kavetski, D., Franks, S., and Thyer, M. 共2006兲. “Towards a Bayesian total error analysis of conceptual rainfall-runoff models: Characterizing model error using storm-dependent parameters.” J. Hydrol., 331, 161–177. Lall, U. 共1995兲. “Recent advances in nonparametric function estimation: Hydrologic applications.” Rev. Geophys., 33共S兲, 1093–1102. Madsen, H. 共2000兲. “Automatic calibration of a conceptual rainfall-runoff model using multiple objectives.” J. Hydrol., 235, 276–288. McKay, M. D., Beckman, R. J., and Conover, W. J. 共1979兲. “A comparison of three methods for selecting values of input variables in the analysis of output from a computer code.” Technometrics, 21共2兲, 239–245. Melching, C. S., and Bauwens, W. 共2001兲. “Uncertainty in coupled nonpoint source and stream water-quality models.” J. Water Resour. Plann. Manage., 127共6兲, 403–413. Muleta, M. K., and Nicklow, J. W. 共2005兲. “Sensitivity and uncertainty analysis coupled with automatic calibration for a distributed watershed model.” J. Hydrol., 306, 127–145. National Research Council 共NRC兲. 共2001a兲. Assessing the TMDL approach to water quality management, National Academy Press, Washington, D.C. National Research Council 共NRC兲. 共2001b兲. Envisioning the agenda for water resources research in the twenty-first century, National Academy Press, Washington, D.C. National Research Council 共NRC兲. 共2004兲. Confronting the nation’s water problems: The role of research, National Academies Press, Washington, D.C. Reckhow, K. H. 共1994a兲. “Forum: Importance of scientific uncertainty in

596 / JOURNAL OF HYDROLOGIC ENGINEERING © ASCE / JULY 2008

decision making.” Environ. Manage. (N.Y.), 18共2兲, 161–166. Reckhow, K. H. 共1994b兲. “Water quality simulation modeling uncertainty analysis for risk assessment and decision making.” Ecol. Modell., 72, 1–20. Refsgaard, J. C., Henriksen, H. J., Harrar, W. G., Scholten, H., and Kassahun, A. 共2005兲. “Quality assurance in model based water management—Review of existing practice and outline of new approaches.” Environ. Modell. Software, 20, 1201–1215. Samanta, S., Mackay, D. S., Clayton, M. K., Kruger, E. L., and Ewers, B. E. 共2007兲. “Bayesian analysis for uncertainty estimation of a canopy transpiration model.” Water Resour. Res., 43, 1. Santhi, C., Arnold, J. G., Williams, J. R., Dugas, W. A., Srinivasan, R., and Hauck, L. M. 共2001兲. “Validation of the SWAT model on a large river basin with point and nonpoint sources.” J. Am. Water Resour. Assoc., 37共5兲, 1169–1188. Schulz, K., Beven, K., and Huwe, B. 共1999兲. “Equifinality and the problem of robust calibration in nitrogen budget simulations.” Soil Sci. Soc. Am. J., 63, 1934–1941. Shirmohammadi, A., et al. 共2006兲. “Uncertainty in TMDL models.” Trans. ASABE, 49共4兲, 1033–1049. Shumway, R. H. 共1988兲. Applied statistical time series analysis, PrenticeHall, Englewood Cliffs, N.J. Silverman, B. W. 共1986兲. Density estimation for statistics and data analysis, Chapman and Hall, New York. Singh, V. P., and Woolhiser, D. 共2002兲. “Mathematical modeling of watershed hydrology.” J. Hydrol. Eng., 7共4兲, 270–292. Sorooshian, S., and Dracup, J. A. 共1980兲. “Stochastic parameter estimation procedures for hydrologic rainfall-runoff models: Correlated and heteroscedastic error cases.” Water Resour. Res., 16共2兲, 430–442. Spear, R. C., and Hornberger, G. M. 共1980兲. “Eutrophication in peel inlet—II. Identification of critical uncertainties via generalized sensitivity analysis.” Water Res., 14, 43–49. van der Perk, M., and Bierkens, M. F. P. 共1997兲. “The identifiability of parameters in a water quality model of the Biebrza, Poland.” J. Hydrol., 200, 307–322. Vrugt, J. A., Gupta, H. V., Bastidas, L. A., Bouten, W., and Sorooshian, S. 共2003a兲. “Effective and efficient algorithm for multiobjective optimization of hydrologic models.” Water Resour. Res., 39共8兲, 1214. Vrugt, J. A., Gupta, H. V., Bouten, W., and Sorooshian, S. 共2003b兲. “A shuffled complex evolution metropolis algorithm for optimization and uncertainty assessment of hydrologic model parameters.” Water Resour. Res., 39共8兲, 1201. Wagener, T., Boyle, D. P., Lees, M. J., Wheater, H. S., Gupta, H. V., and Sorooshian, S. 共2001兲. “A framework for development and application of hydrologic models.” Hydrology Earth Syst. Sci., 5共1兲, 13–26. Zhang, H. X., and Yu, S. L. 共2004兲. “Applying the first-order error analysis in determining the margin of safety for total maximum daily load computations.” J. Environ. Eng., 130共6兲, 664–673.