GIS BASED STOCHASTIC MODELING OF GROUNDWATER SYSTEMS USING MONTE CARLO SIMULATION By Dipa Dey A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Civil Engineering 2010 ABSTRACT GIS BASED STOCHASTIC MODELING OF GROUNDWATER SYSTEMS UISNG MONTE CARLO SIMULATION By Dipa Dey Stochastic modeling of subsurface flow and transport has become a subject of wide interest and intensive research for last few decades and results evolution of many stochastic theories. These theories, however, have had relatively little impact on practical groundwater modeling. In a recent forum on stochastic subsurface hydrology: from theory to application, a number of leading experts in stochastic modeling stress that data limitation, the assumptions of linearization, stationarity, Gaussianity, and excessive computations required are the bottlenecks in practical stochastic modeling. These bottlenecks must be removed or substantially relaxed before stochastic modeling methods can be routinely applied in practice. Motivated by these critical assessments, this research addresses the issue of Gaussianity and issues of data limitations in stochastic modeling. The issue of Gaussianity is addressed by Monte Carlo Simulation (MCS). Data limitations issues are addressed using new source of Geographic Information Systems (GIS) database. My first application considers a comprehensive study of synthetic scenarios to investigate the probabilistic structure of basic hydrogeological variables such as hydraulic head, groundwater velocity, concentration, seepage flux and solute flux. Results indicate that the statistical structure of groundwater systems in general non-Gaussian, nonstationary and anisotropic. Some critical state variables are extremely complex, with the probability distribution varying rapidly with locations and directions even for very weak heterogeneity. This study concludes that we can not always use variance as a good measure of uncertainty. Actual probability distribution from more accurate and generalized method such as MCS is a better way to characterize the structure of hydrogeological variables. This research represents the first systematic analysis of the probabilistic structure of basic hydrogeological variables and findings from this work have significant implications on theoretical and practical stochastic subsurface hydrology. My second application involves probabilistic delineation of well capture zones. In this paper, we explore the use of a recently-developed statewide GIS database in Michigan. We are particularly interested in exploring if the relatively crude, but detailed datasets can be used to characterize aquifer heterogeneity with sufficient details to enable practical stochastic modeling. We consider three approaches such as deterministic, stochastic macrodispersion, stochastic Monte Carlo, to delineate well capture zones, each representing a different way to conceptualize aquifer heterogeneity. Results show deterministic approach that accounts only for the effect of large-scale trend leads to capture zones that are significantly smaller than its stochastic counterpart. Stochastic Monte Carlo approach that models large-scale trends deterministically and small-scale heterogeneity as random field provides a probability map of well capture zone which is useful for risk-based decision making processes. Stochastic macrodispersion approach that models large-scale trends deterministically and small-scale heterogeneity as effective macrodispersion, provides a computationally efficient alternative to delineate well capture zones. A probability map of well capture zone has important implications for environmental policy on source water protection, risk management and sampling design. DEDICATION TO MY FAMILY iv ACKNOWLEDGMENTS I would like to express my deep sense of gratitude to my advisor Dr. Shu-Guang Li, for his invaluable support, patience and encouragement throughout my studies in Michigan State University. I wish to extend my appreciation to Dr. Mantha Phanikumar, Dr. Alexandra Kravchenko, Dr. Roger Wallace, and Dr. Howard Reeves for serving in my dissertation committee and for valuable inputs in this study. I am thankful to Dr. Huasheng Liao for his invaluable suggestions during my research work. I am also thankful to my research team members, Prasanna Sampath, Hassan Abbas, and Mehmet Otzan. Last, but not the least, I would like to thank my husband Mr. Rajen Ghosh for his understanding and love during the past few years. His support and encouragement was in the end what made this dissertation possible. My family members receive my deepest gratitude and love for their support during my studies. v TABLE OF CONTENTS LIST OF TABLES……………………………………………………………………… vii LIST OF FIGURES…………………………………………………………………...... viii PAPER1: PROBABILISTIC STRUCTURE OF HYDROGEOLOGIC VARIABLES……………………………………………………………..... 1.1 ABSTRACT……………………………………………………………….. 1.2 INTRODUCTION…………………………………………………………. 1.3 OBJECTIVE……………………………………………………………….. 1.4 APPROACH……………………………………………………………….. 1.5 PROBLEM DESIGN………………………………………………………. 1.6 RESULTS AND DISCUSSION…………………………………………… 1.7 CONCLUSIONS…………………………………………………………… APPENDIX………………………………………………………………… REFERENCES…………………………………………………………….. 1 1 3 5 7 8 14 38 41 65 PAPER2: PROBABILISTIC WELL CAPTURE ZONE DELINEATION USING A STATEWIDE GROUNDWATER DATABASE SYSTEM………………. 2.1 ABSTRACT……………………………………………………………….. 2.2 MOTIVATION AND OBJECTIVE……………………………………….. 2.3 DATA INTENSIVE STATEWIDE MODLEING SYSTEM……………… 2.4 RESEARCH APPROACHES……………………………………………… 2.5 APPLICATION EXAMPLES……………………………………………... 2.6 RESULST AND DISCSSION……………………………………………... 2.7 SUMMARY AND CONCLUSIONS……………………………………… REFERENCES…………………………………………………………….. 70 70 72 76 78 82 84 108 111 vi LIST OF TABLES Table 1.1 Parameter Definitions for four cases……………………………………………10 Table 1.2 Statistical parameters (Skewness and Kurtosis) of Head, X-Velocity and Yvelocity for low and high lnK Variances……………………………………….41 Table 1.3 Statistical parameters (Skewness and Kurtosis) of Head, X-Velocity and YVelocity for Nonstationary with pumping wells………………………………..42 Table 1.4 Statistical parameters (Skewness and Kurtosis) of Head, for Nonstationary with pumping wells having large scale heterogeneity………………………………..42 Table 1.5 Statistical parameters of Head in the sub model domain for nonstationary case with pumping well………………………………………………………………42 Table 1.6 Statistical parameters (Skewness and Kurtosis) of Concentration For low variance of lnK and for nonstationary case with pumping wells………43 Table 1.7 Statistical parameters (Skewness and Kurtosis)of seepage and solute fluxes For low and high lnK variances………………………………………………...43 Table 2.1 Geostatistical Parameters for three sites………………………………………...84 Table 2.2 Parameters definitions for three sites…………………………………………...95 vii LIST OF FIGURES Figure 1.1 Model domain showing locations of monitoring wells, polylines (each polyline is a set of two polylines as shown in extreme right of the figure), pumping wells, plume source and different conductivity zones, Zone 1 (lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day). Zone 1 and Zone 2 are shown in shaded areas For interpretation of the references to color in this and all other figures, the reader is referred to the electronic version of this thesis (or dissertation).……………………. 9 Figure 1.2 Probability distributions of hydraulic head for stationary cases (Top: lnK variance =0.1, Bottom: lnK variance =2.0). Black circles are monitoring wells and thin solid lines are head contours in the model domain…………. 16 Figure 1.3 Probability distributions of hydraulic head for Nonstationary case without pumping wells (Zone 1(lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day)). Black circles are monitoring wells and thin solid lines are head contours in the model domain………………………………. 17 Figure 1.4a Probability distributions of hydraulic head for Nonstationary case with pumping wells (Zone 1 (lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day)). Black circles are monitoring wells and thin solid lines are head contours in the model domain………………………………. 18 Figure 1.4b Probability distributions of hydraulic head observed close to the pumping well. Black circles are monitoring wells and thin solid lines are the head contours in the model domain……………………………………………… 19 Figure 1.4c Probability distributions of hydraulic head for Nonstationary case with pumping wells having large correlation scale (Zone 1 (lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day)). Black circles are monitoring wells and thin solid lines are head contours in the model domain……………………………………………………………………… 20 Figure 1.5a Probability distributions of X-velocity for stationary cases (Top: lnK variance =0.1, Bottom: lnK variance =2.0). Black circles are monitoring wells and thin solid lines are head contours in the model domain…………. 23 Figure 1.5b Probability distributions of Y-velocity for stationary cases (Top: lnK variance =0.1. Bottom: lnK variance =2.0). Black circles are monitoring wells and thin solid lines are head contours in the model domain…………. 24 viii Figure 1.6a Probability distributions of X-Velocity for Nonstationary case without pumping wells (Zone 1 (lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day)). Black circles are monitoring wells and thin solid lines are head contours in the model domain………………………………. 25 Figure 1.6b Probability distributions of Y-Velocity for Nonstationary case without pumping wells (Zone 1 (lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day)). Black circles are monitoring wells and thin solid lines are head contours in the model domain………………………………. 26 Figure 1.7a Probability distributions of X-Velocity for Nonstationary case with pumping wells (Zone 1 (lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day)). Black circles are monitoring wells and thin solid lines are head contours in the model domain………………………………. 27 Figure 1.7b Probability distributions of Y-Velocity for Nonstationary case with pumping wells (Zone 1 (lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day)). Black circles are monitoring wells and thin solid lines are head contours in the model domain………………………………. 28 Figure 1.8 Probability distributions of concentration for Stationary case. Black circles are monitoring wells and the color contours are the concentration………… 30 Figure 1.9 Probability distributions of concentration for Nonstationary case (Zone 1 (lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day)). Black circles are monitoring wells and the color contours are the concentration…………………………………………………………… 31 Figure 1.10 Probability distributions of seepage flux for stationary cases at each polylines. Thick solid lines represent the length of the polylines and thin solid lines are the head contours in the model domain. Top: lnK variance =0.1, Bottom: lnK variance =2.0…………………………………………… 33 Figure 1.11 Probability distributions of seepage flux at each polylines for nonstationary case with pumping wells (Zone 1 (lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day)). Thick solid lines represent the length of the polylines and thin solid lines are the head contours in the model domain……………………………………………… ix 34 Figure 1.12 Figure 1.13 Figure 1.14 Figure 1.15 Figure 1.16 Probability distributions of solute flux at each polylines for Stationary case. Thick solid lines represent the length of the polylines and the color contours are the concentration……………………………………………... 36 Probability distributions of solute flux at each polylines for Stationary case. Thick solid lines represent the length of the polylines and the color contours are the concentration……………………………………………... 37 Q-Q plot of Hydraulic head data versus normal distribution for low variance and high variance cases. MW3: located near the west boundary, MW28: Located near the east boundary, MW13: located approximately middle of the model domain……………………………………………….. 44 Q-Q plot of X-Velocity data versus normal distribution for low and high lnK variance………………………………………………………………... 50 Q-Q plot of Y-Velocity data versus normal distribution for low and high lnK variance………………………………………………………………... 52 Figure 1.17 Q-Q plot of concentration data versus normal distribution………………… 54 Figure 1.18 Q-Q plot of seepage flux data versus normal distribution for low and high lnK variance………………………………………………………………... 56 Q-Q plot of solute flux data versus normal distribution for low and high lnK variance………………………………………………………………... 60 Figure 2.1 Approaches used to delineate well capture zones………………………….. 81 Figure 2.2 Sites located in Michigan showing GIS features such as wells, streams, lakes, city, Type 1 wells and underground storage tanks…………………... 83 Figure 2.3 Variogram modeling at regional scale for log hydraulic conductivity for three sites…………………………………………………………………… 85 Figure 2.4 Variogram modeling for detrended log hydraulic conductivity for three sites…………………………………………………………………………. 86 Figure 2.5 Kriged hydraulic conductivity at regional scale. The outer boundary is the extent of regional model and inner boundary is the location of local model. Figure 1.19 89 Figure 2.6 A realization of residual hydraulic conductivity (Bottom) and large-scale hydraulic conductivity (Top) for the local scale model area……………...... 91 Figure 2.7 Site1: Head field for regional model (bottom) and local model (top)……… 97 x Figure 2.8 Site2: Head field for regional model (bottom) and local model (top)……… 98 Figure 2.9 Site3: Head field for regional model (bottom) and local model (top)……… 99 Figure 2.10 Comparison of predicted head and observed head (Static water level) for three sites from water wells in the regional model. Dots: Static water level from well data, Triangles: Mean static water level from well data………… 100 Figure 2.11 Site 1: 10 years capture zone from three approaches. …………………….. 104 Figure 2.12 Site 2: 10 years capture zone from three approaches. …………………….. 105 Figure 2.13 Site 3: 10 years capture zone from three approaches. …………………….. 106 xi PAPER 1 PROBABILISTIC STRUCTURE OF HYDROGEOLOGIC VARIABLES Dipa Dey and Huasheng Liao Department of Civil and Environmental Engineering Michigan State University, East Lansing, Michigan 48824 1.1 ABSTRACT Most stochastic groundwater analyses rely on the critical assumption of Gaussianity for hydrogeological variables involved, although little is known systematically on their probabilistic structure. The conventional wisdom is – as long as the input source variability is not too strong and Gaussian, most dependent hydrogeological variables will be approximately normally distributed and a first-order method can be used to compute the first two statistical moments which can then be used to characterize the probabilistic distributions. This assumption is essentially the foundation of stochastic perturbation theory and the recent stochastic revolution that produced a huge body of literature on perturbation methods in the past three decades. Many researchers, however, question the validity of this assumption, stressing that the reliance on Gaussianity is one of the reasons stochastic methods are so difficult to apply in practice. In this study, we investigate systematically the probabilistic structure of basic hydrogeological variables by performing an integrated analysis of groundwater flow and solute transport in heterogeneous media. In particular, we investigate the probabilistic distribution of hydraulic head, velocity, solute concentration, seepage flux, and solute flux, assuming hydraulic conductivity to be the source of uncertainty. We also investigate the effect of spatial integration on the probabilistic structure of seepage and solute fluxes. We consider both statistically uniform 1 and nonuniform flows for a range of situations, including weakly and strongly heterogeneous media, stationary and nonstationary media, and with and without aquifer stresses. We use the Monte Carlo approach to solve the stochastic groundwater flow and transport equations and make use of a recently-developed stochastic modeling software system. This software system is based on a new computing paradigm that enables real-time, parallel visual analysis and modeling of large numbers of realizations without storing them in memory or on disk, eliminating an implementation bottleneck in stochastic modeling. We observe and analyze the simulated probability distributions of interested variables throughout the computational domain at representative locations and of spatially integrated fluxes across transects of different lengths. We assess the Gaussianity of these probability distributions both visually and using the Q-Q plot technique. We also assess the feasibility to fit the simulated probabilities with empirical distribution under general conditions. The following key findings emerge from this research: • The statistical structure of a groundwater system under general, nonstationary condition is complex and cannot be characterized by a few simple statistical moments. Even when lnK is assumed to be Gaussian, the dependent hydrogeological variables can be strongly non-Gaussian, nonstationary, and for those that are directional, anisotropic. • The closest to Gaussian, among all variables simulated, is the probability distribution of hydraulic head, when the scale of heterogeneity is much smaller than the domain size (e.g., 1/33 of the problem scale). • The probability distribution of hydraulic head in response to larger scale heterogeneity (e.g., greater 1/10 of the problem scale) is often non-Gaussian and nonstationary and cannot be accurately fitted with a single functional distribution. 2 • The probability distribution of groundwater velocity is direction-dependent. The probabilistic structure of longitudinal velocity is strongly skewed. • The probability distribution of transverse velocity is approximately Gaussian or symmetrically distributed, even in the presence of strong heterogeneity. • Seepage flux across a transect – representing integration of strongly skewed velocity, is approximately Gaussian, even in the presence of strong heterogeneity. • Probability distribution of solute concentration is far from Gaussian and strongly nonstationary when the source size is relatively small (e.g., less than 10 lnK correlation scales), even for aquifers with very weak heterogeneity. The complex concentration probability distribution cannot be fitted with a single functional distribution. • Solute flux across a transect, unlike its seepage flux counterpart, is strongly non-Gaussian and location-dependent. This research represents the first comprehensive analysis of the probabilistic structure of basic hydrogeological variables and findings from this work have significant implications on theoretical and practical stochastic subsurface hydrology. 1.2 INTRODUCTION Most stochastic groundwater analyses rely on the critical assumption of Gaussianity. The conventional wisdom is – as long as the source variability is not too strong and is Gaussian, most dependent hydrogeological variables would be approximately normally distributed. This is essentially the foundation of stochastic perturbation theory and the huge body of literature on 3 stochastic perturbation methods for flow and transport developed in the past three decades [Dagan, 1989; Gelhar, 1993; Cushman, 1997; Neuman, 1997; Zhang, 2001; Rubin, 2003]. Many researchers, however, question the validity of this assumption, stressing that the reliance on Gaussianity is one of the major reasons why stochastic methods are so difficult to apply in practice [Bellin et al., 1992a; Zhang and Neuman, 1995a-d; Hamed et al., 1996; Saladin and Fiorotto, 1998; Guadagnini and Neuman, 1999a,b; Zhang, 2001; Lu et al., 2002; Dagan, 2002; Li et al., 2003; Winter et al., 2003; Zhang and Zhang, 2004; Ginn, 2004; Neuman, 2004]. In a recent forum on stochastic subsurface hydrology: from theory to application, a number of leading experts in stochastic modeling stress that data limitation, the assumptions of linearization, stationarity, Gaussianity, and excessive computations required are the bottlenecks in practical stochastic modeling [Zhang and Zhang, 2004; Molz, 2004; Molz et al., 2004; Ginn, 2004; Winter, 2004; Neuman, 2004; Rubin, 2004; Dagan, 2004] Experience from our own applications also shows that the observed variances of hydraulic head from linearized stochastic theories often cannot explain much of the observed variability in practical modeling. In fact, the observed variance in response to random, small-scale heterogeneity is almost negligible as compared to the uncertainty caused by errors in the representation of large scale variability, model conceptualization, and boundary conditions. Gelhar [1993] stressed that the most significant effect of small-scale heterogeneity is on groundwater velocity and solute transport [Gelhar et al., 1979; Kapoor and Kitanidis, 1997; Salandin and Fiorotto, 2000; Darvini and Salandin, 2006; Morales-Casique et al., 2006a, 2006b]. 4 In particular, Gelhar stressed that the observed variances from stochastic theories in groundwater velocity and solute concentration are much larger than that in hydraulic head. This strong sensitivity of solute transport to very small scale heterogeneity is in fact what attracted the recent explosion of research in stochastic subsurface hydrology. The significantly increased variance, however, is also what makes the application of stochastic perturbation theories problematic Nowak et al. [2008] shows hydraulic heads in the presence of large scale heterogeneity can be non-Gaussian and approximated as beta distribution for statistically uniform flow between two constant head boundaries. Engler et al. [2006] shows that velocity in randomly heterogeneous media is approximately lognormal, with the skewness increasing with the variance of log hydraulic conductivity. Fiorotto and Caroni [2002], Caroni and Fiorotto [2005], Bellin and Tonina [2007] and Schwede et al. [2008] have studied the statistical distribution of concentration in a random heterogeneous medium and concluded that the probability density function of concentration follows approximately a beta distribution. 1.3 OBJECTIVE In this study, we systematically investigate the probabilistic structure of basic hydrogeological variables. In particular, we investigate the probabilistic distributions of a complete set of basic dependent state variables - hydraulic head, velocity, solute concentration, and seepage flux and solute flux at a point and across a polyline. We assume log hydraulic conductivity is the only source of uncertainty and is normally distributed. The paper focuses on understanding how the nonlinear flow and transport equations transform the probability distributions of hydrogeological variables. This paper builds on several recent research that highlights the need to systematically 5 understand the probabilistic structure of hydrological variables [Engler et al., 2006; Fiorotto and Caroni, 2002; Caroni and Fiorotto, 2005; Bellin and Tonina, 2007; Nowak et al., 2008; Schwede et al., 2008]. Engler et al. [2006] stressed that with increasing variance of log hydraulic conductivity, the skewness and kurtosis of velocity components increases. This implies ignoring the actual shapes of the distributions may have severe consequences on the stochastic flow and transport theory. The paper investigates the following important questions: 1. What is the probabilistic structure of commonly encountered hydrogeological variables in realistically complex groundwater systems? 2. Can the predicted variances of hydrogeological variables be used as an effective measure of uncertainty in groundwater modeling? 3. How does the scale of heterogeneity impact the probability distribution? How do the scale of variability and complex sources and sinks interact to control the probability distribution? 4. Is the head probability distribution approximately Gaussian? Under what conditions it can be assumed to be Gaussian? 5. What is the probability distribution of groundwater velocity and solute concentration? How do their probability distributions vary in space in a nonstationary flow system? 6. Does the probability distribution of groundwater velocity vary with direction? 6 7. Is aggregated probability distribution of seepage flux across a polyline more Gaussian? How does the scale of integration affect the probability distributions of seepage flux? 8. What is the probability distribution of solute fluxes? 1.4 APPROACH To address these questions systematically, we perform a comprehensive Monte Carlo simulation exercise involving integrated flow and transport modeling. This simulation exercise is conceptually straightforward, but the results, as we will show further on, have significant implications and provide very useful insights for practical stochastic groundwater modeling. We consider a range of modeling scenarios in this exercise, including • weakly and strongly heterogeneous media; • stationary and nonstationary media; • with and without external stresses. We make use of a recently-developed modeling software system, Interactive Groundwater (IGW) [Li and Liu, 2003, 2006; Li et al., 2006] to solve the stochastic groundwater flow and transport equations. This software system differs from a traditional one in that it allows real-time, parallel analysis and visualization of model statistics and probabilities, enabling simulating large numbers of realizations for a range of modeling scenarios. The software also allows stochastic simulations without simultaneously storing them in memory or on disk, eliminating an implementation bottleneck in stochastic modeling. 7 We monitor, analyze, and compare the probability distributions of interested key state variables at representative locations throughout the simulation domain. We assess the Gaussianity of these probability distributions both visually and using the Q-Q plot technique [Chambers et al., 1983] We also assess the feasibility to fit the simulated probabilities with a general empirical distribution. In a Q-Q (Quantile-Quantile) plot, the quantiles of dependent variable from the Monte Carlo simulation are plotted against the quantiles of a standard normal distribution. A straight line implies that the distribution is a normal distribution. Any deviation from the straight line is the indication of deviation from the normality. 1.5 PROBLEM DESIGN We consider two dimensional steady state flows in confined aquifer. Figure 1.1 shows the model domain, boundary conditions, aquifer stresses, and monitoring network. We develop numerical models with no flow boundary conditions along north and south boundaries and constant heads along east and west boundaries. A constant head of 2 m is assigned to west boundary and 0 (zero) m to east boundary. The domain size is 1000 m x 750 m. A source of constant recharge of 0.001 m/day, two pumping wells, and statistically nonstationary conductivity zones are also introduced in some modeling scenarios. A contaminant plume is injected from the west side of the model. The plume source is continuous with concentration of 100 ppm (mg/L) spanning across twice the length of correlation scale of log hydraulic conductivity. 8 No Flow Boundary 1 MW1 MW 6 MW11 MW16 MW21 MW26 Constant Head = 0 m σ 2 = 0 .5 Y λ=2m K = 2 m/d Source Constant Head = 2m Zone 1 MW5 MW10 MW15 MW20 MW25 MW30 Zone 3 Zone 2 σ 2 = 1. 5 Y 2 σ Y = 3 .0 λ =10 m K = 20 m/d 0 2 λ = 20 m K = 200 m/d No Flow Boundary 250 500 m Pumping Wells σ2 = Y Monitoring Wells lnK variance Figure1.1: Model domain showing locations of monitoring wells, polylines (each polyline is a set of two polylines as shown in extreme right of the figure) , pumping wells, plume source and different conductivity zones, Zone 1 (lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day). Zone 1 and Zone 2 are shown in shaded areas. For interpretation of the references to color in this and all other figures, the reader is referred to the electronic version of this thesis (or dissertation). 9 To investigate the probabilistic structure of hydrogeological variables, we consider four different simulation scenarios, as shown in Table 1.1. Table 1.1: Parameter Definitions for four cases Case 1, Stationary flow, low lnK variance Case 2, Stationary flow, high lnK variance Case 3, nonstationary flow, nonstationary lnK variances Case 4, nonstationary flow, with pumping Domain length (m) 1000 x 750 1000 x 750 1000 x 750 1000 x 750 West boundary condition Constant Head Constant Head Constant Head Constant Head East boundary condition Constant Head Constant Head Constant Head Constant Head North boundary condition No Flow No Flow No Flow No Flow South boundary condition No Flow No Flow No Flow No Flow Global Recharge(m/day) 0.001 0.001 0.001 0.001 Number of wells 0 0 0 2 Geometric mean conductivity (KG) (m/day) 100 100 2,20,200 2,20,200 LnK correlation scales (λ)(m) 30 30 2, 10, 20 2, 10, 20 Ln K Variance 0.1 2 0.5, 1.5, 3 0.5, 1.5, 3 2 ( σY ) 10 The first two cases consider a stationary system with a variance of 0.1 (weakly heterogeneous) and 2.0 (strongly heterogeneous) for log hydraulic conductivity, respectively. The third and fourth cases involve nonstationary systems with three different hydraulic conductivity zones and pumping wells, respectively. The lnK variances used reflect the range of variability observed in real-world sites. The well known Cape Cod site in Massachusetts is relatively homogeneous and has lnK variance of 0.24. The popular Columbus research site in Mississippi is strongly heterogeneous and has lnK variance of 4.5 [Gelhar, 1993] Monitoring wells are distributed in a grid format of 5 rows and 6 columns (total of 30) to observe the probability distributions of hydraulic head, longitudinal velocity, transverse velocity and concentration (Figure 1.1). Groups of polylines at seven different locations perpendicular to the flow direction are used to observe the probability distributions of integrated seepage and solute fluxes. Within each group, polylines of two different lengths [see Figure 1.1, extreme right numbered as 1 and 2] are used to investigate the effect of scale of integration on the shape of probability distributions of seepage and solute fluxes. Monte Carlo Simulation: The MC simulation is carried out in three steps: Random Field Generation: This step creates realizations of random hydraulic conductivity. The random field of hydraulic conductivity is generated using a Fast Fourier Transform-based spectral algorithm. The exponential covariance model [Gelhar and Axness, 1983; Gelhar, 1986; 11 Dagan, 1994; Ni and Li, 2005, 2006] is used to describe the spatial structure of hydraulic conductivity. Numerical Solution: For a given realization, the governing equations become deterministic and are solved by central finite difference methods for groundwater flow and the modified methods of characteristics for solute transport. Discretization cell sizes are designed such that the scale of heterogeneity can be resolved [Li et al. 2003, 2004]. Statistical Postprocessing: Once the realizations are solved, the solution of flow and transport quantities such as hydraulic head, velocity, concentration and fluxes are processed to compute statistical moments and probability distributions. The statistical moments and probability distribution of hydraulic head at points of interest are computed by 1 N h (x) h(x) = s N ∑ n n =1 1 N [h (x) − h(x) ] 2 Var [h(x)]s = ∑ n s N n =1 (1.1) (1.2) Here N is the total number of realizations and hn (x) is the head at x from nth realization. The covariance of head at the points x and x1 is computed by 12 1 N h (x)h (x ) − h(x) h(x ) Cov[h(x),h(x1 )]s = n 1 1 s s N ∑ n n =1 (1.3) The probability distribution is computed by: p(h j ) = n j /N, j = 1,2, ..., M (1.4) where nj is the total number times the simulated head values fall into the interval [h j , h j + 1 ] , M is the number of interval the entire range of simulated heads are divided into, and M ∑ nj =N j =1 hj=hmin+ j (hmax-hmin)/M, The symbols hmin and hmax are respectively the minimum and maximum head values at the monitoring point. The accuracy of Monte Carlo simulation depends on the number of realizations. The number of realizations is selected to ensure that the observed probability distributions stabilize. The number of realizations used in this simulation exercise is generally in the range of 5000- 6000. 13 1.6 RESULTS AND DISCUSSION Modeling results are summarized in Figures 1.2 through 1.13. Selected Q-Q plots are shown in Figures 1.14 through 1.19. Additional information on the distribution statistics (skewness and kurtosis) are presented in Table 1.2 -1.7 in the Appendix. The results from Monte Carlo simulation clearly show that probability distribution of hydrogeological variables is, in general, non-Gaussian, nonstationary, and cannot be described by the first two statistical moments or a predefined distribution function based on these statistical moments. The probability distribution can vary structurally from variable to variable and location to location. The probability distributions of solute concentration and concentration fluxes are particularly complex, departing sharply from Gaussian distribution and varying strongly and rapidly with locations even for weakly heterogeneous conductivity fields. The following is a variable by variable discussion of the probability structure for all scenarios considered. Hydraulic Head Figures 1.2-1.4 show the observed probability distribution of hydraulic head for different cases. The results show that the probability distribution of hydraulic head, among all dependent variables, is the closest to Gaussian. The exceptions are in areas close to pumping well and boundary conditions. Figures 1.2 to Figure 1.4a show that probability distributions of hydraulic head when the scale of lnK heterogeneity is 1/33 of the simulation domain. The simulated probability distributions at 14 most locations are approximately Gaussian. The Q-Q plot in Figure 1.14 shows a straight line, indicating normality. The gaussianity of the head can be attributed to the fact that head at steady state depends on log conductivity, not conductivity itself, and the transformation between lnK fluctuation and head fluctuation is almost linear in areas where mean hydraulic gradient is small, especially when the scale of heterogeneity is also relatively small. For the high variance case, the probability distributions of hydraulic head at locations near the constant head boundaries (e.g., MW3, MW28) are slightly skewed towards the constant head values specified on the boundaries. The Q-Q plot of hydraulic head (see Appendix, Figure 1.14) at these locations shows a departure from the straight line, indicating local non-Gaussianity. Figure 1.3 shows the presence of the three conductivity zones alters the hydraulic gradient, head variance, and probability distribution, but not much in the shape of the probability distributions. The Q-Q plots for these probability distributions still show a straight line, indicating normality. Figures 1.4ab shows pumping can impact significantly the shape of the probability distribution of hydraulic head, especially in the proximity of the pumping well. The probability distributions of head become skewed in the area where there is significant drawdown. Figure 1.4b shows that the probability distributions for hydraulic head close to the pumping well are also strongly nonstationary in space. Figure 1.4c shows the probability distribution of hydraulic head when the scale of heterogeneity is increased to 1/10 of the simulation domain. Large-scale heterogeneity results clearly more 15 non-Gaussian behavior, which in turn may result a non-Gaussian distribution on other dependent variables. The shape of probability distributions becomes more complex and is influenced more by boundary conditions, nonstationary conductivity distribution, and the presence of sources and sinks. 20 MW13 10 0 1.54 1.61 10 0 1.68 1.12 Y 10 1.19 0.30 Head (m) 0.35 0.40 Source Constant Head= 0 m Head (m) 6 2 0 1.20 1.50 1.80 Head (m) 4 3 2 1 0 6 MW13 Density 4 MW3 Density Density σ 2 = 0.1 0 1.05 Head (m) Constant Head= 2 m 20 MW28 Density Density Density 20 MW3 100 4 σ 2 = 2 .0 Y 2 0 0.70 1.05 1.40 Head (m) 0 MW28 0.20 0.40 0.60 Head (m) 200 m Figure 1.2: Probability distributions of hydraulic head for stationary cases (Top: lnK variance =0.1, Bottom: lnK variance =2.0). Black circles are monitoring wells and thin solid lines are head contours in the model domain. 16 5 0 1.50 0 1.60 1.70 Head (m) 100 MW13 15 Density Density 10 MW3 8 6 4 2 0 Constant Head=0m Zone 3 Zone 2 Zone 1 Source Constant Head=2m Density 15 MW28 10 5 0 0.60 0.75 0.90 Head (m) 0.30 0.40 0.50 Head (m) 200 m Figure 1.3: Probability distributions of hydraulic head for Nonstationary case without pumping wells (Zone 1(lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day)). Black circles are monitoring wells and thin solid lines are head contours in the model domain. 17 3.0 MW3 Density Density 2.0 1.5 1.0 0.5 0.0 MW23 2.0 1.0 0.0 0.60 1.20 1.80 -2.00 0.00 Head (m) 2.0 1.5 1.0 0.5 0.0 3.0 D ns i y e t 6 4 MW 13 2 0 - 0. 0 2 0 00 . 0. 0 2 MW13 -1.00 0.00 1.00 MW18 100 MW28 2.0 1.0 0.0 -2.00 Head (m) 0 Density H ed m ) a ( Density Density 2.0 1.5 1.0 0.5 0.0 Constant Head =0m Zone 2 Source Constant Head =2m Zone 1 Zone 3 3 Zone Head (m) -1.00 -1.00 0.00 Head (m) -2.00 -1.00 0.00 Head (m) 200 m Figure 1.4a: Probability distributions of hydraulic head for Nonstationary case with pumping wells (Zone 1 (lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day)). Black circles are monitoring wells and thin solid lines are head contours in the model domain. 18 Density 15 MW9 15 Density 10 5 0 -0.60 Density 6 -0.40 MW7 10 5 0 -0.20 -0.60 Head (m) MW2 -0.40 -0.20 Head (m) 4 Prescribed Head 2 -1.00 -0.50 Head (m) Density 30 MW20 Prescribed Head -1.50 Prescribed Head 0 0 20 50 100 m Prescribed Head 10 0 -0.24 Head (m) -0.18 Density -0.30 Sub model domain 40 30 20 10 0 MW25 -0.25 -0.20 -0.15 Head (m) Figure 1.4b: Probability distributions of hydraulic head observed close to the pumping well. Black circles are monitoring wells and thin solid lines are the head contours in the model domain. 19 3.0 MW3 Density Density 2.0 1.5 1.0 0.5 0.0 0.60 1.20 MW23 2.0 1.0 0.0 1.80 -2.00 Head (m) -1.00 0.00 2.0 1.5 1.0 0.5 0.0 3.0 De ns i t y 6 4 2 M W 13 0 -1.00 0.00 1.00 0. 0 0 0. 20 He a d ( m ) MW18 100 MW28 2.0 1.0 0.0 -2.00 Head (m) 0 Density - 0. 0 2 MW13 Density 2.0 1.5 1.0 0.5 0.0 Constant Head =0m Zone 2 Source Density Constant Head =2m Zone 1 Zone 3 Head (m) -1.00 0.00 Head (m) -2.00 -1.00 0.00 Head (m) 200 m Figure 1.4c: Probability distributions of hydraulic head for Nonstationary case with pumping wells having large correlation scale (Zone 1 (lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day)). Black circles are monitoring wells and thin solid lines are head contours in the model domain. 20 Velocity Figures 1.5-1.7 show the observed probability distribution of groundwater velocity for different cases. It is interesting to note that probability distributions of longitudinal velocity and transverse velocity are structurally different. Figure 1.5a shows that the probability distribution of longitudinal velocity at low lnK variance is slightly asymmetric, but becomes strongly skewed at the higher lnK variance of 2.0. The Q-Q plot of longitudinal velocity (see Appendix, Figure 1.15) indicates that even at low variance the quantiles of longitudinal velocity deviates from the straight line and, as the lnK variance increases, the deviation becomes more prominent. The skewness in the velocity probability distribution can be directly attributed to the fact that velocity is proportional to hydraulic conductivity and its variability is dominated by the variability in hydraulic conductivity. Since hydraulic conductivity is assumed to be lognormal, the longitudinal velocity is expected to be strongly skewed. Figure 1.5b shows that the probability distribution of transverse velocity is interestingly always symmetrically distributed, even in the presence of strong heterogeneity. At a low lnK variance, the corresponding Q-Q plot (see Appendix, Figure 1.16) closely matches with the theoretical straight line, indicating normality. 21 For strongly heterogeneous media at high lnK variance of 2.0, the Q-Q plot shows deviation from the theoretical straight line, implying that the distribution, though symmetrical, is nonGaussian (see Appendix, Figure 1.16). Figures 1.6ab and 1.7ab show that the probability distributions of velocity in a particular direction in a nonstationary flow system are strongly nonstationary. As the direction of flow changes the shapes of the probability distributions in X-velocity and Y-velocity change significantly. However, if we focus our attention on the longitudinal and transverse velocities, their probability distributions change little with location and are not particularly sensitive to boundary conditions, aquifer stresses and other sources and sinks. 22 4 3 2 1 0 1.20 0.40 0.80 1.20 4 3 2 1 0 MW28 σ 2 = 0.1 Y 0.40 1.20 X-Velocity (m/d) Constant Head= 0 m X-Velocity (m/d) 0.80 Source 1.0 1.5 MW3 0.5 0.0 1.5 MW13 0.5 0.0 2.00 4.00 6.00 X-Velocity (m/d) 0 1.0 Density 1.5 Density Constant Head= 2 m X-Velocity (m/d) Density MW13 Density Density Density 4 MW3 3 2 1 0 0.40 0.80 100 0.00 2.40 4.80 X-Velocity (m/d) 1.0 0.5 MW28 σ 2 = 2 .0 Y 0.0 1.60 3.20 4.80 X-Velocity (m/d) 200 m Figure 1.5a: Probability distributions of X-velocity for stationary cases (Top: lnK variance =0.1, Bottom: lnK variance =2.0). Black circles are monitoring wells and thin solid lines are head contours in the model domain. 23 0.20 0.20 Density Constant Head= 0 m 1.0 0.0 3.0 1.0 0.0 -2.00 0.00 2.00 Y-Velocity (m/d) 0 MW13 2.0 Density 2.0 3.0 MW3 Density 3.0 0.1 Y-Velocity (m/d) 8 6 MW13 4 2 0 -0.20 0.00 Source Constant Head= 2 m Y-Velocity (m/d) Density 2 Y 8 σ = 6 MW28 4 2 0 -0.20 0.00 0.20 Y-Velocity (m/d) Density Density 8 6 MW3 4 2 0 -0.20 0.00 100 -2.00 0.00 2.00 Y-Velocity (m/d) 2.0 MW28 σ 2 = 2 .0 Y 1.0 0.0 -2.00 0.00 2.00 Y-Velocity (m/d) 200 m Figure 1.5b: Probability distributions of Y-velocity for stationary cases (Top: lnK variance =0.1. Bottom: lnK variance =2.0). Black circles are monitoring wells and thin solid lines are head contours in the model domain. 24 Zone 3 Zone 2 0.03 X-Velocity (m/d) 100 0.05 8 6 4 2 0 3.0 Density Density MW3 0.01 0 Constant Head=0m Zone 1 Source Constant Head=2m Density 80 60 40 20 0 MW13 MW28 2.0 1.0 0.0 0.24 0.48 0.72 X-Velocity (m/d) 0.00 1.00 2.00 X-Velocity (m/d) 200 m Figure 1.6a: Probability distributions of X-Velocity for Nonstationary case without pumping wells (Zone 1 (lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day)). Black circles are monitoring wells and thin solid lines are head contours in the model domain. 25 Zone 3 Zone 2 0 -0.01 0.00 Y-Velocity (m/d) MW13 15 Density Density Density 50 -0.03 10 5 0 100 MW28 10 5 0 0.00 0.24 Y-Velocity (m/d) 0 Constant Head=0m Zone 1 Constant Head=2m Source 15 100 MW3 -1.60 0.00 1.60 Y-Velocity (m/d) 200 m Figure 1.6b: Probability distributions of Y-Velocity for Nonstationary case without pumping wells (Zone 1 (lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day)). Black circles are monitoring wells and thin solid lines are head contours in the model domain. 26 2.0 Density Density 60 MW3 40 20 0 MW18 1.0 0.0 0.03 0.06 0.09 0.00 4.00 X-Velocity (m/d) 15 0.50 1.00 X-Velocity (m/d) 0 100 MW25 10 Density MW13 Density Density 4 3 2 1 0 Constant Head =0m Zone 2 Source Constant Head =2m Zone 1 Zone 3 X-Velocity (m/d) 2.00 5 0 -0.30 0.00 0.30 X-Velocity (m/d) 4 3 2 1 0 MW28 -2.00 -1.00 0.00 X-Velocity (m/d) 200 m Figure 1.7a: Probability distributions of X-Velocity for Nonstationary case with pumping wells (Zone 1 (lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day)). Black circles are monitoring wells and thin solid lines are head contours in the model domain. 27 60 MW28 Density Density 60 40 20 0 20 0 -0.04 -0.02 MW18 40 0.00 0.00 2.00 4.00 Y-Velocity (m/d) 0.00 0.46 Y-Velocity (m/d) 8 6 4 2 0 6 MW25 100 MW28 4 2 0 0.30 0.60 Y-Velocity (m/d) 0 Constant Head =0m Zone 2 MW13 Density 8 6 4 2 0 -0.46 Density Density Source Constant Head =2m Zone 1 Zone 3 Y-Velocity (m/d) -1.00 0.00 1.00 Y-Velocity (m/d) 200 m Figure 1.7b: Probability distributions of Y-Velocity for Nonstationary case with pumping wells (Zone 1 (lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day)). Black circles are monitoring wells and thin solid lines are head contours in the model domain. 28 Concentration Figures 1.8-1.9 show the observed probability distribution of solute concentration for both stationary and nonstationary flow cases. The results show that concentration probability distributions are the most complex, strongly non-Gaussian and nonstationary in both space and time even for lnK variance of 0.1. The Q-Q plot of concentration (see Appendix, Figure 1.17) show that simulated concentration probability distribution deviates substantially from the theoretical normal distribution. The results also show that the complex concentration probability structure can be best characterized based on its relative location to the mean plume. In particular, both for uniform and non-uniform flow situations, the results show that the concentration probability distributions at locations close to the source concentration (the specified source concentration is 100 ppm) are strongly skewed right. The probability distribution on the edge of the mean plume is strongly skewed toward the left or zero. The distributions at locations within the mean plume represent a transition between the two extreme distributions. Because of the extremely skewed nature in concentration probability, variance provides little information on concentration uncertainty and it is impractical fit accurately the highly nonstationary probability distribution using a predefined function. 29 0.030 MW18 Density MW 12 Density Density 0.20 0.15 0.10 0.05 0.00 0.020 0.010 0.000 20 40 60 Concentration (ppm) 30 60 90 Concentration (ppm) 0.030 MW23 0.020 0.010 0.000 20 40 60 Concentration (ppm) σ 2 = 0.1 MW3 MW14 0.03 0.5 0.0 97 98 99 10 20 Concentration (ppm) MW19 Constant Head =0m 1.0 Density Density 4 3 2 1 0 Density Source Constant Head =2m Y 0.02 0.01 0.00 20 40 60 Concentration (ppm) Concentration (ppm) 0 100 200 m Figure 1.8: Probability distributions of concentration for Stationary case. Black circles are monitoring wells and the color contours are the concentration. 30 0.02 0.01 0.00 40 30 MW8 20 10 0 10 20 30 Concentration (ppm) 1 Density MW3 Density Density 0.03 2.0 1.5 1.0 0.5 0.0 MW13 10 20 30 Concentration (ppm) 2 Density 4 2 0 6 12 Concentration (ppm) 0 Constant Head =0m Zone 2 MW5 100 0.020 0.015 0.010 0.005 0.000 0.15 MW4 Density 6 Density Zone 1 Source Constant Head =2m Zone 3 Concentration (ppm) MW9 0.10 0.05 0.00 30 60 90 Concentration (ppm) 20 40 60 Concentration (ppm) 200 m Figure 1.9: Probability distributions of concentration for Nonstationary case (Zone 1 (lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day)). Black circles are monitoring wells and the color contours are the concentration. 31 Seepage flux across a polyline Figures 1.10-1.11 present the observed probability distribution of seepage flux for different cases. The results show that the probability distributions of seepage flux is approximately Gaussian at low variance when it is integrated along a polyline. As the variance of lnK increases, the probability distributions of seepage flux become more non-Gaussian, but its skewness is less extreme than that in the original velocity. Spatial integration has the effect of reducing variability, rendering the transformation from lnK fluctuation to flux fluctuation more linear. Figure 1.18 shows that the seepage flux quantiles deviate from the theoretical straight line at high variance, indicating non-Gaussianity. Overall, the shapes of the seepage flux distributions depend on the scale of integration along the polyline and on degree of heterogeneity. The probability distribution of the seepage flux becomes more Gaussian as the scale of integration increases or as lnK variance decreases. 32 0.002 0.001 0.000 4800 5200 5600 100 3 Y 200 300 3 Total Seepage Flux (m /d) Source σ 2 = 2 .0 0.0.0006 Density Y Density Constant Head= 2 m Total Seepage Flux (m /d) σ 2 = 0.1 0.020 PL2 0.015 0.010 0.005 0.000 PL1 0.0.0004 0.0.0002 0.0000 4000 6000 8000 Total Seepage Flux (m /d) 100 PL2 0.0.004 0.0.002 0.000 3 0 0.0.006 Constant Head= 0 m Density Density 0.003 PL1 0 500 3 Total Seepage Flux (m /d) 200 m Figure 1.10: Probability distributions of seepage flux for stationary cases at each polylines. Thick solid lines represent the length of the polylines and thin solid lines are the head contours in the model domain. Top: lnK variance =0.1, Bottom: lnK variance =2.0. 33 Density Density 0.003 0.20 PL2 0.15 0.10 0.05 0.00 10 0.001 0.000 15 -900 -600 -300 20 3 3 Total Seepage Flux (m /d) PL6 70 140 210 PL8 300 600 3 Total Seepage Flux (m /d) 100 0.008 0.006 0.004 0.002 0.000 0 3 0 Constant Head= 0m Zone3 Zone2 0.020 0.015 0.010 0.005 0.000 Density Density Source Zone1 Total Seepage Flux (m /d) Constant Head= 2m PL10 0.002 Total Seepage Flux (m /d) 200 m Figure 1.11: Probability distributions of seepage flux at each polylines for nonstationary case with pumping wells (Zone 1 (lnK variance =0.5, λ= 2m, K= 2m/day), Zone 2 (lnK variance =3.0, λ= 20m, K= 200m/day) and Zone3 (lnK variance =1.5, λ= 10m, K= 20m/day)). Thick solid lines represent the length of the polylines and thin solid lines are the head contours in the model domain. 34 Solute Flux across a Polyline Figures 1.12-1.13 present the observed probability distribution of solute flux for different cases. The results show that the probability distribution of integrated solute flux is non-Gaussian except when lnK variance is very small. At a low lnK variance of 0.1, solute flux is slightly skewed, and the Q-Q plot in Figure 1.19 shows a slight departure from the theoretical straight line of normality. However, as lnK variance increases, the probability distributions in solute flux become strongly skewed, though less extreme than the original concentration probability distribution because of spatial integration. In general, the skewness of the distribution for solute flux across a polyline increases as the length of the polyline decreases and the degree of heterogeneity increases (see Appendix, Figure 1.19). Because of the strongly skewed nature of solute flux, variance in its distribution provides little information on uncertainty and it is impractical fit accurately the highly nonstationary probability distribution using a predefined function. 35 30000 Density Density 0.00020 PL2 0.00015 0.00010 0.00005 0.00000 10000 20000 Total Solute Flux (kg/d) 0.00020 0.00015 0.00010 0.00005 0.00000 PL6 8000 16000 24000 Density Y 0.00008 0.00006 0.00004 0.00002 0.00000 PL5 Constant Head= 0m σ 2 = 0.1 Source Constant Head= 2m Total Solute Flux (kg/d) 40000 60000 80000 0 100 200 m Density Total Solute Flux (kg/d) 0.00008 0.00006 0.00004 0.00002 0.00000 PL1 40000 60000 80000 Total Solute Flux (kg/d) Figure 1.12: Probability distributions of solute flux at each polylines for Stationary case. Thick solid lines represent the length of the polylines and the color contours are the concentration. 36 Density 0.00006 PL2 0.00004 0.00002 0 60000 Density 0.00000 120000 Total Solute Flux (kg/d) 0.00020 0.00015 0.00010 0.00005 0.00000 PL6 0 40000 Total Solute Flux (kg/d) Constant Head=0m Y Source Constant Head=2m σ 2 = 2 .0 Density 0.000015 PL5 0.000010 0.000005 0.000000 80000 160000 240000 100 200 m 0.000015 Density 0 PL1 Total Solute Flux (kg/d) 0.000010 0.000005 0.000000 100000 200000 300000 Total Solute Flux (kg/d) Figure 1.13: Probability distributions of solute flux at each polylines for Stationary case. Thick solid lines represent the length of the polylines and the color contours are the concentration. 37 1.7 CONCLUSIONS The following key findings emerge from this research: • The statistical structure of a groundwater system under general, nonstationary condition is complex and cannot be characterized by a few simple statistical moments. Even when lnK is assumed to be Gaussian, the dependent hydrogeological variables can be strongly nongaussian, nonstationary, and for those that are directional, anisotropic. • The closest to Gaussian, among all variables simulated, is the probability distribution of hydraulic head, when the scale of heterogeneity is much smaller than the domain size (e.g., 1/33of the problem scale). • The probability distribution of hydraulic head in response to larger scale heterogeneity (e.g., greater 1/10 of the problem scale) is often non-Gaussian and nonstationary and cannot be accurately fitted with a single functional distribution. • The probability distribution of groundwater velocity is direction-dependent. The probabilistic structure of longitudinal velocity is strongly skewed. • The probability distribution of transverse velocity is approximately Gaussian or symmetrically distributed, even in the presence of strong heterogeneity. • Seepage flux across a transect – representing integration of strongly skewed velocity, is approximately Gaussian, even in the presence of strong heterogeneity. • Probability distribution of solute concentration is far from Gaussian and strongly nonstationary when the source size is relatively small (e.g., less than 10 lnK correlation 38 scales), even for aquifers with very weak heterogeneity. The complex concentration probability distribution cannot be fitted with a single functional distribution. • Solute flux across a transect, unlike its seepage flux counterpart, is strongly non-Gaussian and location-dependent. This research represents the first comprehensive analysis of the probabilistic structure of basic hydrogeological variables and findings from this work have significant implications on theoretical and practical stochastic subsurface hydrology. 39 APPENDIX 40 APPENDIX STATISTICAL PARAMETERS OF GROUNDWATER VARIABLES Table 1.2: Statistical parameters (Skewness and Kurtosis)of Head, X-Velocity and Yvelocity for low and high lnK Variances Head lnK variance=0.1 lnK variance=2.0 Nonstationary without pumping well MW3 Skewness -0.04 -0.54 -0.16 Kurtosis -0.05 0.45 -0.09 MW13 Skewness -0.04 -0.06 0.11 Kurtosis -0.01 -0.17 0.11 MW28 Skewness 0.04 0.54 0.12 Kurtosis 0.05 0.45 0.02 0.29 10.10 5.08 0.51 2.32 2.09 0.49 9.76 6.31 0.42 6.66 4.23 -0.04 0.25 0.07 0.29 10.77 9.59 X-Velocity lnK variance=0.1 lnK variance=2.0 Nonstationary without pumping well 0.51 2.32 0.64 0.49 9.76 0.67 0.46 2.38 1.78 Y-Velocity lnK variance=0.1 lnK variance=2.0 Nonstationary without pumping well 0.04 -0.25 -0.50 0.29 10.77 0.60 -0.00 0.17 1.25 41 Table 1.3 : Statistical parameters (Skewness and Kurtosis)of Head, X-Velocity and YVelocity for Nonstationary with pumping wells Head X-Velocity Y-Velocity Skewness MW3 MW13 MW18 MW23 MW28 Kurtosis Skewness Kurtosis Skewness -0.18 0.19 -0.01 -0.67 -0.23 -0.07 0.13 0.22 1.59 0.07 0.70 1.65 2.09 -1.96 -2.31 0.95 4.34 6.50 5.47 8.10 -0.51 1.60 1.86 1.78 1.23 Table 1.4: Statistical parameters (Skewness and Kurtosis)of Head, for Nonstationary with pumping wells having large scale heterogeneity MW3 MW13 MW18 MW23 MW28 Head Skewness -0.67 -0.12 -0.50 -0.77 -1.11 Kurtosis 0.40 0.88 2.59 3.37 4.73 Table 1.5: Statistical parameters of Head in the sub model domain for nonstationary case with pumping well Skewness Kurtosis MW2 -5.76 49.12 MW7 -1.45 4.45 42 MW9 -1.39 5.40 MW20 -0.64 0.39 MW23 -1.23 4.09 Kurtos is 0.64 6.26 10.48 7.74 10.17 Table 1.6 : Statistical parameters (Skewness and Kurtosis) of Concentration For low variance of lnK and for nonstationary case with pumping wells lnK variance=0.1 MW5 MW12 MW14 MW18 MW19 MW23 Skewness -6.20 1.57 3.63 -0.84 2.18 2.18 Non stationary with pumping well Kurtosis 57.17 2.28 19.27 -0.04 5.95 4.95 Skewness 2.02 -0.07 9.92 13.32 1.11 5.08 MW3 MW4 MW5 MW8 MW9 MW13 Kurtosis 6.03 -1.04 147.59 276.68 0.69 33.63 Table 1.7: Statistical parameters (Skewness and Kurtosis)of seepage and solute fluxes For low and high lnK variances Seepage Flux Skewness ( lnK variance=0.1) Kurtosis ( lnK variance=0.1) Skewness ( lnK variance=2.0) Kurtosis ( lnK variance=2.0) Solute Flux PL1 0.06 -0.12 0.36 0.16 PL1 0.26 0.08 1.18 2.30 43 PL2 0.40 0.31 1.86 5.68 PL2 0.40 0.29 1.92 6.20 PL5 0.23 0.19 1.01 1.46 PL6 0.04 0.21 2.22 8.15 QUANTILES-QUANTILES PLOTS Q-Q Plot of Hydraulic Head versus Standard Normal Distribution MW 3 Quantiles of Head lnK variance= 0.1 Theoretical Data Standard Normal Quantiles Figure 1.14: Q-Q plot of Hydraulic head data versus normal distribution for low variance and high variance cases. MW3: located near the west boundary, MW28: Located near the east boundary, MW13: located approximately middle of the model domain. 44 Q-Q Plot of Hydraulic Head versus Standard Normal Distribution MW 3 Quantiles of Head lnK variance= 2.0 Theoretical Data Standard Normal Quantiles Figure 1.14 (cont’d). 45 Q-Q Plot of Hydraulic Head versus Standard Normal Distribution MW 13 Quantiles of Head lnK variance= 0.1 Theoretical Data Standard Normal Quantiles Figure 1.14 (cont’d). 46 Q-Q Plot of Hydraulic Head versus Standard Normal Distribution MW13 Quantiles of Head lnK variance= 2.0 Theoretical Data Standard Normal Quantiles Figure 1.14 (cont’d). 47 Q-Q Plot of Hydraulic Head versus Standard Normal Distribution MW28 Quantiles of Head lnK variance=0.1 Theoretical Data Standard Normal Quantiles Figure 1.14 (cont’d). 48 Q-Q Plot of Hydraulic Head versus Standard Normal Distribution MW 28 Quantiles of Head lnK variance= 2.0 Theoretical Data Standard Normal Quantiles Figure 1.14 (cont’d). 49 Q-Q Plot of X-Velocity versus Standard Normal Distribution MW 3 Quantiles of X-Velocity lnK variance= 0.1 Theoretical Data Standard Normal Quantiles Figure 1.15: Q-Q plot of X-Velocity data versus normal distribution for low and high lnK variance 50 Q-Q Plot of X-Velocity versus Standard Normal Distribution MW 3 Quantiles of X-Velocity lnK variance= 2.0 Theoretical Data Standard Normal Quantiles Figure 1.15 (cont’d). 51 Q-Q Plot of Y-Velocity versus Standard Normal Distribution MW 3 Quantiles of Y-Velocity lnK variance= 0.1 Theoretical Data Standard Normal Quantiles Figure 1.16: Q-Q plot of Y-Velocity data versus normal distribution for low and high lnK variance 52 Q-Q Plot of Y-Velocity versus Standard Normal Distribution MW 3 Quantiles of Y-Velocity lnK variance= 2.0 Theoretical Data Standard Normal Quantiles Figure 1.16 (cont’d). 53 Q-Q Plot of Concentration versus Standard Normal Distribution MW 3 Quantiles of Concentration lnK variance= 0.1 Theoretical Data Standard Normal Quantiles Figure 1.17: Q-Q plot of concentration data versus normal distribution 54 Q-Q Plot of Concentration versus Standard Normal Distribution MW 12 Quantiles of Concentration lnK variance= 0.1 Theoretical Data Standard Normal Quantiles Figure 1.17 (cont’d). 55 Q-Q Plot of Seepage Flux versus Standard Normal Distribution PL1 Quantiles of Seepage Flux lnK variance= 0.1 Theoretical Data Standard Normal Quantiles Figure 1.18: Q-Q plot of seepage flux data versus normal distribution for low and high lnK variance 56 Q-Q Plot of Seepage Flux versus Standard Normal Distribution PL1 Quantiles of Seepage Flux lnK variance= 2.0 Theoretical Data Standard Normal Quantiles Figure 1.18 (cont’d). 57 Q-Q Plot of Seepage Flux versus Standard Normal Distribution PL2 Quantiles of Seepage Flux lnK variance= 0.1 Theoretical Data Standard Normal Quantiles Figure 1.18 (cont’d) 58 Q-Q Plot of Seepage Flux versus Standard Normal Distribution PL2 Quantiles of Seepage Flux lnK variance= 2.0 Theoretical Data Standard Normal Quantiles Figure 1.18 (cont’d). 59 Q-Q Plot of Solute Flux versus Standard Normal Distribution x 4 10 PL1 Quantiles of Solute Flux lnK variance= 0.1 Theoretical Data Standard Normal Quantiles Figure 1.19: Q-Q plot of solute flux data versus normal distribution for low and high lnK variance 60 Q-Q Plot of Solute Flux versus Standard Normal Distribution X 10 5 PL1 Quantiles of Solute Flux lnK variance= 2.0 Theoretical Data Standard Normal Quantiles Figure 1.19 (cont’d). 61 Q-Q Plot of Solute Flux versus Standard Normal Distribution X 10 4 PL2 Quantiles of Solute Flux lnK variance= 0.1 Theoretical Data Standard Normal Quantiles Figure 1.19 (cont’d). 62 Q-Q Plot of Solute Flux versus Standard Normal Distribution X 4 10 PL2 Quantiles of Solute Flux lnK variance= 2.0 Theoretical Data Standard Normal Quantiles Figure 1.19 (cont’d). 63 REFERENCES 64 REFERENCES 1. Ballio, F., and A. Guadagnini (2004), Convergence assessment of numerical Monte Carlo simulations in groundwater hydrology, Water Resour. Res., 40, W04603, doi:10.1029/2003WR002876. 2. Bellin, A., and D. Tonina (2007), Probability density function of nonreactive solute concentration in heterogeneous porous formations, J. Contam. Hydrol., 94(1 – 2), 109 – 125. 3. Bellin, A., P. Salandin, and A. Rinaldo (1992a), Simulation of dispersion in heterogeneous porous formations: Statistics first-order theories, convergence of computations, Water Resour. Res., 28(9), 2211-2227. 4. Carle, S.F., and G.E. Fogg (1996), Transition probability-based indicator geostatistics, Math Geol., 28(4), 453–476. 5. Caroni, E., and V. Fiorotto (2005), Analysis of concentration as sampled in natural aquifers, Transp. Porous Media, 59(1), 19–45. 6. Chambers, J., W. Cleveland, B. Kleiner, P. Tukey (1983), Graphical Methods for Data Analysis, Wadsworth. 7. Cushman, J. H. (1997), The Physics of Fluids in Hierarchical Porous Media, Angstroms to Miles, Kluwer Acad. Publ.,Dordrecht. 8. Dagan, G. (1989), Flow and Transport in Porous Formations, Springer-Verlag, New York 9. Dagan, G. (2002), An overview of stochastic modeling of groundwater floe and transport: from theory to applications, EOS, 83(53). 10. Dagan, G.(2004), On application of stochastic modeling of groundwater flow and transport, Stochastic Environmental Research and Risk assessment, 18(4),266-267. 11. Dagan, G., (1994), An exact nonlinear correction to transverse macrodispersivity for transport in heterogeneous formations, Water Resour. Res., 30(10), 2699-2705. 12. Darvini, G., and P. Salandin (2006), Nonstationary flow and non ergodic transport in random porous media, Water Resour. Res., 42, W12409, doi:10.1029/2005WR004846. 13. Darvini, G., and P.Salandin (2006), Nonstationary flow and nonergodic transport in random porous media, Water Resour. Res.42, W12409, doi:10.1029/2005WR004846. 14. de Marsily, G., F. Delay, J. Goncxalve`s, P. Renard, V. Teles, and S. Violette (2005), 65 Dealing with spatial heterogeneity, Hydrogeology Journal 13(1), 161–183. 15. Englert, A., J. Vanderborght, and H. Vereecken (2006), Prediction of velocity statistics in three-dimensional multi-Gaussian hydraulic conductivity fields, Water Resour. Res., 42, W03418, doi:10.1029/2005WR004014. 16. Fiorotto,V., and E. Caroni (2002), Solute concentration statistics in heterogeneous aquifers for finite Peclet values, Transport Porous Media, 48, 331–351. 17. Franssen, H.J.W.M.H., and J.J. Go`mez-Herna`ndez (2002), 3D inverse modeling of groundwater flow at a fractured site using a stochastic continuum model with multiple statistical populations, Stochastic Environmental Research and Risk Assessment 16(2), 155–174. 18. Gelhar, L. W., (1986), Stochastic subsurface hydrology: from theory to applications, Water Resour. Res., 22(9),135S. 19. Gelhar, L. W., A. L. Gutiahr, and R. L. Naff (1979), Stochastic analysis of macrodispersion in a stratified aquifer, Water Resour. Res., 15, 1387–1397. 20. Gelhar, L. W., and C. L. Axness (1983), Three-dimensional stochastic analysis of macrodispersion in aquifers, Water Resour. Res., 19(1), 161-180. 21. Gelhar, L.W. (1993), Stochastic subsurface hydrology, Englewood Cliffs, NJ: PrenticeHall. 22. Ginn, T.R. (2004), On the application of stochastic approaches in hydrogeology, Stochastic Environmental Research and Risk assessment, 18(4), 282-284. 23. Guadagnini, A., and S.P. Neuman (1999a), Nonlocal and localized analysis of conditional mean steady state flow in bounded randomly nonuniform domains 1. Theory and computational approach, Water Resour. Res.35 (10), 2999-3018. 24. Guadagnini, A., and S.P. Neuman (1999b), Nonlocal and localized analysis of conditional mean steady state flow in bounded randomly nonuniform domains 2. Computational examples, Water Resour. Res.35 (10), 3019-3039. 25. Hamed, M.M., P.B. Bedient, and C.N. Dawson (1996), Probabilistic modeling of aquifer heterogeneity using reliability methods, Adv. Water Resour., 19(5), 277-295. 26. Kapoor, V., and P. K. Kitanidis (1997), Advection-diffusion in spatially random flows: Formulation of concentration covariance, Stochastic Hydrol.Hydraul., 11, 397–422. 27. Li, L., A. Hamdi, Tchelepi, D. Zhang (2003), Perturbation-based moment equation approach for flow in heterogeneous porous media: applicability range and analysis of high-order terms, Journal of Computational Physics, 188, 296-317. 66 28. Li, S.G., and Q. Liu (2003), Interactive Ground Water (IGW): An Innovative Digital Laboratory for Groundwater Education and Research, Computer Applications in Engineering Education, 11(4), 179-202. 29. Li, S.G., and Q. Liu (2006), A Real-time, Computational Steering Environment for Integrated Groundwater Modeling, Journal of Ground Water, 44(5), 758-763. 30. Li, S.G., D. McLaughlin, H.S. Liao (2003), A computationally practical method for stochastic groundwater modeling, Adv. Water Resour., 26(11), 1137-1148. 31. Li, S.G., D.B. McLaughlin, and H.S.Liao (2003), A computationally practical method for stochastic groundwater modeling, Adv. Water Resour., 26(11), 1137-1148. 32. Li, S.G., H.S. Liao, and C.F. Ni (2004), Stochastic modeling of complex nonstationary groundwater systems, Adv. Water Resour., 27(11), 1087-1104. 33. Li, S.G., Q. Liu, and S. Afshari (2006), An Object-Oriented Hierarchical Patch Dynamics Paradigm (HPDP) for Groundwater Modeling, Environmental Modeling and Software, 21(5), 601-758. 34. Lu, Z., and D. Zhang (2002), On stochastic modeling of flow in multimodal heterogeneous formations, Water Resour. Res. 38(10):10.1029/2001WR001026. 35. Lu, Z., and D. Zhang (2003), On importance sampling Monte Carlo approach to uncertainty analysis for flow and transport in porous media, Adv. Water Resour., 26, 1177-1188. 36. Lu, Z., S.P. Neuman, A. Guadagnini, D.M. Tartakovsky (2002), Conditional moment analysis of steady state unsaturated flow in bounded, randomly heterogeneous soils, Water Resour. Res. 38(4): 10.1029/2001WR000278. 37. Molz, F. (2004), A rational role for stochastic concepts in subsurface hydrology: a personal perspective, Stochastic Environmental Research and Risk assessment, 18(4), 278-279. 38. Morales-Casique, E., S. P. Neuman, and A. Guadagnini (2006a), Non-local and localized analyses of non-reactive solute transport in bounded randomly heterogeneous porous media: Theoretical framework, Adv. Water Resour., 29(8), 1238– 1255. 39. Morales-Casique, E., S. P. Neuman, and A. Guadagnini (2006b), Nonlocal and localized analyses of nonreactive solute transport in bounded randomly heterogeneous porous media: Computational analysis, Adv. Water Resour., 29(9), 1399–1418. 40. Neuman, S. P. (1997), Stochastic approach to subsurface flow and transport: A view to the future, in Subsurface Flow and Transport: A Stochastic Approach, edited by G. 67 Dagan and S. P. Neuman, pp. 231-241, Cambridge Univ. Press, New York. 41. Neuman, S.P. (2004), Stochastic groundwater models in practice, Stochastic Environmental Research and Risk Assessment,18(4),268-270. 42. Ni, C.F. and S.G. Li (2006), Modeling Groundwater Velocity Uncertainty in Complex Composite Media, Adv. in Water Resources, 29(12), 1866-1875. 43. Ni, C.F., and S.G. Li (2005), Simple closed form formulas for predicting groundwater flow model uncertainty in complex, heterogeneous trending media, Water Resour. Res., 41, W11503, doi:10.1029/2005WR004143. 44. Nowak, W., R. L. Schwede, O. A. Cirpka, and I. Neuweiler (2008), Probability density functions of hydraulic head and velocity in three-dimensional heterogeneous porous media, Water Resour. Res., 44, W08452, doi:10.1029/2007WR006383. 45. Renard, P. (2007), Stochastic hydrogeology: What Professionals really need? Ground Water, 45(5), 531-541. 46. Rubin, Y. (2003), Applied Stochastic Hydrogeology, Oxford Univ. Press, New York 47. Rubin, Y. (2004), Stochastic hydrogeology-challenges and misconceptions, Stochastic Environmental Research and Risk assessment, 18(4), 280-281. 48. Sahuiquillo, A., J.E. Capilla, J.J. Go`mez-Herna`ndez, J. Andreu (1992),Conditional simulation of transmissivity fields honoring piezometric head data, Hydraulic engineering software IV, Fluid flow modeling, Vol. II, Blair WR, Cabrera, E (eds) Elsevier Oxford,201–214. 49. Salandin, P. and V. Fiorotto (1998), Solute Transport in highly heterogeneous aquifers, Water Resour. Res.34 (5), 949-961. 50. Salandin, P., and V.Fiorotto (2000), Dispersion tensor evaluation in heterogeneous media for finite Peclet values, Water Resour. Res., 36, 1449–1455. 51. Schwede, R.L., O.A. Cirpka, W. Nowak, I. Neuweiler (2008), Impact of sampling volume on the probability density function of steady state concentration, Water Resour. Res., 44, W12433, doi:10.1029/2007WR006668. 52. Winter, C. L., D.M. Tartakovsky, and A.Guadagnini (2003), Moment Differential Equations for Flow in Highly Heterogeneous Porous Media. Surveys in Geophysics, 24, 81-106. 53. Winter, C.L. (2004), Stochastic hydrology: Practical alternative exist, Stochastic Environmental Research and Risk Assessment, 18(4), 271-273. 68 54. Winter, C.L., and D.M. Tartakovsky (2002), Groundwater flow in heterogeneous composite aquifers, Water Resour Res. 38(8):10.1029/2001WR000450. 55. Zhang, D. (2001), Stochastic Methods for Flow in Porous Media: Coping With Uncertainties, Academic Press, San Diego, CA. 56. Zhang, D., and S.P. Neuman (1995a), Eulerian-Lagrangian analysis of transport conditioned on hydraulic data: 1. Analytical-numerical approach, Water Resour. Res.31 (1), 39-51. 57. Zhang, D., and S.P. Neuman (1995b), Eulerian-Lagrangian analysis of transport conditioned on hydraulic data: 2. Effect of log trasmissivity and hydraulic head measurements, Water Resour. Res.31 (1), 53-63. 58. Zhang, D., and S.P. Neuman (1995c), Eulerian-Lagrangian analysis of transport conditioned on hydraulic data: 3. Spatial moments, travel time distribution, mass flow rate and cumulative release across a Compliance surface, Water Resour. Res.31 (1), 6575. 59. Zhang, D., and S.P. Neuman (1995d), Eulerian-Lagrangian analysis of transport conditioned on hydraulic data: 4. Uncertain initial plume state and Non-Gaussian velocities, Analytical-numerical approach, Water Resour. Res.31 (1), 77-88. 60. Zhang, Y.-K., and D. Zhang (2004), Forum: The state of stochastic hydrology, Stochastic Environmental Research and Risk assessment, 18(4), 265. 61. Zheng, C., and S.M. Gorelick (2003), Analysis of solute transport in flow fields influenced by preferential flow paths at the decimeter scale, Ground Water, 41, 142–155. 69 PAPER 2 PROBABILISTIC WELL CAPTURE ZONE DELINEATION USING A STATEWIDE GROUNDWATER DATABASE SYSTEM Dipa Dey and Huasheng Liao Department of Civil and Environmental Engineering Michigan State University, East Lansing, Michigan 48824 2.1 ABSTRACT It is widely known that data limitation is a major practical bottleneck in stochastic groundwater modeling. An expert panel at a recent forum on the state of practice in stochastic subsurface modeling and hydrology suggested that new measurement technologies or data sources of much better resolution are needed if we are to enable routine application of stochastic methods in real world investigations. In this paper, we explore the use of a recently-developed statewide database in Michigan for stochastic groundwater modeling and apply it to probabilistic capture zone delineation. We are particularly interested in testing if the relatively crude, but detailed datasets can be used to characterize aquifer heterogeneity with sufficient details to enable practical stochastic modeling in the context of source water protection. This research represents an extension of our earlier work that uses this data source for basic deterministic modeling and cost effective source water management in Michigan. Our earlier research clearly shows that this new data source, despite its significant uncertainty, is effective in modeling regional mean flow patterns, directions and particularly source water protection areas. Our earlier research also shows that capture zone delineation based on regional patterns without taking into account smaller scale heterogeneity 70 can significantly underestimate the size of capture zone. As a pragmatic way to improve the delineation, managers and practitioners often manually draw a “grey area” around “skinny” capture zones to account for the effect of un-modeled heterogeneity. We attempt in this research to quantify these grey areas probabilistically in a systematic manner. In particular, we consider the following four approaches to delineate well capture zones, each representing a different way to conceptualize aquifer heterogeneity. • Deterministic approach: modeling large scale heterogeneity as a deterministic trend, ignoring the uncertainty in large-scale heterogeneity and the effect of small scale heterogeneity; • Stochastic unconditional method: modeling heterogeneity as a random field using Monte Carlo simulation; • Stochastic semiconditional method: modeling large scale heterogeneity as a deterministic trend and smaller scale heterogeneity as a random field using unconditional Monte Carlo simulation; • Stochastic macrodispersion method: modeling large-scale heterogeneity as a deterministic trend and small scale heterogeneity deterministically using effective macrodispersion. These approaches were applied to several selected sites with distinct characteristics to delineate well capture zones. The results are systematically analyzed and compared. The following findings emerge from this research: • A significant portion of the variability in the statewide conductivity dataset is uncorrelated noise and this data noise can be quantified through variogram modeling. • Despite the data noise, two scales of natural variability in conductivity can be identified; 71 • Kriging based on the large-scale variograms, with small scale variability and data noises represented as nuggets, can be used to delineate large scale trends. • Variogram modeling of detrended conductivity data can be used to identify smaller scale heterogeneity. The residual variograms are significantly noisier, but statistically meaningful structures can be detected; • Deterministic delineation that accounts only for the effect of large scale trend leads to capture zones that are significantly smaller than its stochastic counterpart, especially in areas where regional gradient is relatively strong; • Stochastic Monte Carlo simulation provides a map of well capture probability. Combined with the statewide database, MC provides a useful, systematic tool for risk-based decision making; • Stochastic macrodispersion method provides a computationally efficient alternative to delineate wellhead protection area in a way that accounts for the effect of both small and large scale heterogeneity. Future research should focus on 3D stochastic modeling using 3D lithology and calibrating the conductivity to statewide static water level dataset. 2.2 MOTIVATION AND OBJECTIVE It is widely known that data limitation is a major practical bottleneck in stochastic groundwater modeling. An expert panel at a recent forum on the state of practice in stochastic subsurface modeling and hydrology suggested that new measurement technologies or data sources of much better resolution are needed if we are to enable routine application of stochastic methods in real 72 world investigations [Dagan, 2002; Harrar et al., 2003; Zheng and Gorelick, 2003; Molz et al., 2003; Zhang and Zhang 2004; Winter, 2004; Dagan 2004; Neuman, 2004; Molz, 2004]. In this paper, we explore the use of a recently-developed statewide database system in Michigan [GWIM, 2006; Wellogic, 2006; Map Image viewer, RS&GIS-MSU] for stochastic groundwater modeling and apply it to probabilistic well capture zone delineation. We are particularly interested in exploring if the relatively crude, but detailed water well records can be used to characterize aquifer heterogeneity with sufficient details to enable practical stochastic modeling in the context of source water protection. This research represents an extension of earlier work that used this database for basic deterministic modeling and improved source water management [Simard, 2007; Li et al., 2009; Oztan, 2010]. This earlier research clearly showed that this new data source, despite its significant uncertainty, is effective in modeling regional mean flow patterns, flow directions and particularly source water protection areas. This earlier research also showed that capture zone delineation based on regional patterns without taking into account smaller scale heterogeneity can significantly underestimate the size of well capture zone. The significant effect of heterogeneity on groundwater flow and transport has been widely recognized and an area of intense research in past decades [Dagan, 1989; Gelhar 1993; Rubin, 2003]. Variability in aquifer properties translates into uncertainty in groundwater seepage velocity, which in turn causes uncertainty in contaminant transport. A number of researchers recently stressed that spatial variability of aquifer properties (e.g., hydraulic conductivity) is one of the main sources of uncertainty in capture zone delineation [e.g., Varljen and Shafer, 1991; Wheater et al., 1998; van Leeuwen et al., 1998; Feyen et al., 2001; Stauffer et al., 2002; Feyen et al., 2003; Hendricks Franssen et al., 2004b; Stauffer et al., 2005] and have used stochastic 73 approaches to model well capture zone areas [Varljen and Shafer, 1991; Bair et al., 1991; Vassolo et al., 1998; van Leeuwen et al., 2000; Feyen et al., 2001]. A stochastic approach simulates uncertain, heterogeneous input parameter fields as random fields and allows modeling not only the area of well capture but also the probability of capture, providing a systematic, riskbased framework for decision making [e.g., Franzetti and Guadagnini, 1996; Guadagnini and Franzetti, 1999; Riva et al., 1999; Delhomme, 1978, 1979; Varljen and Shafer, 1991; van Leeuwen et al., 1998, 1999; Starn et al., 2000; Van Leeuwen et al., 2000; Riva et al., 2006]. The Achilles’ heel of stochastic approaches, however, is their intensive data requirements [Dagan, 2002; Harrar et al., 2003; Zhang and Zhang 2004, Winter, 2004; Dagan 2004; Neuman, 2004; Molz, 2004]. Unlike deterministic modeling, stochastic approaches require not only characterizing the means of input fields but also modeling their variogram structures to estimate the geostatistical parameters. Datasets allowing such estimation are often very difficult to obtain in practice. A literature review on stochastic groundwater modeling over the years shows that stochastic approaches are rarely used outside of research applications [Harrar et al., 2003]. In particular, the review shows current efforts in stochastic capture zone modeling suffer from the following limitations: 1. The vast majority of the stochastic modeling studies focus on hypothetical sites; assuming known statistical structure of small-scale heterogeneity [e.g., Varljen and Shafer 1991; Cole 1995; Franzetti and Guadagnini, 1996; Kinzelbach et al., 1996; Cole and Silliman, 1996; van Leeuwen et al.,1998; Hendricks Franssen et al., 2004b]. 74 2. Among the relatively few site-specific, real-world applications, the datasets used are often highly limited - too sparse to reliably characterize real world heterogeneity [e.g., Johanson, 1992; Vassolo et al., 1998; Riva et al., 2006; Theodossiou et al., 2005]. Because of this, as Gelhar [1993], stressed there is often significant uncertainty about predicted uncertainty in stochastic modeling. 3. While the original motivation of stochastic modeling was to enable better coping with complex realities, most of the newly developed stochastic methods (e.g., perturbation methods) apply only for simple, statistically uniform processes, free of complex sources and sinks, boundary conditions, and nonstationarities [Li and McLaughlin, 1993, 1995; Li, Liao and Ni, 2004]. 4. Some of the stochastic methods developed only allow representing uncertain aquifer properties or stresses as random constants [e.g., Evers and Lerner, 1998, Kinzelbach et al., 1996]. Future direction of stochastic subsurface hydrology research should clearly focus on the development of approaches, data sources, or measurement technologies that allow modeling problems of realistic sizes and complexities [Zhang and Zhang, 2004; Ginn, 2004; Neuman, 2004]. Stochastic methods need data of good resolution and must allow modeling both largescale (systematic) and small-scale (random) heterogeneity and stresses before they can become routine tools for practical applications. 75 2.3 DATA INTENSIVE STATEWIDE MODELING SYSTEM Over the past decade, the Michigan Department of Natural Resources and Environment (DNRE) has engaged in a major statewide data integration effort, capitalizing on the recent dramatic developments in remote sensing, Global Position System (GPS), and Geographic Information System (GIS) technologies. In this process, the DNRE has created a statewide GIS-based groundwater database system [Wellogic, 2006; GWIM, 2006; Map Image viewer, RS&GIS-MSU] that can be used to model regional scale processes and potentially local flow systems. The network of databases contains virtually “all” data needed for typical groundwater flow simulations [Li et al., 2009]. In particular, the databases contain data for defining model areas (e.g., watershed delineations, the Great Lakes, streams, and inland lakes), aquifer elevations (e.g., digital elevation model - DEM, lithologies), aquifer properties (e.g., lithologies, specific capacities), aquifer heterogeneity (lithologies, glacial land systems), aquifer stresses (e.g., estimated recharge, streams, lakes, wetlands, drains, and wells), surface water levels (e.g., DEM), contamination sources (e.g., industrial sites, landfills, leaky underground storage tanks), as well as, data for model calibration (e.g., static water levels, estimated baseflows, aquifer test hydrographs, and hydrographs from long term USGS monitoring wells). To enable systematic use of this vast “new” data source, Li et al. [2009] have recently developed an intelligent link between the databases and a real-time hierarchical groundwater modeling system [Li and Liu, 2006, 2008; Li et al., 2009]. The result is a data-intensive, statewide platform that can be used to model groundwater flow virtually anywhere (the areas covered by the database) in Michigan’s aquifers. The modeling system is “live-linked” in a hierarchical fashion to Michigan’s streams, lakes, geology, and water wells. 76 The hierarchical modeling system simulates state’s two uppermost aquifer layers – the glacial drift aquifer and/or the bedrock aquifer. Aquifer elevations are defined based on approximately half a million of well logs, including logs from water wells, as well as, oil and gas wells. Hydraulic properties are estimated using more than 300,000 water well records, including test pumping information and lithologic descriptions. Natural recharge to the glacial aquifer is estimated using data on land use, soil conditions, watershed characteristics, and observed stream flow hydrographs at more than 400 gaging stations statewide [GWIM,2006]. The Great Lakes are represented in the glacial aquifer as prescribed heads equal to the long term average of observed lake levels. Large inland lakes and rivers are represented as head-dependent fluxes. Wetlands, ponds, small lakes, and small streams are treated by default as drains. A default leakance value is assigned to all surface water bodies based on their size class or stream order. A steady state water level elevation is assigned to all lakes, wetlands, and stream arcs based on the 10 to 30 m statewide DEM. All default representations and parameter values in the modeling system can be interactively modified. The estimation of hydraulic conductivity in the database is briefly presented below. Hydraulic Conductivity: The horizontal hydraulic conductivity in the database was calculated from water well records using the following equation [Freeze and Cherry, 1979]: n ∑ Ki Bi K h = i =1 B 77 (2.1) where, Kh is the equivalent horizontal hydraulic conductivity averaged over saturated aquifer thickness, n is the total number of layers in the unit, Ki is the hydraulic conductivity of lithologic layer i, Bi is the thickness of the layer i, B is the total thickness of the saturated unit. Since hydraulic conductivity depends on not only local lithologies but also the regional geological processes that produce them (e.g., sand in a lacustrine deposition environment may be less permeable than that in a glacial outwash environment), a range of three typical hydraulic conductivity values (minimum, mean, maximum) were defined for each lithological layer. One of the three values was assigned for the corresponding lithology, depending on the glacial land system that the water wells were drilled in [GWIM, 2006; State of Michigan, 2006]. 2.4 RESEARCH APPROACHES In this paper, we investigate how effectively the Michigan statewide groundwater database can be used to probabilistically model Michigan’s groundwater systems. We consider three different approaches to model local groundwater flow and to delineate well capture zones, each representing a different way of conceptualizing aquifer heterogeneity. These approaches are briefly described below: 1. Deterministic Approach: This is our baseline approach. We model large-scale heterogeneity as deterministic trends and the effects of small-scale heterogeneity are ignored. The approach involves solving the groundwater flow equation and performing reverse particle tracking to delineate well capture zone. 78 2. Stochastic Monte Carlo Approach: This approach builds on the deterministic approach and we model both large-scale heterogeneity as deterministic trend and small-scale heterogeneity as a random field using Monte Carlo simulation (MC). We choose to use MC over the newer stochastic perturbation methods because of its conceptual simplicity, general applicability, and ability to fully quantify uncertainty in model output. MC involves directly generating many likely realizations of aquifer flow and particle tracking models in such a manner that they reflect the observed parameter uncertainty. The results are subsequently analyzed in a statistical manner to quantify the uncertainty inherent in the expected result. Combining with Michigan’s statewide databases, MC represents a general and potentially very powerful tool for risk-based groundwater investigations. 3. Stochastic Macrodispersion Approach: This approach aims at significantly reducing the computational cost in practical implementation of stochastic capture zone models. The approach involves modeling groundwater flow and performing random-walk based reverse particle tracking to delineate well capture zones. Random-walk based particle tracking allows incorporating the effect of macrodispersion. The longitudinal macrodispersivity is estimated from the following equation [Gelhar and Axness, 1983; Gelhar et al., 1992]: AL = σ 2 λ lnK (2.2) where λ is the scale of small-scale variability and obtained from the variaogram of fluctuation of 2 log hydraulic conductivity, σ lnK represents the variance of small-scale heterogeneity in log 79 hydraulic conductivity. Transverse dispersivity is assumed to be one third of the longitudinal dispersivity [EPA, 1986]. The macrodispersion approach was originally developed for the forward solution of solute transport equation in the presence of random, small-scale heterogeneity. We extend the approach to model capture zone “macrodispersion” due to random, small-scale heterogeneity in hydraulic conductivity. We are interested in testing if macrodispersion approach can approximately reproduce the ensemble of likely capture zone realizations from Monte Carlo simulation. For all three modeling approaches, we begin with regional scale flow modeling. A regional scale model simulates large-scale processes controlled by major streams, lakes, complex topographies and stratigraphies and provides boundary conditions for local scale models. These three approaches discussed above are summarized graphically in Figure 2.1 80 GIS DATA Water well records, Digital elevation model, National Hydrologic datasets, Glacial and Bedrock Geology REGIONAL SCALE MODEL (Simulating regional flow process and providing boundary conditions to local scale model) LOCAL SCALE MODEL (Simulating local flow process and performing reverse particle tracking) DETERMINISTIC APPROACH (Model large scale heterogeneity as a deterministic trend and ignore small scale heterogeneity) Kriging to interpolate large scale hydraulic conductivity Groundwater flow modeling Reverse particle tracking Regional variogram modeling Kriging to interpolate large scale lnK Kriging to interpolate large scale lnK Smoothing to remove noise due to data clustering Smoothing to remove noise due to data clustering Detrending lnK data to obtain small scale fluctuation Detrending lnK data to obtain small scale fluctuation Local scale residual variogram modeling Local scale residual variogram modeling Random field generation Computing macrodispersivity from stochastic macroDispersion theory Groundwater flow modeling Well capture zone delineation STOCHASTIC MACRO DISPERSION APPROACH (Model large scale heterogeneity as deterministic trend and small scale heterogeneity implicitly using effective macrodispersion) Regional variogram modeling Reverse particle tracking Well capture zone delineation Realization Loop Variogram modeling STOCHASTIC MONTE CARLO APPROACH (Model large scale heterogeneity as a deterministic trend and small scale heterogeneity as a random field) Reverse particle tracking and random walk Well capture zone delineation Visualization Statistical/probabilistic analysis and visualization Visualization Figure 2.1: Approaches used to delineate well capture zones 81 2.5 APPLICATION EXAMPLES In this section, we test the utility of the statewide database and modeling system to determine probabilistic well capture zone through a systematic case study of three real sites. Our sites are selected based on data availability, data resolution, and complexity. Figure 2.2 presents a map showing the sites, key hydrological features, pumping wells to be delineated, data distributions from water well records. We investigate the following questions in this evaluation: 1. How useful is the statewide database for practical stochastic groundwater modeling? 2. Can the relatively crude statewide conductivity dataset be used to characterize aquifer heterogeneity with sufficient details to enable probabilistic wellhead delineation? 3. Can we quantify the data uncertainty? How can we separate natural variability from data noise? 4. How can the different database components [e.g., DEM, National Hydrological Datasets (NHDs), glacial and bedrock geology, and water well records) be combined to enable modeling stochastically complex groundwater systems? Our systematic evaluation proceeds through the following 4 steps: 1. Characterization of heterogeneity at regional and local scales, 2. Regional scale and local scale deterministic flow modeling, 3. Capture zone modeling at local scale using reverse particle tracking, and 4. Systematic comparative analysis. 82 Site 1- County: Van Buren, City: Mattawan Site 3- County: Washtenaw, City: Saline Well Delineated Well Delineated 0 3,500 7,000 m 0 10,000 20,000 m Site 2- County: Van Buren, City: Lawton Type 1 Wells Well Delineated Wells UST 0 150,000 300,000 m Streams Lakes City 0 8,000 16,000 m Figure 2.2: Sites located in Michigan showing GIS features such as wells, streams, lakes, city, Type 1 wells and underground storage tanks(UST) (Site1-County:VanBuren, City: Mattawan; Site2- County-Van Buren, City: Lawton; Site3- County: Washtenaw, City: Saline) 83 2.6 RESULTS AND DISCUSSION This section begins with the discussion of geostatistical site characterization using variogram modeling, at regional and local scales, followed by discussion of deterministic flow modeling at regional and local scale. Finally capture zone modeling at local scale using deterministic approach, stochastic Monte Carlo approach and stochastic macrodispersion approach are discussed and compared. Characterization of Heterogeneity Figure 2.3 presents variogram modeling at regional scale for log hydraulic conductivity. Figure 2.4 shows variogram modeling for the detrended log hydraulic conductivity. Table 2.1 presents a summary of the parameters of regional and local variogram models. Table 2.1: Geostatistical Parameters for three sites Site 1 Site 2 Site3 Variogram Information For Log Hydraulic Conductivity (Large-scale) Theoretical Model Gaussian Exponential Exponential Nugget 0.26 0.3 0.8 Range (m) 3581 6185 5796 Resolved variability 0.14 0.14 1.2 Variogram Information For Log Hydraulic Conductivity Fluctuation (Small-scale) Theoretical Model Exponential Exponential Exponential Nugget 0.11 0.17 0.43 Range (m) 120 120 229 Resolved variability 0.17 0.16 0.56 84 0.5 (a) Site 1: Log K Variogram Variogram 0.4 0.3 Sample Variogram Theoretical Model 0.2 0.1 0.0 0 2000 4000 6000 8000 10000 Lag Distance (m) Variogram 0.6 (b) Site2: Log K Variogram 0.4 Sample Variogram Theoretical Model 0.2 0.0 0 5000 10000 15000 Lag Distance (m) 3.0 Variogram (c) Site3: Log K Variogram 2.0 Sample Variogram Theoretical Model 1.0 0.0 0 5000 10000 Lag Distance (m) 15000 Figure 2.3: Variogram modeling at regional scale for log hydraulic conductivity for three sites 85 0.4 (a) Site 1: Detrended Log K Variogram Variogram 0.3 0.2 Sample Variogram Theoretical Model 0.1 0.0 0 100 200 300 400 500 600 Lag Distance (m) Variogram 0.6 (b) Site 2: Detrended Log K Variogram 0.4 0.2 Sample Variogram Theoretical Model 0.0 0 100 200 300 400 500 Lag Distance (m) Variogram 1.5 (c) Site 3: Detrended Log K Variogram 1.0 Sample Variogram Theoretical Model 0.5 0.0 0 100 200 300 400 500 600 Lag Distance (m) Figure 2.4: Variogram modeling for detrended log hydraulic conductivity for three sites 86 The results and our experience gained through the modeling process show: • Variogram modeling process for hydraulic conductivity using data from the statewide database is much more robust (statistically more meaningful) than that based on limited measured data from a site-specific hydrogeological study. • Despite data uncertainty, the variogram modeled show a clear “spatial structure”. • Variogram modeling allows identifying three-disparate scales of variability: uncorrelated data errors/noise, small-scale variability, and large scale variability. • Variogram modeling at regional scale allow separating large scale variability from small scale variability and uncorrelated data errors. • Kriging based on the large-scale variogram, with small scale variability and data noises represented as nuggets, can be used to delineate large scale trends. The large scale variability for the three sites has a characteristic length scale on the order of 3000m, reflecting variability associated with different glacial land systems. • The large-scale variation from regional Kriging can be used to detrend hydraulic conductivity data. • The residual variograms are significantly noisier than the regional variogram, but statistically meaningful structures can be detected. • Variogram modeling of detrended log conductivity data allow identifying/ quantifying smaller spatial variability, and uncorrelated data noise. • Variances identified for natural variability in the local models are between 0.17 and 0.52. This is relatively low because hydraulic conductivity analyzed is depth averaged. Much of the vertical variability is already averaged out. 87 • The small scale heterogeneities identified for the three sites have a scale on the order of 100 to 200m. • Data noise or nugget at local scale is high, accounting for up to 40-50% of the total data variability. • The high nugget represents data noise caused by location uncertainty associated with approximate address matching, lithologic description uncertainty, well depth variability, and reporting errors. • Uncertainty in the assumed typical values is another important source of uncertainty and is not modeled in this study; • DNRE’s new, web-based water well data collection system with GPS-based input, realtime quality control and validation check, should significantly improve our practical ability to characterize aquifer heterogeneity in the future. Figure 2.5 shows the kriging of conductivity distribution at regional scale and Figure 2.6 shows a realization of large-scale and detrended log conductivity at local scale. 88 (a) Site 1 K (m/day) 13-20 21-25 26-31 32-38 39-50 0 3,500 7,000 m (b) Site 2 K (m/day) 10-18 19-23 24-28 29-33 34-41 0 8,000 16,000 m Figure 2.5: Kriged hydraulic conductivity at regional scale. The outer boundary is the extent of regional model and inner boundary is the location of local model. 89 (c) Site 3 K (m/day) 0-4 5-8 9-12 13-16 17-23 0 10,000 20,000 m Figure 2.5 (cont’d). 90 Large-scale K (a) Site 1 K (m/day) 15-18 19-21 22-24 25-27 28-31 0 1,000 2,000 m Small-scale K K (m/day) 15-18 19-21 22-24 25-27 28-31 0 1,000 2,000 m Figure 2.6: A realization of residual hydraulic conductivity (Bottom) and large-scale hydraulic conductivity (Top) for the local scale model area 91 Large-scale K (b) Site 2 K (m/day) 20-22 23-25 26-28 29-32 33-36 0 1,500 3,000 m Small-scale K K (m/day) 20-22 23-25 26-28 29-32 33-36 0 1,500 3,000 m Figure 2.6 (cont’d). 92 Large-scale K (c) Site 3 K (m/day) 6-10 11-12 13-13 14-15 16-19 0 2,000 4,000 m Small-scale K K (m/day) 6-10 11-12 13-13 14-15 16-19 0 2,000 4,000 m Figure 2.6 (cont’d). 93 A few comments are in order: • Mean conductivity at regional scale shows a spatial pattern that is much more complex than what most site-specific hydrogeological study can provide. • For a particular site, the simulated complex pattern is a function of (or is controlled by) only a few assigned typical conductivity values (e.g., for sand, silt, clay, gravel) • Interpolated or simulated conductivity distribution does not go through data values because of the significant data uncertainty or high nugget. • Large scale heterogeneity in conductivity is represented as trends deterministically in the model. • The fluctuations around the trends are represented probabilistically in the local models. Regional and Local Deterministic Flow Modeling Figures 2.7 -2.9 show the simulated head distribution from the regional and local flow models for different sites. Figures 2.10 show a comparison of simulated head and static water levels from water wells in the regional model. Table 2.1 shows key inputs and numerical parameters used in local and regional scale models. The regional scale models are simulated deterministically using no flow boundary conditions. The grid sizes and other input parameters are listed in Table 2.2. 94 Table 2.2: Parameters definitions for three sites Site1 Site2 Site3 Boundary Condition of Regional Model No flow No flow No flow Grid Design of Regional Model (Nx, Ny) 359,287 350,254 328,330 Prescribed Prescribed Prescribed Head Head Head Grid Design of Local Model (Nx, Ny) 233,181 343,253 262,256 Effective Porosity 0.15 0.15 0.15 1,2,5,10 1,2,5,10 1,2,5,10 No. of Pumping Wells 2 2 2 Pumping Capacity (GPM) 700,400 950, 750 700,525 No. of particles released in each well 2000 2000 2000 Boundary Condition of Local model Leakance (1/day) for first, second, third, fourth order streams The results and our experience gained in the modeling process shows: • Hierarchical modeling provides a systematic framework to incorporate large scale processes in local scale modeling and to allow making maximum use of data available at different spatial scales. • Data intensive modeling allows capturing dominant large scale groundwater processes largely controlled by surface water systems, topographies, and natural recharge. The high resolution DEM and NHDs and other derived GIS layers (e.g., estimated recharge) significantly improve our practical ability to model complex groundwater systems. 95 • Even without calibration, simulated heads from data intensive modeling often compares very well with observed static water levels from water well records, especially in areas glacial layer is not too thick. • The significant spread in the model-data comparison does not mean that the simulated heads are not accurate but that data being compared to is noisy. The spread is the result of significant data uncertainty reflecting temporal variability, vertical variability, location errors, and “driller variability”. Comparison between the model and static water levels from water well records must be interpreted with care. 96 Local Model Head (m) 243 - 245 246 - 247 248 - 248 249 - 250 251 - 252 253 - 254 255 - 255 0 1,000 256 - 257 2,000 m 258 - 258 259 - 261 Head (m) 211 - 217 Legend 218 - 223 224 - 228 Type 1 Wells 229 - 233 Wells 234 - 239 240 - 244 UST 245 - 250 Streams 251 - 256 Lakes 257 - 261 0 3,500 7,000 m Regional Model City 262 - 265 Figure 2.7: Site1: Head field for regional model (bottom) and local model (top) 97 Parent Model Local Model Head (m) Legend 233 - 236 Type 1 Wells 237 - 239 240 - 241 Wells 242 - 244 UST 245 - 247 Streams 248 - 249 Lakes 250 - 252 253 - 255 City 0 1,500 3,000 m 256 - 257 258 - 260 Head (m) 202 - 209 210 - 215 216 - 222 223 - 227 228 - 233 234 - 241 242 - 249 250 - 255 256 - 260 261 - 265 0 8,000 16,000 m Regional Model Figure 2.8: Site2: Head field for regional model (bottom) and local model (top) 98 Local Model Head (m) 222 - 231 232 - 236 237 - 240 241 - 244 245 - 248 249 - 252 253 - 256 257 - 260 261 - 265 0 2,000 4,000 m 266 - 271 Legend Head (m) 189 - 207 Type 1 Wells 208 - 216 Wells 217 - 224 UST 225 - 233 Streams 234 - 241 242 - 250 251 - 258 Lakes City 259 - 266 267 - 274 0 10,000 20,000 Regional Model m 275 - 284 Figure 2.9: Site3: Head field for regional model (bottom) and local model (top) 99 Predicted Head (m) (a) Site 1 Observed Head (m) Figure 2.10: Comparison of predicted head and observed head (Static water level) for three sites from water wells in the regional model. Dots: Static water level from well data, Triangles: Mean static water level from well data 100 Predicted Head (m) (b) Site 2 Observed Head (m) Figure 2.10 (cont’d). 101 Predicted Head (m) (c) Site 3 Observed Head (m) Figure 2.10 (cont’d). 102 Capture Zone Modeling: Finally, Figures 2.11- 2.13 show well capture zone delineation by 3 different methods – the deterministic, stochastic Monte Carlo and stochastic macrodisperion. Deterministic method results a single envelope for the extent of capture zones. The capture zone boundary is aligned with mean flow. For the sites simulated, the predicted capture zone is relatively narrow because of low pumping and strong regional flow. Unmodeled small scale processes can potentially causes significant errors in the capture zone distribution. The stochastic macrodispersion approach is similar to the deterministic, resulting in a single envelope for the simulated capture zone. But the approach allows modeling the dispersive effects of small-scale heterogeneity. It is interesting to notice that capture zone from stochastic macrodispersion modeling is noticeably larger than zone of uncertainty from deterministic modeling. The extent of the capture zone is lengthier, wider, and more conservative. The difference between the deterministic and macrodispersion-based delineation is the direct result of unmodeled small-scale processes and is controlled by the spatial correlation structure (scale, variance, and nugget) of the detrended hydraulic conductivity. 103 (a) Deterministic Approach (Grey color shows 10 years capture zone) 0 1,000 (b)Stochastic Macrodispersion Approach (Grey color shows 10 years capture zone) 2,000 m 0 1,000 2,000 m Legend (b)Stochastic Monte Carlo Approach (Grey colors show probabilities of 10 years ensemble capture zone) Type 1 Wells Probability (%) Wells Underground Storage Tanks City < 20 20 < − ≤ 40 40 < − ≤ 60 60 < − ≤ 80 > 80 0 1,000 2,000 m Figure 2.11: Site 1: 10 years capture zone from three approaches. 104 (a) Deterministic Approach (Grey color shows 10 years capture zone) (b)Stochastic Macrodispersion Approach (Grey color shows 10 years capture zone) 0 0 1,500 3,000 m 1,500 3,000 m Legend (b)Stochastic Monte Carlo Approach (Grey colors show probabilities of 10 years ensemble capture zone) Type 1 Wells Wells Underground Storage Tanks City Probability (%) < 20 20 < − ≤ 40 40 < − ≤ 60 60 < − ≤ 80 > 80 0 1,500 3,000 m Figure 2.12: Site 2: 10 years capture zone from three approaches. 105 (a) Deterministic Approach (Grey color shows 10 years capture zone) 0 2,000 4,000 m (b)Stochastic Macrodispersion Approach (Grey color shows 10 years capture zone) 0 2,000 4,000 m (b)Stochastic Monte Carlo Approach (Grey colors show probabilities of 10 years ensemble capture zone) Legend Type 1 Wells Wells Underground Storage Tanks City Probability (%) < 20 20 < − ≤ 40 40 < − ≤ 60 60 < − ≤ 80 > 80 0 2,000 4,000 m Figure 2.13: Site 3: 10 years capture zone from three approaches. 106 Although macrodispersion concept is widely applied in the forward modeling of contaminant transport in heterogeneous media, we are not aware of any application in the context of probabilistic wellhead protection area modeling. Stochastic Monte Carlo gives an ensemble of capture zones. Simulation of a large number of likely capture zone realizations allows computing the likelihood of being captured for a water particle in a discrete modeling cell. The occurrence of a certain location to fall in the capture zone is presented by probability map. In Figure 2.11-2.13, dark grey color represents the area that will “certainly” be captured in 10 years, despite data uncertainty. The light grey color represents the area that is unlikely to be captured with a probability of less than 20%. The other colors represent areas that may be captured, with different levels of probability. We can clearly see how uncertainty in hydraulic conductivity results in uncertainty in the boundaries of the capture zone. Stochastic modeling with its view of complex natural processes as spatial random variables is an ideal basis to build a risk assessment approach. Such an approach enables the concept of risk to be based on the principles of uncertainty that relate to the variability of the media. Thus uncertainty as to whether a point lies within the true capture zone can be estimated. Planning and protection considerations can be implemented according to the level of probability that a particular location lies within a capture zone in combination with the associated risk of the site polluting the aquifer. Stochastic modeling therefore represents a more reasonable solution for protecting groundwater supplies and ensuring sustainable future resources. 107 Limitations and Future Work This case study has also identified a number of limitations in the database and modeling system. The most significant relevant to stochastic modeling are: The stochastic models presented are two dimensional. Aquifer heterogeneity delineated is depth averaged. Future research should focus on delineation of 3D heterogeneity and 3D stochastic modeling, taking advantage of the highly valuable, high resolution 3D lithologies available in the rapidly expanding water well database. The statewide conductivity dataset was derived based on assumed typical values for each lithologic facies. These typical values can be uncertain. Future research should consider calibrating these typical values to observed static water levels. 2.7 SUMMARY AND CONCLUSIONS The incomplete knowledge of essential parameters and limited data availability results uncertainty in the prediction of well capture zones. In particular, the location of data points is often restricted to specific regions within the aquifer due to economic and logistic reasons. Furthermore, experimental data are always corrupted to some degree by measurement and by interpretive errors. Consequently, the location of the protection zones can only be defined in statistical manner and should therefore be represented using a probabilistic capture zone map. In this study, we have demonstrated the usefulness of the Michigan statewide groundwater database in stochastic simulation of well capture zones. We have compared results of capture zone delineation using the traditional deterministic method, macrodispersion method, Monte Carlo simulation method. 108 The following findings emerge from this research: • A significant portion of the variability in the statewide conductivity dataset is uncorrelated noise and this data noise can be quantified through variogram modeling. • Despite the data noise, two scales of natural variability in conductivity can be identified; • Kriging based on the large-scale variograms, with small scale variability and data noises represented as nuggets, can be used to delineate large scale trends. • Variogram modeling of detrended conductivity data can be used to identify smaller scale heterogeneity. The residual variograms are significantly noisier, but statistically meaningful structures can be detected; • Deterministic delineation that accounts only for the effect of large scale trend leads to capture zones that are significantly smaller than its stochastic counterpart, especially in areas where regional gradient is relatively strong; • Stochastic Monte Carlo simulation provides a map of well capture probability. Combined with the statewide database, MC provides a useful, systematic tool for risk-based decision making; • Stochastic macrodispersion method provides a computationally efficient alternative to delineate wellhead protection area in a way that accounts for the effect of both small and large scale heterogeneity. Future research should focus on 3D stochastic modeling using 3D lithology and calibrating the conductivity to statewide static water level dataset. 109 REFERENCES 110 REFERENCES 1. Bair, E. S., C. M. Safreed, and E. A. Stasny (1991), A Monte Carlo-based approach for determining travel time-related capture zones of wells using convex hulls as confidence regions, Ground Water, 29, 849-855. 2. Cole, B. E., and S. E. Silliman (1996), Capture zones for passive wells in heterogeneous unconfined aquifers, Ground Water, 35(1), 92– 98. 3. Dagan, G. (1989), Flow and transport in porous formations, Springer-Verlag, Berlin. 4. Dagan, G. (2002), An overview of stochastic modeling of groundwater flow and transport: from theory to applications, EOS, 83(53). 5. Dagan, G. (2004), On application of stochastic modeling of groundwater flow and transport, Stochastic Environmental Research and Risk assessment, 18(4),266-267. 6. Delhomme, J.P. (1978), Kriging in the hydrosciences, Adv. Water Resour., 1(5), 251-266 7. Delhomme, J.P. (1979), Spatial variability and uncertainty in groundwater flow parameters: A geostatistical approach, Water Resour. Res., 15, 269-280. 8. Evers, S., and D. N. Lerner (1998), How uncertain is our estimate of a wellhead protection zone? Ground Water, 36, 49-57. 9. Feyen, L., K.J. Beven, F.D. Smedt, J. Freer (2001), Stochastic capture zones delineation within the glue-methodology: conditioning on head observations, Water Resour. Res., 37 (3), 625–638 10. Feyen, L., P.J. Ribeiro Jr., F. De Smedt, P.J. Diggle (2003), Stochastic delineation of capture zones: classical versus Bayesian approach, Journal of Hydrol., 281, 313-324 11. Franzetti, S., and A. Guadagnini (1996), Probabilistic estimation of well catchments in heterogeneous aquifers, J. Hydrol., 174, 149-171. 12. Freeze, R.A., and J.A. Cherry (1979), Groundwater, Englewood Cliffs, New Jersey, Prentice-Hall. 13. Gelhar, L.W. (1993), Stochastic subsurface hydrology, Englewood Cliffs, NJ: PrenticeHall. 14. Gelhar, L.W., and C.L. Axness (1983), Three-dimensional stochastic analysis of macrodispersion in aquifers, Water Resour. Res., 19 (1), 161–180. 111 15. Gelhar, L.W., C.Welty, and K.R. Rehfeldt (1992), A critical review of data on field-scale dispersion in aquifers, Water Resour. Res., 28(7), 1955-1974. 16. Ginn, T.R. (2004), On the application of stochastic approaches in hydrogeology, Stochastic Environmental Research and Risk assessment, 18(4), 282-284. 17. Guadagnini, A., and S. Franzetti (1999), Time-related capture zones for contaminants in randomly heterogeneous formations, Ground Water, 37(2), 253– 260. 18. Guadagnini, A., U. Maione, A. Martinoli, and M.G.Tanda (1992), Pollutant Extraction by Pumping Wells: Sensitivity Analysis and Influence of Hydrodinamic Dispersion. Excerpta, 6,167-184. 19. Harrar,W. G., T. O. Sonnenborg, H.J. Henriksen (2003),Capture zone, travel time, and solute-transport predictions using inverse modeling and different geological models, Hydrology Journal, 11, 536-548. 20. Hendricks Franssen, H. J., F. Stauffer, and W. Kinzelbach (2004b), Influence of uncertainty of mean trasmissivity, trasmissivity variogram and boundary conditions on estimation of well capture zones, GeoEnv IV-Geostatistics for Environmental Applications, 13(4) 223–234. 21. Johanson, M.G. (1992), delineation of time related capture zones with estimates of uncertainty using conditional simulation of hydraulic conductivity and numerical modeling, M.S.thesis, Department of Geology, University of New Orleans, Louisiana 22. Kinzelbach, W., S. Vassolo, and G.M. Li (1996), Determination of capture zones of wells by Monte Carlo simulation, Proceedings of the Model CARE 96 Conference, IAHS Publ., no. 237. 23. Li, S.G., and D.B. McLaughlin (1993), Reply to a discussion on A Stochastic Theory for Irregular Stream Modeling, Part 1: Flow Resistance, Journal of Hydraulic Engineering, ASCE, 119, (11). 24. Li, S.G., H.S. Liao, S. Afshari, M. Oztan, H. Abbas, R. Mandle (2009), A GIS-enabled Hierarchical Modeling Patch Dynamics Paradigm for Modeling Complex Groundwater Systems Across Multiple Scales. Invited Book Chapter. In Modeling of Pollutants in Complex Environmental Systems. Grady Hanrahan (Ed), ILM Publications. In press. 25. Li, S.G., and D.B. McLaughlin (1995), Using the nonstationary spectral method to analyze flow through heterogeneous trending media, Water Resour. Res., 31(3), 541–551. 26. Li, S. G., and Q. Liu (2006), A real-time, Computational Steering Environment for Integrated Groundwater Modeling, Journal of Ground Water, 44(5),758-763. 27. Li, S.G., and Q. Liu (2008), A New Paradigm for Groundwater Modeling, Studies in 112 Computational Intelligence (SCI) 79, 19–41. 28. Li, S.G., H.S. Liao, and C.F. Ni (2004), Stochastic modeling of complex nonstationary groundwater systems, Adv. Water Resour., 27(11), 1087-1104. 29. MDEQ (Michigan Department of Environmental Quality), (2006), Wellogic System. Available at http://www.deq.state.mi.us/wellogic/mail.html. 30. Michigan Map Image Viewer, http://www.rsgis.msu.edu/mapimage/index.htm 31. Molz, F. (2004), A rational role for stochastic concepts in subsurface hydrology: a personal perspective, Stochastic Environmental Research and Risk assessment, 18(4), 278-279. 32. Molz, F.J, C.L. Dinwiddie, J.L. Wilson(2003), A physical basis for calculating instrument spatial weighting functions in homogeneous systems, Water Resour. Res., 39, 1-4. 33. Neuman, S.P. (2004), Stochastic groundwater models in practice, Stochastic Environmental Research and Risk Assessment,18(4),268-270. 34. Oztan, M. (2010), GIS-Enabled Modeling of Michigan's Groundwater Systems, PhD Thesis, Michigan State University, East Lansing, MI. 35. Rajaram, H., and D. McLaughlin (1990), Identification of large-scale spatial trends in hydrologic data, Water Resour. Res., 26(10), 2411-2423. 36. Riva, M., A. Guadagnini, and F. Ballio (1999), Time-related capture zones for radial flow in two dimensional randomly heterogeneous media, Stochastic Environmental Research and Risk Assessment, 13(3), 217–230. 37. Riva, M., L.Guadagnini, A. Guadagnini, T. Ptak, E. Martac (2006), Probabilistic study of well capture zones distribution at the Lauswiesen field site, Journal of Contaminant Hydrology, 88, 92-118. 38. Rubin, Y. (2003), Applied Stochastic Hydrogeology, Oxford Univ. Press, New York. 39. Rubin, Y. (2004), Stochastic hydrogeology-challenges and misconceptions, Stochastic Environmental Research and Risk assessment, 18(4), 280-281. 40. Simard, A. (2007), Predicting groundwater flow and transport using Michigan’s Statewide Wellogic Database, PhD Thesis, Michigan State University, East Lansing, MI. 41. State of Michigan Act 148 (2006), Groundwater Inventory and Map (GWIM) Project. Technical Report. 42. Stauffer, F., A. Guadagnini, A. Butler, H.J.H. Franssen, N.V. de. Wiel, M. Bakr, M. 113 Riva, L. Guadagnini (2005), Delineation of source protection zones Using statistical methods, Water Resour. Manage., 19,163–185. 43. Stauffer, F., S. Attinger, S. Zimmermann, and W. Kinzelbach (2002), Uncertainty estimation of well catchments in heterogeneous aquifers, Water Resour. Res., 38(11), 1238, doi:10.1029/2001WR000819. 44. Theodossiou, N., P. Latinopoulos and E. Fotopoulou (2005), Application of Monte Carlo analysis in the delineation of well head protection areas, Proceedings of the 9th International Conference on Environmental Science and Technology, Rhodes Island, Greece, 1453-1458 45. U.S. Environmental Protection Agency (1986), Background Document for the GroundWater Screening Procedure to Support 40 CFR Part 269 — Land Disposal, EPA/530SW-86-047. 46. Uffink, J.M.G. (1990), Analysis of Dispersion by the Random Walk Method, Ph. D. Thesis, Technical University of Delft. 47. van der Hoek, C.J. (1992), Contamination of a well in a uniform background flow. Stochastic Hydrology and Hydraulics, 6, 191-208. 48. van Herwaarden, O.A. (1994), Spread of Pollution by Dispersive Groundwater Flow. SIAM J. Appl. Math., 54(1), 26-41. 49. van Herwaarden, O.A., and J. Grasman (1991), Dispersive Groundwater Flow and Pollution, Math.Model. and Meth. Appl. Sci., 1, 61-81. 50. van Kooten, J.J.A. (1994), Groundwater contaminant transport including adsorption ,and first order decay. Stochastic Hydrology and Hydraulics, 8,185-205. 51. van Kooten, J.J.A. (1995), An Asymptotic Method for Predicting the Contamination of a Pumping Well, Adv. Water Resour., 18(5), 295-313. 52. van Leeuwen, M. (2000), Stochastic determination of Well Capture Zones Conditioned on Transmissivity Data, Ph.D. thesis, Imperial College of Science, Technology and Medicine, London. 53. van Leeuwen, M., C.B.M. te Stroet, A.P. Butler, J.A.Tompkins (1998), Stochastic determination of well capture zones, Water Resour. Res., 34, 2215-2223 54. van Leeuwen, M., C.B.M. te Stroet, A.P. Butler, J.A.Tompkins (1999), Stochastic determination of the Wierden (Netharlands) capture zones, Ground water, 37(1),8-17 55. van Leeuwen, M., C.B.M. te Stroet, A.P. Butler, J.A.Tompkins (2000), Stochastic determination of well capture zones conditioned on regular grids of transmissivity 114 measurements, Water Resour. Res., 36 (4), 949–957. 56. Varljen, M.D., J.M. Shafer (1991), Assessment of uncertainty in time-related capture zones using conditional simulation of hydraulic conductivity, Ground Water, 29, 737– 748. 57. Vassolo, S. W. Kinzelbach, W. Schaefer (1998), Determination of a well head zone by stochastic inverse modeling, Journal of Hydrol., 206,268–280. 58. Wheater, H.S, J.A. Tompkins, M. van Leeuwen, A.P. Butler (1998), Uncertainty in groundwater flow and transport modeling – a stochastic analysis of well –protection zones, Hydrogeological processes, 14, 2019-2029. 59. Winter, C.L. (2004), Stochastic hydrology: Practical alternative exist, Stochastic Environmental Research and Risk Assessment, 18(4), 271-273. 60. Zhang, Y.-K., and D. Zhang (2004), Forum: The state of stochastic hydrology, Stochastic Environmental Research and Risk assessment, 18(4), 265. 61. Zheng, C., and S.M. Gorelick (2003), Analysis of solute transport in flow fields influenced by preferential flow paths at the decimeter scale, Ground Water, 41, 142–155. 115