CONTAMINANT INTRUSION IN WATER DISTRIBUTION SYSTEMS: ADVANCED MODELLING APPROACHES by Saheb Mansour Rezaei Fumani B.A.Sc., University of Tehran, Iran, 2006 M.A.Sc., Sharif University of Technology, Iran, 2009 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE COLLEGE OF GRADUATE STUDIES (Civil Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Okanagan) March 2013 © Saheb Mansour Rezaei Fumani, 2013 i Abstract Since exposure to contaminants may have direct adverse impacts on public health, contaminant intrusion has been recognized as one of the top priority in drinking water supply research. Three components must exist to cause contaminant intrusion into a water distribution system. These include the availability of source(s) of contaminant(s) around a water distribution system, the existence of driving forces (low/negative pressure) to make a contaminant enter into a water distribution system, and the presence of pathway(s) through which contaminant(s) intrude into a water distribution system (WDS). Exposure assessment is the most challenging part as location of contaminant intrusion, rate of intrusion, and the fate of contaminants within WDS need to be estimated accurately. In this dissertation, first, common uncertainty analysis techniques are discussed in the context of conservativeness, execution time, ease of formulation, and complexity. Second, a fuzzyrule based model has been developed to identify contaminant intrusion potential in a WDS. The potential of contaminant intrusion has been determined by integrating the potentials for contaminant sources existence, driving forces, and pathways. Third, a novel ingress model has been developed for more realistic estimation of intrusion rate by taking into account the effects of surrounding soil on intrusion rate. Coupled with an Eulerian-based transient hydraulic model, a Lagrangian transient water quality model is developed to predict the fate of the contaminant throughout a WDS. The proposed models are applied to case studies available in the literature to investigate the applicability of the models. The proposed models enhance the reliability and safety of WDSs by improving the prediction ability of the existing modelling tools. ii Preface The core contribution of this thesis has been structured in publication format. The materials of this dissertation are essentially based on already published journal papers as well as those submitted (under review) for publications. Dr. Bahman Naser supervised this dissertation and is a co-author of all publications. Uncertainty analysis in subsurface contaminant transport described in Chapter 3, was published as a journal article (Mansour-Rezaei et al., 2011a) and a conference proceeding (Mansour-Rezaei et al., 2010). Collaborating with Dr. Rehan Sadiq, I was responsible for developing the methodology, writing computer program, simulations and writing the articles. Dr. Sadiq provided review and Dr. Naser finalized the article. Potential of contaminant intrusion in water distribution systems, described in Chapter 4, has been sent for publication as a journal article (Mansour-Rezaei et al., 2011b). This chapter is also based on the collaboration with Dr. Rehan Sadiq. I developed the model, simulations, and writing the article. While Dr. Sadiq reviewed the paper and provided some comments, Dr. Naser finalized it. The ingress model, described in Chapter 5, was published as a journal article (MansourRezaei and Naser, 2012a). Also a part of that was published in a conference proceeding (Mansour-Rezaei and Naser, 2012c). I was responsible for developing methodology, design of experiments, data analysis, and writing the article. Dr. Naser finalized the work by his thorough comments. Water distribution system modeling, Chapter 6, has been accepted for publication as a journal article and is to be published in June 2013 (Mansour-Rezaei et al., 2012b). Also, parts of this chapter were published in two conference proceeding (Mansour-Rezaei and Naser, 2012d and Mansour-Rezaei et al., 2011c). This chapter is based on the collaboration work with Mr. Ahmad Malekpour and Dr. Brayan Karney. I designed the algorithm, programmed the water quality model and wrote the article. Mr. Malekpour provided the hydraulic model. Dr. Karney reviewed the article and provided some comments, while Dr. Naser finalized it. iii Publications Journal Papers Mansour-Rezaei, S., Naser, Gh., and Sadiq, R., (2011a). "A Comparison of Various Uncertainty Models: an Example of Subsurface Contaminant Transport." Journal of Hydro-Environment Research (http://dx.doi.org/10.1016/j.jher.2011.09.001). Mansour-Rezaei, S., Naser, Gh., and Sadiq, R., (2011b), “Potential of Contaminant Intrusion and Propagation in Water Distribution Systems: A Fuzzy-Based Approach”. Submitted for publication in Urban Water Journal, Under Review, Manuscript ID NURW-20110085. Mansour-Rezaei, S. and Naser, Gh., (2013), "Contaminant Intrusion in Water Distribution Systems: An Ingress Model", Journal-AWWA, Vol. 105, No. 1. Mansour-Rezaei, S., Naser, Gh., Malekpour, A., and Karney B. W., (2013), “Contaminant Intrusion in Water Distribution Systems”, Accepted for publication in Journal-AWWA, Paper No. 08232012-JAWWA0104, To be published in June 2013. Peer-Reviewed Conference Proceedings Mansour-Rezaei, S., Naser, Gh., Sadiq, R., (2010), “Model-Based Uncertainty and Its Propagation in Water Resources Systems”, Water Distribution System Analysis, WDSA, Tucson, Arizona, USA, September 12-15. Mansour-Rezaei, S., Naser, Gh., Sadiq, R., (2011c), “Uncertainty Analysis of Contaminant Propagation in Water Distribution Systems”, Canadian Society for Civil Engineering, CSCE, Ottawa, Ontario, Canada, June 14-17. Mansour-Rezaei, S., Naser, Gh., (2012a), “Contaminant Intrusion in Water Distribution Systems: An Ingress Model”, Environmental and Water Resources Institute, EWRI, Albuquerque, New Mexico, USA, May 20-24. Mansour-Rezaei, S., Naser, Gh., (2012b), “A Lagrangian-Transient Water Quality Model for Contaminant Intrusion in Water Distribution Systems”, Water Distribution System Analysis, WDSA, Adelaide, South Australia, September 24-27. iv Table of Contents Abstract ..................................................................................................................................... ii Preface...................................................................................................................................... iii Publications .............................................................................................................................. iv Table of Contents ...................................................................................................................... v List of Tables ......................................................................................................................... viii List of Figures ........................................................................................................................... x Acknowledgments.................................................................................................................. xiii Dedication .............................................................................................................................. xiv Chapter 1: Introduction ............................................................................................................ .1 1.1 Significance ...................................................................................................................... 1 1.2 Research Objectives ......................................................................................................... 2 1.3 Structure of Thesis............................................................................................................ 4 Chapter 2: Literature Review .................................................................................................... 6 2.1 Contaminant Intrusion Components................................................................................. 6 2.1.1 Contaminant source ................................................................................................... 7 2.1.2 Driving Force............................................................................................................. 9 2.1.3 Pathway ..................................................................................................................... 9 2.2 Uncertainty ..................................................................................................................... 10 2.2.1 Uncertainty Analysis in Water Resources Engineering .......................................... 11 2.3 Likelihood/Potential of Contaminant Intrusion .............................................................. 12 2.4 Ingress Model ................................................................................................................. 14 2.5 Water Quality Model ...................................................................................................... 15 Chapter 3: Uncertainty Analysis in Subsurface Contaminant Transport ............................... 18 3.1 Contaminant Transport in Groundwater......................................................................... 18 3.1.1 Governing Equation................................................................................................. 19 3.1.2 Boundary Conditions ............................................................................................... 20 v 3.1.3 Solution Approaches ............................................................................................... 20 3.1.4 Case Study ............................................................................................................... 21 3.2 Uncertainty Models ........................................................................................................ 22 3.2.1 Fuzzy Set Theory ..................................................................................................... 22 3.2.2 Monte Carlo Simulations ......................................................................................... 26 3.2.3 Probability Box ........................................................................................................ 33 3.3 Comparison of Uncertainty Models ............................................................................... 34 Chapter 4: Potential of Contaminant Intrusion in Water Distribution Systems...................... 38 4.1 Uncertainty Analysis for Contaminant Intrusion ........................................................... 39 4.2 Contaminant Intrusion Potential..................................................................................... 43 4.2.1 Potential for Contaminant Source ........................................................................... 43 4.2.2 Potential for Driving Force ...................................................................................... 44 4.2.3 Potential of Pathway ................................................................................................ 47 4.2.4 Intrusion potential .................................................................................................... 48 4.3 Case Study ...................................................................................................................... 50 4.4 Results and Discussion ................................................................................................... 51 Chapter 5: Ingress Model ........................................................................................................ 56 5.1 Modelling Approach ...................................................................................................... 57 5.1.1 Model Assumptions ................................................................................................. 57 5.1.2 Governing Equations ............................................................................................... 58 5.2 Effective Factors ............................................................................................................ 60 5.3 Dimensional Analysis .................................................................................................... 65 5.4 Test Cases ...................................................................................................................... 66 5.5 Results and Discussion .................................................................................................. 67 Chapter 6: Lagrangian-Eulerian Transient Water Quality Model ......................................... 75 6.1 Lagrangian-Eulerian Transient Model .......................................................................... 76 6.1.1 Eulerian-Based Hydraulic Model ............................................................................ 76 6.1.2 Ingress Model .......................................................................................................... 77 6.1.3 Lagrangian-Based Contaminant Transport Model .................................................. 78 6.2 Hazard and Vulnerability indices .................................................................................. 81 vi 6.3 Model Applications, Results and Discussion ................................................................ 82 Chapter 7: Conclusion, Contributions, and Recommendations for Future Research .......... 100 7.1 Conclusion ................................................................................................................... 100 7.2 Contributions ............................................................................................................... 101 7.4 Limitation .................................................................................................................... 101 7.4 Recommendations for Future Research ...................................................................... 102 Reference ............................................................................................................................. 103 vii List of Tables Table 2.1- Recent examples of uncertainty analyses in water resources engineering ............ 13 Table 2.2- Summery of the existing water quality models ..................................................... 17 Table 3.1- Required data in different uncertainty models ...................................................... 29 Table 3.2- Matrix of interval focal elements .......................................................................... 29 Table 3.3– 5th, 50th and 95th percentiles obtained in different uncertainty models ................. 35 Table 3.4– Comparison of uncertainty models based on various qualitative attributes ......... 35 Table 4.1- Fuzzy linguistic sets for causal and effect concepts .............................................. 45 Table 4.2- Rule set for potential of contaminant source ......................................................... 46 Table 4.3- Potential of contaminant intrusion rule set ............................................................ 49 Table 4.4- Characteristics of the pipes in the WDS case study .............................................. 51 Table 4.5- Potential of contaminant intrusion in the WDS ..................................................... 55 Table 5.1- (-1) and (+1) levels of factors ................................................................................ 61 Table 5.2- Two level factorial design matrix with results ...................................................... 62 Table 5.3- Characteristics of tests were conducted for each constant porosity (ϵ) ................. 68 Table 5.4- Geometry of the developed models ....................................................................... 68 Table 5.5- The optimized coefficient in Equation (5.12)........................................................ 70 Table 5.6- The constant values in Equation (5.14) ................................................................. 71 Table 5.7- Characteristics of tests were conducted for validation of proposed model ........... 72 Table 6.1- Pipes characteristics for case study 1 .................................................................... 83 Table 6.2- Hazard and vulnerability evaluation of case study 1 ............................................. 87 Table 6.3- Pipes characteristics for case study 2 .................................................................... 90 Table 6.4- Hazard and vulnerability evaluation for case study 2 ........................................... 93 Table 6.5- The lowest, most likely and highest demand at junctions ..................................... 96 viii Table 6.6- Fuzzy analysis for maximum contaminant concentrations at junctions 4 and 5 ... 97 ix List of Figures Fig. 1.1- Proposed framework for the dissertation ................................................................... 4 Fig. 2.1- A conceptual structure for contaminant intrusion in a WDS ..................................... 8 Fig. 3.1- Continuous contaminant leaking into shallow aquifer- Observation location at x =1220 m at two different times (400 days and 800 days) ................... 20 Fig. 3.2- Example of a triangular fuzzy set ............................................................................. 23 Fig. 3.3- Analytical and numerical results at location of x =1220 m after 400 and 800 days ................................................................................................................... 25 Fig. 3.4- Upper, lower and most likely bounds of concentration at x = 1220 m (Computed by finite-difference method) ................................................................. 26 Fig. 3.5- Histogram obtained from Monte Carlo simulation (1000 iterations) ....................... 28 Fig. 3.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) ................................. 28 Fig. 3.7- Procedural flowchart for 2-D Monte Carlo simulation ............................................ 31 Fig. 3.8- A comparison between 2-D Monte Carlo simulation and fuzzy set results ............. 32 Fig. 3.9- A comparison of P-Box and fuzzy set results at location of x =1220 m and at time t = 400 days .................................................................................................. 34 Fig. 4.1- Two causal concepts connecting to an effect concept (Sadiq et al., 2009) .............. 40 Fig. 4.2- Fuzzy rule-based model: making inference using two ‘causal factors’ (Sadiq et al., 2009) ................................................................................................... 42 Fig. 4.3- Peseudocolor plot of the proposed fuzzy rule-based model ..................................... 45 Fig. 4.4- Contamination zone underneath a contaminant source (Vairavamoorthy et al., 2007a) ................................................................................ 45 x Fig. 4.5- Fuzzy rule-based model for potential of contaminant intrusion .............................. 48 Fig. 4.6- Membership functions of the inputs ......................................................................... 50 Fig. 4.7- Membership functions of the output ........................................................................ 50 Fig. 4.8- Case study: WDS layout .......................................................................................... 52 Fig. 4.9- Variations of head with time for junctions 1, 2, 3, 4, and 5 ..................................... 53 Fig. 4.10- Variations of (dext+dint)/t with pipe age for pipe 5 (top) and pipe 9 (bottom) ........ 54 Fig. 5.1- A two-dimensional (2D) finite element model to simulate fluid flow through porous media ............................................................................................................ 60 Fig. 5.2- Normal probability plot of factors and interactions effects...................................... 63 Fig. 5.3- Main effects of factors.............................................................................................. 64 Fig. 5.4- Interactions of P, ϵ, Ds, Do, Dp, h on intrusion. The straight black line indicates the effect of low values (-1) and that a red dashed line indicates the effect of high values (1) ..................................................................................... 65 Fig. 5.5- Pressure contour (Pa) around the pipe during occurrence of negative pressure ...... 69 Fig. 5.6- Velocity contour (m/s) (a) and velocity vector (b) at leak-point.............................. 69 Fig. 5.7- Transient pressure induced by valve closure (Collins et al., 2011).......................... 73 Fig. 6.1- Flowchart of the Lagrangian-based transient WQM ................................................ 76 Fig. 6.2- Lagrangian contaminant transport model ................................................................. 79 Fig. 6.3- Contaminant parcel reaching point s at time t1 ........................................................ 81 Fig. 6.4- Schematic view of the network for case study 1 ...................................................... 83 Fig. 6.5- Time variations of head at junctions 1, 2, 3, 4, and 5 of case study......................... 85 Fig. 6.6- Variations of head and concentration at junction 3 of case study 1 ......................... 86 Fig. 6.7- Variations of concentration at junctions 1, 4, and 5 due to the intrusion at junction 3 for case study 1 ....................................................................................... 88 xi Fig. 6.8- Schematic view of the water network for case study 2 ............................................ 91 Fig. 6.9- Variations of head at junctions 2, 6, 11, and 14 in case study 2 .............................. 94 Fig. 6.10- Variations of concentration at junctions 2, 5, 9, 11, and 14 for case study 2......... 95 Fig. 6.11- Time variations of contaminant concentration at junction 4 (top) and 5 (bottom) due to intrusion at junction 1 .................................................................. 98 Fig. 6.12- Time variations of contaminant concentration at junction 4 (top) and 5 (bottom)- 26 numbers of simulations at α-level=0 ................................................. 98 xii Acknowledgments I wish to express my gratitude to all those who provided me with the opportunity to complete my dissertation. I would like to thank my supervisor. Dr. Bahman Naser. His friendliness, advise, constant support, inspiration, kind guidance, and patience are a key basis for my PhD. Thank you for everything you have done for me. I would like to thank the members of my doctoral committee; Dr. Rehan Sadiq, Dr. Kasun Hewage, Dr. Abbas Milani, and Dr. Yves Filion for their accessibility, constructive suggestions, and recommendations. My deepest thanks go to my family. I do not have suitable words to fully describe their endless support through my life. Their support and encouragement made this dissertation possible. I also appreciate all financial support for my PhD study, including the Okanagan Basin Water Board, NSERC General Revenue Fund, University of British Columbia PhD Tuition Fee Okanagan Award, and the University of British Columbia International Partial Tuition Okanagan Scholarship. xiii Dedication This thesis is dedicated to the memory of my Dad, Karim Mansour Rezaei Fumani, who passed away during the time I was writing my thesis. He will always be loved and shall forever be alive in my heart and mind. xiv Chapter 1: Introduction “ … Water is essential for all dimensions of life. Over the past few decades, use of water has increased, and in many places water availability is falling to crisis levels. More than eighty countries, with forty percent of the world’s population, are already facing water shortages, while by year 2020 the world’s population will double. The costs of water infrastructure have risen dramatically. The quality of water in rivers and underground has deteriorated, due to pollution by waste and contaminants from cities, industry and agriculture. Ecosystems are being destroyed, sometimes permanently. Over one billion people lack safe water, and three billion lack sanitation; eighty per cent of infectious diseases are waterborne, killing millions of children each year …”. World Bank Institute, Water Policy Reform Program, November 1999 As the above quotation implies, fresh-drinking water reserves are scarce from both quantity as well as quality points of view. Thus, ensuring a safe and reliable water distribution system (WDS) to deliver water from trusted sources to the consumers is a top priority for all water authorities. This exigency highlights the quest for a comprehensive understanding of the hydraulics of the flow and the safety of the delivered water. Such understanding would be useful not only for designing and operating a system, but also for serving as an early warning mechanism for any possible mishaps. 1.1 Significance A WDS is a complex and dynamic environment that over its life-time, water quality degradation may occur due to variety of physical, chemical, and biological processes. Water quality disturbances may occur in a WDS anywhere from water source(s) to transmission main(s) and from treatment plant(s) to the users. Although many municipalities monitor water quality at sources and treatment plants (National Research Council, NRC, 2004), outbreaks sometimes occur. Payment et al. (1997) reported an excess of gastrointestinal illnesses in 14–40% of consumers of drinking water from Canadian WDSs. Schuster et al. (2005) reported 288 disease outbreaks resulting from Canadian WDSs between 1974- 2001. 1 While at least 8000 cases of illness were linked to these outbreaks, the real number could be 10 to 1000 times more depending upon the severity of the symptoms (http://www.ec.gc.ca/inre-nwri/default.asp?lang=En&n=235D11EB-1&offset=2&toc=show). Research indicated viruses in water-samples collected from household taps (Borcharddt et al., 2009) possibly entered the system at distribution level (Besner et al., 2011). The 1993outbreak in Milwaukee resulted in 69 deaths and over 400,000 cases of illness at an estimated total cost of $96.2 million (Corso et al., 2003). These, along with the infection outbreaks in Walkerton in 2000 and North Battleford in 2001, suggest that the consequences of an outbreak could be catastrophic even if it does not occur frequently. The former resulted in 7 deaths, 2300 sick, and $64.5 million cost, while the later made 40–50% of residents become sick (Woo and Vicente, 2003). Water quality failures are likely to occur more frequently due to the continued deterioration of ageing WDSs and the impacts of climate change and urbanization on water resources and emerging contaminants. Additionally, an unrealistic model of a WDS could cause delay in detection of hydraulic and water quality failure (if detected at all) by providing inaccurate estimates of flow parameters. This delay could eventually result in more economic loss, human illness or death. 1.2 Research Objectives Contaminant intrusion in a WDS is an important mechanism that may result in the system’s failure from water quality perspective. Three prerequisite conditions of contaminant intrusion are contaminant sources, driving forces (low/negative pressure), and pathways. Different contaminants may exist immediately external to WDS. A contaminant source near a WDS can be the contaminated soil/water containing the chemical and/or microbial contaminants. The service connections or structurally deficient locations (such as submerged air valves or leak-points) provide pathways through which a contaminant may intrude a WDS. Driving forces can be induced by an increase in external pressure or a decrease in internal pressure. As any WDS has leakage points that can provide entry pathways for intrusion under low/negative pressure conditions, a WDS can be vulnerable to intrusion and subsequently drinking water contamination. Thus, the prime objective of this dissertation is to enhance the reliability of WDSs by improving the prediction ability of the existing modelling tools. This will ensure safe and 2 reliable supply of drinking water that eventually increases consumers’ confidence. This will be achieved through proposing models to identify, characterize, and quantify potential of contaminant intrusion into a WDS and propagation of the contaminant spatially and temporally. Such models will assist the authorities in making informed decisions that will help reducing the likelihood of a contaminant intrusion occurrence and the related consequences (in the case of occurrence). While Fig. 1.1 indicates the framework for this dissertation, the main objective will be gathered through a set of sub-objectives including: 1. Potential for Contaminant Intrusion: A conceptual model is proposed to identify susceptible locations of WDS based on the potential for contaminant intrusion. The goal is not to develop a sophisticated and comprehensive tool, but rather it is to create a framework to model systems for which only vague information is available. 2. Ingress Model: This is to estimate the rate of a contaminant intrusion to a WDS more realistically. A model is created to consider the effects of surrounding soil (including porosity and size of soil particles) for more realistic estimate of intrusion rates into a pipe through a leak-point. 3. WDS Modeling: This is to investigate the fate of a contaminant in a WDS by developing a transient water quality model (WQM). The WQM will be obtained by coupling a Lagrangian contaminant transport model with an Eulerian transient hydraulic model (for detecting junctions of a WDS experiencing negative pressures). 4. Case Studies: The proposed models are employed to simulate intrusion in two wellknown WDSs reported in the literature. Hazard and vulnerability of the systems will be evaluated to show the applicability of the proposed WQM. 3 Potential for Driving Force Objective 1 Potential for Contaminant Intrusion Potential for Contaminant Source Potential for Pathway Objective 2 Ingress Model Objective 3 Water Distribution System Modelling Objective 4 Case Study Rate of Intrusion Hydraulic Model Contaminant Fate Hazard Analysis Vulnerability Analysis Fig. 1.1 – Proposed framework for the dissertation. 1.3 Structure of Thesis The structure of this dissertation is generally based on published/ submitted journal and conference papers. It includes four chapters with a small degree of overlap. • Chapter 2 reviews the available literature related to the contaminant intrusion. • Chapter 3 discusses and compares the commonly-used uncertainty analysis techniques including fuzzy set theory, 1-D and 2-D Monte Carlo simulations, and P-Box method. • Chapter 4 applies a fuzzy-logic based model to identify potential of contaminant intrusion in a WDS. The potential of existence of contaminant sources, driving forces, 4 and pathway existence are integrated to determine the potential of a WDS to a contaminant intrusion. • Chapter 5 proposes an ingress model to estimate intrusion flow rate more realistically. A model is proposed to consider the characteristics of the effects of surrounding soil on intrusion flow rates by conducting computer experiments, factorial design, dimensional analysis, and regression analysis. • Chapter 6 develops a Lagrangian-based transient WQM to simulate the fate of an intruded contaminant along a WDS. A Lagrangian contaminant transport model is coupled with an Eulerian transient hydraulic model to accurately simulate propagation of intruded contaminant. • Chapter 7 concludes the dissertation by discussing the key results and contributions. The Chapter also suggests some improvements for future work. 5 Chapter 2: Literature Review This chapter reviews contaminant intrusion and its required prerequisite components including contaminant source, driving force, and pathway. The chapter discusses the techniques commonly used in water resources engineering for uncertainty analysis, and likelihood of contaminant intrusion in a WDS. The chapter also introduces the existing ingress and water quality models for the simulation of contaminant intrusion and its propagation along a WDS. 2.1 Contaminant Intrusion Components The American Water Works Association (AWWA) highlighted issues of metal accumulation and the risk of release from pipes, contaminant intrusion, nitrification and its control, and disinfection-by-products (DBP) as water quality concerns in future (MacPhee, 2005). Contaminant intrusion is one of the most common causes of water quality failure in a WDS. Because exposure to a contaminant intrusion can have directly adverse impacts on public health, the Research and Information Collection Partnership consisting of United States Environmental Protection Agency (US EPA) and Water Research Foundation (WRF) has recognized contaminant intrusion as one of the ten top priority research areas (http://www.waterrf.org/Research/ResearchPrograms/Partnership/Pages/RICP.aspx). While recognized intrusion events are not well documented (Besner et al., 2011), a low/negative transient pressure inside a pipe allows the entry of contaminants to the pipe if pathway(s) and contaminant source(s) exist. Geldreich et al. (1992), Payment et al. (1997) and Blackburn et al. (2004) indicated outbreaks possibly related to intrusion at WDS level. Suggesting contaminant intrusion as a possible mechanism of WDS contamination, NRC (2006) highlighted issues of low disinfectant residual and transient pressure in terms of their potential health risk for the observed gastrointestinal illnesses in the tap-water consumers. Kirmeyer et al. (2001) provided a detailed list of causes for low/negative pressure events in the normal operation of a WDS. Pressure monitoring studies (Kirmeyer et al., 2001; Gullick et al., 2004 and 2005; Hooper et al., 2006; and Besner et al., 2009) indicated low/negative pressures in WDSs. Besner et al. (2011) discussed two types of transient low/negative pressure events in a WDS: 1) transient events lasting from milliseconds to minutes and 2) 6 sustained events with durations in the range of few minutes to even hours. Events such as the United States northeast blackout of 14/9/2003 and watermain breaks and subsequent repair site isolation can cause sustained events (Besner et al., 2011). Lindley (2001) indicated adverse pressure gradient as a susceptibility condition for contaminant intrusion in 62.8% of the outbreaks in USA during 1971-1998. Although he did not claim it, an adverse pressure gradient can be attributed to transient flow. Three components must exist to cause contaminant intrusion into a WDS. The first is the availability of source(s) of contaminant(s) in the proximity of a WDS. The second is the existence of a driving force (low/negative pressure) to make a contaminant enter a WDS. The third is pathway for a contaminant to intrude into a WDS. Fig. 2.1 shows a detailed conceptual structure describing the three essential elements of a contaminant intrusion. The following sections review the literature related to the three components of contaminant intrusion. 2.1.1. Contaminant Source Contaminant sources in urban areas include sewer lines, garbage disposal areas, leaking underground petroleum storage units and so on. Chemical and microbial contaminants are often found around buried water mains. Schuster et al. (2005) analyzed 288 waterborne outbreaks in Canada occurring during 1974-2001. They identified a variety of pathogens in the outbreaks. They also indicated severe weather, close proximity to animal population, treatment system malfunctions, and poor maintenance and treatment practices as the contributing factors to these outbreaks. Karim et al. (2003) reported microbial contaminants in 66 soil and water samples collected from eight WDS-utilities in six states in USA. They observed different types of microbial contaminants including total coliforms, fecal coliforms, Clostridium, Bacillus, viruses, and Coliphage. 7 Fig. 2.1- A conceptual structure for contaminant intrusion in a WDS. 8 2.1.2. Driving Force A driving force (e.g., low/negative pressure) can be induced by a pump operation, a pipe replacement, a valve closure/opening, or demand variations. Lindley (2001) illustrated the relative frequency of deficiencies attributed to the waterborne outbreaks in the USA. In 62.8% of the events during 1971-1998, susceptibility conditions met adverse pressure. Besner et al. (2011) monitored the frequency and magnitude of negative pressures at 19 sites of the WDS in a suburban area of Montreal (Canada) from June 2006 to November 2007. They recorded 36 negative pressure events in the system; 55% of them lasted less than 3 minutes. Gullick et al. (2004) also carried out pressure monitoring for 43 sites in 8 WDSs. Over a period of ~4640 days, they reported 21 negative pressures that lasted less than 3 minutes mainly caused by sudden pump shutdowns. Gullick et al. (2004) recorded negative pressures as low as -10 meters, while Besner et al. (2011) reported a minimum pressure of -5 meters. In addition to the changes in operating conditions and boundary devices, demand variations in a WDS may also lead to a low/negative pressure event. In 1980, in northern Georgia, 1500 people were affected by an outbreak of viral gastroenteritis. The evidence indicated that the low pressure in the WDS during the peak demand caused contaminant intrusion from an industrial system (Kaplan et al., 1982). Kirmeyer et al. (2001) observed pressure fluctuations possibly induced due to the changes in water use patterns (LeChevallier et al., 2003). Besner et al. (2011) also reported some negative pressure events lasted from 13 to 1719 seconds caused by unknown sources (possibly due to rapid demand variations). 2.1.3. Pathway A contaminant may intrude into a WDS through a variety of pathways including submerged air valves, leak points, repair and installations, faulty seals, joints, and service connections (Thomason and Wang, 2009). Internal/external corrosion reduces the hydraulic capacity of the pipe and results in mechanical failure (e.g., pipe breaks). Approximately 75% of WDSs in the USA consist of various forms of ferrous pipes (Thomason and Wang, 2009). About 50% of the pipes in Canada are old gray cast-iron (Rajani and McDonald, 1995). Every day, an average of 700 water main breaks occurs in Canada and the USA (http://www.watermainbreakclock.com/breakclock.htm). Pipe 9 breakage occurs when the total internal and external stresses acting on a pipe is more than the pipe strength, which is reduced by internal and external corrosion (Yaminighaeshi, 2009). Sadiq et al. (2004) studied the failure of gray cast iron (CI) pipes due to the external corrosion. Parameters of external corrosion model, stresses acting on a pipe, and mechanical characteristics of the pipe were considered as uncertain parameters. Yamini and Lence (2010) defined the rate of internal corrosion of a CI pipe as a function of residual chlorine concentration, the chlorine decay constant, and velocity of the flowing water in a pipe segment with a given length and diameters. They conducted a probabilistic analysis to evaluate uncertainties related to the input parameters of the model. 2.2 Uncertainty Uncertainties can arise from various sources including 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 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 10 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). 2.2.1 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 2.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 2.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 twodimensional (1-D and 2-D) Monte Carlo simulations, Bayesian and Probability-Box (PBox) methods. There is still no clear and direct study comparing the performance and robustness of these techniques. This dissertation provides a comparison among commonly used uncertainty analysis techniques (i.e. fuzzy set theory, Monte Carlo simulations, Bayesian and Probability-Box (P-Box) methods). 11 2.3 Likelihood/Potential of Contaminant Intrusion Contaminant intrusion is a complex phenomenon for which the data are generally uncertain (Sadiq et al., 2007). It is very challenging (if not impossible) and often expensive to reliably and accurately provide all the required data for a real-life system. Uncertainties related to the contaminant intrusion events can arise from data, structural, and parameter uncertainties. For example, to accurately estimate the contaminant potential around a WDS a model is to be employed for soil around the system. However, choosing the right model is not a straightforward task. The early question to select an appropriate model is about what type of model is the best representative of reality? A decision has to be made among variety of models including three-, two- and one-dimensional (3-D, 2-D, and 1-D) models, saturated or unsaturated zone, continuous or instantaneous contaminant release, and so on. One should have a reliable set of detailed information about the study area (e.g. ground water level, aquifer depth, boundary condition, soil permeability, and dispersivity and so). As well, the tradeoff between the simplicity and accuracy is another criterion for selecting an appropriate model. Even if an appropriate model is selected for groundwater flow through soil, there are still new questions to be answered with regards to the value of the model parameters. To define the likelihood/potential of a water quality failure due to a contaminant intrusion event, different methods have been proposed in the literature including evidential reasoning (Sadiq et al. 2006), and fault tree analysis (Sadiq et al., 2008). Evidential reasoning, which is based on Dempster-Shafer theory (DST) is useful to incorporate both aleatory and epistemic uncertainties and it is also able to handle ambiguous and conflicting information in the inference mechanism (Sadiq et al., 2006). Sadiq et al. (2008) evaluated the risk of water quality failure in WDSs by using fault tree analysis. Fault tree was composed of some basic events and the logical relationships (and and or) among them. This dissertation estimates potential of intrusion in cast-iron WDSs by employing a fuzzy-rule based model (to identify the potential of contaminant source), internal and external corrosion models (to assess the potential for pathways), and a transient hydraulic model (to study the potential for driving forces). 12 Table 2.1- Recent examples of uncertainty analyses in water resources engineering. Group Analytical Techniques (Tung 2004) Techniques Target Uncertain Parameters • Hydrologic Parameters • Hydraulic Parameters • Derived Distribution Technique • Fourier Transform Technique • Laplace and Exponential Transform Techniques • Mellin Transform Technique • Water Resources • Climate Change • First-Order Variance Estimation Method (Tung, 2004). • Probabilistic Point Estimation Methods (Tung 2004). • Water Resources • Climate Change • Hydrologic Parameters • Hydraulic Parameters • Ground water flow (Balakrishnan et al. 2005) • Uranium decay chain transport in the groundwater (Nair et al. 2009) • Hydraulic conductivity and recharge rate field. • Velocity, dispersivity, porosity, etc. • Hydraulic conductivity • Recharge rate • Biotic and abiotic reaction parameters. Stochastic Response Surface Method (SRSM). High Dimensional Model Representation (HDMR) Approximation Techniques Application Generalized Likelihood Uncertainty Estimation (GLUE). Parameter Solution (ParaSol). Sequential Uncertainty Fitting algorithm. Bayesian inference based on Markov chain Monte Carlo (MCMC). • 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 et al. 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). • 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. 13 Table 2.1 Cont’d. Recent examples of uncertainty analyses in water resources engineering. Group Techniques Application Bayesian inference based on importance sampling (IS) • 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) • Permeability, porosity, hydraulic conductivity, dispersivity 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) Fuzzy Set Theory Approximation Techniques Monte Carlo Simulation Dempster-Shafer Structure (P-Box) • Risk assessment of contaminated land • • Target Uncertain Parameters Outlet discharge Parameters Hydraulic conductivity and porosity. Velocity and dispersivity Hydraulic conductivity 2.4 Ingress Model Exposure assessment is the most challenging part in contaminant intrusion risk assessment (LeChevallier et al., 2011). For exposure assessment, location and duration of contaminant intrusion, contaminant volume, pathogen concentration, and the fate of intruded contaminants along a WDS should be estimated properly. Duration of an intrusion could be small; however, even a small volume of intrusion may have severe consequences on public health. Microbial contaminants even with dilution can cause infections. Some microbes, even with a single organism, could cause an infection (LeChevallier et al., 2003). Thomason and Wang (2009) described different structural failure modes for a cast iron WDS including circumferential breaks, longitudinal breaks, spiral breaks, and through hole, and so on. Although, due to its simplicity, an orifice equation is often suggested for estimating the intrusion rate into a WDS through a leak-point, it may not be necessarily accurate for different shapes and sizes of orifices (Besner et al., 2011). The 14 orifice equation: Q = CDAo(∆P)0.5, assumed the intrusion rate to be directly proportional to the orifice area (Ao) and square root of pressure difference across the orifice (∆P). Besner et al. (2011) argued that due to the pipe material behavior, orifice equation does not necessarily estimate intrusion rates accurately; the orifice area can be changed by variation of pressure. As well, other research has indicated that the exponent value may differ from the theoretical value of 0.5 as the effective leak area may, at least in some cases, be pressure dependent (Lambert, 2001 and Greyvenstein and van Zyl, 2007). Moreover, in the orifice equation, the effects of the surrounding soil and type of flow through porous media (laminar and/or turbulent) are not taken into account, indeed the soil hydraulic conductivity may have significant impacts on the intrusion rate (Collins et al., 2010 and 2011 and Besner et al., 2011). Collins et al. (2010 and 2011) evaluated numerically and experimentally how the soil around a buried pipe affects the intrusion rate. Although they did not propose any model to estimate intrusion rate, the comparison between the results predicted by the orifice equation and those computed by the numerical model showed that the existence of surrounding soil decreases the intrusion rate by up to 50% (Collins et al. 2010). Moreover, it was observed that the presence of an orifice decreases the magnitude of low/negative pressure (Collins et al., 2011). To have more accurate estimation of intrusion rates, this dissertation provides an ingress model considering the surrounding soil properties (porosity and average diameter of soil particles). 2.5 Water Quality Model Without knowledge of true conditions of flow in a WDS, any estimate of an outbreak or a contaminant fate in the system is unreliable regardless of how/where/why it enters the system. The available water quality models (WQMs) largely ignore inertia and compressibility effects and assume either steady or quasi-steady conditions for flow. EPANET (Rossman, 2000) is a typical example. A real WDS may undergo transient conditions in which both fluid inertia and compressibility effects are significant. Steady or quasi-steady models provide an overall understanding of flow in a WDS. Continuous demand variations along with operations of local devices (valves, pumps, etc.) create transients that cannot be captured by steady or quasi-steady models. A realistic model can more reliably define what operating strategy must be undertaken to ensure 15 acceptable water quality in an efficient manner. Research indicated enhanced survival of pathogenic microbes in biofilms (Percival and Walker, 1999; Skraber et al., 2005). The sudden changes in flow create aggressive shear forces in a WDS such that the materials contributing to water quality are removed from pipe-wall or biofilms and mobilized with the flow (Boxall et al., 2003) posing health risks to consumers. Using 1-D assumptions for flow and lumping chlorine decay at bulk- and wall-region together, Fernandes and Karney (2004) studied transient conditions for hydraulics and chlorine transport in a pipe. Given that a 1-D model cannot capture the radial variations of concentration, Naser and Karney (2007 and 2008) developed/compared transient 1-D and 2-D water quality models for flow in a pipe. The 2-D simulation did exhibit more realistic results than its 1-D counterpart when a decay process occurred at pipe wall. Interestingly, however, there were no significant differences between the models when the process occurred in a bulk region. The advantages in flow and water quality representations in a 2-D model were severely compromised due to computational requirements and simulation time. This makes a 2-D simulation computationallyimpractical for analyzing a real-life WDS. Fernandes and Karney (2004) reported a method of characteristic (MOC) as a computationally efficient model for transient hydraulic analysis of a WDS. The literature provides Eulerian, Lagrangian and mixed Eulerian–Lagrangian models for modeling contaminant transport in pipes (Rossman 2000; Basha and Malaeb, 2007; Andrade et al., 2010). Table 2.2 shows the existing WQMs in the literature. The prime assumption in all currently-available Lagrangian or mixed Eulerian–Lagrangian models is that flow velocity remains unchanged (i.e., steady-state hydraulic conditions). Moreover, comparing Eulerian and Lagrangian methods, Rossman and Boulos (1996) argued that Lagrangian methods are more computationally efficient in terms of CPU simulation time and memory usage. They also indicated that Eulerian methods were subject to numerical dissipation and dispersion errors (i.e., amplitude and phase shift error). However, Lagrangian methods were free of numerical dissipation and dispersion errors. To accurately simulate propagation of intruded contaminant, this dissertation provides a Lagrangian-Eulerian transient WQM by coupling an Eulerian transient hydraulic model with a Lagrangian contaminant transport model. 16 Table 2.2 - Summery of the existing water quality models. Contaminant Transport Model Numerical scheme Governing Equation Hydraulic Model Authors • Lagrangian Event-Driven Method Advection-Reaction Equation Steady-State Modelling Boulos et al. (1994) • • • • Eulerian Finite-Difference Method Eulerian Discrete-Volume Method Lagrangian Time-Driven Method Lagrangian Event-Driven Method Advection-Reaction Equation Steady-State Modelling Rossman, and Boulos (1996) • Eulerian-Lagrangian Method Advection-Diffusion -Reaction Equation Steady-State Modelling Tzatchkov et al. (2002) • Eulerian Finite-Difference Method Advection-Reaction Equation Transient Analysis Fernandes and Karney (2004) • Lagrangian Modified Event-Driven Method Advection-Reaction Equation Steady-State Modelling Munavalli and Kumar (2004) • Eulerian- Lagrangian Method Advection-Diffusion -Reaction Equation Steady-State Modelling Basha and Malaeb (2007) Summery This chapter reviewed the available literature on contaminant intrusion. The chapter also introduced commonly used uncertainty analysis techniques, likelihood/potential of contaminant intrusion, ingress models, and water quality models. 17 Chapter 3: Uncertainty Analysis in Subsurface Contaminant Transport A version of this chapter has been published in Journal of Hydro-environment Research. Mansour-Rezaei, S., Naser, Gh., and Sadiq, R. (2011). "A Comparison of Various Uncertainty Models: an Example of Subsurface Contaminant Transport". (http://dx.doi.org/10.1016/j.jher.2011.09.001) Preface 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-chemical, biological and socio-economical changes in the surrounding environment, tremendous efforts have 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 limit human ability to make correct predictions. Therefore, proper treatment of uncertainties is required for modelling hydrosystems. This chapter aims at comparing the results of four common uncertainty models using an example of contaminant transport in groundwater. An advection–dispersion equation (ADE) is employed to describe the transport of a contaminant in groundwater. For simplicity, two parameters - dispersion coefficient and velocity - are 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 are investigated. 3.1 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 18 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 modelling 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. 3.1.1 Governing Equation Fig. 3.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 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: ∂C ∂ 2C ∂C = D 2 −V ∂t ∂x ∂x 3.1 Here, C is the contaminant concentration (ML-3); V is the flow velocity (LT-1) in xdirection; and D is the dispersion coefficient (L2T-1). 19 Fig. 3.1- Continuous contaminant leaking into shallow aquifer- Observation location at x =1220 m at two different times (400 days and 800 days) 3.1.2 Boundary Conditions Regardless of a solution approach, the Equation (3.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 step-wise continuous release of a contaminant with constant concentration Co at x = 0. • Outlet Boundary Condition: This is considered as a no-flux boundary (i.e. ∂C / ∂x = 0 ) at x = l x in which lx is some characteristic length-scale representing the size of the contaminant plume. 3.1.3 Solution Approaches Solving differential equations such as (3.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. 3.1.3.1 Numerical Solution Several approaches have been traditionally used to numerically solve Equation (3.1). These include (but are not limited to) finite element method, finite difference method, 20 finite volume method, and boundary volume method. A central differencing scheme of an implicit finite difference method is used in this study to numerically integrate Equation (3.1). The method is approximated as: [ Cit +1 − C it ] [ C t +1 − 2C it +1 + Cit−+11 + C it+1 − 2Cit + Cit−1 ] = D × i +1 ∆t 2 ∆x 2 [ C t +1 − C it−+11 + Cit+1 − Cit−1 ] − V × i +1 4 ∆x 3.2 where ∆t and ∆x 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 the Equation (3.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. 3.1.3.2 Analytical Solution Given the conditions at inlet/outlet boundaries introduced in section 3.1.2, Equation (3.1) can also be solved analytically providing the following general solution (Dou et al. 1997): C 1 x − Vt 1 x + Vt Vx = erfc + exp erfc C0 2 D 2 Dt 2 2 Dt 3.3 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.3) to verify the validity of the proposed numerical model by a direct comparison of the results. 3.1.4 Case Study This research applies the 1-D contaminant transport equation in a uniform and hydraulically steady groundwater flow (Fig. 3.1) initially introduced by Duo 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 (∆x) of 15.25 m. A time step (∆t) of one day is also applied to numerically 21 integrate the Equation (3.2). Additional information about the case study can be found in Duo 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.2 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. 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 similar to Duo et al. (1997), and for simplicity purpose, this research assumes velocity and dispersion coefficient are mutually independent parameters, although empirical evidence indicated dependency between velocity and dispersion coefficient (Klotz 1980). 3.2.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) = {( x, µ A ( x)), x ∈ X , µ A ( x) ∈[0,1]} 3.4 where X = {x} is a universal set of elements, A(x) is a fuzzy set of X, and µA(x) is the membership grade of x in A. The parameter µA(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 µA(x) = 1. The fuzzy set A is convex if for every real number a, b, and c, µ A (c) ≥ min[µ A (a), µ A (b)] given that a<b<c. This means that the membership function may increase or decrease. The αlevel cut, described as Aα = {x µ A (x) ≥ α }, consists of all components of X whose membership grade is greater than or equal to α. Fig. 3.2 shows the concepts of the α- 22 level cut and support in a fuzzy set. The support of the fuzzy set A is the set Supp(A)={x| µA(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. 3.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 3.1. Membership Grade Most Credible Value 1 0.75 α-Level cut 0.5 0.25 0 Lowest possible value Highest possible value Support Parameter Value Fig. 3.2- Example of a triangular fuzzy set 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 modelling 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 α-level cut concept and interval analysis technique. This method is easy to implement and has been successfully applied 23 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 x1 , x2 ,..., xn inputs through y = f ( x1 , x2 ,..., xn ) . 2. Discretize the range of membership grade [0, 1] into a finite number of values α1 , α 2 ,...,α M . The number of discretization depends on the degree of accuracy in approximation. 3. For each membership value α j , find the corresponding interval for Ai in xi for i = 1, 2, ..., n . These are the supports of the α j -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 and combine them into an n-array. There are 2 n combinations for the vector ( x1 , x2 ,..., xn ) which are obtained from 2 n distinct permutations. 5. Evaluate f ( x1 , x2 ,..., xn ) for each of 2 n vectors and obtain y1 , y2 ,..., yn ( 2 n values for y ). [∧ y k , ∨ y k ] is the desired interval for y , where k k ∧, ∨ show the maximum k k and minimum values of y k , respectively. These points define the support of α j - cut. 6. Repeat the process for other α -level combinations to obtain additional [∧ y k , ∨ y k ] k k After substituting fuzzy velocity and dispersion coefficient parameters in the Equations (3.2) and (3.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 α-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.3 illustrates the membership functions of concentration for both analytical and numerical 24 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. 3.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. 3.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 will help the engineers or other authorities to make more informed decisions when required. 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.3- Analytical fuzzy and numerical fuzzy results at location of x =1220 m after 400 and 800 days It should be emphasized that this research sets the output of α-level in the fuzzy set theory by searching for the upper and lower bounds of the corresponding output range. 25 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. 1.2 Upper Bound 1 Lower Bound C/C0 0.8 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. 3.4- Upper, lower and most likely bounds of concentration at x = 1220 m (Computed by finite-difference method) 3.2.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.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. 26 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 3.1. 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 (3.2) to compute concentration. The procedure is repeated for a large number of times. Fig. 3.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. 3.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. 3.5 was multiplied by a factor of “5.2” such that the maximum relative frequency equals to unity (Guyoanet et al. 1999). Fig. 3.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 27 Fig. 3.5- Histogram obtained from Monte Carlo simulation at t = 400 days (1000 iterations) Fig. 3.6- A comparison between Monte Carlo simulations and fuzzy set results at t = 400 days (To make a fair comparison, the relative frequencies were multiplied by a factor such that the maximum relative frequency equals 1) 28 Table 3.1- Required data in different uncertainty models Velocity (mday-1) Low (L) = 1.8 Most Likely (M) = 2.4 High (H) = 3.1 Low (L) = 1.8 Most Likely (M) = 2.4 High (H) = 3.1 Dispersion coefficient (m2day-1) Low (L) = 3.45 Most Likely (M) = 36.72 High (H) = 61.38 Low (L) = 3.45 Most Likely (M) = 36.72 High (H) = 61.38 2-D Monte Carlo Low (L) (µ,σ) = (1.8, 0.1) Most Likely (M) (µ,σ)= (2.4, 0.1) High (H) (µ,σ) = (3.1, 0.1) Low (L) (µ,σ) = (3.45, 1) Most Likely (M) (µ,σ)= (36.72, 7) High (H) (µ,σ) = (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 Method Fuzzy set theory 1-D Monte Carlo Table 3.2- Matrix of interval focal elements Dispersion coefficient: [1.80-3.10] 0.20 [1.95-2.92] 0.20 Velocity [2.10-2.75] 0.20 [2.25-2.57] 0.20 [2.40-2.40] 0.20 [3.45-61.38]* 0.2 [11.77-55.22] 0.20 [0.000-0.570]** 0.04 [0.000-0.576] 0.04 [0.000-0.445] 0.04 [0.000-0.438] 0.04 [0.000-0.326] 0.04 [0.000-0.313] 0.04 [0.000-0.222] 0.04 [0.002-0.207] 0.04 [0.000-0.141] 0.04 [0.007-0.126] 0.04 [20.09-49.05] 0.20 [0.000-0.574] 0.04 [0.001-0.430] 0.04 [0.002-0.299] 0.04 [0.009-0.191] 0.04 [0.026-0.111] 0.04 [28.40-42.89] 0.20 [0.001-0.571] 0.04 [0.003-0.421] 0.04 [0.008-0.283] 0.04 [0.022-0.172] 0.04 [0.051-0.094] 0.04 [36.72-36.72] 0.2 [0.003-0.570] 0.04 [0.007-0.409] 0.04 [0.017-0.264] 0.04 [0.038-0.151] 0.04 [0.076-0.076] 0.04 * Focal elements and probability mass. ** The maximum and minimum values of contaminant concentration and probability mass 3.2.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. 3.7 provides the 29 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 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~ µ, σ) were used to describe epistemic uncertainty for low (L), most likely (M), and high (H) values of triangular distribution of model parameters. Table 3.1 provides the mean and variance of L, M, and H values for velocity and dispersion coefficient. 2- 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. 3- 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. 4- Step 3 was repeated for sufficient number of times (say, 1000 iterations). Similar to 1-D Monte Carlo, 1000 iterations was found to be sufficient to obtain a smooth and consistent CDF curve. The required measures of aleatory uncertainty were obtained in this step. 5- 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 timeconsuming and computationally demanding. 30 Check for no.of simulation (if n<N2) No yes Sample Aleatory variables V and D are implemented in Model CDF 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<N1) out Put yes Random (M) No No if M>L if H>M Sample Epistemic variables yes yes Triangular Shape PDF Fig. 3.7-Procedural flowchart for 2-D Monte Carlo simulation 31 The results for the contaminant concentration are provided through a family of curves on Fig. 3.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 of PDF's describes epistemic uncertainty. Fig. 3.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, 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 decision-making point of view as it may have substantial environmental impacts. Fig. 3.8- A comparison between 2-D Monte Carlo simulation (red lines) and fuzzy set results (black line) 32 3.2.3 Probability Box Using the Dempster-Shafer structure (Yager 1986), a probability box (P-Box) provides a comprehensive way to describe and propagate uncertainties in a complex modelling 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 3.1. 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.2 shows how P-Box method can be used to determine uncertainties in contaminant transport modelling (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 of the matrix in Table 3.2, respectively. Suppose U(i, j) is an arbitrary cell inside the matrix on Table 3.2 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) shows 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 (3.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.2) are corresponding with lower 33 and upper values of contaminant concentration computed by fuzzy set theory at five different α-levels: 0, 0.25, 0.5, 0.75 and 1. Fig. 3.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. 3.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 distribution is large. Fig. 3.9- A comparison of P-Box and fuzzy set results at location of x =1220 m and at time t = 400 days 3.3 Comparison of Uncertainty Models To have more illustrative comparison between the four different uncertainty models, Table 3.3 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 34 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 concentration, and the 95th percentile exceeds all but 5 percent of the values. Table 3.3 – 5th, 50th and 95th percentiles obtained in different uncertainty models Method 5th (%) 50th (%) 95th (%) 1-D Monte Carlo - 0.065 0.354 2-D Monte Carlo (L, H) Fuzzy Set Theory (L, M, H) P-Box (L, H) (0.000, 0.046,0.092) (0.000, 0.094) (0.014, 0.127) (0.002, 0.151,0.299) (0.003, 0.306) (0.165, 0.546) (0.062, 0.298,0.535) (0.076, 0.576) L: low, M: most likely, H: high. Table 3.4 –Comparison of uncertainty models based on various qualitative attributesThe numbers in the table shows the rank of models. Qualitative attributes Fuzzy Set Theory Conservativeness Computing Time Ease of formulation Interpretation Complexity Use in literature 2 1 1 1 2 1-D Monte Carlo 4 2 3 4 1 2-D Monte Carlo 3 3 4 3 3 Probability Box 1 1 2 2 4 Table 3.4 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 are 53.67, 708.34, 1.079E5, and 182.59 seconds, respectively. While the 2-D Monte Carlo simulation is the most time consuming method, the P-Box approach as well as fuzzy set are 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 35 P-Box methods to model uncertainty related to parameters of hydraulic model is easier than Monte Carlo simulations, which need to select random parameters. While it is frequently used by researchers to model uncertainty, but Monte Carlo simulations are much more complex than P-Box method and fuzzy set theory. Moreover, P-Box method and fuzzy set theory provide more conservative upper and lower bands of uncertainty. Summary In this chapter, four most common uncertainty analysis techniques were compared using an example of 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. Fuzzy set theory was used to propagate uncertainty associated with the velocity and dispersion coefficient. The comparison of the results of 1-D Monte Carlo simulation and those of the fuzzy set theory showed that the later was more conservative. The application of a 2-D 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 also indicated a good agreement between P-Box and fuzzy set results. This chapter aimed at selecting an appropriate technique to model uncertainty in contaminant intrusion events. The results of this chapter showed the abilities of the introduced techniques to model uncertainty related to the parameter and data. However, further investigations showed that in contaminant intrusion events, in addition to the parameter and data uncertainties, structural uncertainty exists. Structural uncertainty refers to a situation that we are uncertain about the form of the model in-use. For example, to determine the potential of contaminant around a WDS, contaminant propagation model needs to be selected. However, choosing the right model is not often a straightforward task because it depends on the conditions of the soil and groundwater level (saturated/unsaturated). In other words to model uncertainty in contaminant intrusion events the selected model should be able to estimate structural uncertainties as 36 well as data and parameter uncertainties. Fuzzy rule based model can provide a facility to model structure, data and parameters uncertainties. In addition to its flexibility and simplicity, using the fuzzy rules, the rigorous mathematical complexities can be circumvented. A fuzzy rule based model can handle complex and vague situations, and can describe nonlinear relationships through simple rules. 37 Chapter 4: Potential of Contaminant Intrusion in Water Distribution Systems A version of this chapter has been submitted for publication in Urban Water Journal. Mansour-Rezaei, S., Naser, Gh., and Sadiq, R. (2011). "Potential of Contaminant Intrusion in Water Distribution Systems". Preface Contaminant intrusion in a WDS has three main prerequisites including the availability of (a) a contaminant source in the proximity of the system, (b) a driving force (low/negative pressure), and (c) a pathway. A driving force (e.g., low/negative pressure) is induced by a pump operation, a pipe replacement, a valve closure/opening, or demand variations. A contaminant may intrude a WDS through a variety of pathways including submerged air valves, leak points, repair and installations, faulty seals, joints, and service connections. Pipe corrosion can also result in mechanical failure or pipe break. Contaminant intrusion is a complex phenomenon for which the data are generally uncertain. It is very challenging (if not impossible) and often expensive to reliably and accurately provide all the required data for a real-life WDS. Uncertainties related to contaminant intrusion events can arise from data, structural, and parameter uncertainties. For example, to accurately estimate the contaminant potential around a WDS a model should be employed for the conditions of the soil around the system. However, choosing the right model is not a straightforward task. Even if an appropriate model is selected for groundwater flow and soil condition, there are still questions to be answered with regards to the value of the model parameters. Fuzzy logic considers uncertainties related to a model structure and/or parameters. In addition to its flexibility and simplicity, using the fuzzy rules, the rigorous mathematical modeling can be circumvented. Studying a WDS with cast-iron pipes, this chapter developed a model to determine the potential of contaminant intrusion at each junction of the system. To estimate the potentials of existence of contaminant sources, a fuzzy-rule model was created. To estimate the potentials of driving force, a transient hydraulic model was proposed. To estimate the potential of existence of pathway, a pit corrosion model for cast iron pipes 38 was developed. A fuzzy-rule model was then employed to integrate the potentials of existence of contaminant sources, driving forces, and pathways existence and to determine the potential of a WDS to contaminant intrusion. 4.1 Uncertainty Analysis for Contaminant Intrusion Over the last four decades, the practical results of fuzzy system for management of uncertainties have led to its general acceptance in various engineering disciplines (Sadiq et al., 2009). The required data to model contaminant intrusion in a WDS are generally incomplete, non-specific, and uncertain (Deng et al., 2010). Thus, to make inference about the potential of contaminant intrusion, this research employed a fuzzy-rule-based model. A fuzzy-rule-based model can handle complex and vague situations, and can describe nonlinear relationships through simple rules. Since a fuzzy-rule-based model is based on qualitative descriptions used in natural language, it is easy to use (Jang et. al., 1998). In a fuzzy rule-based model, fuzzy if-then rules are used to represent the relationships between variables (Sadiq et al., 2009): if antecedent proposition then consequent proposition. The antecedent proposition defined by a fuzzy proposition of the form ‘X is A’ where X and A are a linguistic variable and a constant, respectively. Depends on the degree of similarity between X and A, a real number between zero and one is defined as a proposition’s membership value (Mamdani, 1977). As example of if-then rules is: Ri: If X is Ai then Y is Bj i = 1, 2, … , L; j = 1, 2, …, N 4.1 where Ri is the rule number i, X is the antecedent (input) linguistic (fuzzy) variable and Ai is a fuzzy subset corresponding to an antecedent linguistic constant. Similarly, Y is the consequent (output) linguistic (fuzzy) variable and Bj is a fuzzy subset corresponding to a consequent linguistic constant. The full relationship between X and Y according to rule i can be computed by Mamdani’s implications, which is the most commonly used implications in engineering applications of fuzzy rule-based model (Sadiq et al., 2009). In Mamdani’s method, conjunction A ∧ B is computed by a minimum (and type t-norm or conjunctive) operator. For a two-input single-output fuzzy rule-based model, Mamdani’s inference algorithm is described as: 39 i = 1, 2, …, L 4.2 Ri,j: If X1 is Ai and X2 is Cj then y is Bk; j = 1, 2, …, M k = 1, 2, …, N In this model, the antecedent proposition is obtained as the Cartesian product of fuzzy sets A and C. Hence, the degree of fulfillment is given by: { } βi , j = max µ A ( x1 ) ∧ µ A ' ( x1 ) ∧ max µC ( x2 ) ∧ µC ' ( x2 ) X X 1 i j 2 4.3 where A’ and C’ are input fuzzy numbers (or singletons), which are mapped on set A and C, respectively. The output fuzzy set of the linguistic model becomes: µ B ( y) = max βi, j ∧ µ B ( y) i =1,2,..., L k 4.4 j =1,2,.., M Please note that the Mamdani’s algorithm can also be extended to multiple-outputs model, which is a set of multi-inputs single-output models. A AND B C Fig. 4.1- Two causal concepts connecting to an effect concept (Sadiq et al., 2009). Consider an effect node B, which is connected by two causal concepts A and C. Fig. 4.1 shows a graphical representation of the causal relationships. The AND action indicates that both A and C are simultaneously required for B to occur. Fig. 4.2 provides additional details of this two-input and one-output model. The process is shown in three distinct steps, namely, fuzzification, inference (a rule base and an inference engine) and defuzzification. Assume that causal concepts A and C are activated at levels of A’ = 0.4 and C’ = 0.6, respectively. The rule set consists of 6 rules (3 × 2) and input activation signals (A’) and (C’) fire the first 4 rules to determine output B’, which is defined over the universe of discourse Y. The defuzzification step provides a discrete (crisp) value of 40 an output of B (i.e., B’). The crisp value approximates the deterministic characteristics of the fuzzy reasoning process based on the output fuzzy set µBk(y). Converting the uncertainty to the deterministic values helps to take an applicable action when solving real-world problems. This research employed centeroid method for defuzzification step. The centeroid method calculates the mass center of the results to provide a crisp value. 41 Step 1: Fuzzification A’ = 0.4 A2 A1 Step 2: Inference C’ = 0.6 A3 B1 C2 C1 B2 B3 β2,1 β2,2 B4 0.7 0.6 0.4 y 0.3 β1,1 β1,2 Y If and x1 x2 Then: Y Rule-base (Ri,j): Antecedent 1 = Ai; i = 1,…, L; L = 3 Antecedent 2 = Cj; j = 1,…, M; M = 2 Consequent = Bk; k = 1,…, N; N = 4 Total number of rules = 3 × 2 = 6 Rule 1 Rule 2 Rule 3 Rule 4 Rule 5 Rule 6 R1,1 : If R1,2 : If R2,1 : If R2,2 : If R3,1 : If R3,2 : If B’ Step 3: Defuzzification A1 and C1 then B1 A1 and C2 then B2 A2 and C1 then B2 A2 and C2 then B3 A3 and C1 then B3 A3 and C2 then B4 Defuzzified value yo yo Defuz. (B’) • x1 has support (membership >0) in A1 and A2, and x2 has support in C1 and C2, consequently only the first 4 rules are “fired” (or activated). • Rule 1: µA1(x1) = 0.3 , : µC1(x1) = 0.4 ⇒ β1,1 = 0.3 ∧ 0.4 = 0.3 ⇒ µB1(y) = 0.3 • Rule 2: µA1(x1) = 0.3 , : µC2(x1) = 0.6 ⇒ β1,2 = 0.3 ∧ 0.6 = 0.3 ⇒ µB2(y) = 0.3 • Rule 3: µA2(x1) = 0.7 , : µC1(x1) = 0.4 ⇒ β2,1 = 0.7 ∧ 0.4 = 0.4 ⇒ µB2(y) = 0.4 • Rule 4: µA2(x1) = 0.7 , : µC2(x1) = 0.6 • For every subset Bk choose maximum membership, i.e., µB2(y) = 0.4 • Defuzzification using centroid method, i.e., Defz. (B’)=0.4 ⇒ β2,2 = 0.7 ∧ 0.6 = 0.6 ⇒ µB3(y) = 0.6 Fig. 4.2- Fuzzy rule-based model: making inference using two ‘causal factors’ (Sadiq et al., 2009). 42 4.2 Contaminant Intrusion Potential 4.2.1 Potential for Contaminant Source The first prerequisite for contaminant intrusion is the availability of a contaminant source. In a field study in the province of Quebec (Canada), Ebacher et al. (2011) indicated that some parts of the WDS were below the groundwater table. They reported that the groundwater depth varies between 0.49 m and 3.18 m below the ground surface, while the minimum pipe burial depth based on provincial requirements was 2.2 m in the study area. This research assumed two scenarios for contaminant around a WDS. Scenario 1 assumed the groundwater table to be above a pipe centerline and thus a contaminant reaches the WDS by groundwater flow. Scenario 2 examined a WDS located in an area contaminated by a pollutant leaking from a source (e.g., sewage lines or so). For a case of the first scenario, this research assumed potential of a contaminant source (S) at a junction of a WDS is affected by two causal concepts: 1) distance between the contaminant source and the junction (R) and 2) the angle between groundwater flow direction and the straight line connecting the source to the junction (θ). MATLAB fuzzy logic tool-box (Jang et al., 1998) was employed to estimate the potential of a contaminant source at each junction of a WDS. A four-level linguistic scale was assumed for R and θ that includes low (L), medium (M), high (H), and very high (VH). Similarly, a three-level scale was assumed for S including low (L), medium (M), and High (H). The rule set consists of 4×4=16 rules. Table 4.1 indicates the break points of membership functions. The trapezoidal and triangular membership functions have four and three break points, respectively. Table 4.2 shows the potential of contaminant rule set. The full relationship between R, θ and S according to rule i was computed using Mamdani’s inference method (Mamdani, 1977). For defuzzification, the centeroid method was used, which calculated the center of the area under the curve. The crisp value approximated the deterministic characteristics of the fuzzy reasoning process based on the output fuzzy set µC. Fig. 4.3 illustrates the Peseudocolor plot of the proposed fuzzy rule-based model for the potential of contaminant. The horizontal axis shows R and vertical axis indicates θ; the brighter color refers to the higher potential of the contaminant source (S). 43 For a case of the second scenario, contaminant seeps from a pollutant source and develops a contaminated zone underneath. In this case, the research assumed the potential of contaminant source around the WDS to be proportional to the distance between the contaminant source and the WDS. Vairavamoorthy et al. (2007a and 2007b) conducted analytical studies of seepage from the pollutant sources in saturated and unsaturated soil (e.g., ditch, open drains, and sewer pipes). They showed that the width of contaminated zone under a leaking pollutant source approaches to B+2H where B is the width of contaminant source and H is contaminated liquid depth in the source (Fig. 4.4). Thus, for the parts of the WDS located in vicinity of pollutant source (but not exposed to groundwater flow), the potential of a contaminant source (S) was estimated by: for X ≤ ( B + 2 H ) / 2 1 S = ( B + 2 H ) / 2 X for X > ( B + 2 H ) / 2 4.5 where X is horizontal distance between the center of contaminant source and the junction. 4.2.2 Potential of Deriving Force The second prerequisite for contaminant intrusion is a driving force. To detect low/negative pressure induced under rapidly changing conditions, transient hydraulic analysis is required (Boulos et al., 2005 and Fernandes and Karney, 2004). EPANET (Rossman, 2000) and other existing hydraulic models ignore inertia and compressibility effects and assume either steady or quasi-steady conditions for flow. Fernandes and Karney (2004) showed that steady and quasi-steady assumptions for a real-life WDS must be reconsidered. Conducting 1-D and 2-D transient modeling, Naser and Karney (2007 and 2008) showed that the advantages in flow and water quality representations in a 2-D model were drastically compromised due to computational requirements and simulation time. This makes a 2-D simulation computationally-impractical for analyzing a real-life WDS. Thus, this research developed a 1-D transient hydraulic model to detect locations in a WDS with low/negative pressure and to estimate the magnitude of pressure. 44 Table 4.1- Fuzzy linguistic sets for causal and effect concepts L M R Trapezoidal MF (0, 0, 10, 20) Θ Trapezoidal MF (0, 0, 15, 30) S Triangular MF (0, 0, 0.4) (10, 20, 50, 70) (15, 30, 45, 60) (0.1, 0.5, 0.9) H (50, 70, 100, 120) (45, 60, 75, 90) (0.6, 1, 1) VH (100, 120, 250, 250) (75, 90,180,180) Fig. 4.3- Peseudocolor plot of the proposed fuzzy rule-based model. Fig. 4.4- Contamination zone underneath a contaminant source (Vairavamoorthy et al., 2007a). 45 Table 4.2- Rule set for potential of contaminant source. Rule R Outcome 1 L AND L THEN H 2 L AND M THEN H 3 L AND H THEN H 4 L AND VH THEN L 5 M AND L THEN H 6 M AND M THEN M 7 M AND H THEN M 8 M AND VH THEN L 9 H AND L THEN M 10 H AND M THEN L 11 H AND H THEN L 12 H AND VH THEN L 13 VH AND L THEN L 14 VH AND M THEN L 15 VH AND H THEN L 16 VH AND VH THEN L θ A transient flow is modeled by continuity and momentum equations (Boulos et al., 2005 and Larock et al., 2001): ∂H a 2 ∂V + =0 ∂t g ∂x 4.6 ∂V ∂H fV V +g + =0 ∂t ∂x 2d 4.7 where V is flow velocity (m/s), H is piezometric head (m), d is internal pipe diameter (m), g is acceleration (m/s2) due to gravity, f is Darcy-Weisbach friction factor, a is acoustic wave speed (m/s), x is the distance (m) along the pipe centerline, and t is time (second). To solve the governing equations both Method of Characteristics (MOC) and Wave Characteristic Method (WCM) can be used. The former is an example of Eulerian method, and the later is an example of Lagrangian method. Wood et al. (2005) showed that both approaches always lead to identical results when same data and model are used. This research used a MOC-based transient hydraulic model to detect low/negative 46 pressure in a WDS. Larock et al. (2000) and Boulos et al. (2005) provided further details about MOC. 4.2.3 Potential of Pathway Intrusion can occur through a several pathways including submerged air valves, leak points, repair and installations, faulty seals and joints, service connections, and so on. Given that 50% of the water pipes in Canada are old corroded gray cast-iron; this research estimated the potential of a pipe break/failure due to its internal/external corrosion. Assuming a high corrosion rate during early pipe ages followed by a stabilizing rate at a certain value during the self-inhibiting process, Rajani and Markar (2000) proposed a two-phase model to study the external corrosion of a pipe-wall: ( d ext = aT + b 1 − e − cT ) 4.8 where, dext is external corrosion depth (mm), a is final pitting rate constant (mm/year), b is pitting depth scaling constant (mm), c is corrosion rate inhibition factor (year-1), and T is exposure time (year). To determine internal corrosion depth, Yamini and Lence (2010) employed a linear corrosion model of: d int = C r T 4.9 where, T is exposure time (year) representing the pipe age and Cr is internal corrosion rate (mm/year). Yamini and Lence (2010) related the later to the chlorine decay rate by using the 0th-order model of Kiene et al. (1998): C C r = C 0 1 − t C0 DV K L 4.10 where, V is flow velocity (m/s), L is effective pipe length (m), Ct is chlorine concentration at time t, D is pipe diameter (m), K is theoretical decay constant, and C0 is initial chlorine concentration (mg/L). It should be indicated that Yamini and Lence (2010) assumed a uniform internal corrosion along the entire length of a pipe with constant rates; assumptions that are far from reality. However, it is adopted in this research due to its simplicity and ease of application. Moreover, this model is the only existing model in the literature for internal corrosion of a pipe. This research employed the ratio of the total corrosion depth (dext+dint) over the initial pipe thickness (t) at each 47 junction to estimate the potential of a pipe failure/pathway. Please note that dext and dint vary with time and so does (dext+dint)/t. 4.2.4 Intrusion Potential The potential of contaminant intrusion at each junction of a WDS was determined by applying a three-input one-output fuzzy rule-based model sketched in Fig. 4.5. The three inputs were “potential of contaminant”, “potential of pathway”, and “potential of low/negative pressure”, which were defined in Sections 4.2.1, 4.2.2, and 4.2.3. The output of the fuzzy rule-based model was “the potential of contaminant intrusion”. The domain of each input was divided into three sub-domains (i.e. “low”, “medium” and “high”). Thus, a total of 33 = 27 rules formed (Table 4.3). Fig. 4.6 and Fig. 4.7 shows the membership functions of inputs and output, respectively. The universe of discourse for “potential of contaminant”, “potential of pathway”, and “potential of low/negative pressure” were described over the continuous intervals [0, 1], [0, 1], and [-10, 50], respectively. Note that the first and second input parameters are dimensionless, while the unit of the third input parameter is metric. It should be emphasized that the universe of discourse was the scale on which each input is measured/observed. The values of three input parameters were mapped over their universe of discourse. The membership functions were evaluated with respect to the linguistic constants low, medium and high (L, M, and H). Then, the proposed rule set was used to make an inference at each junction. Section 34.1 described the inference mechanism. Fig. 4.5- Fuzzy rule-based model for potential of contaminant intrusion. 48 Table 4.3- Potential of contaminant intrusion rule set. Rule Potential of Contaminant 1 L Potential of Low/Negative Pressure AND L Potential of Pathway AND Outcome L THEN QL QL L 2 3 L L AND AND M H AND AND L L THEN THEN 4 5 L L AND AND L M AND AND M M THEN THEN L L 6 7 L L AND AND H L AND AND M H THEN THEN M M 8 9 L L AND AND M H AND AND H H THEN THEN H H 10 11 M M AND AND L M AND AND L L THEN THEN QL L 12 13 M M AND AND H L AND AND L M THEN THEN M L 14 15 M M AND AND M H AND AND M M THEN THEN M H 16 17 M M AND AND L M AND AND H H THEN THEN M H 18 19 M H AND AND H L AND AND H L THEN THEN QH L 20 21 H H AND AND M H AND AND L L THEN THEN L M 22 23 H H AND AND L M AND AND M M THEN THEN M H 24 25 H H AND AND H L AND AND M H THEN THEN H H 26 27 H H THEN AND M AND H THEN AND H AND H QL-quiet low; L-low; M-medium; H-high; QH-quiet high QH QH 49 Fig. 4.6- Membership functions of the inputs. Fig. 4.7- Membership functions of the output. 4.3 Case Study Adopted from literature (Boulos et al., 2005; Wood et al., 2005; and Jung et al., 2009), Fig. 4.8 depicts a schematic plan-view of a WDS considered in this study. On the figure, the number in a circle indicates the junction, while that in a square represents the pipe. The arrows in the figure show the direction of initial steady-state flow in a pipe. The system consists of nine cast-iron pipes, five junctions, a valve, and a reservoir with constant head of 80 m. Table 4.4 shows further information for each pipe including diameter, length, friction factor, wave speed, thickness, and age. As Table 4.4 shows the system has large-diameter pipes. This could be a skeletonized network showing only the big trunk mains of the system. The junctions were located at the same elevations. The demands at junctions 1, 2, 3, 4, and 5 were 0.6, 0.5, 0.5, 0.6 and 0.5 cms, respectively. The demand at the valve located at the downstream end of the pipe 7 was 0.5 cms. The system operated for 50 seconds steady-state condition. The steady-state discharges in the pipes 1, 2, 3, 4, 5, 6, 7, 8, and 9, were 3.2, 1.42, 1.28, 0.52, 0.41, 0.84, 0.5, 0.26, and 0.25 cms, respectively, while the pressure head at junctions 1, 2, 3, 4, and 5 were 41.4, 62.5, 49.7, 38.1, and 38.8 m, respectively. A transient scenario was generated through a 50 valve closure. The valve was closed linearly from the fully open position to completely shut within one second. The dimensionless head-loss coefficient of the valve (KL) for 100, 80, 60, 40, and 20 percent opening was assumed 0.19, 1.15, 5.6, and 31.9, respectively. Fig. 4.8 indicates two zones assumed around the WDS. The groundwater depth in Zone 1 is lower than the pipe burial depth, while it is higher than the pipe burial depth in Zone 2. As Fig. 4.8 shows, the origin of the coordinate system located at node 1. Groundwater in Zone 1 is contaminated with a pollutant (S1) leaking at a constant rate at a point at X = -47 m and Y = 17 m Zone 1 (Fig. 4.8). The direction of regional groundwater flow with X-direction was also assumed to be 20º. Another leaking contaminant source (S2) with B=5 m and H=2 m located at a point at X = 337 m and Y = 332 m in Zone 2 (Fig. 4.8). Table 4.4- Characteristics of the pipes in the WDS case study. Pipe Number 1 2 3 4 5 6 7 8 9 Diameter (mm) Length (m) Friction Factor (f) Wave Speed (m/s) Thickness (mm) Age (year) 900 610 0.02 1000 35 10 750 914 0.02 1000 32 15 600 610 0.02 1000 28 20 450 457 0.02 1000 25 20 450 549 0.02 1000 25 40 750 671 0.02 1000 32 15 900 610 0.02 1000 35 50 600 457 0.02 1000 28 40 450 488 0.02 1000 25 40 4.4 Results and Discussion Applying the three-input one-output fuzzy rule-based model (Section 4.2), this section makes inference about the potential of contaminant intrusion by bringing together the potentials of the three essential elements of intrusion (i.e., contaminant source, driving force, and pathway). Potential of Contaminant Source: Using R and θ for junctions located in the zone 1 and employing fuzzy rule-based model (Section 4.2.1), the potential of contaminant source at junctions 1, and 4 were estimated at 0.87, and 0.14, respectively. As well, using the ratio of (B+2H)/2X for junctions located in the zone 2 (Section 4.2.1), the potential of contaminant source at junctions 2, 3, and 5 were estimated at 0, 0.23, and 0.01, respectively. The higher potential indicates the higher possibility for contaminant existence at a junction. Thus, junctions 1 and 3 had the highest potential of contaminant. 51 1 S2 4 3 2 5 2 5 8 Zone 2 Y 3 9 Zone 1 S1 X 1 4 6 7 Regional groundwater flow direction Fig. 4.8- Case study: WDS layout. Potential of Driving Force: The variations of pressure head during transient events were simulated by employing the transient hydraulic model proposed in Section 4.2.2. Fig. 4.9 indicates the results of pressure head variations in the system during valve closure. The results indicated that the lowest pressure during transient event at junctions 1, 2, 3, 4, and 5 were -3, 23.8, -6.3, -6.8, and 1.1 m, respectively. Thus, junctions 1, 3, and 4 experienced negative pressures. The pressure at junction 1 dropped from 41.4 to 3 m at time 58.9, for junction 3 dropped from 49.7 m to -6.3 m at time 59.8 s, and for junction 4 the drop was from 38.1 m to -6.8 m at time 58.6 s. The system reached the new steady-state hydraulic conditions after about 150 seconds. Note that the new steady state discharges after the transient event in pipes 1, 2, 3, 4, 5, 6, 7, 8, and 9 were 2.07, 1.182, 1.018, 0.381, 0.301, 0.542, 0.00, 0.059 and 0.178 cms, respectively. 52 Fig. 4.9- Variations of head with time for junctions 1, 2, 3, 4, and 5. Potential of Pathway: To define the potential of a pathway, the ratio of (dext+dint)/t proposed in Section 4.2.3 was considered. The parameters a, b, and c in the external corrosion model were assumed to be 0.009 mm/year, 6.27 mm, and 0.1 year-1, respectively (Sadiq et al., 2004), while C0, K, Ct/C0, and L in internal corrosion model were set at 1 mg/L, 1.01x10-6, 0.94, and 6 m, respectively (Yamini and Lence, 2010). Using the age and pipe wall thickness given in Table 4.4, the potential of pathway was calculated for each junction. Fig. 4.10 illustrates the variations of (dext+dint)/t with pipe age for pipes 5 and 9 as an example of results. The results show that pipe 5 corrodes faster than pipe 9. As Error! Reference source not found. suggests the external corrosion depths (dext) for both pipes were identical. However, the internal corrosion rate (Cr) for pipe 5 was higher than that for pipe 9 because of the higher flow velocity and larger pipe diameter (see 53 ). It should be noted that the wall thicknesses for both pipes were identical (25 mm). Considering the pipes ages as indicated in Table 4.4, (dext+dint)/t for junctions 1, 2, 3, 4, and 5 were obtained as 0.8, 0.5, 0.5, 0.4, and 0.55, respectively. The potential of pathway for each junction was assumed to be the maximum of the potential of pathway of the pipes joining to the junctions. For example, the potential of pathway for pipes 3 and 5 were 0.75 and 0.8, respectively. Thus, the potential of pathway for junction 1 was 0.8. Fig. 4.10-Variations of (dext+dint)/t with pipe age for pipe 5 (top) and pipe 9 (bottom). Potential of Intrusion: To define the potential of contaminant intrusion at each junction of the WDS, the three-input one-output fuzzy rule-based model proposed in Section 4.2.4 was used. Table 4.5 indicates the potential of contaminant intrusion at each junction of the WDS computed by fuzzy rule-based model. The table also shows the potential of contaminant, low/negative pressure, and the pathway for each junction of the WDS. The results showed that junction 1 has the highest potential for contaminant 54 intrusion. Please note that the potential of contaminant intrusion shows the possibility of an intrusion event, which is different from the probability. Table 4.5- Potential of contaminant intrusion in the WDS Junction 1 2 3 4 5 Potential of Contaminant 0.87 0.00 0.23 0.14 0.01 Potential of Low/Negative Pressure (m) -3.0 23.8 -6.3 -6.8 1.1 Potential of Pathway 0.8 0.5 0.5 0.4 0.55 Potential of Contaminant Intrusion 0.92 0.25 0.67 0.52 0.52 Summary In this chapter, a fuzzy rule-based model was employed to define the susceptible locations of a WDS to contaminant intrusion. The potential of a contaminant source around the system, the potential of driving force, and the potential of pathway due to pipe corrosion were integrated to make inference on possibility of the contaminant intrusion to the system. The results showed that fuzzy rule-based model can prioritize junctions of a real-life WDS based on the potential risk of contaminant intrusion using very vague and uncertain data. Such information will help engineers and water authorities to study and understand hazard and vulnerability of a WDS to contaminant intrusion. Such understanding can eventually help them to make more informed decisions on selecting the most efficient operating strategy to mitigate the risk. To study contaminant intrusion and its propagation along a WDS, the proposed model can be coupled with an ingress model (for estimating rate of contaminant entering a system during an intrusion event) and a water quality model (for tracking the fate of a contaminant along a system). Such a comprehensive study will provide a valuable tool for water authorities to eventually conduct an exposure and risk assessment if an intrusion occurs. 55 Chapter 5: Ingress Model A version of this chapter has been accepted for publication in Journal American Water Work Association (AWWA). Mansour-Rezaei, S., Naser, Gh. (2013). "Contaminant Intrusion in Water Distribution Systems: An Ingress Model". Preface In a water distribution system, a low/negative pressure event may occur under normal operation of the system that may lead to contaminant intrusion. To assess the risk associated with contaminant intrusion, exposure assessment is required. While a model was proposed in Chapter 4 to identify susceptible locations of WDS to contaminant intrusion, in this chapter a model is proposed to estimate intrusion flow rate. Although, due to its simplicity, an orifice equation is often suggested for estimating the intrusion rate into a WDS through a leak-point, it may not be necessarily accurate for different shapes and sizes of orifices. The orifice equation; Q = CDAo(∆P)0.5, assumed the intrusion rate to be directly proportional to the orifice area (Ao) and square root of pressure difference across the orifice (∆P). In the orifice equation, the effects of the surrounding soil and type of flow through porous media (laminar and/or turbulent) are not taken into account. Indeed, the soil hydraulic conductivity may have significant impacts on the intrusion rate. To have more accurate intrusion exposure assessment, the impacts of surrounding soil on intrusion rate must be considered. In this chapter a model is proposed to make a more realistic estimation of intrusion rate. A two-level factorial design is used to screen the effect of pipe internal pressure, soil porosity, average diameter of soil particles, size of the leak-point, and pipe burial depth. The Bukingham π-theorem is applied to group the influential factors into a set of dimensionless π-parameters. The relationship among the π-parameters is obtained by using the numerical data generated by a twodimensional finite element model for flow in a saturated zone surrounding a pipe. A regression analysis is then applied to find the best fit for the data. 56 5.1 Modelling Approach 5.1.1 Model Assumptions Intrusion into a WDS is a complex phenomenon; saturation of surrounding soil can change temporally and spatially. The induced low/negative pressure has a transient nature; the low/negative pressure also has an interaction with a leak-point. Thus, to simulate intrusion rate into a WDS, this research employed three simplifications. In reality the saturation of surrounding soil depends on groundwater level, leakage flow rate (before intrusion events), and physical characteristics of the soil. Given that water leaks outside a pipe for a while through the leak-points under positive internal pressure conditions prior to the low/negative pressure incident, the soil around a leak-point may have high saturation. Thus, this research assumed the surrounding soil fully saturated. Although this simplifying assumption may not accurately represent reality, it leads to more conservative results (i.e. greatest intrusion). As the low/negative pressure inside a pipe has impacts on a relatively small region of the soil around a leak-point, this research assumed the soil around the leak-point to be homogeneous and isotropic. It should be mentioned that this research neglected the effect of washout and void formation surrounding the leak point. Besner et al. (2011) reported two types of low/negative pressures that may be observed in a WDS. These include transient and sustained low/negative pressures. A transient low/negative pressure induced by rapid changes in flow velocity may last from milliseconds to minutes. A sustained low/negative pressure induced by operation/maintenance activities lasts in the range of minutes to hours. Although low/negative pressure may vary rapidly with time (transient low/negative pressure) (Collins et al., 2011; Ebacher et al., 2011 and Besner et al., 2011), this research modeled an intrusion event occurred under a steady-state low/negative pressure scenario. The proposed model in this research, however, can also be employed to estimate intrusion rate under a transient low/negative pressure event. This is explored later. Although Collins et al. (2011) argued a decrease in the magnitude of transient low/negative pressure by an increase in the intrusion rate; in this work, such effects were ignored for the purpose of simplicity. Moreover, the research ignored the role of a leak-point as a pressure relief in WDS during to the transient. 57 5.1.2 Governing Equations Momentum Equation. To describe fluid flow through a porous media, the Darcy equation is often suggested (Huang and Ayoub, 2008; Collins et al., 2010): −∇p = µ k V 5.1 where ∇p is pressure gradient vector, µ is dynamic viscosity (ML-1T-1) of the fluid, k is permeability (L2) of porous media, and V=q/A is superficial velocity (LT-1), q is volumetric flow rate (L3T-1), and A is surface area (L2) of the porous media perpendicular to the flow direction. The Darcy equation assumes a linear relation between pressure loss and fluid flow (Huang and Ayoub, 2008 and Collins et al., 2010). Fluid flow through a porous media depends on both permeability and the pressure gradient; however, due to high pressure gradients, non laminar flow is expected in most surrounding soil during leakage/intrusion (Van Zyl and Clayton, 2007 and Collins et al., 2010). The Darcy equation is valid for laminar flows when the Reynolds number is smaller than unity. However, for situations including higher flow velocities, flow is generally believed to follow the Forchheimer equation (Macdonald et al., 1979; Vafai, 2005; Huang and Ayoub, 2008; and Collins et al., 2010). Thus, this research employed the Forchheimer equation (Collins et al., 2010) to simulate fluid flow through the porous media during intrusion events: µ −∇p = V + b ρ V V k 5.2 where ρ is density (ML-3) of the fluid, and b is Forchheimer constant (i.e., a properties of the porous media). The second term in the right side of the Forchheimer equation is to consider pressure drop due to the nonlinear effects in high-velocity flows. Using Ergun’s empirical expression, permeability can be related to physical characteristics of soil (Macdonald et al., 1979 and Vafai, 2005) as: K = ε 3 Ds2 /150(1 − ε ) 2 5.3 where Ds is the average diameter (L) of soil particles and ε is the soil porosity. The Forchheimer constant b can be determined by (Macdonald et al., 1979 and Vafai, 2005): b = 1.75 (1 − ε ) 1 ε 3 Ds2 5.4 58 It should also be noted that Forchheimer equation approaches to Darcy equation when the flow is laminar (Faybishenko, 2005). Continuity Equation. Continuity equation for an incompressible fluid flow in a porous media is: ∇.V = 0 5.5 The pore pressure p and velocity field V were determined by numerical integration of the momentum and continuity Equations (5.2 and 5.5) along with appropriate boundary conditions for a saturated zone surrounding a leak-point. Boundary Conditions. Considering a Cartesian computational domain around a leakpipe (Fig. 5.1), this research assumed a no-flow boundary condition for the nodes located at the left, right, and bottom of the domain. However, for the nodes at the top (on the ground-level), the pore pressure was set at zero (to make sure the soil always stay saturated). It should be noted that hydrostatic-pressure condition could also be applied for the left, right, and bottom boundaries. However, no flow conditions at these boundaries lead to hydrostatic pressure at the corresponding nodes if the computational domain is larger than the leak point. This will be discussed with further details in the “Results and Discussion” section. For the nodes located on the leak-point, the pore pressure was set as the pipes’ internal pressure. Notice that due to the gravity effect, locating the leak point at the top of the pipe leads to more intrusion flow rate in compare to the case the leak point located at the bottom of that. Intrusion is, in its nature, a three-dimensional (3D) event. In fact, the 3D analyses of Collins et al. (2010) showed that a contaminant can even migrate from locations below a pipe and intrude into the pipe. The three dimensional (3D) intrusion phenomena can be studied by a two-dimensional (2D) model in two ways including 2D in plan- or surface-view and 2D in section- or depth-view. In fact, a 2D model (depth-view) may lead to more conservative results (i.e. the higher intrusion rate), while it is simpler. Thus, this research employed a 2D model for simplicity purpose. Thus, intrusion rate was determined by the average velocity perpendicular to the leak-point multiplied by leak-point cross sectional area. To simulate fluid flow through a porous media, this research applied the 2D finite element model established in the software Abaqus (Abaqus Analysis User's Manual, 2010). For numerical simulation, quad-shape finite elements were used, while their displacements were restrained along the coordinate 59 axes. The preliminary simulation results showed that the pressure gradient is very high only at the vicinity of leak-point. Thus, to improve accuracy and computational efficiency (e.g., smaller simulation time), non-uniform elements were applied. Moreover, finer elements were chosen at the vicinity of the leak-point, while coarser elements were employed at further distances from the pipe (Fig. 5.1). Although nonuniform elements may lead to the distortion of some elements, the mesh verification showed that no error occurred due to non-uniformity. Abaqus Analysis User's Manual (2010) provides additional details about mesh generation/verification. Fig.5.1-A two-dimensional (2D) finite element model to simulate fluid flow through porous media. 5.2 Effective Factors In a complex system, conducting comprehensive experiments to consider all the effective parameters is quite challenging (if possible at all). A practical and efficient approach to collect required information is to generate experiments by computer simulations (Fang, 2006 and Motamedi and Milani, 2010). Three components of a typical experiment are factors (inputs), responses (outputs), and levels (settings of 60 factors). Modellers may screen factors effects (sensitivity analysis) for many reasons (Hamby, 1994) such as the need to detect factors with the highest contributions on the responses or to detect insignificant factors which can be eliminated from a model or experiment. Factorial analysis or design offers a simple but efficient way to estimate the main and interaction effects of factors on a system’s response (Montgomery, 2008). The main effect is the effect of one factor on the response ignoring the effects of all other factors. An interaction effect is the effect of combination of factors on the response. In full factorial design, all possible combinations of the levels of factors are investigated (Fang, 2006 and Montgomery, 2008). A two-level factorial design (e.g., 2k factorial design with k refers to the number of factors) is a special case of factorial design in which each factor is investigated only at two levels. A two-level factorial design assumes linear variations of response over the factors levels. This research applied Minitab software (Ryan and Joiner, 2001) to perform a two-level factorial design. Table 5.1- (-1) and (+1) levels of factors Levels Factors ε -1 1 P (kPa) Ds (m) Do (m) Dp (m) h (m) -90 0.25 0.001 0.001 0.2 1 0 0.4 0.01 0.01 0.6 3 Six factors (k = 6) that are expected to play key roles in intrusion events are a pipe internal-pressure (P), soil porosity (ε), average diameter of soil particles (Ds), leak-point diameter (Do) or cross sectional area (Ao), pipe diameter (Dp), and pipe external pressure or its burial depth (h). The latter is the distance between ground surface and the pipe centerline. This research assumed high and low levels for the influential factors listed in Table 5.1. The low and high values for ε, Ds, and Do are similar to those selected by Collins et al. (2010). The low value of P was assumed such that no cavitation occurs in the pipe system (-90 kPa). The low (high) values of Dp and h were also assumed to be 0.2 (0.6) m and 1 (3) m, respectively. For simplicity, high and low levels of factor were represented by 1 and -1, respectively. Sixty four (26=64) tests (i.e., computer simulations) were carried out to consider all the combinations of two levels of these six factors. Table 5.2 shows the experiments’ matrix with corresponding response 61 computed by the numerical model. To calculate intrusion rate or discharge (Q), the average velocity over the leak-point (V) was multiplied by corresponding leak-point cross sectional area (Ao). Table 5.2- Two level factorial design matrix with results Run No. P(kPa) ε Ds(m) Do(m) Dp(m) h(m) Q(m3/s) 1 1 -1 1 1 -1 1 6.90E-05 2 1 1 1 -1 1 1 4.70E-06 3 -1 -1 1 -1 -1 -1 3.54E-06 4 -1 1 -1 1 1 1 8.97E-05 5 -1 1 1 1 1 -1 2.78E-04 6 1 1 -1 1 -1 1 4.03E-05 7 1 -1 1 1 1 -1 3.15E-05 8 -1 1 -1 -1 1 -1 2.63E-06 9 -1 -1 -1 -1 1 1 1.20E-06 10 -1 -1 1 -1 1 1 4.32E-06 11 1 -1 1 1 1 1 6.51E-05 12 -1 1 -1 1 -1 1 9.24E-05 13 -1 -1 1 -1 -1 1 4.34E-06 14 -1 1 1 -1 -1 1 9.92E-06 15 -1 1 -1 1 1 -1 7.95E-05 16 1 1 -1 1 -1 -1 2.00E-05 17 -1 1 1 1 -1 1 3.26E-04 18 1 1 -1 -1 1 1 1.35E-06 19 -1 -1 1 1 -1 1 1.42E-04 20 1 -1 -1 1 -1 1 1.19E-05 21 1 -1 -1 -1 1 1 4.96E-07 22 -1 1 -1 1 -1 -1 8.44E-05 23 1 1 1 -1 -1 -1 4.86E-06 24 -1 1 1 -1 -1 -1 8.10E-06 25 1 1 -1 -1 -1 -1 6.61E-07 26 1 1 1 -1 1 -1 2.30E-06 27 -1 -1 -1 1 1 1 3.20E-05 28 1 -1 1 -1 1 1 2.05E-06 29 1 -1 1 -1 -1 1 2.12E-06 30 -1 1 1 1 -1 -1 2.94E-04 31 -1 -1 -1 1 -1 -1 3.08E-05 32 1 -1 -1 -1 -1 1 5.21E-07 Run No. P(kPa) ε Ds(m) Do(m) Dp(m) h(m) Q(m3/s) 33 1 1 -1 1 1 1 3.77E-05 34 -1 1 1 1 1 1 3.17E-04 35 1 1 1 1 -1 -1 8.73E-05 36 -1 -1 1 1 1 1 1.38E-04 37 1 -1 -1 1 -1 -1 5.09E-06 38 1 1 1 1 1 -1 7.33E-05 39 -1 1 -1 -1 -1 -1 2.45E-06 40 1 1 1 1 1 1 1.51E-04 41 1 -1 1 -1 -1 -1 1.05E-06 42 -1 -1 1 -1 1 -1 3.80E-06 43 1 -1 -1 1 1 -1 3.86E-06 44 -1 1 -1 -1 -1 1 3.00E-06 45 -1 1 1 -1 1 -1 8.69E-06 46 -1 -1 -1 -1 1 -1 1.07E-06 47 1 1 1 1 -1 1 1.59E-04 48 -1 1 -1 -1 1 1 2.99E-06 49 -1 -1 -1 1 -1 1 3.31E-05 50 -1 1 1 -1 1 1 9.88E-06 51 1 -1 -1 1 1 1 1.09E-05 52 1 1 -1 -1 -1 1 1.40E-06 53 1 -1 1 1 -1 -1 3.76E-05 54 -1 -1 -1 1 1 -1 2.89E-05 55 1 -1 -1 -1 1 -1 2.03E-07 56 -1 -1 1 1 1 -1 1.21E-04 57 1 -1 1 -1 1 -1 1.00E-06 58 -1 -1 -1 -1 -1 1 1.21E-06 59 1 1 -1 1 1 -1 1.60E-05 60 1 -1 -1 -1 -1 -1 2.22E-07 61 -1 -1 -1 -1 -1 -1 9.90E-07 62 1 1 -1 -1 1 -1 6.22E-07 63 -1 -1 1 1 -1 -1 1.28E-04 64 1 1 1 -1 -1 1 2.41E-06 Using factorial design, the main and interaction effects of factors were estimated by (Jafari and Babadagli, 2008): 1 −1 E f = Rave − Rave 5.6 1 −1 where Ef is the effect of a factor, Rave and Rave are the average responses at levels 1 and -1, respectively. For example, the effect of P on intrusion was calculated as: 62 Ef = (6.9 x 10-5+4.7x 10-6+...)/32-(3.54 x 10-6+8.97x 10-5+...) /32= 4.49x 10-5 5.7 Note that the numbers in the first parenthesis in Equation (5.7) are discharges when P was at level 1, while those in the second parenthesis are discharges when P was at level -1. Table 5.2 provides additional details. Note that in 32 experiments P was at level -1; while it was at level 1 in the other 32 experiments. Fig. 5.2- Normal probability plot of factors and interactions effects To determine interaction effects, all possible combinations of factors were assessed including combinations of two factors (e.g. P x ε), three factors (e.g. P x ε x Ds), four factors (e.g. P x ε x Ds x Do), five factors (e.g. P x ε x Ds x Do x Dp) and six factors (e.g. P x ε x Ds x Do x Dp x h). The levels of interactions were calculated using the levels of main factors. For example, if the level of P was 1 and the level of ε was -1, the level of P x ε was determined as 1 x -1 = -1. Fig. 5.2 illustrates the normal probability plot of the effects. In factorial design, normal probability plot is employed to test whether a data set is normally distributed or not. In a normal probability plot, the data set is graphed against theoretical normal distribution. The data forming a straight line follows normal distribution. However, those data point lying far from straight line show a departure from normal distribution. In Fig. 5.2, the factors or interactions of factors having a large effect lie far from the normal line, while those having negligible effects lie along the normal line (Fernandez et al., 2003 and Montgomery, 2008). Fig. 5.2 63 clearly highlights that P, ε, Ds, Do, h and some of their interactions (P x Do, P x Ds, and Do x Ds and so on) are located far from the normal line. Thus, they have significant effects on intrusion. However, the effect of Dp is negligible. Fig. 5.3 shows the main effects of factors. The figure illustrates that the effect of P is negative meaning that if P increases, intrusion will decrease. Fig. 5.3 also shows that increasing the level of Do from -1 to 1 leads to the largest variations in response. On the figure, the numbers on the top shows the values of the factors at levels -1 and 1. The slopes show changes in intrusion discharge to changes in factors. This indicates Do as the most significant factor on the response. As well, the results show that the effect of pipe diameter Dp is negligible because increasing the level of Dp from -1 to 1 does not have any significant impact on the response. Fig. 5.2 also shows that the effects of P, ε, Ds, and h are almost equally significant. Fig. 5.3- Main effects of factors Interaction among factors occurs when the effect of one factor on the response depends on the level of another factor (Fernandez et al., 2003 and Montgomery, 2008). Fig. 5.4 shows the interaction plot for the responses. Fig. 5.4-1 indicates the interaction between P and ε. In this figure, the red line shows the variation of intrusion with ε when P is at 1 level, while the black line shows the variation of intrusion with ε when P is at -1 level. Fig. 5.4-1 shows that intrusion increases when ε increases from level -1 (0.25) to 1 (0.5). However, the increase is higher when P is at level -1 in comparison with the case when P is at level 1. This observation can be interpreted as an interaction exists between ε and P. Similarly, Fig. 5.4-2, 5.4-3, 5.4-6, 5.4-7, and 5.4-10 show significant interactions among P, ε, Ds, and Do; facts that are also observed in Fig. 5.2. As Figs. 5.2 64 and 5.3 show the main effect of Dp is negligible. Figs. 5.2 and 5.4 also show the interactions of Dp with the other factors are negligible. Thus, Dp is removed from the set of significant contributing factor. 5.3 Dimensional Analysis Dimensional analysis is a useful approach for studying physical phenomena for which the relations among factors are not completely known (Misic et al., 2010). Dimensional analysis reduces the number of factors, and allows one to conduct fewer/selective experiments to investigate the relation among the contributing factors (Streeter and Wylie, 1998 and Misic et al., 2010). Dimensional analysis applies the Buckingham πtheorem to develop a set of dimensionless parameters among the contributing factors in a phenomenon (Streeter and Wylie, 1998 and Miragliotta, 2010). The theorem states that a physical phenomenon involving n factors can be rearranged into n - m independent dimensionless parameters for which m is the number of fundamental dimensions involved in the phenomenon (mass, length, time, etc.). Fig. 5.4- Interactions of P, ϵ, Ds, Do, Dp, h on intrusion. The straight black line indicates the effect of low values (-1) and that a red dashed line indicates the effect of high values (1). 65 The previous section concluded that internal pipe pressure (P), soil porosity (ε), average diameters of soil particles (Ds), leak-point diameter (Do) or cross sectional area (Ao), and pipe burial depth (h) were the most significant factors contributing to the intrusion rate (Q). Moreover, density and dynamic viscosity of contaminated water (ρ and µ) around a pipe as well as gravitational acceleration (g) are other essential factors for an intrusion event (Equation 5.1). Thus, a relation among factors can be expressed as: F1 (Q, ε , Ds , A0 , ∆P, ρ , g , µ ) = 0 5.8 where ∆P is the pressure difference between inside and outside a pipe equal to ρgh-P. Note that Equation (5.8) consists of eight factors with three repeating fundamental dimensions of time, length, and mass. Employing the Buckingham π-theorem, the following five dimensionless π-parameters can be identified: π1 = ρ Q / A0 Ds A ρ gDs Q , π 2 = 02 , π 3 = ,π 4 = 2 ,π 5 = ε µ Ds ∆P Ds ∆P / ρ 5.9 Note that π1 is the Reynolds number showing the ratio of inertia force to viscose force. The Reynolds number plays a key role in flow through a porous media (Khajeh and Maijer, 2010 and Khajeh and Maijer, 2011). Using dimensionless numbers, F1 was reformed as: F1 (π1 , π 2 , π 3 , π 4 , π 5 ) = 0 5.10 The functionality F1 can be rewritten as F2 for each constant value of π 5 as: π 4 = F2 (π1 , π 2 , π 3 ) for each π 5 5.11 To yield additional information about F2, experimental data is required. This research, however, employed the numerical model to generate a set of data and to ultimately identify the functionality F2. The next section discusses the test cases followed by the results. 5.4 Test Cases Numerical simulations were conducted for three different values of soil porosity ε (0.25, 0.3, and 0.4). For each ε, twenty four computer-tests were performed (Table 5.3), which were a full combinations of different values of internal pressure P (-90 and 0 kPa), soil particles size Ds (1 and 10 mm), leak-point of circular shape with diameter Do (1, 5, and 10 mm), and the pipe burial depth h (1 and 3 m). For all tests, the pipe 66 diameter (Dp) was assumed unchanged and equal to 200 mm since its effect was negligible. To conduct tests, six models indicated in Table 5.4 were made. On the Table, column 5 (a × b) shows the dimension of computational domain around the leak-pipe; and that was assumed to be 4 × 4 m2 when burial depth was 3 m and 2 × 2 m2 when burial depth was 1 m. Please note that the effect of mesh refining on the results was also evaluated. Instead of 50 mm in M4 (Table 5.4), dl was assumed 25 mm. Comparison showed that changes in the results due to the refinement were negligible. 5.5 Results and Discussions Fig. 5.5 indicates the simulations results as contours for pore pressure (Pa). Fig. 5.5-a shows that pore pressure was hydrostatic at points far from the leak-point. On the other hand, Fig. 5.5-b indicates that pore pressure near the leak-point was affected by negative pressure inside the pipe. A very high pressure gradient was observed around the leakpoint (i.e., the pressure increased from -90 to 0 kPa over 1.5 cm distances from the leakpoint). Fig. 5.6-a highlights the contours for velocity (m/s) next to the leak-point. Also, Fig. 5.6-b shows the corresponding velocity vector at the leak-point. As the figure indicates, velocity decreased by increasing the distance from the edges of the leak-point. The maximum velocity was 5.06 m/s happening at the edge of the leak-point, while a minimum velocity of 1.92 m/s occurred at the leak-point center. Calculating the average velocity normal to the leak point and knowing leak-point cross sectional area, average intrusion rate (Q) was calculated as 4.51x10-5 m3/s. Simulations were repeated for other input data (Table 5.3 and 5.4) to produce further data point. After all the simulations were completed, a regression analysis was carried out to find the relationship among the contributing dimensionless parameters (the function F2). The analysis revealed the following relation: π 4 = x4 (π 1 ) x1 (π 2 ) x 2 (π 3 ) x 3 for each π 5 5.12 To find the coefficients of x1, x2, x3, and x4, an optimization process was carried out by using the Solver Toolbox of Microsoft Excel (Table 5.5). Microsoft Excel Solver uses the generalized reduced gradient (GRG) algorithm for optimizing nonlinear problems. GRG algorithm uses the Newton-Raphson iterations during the search (Arora, 2004). 67 Table 5.3- Characteristics of tests conducted for each constant porosity (ε). Test No. P(Kpa) Ds(m) Do(mm) h(m) 1 -90 0.01 10 3 2 -90 0.01 10 1 3 -90 0.01 1 3 4 -90 0.01 1 1 5 -90 0.01 5 3 6 -90 0.01 5 1 7 -90 0.001 10 3 8 -90 0.001 10 1 9 -90 0.001 1 3 10 -90 0.001 1 1 11 -90 0.001 5 3 12 -90 0.001 5 1 13 0 0.01 10 3 14 0 0.01 10 1 15 0 0.01 1 3 16 0 0.01 1 1 17 0 0.01 5 3 18 0 0.01 5 1 19 0 0.001 10 3 20 0 0.001 10 1 21 0 0.001 1 3 22 0 0.001 1 1 23 0 0.001 5 3 24 0 0.001 5 1 Table 5.4- Geometry of the developed models dl at leak-point Model No. Dp (m) h (m) Do (mm) a × b* (m) dl (mm)** M1 0.2 3 10 4 ×4 50 0.05 M2 0.2 1 10 2 ×2 50 0.1 M3 0.2 3 1 4 ×4 50 0.05 M4 0.2 1 1 2 ×2 50 0.1 M5 0.2 3 5 4 ×4 50 0.05 M6 0.2 1 5 2 ×2 The dimension of computational domain around a pipe * ** (mm) 50 0.1 The length of elements 68 a) b) Fig. 5.5- Pressure contour (Pa) around the pipe during occurrence of negative pressure a) b) Fig. 5.6- Velocity contour (m/s) (a) and velocity vector (b) at leak-point. 69 The objective function of the optimization process was to maximize the coefficient of correlation (R2): Maximize R2 for each π 5 5.13 π i 4 − π i 4* / π i 4 ≤ δ Subject to: for i=1,.., 24; for each π 5 where π i 4* is the predicted value of π i 4 computed by Equation (5.12). The constraints in Equation (5.13) keep the relative error less than δ for all the predictions. In fact, both maximum value of R2 and low values of relative error are desired during the optimization process. The parameter δ equal to 20 percent was the smallest value of δ for which there was a feasible solution for optimization. The maximum R2 for π 5 =0.25, 0.3, and 0.4 were found 0.996, 0.993, and 0.995, respectively. To find more explicit relationship for intrusion rate, Equation (5.12) was rewritten as: Q = CDsm1 Aom 2 ∆P m 3 5.14 Table 5.6 shows the values of C, m1, m2, and m3 for different values of ε. Moreover, the orifice equation is: Q = CAo ∆P 0.5 5.15 Comparing Equation (5.14) and (5.15), it is clear that the exponent of ∆P in Equation (5.14) is (about 0.7) higher than that in the orifice equation (0.5). However, the effect of Ao in Equation (5.14) (about 2/3) is lower than that in the orifice equation. Table 5.6 also shows that the effect of constant C decreases by decreasing ε. On contrary, the effect of the average diameter of soil particles (Ds) on intrusion increases by decreasing ε. Note that ε is an external variable in the Equation (5.14). Thus, the values of C, m1, m2, and m3 are to be calculated for every ε. Table 5.5- The optimized coefficient in Equation (5.12) ε x1 x2 x3 x4 0.4 0.159 0.796 -0.076 0.046 0.3 0.177 0.796 -0.064 0.025 0.25 0.206 0.805 -0.053 0.016 70 Table 5.6- The constant values in Equation (5.14) ε C m1 m2 m3 0.4 0.00236 0.584 0.758 0.684 0.3 0.00165 0.633 0.753 0.685 0.25 0.00137 0.685 0.753 0.697 To further verify the predictability of the proposed model, 16 additional hypothetical intrusion-tests were conducted on a 500-mm-diameter pipe buried at 2 m depth. Table 5.7 provides additional details for the tests. Computational domain was a square of size 4 m. Two different circular leak-points were assumed at the top of the pipe with diameters equal to 3 and 6 mm. Different negative pressures were studied including -90, -70, -50, and -30 kPa. Four different soils were examined. Table 5.7 also indicates further information about the porosity and average diameter of the soil particles. The intrusion rates computed by the numerical model were compared with those by Equation (5.14). The comparison indicated relative errors in the range of 0.4 to 13 percent (Table 5.7). The Table also shows that the relative error decreases by decreasing soil particle diameter and therefore by soil permeability (e.g., tests 1 and 5). Also, orifice equation (Equation 5.15) was used to estimate the intrusion rate in each test in Table 5.7. The results showed that the intrusion rate computed by Equation (5.15) is two to sixteen times higher that obtained by numerical model. 71 Table 5.7- Characteristics of tests were conducted for validation of proposed model Test No. P(Kpa) ε Ds(m) Do(m) h(m) Qs (m3/s)* Qp (m3/s)** %Error 1 2 -90 -70 0.4 0.4 0.01 0.01 3 3 2 2 5.02237E-05 4.5281E-05 5.66759E-05 4.93763E-05 12.8 9.0 3 4 -50 -30 0.4 0.4 0.01 0.01 3 3 2 2 3.97364E-05 3.32846E-05 4.15387E-05 3.29445E-05 4.5 1.0 5 6 -90 -70 0.4 0.4 0.001 0.001 3 3 2 2 1.48535E-05 1.32935E-05 1.47805E-05 1.28768E-05 0.5 3.1 7 8 -50 -30 0.4 0.4 0.001 0.001 3 3 2 2 1.1547E-05 9.51232E-06 1.08328E-05 8.59159E-06 6.1 9.6 9 10 -90 -70 0.3 0.3 0.01 0.01 6 6 2 2 8.48817E-05 7.65394E-05 9.5793E-05 8.34458E-05 12.8 9.0 11 12 -50 -30 0.3 0.3 0.01 0.01 6 6 2 2 6.71118E-05 5.61583E-05 7.01903E-05 5.56577E-05 4.5 0.9 13 14 -90 -70 0.3 0.3 0.001 0.001 6 6 2 2 2.25515E-05 1.99403E-05 2.2326E-05 1.94483E-05 1.0 2.4 1.70577E-05 1.36665E-05 1.63589E-05 1.29719E-05 4.1 5.1 15 -50 0.3 0.001 6 2 16 -30 0.3 0.001 6 2 * The intrusion rates computed by the numerical model ** The intrusion rates computed by the proposed model It should be emphasized that this research assumed an intrusion event occurring under steady low/negative pressure. However, a low/negative pressure event in a real-life WDS can be transient in nature. Thus, to estimate intrusion rate in a real-life WDS, either a transient analysis or an extended period simulation (EPS) must be conducted. This is explained by using the experimental data reported by Collins et al. (2011). They employed a 50 mm diameter pipe with an orifice or leak-point (6.35 mm diameter) to study intrusion. To study surrounding soil effects on intrusion, a measuring cylinder containing 5.5-9.5 mm gravel was installed on the top of the orifice. The height of soil column was 290 mm. To keep the soil saturated, the cylinder filled with water (Collins et al., 2011). While the initial pipe pressure was 20 m, a transient low/negative pressure was induced in the pipe by a valve closure scenario. Fig. 5.7 indicates that the pressure eventually reached a new steady state condition of zero pressure after an initial drop of 10 m. Collins et al. (2011) measured 113 ml for the total volume of the intruded water during the transient low/negative pressure (Fig. 5.7). Further details about the experiment and soil properties can be found in Collins et al. (2011). Given the soil porosity, the proposed ingress model, Equation (5.14), was used to estimate intrusion volume. In this regard, the transient pressure given in Fig. 5.7 was divided into some 72 time intervals through which the low/negative pressure was assumed to be steady. Then, Equation (5.14) was used to estimate the intrusion flow rate in each time interval. The total volume of intrusion was ultimately determined by the summation of the multiplication of each intrusion flow rate and its corresponding time interval over the intrusion duration. This results an intrusion volume of 101mL. The relative error between the results of Equation (5.14) and the experimental measurement by Collins et al. (2011) was around 11% indicating an acceptable agreement between both sets of the results. Fig. 5.7- Transient pressure induced by valve closure (Collins et al., 2011) Summery In this chapter, a model was proposed to study the effect of surrounding soil on the intrusion rate. To develop the model, computer simulations were conducted instead of lab experiments. A two-factorial design technique was employed to screen the effect of contributing factors. The results showed that the effect of pipe diameter was negligible, while the leak-point diameter was the most influential factor. As well, the effect of pipe internal pressure, soil porosity, and average diameter of soil particles were found to be almost equally significant. The proposed model was similar to orifice equation; however, the characteristics of the surrounding soil including the porosity and average diameter of soil particles were implemented. The contributions of factors in the proposed model were different compare to those in orifice equation. The proposed model was verified using a new set of numerical simulation data. The results showed that the relative errors between 73 numerical and predicted values were in an acceptable range. The applicability of the proposed model to estimate intrusion rate during transient low/negative pressure was shown by comparing the results obtained from the proposed ingress model coupled with EPS technique and results from an earlier experimental study in the literature. 74 Chapter 6: Lagrangian-Eulerian Transient Water Quality Model A version of this chapter was accepted for publication in the Journal of American Water Work Association (AWWA). Mansour-Rezaei, S., Naser, Gh., Malekpour, A., Karney, B., W. (2012). "Contaminant Intrusion in Water Distribution Systems – A LagrangianEulerian Transient Water Quality Model", Paper No. 08232012-JAWWA0104, to be published in June 2013. Preface The available water quality models largely ignore inertia and compressibility effects and assume either steady or quasi-steady conditions for flow. EPANET is a typical example. Yet a real WDS may undergo transient conditions in which both fluid inertia and compressibility effects are significant. Steady or quasi-steady models provide an overall understanding of flow in a WDS. However, continuous demand variations along with operations of local devices (valves, pumps, etc.) create transients that cannot be captured by steady or quasi-steady models. A realistic model can more reliably define what operating strategy must be undertaken to ensure acceptable water quality in an efficient manner. In this chapter, a WQM is proposed to study contaminant intrusion and its propagation along a WDS by coupling an Eulerian-based transient hydraulic model (for detecting junctions of a WDS experiencing negative pressures), an ingress model (for estimating rate of contaminant entering a system during an intrusion event), and a Lagrangianbased contaminant transport model (for tracking the fate of a contaminant along a system). Compare to Eulerian-based contaminant transport model, Lagrangian methods are more computationally efficient in terms of CPU simulation time and memory usage. Lagrangian methods are also free of numerical dissipation and dispersion errors. To show the applicability of the proposed WQM, the model is employed to simulate intrusion in two well-known WDSs reported in the literature. Investigating the hazard and vulnerability of the systems to contaminant intrusion are investigated, the fate of intruded contaminant under water demand uncertainty are evaluated. 75 6.1 Lagrangian-Eulerian Transient Model Fig. 6.1 presents a flowchart for the proposed model. The model consists of an Eulerianbased transient hydraulic model, an ingress model, and a Lagrangian-based contaminant transport model. The following sections discuss these models. Start of Lagrangian–Transient WQM Transient Hydraulic Model Update V and H at Junctions Update V and H at Pipes Contaminant Intrusion Model Detect Negative Pressure at Junction i Estimate Intrusion Rate at Junction i Q=AC0 |Hi| 0.5 Lagrangian Contaminant Transport Model Update C at Junctions Update C at Pipes End of Lagrangian–Transient WQM Fig. 6.1 - Flowchart of the Lagrangian-based transient WQM. 6.1.1 Eulerian-Based Hydraulic Model Applying an Eulerian model, this research determined the hydraulic conditions of a flow (under a transient event) through the solutions of a hyperbolic set of partial differential equations of continuity and momentum (Boulos et al., 2005 and Larock et al., 2000): ∂H a 2 ∂V + =0 ∂t g ∂x ∂V ∂H fV V +g + =0 ∂t ∂x 2d 6.1 6.2 where V is flow velocity (LT-1), H is piezometric head (L), d is internal pipe diameter (L), g is acceleration (LT-2) due to gravity, f is Darcy-Weisbach friction factor, a is acoustic wave speed (LT-1), x is the distance (L) along the pipe centerline, and t is time 76 (T). Note that a is a function of fluid compressibility, pipe elasticity and its diameter and thickness. While Larock et al. (2000) and Boulos et al. (2005) provided further details about MOC, a fixed-grid MOC replaces the partial differential Equations of 6.1 and 6.2 by an equivalent set of ordinary differential equations (ODE). Using a finite difference scheme, the resulting equations were then integrated along two characteristic lines: dx dt = ±a 6.3 A key to note is that ∆t in the MOC solution must be common to all the pipes. Indeed, a common ∆t forces the characteristic lines to meet at the pipes’ junctions. As Equation 6.3 suggests, MOC provides satisfactory results only if the Courant number (i.e., the ratio of wave speed over numerical speed) is one. Moreover, Karney and Ghidaoui (1997) argued that the Courant numbers less than one increase numerical dissipation. It is, however, impossible to satisfy Courant number one for all the pipes in a real-life WDS. Thus, either the wave speed or pipe length is adjusted or an interpolation technique is used (Fernandes and Karney 2004; Karney and Ghidaoui, 1997). Typical boundary devices in a WDS may include reservoirs and tanks, pumps, valves, and air chambers. The boundary devices were implemented in the hydraulic model as externally imposed conditions on velocity and/or pressure head. Chaudhry (1979) provides additional details on modeling boundary devices. 6.1.2 Ingress Model Intrusion may occur when the internal piezometric head (Hi) is lower than the external piezometric head. The intrusion flow rate at a typical junction i of a WDS can be estimated by using either the original orifice equation (Besner et al., 2011 and McInnis, 2004) or the equation proposed by Mansour-Rezaei and Naser (2012). The model of Mansour-Rezaei and Naser (2012) can provide a more realistic estimate for intrusion rates since it takes into account the characteristics of soil surrounding a pipe. However, this research employed the original orifice equation due to a lack of detailed information about the surrounding environment for two WDSs cases. The original orifice equation can be stated as: Qi* = CD A 2 g ∆H i 6.4 77 where Qi* is intrusion flow rate (L3T-1) at junction i, CD is the dimensionless discharge coefficient, A is orifice area (L2), g is acceleration (LT-2) due to gravity, and ∆Hi is the piezometric head difference across the junction i. 6.1.3 Lagrangian-Based Contaminant Transport Model The transport of a conservative contaminant is modeled by an advection-diffusion equation. Due to the high velocities in most pipes, advection is often the dominate transport mechanism (Boulos et al., 2005 and Naser and Karney, 2007 and 2008). Thus, ignoring diffusion, the transport of a conservative contaminant was modeled by: ∂C ∂C +V =0 ∂t ∂x 6.5 where C is contaminant concentration (ML-3). Note that the contaminant was assumed to be conservative for the sake of simplicity, although the model can be easily modified for a nonconservative contaminant. The research considered no changes in contaminant fate at a partially- or fully-open valve and operating pump. A fully-closed valve or a pump in off-mode was modeled as a no-flux boundary. Suppose a typical junction i (of a WDS) was contaminated during n successive hydraulic time steps. The junction i is contaminated when a pollutant entered the junction from outside the pipe or reached the junction from other pipes joining it. Assuming complete and instantaneous nodal mixing, the contaminant concentration at junction i was calculated by: m ∑Q C j Ci = j + Qi*Ci* j =1 m ∑ Q j + Qi* 6.6 j =1 where Ci is concentration (ML-3) leaving junction i, Qj is inflow discharge (L3T-1) to the junction i by pipe j, Cj is concentration (ML-3) at the end of pipe j, Qi* is intrusion flow rate (L3T-1), Ci* is concentration (ML-3) of intruded contaminant at junction i, and m is the number of pipes joining the junction i. Complete and instantaneous nodal mixing is often accepted in WDSs modelling for the purpose of simplicity, although Ho (2008) confirmed that a complete-mixing assumption is not able to accurately capture the physical mixing process at the junctions. 78 Lagrangian Contaminant Transport Model Start Detect Contaminated Junction i Detect Pipe k at Downstream of Junction i Calculate VCPk and CCPk Determine the Parcel Location (s) in Pipe k Yes ΣVs (tsw) >VCPk /Ak ∆t No Cs (tsw) = CCPk Cs (tsw) =0 ΣVs (tsw) > ∆x/∆t Yes No Cs+1 (tsw) = CCPk Update the Parcel Location (s) in Pipe k Lagrangian Contaminant Transport Model End Fig. 6.2 - Lagrangian contaminant transport model. Considering the contaminated water as a contaminant parcel moving along a pipe in a WDS, this research used a Lagrangian model to solve Equation 6.5. Fig. 6.2 indicates the steps for the Lagrangian-based transport model for a contaminant parcel along a WDS. Thus, the concepts of volume of a contaminant parcel (VCP) and the concentration of a contaminant parcel (CCP) were defined to estimate the fate of a contaminant parcel along a WDS. During the contamination of junction i, a contaminant parcel leaves the junction through pipe k. The VCP and CCP were defined as: 79 n VCPk = ∆tAk ∑Vk (t p ) For t1 < tp < tn 6.7 p =1 n ∑ C (t i CCPk = p ) p =1 6.8 n where VCPk is the volume (L3) of contaminant parcel moving along pipe k, Vk (tp) is the velocity (LT-1) at the beginning of pipe k, tp is a time (T) at which contaminant exists at the upstream junction i of pipe k, Ak is cross sectional area (L2) of pipe k, ∆t is the time step determined by the hydraulic model, Ci (tp) is contaminant concentration (ML-3) at junction i calculated by Equation 6.6 and n is the number of time steps that a contaminant exists at junction i. Fig. 6.3 shows a typical pipe k and its upstream and downstream junctions i and r. The figure also indicates the internal numerical nodes assigned by the hydraulic model. If the contaminant parcel in pipe k reaches the internal node s at time ts1, at this time contaminant concentration at point s changes from zero to CCPk. To determine the contaminant concentration at node s at time tsw, ΣVs(tsw) was calculated where tsw shows the time equal to wth hydraulic time step after reaching contaminant to the internal node s. When this summation exceeds VCPk/(Ak∆t), the contaminant parcel moves away from node s. Therefore, the concentration at node s changes from CCPk to zero. When this summation exceeds ∆x/∆t, the contaminant parcel reaches the node s+1 and the concentration at node s+1 changes from zero to CCPk. As the time advances, the contaminant parcel reaches the junction r downstream of pipe k. At this time, the contaminant parcel may be divided into new parcels depending on the pipes’ flow rates leaving the junction. While a parcel leaves the system through nodal demand at r, the remaining parcels travel along the other pipes leaving junction r. Using Equation 6.7 and 6.8, VCP and CCP can be calculated for the pipes taking water from junction r. It should be noted that when the duration of a contaminant intrusion is relatively long, several contaminant parcels (rather than one) are created at the beginning of each pipe. After this, the fate of the contaminant parcels is simulated using the proposed algorithm. 80 Pipe k Junction i Contaminant parcel s-1 s Junction r s+1 Fig. 6.3 - Contaminant parcel reaching point s at time t1. 6.2 Hazard and Vulnerability Indices When a contaminant intrusion occurs in a WDS, hazardous and vulnerable locations of the system must be detected in order to manage the contaminant threat. This research adopted and employed the hazard and vulnerability indices initially defined by Xin et al. (2012). Hazard indices of a junction were calculated by evaluating the impacts of intrusion events on the downstream junctions. “Cumulative mass (CM)” and “peak concentration (PC)” were the two hazard indices. Assuming intrusion occurred at junction k, cumulative mass (CMk) and peak concentration (PCk) delivered to the other junctions were estimated by: N ,i ≠ k T CM k = ∑ ∑C i =1 ij Q ij ∆t 6.9 j =0 N ,i ≠ k PC k = ∑ i =1 max (C ij Q ij ) j =1→T 6.10 where Cij and Qij are contaminant concentration (ML-3) and demand (L3T-1) at junction i at time j, T is response time, ∆t is time step, and N is the total number of junctions. Vulnerability indices of a junction were calculated by evaluating the impacts of intrusion events on this junction when intrusion occurred at all the other junctions. “Total cumulative mass (TCM)” and “total peak concentration (TPC)” were the two vulnerability indices used in this research. Considering all intrusion scenarios, total cumulative mass (TCMi) and total peak concentration (TPCi) delivered to junction i were estimated by: N ,k ≠i T TCM i = ∑ ∑C jk Q jk ∆t 6.11 k =1 j =0 N , k ≠i TPCi = ∑ max C k =1 j =1→T jk Q jk 6.12 81 To allow comparison among the results, the hazard and vulnerability indices were normalized using the minimum and maximum values calculated for a system by (Xin et al., 2012): G[ Index]k = [ Index]k − min [ Index]k k =1→ N max[ Index]k − min [ Index]k k =1→ N 6.13 k =1→ N where Gk is the normalized value of an index calculated for junction k. Gk equal to 1 shows junction k is the most hazardous/vulnerable junction of the system. 6.3 Model Applications, Results and Discussion Case Study 1 The proposed hydraulic and water quality model was first employed to examine the flow conditions in the WDS initially introduced by Boulos et al. (2005) and later used by Wood et al. (2005) and Jung et al. (2009) as a case study. Fig. 6.4 indicates the schematic view of the WDS. On the figure, the number in a circle indicates the junction number, while that in a square represents the pipe number. Also, the arrows in the figure show the direction of pipe flows before valve closure. The system consists of nine pipes, five junctions, one constant-head reservoir, and a valve. Table 6.1 indicates the diameter, length, friction factor, and wave speed of each pipe. All junctions are located at the same elevation. The reservoir head was 80 m. The demands at junctions 1, 2, 3, 4, and 5 were 0.6, 0.5, 0.5, 0.6 and 0.5 cms, respectively. The demand at the valve located at the downstream end of the pipe 7 was 0.5 cms. The system operated for 50 s resulted in the steady state discharges of 3.2, 1.42, 1.28, 0.52, 0.41, 0.84, 0.5, 0.26, and 0.25 cms in the pipes 1, 2, 3, 4, 5, 6, 7, 8, and 9, respectively. A transient scenario was created in the network by gradual closure of the valve over one second. The local-loss coefficient of the valve (KL) for 100, 80, 60, 40, and 20 percent opening was assumed to be 0.19, 1.15, 5.6, and 31.9, respectively. Given the large demands, pipe flows and large diameters in the test network, it seems that the network is a skeletonised network and only large trunks were considered. 82 1 4 3 5 2 2 5 8 9 3 1 4 6 7 Fig. 6.4 - Schematic view of the network for case study 1. To find ∆t and ∆x for the purpose of hydraulic modelling, the minimum required number of sections (N) for the shortest pipe was defined. Then, ∆x was set as the pipe length divided by N. By knowing a and ∆x and assuming Courant number to be unit, ∆t was calculated for the shortest pipe. The calculated ∆t was considered as the common (fixed) ∆t for the hydraulic modelling. After defining the fixed ∆t, the number of division for the other pipes were calculated such that the Courant number to be as close as to 1. Table 6.1 - Pipes characteristics for case study 1. Pipe Number Diameter (mm) Length (m) Friction factor Wave speed (a) (m/s) 1 900 610 0.02 1000 2 750 914 0.02 1000 3 600 610 0.02 1000 4 450 457 0.02 1000 5 450 549 0.02 1000 6 750 671 0.02 1000 7 900 610 0.02 1000 8 600 457 0.02 1000 9 450 488 0.02 1000 Fig. 6.5 shows the variations of piezometric head (H) with time at junctions 1, 2, 3, 4, and 5. The figure clearly indicates that the negative pressures at junctions 1, 3, and 4 resulted from the valve closure. The pressure at junction 1 dropped from 41.4 to -3 m at time 58.9, for junction 3 dropped from 49.7 m to -6.3 m at time 59.8 s, and for junction 83 4 the drop was from 38.1 m to -6.8 m at time 58.6 s. The duration of negative pressure at junction 1 was 2.2 s, for junction 3 was 0.9 second, and for junction 4 was 0.8 s. The WDS reached a new steady-state hydraulic condition after around 100 s. The study indicated 2.07, 1.18, 1.02, 0.38, 0.3, 0.54, 0.0, 0.06 and 0.18 cms for the new steady state discharges at pipes 1, 2, 3, 4, 5, 6, 7, 8, and 9, respectively. In the presence of a pathway and a contaminant source, contaminant may intrude the system at junctions 1, 3 and 4 due to negative pressures. For simplicity, this research assumed that the pathway and contaminant source existed only in the vicinity of junction 3. Therefore, intrusion occurred at this junction. Knowing the magnitude of negative pressure, Equation 6.4 was used to estimate intrusion rate at junction 3 at each time. In Equation 6.4, C0A was assumed equal to 0.01. Equation 6.6 was used to calculate contaminant concentration at junction 3 during the intrusion while external contaminant concentration was assumed 30 mg/m3. Fig. 6.6 shows the time variations of contaminant concentration at junction 3 as well as the piezometric head at this junction. The figure clearly shows two facts; 1) intrusion occurs during the period of negative pressure at junction 3 and 2) the lower negative pressure causes more contaminant intrusion. Due to contaminant intrusion at junction 3, two parcels of contaminant were created at beginning of pipes 4 and 5 traveling with the flow along the two pipes. Using Equation 6.7 and 5.8, VCP4, CCP4, VCP5 and CCP5 were calculated at 0.36 m3, 0.49 mg/m3, 0.30 m3 and 0.49 mg/m3, respectively. VCP4 was larger than VCP5. This was expected since pipe 4 supplied more water than pipe 5. Note also that CCP5 (or CCP4) is the average of contaminant concentration at junction 3 during the intrusion event (Fig. 6.6). 84 Fig. 6.5 - Time variations of head at junctions 1, 2, 3, 4, and 5 of case study Fig. 6.7 shows the time variations of contaminant concentrations at junctions 1, 4, and 5. Note that VCP4 reached junction 5 at time 249.9 s and moved away from this junction at time 250.8 s. The contaminant concentration at junction 5 during this time period was 0.337 mg/m3. Due to the contamination of junction 1, a parcel of contaminant was created at the beginning of pipe 8 with VCP8 and CCP8 as 0.055 mg/m3 and 0.34 m3. This parcel reached junction 4 at time 2472.8 s and moved away from this junction at time 2473.7 s. The contaminant concentration at junction 4 during this time period was 0.033 mg/m3. VCP5 reached junction 1 at time 348.5 s and moved away from this junction at time 349.5 s. The contaminant concentration at junction 1 during this time period was 0.113 mg/m3. Due to the contamination of junction 5, two parcels of contaminant created at the beginning of pipes 6 and 9 with VCP6 and VCP9 of 0.54 and 0.18 m3 and CCP6 and CCP9 of 0.113 mg/m3. VCP6 reached junction 4 at time 895.9 s 85 and moved away from the system at time 896.9 s. The contaminant concentration at junction 4 during this time period was 0.10 mg/m3. VCP9 reached junction 5 at time 786.5 s and moved from this junction at time 787.5 s. The contaminant concentration at junction 5 during this time period was 0.036 mg/m3. Another parcel of contaminant created at the beginning of pipe 8 with VCP8 and CCP8 equal to 0.06 m3 and 0.036 mg/m3, respectively. This parcel of contaminant reached junction 4 at time 3007.8 s and moved away from the system at time 3008 s. The contaminant concentration at junction 4 during this time period was a mere 0.003 mg/m3 almost certainly associated with the short duration of the occurrence of low/negative pressure. Fig. 6.6 - Variations of head and concentration at junction 3 of case study 1. To demonstrate the applicability of the proposed model to evaluate safety of a WDS in the cases of contaminant intrusion events, hazard and vulnerability were evaluated for this case study. This research employed hazard and vulnerable indices defined by Equation 6.9 to 5.13 initially introduced by Xin et al. (2012). All possible intrusion scenarios were evaluated using the proposed WQM (all the junctions experience negative pressure including junction 1, 3, and 4). Table 6.2 shows the results for G(CM), G(CP), G(TCM), and G(TPC). The results showed that junction 1 and 3 were two hazardous junctions of the system in the case of contaminant intrusion. Although junction 4 experienced negative pressures, the hazard indices for this junction were zero 86 since discharge at the valve located at the downstream of junction 4 is zero after valve closure (Qij is zero in Equation 6.9 and 5.10). According to Table 6.2, G(CM) index indicates junction 1 as the most hazardous junction, while G(CP) index defines junction 3 as the most hazardous junction. It means cumulative mass delivered to junctions located at downstream of junction 1 was higher than the peak concentration delivered to those junctions. Note that the magnitude of negative pressure at junction 1 (-3m) was greater than that at junction 3 (-6.3 m). Thus, concentration of a contaminant parcel created at the beginning of junction 1 was expected to be lower than that at junction 3. Also, total discharge coming to junction 1 from pipe 3 (0.54 cms) and pipe 5 (1.02 cms) was higher than that coming to junction 3 from pipe 2 (1.02 cms) causing more dilution of the intruded contaminant at junction 1 compared to that at junction 3. The G(CM) at junction 1 was higher than that of junction 3 due to the fact that the duration of negative pressure (intrusion) at junction 1 (2.2 s) was higher than that at junction 3 (0.9 s). The last two columns in Table 6.2 show vulnerability indices. Junction 5 was the most vulnerable junction of the system because it was located right after both hazardous junctions 1 and 3. Junction 2 was not vulnerable as it was not affected by intrusion events. Although junction 3 was a hazardous junction, it was not vulnerable to intrusion. The vulnerability of junction 1 was less than that at junction 4. Junction 1 was affected by intrusion at junction 3, while junction 4 was affected by intrusion at both junctions 1 and 3. Table 6.2 - Hazard and vulnerability evaluation of case study 1. Junction number 1 2 3 4 5 Hazard Indices G(CM) G(CP) 1 0 0.19 0 0 0.23 0 1 0 0 Vulnerability indices G(TCM) G(TPC) 0.31 0 0 0.93 1 0.37 0 0 0.58 1 87 Fig. 6.7 - Variations of concentration at junctions 1, 4, and 5 due to the intrusion at junction 3 for case study 1. 88 Case Study 2 The proposed WQM was also used to examine a WDS studied by Boulos et al. (2005). Fig. 6.8 indicates the schematic view of the system with the arrows showing the direction of pipe flows before the pump trip. The system consists of twenty one pipes, sixteen junctions, one constant-head reservoir, two tanks, and a pump. Table 6.3 indicates the characteristics of the pipes including diameter, length, friction factor, and wave speed. All junctions were located at the same elevation. The reservoir head was 70 m, while the heads at both tanks were 70 m. The demand at junctions 1, 2, 5, 6, 8, 9, 10, 11, and 13 were 0.025 cms, while demand was equal to zero at junctions 3, 4, 7, 12, 14, 15 and 16. After 200 s of steady state operation, a transient scenario was created in the network though the failure of the supply the pump. Table 6.3 also indicates the steady state discharges in the pipes before and after the pump failure. The negative discharges showed that the direction of flow reversed after pump failure. The The proposed transient hydraulic model was used to simulate flow condition in the network. For hydraulic modelling, ∆t and ∆x were equal to 5.16×10-2 s and 62 m, respectively. Fig. 6.9 shows time variations of piezometric head at these junctions. The figure showed negative pressures in junctions 2, 6, 11, and 14 after pump failure. The piezometric head at junction 2 dropped from 28.7 m to -0.9 m at time 205.7 s and then after some oscillations it dropped to -1.7 m another time at 209.4 s. The piezometric head at junction 6 dropped from 31.8 m to -1.9 m at time 207.2 s. At junction 11, pressure dropped from 31.3 m to -1.4 m at time 207 s. At junction 14, pressure dropped from 29.2 m to -5.3 m at time 207 s. The duration of negative pressure at junctions 2, 6, 11, and 14 was 4.0, 0.4, 0.3, 0.7, and 0.8 s, respectively. Note that the duration of negative pressure at junction 2 was longer than that for the other junctions. 89 Table 6.3 - Pipes characteristics for case study 2. Pipe Diameter Number (mm) 1 150 2 250 3 300 4 200 5 200 6 150 7 203 8 200 9 150 10 200 11 150 12 200 13 200 14 200 15 250 16 200 17 200 18 200 19 200 20 150 21 300 *discharge before pump failure Wave speed (a) Length Friction (m/s) (m) factor 867 0.02 1200 628 0.025 1200 62 0.025 1200 712 0.015 1200 1005 0.015 1200 604 0.02 1200 209 0.015 1200 803 0.015 1200 956 0.02 1200 502 0.015 1200 854 0.02 1200 446 0.015 1200 731 0.015 1200 777 0.015 1200 534 0.025 1200 467 0.015 1200 896 0.015 1200 592 0.015 1200 1051 0.015 1200 387 0.02 1200 79 0.025 1200 ** discharge after pump failure Q1 * (m3/s) 0.024 0.025 0.1 0.075 0.074 0 0.049 0.076 0 0.04 0 0.026 0.05 0.024 0.051 0.035 0.011 0.076 0.01 0.025 0.1 Q2** (m3/s) -0.01 0.025 0 0.114 0.079 0 0.054 0.071 0 0.039 0 -0.04 0.111 -0.01 -0.015 0.036 0.007 0.071 0.011 0.025 0 90 12 5 20 6 3 11 14 19 16 2 Tank 1 (70 m) 13 2 11 7 5 4 9 17 10 6 1 Reservoir 14 (27m) 18 15 8 7 16 9 Pump 1 4 13 12 8 10 15 21 3 Tank 2 (70m) Fig. 6.8 - Schematic view of the water network for case study 2. 91 Assuming the existence of a pathway and contaminant source only in the vicinity of junction 6, two parcels of contaminated water were created at the beginning of pipes 10 and 17 traveling with the flow along the two pipes. VCP10, CCP10, VCP17 and CCP17 were calculated at 0.02 m3, 3.98 mg/m3, 0.005 m3 and 3.98 mg/m3, respectively. Fig. 6.10 shows time variations of contaminant concentration at junctions 2, 5, 9, 11, and 14. The contaminant parcel moving along pipe 10, reached junction 2 at time 611.1 s and moved away from this junction at time 611.5 s. The contaminant concentration at junction 2 during this time period was 3.11 mg/m3. Note that the concentration of contaminant parcel at pipe 10 was 3.98 mg/m3, while the concentration at junction 2 was 3.11 mg/m3. Clean water coming to junction 2 by pipe 19 diluted the contaminant. When contaminated water reached junction 2, a parcel of contaminant was created at the beginning of pipe 2 with VCP2 and CCP2 as 0.01 m3 and 3.11 mg/m3. This contaminant parcel reached junction 9 at time 1844.2 s and removed from the junction at time 1844.5 s. The contaminant parcel moving along pipe 17 (VCP17) reached junction 11 at time 4313.6 s and moved away from this junction at time 4314.2 s. The contaminant concentration at junction 11 during this time period was 0.45 mg/m3. Note that the concentration of the contaminant at junction 11 was diluted by the clean water coming to the junction through pipe 7. When contaminated water reached junction 11, a parcel of contaminant was created at the beginning of pipe 16. VCP16 and CCP16 were equal to 0.02 m3 and 0.45 mg/m3, respectively. This contaminant parcel reached junction 14 at time 4723.4 s and was removed from the junction at time 4724 s. The contaminant concentration at junction 11 during this time period was 0.45 mg/m3. Two parcels of contaminants was created and moved along pipe 20 and 19. VCP20, CCP20, VCP19 and CCP19 were 0.02 m3, 0.45 mg/m3, 0.007 m3, and 0.45 mg/m3, respectively. Contamination of junction 5 occurred between time 4997.7 and 4998.2 s, while contaminant concentration during this time period was 0.45 mg/m3. Contamination of junction 2 occurred between time 7762.7 and 7763.3 s, while contaminant concentration during this time period was 9.76×10-2 mg/m3. The last parcel of contaminant was created at the beginning of pipe 2 with VCP2 and CCP2 of 0.02 m3 9.76×10-2 mg/m3, respectively. This contaminant parcel reached junction 9 at time 8998.2 s and was removed from the junction at time 8998.7 s. Note that the contaminant never reached junction 3 because the discharge at pipe 11 was zero (Table 6.3). 92 Table 6.4 - Hazard and vulnerability evaluation for case study 2. Junction number 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Hazard Indices G(CM) G(CP) 0 0.62 0 0 0 0.30 0 0 0 0 0.20 0 0 1 0 0 0 0.04 0 0 0 0.51 0 0 0 0 0.41 0 0 1 0 0 Vulnerability indices G(TCM) G(TPC) 0 0.33 0 0 1 0 0 0 1 0 0.04 0 0 0 0 0 0 0.43 0 0 1 0 0 0 0.47 0 0.03 0 0 0 0 0 Vulnerability and hazard were evaluated for this case study. Table 6.4 shows the results for G(CM), G(CP), G(TCM), and G(TPC). The table shows junctions 2, 6, 11, and 14 as the hazardous junctions in the WDS. This was expected as these junctions were the only junctions of the system that experienced negative pressures. The table also shows junction 14 to be the most hazardous junction based on both G(CM) and G(CP) indices. Hydraulic model simulations showed that both duration and magnitude of negative pressure at junction 2 were highest at this junction. Thus, concentration peak as well as cumulative mass delivered to the downstream junction was very high. Junction 2 was ranked as the second most hazardous junction based on G(CM) index; however, it was ranked as the least hazardous junction based on G(CP) index. The high value of G(CM) at junction 2 implies that a considerable amount of contaminant reached the downstream junctions if intrusion occurred at this junction. Note that duration of negative pressure at junction 2 was longer that that at the other junctions. However, the low value of G(CP) shows that the concentration of contaminant reaching the downstream junctions was not very high. The reason was the dilution at junction 2 due to the large incoming discharge from pipe 19 (see Table 6.3). Generally, junctions located downstream of the hazardous junctions were vulnerable. As Table 6.4 shows based on G(TCM) and G(TPC) indices, junctions 2, 5, 9, and 11 were vulnerable. Junctions 14 and 3 were expected to be vulnerable but due to zero demand at these junctions, vulnerability indices were zero. Table 6.4 shows junction 5 to be the most vulnerable junction based on both G(TCM) 93 and G(TCP) indices. This results was expected since junction 5 was located downstream of junction 14 (e.g., the most hazardous junction). As well, G(TCM) was higher than G(TPC) for junction 9. Considering the given explanation about hazard indices of junction 2 located at the upstream of junction 9 reveals the reason of such results. Fig. 6.9 - Variations of head at junctions 2, 6, 11, and 14 in case study 2. 94 Fig. 6.10 - Variations of concentration at junctions 2, 5, 9, 11, and 14 for case study 2. 95 Case Study 3 In this part, the proposed WQM was used to evaluate the fate of intruded contaminant under uncertainty. Considering water demand as uncertain parameter, the WDS introduced in Fig. 6.8 was studied. To propagate uncertainty related to water demand, fuzzy set theory (Chapter 3) was employed. It was assumed intrusion may occur at junction 1 where pathway and contaminant source exist. An external contaminant around the junction was set at 90 mg/m3. Table 6.5 - The lowest, most likely and highest demand at junctions. Junction Number 1 2 3 4 5 valve Highest Value (cms) 0.55 0.45 0.45 0.55 0.45 0.45 Most Likely Value (cms) 0.6 0.5 0.5 0.6 0.5 0.5 Lowest Value (cms) 0.65 0.55 0.55 0.65 0.55 0.55 Due to occurrence of intrusion at junction 1, junctions 4 and 5 were contaminated. Fig. 6.11 shows the results of the WQM for temporal variations of concentration at junctions 4 and 5. To increase credibility of the results, fuzzy set theory was employed to propagate uncertainties associated with demands. Table 6.5 indicates the lowest (L), most likely (M), and the highest (H) values of demands for each junction of the WDS. After substituting the fuzzy demands into the water quality model, the vertex method (Section 3.2.1) was applied to compute the contaminant concentration at each junction of the WDS. It should be noted that vertex method assumes independent fuzzy variables. While considering a large number of α-levels increases accuracy of the results, this research applied the fuzzy analyses with three α-levels of 0, 0.5, and 1 (step 2 in Section 3.2.1) to reduce computational burden. At α-level=0, on each junction of the WDS, the demands had a higher and lower values (step 3 in Section 3.2.1). As a result, 26 numbers of demands combinations were considered (step 4 in Section 3.2.1). These combinations were implemented to the water quality model to compute the higher and lower values of concentrations at α-level=0. Fig. 6.12 shows the results of 26 simulations as the variation of concentration with time at junctions 4 and 5. The results showed that at α-level=0, intrusion did not occur in 53% of the cases. However, when intrusion occurred, the consumers at junctions 4 and 5 might be exposed to the contaminated water with maximum concentration of 0.74 and 0.27 mg/m3, respectively 96 (step 5 in Section 3.2.1). Similarly, at α-level=0.5, 26 demand values were implemented to the water quality model to compute the higher and lower values of concentration (step 6 in Section 3.2.1). The results showed that at α-level=0.5, intrusion did not occur in 11% of the cases. However, when intrusion occurred, the consumers at junctions 4 and 5 were exposed to the contaminated water with maximum concentration of 0.53 and 0.19 mg/m3, respectively. At α-level=1, each junction of the WDS had a single value as the most likely value of the demand. Therefore, only one simulation was carried out to compute the most likely value of the contaminant concentration at each junction of the WDS (step 6 in Section 3.2.1). The results revealed that at α-level=1 the consumers at junctions 4 and 5 were exposed to the contaminated water with concentration of 0.22 and 0.077 mg/m3, respectively. These values were the most likely values of contaminant concentration at junctions 4 and 5. Table 6.6 summarizes the possibility of contamination of junctions 4 and 5 due to the intrusion at junction 1. The lowest possible values of contaminant concentrations at junctions 4 and 5 were zero; the most likely values of contaminant concentrations at junctions 4 and 5 were 0.53 and 0.19 mg/m3, respectively; the highest possible values of contaminant concentrations at junctions 4 and 5 were 0.74 and 0.27 mg/m3, respectively. It was observed that increasing uncertainties in demands (decreasing α-level) caused more severe negative pressures in the system that could potentially lead to more sever contamination at junctions 4 and 5. Table 6.6- Fuzzy analysis for maximum contaminant concentrations at junctions 4 and 5. α-level Not intrusion Intrusion Maximum concentration at junction 4 (mg/m3) Maximum concentration at junction 5 (mg/m3) 1 0% 100% 0.22 0.077 0.5 11% 89% 0.53 0.19 0 53% 47% 0.74 0.27 97 Fig. 6.11-Time variations of contaminant concentration at junction 4 (top) and 5 (bottom) due to intrusion at junction 1. Fig. 6.12- Time variations of contaminant concentration at junction 4 (top) and 5 (bottom)- 26 numbers of simulations at α-level=0. 98 Summery To study a contaminant intrusion event in a WDS, this research coupled a Lagrangianbased contaminant transport model to an Eulerian-based transient hydraulic model. The hydraulic model properly detected negative pressures induced by valve closure or pump failure scenarios. While the WQM enjoyed all the benefits of Lagrangian solutions (no numerical dissipation and dispersion errors), it tracked the fate of an intruded contaminant along a system. The combined model was applied to study contaminant intrusion in two WDSs reported in the literature. The results were then used for a hazard and vulnerably evaluation of the systems. The proposed WQM provided a tool to evaluate the vulnerably of a WDS in contaminant intrusion events. Such tool will help engineers and water authorities to 1) study hazard and vulnerability of a WDS to contaminant intrusion, 2) estimate risk associated with a contaminant intrusion (the combination of hazard and vulnerability evaluations can be used to indicate the concept of risk), and 3) make more informed decisions in selecting the most efficient operating strategy to mitigate the risk. 99 Chapter 7: Conclusion, Contributions, and Recommendations for Future Research 7.1 Conclusion The proposed models in this dissertation allows water authorities and designers to: 1) verify potential of contaminant intrusion in a WDS, 2) identify susceptible locations in a system to contaminant intrusion, 3) estimate intrusion flow rate during transient events, and 4) determine the propagation of contaminant (in the case of intrusion occurrence) throughout a WDS. This will eventually increase safety and reliability of a potable water supply system. The following list summarizes the most significant outcomes of the current dissertation: 1. Uncertainty Analysis: Performance of commonly used uncertainty analysis techniques including fuzzy set theory, 1-D and 2-D Monte Carlo simulations, and PBox method were compared. 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 as well as fuzzy set is comparatively much less than Monte Carlo simulations. 2. Potential for Contaminant Intrusion: Studying a WDS with cast-iron pipes, this research developed a fuzzy rule-based model to determine the potential of contaminant intrusion at each junction of the WDS. The fuzzy rule-based model was able to prioritize the junctions of a WDS with respect to contaminant intrusion potential using vague and uncertain data. Such information will help engineers and water authorities to study susceptible locations of a WDS to contaminant intrusion. 3. An Ingress Model: A model was developed for more realistic estimate of intrusion rates into a pipe through a leak-point. The model was similar to orifice equation; however, the characteristics of the surrounding soil including the porosity and 100 average diameter of soil particles were incorporated. Employing the developed ingress model can help to perform more accurate exposure and risk assessment. . 4. WDS Modeling: A Lagrangian-based contaminant transport model was coupled to an Eulerian-based transient hydraulic model to study contaminant intrusion and its propagation in a water distribution system. The transient hydraulic model properly detected negative pressures in the WDSs. The Lagrangian-based contaminant transport model also reduced simulation time and avoided numerical dissipation and dispersion errors. The proposed WQM provided a tool to evaluate the vulnerably of a WDS in contaminant intrusion events. 7.2 Contributions This dissertation proposed a novel approach to study the contaminant intrusion to a WDS. Moreover, the key contributions of this dissertation are: • Proposing a fuzzy-rule model to predict intrusion by integrating information related to contaminant source, driving force, and pathways in a cast-iron pipe. Using the proposed model, susceptible locations of a WDS to contaminant intrusion can be detected and prioritized. • Developing an ingress model incorporating the effects of the surrounding soil (i.e. soil porosity and particles diameters) to estimate intrusion flow rate. The proposed model was first ingress model considering the surrounding soil effects on intrusion flow rate. • Creating a WQM coupling transient hydraulic model with a Lagrangian contaminant transport model to simulate contaminant intrusion propagation. The proposed WQM was first Lagranigian-Eulerian transient WQM to simulate intrusion propagation. 7.3 Limitations The models developed herein to identify potential of contaminant intrusion, estimate intrusion flow rates, and simulate contaminant propagation have few limitations. The proposed model to identify susceptible locations of a WDS to contaminant intrusion used internal/external corrosion models. The application of this model is limited to the cast-iron WDSs. The proposed ingress model in this research neglects the interaction 101 between leak points and low/negative pressures. However, the magnitude of low/negative pressures in WDSs may reduce due to the presence of leak point, or occurrence of intrusion. The proposed WQM is able to simulate contaminant propagation in intrusion events; however, the model cannot be used to study the fate of a continuous contaminant release in WDSs. Also, the proposed transient WQM does not take into account mass transport due to diffusion/dispersion, which is important in lowflow or dead-end pipes. 7.4 Recommendations for Future Research • In this research, uncertainty analysis techniques were studied, while the uncertain parameters were assumed independent. However, in risk assessment models, there often exist dependencies among the parameters that can have profound impacts on risk being assessed. This needs to be studied in further details in future. • A more comprehensive model is required to estimate potential of contaminant intrusion in a WDS with visco-elastic pipes. • A more realistic ingress model is required to take into account the interaction between leak points and low/negative pressures. A comprehensive investigation is also required to estimate the cross sectional area and shape of leak points in real-life WDSs. • Although the proposed transient WQM was able to track the fate of a few contaminant parcels along a system, it has two obvious modelling shortcomings. First, the model cannot study the fate of a continuous contaminant release. Second, the proposed WQM does not take into account mass transport due to diffusion/dispersion. These need to be studied in more details in future. • A more comprehensive studies need to be conducted to detect contributions of different sources of uncertainties on contaminant intrusion including pipes frictions factors, external contaminant concentration, duration of transient events, and so on. 102 References Andrade, M., Rojano, F., Romero-Gomez, P., and Choi, C. (2010). "Integrated Water Quality Modeling of Water Distribution Systems." Proc. 10th Annual Water Distribution System Analysis Symp, 293-299. Arora, J. S. (2004). Introduction to optimum design. Academic Press. 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(1-4), 204. Basha, H., and Malaeb, L. (2007). "Eulerian–Lagrangian Method for Constituent Transport in Water Distribution Networks." J. Hydraul. Eng., 133 1155. Besner, M. C., Prévost, M., and Regli, S. (2011). "Assessing the public health risk of microbial intrusion events in distribution systems: Conceptual model, available data, and challenges." Water Res., 45(3), 961-979. Besner, M. C., Broséus, R., Lavoie, J., Giovanni, G. D., Payment, P., and Prévost, M. (2009). "Pressure monitoring and characterization of external sources of contamination at the site of the Payment drinking water epidemiological studies." Environ.Sci.Technol., 44(1), 269-277. Blackburn, B. G., Craun, G. F., Yoder, J. S., Hill, V., Calderon, R. L., Chen, N., Lee, S. H., Levy, D. A., and Beach, M. J. (2004). "Surveillance for waterborne-disease outbreaks associated with drinking water-United States, 2001-2002." MMWR Surveill Summ, 53(8), 23-45. Borchardt M., Spencer S., Kieke B., Lambertini E., Loge F. (2009). "Do Water Distribution Systems Contribute to Acute Gastrointestinal Illness Incidence?" Proceedings of AWWA-Water Quality Technology Conference, Seattle, WA, USA, November 15–18. Boulos, P. F., Karney, B. W., Wood, D. J., and Lingireddy, S. (2005). "Hydraulic transient guidelines for protecting water distribution systems." J.Am.Water Works Assoc., 97(4), 111-124. Boulos, P. F., Altman, T., Jarrige, P. A., and Collevati, F. (1994). "An event-driven method for modelling contaminant propagation in water networks", Appl. Math. Model., 18(2), 84-92. 103 Boxall, J., Skipworth, P., and Saul, A. (2003). "Aggressive flushing for discolouration event mitigation in water distribution networks." Water Supply, 3(1-2), 179-186. Chaudhry, M. H. (1987). "Applied hydraulic transients." Second Edition.Van Nostrand Reinhold Co.New York.1987.521, . Collins, R., Beck, S., and Boxall, J. (2011). "Intrusion into water distribution systems through leaks and orifices: Initial experimental results." 11th Computing and Control in the Water Industry 2011. Collins, R., Besner, M. C., Beck, S., Karney, B., and Boxall, J. (2010). "Intrusion Modelling and the Effect of Ground Water Conditions." Proc. 12th Annual Water Distribution Systems Analysis Conference (WDSA 2010), Tucson (AZ). Corso, P. S., Kramer, M. H., Blair, K. A., Addiss, D. G., Davis, J. P., and Haddix, A. C. (2003). "Costs of Illness in the 1993 Waterborne Cryptosporidium Outbreak, Milwaukee, Wisconsin." Emerging Infectious Diseases, 9(4), 426. Deng, Y., Jiang, W., and Sadiq, R. (2010). "Modelling contaminant intrusion in water distribution networks: a new similarity-based DST method." Expert Syst.Appl.. 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(1-2), 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. Ebacher, G., Besner, M. C., Lavoie, J., Jung, B., Karney, B., and Prévost, M. (2011). "Transient Modeling of a Full-Scale Distribution System: Comparison with Field Data." J.Water Resour.Plann.Manage., 137 173. Elfeki, A. M. M. (2006). "Reducing concentration uncertainty using the coupled markov chain approach." Journal of Hydrology, 317(1-2), 1. Esteban Fernandez, J., del Rocio Fernandez, M., Vijande Diaz, R., and Tucho Navarro, R. (2003). "Abrasive wear analysis using factorial experiment design". Wear, 255(1-6), 38-43. Fang, K. (2006). Design and modeling for computer experiments. CRC Press. 104 Faybishenko, B. (2005). Dynamics of fluids and transport in fractured rock. Amer Geophysical Union. Fernandes, C., and Karney, B. (2004). "Modelling the advection equation under water hammer conditions", Urban Water Journal, 1(2), 97-112. Ferson, S., Nelsen, R. B., Hajagos, J., Berleant, D. J., Zhang, J., Tucker, W. T., et al. (2004). "Dependence in probabilistic modeling, dempster-shafer 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(3-4), 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(3-4), 328. Geldreich, E. E., Fox, K. R., Goodrich, J. A., Rice, E. W., Clark, R. M., and Swerdlow, D. L. (1992). "Searching for a water supply connection in the Cabool, Missouri disease outbreak of Escherichia coli 0157: H7." Water Res., 26(8), 1127-1137. Greyvenstein, B., van Zyl, J.E. (2007). "An experimental investigation into the pressure-leakage relationship of some failed water pipes." Journal of Water Supply: Research and Technology-Aqua, 56 (2), 117e124. Gullick, R. W., LeChevallier, M. W., Case, J., Wood, D. J., Funk, J. E., and Friedman, M. J. (2005). "Application of pressure monitoring and modelling to detect and minimize low pressure events in distribution systems." Journal of Water Supply: Research and Technology.AQUA, 54(2), 65-81. Gullick, R. W., LeChevallier, M. W., Svindland, R. C., and Friedman, M. J. (2004). "Occurrence of transient: Low and negative pressures in distribution systems." J.Am.Water Works Assoc., 96(11), 52-66. 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. Hamby, D. (1994). "A review of techniques for parameter sensitivity analysis of environmental models." Environ. Monit. Assess., 32(2), 135-154. 105 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(1-2), 89. Ho, C. K. (2008). "Solute Mixing Models for Water-Distribution Pipe Networks", J. Hydraul. Eng., 134(9), 1236-1244. Hooper S.M., Moe C.L., Uber J.G., Nilson K.A., (2006). "Assessment of Microbiological Water Quality after Low pressure Events in a Distribution System", 8th Annual Water Distribution Systems Analysis Symposium, Cincinnati, Ohio, USA, August 27–30. Huang, H., and Ayoub, J. (2008). "Applicability of the Forchheimer equation for nonDarcy flow in porous media." SPE Journal, 13(1), 112-122. Jang, J. S. R., Gulley, N., and MathWorks, I. (1998). "Fuzzy Logic Toolbox: For Use with MATLAB: User's Guide." Jafari, A., and Babadagli, T. (2008). "A sensitivity analysis for effective parameters on fracture network permeability." SPE Western Regional and Pacific Section AAPG Joint Meeting. Jung, B. S., Boulos, P. F., Wood, D. J., and Bros, C. M. (2009). "A Lagrangian wave characteristic method for simulating transient water column separation", J.Am.Water Works Assoc., 101(6), 64-73. Kaplan, J. E., Goodman, R. A., Schonberger, L. B., Lippy, E. C., and Gary, G. W. (1982). "Gastroenteritis due to Norwalk virus: an outbreak associated with a municipal water system." J.Infect.Dis., 146(2), 190-197. 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. Karney, B. W., and Ghidaoui, M. S. (1997). "Flexible discretization algorithm for fixedgrid MOC in pipelines", J. Hydraul. Eng., 123(11), 1004-1011. Karim M.R., Abbaszadegan M., LeChevallier M., (2003). "Potential for Pathogen Intrusion during Pressure Transients", J.Am.Water Works Assoc., 95(5), 134– 146. Khajeh, E., and Maijer, D. M. (2011). "Permeability of dual structured hypoeutectic aluminum alloys." Acta Materialia. 106 Khajeh, E., and Maijer, D. M. (2010). "Physical and numerical characterization of the near-eutectic permeability of aluminum-copper alloys." Acta Materialia, 58(19), 6334-6344. Kiene, L., Lu, W., and Levi, Y. (1998). "Relative importance of the phenomena responsible for chlorine decay in drinking water distribution systems." Water Science and Technology, 38(6), 219-227. Kirmeyer G.J., Friedman M., Martel K., Howie D., LeChevallier M., Abbaszadegan M., Karim M., Funk J., Harbour J., (2001). "Pathogen Intrusion into Distribution System", AwwaRF, Denver, USA. Klotz, D. (1980). "Dispersivity and velocity relationship from laboratory and field experiments." Journal of Hydrology, 45, 169-184. Larock, B. E., Jeppson, R. W., and Watters, G. Z. (2000). “Hydraulics of Pipeline Systems”, CRC Press LLC., 537p. Lambert, A., (2001). "What do we know about pressure-leakage relationships in distribution systems." In: Proceedings of IWA System approach to leakage control and water distribution systems management, Brno, Czech Republic. LeChevallier, M. W.,Xu, M., and Yang, J., (2011). "Managing Distribution System Low Transient Pressures for Water Quality." Journal-American Water Works Association, [Project #4152]. LeChevallier, M. W., Gullick, R. W., Karim, M. R., Friedman, M., and Funk, J. E., (2003). "The potential for health risks from intrusion of contaminants into the distribution system from pressure transients." Journal of Water and Health, 1(1), 3-14. Li, J. (2003). "Development of an inexact environmental modeling system for the management of petroleum-contaminated sites." PhD. Thesis, University of Regina. Lindley T.R., (2001). "A Framework to Protect Water Distribution Systems against Potential Intrusions.", M.A.Sc. Thesis, University of Cincinnati, Cincinnati, 102p. Macdonald, I., El-Sayed, M., Mow, K., and Dullien, F. (1979). "Flow through porous media-the Ergun equation revisited." Industrial & Engineering Chemistry Fundamentals, 18(3), 199-208. 107 MacPhee M.J., (2005). "Distribution System Water Quality Challenges in the 21st Century: A Strategic Guide.", AWWA, USA, 196p. Mamdani, E. H., (1977). "Application of fuzzy logic to approximate reasoning using linguistic synthesis." IEEE Trans.Comput., 1182-1191. Mannina, G., and Viviani, G. (2010). "An urban drainage stormwater quality model: Model development and uncertainty quantification." Journal of Hydrology, 381(3-4), 248-265. Mansour-Rezaei, S., and Naser, G. (2012). "Contaminant Intrusion in Water Distribution Systems: An Ingress Model." World Environmental and Water Resources Congress 2012@ sCrossing Boundaries, ASCE, 3002-3010. Masky, S. (2004). Modelling uncertainty in flood forecasting systems, Aa Balkema. McInnis, D. (2004). "A relative-risk framework for evaluating transient pathogen intrusion in distribution systems." Urban Water Journal, 1(2), 113-127. Miragliotta, G. (2010). "The power of dimensional analysis in production systems design." Int J Prod Econ. Misic, T., Najdanovic-Lukic, M., and Nesic, L. (2010). "Dimensional analysis in physics and the Buckingham theorem." European Journal of Physics, 31 893. 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. Montgomery, D. C. (2008). Design and analysis of experiments. John Wiley & Sons Inc. Motamedi, D., and Milani, A. S. (2010). "Improving damage resistance of a composite pole using a computer experiment strategy." Proceedings of the ICE- Engineering and Computational Mechanics, 164(3), 171-181. Munavalli, G., and Kumar, M. (2004). "Modified Lagrangian method for modeling water quality in distribution systems." Water Res., 38(13), 2973-2988. Nair, R., Sunny, F., & Manikandan, S. (2009). "Modelling of decay chain transport in groundwater from uranium tailings ponds." Applied Mathematical Modelling. Naser, G., and Karney, B. W. (2008). "A transient 2-D water quality model for pipeline systems." Journal of Hydraulic Research, 46(4), 516-525. 108 Naser, G., and Karney, B. (2007). "A 2-D transient multicomponent simulation model: Application to pipe wall corrosion." Journal of Hydro-Environment Research, 1(1), 56-69. NRC, (2006). "Drinking Water Distribution Systems: Assessing and Reducing Risks", National Research Council, National Academy Press, Washington DC, USA, 404p. NRC, (2004). "Monitoring Water Quality in the Distribution System", National Guide to Sustainable Municipal Infrastructure, Version No. 1.0, 45p. Payment, P., Siemiatycki, J., Richardson, L., Renaud, G., Franco, E., and Prevost, M. (1997). "A prospective epidemiological study of gastrointestinal health effects due to the consumption of drinking water." Int.J.Environ.Health Res., 7(1), 5-31. Percival, S., and Walker, J. (1999). "Potable water and biofilms: a review of the public health implications." Biofouling, 14(2), 99-115. Rajani, B., and Makar, J., (2000). "A methodology to estimate remaining service life of grey cast iron water mains." Canadian Journal of Civil Engineering, 27(6), 1259-1272. Rajani, B., McDonald, S., and Felio, G. (1995). "Water mains break data on different pipe materials for 1992 and 1993." Report no.A-7019.1, National Research Council of Canada, Ottawa. Rossman, L. A., Center for Environmental Research Information, and National Risk Management Research Laboratory. (2000). "EPANET 2: users manual." Rossman, L. A., and Boulos, P. F. (1996). "Numerical methods for modeling water quality in distribution systems: A comparison", J. Water Resour. Plann. Manage., 122 137. Ryan BF, Joiner BL. (2001). Minitabt handbook 4th ed. Canada: Duxbury Thomson Learning. Sadiq, R., Kleiner, Y., and Rajani, B. (2009). "Proof-of-concept model to predict water quality changes in distribution pipe networks (Q-WARP)". Denver, CO: Water Research Foundation. Sadiq, R., Saint-Martin, E., and Kleiner, Y. (2008). "Predicting risk of water quality failures in distribution networks under uncertainties using fault-tree analysis." Urban Water, 5(4), 287-304. 109 Sadiq, R., Kleiner, Y., and Rajani, B., (2007). "Estimating Risk of Contaminant Intrusion in Distribution Networks Using Fuzzy Rule-Based Modeling." Computational Models of Risks to Infrastructure, 6 318. Sadiq, R., Kleiner, Y., and Rajani, B. (2006). "Estimating risk of contaminant intrusion in water distribution networks using Dempster–Shafer theory of evidence." Civ.Eng.Environ.Syst., 23(3), 129-141. Sadiq, R., Rajani, B., and Kleiner, Y., (2004). "Probabilistic risk analysis of corrosion associated failures in cast iron water mains." Reliab.Eng.Syst.Saf., 86(1), 1-10. Sander, P., Bergback, B., & Oberg, T. (2006). "Uncertain numbers and uncertainty in the selection of input distributions - consequences for a probabilistic risk assessment of contaminated land." Risk Analysis, 26(5), 1363. Sentz, K., and Ferson, S. (2002). "Combination of evidence in Dempster-Shafer theory." Report no.SAND2002, 835. 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(1-2), 30. Schuster, C. J., Ellis, A. G., Robertson, W. J., Charron, D. F., Aramini, J. J., Marshall, B. J., and Medeiros, D. T. (2005). "Infectious disease outbreaks related to drinking water in Canada, 1974-2001." Canadian Journal of Public Health.Revue Canadienne De Sante Publique, 96(4), 254. Simonovic, S. P. (1997). "Risk in sustainable water resources management." IAHS- AISH Publication (240) 3. Skraber, S., Schijven, J., Gantzer, C., and de Roda Husman, A. (2005). "Pathogenic viruses in drinking-water biofilms: a public health risk?" Biofilms, 2(2), 105117. Streeter, V. L., Wylie, E. B., and Bedford, K. W. (1998). Fluid mechanics, WCB. Systèmes, D. (2010). "Abaqus 6.10: Analysis User's Manual." Providence, RI: Dassault Systèmes Simulia Corp. Thomason, J., Wang, L., (2009). "Sate of technology review report on condition assessment of ferrous water transmission and distribution systems." EPA/600/R- 09/055. 110 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. Tzatchkov, V. G., Aldama, A. A., and Arreguin, F. I. (2002). "Advection-dispersionreaction modeling in water distribution networks", J. Water Resour. Plann. Manage., 128 334. Vairavamoorthy, K., Yan, J., and Gorantiwar, S. (2007). "Modelling the risk of contaminant intrusion in water mains." Water Management, 160(2), 123-132. Vairavamoorthy, K., Yan, J., Galgale, H. M., and Gorantiwar, S. D. (2007). "IRAWDS: A GIS-based risk analysis tool for water distribution systems." Environmental Modelling & Software, 22(7), 951-965. Vafai, K. (2005). Handbook of porous media. CRC. 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. Van Zyl, J., and Clayton, C. (2007). "The effect of pressure on leakage in water distribution systems." Water Management, 160(2), 109-114. 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(3-4), 283. Willems, P. (2008). "Quantification and relative comparison of different types of uncertainties in sewer water quality modeling." Water Res., 42(13), 3539-3551. Woo, D. M., and Vicente, K. J., (2003). “Sociotechnical Systems, Risk Management, and Public Health: Comparing the North Battleford and Walkerton Outbreaks”, Reliability Engineering and System Safety, 80(3), 253–269. Wood, D. J., Lingireddy, S., Boulos, P. F., Karney, B. W., and McPherson, D. L. (2005). "Numerical methods for modeling transient flow", J.Am.Water Works Assoc., 97(7), 104-115. 111 Xin, K., Tao, T., Wang, Y., and Liu, S. (2012). "Hazard and vulnerability evaluation of water distribution system in cases of contamination intrusion accidents", Frontiers of Environmental Science & Engineering, 1-10. Yager, R. R. (1986). "Arithmetic and other operations on Dempster-Shafer structures." International Journal of Man-Machine Studies, 25(4), 357-366. Yaminighaeshi, H. (2009). "Probability of failure analysis and condition assessment of cast iron pipes due to internal and external corrosion in water distribution systems." Ph.D., Thesis, The University of British Columbia. Yamini, H., and Lence, B. J., (2010). "Probability of Failure Analysis due to Internal Corrosion in Cast-Iron Pipes." J Infrastruct Syst, 16 73. 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(1-2), 1. Yuan, C., & Druzdzel, M. J. (2006). "Importance sampling algorithms for bayesian networks: Principles and performance." Mathematical and Computer Modelling, 43(9-10), 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(1-2), 73. Zimmermann, H. J. (2001). "Fuzzy set theory--and its applications." Kluwer Academic Pub. 112
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Contaminant intrusion in water distribution systems...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Contaminant intrusion in water distribution systems : advanced modelling approaches Mansour Rezaei Fumani, Saheb 2013
pdf
Page Metadata
Item Metadata
Title | Contaminant intrusion in water distribution systems : advanced modelling approaches |
Creator |
Mansour Rezaei Fumani, Saheb |
Publisher | University of British Columbia |
Date Issued | 2013 |
Description | Since exposure to contaminants may have direct adverse impacts on public health, contaminant intrusion has been recognized as one of the top priority in drinking water supply research. Three components must exist to cause contaminant intrusion into a water distribution system. These include the availability of source(s) of contaminant(s) around a water distribution system, the existence of driving forces (low/negative pressure) to make a contaminant enter into a water distribution system, and the presence of pathway(s) through which contaminant(s) intrude into a water distribution system (WDS). Exposure assessment is the most challenging part as location of contaminant intrusion, rate of intrusion, and the fate of contaminants within WDS need to be estimated accurately. In this dissertation, first, common uncertainty analysis techniques are discussed in the context of conservativeness, execution time, ease of formulation, and complexity. Second, a fuzzyrule based model has been developed to identify contaminant intrusion potential in a WDS. The potential of contaminant intrusion has been determined by integrating the potentials for contaminant sources existence, driving forces, and pathways. Third, a novel ingress model has been developed for more realistic estimation of intrusion rate by taking into account the effects of surrounding soil on intrusion rate. Coupled with an Eulerian-based transient hydraulic model, a Lagrangian transient water quality model is developed to predict the fate of the contaminant throughout a WDS. The proposed models are applied to case studies available in the literature to investigate the applicability of the models. The proposed models enhance the reliability and safety of WDSs by improving the prediction ability of the existing modelling tools. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2013-03-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0073626 |
URI | http://hdl.handle.net/2429/44042 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Engineering, School of (Okanagan) |
Degree Grantor | University of British Columbia |
GraduationDate | 2013-05 |
Campus |
UBCO |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2013_spring_mansour rezaei fumani_saheb.pdf [ 2.06MB ]
- Metadata
- JSON: 24-1.0073626.json
- JSON-LD: 24-1.0073626-ld.json
- RDF/XML (Pretty): 24-1.0073626-rdf.xml
- RDF/JSON: 24-1.0073626-rdf.json
- Turtle: 24-1.0073626-turtle.txt
- N-Triples: 24-1.0073626-rdf-ntriples.txt
- Original Record: 24-1.0073626-source.json
- Full Text
- 24-1.0073626-fulltext.txt
- Citation
- 24-1.0073626.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0073626/manifest