Bayesian model for flow-class dependent distributions of fecal-indicator bacterial concentration in surface waters

water research 44 (2010) 1006–1016

Available at www.sciencedirect.com

journal homepage: www.elsevier.com/locate/watres

Bayesian model for flow-class dependent distributions of fecal-indicator bacterial concentration in surface waters Mary E. Schoen a,b,*, Mitchell J. Small a,b, Jeanne M. VanBriesen b,c a

Department of Engineering and Public Policy, Carnegie Mellon University, Pittsburgh, PA 15213-3890, USA Department of Civil and Environmental Engineering, Carnegie Mellon University, Pittsburgh, PA 15213-3890, USA c Department Biomedical Engineering, Carnegie Mellon University, Pittsburgh, PA 15213-3890, USA b

article info

abstract

Article history:

A Bayesian statistical water quality model is demonstrated to predict fecal-indicator

Received 25 March 2009

bacterial concentrations for waterbodies without sufficient monitoring data for data-

Received in revised form

intensive modeling techniques. Using a truncated bivariate normal likelihood function and

21 September 2009

the readily available observations of flow and bacterial concentration, the Bayesian

Accepted 19 October 2009

approach propagates the uncertainty in the model parameterization to the final predictions

Available online 10 November 2009

of in-stream bacterial concentration. The proposed model captures the variation in the magnitude of bacterial loading between dry and wet conditions by estimating separate sets

Keywords:

of model parameters for different flow conditions, but also has the capability to pool the

Fecal indicator

data among flow conditions. The model can be used in two ways: first, the model specifies

Bayesian model

the percent of time that the recreational season in-stream concentration is not in

Surface water

compliance with the relevant water quality standard, and second, the model estimates the necessary bacterial load reduction for multiple flow conditions to meet the relevant water quality standard. Using an eleven year monitoring record for a site sampled at a monthly frequency on the Youghiogheny River in southwestern Pennsylvania, USA, the model parameters are updated and posterior predictions are generated for each 2-year increment. After six years of sampling, the predicted percent of time that the recreational season instream bacterial concentration is not in compliance with the relevant water quality is 82% with 95% CI(80,85), and the bacterial load reductions required to meet the standard are 14.7 CI(14.6,14.8), 14.5 CI(14.3, 14.6), and 14.0 CI(13.9, 14.2) log10(cfu/day) for the high, normal, and dry flow conditions, respectively. The change in estimated load reduction and percent exceedance resulting from additional monitoring for this site becomes small after six years of sampling, indicating that a decision does not need to be postponed for additional monitoring. Published by Elsevier Ltd.

1.

Introduction

To designate a waterbody for recreational use in the US, the waterbody must be in compliance with the applicable water quality standard for fecal-indicator bacterial (FIB)

concentration, as governed by the U.S. Environmental Protection Agency (U.S. EPA). However, the frequency at which recreational waters are monitored for fecal-indicator bacteria meets the minimum sampling frequency specified by the U.S. EPA of five water quality samples per month in only

* Corresponding author at: 26 W. Martin Luther King Drive, Cincinnati, OH 45268, USA. Tel.: þ1 513 569 7647; fax: þ1 513 569 7464. E-mail address: [email protected] (M.E. Schoen). 0043-1354/$ – see front matter Published by Elsevier Ltd. doi:10.1016/j.watres.2009.10.016

water research 44 (2010) 1006–1016

a few instances in high-profile streams and rivers (see, for example, the National Water Quality Inventory Report to Congress, http://www.epa.gov/305b/). The sparse datasets created by the low-frequency monitoring that occurs at most locations could be used to inform management decisions through a calibrated water quality model; however, existing water quality models are often data intensive and difficult to calibrate given such limited input information (Borah and Bera, 2003). The motivation of the present work is the need for statistical water quality models that can complement more mechanistic approaches when only sparse datasets are available to provide predictions for in-stream bacterial concentrations. Process-based water quality models incorporate knowledge of loading, transport, growth, death, deposition, and mobilization of bacteria in a watershed to predict in-stream FIB concentrations at specific time intervals. To build a process-based model, detailed loading information is required along with an intensive monitoring record for accurate calibration (Benham et al., 2006; Borah and Bera, 2003, 2004; Donigian and Huber, 1991; Horn et al., 2004; Wilkinson et al., 1995). Although these models are useful in exploring management options, the intensive data requirements preclude the selection of process-based models for waterbodies with sparse bacterial monitoring records. Furthermore, the process-based models are difficult to re-calibrate as additional data accumulates as part of a long-term monitoring program. An alternative regression approach has been used to predict bacterial concentrations for specified conditions (Eleria and Vogel, 2005; Ferguson et al., 1996; Francy et al., 2000; Mas and Ahlfeld, 2007) but also requires multiple observable quantities including rainfall, chemical parameters, and streamflow (Eleria and Vogel, 2005). The Bayesian, stochastic models presented below are different in that instead of specifying parameters or explanatory variables to account for the effects of fate and transport processes on the FIB concentration, the models make inferences about the distribution of FIB concentration directly from the water quality observations and can be easily updated as observations accumulate over time. Using a Bayesian modeling approach, the observable quantities of interest y, such as FIB concentration, as well as unobservable model parameters q, such as the mean of concentration, are represented with probability models. Inferences about both y and q are made using the available data and a prior distribution p(q) which represents our uncertainty about the unobservable model parameters (Gelman et al., 2004). Whereas a frequentist approach predicts y using point estimates of q estimated from the observations alone, the resulting predictions of y using a Bayesian approach incorporate the uncertainty from the model parameters in the form of a posterior predictive distribution. For sparse datasets, this approach has the advantage of tracking the changes in the uncertainty of the final predictions of in-stream bacterial concentration from the model parameterization as observations accumulate over time. The main objective of this study is to demonstrate how the proposed Bayesian model and resulting posterior predictive distributions of FIB concentration can be used for water quality decision making at a site with a long-term, sparse

1007

monitoring dataset. The model can be used in two ways: first, the model specifies the percent of time that the recreational season in-stream concentration is not in compliance with the relevant water quality standard, and second, the model estimates the necessary bacterial load reduction to meet the relevant water quality standard. Two Bayesian models are developed and compared for a monitoring station on the Youghiogheny River. The models use the relationship between flow and concentration to predict the possible distribution of in-stream fecal-indicator bacterial concentration over the course of the recreational season. The models capture the variation in the magnitude of bacterial loading between dry and wet conditions by estimating separate sets of model parameters for different flow conditions, but have the capability to pool the data among flow conditions to inform the uncertainty in the concentration of fecal-indicator bacteria for a particular streamflow using a hierarchical structure. The Bayesian model with the best performance measures is subsequently used for water quality assessment and the estimation of required bacterial load reductions.

2. Bayesian approach in water quality modeling The qualities of Bayesian data analysis, including the probabilistic representation of quantities of interest and propagation of uncertainty, make it a powerful tool for water quality planning (Hobbs, 1997). Bayes rule has been applied in the water quality modeling calibration process to select improved estimates of model parameters by formulating posterior distributions of those parameters given additional data (Dilks et al., 1992; Freer et al., 1996; Gronewold et al., 2009; Harrison, 2007; Kanso et al., 2003; Liu et al., 2008; Qian and Reckhow, 2007; Qian et al., 2003). For example, Bayesian Monte Carlo techniques can be used to calibrate the parameters of simple, process-based water quality models by using the likelihood function to measure the acceptability of a parameter set given the observed data (Dilks et al., 1992; Freer et al., 1996; Kanso et al., 2003; Qian et al., 2003; Sohn et al., 2000). This technique allows for the uncertainty in the estimated parameters to propagate through the model. The technique has been extended to calibrate more complex models for eutrophication (Arhonditsis et al., 2008, 2007, 2006) and dissolved oxygen (Harrison, 2007). Bayesian Networks are an alternative approach to the basic Bayesian Monte Carlo techniques and useful for modeling multiple, linked water quality processes (Borsuk et al., 2004; Stiber et al., 1999, 2004; Varis, 1998), but are limited in the ability to represent complex probability structures. Hierarchical Bayesian analysis allows additional predictive power when making posterior inferences for model parameters that are believed to be related in some way. Using a hierarchical structure, the parameters (q1,q2,.qc) for separate sets of data ( y1,y2,.yc) are related through a common population distribution, where each parameter qc is modeled as an independent draw from the population distribution. This structure allows posterior estimates for each qc while pooling data among the separate sets. Pooling data between sets is

1008

water research 44 (2010) 1006–1016

County, PA are used (EPA, 2007). Monthly fecal coliform samples were recorded at this site from 1977 through 1987 over the recreational season of May through September. In total, there are 134 samples (three samples per day at 5–10 min intervals); all values are reported as colony forming units (cfu) per 100 ml. The site is located 250 m downstream from the USGS flow gauge 03083500. Since the Youghiogheny River is large (a 1715 square mile watershed) and relatively uniform, the daily-average flows in cubic feet per second (cfs) were used in a 20-year period-of-record analysis to capture the range of possible streamflows during the recreational seasons of 1967– 1987. Fig. 1 displays scatterplots of log10 flow and fecal coliform concentration for three equal probability flow conditions using the entire dataset.

especially useful for modeling applications that have limited available data for a particular set or sets of interest, such as in cross-sectional predictive studies (Dominici et al., 2000; Goyal et al., 2005; Helser and Lai, 2004; Malve and Qian, 2006; Wikle, 2003). Parameters estimated for data-sparse subpopulations are able to ‘‘borrow strength’’ from subsets with more complete or extensive samples (Borsuk et al., 2001b; Ghosh and Rao, 1994). The hierarchical Bayesian approach has been used in water quality modeling applications for cross-sectional analysis. Cross-sectional studies use multiple independent sites to create predictive relationships between observed, explanatory variables and variables of interest. Malve and Qian used the Finnish Lake monitoring network data in a hierarchical model to calculate the probability of the chlorophyll a concentration exceeding the standard for an individual lake under different nitrogen and phosphorus concentration combinations (2006). Using a hierarchical structure, the parameters of the linear regressions relating chlorophyll a concentration to total phosphorus and total nitrogen for multiple lake sites shared a common population distribution. The hierarchical model, by pooling information among lakes, performed better in performance criteria when compared to non-hierarchical, lake specific approaches (Malve and Qian, 2006). To estimate nutrient load reductions in the Neuse Estuary, a hierarchical approach was used to predict sediment oxygen demand as a function of changes in organic matter loading. Again, the hierarchical structure allowed for site specific predictions for 34 estuarine and coastal systems, but pooled information from the entire set of case studies to estimate the parameters of the sediment oxygen demand model (Borsuk et al., 2001a). Hierarchical Bayesian analysis has potential application beyond cross-sectional eutrophication water quality analysis and is evaluated here as an alternative or compliment to process-based models for prediction of fecal-indicator bacterial concentration using sparse datasets.

3.

4.

Model development

Two models are presented and compared for predicting the distribution of fecal-indicator bacterial concentration over the recreational season. Similar to a regression approach, the models use the relationship between flow and concentration to predict fecal-indicator bacterial concentrations. The relationship between flow and concentration may change based on the type of loading. For example, when the bacterial loading is dominated by point sources during times of low flow, the relationship between flow and concentration may be negative (as flow increases, concentration decreases). However, when bacterial loading is dominated by storm runoff during times of high flow, the correlation may be positive (as flow increases, concentration increases). Since the relationship between flow and concentration changes with the type of loading, separate relationships are modeled for different conditions of flow (i.e., low, normal, and high), and used to estimate the overall distribution of in-stream concentration during the recreational season. The first model M1 is a Bayesian model based on the assumption of a truncated bivariate normal relationship between flow and concentration under different flow conditions. The second model M2, the hierarchical Bayesian model, is similar to M1 with an added level of hierarchy that draws the precision parameters of each flow condition from a common population distribution.

Raw data

To demonstrate the applicability of the proposed Bayesian models, the data from Pennsylvania bacterial monitoring station WQS706 on the Youghiogheny River in Westmoreland

2

3

log10(C)

1

2.8

3.2

3.8

4.0

2.8

3.2

3.8

4.0

2.8

3.2

3.8

4.0

log10( ) Fig. 1 – Scatterplots of recreational season fecal coliform concentration and flow. Scatterplots for each of the three flow conditions, i.e. low (1), normal (2), and high (3), show the change in the range of fecal coliform concentration (log10(C )) observed over changes in flow (log10(Q)) at WQS706.

1009

water research 44 (2010) 1006–1016

The form of the relationship between flow and concentration in both models M1 and M2 is based on the bivariate " normal probability model, ( 1 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiexp  fc ðx; yÞ ¼ 2ð1  r2 Þ 2psx sy ð1  r2 Þ  2     #)  y  my y  my 2 x  mx x  mx  þ (1) 2r sx sx sy sy where x represents the natural log of flow and y represents the natural log of concentration. Since we want to represent the different possible relationships between flow and concentration as the magnitude of flow changes, we divide the bivariate distribution into C truncated sections at flows (x1,x2,.xc1) and estimate r separately for each section. By dividing the flow into sections, or conditions, each with equal probability P, the joint distribution for each condition fT,c(x,y) becomes a singly truncated bivariate normal distribution, expressed here as the product of the marginal distribution of x, fx(x), and the conditional distribution of y, fyjx;c ðyjxÞ, fT;c ðyjxÞ ¼

( 2    1 1 x  mx 1 pffiffiffiffiffiffi exp  pffiffiffiffiffiffi exp  2 syjx;c sx P 2p syjx;c 2p 2 )  y  m 1 yjx;c  2 syjx;c

(2)

where



LNðxci ÞwTN ux ; s2x ; lc ; r

(10)



LN yci wN myjx;ci ; s2yjx;c

(11)

uyjx;ci ¼ ac þ bc  ðxci  xc Þ

(12)

where LN(xci) and LN( yci) are the ith observations of transformed variables LN(Q) and LN(C ) from flow condition c. The parameters of the marginal truncated normal distribution of LN(xci) in Eq. (10) include qx ¼ ðmx ; s2x ; lc ; rc Þ where (lc,rc) specify the left and right truncation points for each condition. The parameters qx are not estimated using the Bayesian approach; instead, these values are constants and entered into the model based on maximum likelihood approximations from the long-term flow record. The parameters of the conditional distribution of LN( yci) in Eq. (11) are updated; the parameters qyjx;ci ¼ ðmyjx;ci ; s2yjx;c Þ include the conditional mean for each observation i in each condition c and the conditional variance in each flow condition. Again, only parameters related to y are subject to Bayesian analysis; therefore, xc is entered into the model as a vector of constants estimated from the long-term flow record. The parameters (ac,bc,sc) are assigned non-informative priors,

uyjx;c ¼ ac þ bc ðx  xc Þ

(3)

ac wNð0:0; 0:01Þ

(13)

r sy;c bc ¼ c sx

(4)

bc wNð0:0; 0:01Þ

(14)

ac ¼ my;c  bc ðmx  xc Þ

(5)



s2yjx;c ¼ s2y;c 1  r2c

syjx;c wUð0:0; 10:0Þ

(15)

(6)

and for a given sample of n specimens in the truncated section, P x ¼ n1 xi =n. In this form, the conditional distribution of y is normally distributed with a standard deviation syjx;c and mean uyjx;c . The truncated bivariate normal distribution is flexible in that when the relationship between flow and concentration is weak with a correlation near zero, the conditional variance and mean of concentration calculated in Eqs. (6) and (3) approach the independent mean and variance of concentration. Solving for Eqs. (3)–(6), the non-conditional mean and standard deviation of y along with the correlation coefficient for each flow condition are equivalent to my;c ¼ ac þ bc ðmc  xc Þ sy;c ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s2yjx;c þ s2x b2c

sx bc rc ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 syjx;c þ s2x b2c

The priors for (ac,bc) are relatively flat normal distributions with variances of 100 while the prior for sc is a wide uniform distribution with an upper limit of 10. The selection of noninformative priors ensures that the posterior inference is dominated by the evidence provided by the data. For watersheds with additional information about water quality, an informative prior may be specified in order to combine the observational data with the additional information, such as modeling output or expert opinion (Arhonditsis et al., 2007; Gronewold et al., 2009; Qian et al., 2005; Wikle, 2003). Additional parameters of interest, my,c and sy,c, are calculated based on those subject to Bayesian inference using Eqs. (7) and (8).

(7)

4.2.

Hierarchical Bayesian model M2

(8) (9)

where mx and sx are the mean and standard deviation of the non-truncated  distribution (Cohen, 1955).

The hierarchical Bayesian model M2 adds an additional hierarchical structure to M1 in order to pool information among flow conditions c. The hierarchical Bayesian model samples the conditional standard deviation syjx;c from a common population distribution, syjx;c wUða; bÞ

4.1.

Bayesian model M1

The Bayesian model M1 estimates the parameters for the singly truncated bivariate normal distribution for each flow condition c as follows:

This structure represents each syjx;c uniform population distribution in information from all flow conditions limits of the uniform distribution in informative priors,

(16) as a sample from the Eq. (16), which pools c. The lower and upper Eq. (16) are given non-

1010

water research 44 (2010) 1006–1016

awUð0:0; 10:0Þ

(17)

bwParð10:0; aÞ

(18)

The lower limit a is modeled as uniform in Eq. (17) with an upper limit of 10. The upper limit b is modeled as Pareto with a scale parameter set at 10 and shape parameter set as a. The use of the Pareto distribution for b ensures that the upper limit is greater than the lower limit a. This additional level of hierarchy is introduced to improve prediction when samples are missing from a flow condition of interest. Graphical depictions of the models along with WinBUGs code are presented in the Appendix.

4.3.

Model assumptions

Both models M1 and M2 share common assumptions related to model structure. First, both models assume that the longterm daily-average distribution of flow is well represented as a lognormal distribution. This structural assumption is typical for period-of-record flow frequency curve analysis, where the flow cumulative distribution function is an average curve based on the period of time for which it is constructed. In a review of statistical approaches to fitting flow frequency curves, Riggs notes that no one distribution gives a perfect fit (1989), but the lognormal distribution is most commonly used in water resources applications (Fennessey and Vogel, 1990). Second, both models assume that the distribution of concentration within each flow condition is lognormal (Ott, 1995; Vogel et al., 2005). Exploratory data analysis of the applied water quality sample in the ‘‘Raw Data’’ section supports this assumption. M2 adds an additional assumption that the hierarchical prior population distribution is uniformly distributed. A posterior predictive check of my,c and sy,c is presented in the ‘‘Model Performance’’ section to lend support to the structural assumptions of M1 and M2. Third, both models divide the distribution of flow into three conditions. Splitting the data into flow conditions involves two considerations. The divisions should be meaningful from a management perspective. For example, the divisions of high, normal, and low flows represent possible different flow-related sources of contamination, being nonpoint, mixed, and point sources. Also, each condition must have enough data for the Bayesian model to converge.

5.

Computation

Analytic solutions for posterior parameter estimates are difficult to specify given the complex nature of the model structure; instead, a Markov chain Monte Carlo (MCMC) simulation method is implemented using the WinBUGs software (Spiegelhalter et al., 1999). The MCMC method simultaneously estimates the distribution of unknown parameters by sampling the parameters from their joint posterior distribution. Three chains with over-dispersed starting points of length 40,000 were used for posterior inference. The burn-in for both models was 20,000 iterations, based on the Gelman and Rubin Diagnostic (Gelman and Rubin, 1992), which left 20,000 iterations for posterior inference. Of the 20,000

iterations, a thinning procedure with a thin rate of 60 reduced the sample down to 1000 iterations for each parameter, removing any significant autocorrelation from the final MCMC chain. The resulting chains of samples for each parameter are used for posterior inference. The posterior predictive distributions of FIB concentration in each flow condition are used to specify (1) the percent of time that the recreational season instream bacterial concentration is not in compliance with the relevant water quality standard, and (2) the 90th percentile load reduction necessary to meet the water quality standard. To calculate the posterior predictive distribution of fecalindicator bacterial concentration, 1000 flow samples are drawn from each of the equal probability flow conditions. Next, a concentration is sampled from each of the 1000 posterior MCMC parameter estimates of myjx;c and syjx;c , conditional on the flow condition c and sampled flow. The distribution of the combined concentration samples from all equal probability flow conditions is the posterior predictive distribution of the recreational season fecal-indicator bacterial concentration. The 90th percentile load reduction is calculated by sampling the 90th percentile conditional concentration from each of the 1000 posterior MCMC parameter estimates of myjx;c and syjx;c . This posterior predictive 90th percentile conditional concentration is then multiplied with flow and the allowable in-stream load for that flow is subtracted. In addition, flow-condition-specific posterior predictive distributions of the recreational season concentration are calculated.

6.

Results

6.1.

Model performance

Table 1 shows the mean M1 and M2 posterior parameter estimates for my;c , syjx;c , rc and sy,c along with 95% credible intervals (CIs). Both models estimate higher average concentration for the high flow condition, mM1 y;3 ¼ 3:402 CIð3:039; 3:553Þ ¼ 3:289 CIð2:920; 3:653Þ, compared to the normal flow and mM2 y;3 ¼ 2:842 CIð2:495; 3:355Þ and condition estimates, mM1 y;2 ¼ 2:929 CIð2:576; 3:286Þ. The correlation coefficients for all mM2 y;2 conditions are small for both models with credible intervals that range from negative to positive. The weak relationships between flow and concentration result in little difference between the estimates of syjx;c and sy,c for a particular flow condition and model. When comparing the sM1 y;c among conditions, the high flow mean posterior estimate of sM1 y;3 ¼ 0:481 CIð0:380; 0:609Þ is smaller than the 95% CIs for the low and normal flow conditions, sM1 y;1 ¼ 0:847 CIð0:610; 1:289Þ and sM1 y;2 ¼ 0:925 CIð0:630; 1:620Þ. This pattern in spread is evident in the raw data scatterplots of Fig. 1. The posterior estimate of sM2 y;3 ¼ 0:640 CIð0:551; 0:761Þ has a credible interval that overlaps with both estimates of and sM2 sM2 y;1 ¼ 0:778 CIð0:599; 1:247Þ y;2 ¼ 0:827 CIð0:604; 1:137Þ. The posterior parameter confidence intervals overlap due to the hierarchical structure that samplessM2 yjx;c from a common population distribution. The posterior distribution of sM2 y;3 M1 shifts toward the mean posterior estimates of sM1 y;1 and sy;2 since the high flow condition has the fewest samples. The

1011

water research 44 (2010) 1006–1016

Table 1 – Posterior parameters estimates for Bayesian model M1 and the hierarchical model M2. Mean posterior parameter estimates for bacterial concentration and correlation coefficient (with flow), mean, standard deviation, and conditional standard deviation (for flow x) are presented in log10(cfu/100 ml) for each of the three flow conditions c, i.e. low (1), normal (2), and high (3), with 95% credible intervals shown in parenthesis. Parameter

rc m y, c s y, c syjx;c

M1

M2

1 n ¼ 43

2 n ¼ 57

3 n ¼ 34

1 n ¼ 43

2 n ¼ 57

3 n ¼ 34

0.253 (0.605, 0.828) 3.011 (2.274, 3.767) 0.847 (0.610, 1.289) 0.717 (0.586, 0.886)

0.101 (0.834, 0.869) 2.842 (2.495, 3.355) 0.925 (0.630, 1.620) 0.728 (0.597, 0.892)

0.227 (0.536, 0.154) 3.402 (3.039, 3.553) 0.481 (0.380, 0.609) 0.459 (0.363, 0.579)

0.311 (0.623, 0.843) 3.034 (2.323, 3.723) 0.778 (0.599, 1.247) 0.652 (0.599, 1.247)

0.084 (0.816, 0.868) 2.929 (2.576, 3.286) 0.827 (0.604, 1.374) 0.654 (0.604, 1.137)

0.016 (0.499, 0.219) 3.289 (2.920, 3.653) 0.640 (0.551, 0.761) 0.619 (0.551, 0.761)

hierarchical structure also results in mean posterior estimates M2 for sM2 y;1 and sy;2 at the lower end of the M1 95% CIs. The differences in prediction between M1 and M2 are illustrated by the posterior predictive check process. Posterior predictive checks provide evidence that the observed data support the key model structural assumptions. The basic idea is that if the model assumptions are valid, the observed data should look plausible under the posterior predictive distribution. Any observed discrepancy can be due to poor model structural assumptions or chance. Plausibility is informed by the Bayesian p-value. The Bayesian p-value is defined as the probability that the replicated data could be more extreme than the observed data. A model is suspect if the Bayesian p-value of the observed data is close to 0 or 1, thereby indicating that the observed pattern would be unlikely in replication of the data if the model were true (Gelman et al., 2004). Figs. 2 and 3 show the posterior predictive checks for my,c and sy,c, the mean and standard deviation of concentration for each flow condition for M1 and M2. The posterior predictive distributions are compared to the maximum likelihood parameter estimates for the truncated bivariate normal distribution using the data, _

m

y;c

_

s

y;c

  sy;c ðxc  ux;c Þ ¼ yc þ rc sx;c vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  u u1  1  s2x;c 1  r2

u c s2x;c  ¼ sy;c u t  s2x;c 1  1  s2

(19)

(20)

x;c

where rc is the correlation coefficient and bars over the parameters (such as rc ) indicate that only observed data from within the flow condition are used in estimation (Cohen, 1955). The posterior predictive distributions, observed maximum likelihood estimates, and Bayesian p-values for my,c and sy,c are shown for log10 transformed variables in Fig. 2 using the Bayesian M1. The Bayesian p-value range of 0.87 and 0.56 indicates that the observed data are plausible under the posterior predictive distributions. Fig. 3 shows similar results for M2 with the Bayesian p-values ranging between 0.77 and 0.56, indicating good model fit for all parameters but sy,3. The standard deviation of concentration in the high flow condition

has a p-value of 1. This result is attributed to the hierarchical structures tendency to homogenize model parameter estimates. A model assumption potentially influential in posterior inference is the selection of the uniform prior for the standard deviation of conditional concentration. An alternative prior for the precision (or inverse of variance) is the gamma distribution. Using the alternative gamma distribution (Gamma(1, 0.001)), the posterior parameter estimates of sy,c are M1 sM1 y;1 ¼ 0:826 CIð0:605; 1:276Þ, sy;2 ¼ 0:896 CIð0:613; 1:512Þ, and M1 sy;3 ¼ 0:460 CIð0:364; 0:589Þ and are smaller than the posterior estimates using the uniform prior. The uniform prior was selected over the gamma so that both models M1 and M2 converge each update.

6.2.

Model comparison and selection

In addition to posterior predictive check results, a second way to compare the alternative model structures is to calculate the Deviance Information Criterion (DIC) as a measure of a model’s out-of-sample predictive power, DIC ¼ D þ pD

(21)

where D is the posterior mean of the deviance (with deviance defined as 2 log pðyjqÞ), and pD is the effective number of parameters. The difference in DIC between M1 and M2 is less than 5.0, suggesting a negligible difference between the predictive powers of the two models. Based on the posterior predictive check results, M1 is selected to illustrate the use of posterior inference in compliance assessment and the calculation of load reductions.

6.3. Effect of additional monitoring on water quality management The Bayesian model M1 is used to evaluate the effect of additional monitoring on the water quality assessment for WQS706 and the estimation of bacterial load reductions. Since both assessment and load estimation depend on the posterior parameter estimates, the relative changes in the model parameters resulting from bi-yearly updates provide a first indication of possible implications for water quality management. Fig. 4 shows how the mean posterior parameter

1012

water research 44 (2010) 1006–1016

3.0

0.72

1.5

2.5

3.5

4.5

2.0 1.0 0.0

1.0

2.0

0.84

0.0

1.0

2.0

0.56

0.0

ƒMCMC(μy,c)

3

3.0

2

3.0

1

1.5

2.5

3.5

4.5

1.5

2.5

3.5

4.5

0.5

1.5

2.5

2

4

6

0.72

0

2

4

6

0.87

0

2

4

6

0.84

0

ƒMCMC(σy,c)

log10(μy,c)

0.5

1.5

2.5

0.5

1.5

2.5

log10(σy,c) Fig. 2 – Posterior predictive checks for Bayesian model M1. The posterior predictive distributions, observed maximum likelihood parameter estimates, and p-values for bacterial concentration mean my;c and standard deviation sy;c are presented for log10(cfu/100 ml) transformed variables with the flow conditions c increasing in magnitude from left to right, i.e. low (1), normal (2), and high (3).

3.0

0.76

0.56

1.5

2.5

3.5

4.5

2.0 0.0

1.0

2.0 1.0 0.0

1.0

2.0

0.77

0.0

ƒMCMC(μy,c)

3

3.0

2

3.0

1

1.5

2.5

3.5

4.5

1.5

2.5

3.5

4.5

0.5

1.5

2.5

6 2

4

1

0

2

4

6

0.68

0

2

4

6

0.63

0

ƒMCMC(σy,c)

log10(μy,c)

0.5

1.5

2.5

0.5

1.5

2.5

log10(σy,c) Fig. 3 – Posterior predictive checks for hierarchical model M2. The posterior predictive distributions, observed maximum likelihood parameter estimates, and p-values for bacterial concentration mean my;c and standard deviation sy;c are presented for log10(cfu/100 ml) transformed variables with the flow conditions c increasing in magnitude from left to right, i.e. low (1), normal (2), and high (3).

1013

6

1.0

water research 44 (2010) 1006–1016

3.0

0.5

2.5

-1.0

2

0.5

0.5

0.0

ρc

2.0 1.5

1.5

log10(σy,c)

4 3

log10(μy,c)

5

C1 C2 C3

2

4

6

8

10

2

4

6

8

10

2

4

6

8

10

Year Fig. 4 – Bi-yearly posterior parameter estimates of concentration for Bayesian model M1. Panel 1 plots mean posterior estimates of bacterial concentration mean log10 ðmy;c Þ, standard deviation log10 ðsy;c Þ, and correlation coefficient rc for each flow condition, i.e. low (C1), normal (C2), and high (C3), given bi-yearly updates of sample data along with 95% CIs. The corresponding Bayes Factors (Bi,iD2) are calculated as 0.0, 33.0, 0.7, 7.7, and 3.0, respectively.

6.4. Effect of additional monitoring on water quality assessment The posterior inferences from M1 are used to specify the percent of time that the in-stream FIB concentration exceeds the contaminant limit set by the relevant water quality standard during the recreational season. The relevant water quality standard is 25 Pa. Code x 93, which targets a fecal coliform concentration of 200 cfu/100 ml. The resulting posterior predictive cumulative distributions (CDFs) of recreational season fecal-indicator bacterial concentration using all data are displayed in Fig. 5 (boldface type) along with the CDFs of bacterial concentration for each flow condition. The concentration limit set by the recreational standard (200 cful/100 ml) is shown as a vertical line. The CDFs for M1 predict that at WQS706, the in-stream FIB concentration exceeds the allowable in-stream limit set by the standard 82.0% of the time during the recreational season, and the high flow condition shows the greatest exceedance at

97.0% of the time. Table 2 presents the median predicted percent exceedances, along with each flow condition, over time with 95% CIs shown in parenthesis. After the addition of monitoring data from years 5 and 6, the median predicted percent exceedance of 79.0% stays relatively stationary and approaches the year 11 prediction of 82.0% CI(80.0, 85.0). The predicted median percent exceedances within flow condition change by less than 2.0% after the addition of monitoring data from years 7 and 8 for the low and high flow conditions and 7.0% for the normal flow condition.

6.5. Effect of additional monitoring on the estimation of bacterial load reduction The median 90th percentile log10 daily load reductions (cfu) for each flow condition are reported in Table 3 along with 95% CIs.

P(log10( ))

estimates (plotted on the y-axis) change with each 2-year update for M1. Since the period-of-record flow distribution is assumed constant over time, the parameters of interest are those of concentration and the correlation between concentration and flow. The 1st panel tracks the mean posterior parameter estimate of log10(my,c) for each condition for bi-yearly updates of monthly monitoring data with lines representing the 95% CIs of the posterior estimates. The mean estimates for log10(my,c) are subject to large shifts with the addition of monitoring data from years 3 and 4. This shift is paired with a dramatic reduction in uncertainty. Panel 2 shows mean posterior estimates for log10(sy,c) with a similar trend at year 4. Panel 3 tracks ry,c over time; unlike the fecal-indicator bacterial concentration mean and standard deviation, the confidence intervals of the correlation coefficient for the low and normal flow conditions’ remain large, even after the fifth sampling period in year 10, based on the variability observed in the data. Taken together, the changes in model parameters over time suggest that 2 years or less of sampling at a monthly frequency may result in misleading predictions for WSQ706.

C1 C2 C3 1

2

3

4

5

6

log10( ) Fig. 5 – Posterior predictive indicator bacterial concentration CDFs for Bayesian model M1. The CDF of bacterial concentration log10(C ) is plotted in bold type along in units log10(cfu/100 ml) with the CDFs for each flow condition, i.e. low (C1), normal (C2), and high (C3). The concentration limit set by the recreational standard (200 cfu/100 ml) is shown as a vertical line.

1014

water research 44 (2010) 1006–1016

Table 2 – Predicted percent exceedance of the recreational standard over time for Bayesian model M1. The predicted median percent of time during the recreational season with bacterial indicator concentrations in exceedance of the recreational standard (200 cfu/100 ml) are presented for each flow condition, i.e. low (1), normal (2), and high (3), for each biyearly update with 95% credible intervals in parenthesis. Condition

Year

1 2 3 Overall

2

4

6

8

10

11

61% (55–66) 41% (36–47) 87% (84–91) 63% (61–67)

58% (54–64) 71% (66–75) 95% (91–97) 75% (72–77)

65% (60–71) 74% (67–78) 97% (95–99) 79% (76–81)

72% (67–76) 72% (67–78) 98% (96–100) 81% (78–83)

70% (64–76) 76% (71–81) 99% (97–100) 82% (80–84)

71% (65–76) 79% (75–83) 97% (95–99) 82% (80–85)

The median 90th percentile load reductions decrease by approximately 1 log unit for the low and high flow conditions with the addition of bacterial monitoring data from years 3 and 4. After the addition of monitoring data from years 5 and 6, the load reduction predictions and uncertainty in the predictions from each flow condition start to approach 14.0 CI(13.9,14.2), 14.5 CI(14.3,14.6), and 14.7 CI(14.6,14.8) log10 (cfu/ day).

7.

a relatively high load reduction estimate for the low flow condition in year 2 of 14.7 CI(14.5, 14.9) log10 (cfu/day) compared to the final year prediction of 14.0 CI(13.9, 14.2) log10 (cfu/day) after additional data provided evidence for a weaker relationship between flow and concentration. Once past the first three updates, the decision variables for WSQ706 remained relatively stationary over time with only a 2.0% percent change in the predicted percent exceedance of the water quality standard and less than one log unit change in the predicted load reduction. In addition, the credible intervals for the predicted decision variables were small, approximately 5.0% for the percent exceedance and less than one log unit for the load reduction. For the management purposes of assessing impairment or specifying load reductions, the observed changes and uncertainties in the decision variables are small and suggest that a decision does not need to be postponed for additional monitoring. Still, the model can be continually updated over the course of a monitoring program after a decision has been made to check for changes in the model parameterization over time. Fluctuations in model parameterization may occur over a monitoring program due to the variation inherent to datasets of small sample size or changes within the watershed. When changes within in a watershed occur over time, such as landuse and climate change, the posterior model predictions that inform management decisions may not be stationary (Milly et al., 2008). For these sites, a more complex model that

Discussion

Using a Bayesian approach, a decision maker can observe changes in the model parameter uncertainty and the ultimate effect on the predicted decision variables to inform when to make a water quality decision over the course of a monitoring program. For WSQ706, the posterior estimates of model parameters showed large fluctuations accompanied by large credible intervals in the first few bi-yearly model updates. These fluctuations, due to differences in the small bi-yearly datasets, suggest that a decision should be postponed until more data can be collected. For example, the linear relationship between flow and bacterial concentration for the low flow condition was strongly positive after only two years of sampling at a monthly frequency. Through the bivariate relationship, the low flow mean bacterial concentration and standard deviation were pulled upward. This resulting in

Table 3 – Predicted 90th percentile bacterial load reductions. The median predicted 90th percentile log10 daily bacterial load reductions necessary to meet the recreational standard (200 cfu/100 ml) are presented in units log10(cfu/day) for each flow condition, i.e. low (1), normal (2), and high (3), for each bi-yearly update along with 95% credible intervals in parenthesis. Condition

1 2 3

Year 2

4

6

8

10

11

14.7 (14.5, 14.9) 14.2 (13.9, 14.5) 16.2 (15.9, 16.6)

13.9 (13.7, 14.0) 14.9 (14.6, 15.1) 14.7 (14.6, 14.8)

14.0 (13.9, 14.2) 14.6 (14.5, 14.8) 14.8 (14.7, 14.8)

14.1 (13.9, 14.2) 14.5 (14.3, 14.7) 14.6 (14.5, 14.7)

14.0 (13.9, 14.2) 14.4 (14.3, 14.6) 14.7 (14.6, 14.8)

14.0 (13.9, 14.2) 14.5 (14.3, 14.6) 14.7 (14.6, 14.8)

water research 44 (2010) 1006–1016

incorporates temporal relationships may be necessary to predict future water quality (Wikle, 2003). Although model performance measures were used to select the ‘‘best’’ model, alternative model structures and priors could result in predictions as good as the model structure presented. For example, the alternative hierarchical model M2 has a comparable DIC to M1 and predicts a similar percent exceedance in year 11; using M1, the predicted exceedance is 82.0% CI(80.0, 85.0) and 82.2% CI(80.0, 85.0) using M2. An alternative, simple model (M3) is presented in the Appendix with an independent normal distribution representing bacterial concentration. The difference in DIC between M1 and M3 is less than 5.0, suggesting a negligible difference between the predictive powers of the two models. Using the simple model M3, the predicted exceedance is slightly higher 84.5% CI(84.2, 84.7), but still within the 95% CI of models M1 and M2.

8.

Conclusions

The proposed stochastic, Bayesian model for prediction of the in-stream bacterial concentration distribution is an alternative to more data-intensive models for sites with sparse monitoring datasets, using only flow and the observed fecal-indicator concentrations during different flow conditions. The model can be easily updated to incorporate additional data over the course of a long-term monitoring program. The posterior predictions for site WQS706 on the Youghiogheny River illustrate the effect of additional monitoring on estimating the daily bacterial load reduction and the predicted percent of time that the recreational season in-stream bacterial concentration is not in compliance with the relevant water quality standard. The change in estimated load reduction and percent exceedance resulting from additional monitoring for this site becomes small over time indicating that a decision does not need to be postponed for additional monitoring. Although the model results are not a substitution for intense monitoring, the Bayesian approach can enable the use of previously-underutilized datasets with an explicit communication of uncertainty for water quality management and decision making.

Acknowledgments Funding for this work was provided by the Center for Water Quality in Urban Environmental Systems (WaterQuest) at Carnegie Mellon University and by Three Rivers Wet Weather, Inc. and by the Pittsburgh Parks Conservancy (PPC). We would like to thank our three reviewers for valuable feedback on the content of this work.

Appendix. Supplementary information Supplementary data associated with this article can be found in the online version at doi:10.1016/j.watres.2009.10.016.

1015

references

Arhonditsis, G.B., Papantou, D., Zhang, W., Perhar, G., Massos, E., Shi, M., 2008. Bayesian calibration of mechanistic aquatic biogeochemical models and benefits for environmental management. Journal of Marine Systems 73 (1–2), 8–30. Arhonditsis, G.B., Qian, S.S., Stow, C.A., Lamon, E.C., Reckhow, K. H., 2007. Eutrophication risk assessment using Bayesian calibration of process-based models: application to a mesotrophic lake. Ecological Modelling 208 (2–4), 215–229. Arhonditsis, G.B., Stow, C.A., Steinberg, L.J., Kenney, M.A., Lathrop, R.C., McBride, S.J., Reckhow, K.H., 2006. Exploring ecological patterns with structural equation modeling and Bayesian analysis. Ecological Modelling 192 (3–4), 385–409. Benham, B.L., Baffaut, C., Zeckoski, R.W., Mankin, K.R., Pachepsky, Y.A., Sadeghi, A.M., Brannan, K.M., Soupir, M.L., Habersack, M.J., 2006. Modeling bacteria fate and transport in watersheds to support TMDLs. In: ASABE Annual International Meeting Portland, Oregon. Borah, D.K., Bera, M., 2003. Watershed-scale hydrologic and nonpoint-source pollution models: review of mathematical bases. Transactions of the ASAE 46 (6), 1553–1566. Borah, D.K., Bera, M., 2004. Watershed-scale hydrologic and nonpoint-source pollution models: review of applications. Transactions of the ASAE 47 (3), 789–803. Borsuk, M.E., Higdon, D., Stow, C.A., Reckhow, K.H., 2001a. A Bayesian hierarchical model to predict benthic oxygen demand from organic matter loading in estuaries and coastal zones. Ecological Modelling 143 (3), 165–181. Borsuk, M.E., Stow, C.A., Luettich, R.A., Paerl, H.W., Pinckney, J.L., 2001b. Modeling oxygen dynamics in an intermittently stratified estuary: estimation of process rates using field data. Estuarine, Coastal and Shelf Science 52 (1), 33–49. Borsuk, M.E., Stow, C.A., Reckhow, K.H., 2004. A Bayesian network of eutrophication models for synthesis, prediction, and uncertainty analysis. Ecological Modelling 173 (2–3), 219–239. Cohen, A.C., 1955. Restriction and selection in samples from bivariate normal distributions. Journal of the American Statistical Association 50 (271), 884–893. Dilks, D.W., Canale, R.P., Meier, P.G., 1992. Development of Bayesian Monte-Carlo techniques for water-quality model uncertainty. Ecological Modelling 62 (1–3), 149–162. Dominici, F., Zeger, S., Samet, J., 2000. Combining evidence on air pollution and daily mortality from the largest 20 US cities: a hierarchical modeling strategy. Journal of the Royal Statistical Society: Series A 163, 263–302. Donigian, A.S., Huber, W.C., 1991. Modeling of nonpoint source water quality in urban and non-urban areas. http://www.epa. gov/waterscience/basins/docs/npsmodel.pdf Available at: Eleria, A., Vogel, R.M., 2005. Predicting fecal coliform bacteria levels in the Charles River, Massachusetts, USA. Journal of the American Water Resources Association 41 (5), 1195–1209. EPA, 2007. STORET. http://www.epa.gov/storet/dw_home.html. Fennessey, N., Vogel, R.M., 1990. Regional flow-duration curves for ungauged sites in Massachusetts. Journal of Water Resources Planning and Management-ASCE 116 (4), 530–549. Ferguson, C.M., Coote, B.G., Ashbolt, N.J., Stevenson, I.M., 1996. Relationships between indicators, pathogens and water quality in an estuarine system. Water Research 30 (9), 2045–2054. Francy, D.S., Helsel, D.R., Nally, R.A., 2000. Occurrence and distribution of microbiological indicators in groundwater and stream water. Water Environment Research 72 (2), 152–161. Freer, J., Beven, K., Ambroise, B., 1996. Bayesian estimation of uncertainty in runoff prediction and the value of data: an application of the GLUE approach. Water Resources Research 32 (7), 2161–2173.

1016

water research 44 (2010) 1006–1016

Gelman, A., Carlin, J.B., Stern, H.S., Rubin, D.B., 2004. Bayesian Data Analysis. Chapman and Hall/CRC, New York. Gelman, A., Rubin, D.B., 1992. Inference from iterative simulation using multiple sequences. Statistical Science 7 (4), 457–472. Ghosh, M., Rao, J.N.K., 1994. Small area estimation: an appraisal. Statistical Science 9 (1), 55–76. Goyal, A., Small, M., von Stackelberg, K., 2005. Estimation of fugitive lead emission rates from secondary lead facilities using hierarchical Bayesian models. Environmental Science and Technology 39, 4929–4937. Gronewold, A.D., Qian, S.S., Wolpert, R.L., Reckhow, K.H., 2009. Calibrating and validating bacterial water quality models: a Bayesian approach. Water Research 43 (10), 2688–2698. Harrison, K.W., 2007. Test application of Bayesian Programming: adaptive water quality management under uncertainty. Advances in Water Resources 30 (3), 606–622. Helser, T., Lai, H., 2004. A Bayesian hierarchical meta-analysis of fish growth: with an example for North American largemouth bass, Micropterus salmoides. Ecological Modeling 178, 399–416. Hobbs, B.F., 1997. Bayesian methods for analyzing climate change and water resource uncertainties. Journal of Environmental Management 49 (1), 53–72. Horn, A.L., Rueda, F.J., Hormann, G., Fohrer, N., 2004. Implementing river water quality modeling issues in mesoscale watershed models for water policy demands – an overview on current concepts, deficits, and future tasks. Physics and Chemistry of the Earth 29 (11–12), 725–737. Kanso, A., Gromaire, M.C., Gaume, E., Tassin, B., Chebbo, G., 2003. Bayesian approach for the calibration of models: application to an urban stormwater pollution model. Water Science and Technology 47 (4), 77–84. Liu, Y., Yang, P., Hu, C., Guo, H., 2008. Water quality modeling for load reduction under uncertainty: a Bayesian approach. Water Research 42 (13), 3305–3314. Malve, O., Qian, S.S., 2006. Estimating nutrients and chlorophyll a relationships in Finnish lakes. Environmental Science and Technology 40 (24), 7848–7853. Mas, D., Ahlfeld, D., 2007. Comparing artificial neural networks and regression models for predicting fecal coliform concentrations. Hydrological Sciences Journal 52 (4), 713–731.

Milly, P.C.D., Betancourt, J., Falkenmark, M., Hirsch, R.M., Kundzewicz, Z.W., Lettenmaier, D.P., Stouffer, R.J., 2008. Stationarity is DEAD: whither water management? Science 319 (5863), 573–574. Ott, W., 1995. Environmental Statistics and Data Analysis. Boca Raton. Qian, S.S., Reckhow, K.H., 2007. Combining model results and monitoring data for water quality assessment. Environmental Science and Technology 41 (14), 5008–5013. Qian, S.S., Reckhow, K.H., Zhai, J., McMahon, G., 2005. Nonlinear regression modeling of nutrient loads in streams: a Bayesian approach. Water Resources Research 41 (7). Qian, S.S., Stow, C.A., Borsuk, M.E., 2003. On Monte Carlo methods for Bayesian inference. Ecological Modelling 159 (2–3), 269–277. Riggs, H.C., 1989. Frequency Curves. United States Government Printing Office. Sohn, M.D., Small, M.J., Pantazidou, M., 2000. Reducing uncertainty in site characterization using Bayes Monte Carlo methods. Journal of Environmental Engineering-ASCE 126 (10), 893–902. Spiegelhalter, D.J., Thomas, A., Best, N.G., 1999. WinBUGS Version 1.2 User Manual. MRC Biostatistics Unit. Stiber, N.A., Pantazidou, M., Small, M.J., 1999. Expert system methodology for evaluating reductive dechlorination at TCE sites. Environmental Science and Technology 33 (17), 3012–3020. Stiber, N.A., Pantazidou, M., Small, M.J., 2004. Embedding expert knowledge in a decision model: evaluating natural attenuation at TCE sites. Journal of Hazardous Materials 110 (1–3), 151–160. Varis, O., 1998. A belief network approach to optimization and parameter estimation: application to resource and environmental management. Artificial Intelligence 101 (1–2), 135–163. Vogel, R.M., Rudolph, B.E., Hooper, R.P., 2005. Probabilistic behavior of water-quality loads. Journal of Environmental Engineering-ASCE 131 (7), 1081–1089. Wikle, C.K., 2003. Hierarchical models in environmental science. International Statistical Review 71 (2), 181–199. Wilkinson, J., Jenkins, A., Wyer, M., Kay, D., 1995. Modelling Feacal Coliform Concentrations in Streams. Institute of Hydrology, Oxfordshire, U.K.