A comparison of various uncertainty models: An example of subsurface contaminant transport

Available online at www.sciencedirect.com

Journal of Hydro-environment Research 6 (2012) 311e321 www.elsevier.com/locate/jher

Research paper

A comparison of various uncertainty models: An example of subsurface contaminant transport Saheb Mansour-Rezaei, Gholamreza Naser*, Rehan Sadiq Okanagan School of Engineering, The University of British Columbia, Kelowna, BC, Canada V1V 1V7 Received 9 April 2010; revised 21 October 2010; accepted 23 September 2011

Abstract Groundwater resources are under increasing threat of contamination and wasteful use in many parts of the world. Groundwater flow and integrated contaminant transport models are commonly used to predict the fate of contaminants in the subsurface environment. However, the lack of reliable data and complexity of the natural environmental systems, the predictions are subjected to large uncertainties. For reliable decision-making, these contaminant transport models are required to explicitly consider associated uncertainties in their parameters. This paper aims to compare the results of four common uncertainty models using an example of contaminant transport in groundwater. The research employed an advectionedispersion equation (ADE) to describe the transport of a contaminant in groundwater. For simplicity, two parameters e dispersion coefficient and velocity e were considered in the uncertainty analysis. Fuzzy set theory, one- and two-dimensional (1-D and 2-D) Monte Carlo simulations, and Probability Box (P-Box) methods were investigated. The cumulative distribution functions generated from these analyses were compared to evaluate the capabilities of these methods. The comparison showed that P-Box method provides a more comprehensive analysis with lesser assumptions as compared to other methods, and also found to be more pragmatic way to describe and propagate uncertainties in complex environmental systems. Furthermore, execution time required to perform uncertainty analysis using P-Box method is comparatively much less than 2-D Monte Carlo simulations. Crown Copyright
2011 Published by Elsevier B.V. on behalf of International Association for Hydro-environment Engineering and Research, Asia Pacific Division. All rights reserved. Keywords: Contaminant transport; Uncertainty analysis; Probability box; Monte Carlo simulations; Fuzzy set theory

1. Introduction The level of uncertainty associated with a system is proportional to its complexity, which arises as a result of vaguely known relationships among various entities, and randomness in the mechanisms governing the domain. The complex systems like environmental, socio-political, engineering, or economic systems, which involve human interventions with vast arrays of inputs and outputs cannot all possibly be captured analytically or controlled in a conventional sense. Water systems are extremely complex and dynamic in nature. To understand the impacts of the physico-

* Corresponding author. Tel.: þ1 250 300 7274; fax: þ1 250 807 9850. E-mail address: [email protected] (G. Naser).

chemical, biological and socio-economical changes in the surrounding environment, tremendous efforts have been devoted to enhance the body of knowledge about various processes occurring in hydrosystems. However, the lack of reliable data due to resource constraints and complexities inherent in the natural environmental systems limits human ability to make correct predictions. Therefore, proper treatment of uncertainties is required for modeling hydrosystems. 1.1. Taxonomy of uncertainties Uncertainties can arise from various sources in water resources engineering. These include data uncertainty, structural uncertainty (raised from the imperfect description of physical reality by a limited number of mathematical relations), and parameters uncertainty (Mannina and Viviani, 2010; Freni

1570-6443/$ - see front matter Crown Copyright
2011 Published by Elsevier B.V. on behalf of International Association for Hydro-environment Engineering and Research, Asia Pacific Division. All rights reserved. doi:10.1016/j.jher.2011.09.001

312

S. Mansour-Rezaei et al. / Journal of Hydro-environment Research 6 (2012) 311e321

et al., 2009; Willems, 2008). Data uncertainties may include (but are not limited to) measurement errors, inconsistency and non-homogeneity of data, data handling and transcription errors, and inadequate representation of data sample due to time and space limitations. Moreover, the uncertainty in the estimation of model parameters implicitly considers the error induced by incorrect model structure, data errors and climatic variation factors. Simonovic (1997) indicated randomness and lack of knowledge as the two major sources of uncertainty. As a result, uncertainty is either an objective fact of the phenomenon under consideration or a subjective impression of human perception (Zimmermann, 2001). Aleatory uncertainty (also known as stochastic or objective uncertainty) results from the fact that a system can behave in random ways. In general, the uncertainties due to inherent randomness of a physical process cannot be eliminated or reduced. An epistemic uncertainty (also known as subjective or ignorance) results from the lack of knowledge about a system. Due to stochastic nature of natural phenomena, an aleatory uncertainty is always associated with the natural systems. On the other hand, an epistemic uncertainty is reducible by data analysis, making additional monitoring, and deepening our understanding and knowledge of the phenomenon. The traditional approach to handle an aleatory uncertainty is probabilistic analysis based on historical data (a frequentist approach). An epistemic uncertainty, on the other hand, has traditionally been addressed through Bayesian approach even though the approach was limiting as it required priori assumptions (Sentz and Ferson, 2002). 1.2. Uncertainty analysis in water resources engineering In water resources engineering, the design quantities and system outputs depend on several system parameters, and not all of them can be quantified with absolute accuracy. Thus, several techniques with different levels of mathematical complexity and data requirements have been reported in the literature for conducting uncertainty analysis in different areas of water resources engineering. Selection of an appropriate technique to be used in a particular problem depends strongly on the nature of the problem, availability of information, resources constraint, model complexity, and type and the desired level of the accuracy and/or reliability of the results. Overall, the uncertainty methods are generally classified into the two groups of analytical and approximation techniques. Table 1 provides a summary of this classification along with several recent applications in water resources engineering. Analytical techniques: These techniques (e.g., first and second order reliability methods) are mathematically less demanding, and therefore can be implemented straightforwardly. Their applications, as indicated in Table 1, are fairly limited due to oversimplification and/or assumptions. Approximation techniques: These techniques are mathematically more complex and have traditionally been applied to a broad range of water resources problems. The most common examples of approximation techniques are fuzzy set theory, one- and two-dimensional (1-D and 2-D) Monte Carlo simulations, Bayesian and Probability Box (P-Box) methods.

The fuzzy set theory, 1-D and 2-D Monte Carlo simulations, and P-Box method have been frequently applied to conduct an uncertainty analysis in water resources area. Although Guyoanet et al. (1999) reported more realistic results for the fuzzy set theory (in human health and environmental risk analyses where response values at the tails are important), there is still no clear and direct study comparing the performance and robustness of these techniques. Thus, the primary objective of this study is to compare the four common approximation techniques (i.e., fuzzy set theory, 1-D and 2-D Monte Carlo simulations, and P-Box method) with regards to their conservativeness, execution time, ease of formulation, and complexity in uncertainty analysis in water resources. The study advances by applying the techniques to a case of contaminant transport in a groundwater flow. Both numerical and analytical solutions of the governing advectionedispersion equation (ADE) are used to compute temporal and spatial variations of a contaminant concentration. Two parameters, namely, velocity and contaminant dispersion coefficient, are identified to perform uncertainty analyses using above four methods. 2. Contaminant transport in groundwater Groundwater resources are vulnerable to a wide variety of hazards that could potentially limit their ability to perform satisfactorily. Groundwater pollution resulting from agriculture, industrial, and waste-disposal activities is a very serious problem that often requires extensive treatments. In case of groundwater pollution, remediation procedures are required to keep contaminant concentrations below threshold limits. Moreover, the remediation procedures are often cost-effective and require careful evaluations of the extent of the pollutions. On the other hand, such evaluations rely on the accuracy of the data for aquifer properties. Vague or imprecise information especially arises in the identification and determination of aquifer properties. Comprehensive collection of aquifer data is extremely expensive due to their large spatial and temporal variations of the effective parameters. The lack of field data and parameters variations may have impacts on the reliability of a model’s results. The diversity of uncertainty sources presents a great challenge to groundwater systems’ planning and management. Therefore, a comprehensive modeling approach is required to 1) consider parameters’ uncertainties, 2) propagate uncertainties throughout the system, and 3) help authorities to wisely make informed decisions at planning, design, and management levels. This is particularly important because conventional deterministic techniques in practice are unable to account for possible variations of system responses. 2.1. Governing equation Fig. 1 shows a typical one-dimensional contaminant transport in a uniform groundwater flow field. The aquifer is assumed to be homogeneous and isotropic. For simplicity, the contaminant is considered conservative; therefore no decay or growth process is considered. Assuming a steady uniform flow

S. Mansour-Rezaei et al. / Journal of Hydro-environment Research 6 (2012) 311e321

313

Table 1 Some recent examples of uncertainty analyses in water resources engineering. Group

Techniques

Application

Target uncertain parameters

Analytical techniques (Tung, 2004)

 Derived distribution technique  Fourier transform technique  Laplace and exponential transform techniques  Mellin transform technique

 Water resources  Climate change

 Hydrologic parameters  Hydraulic parameters

Approximation techniques

 First-order variance estimation method (Tung, 2004)  Probabilistic point estimation methods (Tung, 2004) Stochastic response surface method (SRSM)

 Water resources  Climate change

 Hydrologic parameters  Hydraulic parameters

 Groundwater flow (Balakrishnan et al., 2005)  Uranium decay chain transport in the groundwater (Nair et al., 2009)  Groundwater flow (Balakrishnan et al., 2005)  Dynamics of trace metals, metalloids, radio-nuclides in saturated porous media (Wang et al., 2003)  Distributed watershed models (Yang et al., 2008)  Urban drainage system management (Freni et al., 2009)  Stochastic groundwater flow (Hassan et al., 2008)  Distributed watershed models (Yang et al., 2008)  River basin water quality (van Griensven and Meixner, 2006)  Distributed watershed models (Yang et al., 2008).  To assess the global freshwater availability at a sub-basin level (Schuol et al., 2008)  Distributed watershed models (Yang et al., 2008)  Groundwater flow and mass transport predictions (Fu and Jaime Gomez-Hernandez, 2009)  Groundwater flow and transport models (Elfeki, 2006)  Distributed watershed models (Yang et al., 2008)  Subsurface flow and particles travel time (Yuan and Druzdzel, 2006)  Contaminated groundwater system (Dou et al., 1997)  Random and nonrandom uncertainties in contaminated groundwater system (Zhang et al., 2009)  Contaminant leakage into the soil (Guyoanet et al., 1999)  Microbial-mediated contaminant in aquifer (Mohamed et al., 2006)

 Hydraulic conductivity and recharge rate field  Velocity, dispersion coefficient, porosity, etc.

High dimensional model representation (HDMR)

Generalized likelihood uncertainty estimation (GLUE)

Parameter solution (ParaSol)

Sequential uncertainty fitting algorithm

Bayesian inference based on Markov chain Monte Carlo (MCMC)

Bayesian inference based on importance sampling (IS)

Fuzzy set theory

Monte Carlo simulation

DempstereShafer structure (P-Box)

 Risk assessment of contaminated land

 Hydraulic conductivity  Recharge rate  Biotic and abiotic reaction parameters

 Outlet discharge parameters  Quantity and quality modules  Hydraulic conductivity and recharge

 Outlet discharge parameters  Outlet discharge parameters

 Input parameters of a hydrological model  Outlet discharge

 Outlet discharge parameters  Hydraulic conductivity, piezometric head, dispersion  Geological structure

 Outlet discharge parameters  Hydraulic conductivity and porosity

 Velocity and dispersion coefficient  Hydraulic conductivity

 Permeability, porosity, hydraulic conductivity, dispersion coefficient of aquifer  Microbial biomass concentration and hydraulic conductivity  Concentration in the soil, characteristics of the soil, and variables affecting the long term exposure to cadmium (Sander et al., 2006)

314

S. Mansour-Rezaei et al. / Journal of Hydro-environment Research 6 (2012) 311e321

differencing scheme of an implicit finite-difference method is used in this study to numerically integrate Equation (1). The method is approximated as:
tþ1 
tþ1  t t t Ci  Cti Ciþ1  2Ctþ1 þ Ctþ1 i i1 þ Ciþ1  2Ci þ Ci1 ¼D Dt 2Dx2
tþ1  tþ1 Ciþ1  Ci1 þ Ctiþ1  Cti1 V  4Dx ð2Þ

Fig. 1. Continuous contaminant leaking into shallow aquifer. Observation location at x ¼ 1220 m at two different times (400 days and 800 days).

velocity, the spatial and temporal variations of the contaminant concentration at any specified location and time can be mathematically modeled by only an unsteady 1-D ADE expressed by the following parabolic partial differential equation: vC v2 C vC ¼D 2 V vt vx vx

ð1Þ

Here, C is the contaminant concentration (ML3); V is the flow velocity (LT1) in x-direction; and D is the dispersion coefficient (L2T1). 2.2. Boundary conditions Regardless of a solution approach, Equation (1) must be solved with respect to the boundary conditions. Similar to Dou et al. (1997), this study examines two types of conditions at inlet and outlet boundaries of the flow field:  Inlet Boundary Condition: This includes a stepwise continuous release of a contaminant with constant concentration Co at x ¼ 0.  Outlet Boundary Condition: This is considered as a noflux boundary (i.e.vC=vx ¼ 0) at x ¼ lx in which lx is some characteristic length-scale representing the size of the contaminant plume.

2.3. Solution approaches Solving differential equations such as (1) is a fundamental problem in science and engineering. Occasionally, one can find closed-form analytical solutions for very simple cases. However, numerical approximations are perhaps the only solution methods to more mathematically complex systems. 2.3.1. Numerical solution Several approaches have been traditionally used to numerically solve Equation (1). These include (but are not limited to) finite element method, finite-difference method, finite volume method, and boundary volume method. A central

where Dt and Dx are time and space increments (T and L). The index “i” refer to location in x-direction, while the super index “t” represents time. It should be noted that Equation (2) will eventually result in a system of algebraic three-diagonal linear equations that must be solved for concentration at each node in space and time with regards to the aforementioned boundary conditions. 2.3.2. Analytical solution Given the conditions at inlet/outlet boundaries introduced in Section 2.2, Equation (1) can also be solved analytically providing the following general solution (Dou et al., 1997):       C 1 x  Vt 1 Vx x þ Vt ð3Þ erfc pffiffiffiffiffi ¼ erfc pffiffiffiffiffi þ exp C0 2 2 D 2 Dt 2 Dt Here, “exp” and “erfc” are the exponential and complementary cumulative error functions, respectively. Note that this study applies the analytical approach (i.e., Equation (3)) to verify the validity of the proposed numerical model by a direct comparison of the results. 2.4. Case study This research applies the 1-D contaminant transport equation in a uniform and hydraulically steady groundwater flow (Fig. 1) initially introduced by Dou et al. (1997). A conservative contaminant with concentration of Co ¼ 100 mg/L is constantly released at the inlet section. The characteristic length-scale (i.e., the size of the contaminant plume) is lx ¼ 1525 m. The Cartesian coordinate axes were set up at the release point. Similar to Dou et al. (1997), the 1-D flow domain in this study is discretized into a numerical mesh consisting of 101 equally spaced nodes with spatial increment (Dx) of 15.25 m. A time step (Dt) of one day is also applied to numerically integrate Equation (2). Additional information about the case study can be found in Dou et al. (1997). It should also be emphasized that this study reports the analytical and numerical results for the contaminant concentration at a point located at x ¼ 1220 m (i.e., node number 81) after 400 and 800 days of contaminant release. 3. Uncertainty models Several techniques have been reported in the literature to determine the magnitude of the dispersion coefficient including laboratory and natural gradient field tests (Dou et al.,

S. Mansour-Rezaei et al. / Journal of Hydro-environment Research 6 (2012) 311e321

1997). At field scale, however, it is very difficult (if not impossible) to precisely determine dispersion coefficient because of its scale dependence. On the other hand, applying ideal assumptions and simplifications, groundwater velocity is a highly uncertain parameter as well given that it is often derived from imprecise aquifer properties such as hydraulic conductivity. It should be emphasized that this research assumes velocity and dispersion coefficient are mutually independent parameters, although empirical evidence has suggested that dependency relations exist between velocity and dispersion coefficient (Klotz, 1980). 3.1. Fuzzy set theory Fuzzy sets, introduced by Zadeh (1965), have been widely used in all engineering fields, decision making and control processes (Dou et al., 1997). A fuzzy set is defined as a class of objects with a continuum of grades of membership with values from zero to one (Li, 2003) as follows: AðxÞ ¼ fðx; mA ðxÞÞ; x˛X; mA ðxÞ˛½0; 1g

ð4Þ

where X ¼ {x} is a universal set of elements, A(x) is a fuzzy set of X, and mA(x) is the membership grade of x in A. The parameter mA(x) is a real number in the interval [0, 1] that zero represents non-membership and one represents full membership. A fuzzy set A is normal if and only if at least one x exists such that mA(x) ¼ 1. The fuzzy set A is convex if for every real number a, b, and c, mA ðcÞ  min½mA ðaÞ; mA ðbÞ given that a < b < c. This means that the membership function may increase or decrease. The a-level cut, described as Aa ¼ fxjmA ðxÞ  ag, consists of all components of X whose membership grade is greater than or equal to a. Fig. 2 shows the concepts of the a-level cut and support in a fuzzy set. The support of the fuzzy set A is the set Supp(A) ¼ {xjmA(x) > 0}. The convexity condition ensures that the support is an interval. In fuzzy set theory, membership function plays a vital role. Two commonly used membership functions for characterizing fuzzy numbers are a bounded bell-shape and triangular functions. Due to its simplicity (Li, 2003), this research focuses on triangular membership functions to define the fuzzy input parameters. An example of a triangular fuzzy number (TFN) is

shown in Fig. 2. A TFN can be defined by specifying three values on universe of discourse: the most credible value, the lowest possible value, and the highest possible value. The most credible value is assigned a membership value of one, and number that falls short of the lowest possible value or exceeds the highest possible value will get a membership grade of zero (Dou et al., 1997; Li, 2003). The low (L), most likely (M), and high (H) values for velocity and dispersion coefficient of the case considered in this research are provided in Table 2. To compute the uncertainties of contaminant concentrations in groundwater, Li (2003) proposed a three-step fuzzy methodology that is adopted in this study. The steps are: 1) define input parameters which are known imprecisely, 2) definition of fuzzy parameters and selection of fuzzy operators, and 3) analysis of fuzzy modeling outputs. In this work, a numerical approach, called vertex method is applied for the fuzzy analyses. The vertex method (Dong and Shah, 1987) is based on the a-level cut concept and interval analysis technique. This method is easy to implement and has been successfully applied to solve contaminant transport problems (Dong and Shah, 1987; Dou et al., 1997). The algorithm of this approach is described as following: 1. Consider fuzzy set A1,A2,.,An defined on the real line R, and denoted an element of Ai by xi (where i ¼ 1,2,.,n and n is the total number of uncertain input parameters). Output y can be related to inputs x1,x2,.,xn through y ¼ f(x1,x2,.,xn). 2. Discretize the range of membership grade [0, 1] into a finite number of values a1,a2,.,aM. The number of discretization depends on the degree of accuracy in approximation. 3. For each membership value aj, find the corresponding interval for Ai in xi for i ¼ 1,2,.,n. These are the supports

Table 2 Required data in different uncertainty models. Method

Low (L) ¼ 3.45 Most likely (M) ¼ 36.72 High (H) ¼ 61.38

1-D Monte Carlo

Low (L) ¼ 1.8 Most likely (M) ¼ 2.4 High (H) ¼ 3.1

Low (L) ¼ 3.45 Most likely (M) ¼ 36.72 High (H) ¼ 61.38

2-D Monte Carlo

Low (L) (m,s) ¼ (1.8, 0.1) Most likely (M) (m,s) ¼ (2.4, 0.1) High (H) (m,s) ¼ (3.1, 0.1)

Low (L) (m,s) ¼ (3.45, 1) Most likely (M) (m,s) ¼ (36.72, 7) High (H) (m,s) ¼ (61.38, 7)

P-Box

Low (L) ¼ 1.8 Most likely (M) ¼ 2.4 High (H) ¼ 3.1

Low (L) ¼ 3.45 Most likely (M) ¼ 36.72 High (H) ¼ 61.38

Membership Grade

0.75 α-Level cut

0.25 0

Highe st possible value Support Parameter Value

Fig. 2. Example of a triangular fuzzy set.

Dispersion coefficient (m2 day1)

Low (L) ¼ 1.8 Most likely (M) ¼ 2.4 High (H) ¼ 3.1

1

Lowe st possible value

Velocity (m day1)

Fuzzy set theory

Most Credible Value

0.5

315

S. Mansour-Rezaei et al. / Journal of Hydro-environment Research 6 (2012) 311e321

of the aj-cuts of A1,A2,.,An. Denote the end points of these intervals by ½a1 ; b1 ; ½a2 ; b2 ; .; ½an ; bn . In case a1 equals b1, the interval reduces to a point estimate. 4. Take one end point from each of the intervals. The obtained end points can be combined into an n-array. There are 2n distinct permutations, which gives 2n combinations for the vector ðx1 ; x2 ; .; xn Þ. 5. Evaluate the function f ðx1 ; x2 ; .; xn Þ for each of 2n combinations and obtain 2n values for y. Denote these by y1,y2,.,yn. The desired interval for y is given by ½^ yk ; n yk . These points define the support of aj-cut. k

k

6. Repeat the process for other a-level combinations to obtain additional ½^ yk ; n yk . k

k

After substituting fuzzy velocity and dispersion coefficient parameters in Equations (2) and (3), the vertex method is used to compute contaminant concentration at desired distance and times. Similar to Dou et al. (1997) and Li (2003), we have carried out the analyses at five different a-levels: 0, 0.25, 0.5, 0.75 and 1. By applying the same fuzzy inputs, the fuzzy finite-difference simulation results compared with the corresponding analytical results to evaluate the validity of the numerical model. Fig. 3 illustrates the membership functions of concentration for both analytical and numerical results at a point located at x ¼ 1220 m at two different times (i.e., 400 days and 800 days after release). The illustration clearly indicates that the numerical results are in excellent agreement with the analytical results. Fig. 4 highlights the influence of flow hydraulics and dispersion on contaminant transport. The illustration depicts the front evolution for the concentration through the groundwater and clearly shows a symmetric dispersion of the contaminant around the mean front. Fig. 4 also shows the lower, upper, and most likely bounds of concentration (at the zero-level-cut) as a function of time at a point located at x ¼ 1220 m. In practice, this graphical representation can be used to determine the maximum and minimum amount of contaminant concentration at any given point and time, and help estimating dynamic risk associated with the contaminant. Such information not only will help the engineers to identify the critical points of the system but also 1.2

Membership Function

FDM, T=400 day FDM, T=400 day

1.0

Analytical, T=400 day Analytical, T=400 day

0.8

FDM, T=800 day FDM, T=800 day

0.6

Analytical , T=800 day Analytical , T=800 day

0.4 0.2 0.0 0.0

0.2

0.4

0.6

0.8

1.0

1.2

Contaminat Concentration (C/C0)

Fig. 3. Analytical and numerical results at location of x ¼ 1220 m after 400 and 800 days.

1.2 1 0.8

C/C0

316

Upper Bound Lower Bound Most likley

0.6 0.4 0.2 0 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 850 900

Time ( Day)

Fig. 4. Upper, lower and most likely bounds of concentration at x ¼ 1220 m (computed by finite-difference method).

will help the engineers or other authorities to make more informed decisions when required. It should be emphasized that in this paper the output of alevel sets in the fuzzy set theory by searching for the upper and lower bounds of the corresponding output range. Sufficient accuracy can only be obtained by a larger number of calls of the output function evaluated based on the input data and that increases considerably with the number of input variables. Therefore, the fuzzy set theory may be feasible only in the case of medium size problems with a few input variables. On the other hand, monotonicity of the output function may increase accuracy and therefore help reducing the number of computations required. These need to be investigated in future. 3.2. Monte Carlo simulations Monte Carlo simulation is a classical method that uses a class of computational algorithms that rely on repeated random sampling of input parameters (i.e., velocity and dispersion coefficient) to compute response of the groundwater system (concentration of contaminant). 3.2.1. 1-D Monte Carlo Monte Carlo simulation is conceptually very simple and called a parametric method. It assumes that the model parameters (i.e., velocity and dispersion coefficient) are random variables, which follow a specific probability density function (PDF) (Guyoanet et al., 1999). First step in Monte Carlo simulation is to define PDFs for velocity and dispersion coefficient. To compare the results of Monte Carlo simulation with those of the fuzzy set theory (Guyoanet et al., 1999), we assume triangular PDFs similar to shape of TFNs for both parameters (velocity and dispersion coefficient). The main difference between PDFs and fuzzy membership function is that the area under the PDF is unity. Required data for 1-D Monte Carlo simulation is shown in Table 2. The cumulative distribution functions (CDFs) are first obtained by integrating PDFs. The CDFs are used for random sampling of velocity and dispersion coefficient and are implemented in Equation

S. Mansour-Rezaei et al. / Journal of Hydro-environment Research 6 (2012) 311e321

(2) to compute concentration. The procedure is repeated for a large number of times. Fig. 5 shows the results of the 1-D Monte Carlo simulation as a relative frequency diagram. To obtain a consistent and smooth relative frequency diagram, 1000 iterations were performed. The frequency diagram was obtained by dividing the range of contaminant concentration into 30 classes of equal width and calculating the frequency of each class. As the values of input parameters are chosen randomly during Monte Carlo simulations, the values with low probability have less chance of being randomly selected, especially when the iterations are limited. However, in fuzzy set theory, all possible combinations of parameters values are considered and the maximum and minimum model outputs are directly reflected in the output uncertainty (Masky, 2004). Fig. 6 compares fuzzy set theory results with those of 1-D Monte Carlo simulation at the location of x ¼ 1220 m and at time t ¼ 400 days. Please note that to make such a comparison possible, the relative frequency given in Fig. 5 was multiplied by a factor of “5.2” such that the maximum relative frequency equals to unity (Guyoanet et al., 1999). Fig. 6 clearly indicates the greater spread in case of fuzzy set theory as compared to 1-D Monte Carlo simulation results. This highlights that fuzzy set theory provides more conservative results than the 1-D Monte Carlo simulations. The greater spread can be attributed to the fact that all possible combinations of velocity and dispersion coefficient values are considered in a fuzzy analysis. 3.2.2. 2-D Monte Carlo As mentioned before, the aleatory and epistemic uncertainties are the two major types of uncertainties. The obvious advantage of 2-D over 1-D Monte Carlo simulations is that the former can distinguish between aleatory and epistemic uncertainties. A 2-D Monte Carlo simulation is carried out by sampling epistemic variables in the outer loop and aleatory variables in the inner loop (Karanki et al., 2009). Fig. 7 provides the procedure to implement 2-D Monte Carlo simulation and explained below in stepwise manner (Durga Rao et al., 2007; Karanki et al., 2009): 1. Information regarding the PDFs of the model parameters (i.e., velocity and dispersion coefficient) and the

Fig. 5. Histogram obtained from Monte Carlo simulation (1000 iterations).

317

Fig. 6. A comparison between Monte Carlo simulations and fuzzy set results. (To make a fair comparison, the relative frequencies were multiplied by a factor such that the maximum relative frequency equals 1.)

2.

3.

4.

5.

uncertainty in the parameters of PDFs (epistemic uncertainty in the parameters) were obtained. The current study assumed triangular PDFs for both parameters of the model. Also, normal distributions (N w m, s) were used to describe epistemic uncertainty for low (L), most likely (M), and high (H) values of triangular distribution of model parameters. Table 2 provides the mean and variance of L, M, and H values for velocity and dispersion coefficient. Distributions for PDF parameters of components were first sampled using MATLAB. This action took place in the first loop of the two-loop sampling. The outer loop or first loop focused on epistemic uncertainty, and the second loop or inner loop focused on aleatory uncertainty. Epistemic variables were treated as constants inside the second loop. The sampled low, most likely and high values of velocity and dispersion coefficient were passed on to the second loop to generate triangular PDFs. In the second loop, simulation was carried out by sampling aleatory variables. Step 3 was repeated for sufficient number of times (say, 1000 iterations). Similar to 1-D Monte Carlo, 1000 iterations were found to be sufficient to obtain a smooth and consistent CDF curve. The required measures of aleatory uncertainty were obtained in this step. Sampling was repeated again for epistemic uncertainty and subsequently entered the second loop. The first loop had to be executed for 100 numbers of iterations. This was mainly because large number of iterations for the outer loop was very time-consuming and computationally demanding.

The results for the contaminant concentration are provided through a family of curves in Fig. 8, which shows the probability bounds for contaminant concentration at a point located at x ¼ 1220 m and at time t ¼ 400 days. Each CDF curve indicates the aleatory uncertainty, whereas the spread describes epistemic uncertainty. Fig. 8 also compares the 2-D Monte Carlo simulation results with those of the fuzzy set theory. The comparison shows that fuzzy set results enclose the 2-D Monte Carlo simulation results. As discussed before,

S. Mansour-Rezaei et al. / Journal of Hydro-environment Research 6 (2012) 311e321

V and D are implemented in Model

Check for no.of simulation (if n
No

yes

Sample Aleatory variables

318

Fig. 8. A comparison between 2-D Monte Carlo simulation and fuzzy set results.

CDF

Monte Carlo simulation applies a class of computational algorithms that rely on repeated random sampling of the uncertain parameters. Please also note that the scenarios in a Monte Carlo simulation that can combine low probability parameters values have a very less chance of being randomly selected. Hence, if probabilities are arbitrarily assigned to velocity and dispersion coefficient without being confirmed by a reliable set of field data, possible combinations of the parameter values will be eliminated from the analysis as a result. This could be significant, especially from decisionmaking point of view as it may have substantial environmental impacts.

No

No

if M>L

if H>M

Sample Epistemic variables

yes

yes

Triangular Shape PDF

Random (L)

Random (H)

Normal Dist.

Normal Dist.

Normal Dist.

Mean M SD

Mean H SD

Mean L SD

No

Check for no.of simulation (if n
Fig. 7. Procedural flowchart for 2-D Monte Carlo simulation.

out Put

yes

Random (M)

3.3. Probability box Using the DempstereShafer structure (Yager, 1986), a probability box (P-Box) provides a comprehensive way to describe and propagate uncertainties in a complex modeling environment. A P-Box is a class of distribution functions enclosed by an upper and a lower bound which represent the epistemic uncertainty in a distribution function of a random variable (Ferson et al., 2004). P-Box analysis can comprehensively account for possible deviations arising from uncertainty in distribution parameters, distribution shape or family, inter-variable dependence, and even model structure (Tucker and Ferson, 2003). P-Box is composed of a collection of intervals called focal elements, which paired with a probability mass that depends on the number of discretization. In this research, we assumed that the velocity and dispersion coefficient have triangular distributions with low (L), most likely (M), and high (H) as indicated in Table 2. The distribution functions for velocity and dispersion coefficient were, then, discretized at five levels resulting in a probability mass of 0.2 for each parameter. Table 3 shows how P-Box method can be used to determine uncertainties in contaminant transport modeling (Tucker and Ferson, 2003; Karanki et al., 2009). The matrix cells include an interval and a number representing the focal elements and associated probability mass, respectively. Note that the intervals and associated probability mass of velocity and dispersion coefficient were arrayed along the first raw, and the first column

S. Mansour-Rezaei et al. / Journal of Hydro-environment Research 6 (2012) 311e321

319

Table 3 Matrix of interval focal elements. Dispersion coefficient

Velocity [1.80e3.10] 0.20

[3.45e61.38] 0.2a [11.77e55.22] 0.20 [20.09e49.05] 0.20 [28.40e42.89] 0.20 [36.72e36.72] 0.2 a b

[0.000e0.570] [0.000e0.576] [0.000e0.574] [0.001e0.571] [0.003e0.570]

0.04b 0.04 0.04 0.04 0.04

[1.95e2.92] 0.20

[2.10e2.75] 0.20

[2.25e2.57] 0.20

[2.40e2.40] 0.20

[0.000e0.445] [0.000e0.438] [0.001e0.430] [0.003e0.421] [0.007e0.409]

[0.000e0.326] [0.000e0.313] [0.002e0.299] [0.008e0.283] [0.017e0.264]

[0.000e0.222] [0.002e0.207] [0.009e0.191] [0.022e0.172] [0.038e0.151]

[0.000e0.141] [0.007e0.126] [0.026e0.111] [0.051e0.094] [0.076e0.076]

0.04 0.04 0.04 0.04 0.04

0.04 0.04 0.04 0.04 0.04

0.04 0.04 0.04 0.04 0.04

0.04 0.04 0.04 0.04 0.04

Focal elements and probability mass. The maximum and minimum values of contaminant concentration and probability mass.

of the matrix in Table 3, respectively. Suppose U(i, j ) is an arbitrary cell inside the matrix on Table 3 at row i and column j (except the cells at the first raw and the first column of the matrix). Two ends of each interval in U(i, j ) show the maximum and minimum values of contaminant concentration computed by the corresponding dispersion coefficient and velocity interval at raw i and column j, respectively. To compute each interval in U(i, j ), four permutations based on maximum and minimum values of velocity and dispersion coefficient should be implemented in Equation (2). Because each probability mass of velocity and dispersion coefficient are 0.2, the probability masses associated with the intervals on the U(i, j ) are all 0.04 (¼0.2  0.2). After computing all U(i, j ), the CDF of relative contaminant concentration at specific location and time can be obtained. It should be noted that two ends of intervals of U(i, i) at the diagonal of described matrix (Table 3) are corresponding with lower and upper values of contaminant concentration computed by fuzzy set theory at five different a-levels: 0, 0.25, 0.5, 0.75 and 1. Fig. 9 reveals uncertainty of contaminant transport computed by P-Box method at a point x ¼ 1220 m and at time t ¼ 400 days after the release. Also, in Fig. 9 the P-Box results are compared with those of fuzzy set theory. Although the comparison indicates a fairly good agreement between the two methods, but it also shows that the P-Box results are bounding the fuzzy set results. It should be emphasized that to increase in the numbers of uncertain parameter of the model lead to increase of differences between fuzzy set theory and P-Box results. It should be noted that using P-Boxes can become very imprecise if the distance between the lower and upper distributions is large.

concentration, and the 95th percentile exceeds all but 5 percent of the values. Table 5 ranks and compares the four uncertainty models discussed in this research with respect to various qualitative attributes. The time of simulation for fuzzy set theory, 1-D Monte Carlo, 2-D Monte Carlo, and P-Box is 53.67, 708.34, 1.079E5, and 182.59 s, respectively. While the 2-D Monte Carlo simulation is the most time-consuming method, the PBox approach as well as fuzzy set is more computationally efficient. This is particularly important since computational costs can be a major challenge in engineering applications. Methods have to be developed in a way that they admit assertions about the sensitivity of the output with as few computations as possible. Generally, computational efficiency of the Monte Carlo simulation is amenable to the number of iterations (or sample size). Although a sample size of 1000 appeared sufficient in this research, it can be chosen independently of the number of input variables. It should also bring into attention that a Monte Carlo simulation, in general, enjoys re-sampling, which introduces additional computational efforts. Therefore, it seems feasible to apply a 2-D Monte Carlo simulation only if the problem is of medium size with a small number of random variables. Use of fuzzy set theory and P-Box methods to model uncertainty related to

4. Comparison of uncertainty models To have more illustrative comparison between the four different uncertainty models, Table 4 presents the 5th, 50th, and 95th percentile of 1- and 2-D Monte Carlo simulations [low (L) and high value (H)], fuzzy set theory [low (L), most likely (M), and high value (H)], and P-Box method [low(L) and high value (H)]. It should be noted that there are not any values for 5th percentile of 1- and 2-D Monte Carlo simulations. The 5th percentile value shows the lowest 5 percent of the relative concentration from the rest at specific location and time, the 50th percentile is the same as the median of relative

Fig. 9. A comparison of P-Box and fuzzy set results at location of x ¼ 1220 m and at time t ¼ 400 days.

320

S. Mansour-Rezaei et al. / Journal of Hydro-environment Research 6 (2012) 311e321

Table 4 5th, 50th and 95th percentiles obtained in different uncertainty models. Method

5th (%)

50th (%)

95th (%)

1-D Monte Carlo 2-D Monte Carlo (L, H) Fuzzy set theory (L, M, H) P-Box (L, H)

e e (0.000, 0.046,0.092) (0.000, 0.094)

0.065 (0.014, 0.127) (0.002, 0.151,0.299) (0.003, 0.306)

0.354 (0.165, 0.546) (0.062, 0.298,0.535) (0.076, 0.576)

L: low, M: most likely, H: high.

Table 5 Comparison of uncertainty models based on various qualitative attributes. Qualitative attributes

Fuzzy 1-D Monte 2-D Monte Probability set theory Carlo Carlo box

Conservativeness Computing time Ease of formulation Interpretation complexity Use in literature

2 1 1 1 2

4 2 3 4 1

3 3 4 3 3

1 1 2 2 4

parameters of hydraulic model is easier than Monte Carlo simulations, which need to select random parameters. 1-D Monte Carlo simulations frequently use by researchers to model uncertainty, but interpretation complexity of Monte Carlo simulations is much more than P-Box method and fuzzy set theory, which provide upper and lower bands of probability. 5. Conclusion Water resource systems in general include different types of interconnected components serving broad geographic regions. Therefore, they are at risk of failure due to natural hazards or anthropogenic causes, whether they are unintentional (operational errors and mistakes) or deliberate (terrorist act). This study compared four most common uncertainty models for contaminant transport in groundwater. To compute the concentration of contaminant in an aquifer at any specific place and time, an ADE was employed. Velocity and dispersion coefficient were assumed to be two independent uncertain variables. The finite-difference method was applied for numerical simulation of contaminant transport. Numerical results were compared with analytical solution, which proved the accuracy of numerical model. First, fuzzy set theory was used to propagate uncertainty associated with the velocity and dispersion coefficient. The results of uncertainty assessment, illustrated as concentration membership function, were in good agreement with the analytical results. Second, Monte Carlo simulations were used to model uncertainty. The comparison of results of 1-D Monte Carlo simulation and those of the fuzzy set theory showed that the latter was more conservative. Third, the application of a 2D Monte Carlo simulation revealed that such studies were often computationally inefficient to conduct in practice as it requires a huge number of calculations that considerably increases the simulation time. Finally, P-Box method was used

to compute bounds on response (contaminant concentration) cumulative distribution. The study indicated a good agreement between P-Box and fuzzy set results. In risk assessment models, sometimes there may be a dependency between input parameters. However, due to the lack of information, and sometimes due to modeling difficulties, it is often assumed that parameters are independent. On the other hand, the possible dependency between variables can have profound impacts on risk analysis and that can lead to over or under estimation (Bonissone et al., 2004). This research assumed that the parameters of model are mutually independent. For future work, the authors recommend to implement parametric Copula (Ferson et al., 2004) or T-norm operators to express the effect of dependency on risk assessment of contaminant transport. Acknowledgments The financial support provided by University of British Columbia is greatly acknowledged. References Balakrishnan, S., Roy, A., Ierapetritou, M.G., Flach, G.P., Georgopoulos, P.G., 2005. A comparative assessment of efficient uncertainty analysis techniques for environmental fate and transport models: application to the FACT model. Journal of Hydrology 307 (1e4), 204. Bonissone, P., Goebel, K., Yan, W., 2004. Classifier fusion using triangular norms. Lecture Notes in Computer Science 3077, 154e163. Dong, W., Shah, H.C., 1987. Vertex method for computing functions of fuzzy variables. Fuzzy Sets and Systems 24 (1), 65. Dou, C., Woldt, W., Bogardi, I., Dahab, M., 1997. Numerical solute transport simulation using fuzzy sets approach. Journal of Contaminant Hydrology 27 (1e2), 107. Durga Rao, K., Kushwaha, H.S., Verma, A.K., Srividya, A., 2007. Quantification of epistemic and aleatory uncertainties in level-1 probabilistic safety assessment studies. Reliability Engineering and System Safety 92 (7), 947. Elfeki, A.M.M., 2006. Reducing concentration uncertainty using the coupled Markov chain approach. Journal of Hydrology 317 (1e2), 1. Ferson, S., Nelsen, R.B., Hajagos, J., Berleant, D.J., Zhang, J., Tucker, W.T., et al., 2004. Dependence in Probabilistic Modeling, DempstereShafer Theory, and Probability Bounds Analysis. Sandia National Laboratories, Albuquerque. Freni, G., Mannina, G., Viviani, G., 2009. Uncertainty assessment of an integrated urban drainage model. Journal of Hydrology 373 (3e4), 392. Fu, J., Jaime Gomez-Hernandez, J., 2009. Uncertainty assessment and data worth in groundwater flow and mass transport modeling using a blocking Markov chain Monte Carlo method. Journal of Hydrology 364 (3e4), 328. Guyoanet, D., Come, B., Perrochet, P., Parriaux, A., 1999. Comparing two methods for addressing uncertainty in risk assessments. Journal of Environmental Engineering 125 (7), 660.

S. Mansour-Rezaei et al. / Journal of Hydro-environment Research 6 (2012) 311e321 Hassan, A.E., Bekhit, H.M., Chapman, J.B., 2008. Uncertainty assessment of a stochastic groundwater flow model using GLUE analysis. Journal of Hydrology 362 (1e2), 89. Karanki, D.R., Kushwaha, H.S., Verma, A.K., Ajit, S., 2009. Uncertainty analysis based on probability bounds (P-box) approach in probabilistic safety assessment. Risk Analysis 29 (5), 662. Klotz, D., 1980. Dispersivity and velocity relationship from laboratory and field experiments. Journal of Hydrology 45, 169e184. Li, J., 2003. Development of an Inexact Environmental Modeling System for the Management of Petroleum-contaminated Sites. University of Regina. Mannina, G., Viviani, G., 2010. An urban drainage stormwater quality model: model development and uncertainty quantification. Journal of Hydrology 381 (3e4), 248e265. Masky, S., 2004. Modelling Uncertainty in Flood Forecasting Systems. A.A. Balkema. Mohamed, M.M.A., Hatfield, K., Hassan, A.E., 2006. Monte Carlo evaluation of microbial-mediated contaminant reactions in heterogeneous aquifers. Advances in Water Resources 29 (8), 1123. Nair, R., Sunny, F., Manikandan, S., 2009. Modelling of decay chain transport in groundwater from uranium tailings ponds. Applied Mathematical Modelling. Sander, P., Bergback, B., Oberg, T., 2006. Uncertain numbers and uncertainty in the selection of input distributions e consequences for a probabilistic risk assessment of contaminated land. Risk Analysis 26 (5), 1363. Schuol, J., Abbaspour, K.C., Srinivasan, R., Yang, H., 2008. Estimation of freshwater availability in the West African sub-continent using the SWAT hydrologic model. Journal of Hydrology 352 (1e2), 30. Sentz, K., Ferson, S., 2002. Combination of evidence in DempstereShafer theory. Report no.SAND2002, p. 835. Simonovic, S.P., 1997. Risk in sustainable water resources management. IAHS-AISH Publication, vol. 240, p. 3.

321

Tucker, W.T., Ferson, S., 2003. Probability Bounds Analysis in Environmental Risk Assessment. Applied Biomathematics, Setauket, New York. http:// www.ramas.com/pbawhite.pdf. Tung, Y.K., 2004. Uncertainty and reliability analysis. Water Resources Systems Management Tools 1. van Griensven, A., Meixner, T., 2006. Methods to quantify and identify the sources of uncertainty for river basin water quality models. Water Science and Technology 53 (1), 51. Wang, S., Jaffe, P.R., Li, G., Wang, S.W., Rabitz, H.A., 2003. Simulating bioremediation of uranium-contaminated aquifers; uncertainty assessment of model parameters. Journal of Contaminant Hydrology 64 (3e4), 283. Willems, P., 2008. Quantification and relative comparison of different types of uncertainties in sewer water quality modeling. Water Research 42 (13), 3539e3551. Yager, R.R., 1986. Arithmetic and other operations on DempstereShafer structures. International Journal of ManeMachine Studies 25 (4), 357e366. Yang, J., Reichert, P., Abbaspour, K.C., Xia, J., Yang, H., 2008. Comparing uncertainty analysis techniques for a SWAT application to the Chaohe basin in china. Journal of Hydrology 358 (1e2), 1. Yuan, C., Druzdzel, M.J., 2006. Importance sampling algorithms for Bayesian networks: principles and performance. Mathematical and Computer Modelling 43 (9e10), 1189. Zadeh, L.A., 1965. Fuzzy sets. Information and Control 8 (3), 338. Zhang, K., Li, H., Achari, G., 2009. Fuzzy-stochastic characterization of site uncertainty and variability in groundwater flow and contaminant transport through a heterogeneous aquifer. Journal of Contaminant Hydrology 106 (1e2), 73. Zimmermann, H.J., 2001. Fuzzy Set Theory e and its Applications. Kluwer Academic Pub.