Assessment of uncertainty sources in water quality modeling in the Niagara River

Advances in Water Resources 33 (2010) 493–503

Contents lists available at ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

Assessment of uncertainty sources in water quality modeling in the Niagara River Samuela Franceschini a, Christina W. Tsai b,* a b

Via Loredan 20, Dipartimento di Ingegneria Idraulica, Marittima e Geotecnica, Università di Padova, Padua, Italy 233 Jarvis Hall, Civil, Structural and Environmental Engineering, State University of New York at Buffalo, Buffalo, NY 14260, United States

a r t i c l e

i n f o

Article history: Received 1 February 2009 Received in revised form 2 July 2009 Accepted 1 February 2010 Available online 4 February 2010 Keywords: Uncertainty analysis Water-quality modeling Surface flows Natural rivers

a b s t r a c t This paper presents a framework to quantify the overall variability of the model estimations of Total Polychlorinated Biphenyls (Total PCBs) concentrations in the Niagara River on the basis of the uncertainty of few model parameters and the natural variability embedded in some of the model input variables. The results of the uncertainty analysis are used to understand the importance of stochastic model components and their effect on the overall reliability of the model output and to evaluate multiple sources of uncertainty that might need to be further studied. The uncertainty analysis is performed using a newly developed point estimate method, the Modified Rosenblueth method. The water quality along the Niagara River is simulated by coupling two numerical models the Environmental Fluid Dynamic Code (EFDC) – for the hydrodynamic portion of the study and the Water Quality Analysis and Simulation Program (WASP) – for the fate and transport of contaminants. For the monitoring period from May 1995 to March 1997, the inflow Total PCBs concentration from Lake Erie is the stochastic component that most influences the variability of the modeling results for the simulated concentrations at the exit of the Niagara River. Other significant stochastic components in order are as follows: the suspended sediments concentration, the point source loadings and to a minor degree the atmospheric deposition, the flow and the non-point source loadings. Model results that include estimates of uncertainty provide more comprehensive information about the variability of contaminant concentrations, such as confidence intervals, and, in general offer a better approach to compare model results with measured data. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction An enhanced understanding and better characterization of the transport of contaminants in surface waters has fundamental significance to environmental engineers. Mathematical modeling of such natural phenomena is subject to uncertainty due to inaccurate or oversimplifying representations of the processes in the model, measurement errors, limited availability of data, inherent spatial and temporal randomness of variables and model parameters uncertainty among others. The ability to incorporate the aforementioned uncertainty sources into the modeling of contaminant concentrations is crucial to quantifying the effective risks of exceeding a designated water-quality standard. The variability in data and the uncertainty associated with model parameters in water-quality modeling cannot be avoided. Several techniques are available in uncertainty analysis to account for uncertainty sources in the modeling process. Tung and Yen [26] classified the sources of uncertainty into natural variability and model related uncertainty. The inherent variability of natural sys* Corresponding author. Tel.: +1 7166452114. E-mail addresses: [email protected] (S. Franceschini), ctsai4@eng. buffalo.edu, [email protected] (C.W. Tsai). 0309-1708/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2010.02.001

tems is typical of random phenomena and is intrinsic in the data series used in the model. This type of variability is always present and cannot generally be controlled. The natural variability represents the randomness of geophysical processes such as the frequency and the magnitude of a hydrological event or of spatiotemporal variation of physical variables involved in the transport processes such as in the case of climate conditions. The other source of fluctuations in the model results comes from the uncertainty specific to the mathematical model used to represent the natural system and is often called ‘‘knowledge deficiency” [26]. Such uncertainty may arise from the incomplete understanding of the system being modeled and/or the inability to accurately reproduce the studied processes with mathematical and statistical techniques [11]. The model related uncertainty can be further classified in three main subgroups: model conceptualization uncertainty (the uncertainty deriving from the model formulation, the parameter calibration, and the numerical error), data related uncertainty (the uncertainty due to measurement errors, temporal and spatial sampling resolution, and sampling and data handling errors), and operational uncertainty (the uncertainty linked to the management and operation of the system under study). The main advantage of uncertainty analysis is the ability to statistically quantify the response of the system to input variables and

494

S. Franceschini, C.W. Tsai / Advances in Water Resources 33 (2010) 493–503

parameters that are random or uncertain. As opposed to traditional sensitivity analysis, in which the magnitude of a perturbation to the selected variables is generally arbitrarily chosen, in uncertainty analysis, the selected model variables and parameters are treated as random variables. Thus, the variability of random variables is measured in terms of probability density functions or statistical moments. Furthermore, in uncertainty analysis, the effect of the selected random variables can be evaluated individually or collectively, while in traditional perturbation methods, the sensitivity of the model output to a selected variable is assessed one variable at a time [19]. Uncertainty analysis techniques and their applications have continuously evolved since early attempts to incorporate uncertainty computation in water-quality modeling. Ever since, simulation methods and first-order-second-moment methods have found application in a variety of studies primarily focusing on the evaluation of parameter uncertainty related to the dissolved oxygen concentration and biochemical oxygen demand (e.g. [2,3,10,15,16,19,27]). Jaffe and Ferrara [13] and Yulianti et al. [32] also applied uncertainty analysis to model the behavior of contaminated sediments in the water column. More recently, Preston and Jones [21] assessed the sensitivity of the risk of runoff changes to various sources of uncertainty, including global mean temperature, using uncertainty analysis techniques. In addition, uncertainty analysis has been used as a supporting tool in decision making as discussed by Harrison [12] and Franceschini and Tsai [7] in two applications to water pollution control. In the last two decades, Point Estimate Methods (PEMs) have gained momentum as valuable easy-to-apply techniques that do not require much information on the probability distribution of the stochastic variables. Ünlü [28] discussed the results of the application of one PEM compared to Monte Carlo simulation and First Order Second Moment method in groundwater modeling. Yu et al. [31] applied PEMs to rainfall-runoff studies, Lian and Yen [14] to culvert analysis and Franceschini and Tsai [7] to water-quality management. Tsai and Franceschini [25] have also introduced an improved PEM modified from Rosenblueth [22] to evaluate uncertainty in the computation of settling velocity of particles in a fluid. This study focuses on quantifying the variability of model outputs (contaminant concentrations) due to data variability and model parameter uncertainty. Model errors (or model formulation uncertainty), interpreted as the uncertainty in both the formulation of the problem and the chosen solution method, was not considered in this study. The work presented in this paper is an attempt at ascertaining the degree of influence of several variables and parameters used in modeling Total PCBs concentration in the Niagara River, as summarized in Table 1. The results obtained by this uncertainty analysis can be used to assess the need for further characterization of contaminant transport processes occurring in the river. Ranking the effects of inherently random model variables and parameters is the first step to optimize and improve the modeling effort. The following sections explore in detail the sources of Table 1 Model kinetics and parameters used in the water quality model WASP. Parameter

Total PCBs run

Wind Speed, W Reareation rate coefficient, K a Initial average fraction dissolved Octanol–water partition coefficient, K ow Fraction of organic carbon in solids, focs Henry’s constant, H Molecular weight, Mw Atmospheric concentration, C a Biodegradation rate coefficient, K

6.4 m/s 4.72 I/day 91% 6.101og(Lw/Lo) 0.1 4.90E3 atm-m3/mole 292 g/mole 1E13 g/m3 3.02E3 1/day

data and the expected effects of uncertainty of each stochastic model component used in this work. 2. Study site: the Niagara River The Niagara River is a natural channel between Lake Erie and Lake Ontario and serves as a border between Ontario, Canada and the western region of New York State, USA. The location and its hydraulic features (e.g, Niagara Falls, water diversions for hydropower production) make the Niagara River an important body of water and a significant resource for the western New York region and its economy. For decades, the communities and the industries along its banks have used this natural system as a drinking water supply and an inexpensive source of electric energy. However, past lax waste management practices have left behind soil and groundwater contamination. The incoming concentrations from Lake Erie, polluted seepage and point discharges from wastewater treatment plants and industrial facilities are the major sources of contamination to the Niagara River. In February 1987 Environment Canada, U.S. Environmental Protection Agency (EPA) Region II, the Ontario Ministry of Environment and the New York State Department of Environmental Conservation (NYSDEC) – established the Niagara River Toxic Management Plan (NRTMP) with the intent of reducing toxic chemical inputs to the river, and achieving water-quality standards that are protective of human health and aquatic life. As part of this program, samples are collected at the head (Fort Erie, Canada) and at the mouth (Niagara on the Lake, Canada) of the river on a weekly/bi-weekly basis and analyzed for a suite of 18 priority toxic pollutants. Total Polychlorinated Biphenyls (Total PCBs) is the contaminant modeled in this study. 3. Model development This work proposes a more comprehensive and realistic representation of the fate and transport of contaminants in the Niagara River by coupling deterministic hydrodynamic and water-quality modeling with a probabilistic analysis. It also offers a more rigorous approach to modeling the hydrodynamic and physico-chemical phenomena by coupling a hydrodynamic model and a waterquality model as described in the following sections. 3.1. Hydrodynamic modeling – Environmental Fluid Dynamic Code (EFDC) The hydrodynamic model Environmental Fluid Dynamics Code (EFDC) was used in this study to simulate the movement of water along the Niagara River. The EFDC model was initially developed by Hamrick [9]. The model has been applied to the study of several estuarine systems including the Chesapeake Bay, and to simulations of pathogens, pollutants, and crustaceous larvae fate and transport [24]. The EFDC is configured to solve the three-dimensional, free surface, turbulence averaged equations of motion for a fluid. In this study, the model was run for a four-year monitoring period, starting on April 1, 1993 through March 31, 1997 with a lag period of 30 days between the beginning of the EFDC simulation and the beginning of the water-quality simulation. In other words, EFDC results for the water-quality model were available from May 1, 1993. A 2year data set between April 1993 and March 1995 was used for calibration, while a second 2-year set was used for verification and successively in the uncertainty analysis. The model grid extended from the exit of Lake Erie at Fort Erie up to the exit of the river into Lake Ontario at Niagara on the Lake, as shown in Fig. 1. In total, 130 active cells of varying dimensions

S. Franceschini, C.W. Tsai / Advances in Water Resources 33 (2010) 493–503

495

Fig. 1. EFDC grid super imposed to a satellite picture of the Niagara River (courtesy of NASA Earth Observatory).

followed the general shape of the river. Each cell was vertically composed of two layers. The incoming boundary condition consisted of the incoming flow from Lake Erie at Buffalo as reported by the US Army Corps of Engineers – Detroit District (USACE). At the mouth of the river, we defined an open boundary condition as the standard specification of water surface elevation using combinations of harmonic constituents and time series. The initial depth assigned to each cell represented the initial condition. Flow sinks and sources on the top layers of the grid represented the withdrawal and return flows for hydropower production at the reservoirs indicated in Fig. 1. A control structure based on a weir assumption represented Niagara Falls. The main EFDC modeling coefficients were bottom roughness, bathymetry, and vertical mixing [23]. The bottom roughness coefficient was calibrated against the flow reported at the mouth of the river (Niagara on the Lake, Canada) by the USACE.

specify time-variable exchange coefficients, advective flow loadings, water quality boundary conditions, and other parameters. Vuksanovic et al. [29] indicated that the model is capable of simulating evolution profiles of dissolved and sorbed PCB. The WASP model allows for one-, two- and three-dimensional representation based on a box-modeling approach [30]. Model kinetics and parameters determined the decay and transport of the pollutants in the system. The WASP model considered the processes of volatilization, sorption, and water column biodegradation along with the transport of dissolved and particulate components. The kinetic rates and parameters used to model Total PCBs are summarized in Table 1. Parameters for pollutant decay and sorption were calibrated with respect to the contaminant concentrations measured at Niagara on the Lake, Canada as reported by the NRTMP. The decay rate had very little impact on the overall concentration. The initial and upstream boundary conditions were

3.2. Contaminant fate and transport – Water Quality Analysis Simulation Program (WASP)

Cðx; 0Þ ¼ C 0

The numerical model used for simulating the water quality of the Niagara River is WASP version 7.1. The advantages of using WASP are the flexibility offered by the model and the ability to

and Cð0; tÞ ¼ C in ðtÞ;

where C :0 and C :in: were the upstream measurements for contaminant concentration obtained through the NRTMP. Although the input concentrations were available as monthly mean values, a continuous source was assumed to represent the input from Lake

S. Franceschini, C.W. Tsai / Advances in Water Resources 33 (2010) 493–503

Erie in the WASP model. The numerical model interpolated between incoming concentration values. In the grid representation of the river, the last cell corresponding to the exit boundary from the Niagara River did not coincide with the physical end of the river but was moved out into Lake Ontario. Therefore, the concentration at the downstream boundary was set to zero to assure the presence of a gradient, but at the same time to avoid bringing any influence to the outflow concentration at the physical exit of the river. Incoming boundary concentrations were also applied to the return flows from the reservoirs of the two hydropower plants. These incoming concentrations were computed from the withdrawal concentrations, based on the dimensions of each reservoir, assuming steady state and completely mixed conditions.

Model

Observed

8.000 7.500

Flow (cms)

496

7.000 6.500 6.000 5.500 5.000 0

50

100

150

200

250

300

350

April 1996 - March 1997

3.3. Modified Rosenblueth method

Fig. 2. Comparison of reported flows at Queenstown, NY and model out flows (portion of the validation run).

   g  N m XN  m EðY Þ ¼ ð1  NÞ þ Y þ ðpiþ  gi þ 1ÞY m iþ þ pi Y i i¼1 2 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} variability contributed by m

central value

þ

XN1 XN

Ym ij gij |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} i¼1

j¼iþ1

each random variable X :i:

ð1Þ

variability contributed by the correlated pair ðX :i: ;X :j: Þ

where N is the number of stochastic variables and parameters conP P sidered in the uncertainty analysis, g ¼ i gi , where gi ¼ j gij and 0 0 gij ¼ qij =xiþ xjþ with the exception of gii which by definition is equal to one. The variable Y represents the evaluation of the model function Y ¼ f ðX :1: ; X :2: ; . . . ; X :N: Þ at the mean value of each random variable. The variables Y :iþ: and Y :i: are calculated from the model function, using x:iþ: and x:i: respectively for the variable X :i: and mean values for the other variables X :j: with j–i. The points x:iþ: and x:i: correspond to the discrete points at which the probability mass is concentrated in the two-point representation of the pdf of the variable X :i , while p:iþ: and p:i: are the associated weights. In the last term of Eq. (1), Y :ij: is expressed as Y ij ¼ f ðl1 ; . . . ; li1 ; xiþ ; liþ1 ; . . . ; lj1 ; xjþ ; ljþ1 ; . . . ; lN Þ, which represents the model function evaluated at the mean value of all the variables X :k: , except for k = i and k = j, for which the values of x:iþ: and x:jþ: are used. Details on the derivation of Eq. (1) based on the Taylor’s Series expansion can be found in Tsai and Franceschini [25].

Observed

Calibration

Verification

0.0025

Total PCBs concentration ( µ g/L)

In the uncertainty analysis, the Modified Rosenblueth method was used along with multiple runs of the deterministic model. The Modified Rosenblueth method is a Point Estimate Method based on a two-point representation of the probability density function (pdf) of the stochastic variables, X :i . The mth statistical moment of the dependent variable Y (Eq. (1)) is composed of three terms: a first term based on the mean values of the variables X i , a second term representing the variability due to the pdf of each random variable, and a final term based on the variability due to correlation between pairs of random variables

0.002

0.0015

0.001

0.0005

0 Apr-93

Oct-93

Apr-94

Oct-94

Apr-95

Oct-95

Apr-96

Oct-96

Apr-97

Monitoring Period Fig. 3. Measured concentrations of Total PCBs at the end of the Niagara River and results of the calibration and verification runs.

reported flow values at the exit of the Niagara River, thus indicating the functionality of the model. In the calibration of the water quality phase of the model, some modeling judgment was applied in the utilization of certain loadings (particularly non-point sources) since some reported measurements were not consistent with other data sets and, which may have resulted in over-prediction of concentrations at the end of the Niagara River. The calibration of the deterministic model was performed on a two year data set from May 1993 to March 1995 (Fig. 3). Simulated flows and concentrations at the end of the Niagara River were graphically compared with the respective measured values at the same location. For the verification process, the time series for flow, inflow concentration, suspended sediments, and air and water temperatures for the monitoring period from April 1995 to the end of March 1997 were used.

3.4. Calibration and validation 4. Stochastic components For the calibration of the hydrodynamic model, the optimization software OSTRICH [17] was applied. This optimization software has been created as a model-independent program that allows automation of the processes of model calibration and design optimization. The OSTRICH package includes several algorithms that can be used for weighted least squares calibration. For this study, the Genetic Algorithm was chosen to calibrate the global bottom roughness and the cell related bottom roughness corrections. As shown in Fig. 2, there is a fairly good correspondence between the hydrodynamic modeling results for flow and the

The influences of eight variables (flow, inflow concentrations, suspended sediment concentration, non-point and point source loadings, atmospheric deposition, air temperature, and water temperature) and two model parameters (roughness coefficient and degradation rate) on the variability of estimated concentrations at the end of the Niagara River were investigated. Some details about the stochastic components used in the uncertainty analysis are summarized in Table 2 and discussed in the following paragraphs.

497

S. Franceschini, C.W. Tsai / Advances in Water Resources 33 (2010) 493–503 Table 2 List of random variables and parameters used in this study. Data

Type

Flow Inflow concentrations Point and non-point source loadings Atmospheric deposition Incoming suspended sediment concentration Air temperature Water temperature Degradation rate Roughness coef. a b c d e f g

Variable Variable Variable Variable Variable Variable Variable Parameter Parameter

Source a

USGS Station NRTMPb IHWSc NYSDECd IADNe NRTMP NOAAf DWTPg Calibration Calibration

Collection point

Collection frequency

Buffalo, NY Fort Erie, Canada Along River Lake Erie and Lake Ontario Fort Erie, Canada Buffalo Station Buffalo N/A N/A

Mean Daily Flow Weekly/Bi-weekly Monthly Monthly Weekly/Bi-weekly Monthly Monthly N/A N/A

US geological survey. Niagara River toxic management plan. Inactive hazardous waste sites. New York State Dept. of Environmental Conservation. Integrated atmospheric deposition network. National oceanic atmospheric administration. Drinking water treatment plant.

Limited data were available on the statistical correlation between stochastic components used in this study. Some correlation was expected to exist among model components such as flow and inflow concentration and air and water temperature. However, only the correlation between flow and inflow contaminant concentration was considered. This value was computed from the available data and equaled 0.10. All other correlation coefficients were set to zero in this paper for simplicity. In general, the PEMs and specifically the Modified Rosenblueth method used in this study are less computationally demanding techniques than other uncertainty analysis methods. The number of times the deterministic model is run to evaluate the uncertainty contributed by the stochastic components is controlled by Eq. (1). In this study, for ten stochastic variables and parameters, the deterministic model was run only 22 times. The model function was evaluated once at the mean values and 2N times at the x:þ: and x:: values of each of the considered stochastic variables, where N is the number of stochastic variables. When the correlation is accounted for, the model function is also evaluated ðN :2:  NÞ=2 times, where N, in this case, is the number of correlated variables (which is two in this study). The first run, referred to as the base run, evaluated the model function at the mean values of all uncertain variables and parameters used in the model. Each successive run included one additional stochastic component, as indicated in Table 3. The results of each run were then compared graphically to evaluate the degree of influence of the selected stochastic variables and parameters. The cumulative effects were also evaluated. For the purpose of establishing the order of importance of each stochastic variable

or parameter, the order in which each stochastic component was added was also varied. The figures presented in this paper refer to the order described in Table 3.

4.1. Flow USACE provided the data for flow. The Corps is in charge of computing daily mean Lake Erie outflows and daily mean Lake Ontario inflows from an empirical equation that uses stage data from a gauge station operated by the US Geological Survey in the City of Niagara Falls. The gauge station is located just downstream of the pool that forms below the Niagara Falls. The uncertainty related to the flow included both the natural variability and knowledge deficiency. The latter uncertainty derived from the fact USACE computes both the inflows and the outflows using an empirical formula that relates flow stage to flow rate rather than using on site flow measurements. The statistical properties for the flow distribution were obtained from the data. The mean and variance of the mean daily flow for the selected month were computed according to the following equations.

Zm ¼

n X

Z i =n;

ð2Þ

i¼1

Var½Z m  ¼ Var½Z i =n2 ;

ð3Þ

where Z :m: is the mean daily flow for the selected month; Z :i: is the daily mean flows reported for the selected month and n is the number of days in the month.

Table 3 Order of inclusion of each stochastic model components in the uncertainty analysis. Uncertainty analysis runs Base Run QRun Cin Run SSRun Nps Run Ps Run AD Run A Run W Run Kd Run All Run

Mean values of all the stochastic model components Include one stochastic variable: flow Include two stochastic variables: flow and inflow concentration Include three stochastic variables: flow, inflow concentration and suspended solids Include four stochastic variables: flow, inflow concentration and suspended solids and non-point source loadings Include five stochastic variables: flow, inflow concentration and suspended solids and non-point and point source loadings Include six stochastic variables: flow, inflow concentration and suspended solids, non-point and point source loadings and atmospheric deposition loadings Include seven stochastic variables: flow, inflow concentration and suspended solids, non-point and point source loadings atmospheric deposition loadings and air temperature Include eight stochastic variables: flow, inflow concentration and suspended solids, non-point and point source loadings atmospheric deposition loadings, air and water temperature Include nine stochastic variables: flow, inflow concentration and suspended solids, non-point and point source loadings atmospheric deposition loadings, air and water temperature, biodegradation rate Include all 10 stochastic variables: flow, inflow concentration and suspended solids, non-point and point source loadings atmospheric deposition loadings, air and water temperature, degradation rate and global roughness

498

S. Franceschini, C.W. Tsai / Advances in Water Resources 33 (2010) 493–503

The variability of the flow data was expected to have only a marginal influence on the overall simulated concentration. In general, larger flow would have a more significant effect of flushing the river of both contaminant and suspended sediments with the consequent result of lowering the overall contaminant concentration in the river.

4.4. Loadings Contaminants enter the Niagara River both via point and nonpoint sources. Several contaminant loading sources such as industrial and waste water treatment plant discharges, tributary inflows, surface runoff, contaminated groundwater seepage from inactive hazardous waste sites (IHWS) and subsurface sediments contribute to the pollution of the Niagara River.

4.2. Inflow concentration 4.5. Point sources The Environment Canada manages the Upstream/Downstream monitoring program established with the NRTMP. Through this program, samples are collected and analyzed for both dissolved and suspended sediment adsorbed phases and suspended sediment concentrations on a weekly or bi-weekly basis. This study modeled Total PCBs. Statistical data and distribution properties of the inflow concentration were obtained from the data, based on the values recorded for the chosen monitoring period. The mean inflow concentration and the associated standard deviation for each month in the monitoring period were computed from the data provided through the NRTMP. It was decided not to compute the statistics for any given month from multiple-year datasets, since comparison with recorded concentration values at the end of the Niagara River would have lost significance. All the mean inflow concentrations for the selected month were computed according to the Eqs. (2) and (3), where Z :m: in this case is the mean inflow concentration for the selected month; Z :i: is the inflow concentrations recorded for the selected month and n is the number of reported measurements. The fluctuation in the inflow concentrations included both natural variability and data uncertainty. The latter is related to possible measurement and handling errors. In general, it was expected that the uncertainty of the overall results was largely attributed to the variability of the inflow concentrations. In particular, it was proposed that increases in the inflow concentration induced increases in the simulated concentration at the end of the river.

4.3. Suspended sediment concentration Suspended sediment concentrations are measured through the NRTMP. Samples for suspended sediments are collected by centrifugation on a weekly, or bi-weekly, basis at Fort Erie and Niagara on the Lake over the period 1980–1997. A set of the suspended sediment concentrations was used in the water-quality modeling portion of this study. Statistical properties for suspended sediment concentrations were obtained from the data. The mean suspended concentrations for each month were computed according to Eqs. (2) and (3) where Z :m: is the mean suspended concentration for the selected month, Z :i: is the suspended sediment concentrations recorded for the selected month and n is the number of reported measurements. As for the inflow concentrations, the fluctuation of measured concentrations of suspended sediments combined both natural variability and data handling errors. The variability in the measured suspended sediment concentration influenced the overall contaminant concentration in the river. In fact, it was hypothesized that the fluctuation in the concentration of suspended sediments altered the ratio of suspended to dissolved contaminant concentration. The greater the suspended sediment concentration was, the higher the portion of contaminant in the sorbed phase was. Therefore, due to settling and exiting of the sediments at the end of the Niagara River, the overall concentration of contaminant in the water column was likely to decrease.

In 1981 and 1982 the NYSDEC conducted a comprehensive assessment of the municipal and industrial waste water discharges to the Niagara River. It was established that these discharges contributed up to 95% of all toxic discharges to the Niagara River [6]. The majority of the data for point sources are obtained as discharges from waste water treatment plants and industrial facilities. Other point sources are represented by the tributaries to the Niagara River. Loadings from tributaries are not monitored directly and are not included in this study 4.6. Non-point sources Two reports by Gradient Corporation GeoTrans, Inc. and by TRC Environmental Corporation for USEPA have provided a detailed evaluation of the potential chemical input due to seepage of contaminated groundwater from US IHWS [4]. Estimates of non-point source loadings from IHWS are computed based on averaged concentrations and groundwater flows [6]. Precise and updated estimates of the loadings effectively entering the Niagara River are not available. The TRC study reports only six of the 18 chemicals analyzed by the NRTMP. Data for these six contaminants were used in the past to model their concentrations in the Niagara River and they were used, unaltered, in this study. River bottom sediments are also considered to produce nonpoint loadings. However, due to the rapid currents, the main channel of the Niagara River does not contain substantial deposits of fine-grained cohesive sediments. Other non-point sources include surface runoff. These sources of non-point loadings were not considered in this study due to their limited influence on the overall concentration and the lack of substantial data for the modeling effort. The uncertainty of the loadings data was mostly attributed to knowledge deficiency, in particular the inability to confirm the reported values. However, the randomness of natural phenomena was also part of the variability of the data. The influence of the point and non-point sources on the model estimations was expected to be minimal. This was due to the fact that the magnitude of these loadings was much smaller compared to that of the inflow concentration, which is the factor that most influences the model results. Due to the limited data availability for the monitoring period and the time scale considered in this study, the statistical properties for the loadings were derived by assuming a log-normal distribution and a coefficient of variation (COV) of 25% for the point source loadings and 60% for the non-point source loadings [5,20]. The large COV for non-point source loadings was considered to be representative of the estimated nature of the available data. Non-point source measurements have not been updated in the last two decades. 4.7. Atmospheric deposition The Integrated Atmospheric Deposition Network (IADN) is a binational program coordinated by the Air Quality Research Branch of Environment Canada and the Great Lakes National Program

499

S. Franceschini, C.W. Tsai / Advances in Water Resources 33 (2010) 493–503

The air and water temperatures were taken into account in this study for their importance in the degradation and volatilization processes. The water temperature data are gathered from the Buffalo Water Authority that manages the drinking water treatment plant for the City of Buffalo and from the USACE, Buffalo District. The Buffalo Water Authority intake for the treatment plant is located at the head of the Niagara River. The air temperatures were obtained from two weather stations: one at the Buffalo airport, managed by the National Oceanic Atmospheric Administration (NOAA) and the other in Burlington, Ontario, Canada, managed by Environment Canada. In this study, the temperatures from these two stations were averaged. The uncertainty of the air and water temperatures was attributed to the inherent natural variation and to possible measurement errors. The statistical properties of the temperature distributions were obtained from the data. The mean and variance of the mean daily temperature for the selected months are computed using Eqs. (2) and (3) where Z :m: is the mean daily air or water temperature; Z :i: is the daily mean temperatures recorded for the selected month and n is the number of days in the month. Although the overall effect of the variability of air and water temperature was anticipated to be minimal, these sources of uncertainty were still considered due to their importance in some physico-chemical processes that contaminants undergo while in the water. In particular, in the case of volatilization, it was expected that winter temperature strongly diminished volatilization and consequently reduced this sink of contaminants in the Niagara River. 4.9. Degradation rate The biodegradation rate was obtained through calibration of the water quality portion of the deterministic model (WASP). The biodegradation rate was used as a stochastic component in the uncertainty analysis. It was not expected to have a large influence on the overall results, primarily because Total PCBs are persistent pollutants. A uniform distribution was assumed with minimum and maximum values obtained from the calibration and statistical moments computed from the distribution.

The bottom roughness coefficient was calibrated against the flow reported at the mouth of the river (Niagara on the Lake, Canada) by USACE. The global bottom roughness for the Niagara River was determined to be 0.0137 in the calibration. This value was taken as the mean value of bottom roughness. Statistical data were then obtained from the uniform distribution assuming a COV of 25% [8]. 5. Results Presented in Figs. 4 and 5 are two comparisons of the observed Total PCBs concentrations recorded at Niagara on the Lake with  the simulated mean value of Total PCBs (Fig. 4) modeled using all mean values for the considered stochastic components (as in a deterministic approach) along with the expected value of Y (E[Y], Total PCBs) computed with Eq. (1) for m = 1  the one-standard deviation interval about the expected value E[Y] obtained from the uncertainty analysis based on the considered stochastic model components (Fig. 5). Scatter plots of the observed concentrations of Total PCBs versus the deterministic model output (Mean) and the uncertainty model output (E[Y]) and their respective linear regressions are presented in Fig. 6.

Observed

Deterministic

Uncertainty

0.0025

0.002

Total PCBs (mg/L)

4.8. Air and water temperatures

4.10. Roughness coefficient

0.0015

0.001

0.0005

0 May-95

Jul-95

Oct-95

Jan-96

Apr-96

Jul-96

Oct-96

Jan-97

Monitoring Period

Fig. 4. Measured concentrations of Total PCBs at the end of the Niagara River, results of the ‘‘base run” using mean values for all the uncertainty components and of the expected value E[Y] compute using Eq. (1) for m = 1.

Observed

All Upper

All Lower

0.0025

0.002

Total PCBs (mg/L)

Office of the USEPA. The IADN was initiated in 1990 to measure and report atmospheric concentrations of persistent toxic pollutants in the Great Lakes basin [1]. The data available through the IADN program were used in this study. In particular, for the Niagara River net deposition values were computed as the average of net deposition for Lake Erie and Lake Ontario (which was found to be negative for the period considered in this study indicating volatilization of the compound) and assigned to each segment of the model based on its surface area. Volatilization was also accounted for in the water-quality model WASP. Both forms of volatilization were considered, since overall, the model predictions tended to slightly overestimate the observed data at the end of the river. In addition to the WASP simulated process of volatilization, the model of mass transfer from aqueous to gaseous phases at waterfalls by McLachian et al. [18] was also applied to the cells representing the Niagara Falls. The uncertainty of the air deposition data included the inherently random nature of the processes considered here and the possible measurement and sample handling errors. For the uncertainty part of the study, the COV reported by the IADN program as the ratio between the standard error of detection over mean was used as a statistical property and a log-normal distribution was assumed.

0.0015

0.001

0.0005

0 May-95

Jul-95

Oct-95

Jan-96

Apr-96

Jul-96

Oct-96

Jan-97

Monitoring Period Fig. 5. Measured concentrations of Total PCBs at the end of the Niagara River and a one-standard deviation interval about the expected value E[Y].

500

S. Franceschini, C.W. Tsai / Advances in Water Resources 33 (2010) 493–503 Q

Cin

SS, Nps

Ps

AD, Air, Water, Kd = All

Observed

0.00135

Cin upper

All upper

All lower

Cin lower

0.0025

0.00125

Total PCBs (mg/L)

Mean Total PCBs (mg/L)

0.002 0.00115 0.00105 0.00095 0.00085

0.0015

0.001

0.0005

0.00075 0.00065 0.00055 May-95

0 May-95 Jul-95

Oct-95

Jan-96

Apr-96

Jul-96

Oct-96

Monitoring Period

Fig. 6. Expected value of Total PCBs mean values obtained by adding one uncertainty component at a time.

Mean

E[Y]

Oct-95

Jan-96

Apr-96

Jul-96

Oct-96

Jan-97

Monitoring Period

Fig. 9. Comparison between a one-standard deviation interval about the expected value for flow and inflow concentration uncertainty and a one-standard deviation interval about the expected value for all the uncertainty components.

against the overall one-standard deviation interval computed for the run that includes all the analyzed stochastic components (Figs. 8–11). As for the expected values, the one-standard deviation

Linear (E[Y])

0.0016 0.0014 0.0012

Observed

y = 0.6754x + 0.0002

SS upper

All upper

SS lower

All lower

0.0025

0.001 0.0008

0.002

0.0006 Total PCBs (mg/L)

Observed Concentration (mg/L)

Jul-95

Jan-97

0.0004 0.0002

y = 0.5319x + 0.0003

0 0

0.0002

0.0004

0.0006

0.0008

0.001

0.0012

0.0014

Predicted Deterministic Concentration (mg/L)

0.0015

0.001

0.0005

Fig. 7. Scatter plots of measured concentrations of Total PCBs at the end of the Niagara River, results of the ‘‘base run” using mean values for all the uncertainty components and of the expected value E[Y].

0 May-95

Jul-95

Oct-95

Jan-96

Apr-96

Jul-96

Oct-96

Jan-97

Monitoring Period

Shown in Fig. 7 are the expected values of Total PCBs computed from the uncertainty analysis for each run. The curves for fourth and fifth runs (SS and Nps respectively) are essentially identical. The curves for seventh through tenth runs (AD, Air, Water and Kd) are essentially indistinguishable from the last run that includes all variables (All). The one-standard deviation intervals about the expected value computed for several successive runs are plotted

Fig. 10. Comparison between a one-standard deviation interval about the expected value for flow, inflow concentration and suspended solids uncertainty and a onestandard deviation interval about the expected value for all the uncertainty components.

Observed NPS lower

NPS upper All lower

All upper PS lower

PS upper

0.0025

Observed

Q upper

All upper

Q lower

All lower

0.0025

Total PCBs (mg/L)

0.002

Total PCBs (mg/L)

0.002

0.0015

0.001

0.001

0.0005

0.0005

0 May-95

0.0015

0 May-95

Jul-95

Oct-95

Jan-96

Apr-96

Jul-96

Oct-96

Jan-97

Monitoring Period Jul-95

Oct-95

Jan-96

Apr-96

Jul-96

Oct-96

Jan-97

Monitoring Period

Fig. 8. Comparison between a one-standard deviation interval about the expected value for flow uncertainty only and a one-standard deviation interval about the expected value for all the uncertainty components.

Fig. 11. Comparison between two one-standard deviation intervals about the expected value for flow, inflow concentration and suspended solids uncertainty, one including non-point sources uncertainty only and one including non-point source and point sources uncertainty and a one-standard deviation interval about the expected value for all the uncertainty components.

S. Franceschini, C.W. Tsai / Advances in Water Resources 33 (2010) 493–503

501

Table 4 Order of importance of the stochastic components used in the uncertainty analysis. Order

Stochastic components

Portion of total one-standard deviation interval

1 2 3 4 N/A

Inflow Concentration, Cin Suspended Solids, SS Point Source loadings, PS Flow, Q atmospheric deposition, AD non-point source loadings, NPS Air temperature, air water temperature, water dispersion coefficient, Dx degradation rate, Kd

35–91% 5–40% 2–37% 0–8% No significant contribution

intervals depend on the number and the characteristics of the stochastic variables and parameters considered in the uncertainty analysis. As shown in Fig. 7, the magnitude of the expected values obtained from the uncertainty analysis is altered each time another stochastic component is added, as more uncertainty information is processed. Since the expected value curve of the Total PCBs from the uncertainty analysis is not fixed but may be higher or lower depending on the stochastic components being considered in the uncertainty analysis, the corresponding one-standard deviation interval might overlap or extend outside the computed one-standard deviation interval estimated for the ‘‘All run.” Furthermore, some one-standard deviation intervals are narrower or wider than the interval corresponding to the previous runs with less stochastic components (e.g. one-standard deviation interval for ‘‘SS run” – Fig. 10 and one-standard deviation interval for ‘‘Nps run” – Fig. 11 are narrower and one-standard deviation interval for ‘‘PS” – Fig. 11 is wider). This is due to the particular characteristic of the added stochastic element. For example, in the case of the ‘‘SS run,” the higher suspended sediment concentrations correspond to lower contaminant concentrations at the end of the river since sorbed contaminants are flushed out of the river. The effect of this process on the plotted data is a slightly narrower one-standard deviation interval. However, when evaluated collectively, the single effect of the stochastic components is compensated or smoothed out in the ‘‘All run” one-standard deviation interval. Reported in Table 4 is the order of influence of the stochastic variables and components used in the uncertainty analysis and the associated percentage range of the overall one-standard deviation interval contributed.

6. Discussion Overall, the deterministic representation of the fate and transport of contaminants in the Niagara River captures the physical trend of the measured concentrations of Total PCBs, even though, the magnitude of the estimated concentrations tends to be slightly higher than the measured values, particularly for the 1993/1995 dataset (Fig. 3). At the same time, several data points representing high concentrations are not captured by the model. The deterministic mean and the expected value computed from Eq. (1) are similar. From Fig. 6, it is clear that the deterministic results and the uncertainty analysis results have very little bias since the intercept of their regression lines with the observed concentrations is very close to zero. However, from the same regression lines, it is also evident that predicted mean and predicted expected value do not vary consistently with the observed data. In fact, both slopes of the regression lines differ from one. It is likely that the occasional high concentrations recorded at the end of the Niagara River are the result of poorly characterized loadings or scouring events and/or errors in reporting the data. This source of uncertainty is accounted for in the model, including both the parameter uncertainties and variability of the stochastic components. In fact, when the measured concentrations are compared with the results of the uncertainty analysis (Figs. 4 and 5), the major part of the points is included within an upper and a lower bound, expressed as a

one-standard deviation interval about the expected value. In any given month in the monitoring period, there are between three to five recorded observations. A detailed review of the observed concentrations versus the predicted one-standard deviation interval about the expected value indicates that for any given month, the percentage of observations falling within the interval varies. The percentage of included observations ranges between 25% and 80% with a few months in which none of the observations are included in the interval and one occasion in which all observations are included. These results suggest that the estimated concentrations of Total PCBs are not normally distributed. Overall, however, accounting for the combination of all the possible input scenarios captures both the average trend and the variability of the measured concentrations in this modeling effort. From the results of the uncertainty analysis performed based on the selected stochastic components and available datasets for the monitoring period from May 1995 to March 1997, it is important to note the following: 1. The overall expected value of the Total PCB concentrations computed using the uncertainty analysis depends on the number of stochastic components used in the uncertainty analysis (Fig. 7). The trend of the mean values remains fairly constant for any number of stochastic components used. However, the addition of one extra stochastic component requires the recalculation of the expected value which could be greater or smaller than the previously computed expected value, depending on the degree of skewness of the added stochastic component. 2. The largest contribution to the overall variability of the results is given by the inflow concentrations from Lake Erie (Figs. 8–11 and Table 4). Adding one stochastic component to the model modifies the computed standard deviation by a certain percentage of the total standard deviation calculated using all the considered stochastic components. The portion of overall onestandard deviation contributed by each stochastic component is summarized in Table 4 and was computed as the difference between the widths of the one-standard deviation interval of two successive runs. It is worth pointing out that the order with which the stochastic components were added in each single run was arbitrary. However, changing the order of addition of the stochastic components did not change the order of influence of each component, nor the width of the overall one-standard deviation interval. It did, however, change the portion of the one-standard deviation interval contributed by some of the stochastic components. In the case of inflow concentration, for example, adding Cin as a stochastic variable at the end would indicate a slightly smaller contribution to the overall one-standard deviation interval than adding it at the beginning. This suggests that there exists some functional correlation among components, due to the model structure and dependence among model variables and parameters. This functional correlation might not have an effect on the overall uncertainty, however it remains an interesting point for future exploration. 3. The other stochastic components that bring significant contributions to the overall variability of the model results are indicated in Table 4. The variability of air and water temperature,

502

S. Franceschini, C.W. Tsai / Advances in Water Resources 33 (2010) 493–503

the biodegradation rate and the roughness coefficient do not bring any considerable contribution to the overall uncertainty of the model results. 4. Depending on the order in which stochastic components were added to the uncertainty analysis, the variability of certain components (most notably non-point source, point sources and air deposition) contributes negatively to the overall standard deviation. Although this outcome is counter intuitive, it can be explained as the result of the mathematical derivation of the Modified Rosenblueth method. In fact, when no correlation is used Eq. (1) can be simplified into Eq. (4)

EðY m Þ ¼ ½ð1  NÞY m þ

N X m ½piþ Y m iþ þ pi Y i ;

ð4Þ

i¼1

where the first term refers to the concentration obtained by using all mean values (deterministic mean) and the second term is the sum of the weighted averages of the concentrations Y :iþ: and Y :i: , N is the number of stochastic components considered. If there is no correlation among the terms, the addition of an extra stochastic component, as described in Table 3, implies the addition of a weighted average to the second term in Eq. (4) and the subtraction of a deterministic mean from the first term in Eq. (4). If the weighted average is less than the deterministic mean, then the recomputed moment is less than the value computed in the preceding run. Therefore, the difference between the standard deviations of two successive runs might be negative. 5. The derivation of Eq. (4) above explains also why occasionally the standard deviation computed for certain runs appears to be greater than the overall standard deviation (PS run Fig. 11). In uncertainty analysis, as opposed to the sensitivity analysis, the combined effect of all the stochastic components is evaluated at the same time. Thus the variability of certain components is compensated by other components.

7. Conclusions In this paper, we have coupled a hydrodynamic model and a water-quality model with a recently developed uncertainty analysis method to quantify the overall variability of the estimated Total PCBs concentrations in the Niagara River, based on the uncertainty of model parameters and the natural variability embedded in selected input variables. The results can be attributed in large part to the variability of the inflow concentration from Lake Erie into the Niagara River. Any attempt to limit the PBC loadings from the Niagara River to Lake Ontario should come through better management of the Total PCBs releases and concentrations in the Great Lakes basin upstream of the river, both in terms of dissolved concentrations and incoming suspended sediments. Nevertheless, other sources of variability should not be disregarded. While the influence of the variability of flow and atmospheric deposition cannot be controlled, their variability does not have for the most part a very significant effect on the outflow concentration. On the other hand the variability of the point and non-point source loadings to the Niagara River does have a significant effect. These two stochastic components should be better characterized. For both loadings, but especially for the non-point sources, there is either very limited or obsolete data. A renewed effort should be enacted to assess the effective runoff from IHWS considering that some sites have been remediated. In addition, a study of the presence along the Niagara River of possible new sources of contamination should be considered. This paper has contributed to modeling the fate and transport of toxic contaminants in the Niagara River. Primarily, it provides a framework to account for uncertainty of both input variables and

model parameters as a first step to optimize and improve the modeling effort. Furthermore, simulated concentrations that included uncertainty of model variables and parameters offer a better approach to compare model results with measured data in terms of both the mean trend and variability. Acknowledgements The authors would like to thank Donald J. Williams and the Aquatic Sciences and Reporting Ecosystem Health Division Environment Conservation Branch of the Ontario Region in Canada, the US Army Corps of Engineers, Buffalo District, NY, and the US Coast Guard, Station Niagara, NY for providing essential information for this study. The authors also appreciate the technical support received from Dr. Shawn Mattot on the OSTRICH algorithm. The work reported in this paper is financially supported by the National Science Foundation Graduate Fellowship; the Civil, Structural and Environmental Engineering Department; and the Environment and Society Institute at the State University of New York at Buffalo. This work is also partially supported by the National Science Foundation under grant contract number EAR0510830. References [1] Blanchard P, Audette CV, Hulting ML, Basu I, Brice KA, Chan CH, Dryfhout-Clark H, Froude F, Hites RA, Neilson M. Atmospheric deposition of toxic substances to the great lakes: IADN results through 200, Toronto, Ontario, CA; 2000. [2] Brown LC, Barnwell, TO Jr. The enhanced stream water quality models QUAL2E and QUAL2E-UNCAS: documentation and user manual., Rep. EPA/600/3-87/ 077, US environmental protection agency, Athens, GA; 1987. [3] Chadderton RA, Miller AC, McDonnell AJ. Uncertainty analysis of dissolved oxygen model. J Environ Eng – ASCE 1982;108(5):1003–13. [4] Chang WT. Analysis of contaminant loadings for the niagara river, MS Thesis, SUNY at Buffalo, Buffalo, NY; 1996. [5] Environmental Protection Agency (USEPA). Results of the nationwide urban runoff program. Water Planning Division, Washington, DC; 1983. [6] Franceschini S. Modeling toxic concentrations in the niagara river under uncertainty, MS Thesis, SUNY at Buffalo, Buffalo, NY; 2004. [7] Franceschini S, Tsai C. Incorporating reliability into the definition of the margin of safety in total maximum daily load calculations. J Water Resour Plan Manage – ASCE 2008;134(1):34–44. [8] Gates TK, Al-Zahrani MA. Spatiotemporal stochastic open-channel flow. 2. Simulation experiments. J Hydraul Eng – ASCE 1996;122(11):652–61. [9] Hamrick JM. A three-dimensional environmental fluid dynamics computer code: theoretical and computational aspects., S.R. 317, ed., The College of William and Mary, Virginia Institute of Marine Science; 1992. [10] Han KY, Kim SH, Bae DH. Stochastic water quality analysis using reliability method. J Am Water Resour Assoc 2001;37(3):695–708. [11] Harmel DR, Smith PK. Consideration of measurement uncertainty in the evaluation of goodness-of-fit in hydrologic and water quality modeling. J Hydrol 2007;337(3-4):326–36. [12] Harrison KW. Two-stage decision-making under uncertainty and stochasticity: Bayesian programming. Adv Water Resour 2007;30(3):641–64. [13] Jaffe PR, Ferrara RA. Modeling sediment and water column interactions for hydrophobic pollutants – parameter discrimination and model response to input uncertainty. Water Res 1984;18(9):1169–74. [14] Lian Y, Yen BC. Comparison of risk calculation methods for a culvert. J Hydraul Eng – ASCE 2003;129(2):140–52. [15] Loucks DP, Lynn WR. Probabilistic models for predicting stream quality. Water Resour Res 1966;2(3):593–605. [16] Mailhot A, Villeneuve J-P. Mean-value second-order uncertainty analysis method: application to water quality modelling. Adv Water Resour 2003;26(5):491–9. [17] Mattot LS. Ostrich: an optimization software tool, documentation and user’s guide, Version 1.9. Department of Civil, Structural, and Environmental Engineering, University at Buffalo, Buffalo, NY; 2006. Available from: http:// www.groundwater.buffalo.edu. [18] McLachlan M, Mackay D, Jones PH. A conceptual model of organic chemical volatilization at waterfalls. Environ Sci Technol 1990;24(2):252–7. [19] Melching CS, Bauwens W. Uncertainty in coupled non-point source and stream water-quality models. J Water Resour Plan Manage – ASCE 2001;127(6):403–13. [20] Novotny V. Integrating diffuse/nonpoint pollution control and water body restoration into watershed management. J Am Water Resour Assoc 1999;35(4):717–27. [21] Preston BL, Jones RN. Evaluating sources of uncertainty in Australian runoff projections. Adv Water Resour 2008;31(5):758–75.

S. Franceschini, C.W. Tsai / Advances in Water Resources 33 (2010) 493–503 [22] Rosenblueth E. Point estimates for probability moments. Proc Natl Acad Sci 1975:3812–4. [23] Tetra Tech. Development of the hydrodynamic and water quality models for the Savannah harbor expansion project. USACE – Savannah District Contract No.: DACA65-99-D-0065, Savannah, GA; 2006. [24] TetraTech. User’s Manual for Environmental Fluid Dynamics Code – Hydro Version. US Environmental Protection Agency, Region 4, Atlanta, GA; 2002. [25] Tsai CW, Franceschini S. Evaluation of probabilistic point estimate methods in uncertainty analysis for environmental engineering applications. J Environ Eng – ASCE 2005;131(3):387–95. [26] Tung Y-K, Yen BC. Hydrosystems engineering uncertainty analysis. New York: McGraw-Hill; 2005. [27] Tung YK, Hathhorn WE. Assessment of probability-distribution of dissolvedoxygen deficit. J Environ Eng – ASCE 1988;114(6):1421–35.

503

[28] Ünlü K. Assessing risk of groundwater pollution from land-disposed wastes. J Environ Eng – ASCE 1994;120(6):1578–97. [29] Vuksanovic V, De Smedt F, Van Meerbeeck S. Transport of polychlorinated biphenyls (PCB) in the Scheldt Estuary simulated with the water quality model WASP. J Hydrol 1996;174(1–2):1–18. [30] Wool T, Ambrose B, Martin J, Comer E. Water quality analysis simulation program (WASP) – V 6.0 User’s manual, US Environmental Protection Agency – Region 4, Atlanta, GA; 2002. [31] Yu PS, Yang TC, Chen SJ. Comparison of uncertainty analysis methods for a distributed rainfall-runoff model. J. Hydrol. 2001;244(1–2):43–59. [32] Yulianti JS, Lence BJ, Johnson GV, Takyi AK. Non-point source water quality management under input information uncertainty. J Environ Manage 1999;55(3):199–217.