1 ram ., . .31....92' .53.... 35.3! , . .: ... z: k... . .2. ...W3m.:.v~.. I .- umrz. ”Ma... 23% 3. in... :m.b.n.~....fi.. . .. 3'31 ‘!~.t..b€. I. 16- 5 .I. $13.): )) 9a a: no . . .3... 25.5... v .. . .fi. .anufemuwfinuu .52 .d I .. 2.... ....n1. . .. ‘ l 3 3:0. 59”“: 5.2.3.. . 35...: . 43...... n P ‘ git ‘ t... 91.33....mu r . .(o;s):1.tzgft . 3". 3 .1... it: 3.2%. 2:8,... I..." {vulva w i «5.... .3... “It... . 3 . 3.35:2! » l ’35.: z... 95. 51.23.94! .3953... 1‘12"!!! if. . i .535) I .I ,I’Il I... V 1.3..... vs 2:. .1. n .15??? Hill... , ‘99....) n 1?. I 8.888 9.60? This is to certify that the dissertation entitled STOCHASTIC MODELING OF COMPLEX NONSTATIONARY GROUNDWATER FLOW SYSTEMS presented by Chuen-Fa Ni has been accepted towards fulfillment of the requirements for the PhD. degree in Civil Engineering [:3 ) L3“? \thor PWture / Decemée r. 2,6; 2.00.5?— Date MSU is an Affirmative Action/Equal Opportunity Institution LIBRARY Michigan State University PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE JUN l 9 2007 DEG 39132200? 1‘ 2/05 p:/ClRC/DateDue.indd-p.1 STOCHASTIC MODELING OF COMPLEX NONSTATIONARY GROUNDWATER FLOW SYSTEMS By Chuen-Fa Ni A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Civil and Environmental Engineering 2005 ABSTRACT STOCHASTIC MODELING OF COMPLEX NONSTATIONARY GROUNDWATER FLOW SYSTEMS By Chuen-Fa Ni Although the ”stochastic revolution” has produced an enormous number of theo- retical publications and influenced significantly how we think about the aquifer het- erogeneity, it has had relatively little impact on practical modeling community. A number of recent review articles provided the following reasons and the potential solutions: (1) Stochastic modeling is incompatible with the available measurement technologies. New measurement technologies, and new sources of data of much better resolution to characterize aquifer heterogeneity are urgently needed[9, 10, 59, 87, 57]. (2) Stochastic analytical theories are based on too many restrictive requirements to be practically useful. The assumptions of stationarity, ergodicity, mean uniform flow, gaussian distribution, and small perturbation must be substantially relaxed [9, 47, 42, 43, 62, 59, 87, 19]. (3) Stochastic numerical theories are computationally impractical for most problems of realistic sizes. One must recognize and remove these tough computational bottlenecks before meaningful stochastic modeling applications are possible [9, 59]. Motivated by these critical assessments, this research addressed a number of con- ceptual, computational, and implementation issues in the modeling of subsurface heterogeneity. In particular, this study deveIOped, tested, and implemented the con- ceptually improved, nonstationary stochastic methods for predicting velocity uncer- tainty in two-dimensional flows. An approximate and analytically-based spectral method was presented for predicting velocity variances under mildly nonstationary flow situations. This approximate spectral method (ASM) rely on a linearization of the groundwater flow equation but do not require the common statistical stationar- ity assumptions. To provide general insight into the ASM for mildly nonstationary flows, this study performed intensive numerical experiments to assess the accuracy of ASM under a number of nonstationary situations. The illustrative examples involved nonstationary situations caused by hydraulic conductivity trends, composite media, nonlinear head distributions in unconfined aquifers, transient flows, and deterministic sources and sinks applied in modeling areas. The surprising results in the assessment showed that ASM can reproduce well the solutions of corresponding first-order nu- merical analysis and Monte Carlo simulation except in the proximity of prescribed boundaries and well locations. These regions, however, are limited in 3 to 5As (lnk correlation length) from boundaries and in 5 to 10As from wells. Due to the inherent limitations of the analytically-based ASM, the detail dynam- ics of strongly nonstationary regions such as prescribed head boundaries and strong stresses are not well described. This study further presented a hybrid spectral method (HSM) to predict velocity variances in strongly nonstationary flows. The proposed hybrid method, based on solving stochastic perturbation equations, involves two ma- jor computational steps after solving the mean flow equation. The first step applies the analytically-based ASM to predict the nonstationary variances for the entire mod- eling area. Then the second step employs numerically-based nonstationary spectral method (NSM) to correct the ”regional solution” in localized areas where the vari- ance distribution is considered to be strongly nonstationary. The boundary conditions for the localized numerical solutions are adopted from reginal ASM solutions. This study then illustrated HSM with two steady, two-dimensional, and complex nonsta- tionary flow problems. The results showed that the proposed hybrid method, which takes major advantages of analytical and numerical techniques, can efficiently handle large modeling areas and can accurately predict the detail dynamics around highly nonstationary locations. ACKNOWLEDGMENTS I would like to express my gratitude to my advisor, Dr. Shu—Guang Li, for his support, patience, and encouragement throughout my studies. His technical and editorial advice was essential to the completion of this dissertation and has taught me innumerable lessons and insights on the workings of academic research in general. My thanks also go to the members of my research committee, Dr. Roger Wallace, Dr. Mantha Phanikumar, and Dr. David Hyndman for reading previous drafts of this dissertation and providing many valuable comments that improved the presentation and contents of this dissertation. I am also grateful to my other technical advisors Dr. Huasheng Liao and Dr. Qun Liu. They helped me considerably with programming of FORTRAN and Visual Basic language. The friendship of research team members, Soheil, Andreanne, Dr. Lu, Prasanna, Dr. Lee, and Kyle is much appreciated and has led to many interesting and good-spirited discussions relating to this research. Last, but not least, I would like to thank my wife Hsiao-Ching for her under- standing and love during the past few years. Her support and encouragement was in the end what made this dissertation possible. My parents, my mother-in-law, my brothers, and my sister receive my deepest gratitude and love for their support during my studies. iv TABLE OF CONTENTS LIST OF TABLES vii LIST OF FIGURES viii 1 Introduction 1 1.1 Stochastic Analytic Solutions and Limitations ............. 2 1.2 Stochastic Numerical Solutions and Computational Challenges . . . . 4 1.3 Research Objective and Basic Approaches ............... 10 1.4 Organization of the Dissertation ..................... 12 2 Simple Closed-Form Formulas for Predicting Groundwater Flow Model Uncertainty in Complex, Heterogeneous Trending Media 14 2.1 Abstract .................................. 14 2.2 Introduction ................................ 15 2.3 Problem Formulation ........................... 16 2.4 Spectral Solution ............................. 17 2.5 Approximate Closed Form Solution ................... 18 2.6 Illustrative Examples ........................... 20 2.7 Results ................................... 22 2.8 Conclusions ................................ 24 3 Modeling Groundwater Velocity Uncertainty in Nonstationary Com- posite Porous Media 28 3.1 Abstract .................................. 28 3.2 Introduction ................................ 29 3.3 Problem Formulation ........................... 31 3.4 Spectral Solution in Composite Media ................. 32 3.5 Approximate Spectral Method ...................... 34 3.6 Illustration Examples ........................... 36 3.7 Results and Discussion .......................... 39 3.8 Conclusions ................................ 46 4 Quantifying Flow Uncertainty in Stochastic Groundwater Models 56 4.1 Abstract ............................ , ..... 56 4.2 Introduction ................................ 57 4.3 Problem Formulation ........................... 58 4.3.1 Approximate Perturbation Equations .............. 60 4.3.2 Spectral Solutions (Approximate Spectral Method) ...... 61 4.4 Illustrative Examples and Numerical Consideration .......... 62 4.5 Results and Discussion .......................... 63 4.5.1 Steady state situation ...................... 63 4.5.2 Transient state situation ..................... 67 4.6 Conclusion ................................. 68 5 Prediction of Velocity Uncertainty in Complex, Nonstationary, and Heterogeneous Porous Media - A Hybrid Spectral Method 82 5.1 Abstract .................................. 82 5.2 Introduction ................................ 83 5.3 Problem Formulation ........................... 86 5.3.1 The Nonstationary Spectral Method (NSM) .......... 87 5.3.2 Stochastic Modeling of Moderately Nonstationary Groundwa- ter Systems - Approximate Spectral Method(ASM) ...... 90 5.3.3 Stochastic Modeling of Strongly Nonstationary Groundwater Systems - Hybrid Spectral Method(HSM) ........... 91 5.4 Illustrative Examples and Numerical Consideration .......... 92 5.5 Results and Discussion .......................... 93 5.6 Conclusion ................................. 95 6 Summary 106 References 110 vi 1.1 2.1 3.1 4.1 LIST OF TABLES Some selected stochastic groundwater studies .............. 7 Parameter definitions for the examples ................. 22 Parameter definitions for three test examples .............. 38 The definitions of physical and first-order numerical parameters for six examples .................................. 64 vii 1.1 1.2 1.3 1.4 2.1 2.2 2.3 2.4 LIST OF FIGURES The solution types for different stochastic methods, (a) Monte Carlo method solves a series of realizations, (b) Moment equation method solves the cross covariance function over the modeling domain, (c) Green’s function method solves the Green’s function equation over the modeling domain, (d)Sensitivity derivative method solves the sensitiv- ity derivative function over the modeling domain. (Source: Li et al., 2003) ................................... The comparison of CPU time for different methods, where L repre- sents the domain length and /\ represents the correlation length of log hydraulic conductivity.(Source: Li et al., 2003) ............ The comparison of memory requirement for different methods, where L represents the domain length and /\ represents the correlation length of log hydraulic conductivity.(Source: Li et al., 2003) ......... The conceptual diagram of the Hybrid Spectral Method (HSM). . . . Spatial distribution of the geometric mean conductivity and corre- sponding mean head distribution for the continuous trending case. . . Spatial distribution of the geometric mean conductivity and corre- sponding mean head distribution for the discontinuous trending case. Centerline profile of the predicted velocity standard deviations using the closed-form formulae, nonstationary numerical spectral method, and Monte Carlo simulation for continuous trending case ........ Centerline profile of the predicted velocity standard deviations using the closed-form formulae, nonstationary numerical spectral method, and Monte Carlo simulation for discontinuous trending case. ..... viii 10 11 13 23 24 25 26 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 A realization of nonstationary conductivity distribution plotted along the domain centerline in the horizontal direction for a) example 1, b) example 2, c) example 3 ......................... Spatial distribution of the geometric mean conductivity and the corre- sponding mean head distribution for example 1 ............ Spatial distribution of the geometric mean conductivity and the corre- sponding mean head distribution for example 2 ............ Spatial distribution of the geometric mean conductivity and the corre- sponding mean head distribution for example 3 ............ Spatial distribution of the predicted velocity variances using the closed— form formulas, nonstationary spectral method, and Monte Carlo sim- ulation based on the simple exponential covariance model (example 1) ..................................... Spatial distribution of the predicted velocity variances using the closed- form formulas, nonstationary spectral method, and Monte Carlo sim- ulation based on the simple exponential covariance model (example 2) ...................................... Spatial distribution of the predicted velocity variances using the closed- form formulas, nonstationary spectral method, and Monte Carlo sim- ulation based on the simple exponential covariance model (example 3) ...................................... Centerline profile of the predicted velocity standard deviations using the closed-form formulas, numerical spectral method, and Monte Carlo simulation based on the simple exponential covariance mode1(example 1) ...................................... Centerline profile of the predicted velocity standard deviations using the closed-form formulas, numerical spectral method, and Monte Carlo simulation based on the simple exponential covariance model(example 2) ...................................... Centerline profile of the predicted velocity standard deviations using the closed-form formulas, numerical spectral method, and Monte Carlo simulation based on the simple exponential covariance model(example 3) ...................................... ix 39 40 41 42 47 48 49 50 51 52 3.11 3.12 3.13 4.1 4.2 4.3 4.4 4.5 Centerline profile of the predicted velocity standard deviations using the closed-form formulas, numerical spectral method, and Monte Carlo simulation based on a hole-type covariance model (example 1) Centerline profile of the predicted velocity standard deviations using the closed-form formulas, numerical spectral method, and Monte Carlo simulation based on a hole-type covariance model (example 2) Centerline profile of the predicted velocity standard deviations using the closed-form formulas, numerical spectral method, and Monte Carlo simulation based on a hole-type covariance model (example 3) The conceptual model and solutions for example 1: (a) the conceptual model and the corresponding mean head distribution in an unconfined aquifer, (b) the longitudinal velocity (0“) standard deviation along domain centerline (3:2 = 40m), and (c) the transverse velocity (0,) standard deviation along domain centerline ............... The conceptual model and solutions for example 2: (a) the spatial distribution of the geometric mean conductivity and the corresponding mean head distribution in an unconfined aquifer, (b) the longitudinal velocity (on) standard deviation along domain centerline, and (c) the transverse velocity (0,) standard deviation along domain centerline The conceptual model and solutions for example 3(global recharge): (a) the conceptual model and the selected mean head distribution in an unconfined aquifer, (b) the longitudinal velocity (on) standard de- viation along domain centerline (2:2 = 40m), and (c) the transverse velocity (0,) standard deviation along domain centerline ....... The conceptual model and solutions for example 3 (local recharge): (a) the conceptual model and the selected mean head distribution in an unconfined aquifer, (b) the longitudinal velocity (on) standard de- viation along domain centerline (x2 = 40m), and (c) the transverse velocity (0”) standard deviation along domain centerline ....... The conceptual model and solutions for example 3 (well): (a) the con- ceptual model and the selected mean head distribution in an unconfined aquifer, (b) the longitudinal velocity (on) standard deviation along do- main centerline (11:2 = 40m), and (c) the transverse velocity (0”) stan- dard deviation along domain centerline ................. 53 54 55 69 70 71 72 73 4.6 4.7 4.8 4.9 4.10 4.11 5.1 5.2 5.3 5.4 5.5 The conceptual model and solutions for example 4: (a) the conceptual model and the mean head distribution, (b) the longitudinal velocity (on) standard deviation along profile A-A’, and (c) the transverse ve- locity (0,) standard deviation along profile A-A’ ............ The conceptual model for example 5. The transient flow is driven by the west boundary where the head values are changed with time based on a sinusoidal function. ......................... The modeling results for example 5: (a) the mean head time series at two monitoring wells, (b) the time series of longitudinal velocity (on) standard deviation at two monitoring wells, and (c) the time series of transverse velocity (0,) standard deviation at well locations ..... The conceptual model for example 6. In addition to the transient flow situation, a conductivity trend and multiple source/sinks are applied to the modeling area ............................ The simulation results for example 6: (a) the mean head time series at two monitoring wells, (b) the time series of longitudinal velocity (on) standard deviation at two monitoring wells, and (c) the time series of transverse velocity (0”) standard deviation at well locations ..... The solutions at t = 2.0 (day) for example 6: (a) the mean head distri- bution along profile A-A’, (b) the longitudinal velocity (on) standard deviation along profile A-A’, and (c) the transverse velocity (0,) stan- dard deviation along profile A-A’ .................... The conceptual diagram of the Hybrid Spectral Method (HSM). . . . Spatial distribution of the mean head distributions for pumping well example .................................. Spatial distribution of the geometric mean conductivity and the corre- sponding mean head distributions for the application example Centerline profiles of predicted velocity standard deviations using the approximate spectral method, numerical nonstationary spectral method, hybrid spectral method, and Monte Carlo simulation Spatial distribution of the predicted longitudinal velocity variances us- ing the approximate spectral method, nonstationary spectral method, hybrid spectral method, and Monte Carlo simulation (the correction area near well locations) ......................... xi 74 75 76 77 78 79 95 96 97 98 99 5.6 5.7 5.8 5.9 5.10 Spatial distribution of the predicted transverse velocity variances us- ing the approximate spectral method, nonstationary spectral method, hybrid spectral method, and Monte Carlo simulation (the correction area near well locations) ......................... Spatial distribution of the predicted longitudinal velocity variances us- ing the approximate spectral method, nonstationary spectral method, hybrid spectral method, and Monte Carlo simulation (the correction area near the lake location) ....................... Spatial distribution of the predicted longitudinal velocity variances us- ing the approximate spectral method, nonstationary spectral method, hybrid spectral method, and Monte Carlo simulation (the correction area near the lake location) ....................... The profiles of predicted velocity standard deviations using the ap— proximate spectral method, numerical nonstationary spectral method, hybrid spectral method, and Monte Carlo simulation (x2: 80m) . . . The profiles of predicted velocity standard deviations using the ap- proximate spectral method, numerical nonstationary spectral method, hybrid spectral method, and Monte Carlo simulation (11:2: 110m) xii 100 103 104 CHAPTER 1 Introduction It is now generally agreed that natural subsurface environment is very heterogeneous. Soil heterogeneities, in particular, cause dramatic variation in hydraulic conductivity from point to point within a groundwater formation. Such variation sometimes ap- pear to be random, although many sites also exhibit trends which are related to the sedimentation processes that create stratified deposits of contrasting soils, alluvial fans, deltas, and glacial outwash plains. At first glance it may seem that small-scale variation in hydrogeologic properties should have relatively little effect on the larger- scale movement of contaminants. In reality, small—scale fluctuations in hydraulic conductivity can have significant large—scale consequences, primarily because of the nonlinear relationship which exists between conductivity and concentration [16, 7]. Over the years a number of stochastic techniques have been proposed for ana- lyzing the role of spatial heterogeneity in groundwater flow systems [e.g.[74, 75, 17, 8, 61, 24, 95, 69, 100, 42, 43] ]. Generally, they assume that heterogeneous physical parameters such as hydraulic conductivity are random spatial functions with known statistical properties. It follows that environmental variables which depend on these parameters are also random. So, for example, uncertainty in hydraulic conductiv- ity induces uncertainty in pore velocity which, in turn, induces uncertainty in solute concentration. The flow equations provide a physical basis for relating the moments of random dependent variables (e.g. mean head and velocity, and their associated variances and covariances) to the hydrogeological parameters which are the original source of uncertainty. In practice, however, it is very difficult to derive these mo- ments without making approximations of one kind or another. The art of stochastic groundwater modeling lies in knowing which approximations are most appropriate for a given application. 1.1 Stochastic Analytic Solutions and Limitations One of the most important approximations introduced in stochastic groundwater anal- ysis is the assumption of statistical stationarity. The technical definitions of station- arity are concerned with the behavior of the statistical properties of spatial random function when its time and space origin is shifted [65]. A spatial random function is, for example, wide-sense (or weakly) stationary when its mean and covariance do not depend on absolute location. Most hydrologic and geological variables are not truly stationary since their means can vary over the scales of interest. If these mean variations are small compared to random fluctuations about the mean (i.e. if the mean is nearly constant over many correlation scales), stationarity may be assumed to hold in an approximate or “local” sense. This concept of “local stationarity” is frequently encountered in groundwater applications [17, 95]. When the concept is justified, stationarity assumptions are advantageous both methodologically and conceptually. They enable us to use a variety of analyti- cal techniques to solve stochastic flow and transport problems [see, for example, [18, 8, 95, 69]]. More importantly, they enable us to invoke the ergodic hypothesis, which establishes an equivalence between spatial heterogeneity and uncertainty. The stationarity-based stochastic theories have produced a number of insightful results on field-scale flow and transport and unified the recent stochastic and classical deter- ministic modeling framework [see reviews by [16, 7]]. In particular, the stationarity- based stochastic theories show that “randomly” heterogeneous groundwater forma- tion may be equivalently represented, at field-scale, as a more homogeneous medium characterized by a set of mean effective properties. These effective properties can be related by stochastic theories to the statistical structure of media heterogeneity [17, 8, 95, 69]. Consequently, one can model field—scale flow and transport using a conventional deterministic model as if there were no small-scale heterogeneity as long as soil properties (e.g., hydraulic conductivities) are replaced by their corresponding effective values. As shown in numerous publications dealing with field-scale flow and transport at the well known Bordon site [79, 17] and the Cape Code site [14, 17] and several synthetically-generated heterogeneous sites using large-scale numerical simu- lation [86, 17, 69], stochastic theories accurately predict effective hydraulic conduc- tivities and field-scale mean flow, longitudinal solute macrodispersion, and evolution of spatial moments of the tracer plume, given a reasonable estimate of the statistical and geostatistical parameters of the hydraulic conductivity field. However, groundwater flow at most sites are statistically nonuniform and, in some cases, strongly nonuniform in response to one or more of the following factors: 1) nonstationarity (or vertical and/or horizontal trends) in hydraulic conductivity, 2) internal and external sources/ sinks (e.g. pumping and injection wells), 3) distributed recharge, 4) geologic and hydrologic boundaries and the associated nonstationary stresses imposed, 5) systematically-varying aquifer thickness, and 6) transient mean flow effects. If stochastic modeling is to become a viable tool for real-world ground- water investigations, it must be able to accommodate these nonstationarities and complex sources and sinks in addition to just ”random” small-scale heterogeneity since they are part of almost every realistic applications. 1.2 Stochastic Numerical Solutions and Computa- tional Challenges There are a number of numerical approaches that can be used to analyze general non- stationary flow problems under complex field conditions. These include, for example, Monte Carlo methods [e.g., [74, 75, 20, 21, 61, 37, 28, 29, 95, 69]] and perturbation approaches, such as the moment equation methods [e.g.,[20, 21, 99, 96, 95, 50, 97]], Green’s function methods [e.g.,[8, 61, 83, 84, 85, 24, 25, 95, 69]], and the so called first-order second-moment methods based on Taylor’s expansions [e.g.,[80, 11, 93, 95, 47, 102, 3, 92]]. All of these methods, however, are computationally demanding when applied to flow or transport problems of realistic size [20, 21, 17, 54, 47 , 42, 43, 62, 64], mostly because they all require very fine spatial discretizations in order to resolve (ei- ther directly or indirectly) the underlying small—scale heterogeneous dynamics. The Monte Carlo method has to solve large numbers of numerically-difficult flow equations with highly fluctuating coefficients. The moment method has to solve for at least N coupled large covariance matrices that have a limiting spatial scale similar to that of small-scale heterogeneity and requires 0(N2) words of on-line memory (note the exponential growth with N l), N being the total number of discrete computational nodes. The Green’s function is obtained from a numerical solution of the linearized flow equation, with the forcing term replaced by a Dirac delta function [24, 25].The Green’s function method has to solve N Green’s function equations that are basically equivalent to the covariance equations. The major computational cost involved in the first-order second moment methods lies in the calculation of the so-called sensitivity derivatives. These coefficients can be obtained directly by the sensitivity equation method [94] or by the adjoint state method [80, 81, 28, 40]. Both methods are expen- sive when applied to flow problems. The sensitivity equation method has to solve N sensitivity derivative equations. These equations are also equivalent to the covariance equation equations. The adjoint state method has to solve for N fast-varying adjoint functions and resolve a singular Dirac delta function at every computational node. Note the usual advantage of the adjoint state method over the direct sensitivity equa- tion method is lost when the head sensitivity coefficients are required everywhere, as is the case when we need to quantify the effects of scale interactions on large-scale flow (i.e. need to calculate the closure covariances) and/ or when we need to quantify the uncertainty associated with the mean predictions throughout the computational domain. Figure 1.1 presents the solutions of Monte Carlo simulation and different perturbation approaches when solving a simple uniform flow problem [47]. The mod- eling area is 50 A by 50 A with constant head gradient 0.1 from the west boundary to the east boundary, where A is the correlation length of log hydraulic conductiv- ity. The profile in Figure 1.1 represents the centerline from the west side to the east side of the modeling area. It is clear that under the simple flow condition these methods all need fine grids over entire modeling domains to capture rapidly changing variations and transfer functions. Investigators who introduce classical perturbation approaches as alternative to Monte Carlo method, expect them to be useful for prac- tical applications. However, many researchers fail to realize that the sensitivity of the computational cost to the domain size makes classical perturbation approaches impractical for most realistic problems. The numerical difficulty in calculating the moment, Green’s function, and sensi- tivity equations was pointed out in [47, 42, 43]. The implementation difficulty in the moment equation methods were stressed by [55, 55, 20, 21, 47, 42, 43]. In order to ex- pect accurate solutions by any of these techniques, the discretization required would have to be significantly smaller than the typical scales of variation of the hydraulic A n) V A 0.10 0.08 - , oment at on at : Solving tor Nx cross oov tunes Head-00nd. Cov. Functions 3 Greens Functions Greens Function Met eth:od Solving tor Nx Greens functions g Sensitivity Equations 0.04 - 0.03 - 0.02 - 0.01 - 3: ll ’l Ll eMethod. Solving tor Nx sensitivity derivative functions l'll lll‘ L thin -0.01- -0.02- -0.03- .fl l l I "l’ lill’i -0.04 - -0.05 Figure 1.1. The solution types for different stochastic methods, (a) Monte Carlo method solves a series of realizations, (b) Moment equation method solves the cross covariance function over the modeling domain, (c) Green’s function method solves the Green’s function equation over the modeling domain, (d)Sensitivity derivative method solves the sensitivity derivative function over the modeling domain. (Source: Li et al., 2003) parameter (often order(s) of magnitude smaller than that used in a corresponding deterministic model). As a result, computational effort involved often becomes pro- hibitively expensive, and thus applications of these methods to field situations have been very limited [17 , 9, 47, 42, 43, 59]. Table 1.1 summarizes some important cases that have been demonstrated in the past research. Note that all modeling problems in Table 1.1 are restricted to a relatively small domain size, i.e, the ratio of domain size to the correlation length is very small. Table 1.1: Some selected stochastic groundwater studies. Stochastic Numerical Problem Description Domain Grid No. Source Model Scheme Size(L) MC(Nr) MM FE Flow and transport 44x24 m, 34x23 [20] problem,conditioning on A=2 m conductivity and head SM FE Flow, inverse problems, 560x240 m, 14x16 [80] conditioning on head and A2200 m Ink SM FD Cape Cod site, multi-scale 800x40 m, 120x20 [11] transport problem, A=3.2,3.5 conditioning on tracer data GM,MC FE Flow problem with pumping L—22 A, 120x20 [61] well A—2, 5, 8 N SM,MC FD Flow problem with linear L—8 A, NA [45] conductivity trend SM FE Unsaturated flow, 80x80 cm, 20x20 [39] conditioning on head A=20 cm MC FE Flow, conditioning on head 40x40 20x20, [28] A=2 Nr=400 MM FD Unsaturated flow L=10 A, NA [98] A=10 cm MC FD Saturated flow and 64x40 A, 256x400, [30] transport problem A=l Nr=2400 MC FD Transport problem 5x3x0.1 km, 100x60x50, [53] Ah=250 m, Nr=NA Av=10 m MM,MC FD Flow problem 200x100, NA, [99] A=17.5, 16.6 Nr=NA MC FE Vertical unsaturated flow 24x12.8 m, 80x128, [29] and transport A, =3 Nr=300 Az=0.5 SM FE Flow and transport 40x200 cm, 10x50 , [40] Continued on Next Page... Table 1.1 —- Continued Stochastic Numerical Problem Description Domain Grid No. Source Model Scheme Size(L) MC(Nr) conditioning on head A=40cm GM FE Flow and conditioning 18x8, 180x80 , [24] on head A=1 Nr=4000 MM FD Saturated and unsaturated 120x360, 20x60 , [96] flow A=30 MC FE Vertical unsaturated flow 20x20 m, 50x50 , [38] and transport Ay=5 m N r=50 Az =3 m MC FE Flow problem 12x12 m, 60x60 , [88] A=1 m Nr=5000 MC FE Unsaturated flow problem 4x8A, 20x40 , [49] Nr=3000 MM FE Unsaturated flow problem 40x3A, 40x30 , [82] A=40,4 MC FD Transport problem 53x25.6A, 256x128 , [32] A=1 Nr=2000 NSM,SM, FD Flow problem 50x50A, Ax=0.5-2 , [47] GM,MM, A=1 Nr=20000 MC MC FE Unsaturated flow and 120x360 cm, 20x60 , [50] transport problem A=30 cm Nr=NA TE,MC FD Unsaturated flow problem 10x10A 100x100 , [102] Nr=10000 TE,SM FE Flow inverse problem 400x200 m, 40x20 , [3] A=50 111 MC FD Transport problem 1x2x4 A 128x64x64 , [100] Nr21600 AM FE Plume travel time, 100x50 m, 1000x500 , [6] conditioning on head two zones, and Ink A=2 and 10 m GM,MC FE Flow problem, conditioning 8x4 A 40x20 , [92] on Ink Nr=2500 NSM,MC FD Flow problem 100x100 A, Ax=l,2 A, [43] 200x200 A Nr=20000 MM,MC FE Flow problem 10x10 A 40x40, [97] Nr=5000 MM,MC FE Water-oil two-phase 3x0.96 In, 50x16, [4] problem A=0.3 m Nr=2000 MM FD 'h‘ansport problem 20x10 A, 41x41, [91] A=1 m SM FE Transient Flow, 15x15x15 m, 15x15x15, [103] inverse problem A=1 m The listed cases consider solely Ink as the source of uncertainty. Continued on Next Page. . . Table 1.1 — Continued Stochastic Numerical Problem Description Domain Grid No. Source Model Scheme Size(L) MC(Nr) MM: Moment equation method, MC: Monte Carlo simulation, SM: Sensitivity equation method, GM: Green’s function method, TE: Taylor expansion method, NSM: Nonstationary spectral method, FD: Finite difference numerical scheme, FE: Finite element numerical scheme, Nr: Number of realizations in Monte Carlo simulation. NA: Information is not available. Over the past few years, Li and his co—workers have developed a new stochastic modeling approach, a nonstationary spectral method, for predicting nonstationary flows [45, 47, 42, 43]. This approach is based on an extension of the classical spectral method for stationary flow problems [2, 27, 56, 18, 23, 22, 41, 17] and has several advantages over other methods for analyzing nonstationary groundwater problems. The new approach allows for the first time modeling general field-scale processes and uncertainty without having to resolve numerically the small-scale dynamics, sig- nificantly increasing the size and expanding the range of site characterization prob- lems that can be analyzed with stochastic methods. This nonstationary spectral method combines the best features of analytical and numerical techniques. The sta- tistically stationary small-scale portion of natural variability is described with a com- pact spectral (Fourier) representation while the remaining nonstationary component is described as a larger-scale spatial process driven by mean gradients. The numeri- cal computations focus on departures from stationarity. The division of labor is not prescribed explicitly but is handled naturally with transfer functions that depend on both wave number and location. Discretized numerical computations are performed on spatial grids that are comparable in size to the grids used in traditional deter- ministic modeling applications. The results, as illustrated in Figures 1.2 and 1.3, dramatically decrease the required computational effort to compute statistics such as head and velocity variances. 100 ensitivity and reens Functio ~ A ethod (I) 30 P Monte Carlo 2 . (20,000 simulations) 3 : ._ 60 .. s : a) I Moment Equation E 40 '- ._ i Monte Carlo "‘ : (10,000 simulations) 2 20 ’ O i . Nonstationary Spectral 01 .1111111111111 1111. .11111111111. 0‘”10 20 30 40"L5/0A 60 70 80 90 100 Figure 1.2. The comparison of CPU time for different methods, where L represents the domain length and A represents the correlation length of log hydraulic conduc- tivity.(Source: Li et al., 2003) However, the nonstationary spectral method is still significantly more expensive than deterministic methods, especially for large unsteady problems, primarily because the method still requires solving large numbers of partial differential equations. 1.3 Research Objective and Basic Approaches The objective of the proposed research is to further improve the nonstationary spectral method so that it can be routinely applied in site-specific groundwater investigations. We are particularly concerned with the design of a characterization tool that can adapt to complex and general nonstationary environmental conditions in complex 10 500 From left to right: Greens Function Nonstationary Spectral Method - - - - Monte Carlo Simulation 400 t Sensitivi Derivative Moment quation A I m _ 2- 300 _- V _ 2‘ : 0 200 - E _ a) - Figure 1.3. The comparison of memory requirement for different methods, where L represents the domain length and A represents the correlation length of log hydraulic conductivity.(Source: Li et al., 2003) geometries and trending lithologies. The improved approach is based on the following observations: 0 Most groundwater systems are only weakly to moderately nonstationary ex- cept in localized areas (e.g., around wells, across disconituities, and prescribed boundaries) . 0 There is a significant scale disparity between the mean variation and the small- scale heterogeneity in these weakly or moderately nonstationary areas. 0 One can dramatically simplify the general nonstationary spectral method if the scale disparity is taken advantage of and a highly efficient, closed form solution 11 can be derived. 0 The more intensive nonstationary spectral modeling only needs to be applied in localized areas where the variance dynamics is strongly nonstationary. Given these observations, we propose in this study a hybrid spectral method (HSM) for stochastic groundwater modeling. The improved “hybrid” spectral mod- eling will proceed in two steps: First, the approximate closed form spectral method (ASM) is applied to obtain a “regional” screening level uncertainty analysis. Second, the nonstationary spectral method (NSM) is applied to locally refine or correct the solution where the variance dynamics are rapidly varying. The boundary conditions for the local nonstationary spectral solutions will be based on the “regional” closed form solution. The regional solution and the numerical corrections are both expected to be highly efficient and the computational process will be virtually instantaneous since the former is in closed form and the latter are only applied in very small areas. Figure 5.1 illustrates the concept of hybrid spectral method. 1.4 Organization of the Dissertation This dissertation will be organized in six sections including this introductory chap- ter that reviews where we are in stochastic groundwater modeling, motivates the research, and explains the source of difficulties in applying existing stochastic meth- ods to realistic field problems. Chapter 2 to 5 are written in the form of technical paper. Chapter 2 present the capability of ASM for predicting flow uncertainty in complex, heterogeneous trending media. In this paper, the nonstationary flows were driven solely by the conductivity trends. The small-scale statistical structure of the log conductivity was kept uniform in the modeling areas. In Chapter 3, the ASM was 12 Modeling Area : $13551»; flea—o ' ', 1 ' " " tak-e ‘ ' "I l. _______ I I . . I I | ______ I ASM solution , _ --. WE" I ‘‘‘‘‘‘ I o -------- _. k._- ---—I--:.:.:.='- \. I Wm (NSM): ‘I’ni (ASM) l \_ l g 31 \\ : g. NSM solution 3.1 -\ - 2' \ ‘ ? Well 9'1 . ... Il . \. ' a 0 31 '- - I Apply Nonstatlonary Spectral Method(NSM) \ I 2 g l _ _. Numerical Method ‘\ I :5: g D Apply Approximated Spectral Method(ASM) '\ '9- 5.‘ Analytical Method \. l I \ 'i _ _ _‘" m 060:1" mlosw. .I Figure 1.4. The conceptual diagram of the Hybrid Spectral Method (HSM). employed to predict the flow uncertainty for composite porous media. In the compos- ite porous media, the small-scale statistical structure of the log conductivity became nonuniform. Chapter 4 provides a systematical comparison of ASM for flow under complex flow conditions. These conditions include:(1) groundwater flow with sources and sinks in confined and unconfined aquifers, and (2) transient flow in confined and unconfined aquifer with complex trends and multiple sources and sinks. The general concept and associated applications of the HSM were demonstrated in Chapter 5. Two examples were presented to show the effectiveness and accuracy of the HSM. In the end of this dissertation, I organized the overall conclusion in Chapter 6 for this study. 13 CHAPTER 2 Simple Closed-Form Formulas for Predicting Groundwater Flow Model Uncertainty in Complex, Heterogeneous Trending Media 2.1 Abstract This paper presents approximate, closed form formulas for predicting groundwater velocity variances caused by unmodeled small-scale heterogeneity in hydraulic con- ductivity. These formulas rely on a linearization of the groundwater flow equation but do not require the common statistical stationarity assumptions. The formulas are illustrated with a two-dimensional analysis of steady state flows in complex, multi- dimensional trending media and compared with a first-order numerical analysis and Monte Carlo simulation. This comparison indicates that, despite the simplifications, the closed-form formulas capture the strongly nonstationary distributions of the ve- 14 locity variances and match well with the first order numerical model and the Monte Carlo simulation except in the immediate proximity of prescribed head boundaries and mean conductivity discontinuities. The complex trending examples illustrate that the closed-form formulas have many of the capabilities of a full stochastic numerical model while retaining the convenience of analytical results. 2.2 Introduction It is now widely recognized that hydrogeological properties such as hydraulic con- ductivity vary significantly over a wide range of spatial scales [17]. This variability reflects the aggregate effect of different geological processes acting over extended time periods. It is tempting for analytical reasons to adopt a simplified description of spatial variability which imposes some sort of regular structure on the subsurface environment. Such a description can be deterministic (e.g., the subsurface consists of homogeneous layers) or stochastic (e.g., the subsurface properties are stationary ran- dom fields). Although these simplified descriptions do not capture the true nature of hydrogeologic variability, they may provide reasonable approximations for particular applications. In certain situations, variables such as hydraulic conductivity exhibit local trends that extend over scales comparable to the overall scale of interest in a given problem. Trends are a natural result of the sedimentation processes that create deltas, alluvial fans, and glacial outwash plains [1]. In such cases it seems reasonable to represent hydrogeologic variability as a stationary random field superimposed on a deterministic trend [60, 67]. This paper is concerned with predicting the velocity variances caused by unmodeled small-scale heterogeneity in the presence of systematic trends in mean hydraulic conductivity. 15 Stochastic approaches for analyzing groundwater flow in heterogeneous trending media generally divide into (1) analytical techniques (e.g., [48, 72, 45, 61, 33]), which provide convenient closed-form expressions for quantities such as the velocity vari- ances and effective hydraulic conductivities but depend on restrictive assumptions (e.g., linear trends) and (2) numerical techniques (e.g., [74, 55, 47, 42]), which make fewer assumptions but are more difficult to apply in practice. Here we present sim- ple, closed-form analytical solutions that relax the assumptions required in existing analytical theories of stochastic groundwater flow. This enables us to account for the effect of general trends while retaining the convenience of closed-form results. 2.3 Problem Formulation To illustrate the basic concept we consider steady flow in a heterogeneous multidi- mensional porous medium with a systematic trend in the log hydraulic conductivity. There are many ways to partition a particular log conductivity function into a ”trend” and a ”random fluctuation.” Here we require that the fluctuation have a spatial aver- age of zero and no obvious nonstationarities. The trend should vary more smoothly than the fluctuation. In practice, the trend can be estimated from a sample conductiv- ity function by applying a low-pass or moving window filter [67]. Given these general requirements, we assume that the log hydraulic conductivity is a locally isotropic, ' random field with a known spatially variable ensemble mean equal to the value of the trend at each point in space. We also assume that the fluctuation about the log con— ductivity mean is a wide-sense stationary random field with a known spectral density function. The random log conductivity fluctuation is approximately related to the piezometric head and velocity fluctuations by the following first-order, mean-removed 16 flow equations: 62h’ Bh’ 6f’ aha,“ + iii-(Kla—x; — Ji(x) 6.1:,- x E D, (2 1) , 8h’ 11 = Kg(x) [J,-( ) — Bari] x E D, (2 2) h'(X) = 0 X 6 PD Vh’(x) -n(x) = 0 x E [‘10, where J,(x) = —6}-z/6:1:,- is the mean head gradient and p,(x) = {if/8x,- is the mean log conductivity gradient. The notation Kg(x) is the geometric mean conductiv- ity. These equations are written in Cartesian coordinates, with the vector location symbolized by x and summation implied over repeated indices. The point values of the log hydraulic conductivity fluctuation f’ (x), head fluctuation h’ (x), and velocity fluctuation u§(x),i‘ = 1, 2, 3, are defined throughout the domain D. A homogeneous condition is defined on the specified head boundary I‘D and the specified flux bound- ary [‘10. Note that we consider f is solely the source of uncertainty applied in the aquifer system. 2.4 Spectral Solution Spectral methods offer a particularly convenient way to derive velocity statistics from linearized fluctuation equations. Taking advantage of scale disparity between the mean and fluctuation processes and invoking locally the spectral representation, one can solve (3.1) and (3.2) and obtain the following expression for predicting nonsta- tionary velocity variances in heterogeneous trending media: 0'3..(X) = C (Mt-(X))0§(X)Kg2 (301200, (2-3) 17 where C(“i‘x” Z I: I: (1 ‘ w. 5:11:- (xi) (1’ 02 3:17.th 5”‘“’d‘"‘“’2’ (2.4) Note sf; is the dimensionless spectral density function of f’ (x), 3,, = S f f /a}, 0,2. is the log conductivity variance, (.0,- is the wave number, 022 = w? +1.03, 0,2,1 and 0,2,2 are respectively the longitudinal and transverse velocity variances. The detailed deriva- tion process of (3.9) is very similar to that used in Gelhar (1993) for homogeneous media and is not repeated here. 2.5 Approximate Closed Form Solution Equation (3.9) applies to flows in general, mildly heterogeneous trending media. To obtain explicit results, one must in general evaluate the associated double integrals numerically. In most cases, these evaluations can be quite difficult. Most prior research focused on limited special cases for which (3.10) can be reduced to a form that allows exact, closed-form integration [48, 45, 17]. In this paper, we evaluate (3.10) approximately under general trending conditions. Our approximate solution is based on the following observations: 1. Trends in conductivity influence the variance dynamics through C (11,-(x)) and Kg(x)J(x) (see (3.9)). 2. It is the evaluation of C (p,(x)) for a general, multi-dimensional trend distribu— tion p,(x) that is difficult. The zuj(x)kj(x) term in (3.10) makes the integration analytically intractable. 18 3. It is, however, predominantly K 9(x)J (x) that controls the nonstationary spatial variance dynamics. 4. For most trending situations, change in the mean log conductivity over the characteristic length of small-scale heterogeneity (a correlation scale A) is small (or u§A2 << 1) since the mean conductivity is expected to be much smoother than the fluctuation [17]. To enable variance modeling under general, complex conditions, we propose approx- imating C (u,(x)) via the following Taylor’s expansion-based expression: 00(0) 0000‘» = 0(0) + a”, u,- + z C(O) (2.5) Essentially we suggest that the small, but hard-to—evaluate contribution to variance nonstationarity from C (it, (x)) be ignored relative to the more important contribution from K 9(x)J (x). This assumption may seem quite crude at first sight but proves to be highly effective and makes general, approximate variance modeling in nonstationary trending media possible. Previous studies investigating the efl'ects of trending are all based on full rigorous integration of (3.10) or full solution of (3.1) that is intractable unless the trends are assumed to be of special forms [44, 45, 48, 17]. These highly restrictive assumptions severely limit the practical utilities of the results. Substituting (3.11) into (3.9), we obtain: 0?...(x) = aitxiKitx>J2(x) f: f: (1— ff)? sntwldwtdwg (2.6) Equation (3.12) can be easily integrated in the polar coordinate system. The result is the following explicit expressions: 2 _ 2 2 2 0,,1 (x) — 0.3750f(x)Kg(x)J (x), (2.7) 0'2 (x) = 0.1250}(x)K§(x)J2(x). (2.8) “2 19 These expressions are independent of the specific form of the log conductivity spec- trum or covariance function for the isotrOpic case. Note that (3.13) and (3.14) are of the same form as the equations derived by previous researchers (e.g., [56, 17]) for statistically uniform flow, except that K g and J are now allowed to vary over space as a function of the mean log conductivity gradient. These general equations are simpler than the nonstationary expressions deveIOped by Li and McLaughlin [1995], Loaiciga [1993], and Gelhar [1993] for simple linear trending media. In the following section, we illustrate how these simple closed-form expressions can be used to quantify robustly groundwater velocity uncertainty with surprising accu- racy, even in the presence of complex and strongly nonstationary trending hydraulic conductivity. 2.6 Illustrative Examples Our examples consider two-dimensional steady-state flows in a confined aquifer. The hydraulic conductivity is assumed to exhibit both a randomly varying small-scale fluctuation and a systematic large-scale nonstationarity. The small-scale fluctuation is represented as a random field characterized by a simple exponential and isotropic covariance function with a log conductivity variance of 1.0 and a correlation scale A of 1.0. The large—scale nonstationarity is represented as a deterministic trend. To systematically test the effectiveness of the closed-form solutions, we consider a range of trending conductivity situations, from relatively simple and mild trends to trends that are much stronger, more complex, strongly multi-dimensional, and even discontinuous. For each problem, we first solve the mean deterministic groundwater 20 flow equation without accounting for the small-scale heterogeneity, use the resulting head to evaluate the mean hydraulic gradient, and then substitute it into the ex- plicit expressions in (3.13) and (3.14) to obtain the local variance values. We verify the accuracy of the closed-form variance solutions by comparing them with the cor- responding numerical solutions obtained from the first-order nonstationary spectral method [44, 42, 43] and Monte Carlo simulation (based on 10,000 realizations). In all test examples, the simple closed-form solutions proves to be robust, can capture complex spatial nonstationarities, and match well with the corresponding numerical perturbation solutions and Monte Carlo simulation. We present in this section results from two selected examples involving relatively difficult trending situations. We select these examples because we believe that if a methodology is able to handle such distinctly different mean log conductivity trends, it should be able to handle perhaps most trends that can be practically represented using field data in real-world groundwater modeling. Table 2.6 presents detailed information defining the large-scale mean trends, statistics describing the small-scale heterogeneity, and other inputs used in the examples. The first example involves complex multi-dimensional trends artificially generated by a random field generator with a large correlation scale. We have purposely made the trends more complicated than that may typically be delineated with scattered field data in order to test the robustness of the closed-form formulas. The second example involves discontinuous trends in the mean conductivity. Although we do not expect the closed-form solutions to apply right at the discontinuities where the trending slope u,(x) is infinite, we are interested in determining if and to what degree the closed-form solutions apply overall and away from the discontinuities. The ability for a stochastic method to model discontinuous trends is important since many practical applications require working with zones of distinct materials, sharp geological boundaries, and 21 Table 2.1. Parameter definitions for the examples Continuous trend Discontinuous trend an varTance 1 1 an correlation scale A 1 1 Geometric mean hydraulic A replicate from a random (see Figure 2.7) conductivity field with a correlation scale 40A and variance: 0.5 (see Figure 2.7) Domain length 80A 80A West boundary condition Const Head 100m Const Head 92m East boundary condition Const Head 100m Const Head 92m North boundary condition No flow No flow South boundary condition No flow No How multiple aquifers with different mean conductivities. 2.7 Results Figures 2.7 and 2.7 illustrate the mean conductivity distributions for the trending examples and the distributions of corresponding steady state, mean head distribu- tions. Although the driving head difference in both examples is uniform, the strong trends in the mean log conductivity yield nonuniform flow patterns which in turn cause strong nonstationarity in the velocity variances. Figures 2.7 and 2.7 show the complex, nonuniform distributions of the predicted velocity standard deviations for, respectively, the continuous and discontinuous trend- ing examples. The results are presented in a profile along the domain centerline and obtained using the approximate closed-form formulas, the numerically-oriented non- stationary spectral method, and the Monte Carlo simulation. These plots clearly show that, despite the simplifications and the strong multi-dimensional trending non- stationarities, the simple closed-form solutions reproduce well the corresponding first- order, nonstationary spectral solutions and allow capturing both the spatial structure and the magnitude of the highly nonstationary, complex uncertainty distributions. 22 80 60 ‘11,’ 'A . . I. M 5.0 x20» 8 Constant Head = 100m Head = 92m to «t- b N 20 40 x10» 60 80 Figure 2.1. Spatial distribution of the geometric mean conductivity and corresponding mean head distribution for the continuous trending case. The closed-form predictions of the velocity uncertainty also match reasonably well with those obtained from the Monte Carlo simulation for both examples. For the discontinuous trending case, the results becomes inaccurate at the discontinuities but the inaccuracies appear to be very localized. 23 k1 = 25, k2 = 5, and k3 =125 m/day 60 E s ‘- E II N s... t ‘7? >< I 8 .. a) g I ‘65 8 20 O o 10 20 30 40 50 60 7o 80 x,/A 0 Figure 2.2. Spatial distribution of the geometric mean conductivity and corresponding mean head distribution for the discontinuous trending case. 2.8 Conclusions In this paper, we have developed and demonstrated approximate closed-form formulas to predict velocity variances for flow in heterogeneous porous media with systematic trends in log conductivity. Examples involving nonstationary flows in complex trend- ing media are used to illustrate the approximate methodology. The results reveal that, despite the gross simplifications, the closed-form expressions are highly effective and can reproduce the corresponding numerical, nonstationary spectral and Monte Carlo solutions. The results also show how trends in hydraulic conductivity produce complex structural changes in the spatial distributions in the velocity variances. The closed form formulas make it possible to model the velocity uncertainty in nonsta- 24 Approximate Spectral Method 1.0 - ------------- Nonstationary spectral method 0 Monte Carlo simulation Nr=10000 0.0 7 A A 4 1 . l 1 1 . l A A x 6.0 . :— Approximate Spectral Method 5.0 E- ------------- Nonstationary spectral method 0 Monte Carlo simulation Nr=10000 Figure 2.3. Centerline profile of the predicted velocity standard deviations using the closed-form formulae, nonstationary numerical spectral method, and Monte Carlo simulation for continuous trending case. 25 : Closed-form 2-5 y - ------------- Nonstationary spectral method ' 0 Monte Carlo simulation Nr=10000 601 0| ; Closed-form 2-5 r - ------------- Nonstationary spectral method 0 Monte Carlo simulation Nr=10000 VTY P O "'I Figure 2.4. Centerline profile of the predicted velocity standard deviations using the closed-form formulae, nonstationary numerical spectral method, and Monte Carlo simulation for discontinuous trending case. 26 tionary trending media. The analysis represents a step closer to our ultimate goal to include a systematic uncertainty analysis as a part of routine groundwater modeling. We are currently in the process of extending the approximate methodology to general unsteady flow in both confined and unconfined aquifers in the presence of complex sources and sinks. 27 CHAPTER 3 Modeling Groundwater Velocity Uncertainty in Nonstationary Composite Porous Media 3.1 Abstract In this paper I present approximate, closed-form equations that allow modeling 2D nonstationary flows in statistically inhomogeneous aquifers, including compos- ite aquifers containing multiple zones characterized by diflerent statistical models. The composite representation has the effect of decreasing the variance of deviations from the mean, relaxing the limitation of the small-perturbation assumption. The simple formulas are illustrated with a number of examples and compared with a cor- responding first-order nonstationary numerical analysis and Monte Carlo simulation. The results show that, despite the gross simplifications, the closed-form equations are robust and able to capture complex variance dynamics, reproducing surprisingly well the first-order numerical solutions and the Monte Carlo simulation even in highly 28 nonstationary, variable situations. 3.2 Introduction Probabilistic theories of subsurface flow and transport have had a significant impact on the way we think about uncertainty and heterogeneity. However, they have not had much impact on the way that predictions are generated and reported in practical groundwater modeling studies [101]. One major reason for this significant gap lies in that most existing stochastic theories require that the aquifer of interest be statis- tically homogeneous, the hydraulic gradient statistically uniform, and the deviation from the uniform mean small [47, 42, 43]. A number of recent review articles stressed that if stochastic modeling is to become a viable practical tool, it must be made much more general and flexible [101, 47, 42, 43]. Specifically, a stochastic model must be able to easily incorporate site-specific aquifer structure before it can be routinely applied in practice. It must allow modeling flexible zonations, layering, and general trending as most real-world aquifers exhibit both ”structural” and ”random” variability and the statistics characterizing aquifer heterogeneity can vary (e.g., from region to region and layer to layer) in response to systematic changes in the distribution of aquifer materials. The research, for example by Loaiciga et. al., 1993[48], Li and McLaughlin, 1995[45], Rubin, 1995[68], Indelman and Rubin, 1995,1996[34, 36], Winter and Tartakavosky, 2002[88], and Lu and Zhang, 2002[51], represents first steps in this direction. These studies investigated flows in heterogeneous trending media and pre- sented explicit, special analytical solutions in illustration of the effect of nonstation- arity. derived closed-form solution for flow in special composite media that consist of multiple zones with nonstationary fluctuation statistics. These illustrative solutions 29 provided useful insight into how large-scale nonuniformity and zonations interact with small-scale heterogeneity, although they must be generalized before they can have a significant practical impact. There are a number of numerical approaches that can be used to analyze statis- tically nonuniform flow and transport in more general, statistically inhomogeneous aquifers. These include, for example, Monte Carlo methods [74, 75, 61, 47] and per- turbation techniques, such as the moment equation methods [20, 21, 99, 89, 90, 97] and the so called first-order second-moment methods based on Taylor’s expansions [80, 93, 6]. All of these methods, however, are computationally demanding when ap- plied to flow and transport problems of realistic size [17, 59], mostly because they all require solving large numbers of partial differential equations on very fine spatial dis- cretizations in order to predict the impact of the small-scale heterogeneous dynamics [47, 43]. In this paper we present a general, approximate methodology for modeling two dimensional groundwater systems in complex, composite media. In particular, we apply the popular spectral technique to derive a set of closed-form formulas for pre- dicting nonstationary velocity variances in aquifers that consist of complex zonations characterized by different statistical models. These formulas provide nonstationary predictive capabilities of a stochastic numerical model while retaining the convenience of analytical solutions. We demonstrate the effectiveness of the simple closed-form formulas using a number of examples involving complex material zonations and non- stationary conductivity statistics. 30 3.3 Problem Formulation We can make our presentation and discussion more specific by considering a rela- tively simple problem: steady-state simulation of hydraulic heads and groundwater velocities in a saturated region. Here we assume that the region is composed of multi- ple subregions containing systematically different materials characterized by different statistical models. We further assume that fluctuation in each subregion is locally isotropic and statistically uniform but their statistics can vary from region to region and the means trends vary more smoothly than the fluctuations. In practice, the trend can be estimated from a sample conductivity function by applying a low-pass or moving window filter [67]. The random log conductivity fluctuation is approximately related to the piezo— metric head and velocity fluctuations by the following first-order, mean-removed flow equations: 3271' 8h’ Bf’ 627,63,- + #21305; — J’(x)6_x,- x E D, C“) , , Bh’ u,- = Kg(x) J,(x)f —— x E D, (3.2) 83:,- h’(X) = 0 X 6 PD Vh’(x) - n(x) = 0 x 6 FN, where J,-(x) = ~8fi/8x, is the mean head gradient, p,(x) = Bf/Bzc, is the mean log conductivity gradient, and Kg(x) is the geometric mean conductivity. Both 11,-, K9 are in general spatially variable and may be discontinous across internal subregion boundaries. These equations are written in Cartesian coordinates, with the vector location symbolized by x and summation implied over repeated indices. The point values of the log hydraulic conductivity fluctuation f’ (x), head fluctuation h’ (x), and velocity fluctuation u;(x),z’ = 1, 2, 3, are defined throughout the domain D. A 31 homogeneous condition is defined on the specified head boundary I‘D and the specified flux boundary I‘N. Note that we consider f is solely the source of uncertainty applied in the aquifer system and is locally (within an individual subregion) stationary and globally (across different regions) nonstationary. 3.4 Spectral Solution in Composite Media Spectral methods offer a particularly convenient way to derive groundwater velocity statistics from linearized fluctuation equations such as (3.1) and (3.2). Invoking the spectral representation in each subregion for the locally stationary log conductivity perturbation gives +oo f’(x) = / exp(zk,~a:,-)de(k,x), (3.3) —00 where z = (—1)1/2, k,- is component 2' of the wave number k, and de(k,x) is the random Fourier increment of f’ (x), evaluated at k. The x dependence in the spectral amplitude reflects the fact that it is in generall globally nonstationary and can vary from region to region. The Fourier representation can be viewed as the continuous version of a Fourier series expansion of f (x) The random Fourier increment at a particular wave number is analogous to the random amplitude of one of the terms in the Fourier integral. Stationary Fourier increments within each region satisfy the following orthogonality property [66, 17]: de(k, x)dZ;(k’, x) = 3,,(k, x)6(k — k’)dkdk’ (3.4) where the asterisk superscript represents the complex conjugate, 6(.) is the Dirac delta function, and Sff(k,X) is the spectral density function of the log hydraulic conductivity [66, 17]. 32 [65] shows that the output (e.g., h’, 11.1) of linear transformations such as (3.1) and (3.2) are stationary only if the input (e.g., f’) is stationary and the transformations are spatially invariant. In the problem of interest here, spatial invariance implies that the fluctuation equations (3.1) and (3.2) should have constant coefficients with the boundaries sufficiently distant having no effect on velocity fluctuations in the region of interest. Such spatial invariance requirement is clearly not met because, for heterogeneous composite media, the coefficients n,(x) and J,(x) may both vary with x, [20, 21, ?, 70, 61, 33, 45, 88, 26, e.g.]. Like many investigators [2, 56, 17] and following the pioneering work of [17], here we seek approximate solutions to (3.1) and (3.2). Taking advantage of the scale disparity between the fluctuation process and the mean process, we further assume that the input log-conductivity, the dependent head, and velocity output fluctuation be stationary in a local sense (and away from the immediate proximity of boundaries), so that they also have an approximate spectral representation defined in terms of a location-dependent spectral amplitude +oo +oo h'(x) = [00 exp(zk,—x,~)th(k,x); u3(x) = [00 exp(zk,-:t:,-)dZuj(k,x). (3.5) This local homogeneity assumption implies that we regard [n(x) and J,-(x) in (3.1) and (3.2) as varying slowly in space relative to the characteristic scale of the fluctuation, that is, that they do not change significantly over distance corresponding to the correlation scale of h. Note that 11,71 will be a typical length scale for change in the mean gradient, so that the notion of local statistical homogeneity will be meaningful when the product of it, and the correlation scale is small relative to 1 [17]. By using the local spectral representations (3.5) in (3.1), the spectral amplitude for head in two-dimensional problems is —’Iki Ji (X) th(k, X) : k2 _ zkjuj (x) de(k,x); k2 = k? + kg, (3.6) 33 Substituting (3.5) and (3.6) in (3.2), we obtain the following spectral amplitude of velocity: dZu,(k,x) = Kg(x)J,-(x) (1 — k2 _l:;:::‘j(x)) de(k,x). (3.7) Now, if the :01 coordinate is selected to be aligned with the local mean hydraulic gra- dient (J1(x) = J (x), J2(x) = 0), then the spectral density functions in each geological face for velocity become Sum]. (k, x) = Kg2(x)J2(x) kikl kjkl (1 k2 — 2kmum(x)) (1 k2 + 2kmum(x) Sff(k’x)° (3'8) Integrating in the wave number domain, one obtains the following integral expression for evaluating the velocity variance: of... (x) = 0.. (ulaitle2tle2tx), (3.9) where C()—/+°°/+°° 1— ””1 1— ””1 (kx)dkk u.- ” — -—oo —oo k2~2kjlllj(X) k2+lkjpj(X) 8‘” , 1 2, (3.10) Note an is dimensionless spectral density function, Sff = S f f / Ufa?“ and 032 repre- sent respectively the longitudinal and transverse velocity variances. 3.5 Approximate Spectral Method To obtain explicit results, one must in general evaluate the associated double integrals numerically. In most cases, these evaluations can be quite difficult. Most prior research focuses on some very special cases for which (3.10) can be reduced to a form that allows exact, closed-form integration [8, 17, 95, 69]. In this paper, we evaluate 34 (3.10) approximately under more general conditions. Our approximate solution is based on the following observations: 0 Trends in conductivity influence the variance dynamics through Cu,(p) and Kg(x)J(x) (see (3.9)). o It is the evaluation of Cu, for a general, multi-dimensional trend distribution u,(x) that is difficult. More specifically, it is the presence of the zkj/ij(x) term in (3.10) that makes the integration analytically intractable. o It is, however, predominantly K 9(x)J (x) that controls the nonstationary spatial variance dynamics. 0 For most trending situations, change in the mean log conductivity over the characteristic length of small-scale heterogeneity (a correlation scale A) is small (or u§A2 << 1) since the mean conductivity is expected to be much smoother than the fluctuation [17]. To enable variance modeling under general, complex conditions, we propose approx- imating Cu,(;i) via the following Taylor’s expansion-based expression: 60”,. (O) 3M Cu,(ii) = Cu,(0) + Hj + z Cu,(0) (3.11) Essentially we suggest that the small, but hard-to-evaluate contribution to variance nonstationarity from Cu, be ignored relative to the much more important contribution from Kg(x)J(x). This assumption may seem quite crude but proves to be highly effective and makes general, approximate variance modeling in nonstationary media possible. Previous studies investigating the effects of trending are all based on full integration of (3.9) through (3.10) or full solution of (3.1) that is intractable unless the trends are assumed to be of special forms [44, 17, 34, 35]. These highly restrictive assumptions severely limit the practical utilities of the results. 35 Substituting (3.11) into (3.9), we obtain: 0,2“(x) = a}(x)K§(x)J2(x) [:0 [:0 (1 — k;:1)23ff(k,X)dk1dk2 (3.12) Equation (3.12) can be easily integrated in the polar coordinate system. For the statistically isotropic case, the result is the following simple, explicit expressions, independent of the specific form of the MK spectrum: 03,, (x) = 0.3750§(x)K§(x)J2(x), (3.13) 02 (x) = 0.1250;(X)K§(X)J2(x). (3.14) “2 Note these expressions are of the same form as those derived by many previous re- searchers [56, 17, e.g.] for statistically uniform flow, except that K g and J are now allowed to vary over space as a function of the trending and a} can vary from zone to zone. These general equations are simpler than the nonstationary expressions de- veloped by [45] and [17] for simple linear trending media. In the following section, we illustrate how these simple analytical expressions can be used to quantify robustly and quite accurately groundwater velocity uncertainty in complex, nonstationary hy- draulic conductivity fields. 3.6 Illustration Examples To systematically test the effectiveness, robustness, and generality of the closed-form solutions, we consider a range of nonstationary, composite media configurations and use different statistical models to represent the small scale heterogeneity. Our first example considers a relatively simple situation in which the overall domain of interest contains two embedded, rectangular areas of distinctly different materials. Our sec- ond example is patterned after a real situation involving more complex zonations of 36 irregular shapes delineated based on a map of surficial deposits. Our third example considers fan-shaped deposits of water-transported material (alluvium) that forms at the base of fractured mountain blocks where there is a marked break in slope. The alluvial fan is coarse-grained and very permeable at the mouth and becomes gradually finer-grained towards the edge as it meets with a large surface water body. In all three examples, the overall conductivity distributions are strongly nonsta- tionary both in the mean and fluctuation. The mean conductivity varies within the zones (example 3) and between them (examples 1 through 3) and these large-scale nonstationarities are represented as deterministic trends. The fluctuation statistics in each zone are characterized by an independent statistical model with a zone- dependent variance and correlation scale. It is our opinion that if a methodology is able to predict the groundwater velocity uncertainty in such diflerent composite media configurations, it should be able to handle perhaps most situations that can be realistically represented using field data in real-world groundwater modeling. Table 3.1 presents detailed information defining the aquifer parameters, boundary conditions, and other inputs used in the three examples. Figures 3.1 presents real- izations of the conductivity distributions along the domain center line for example 1 through example 3 respectively. Figures 3.2 through 3.4 illustrate the distributions of the mean conductivity, prescribed stresses, and corresponding steady state mean head. The strong trends and irregular zonations in the mean log hydraulic conduc- tivities yield nonuniform flow patterns which in turn cause strong nonstationarity in the velocity variances. To compute the velocity variances and demonstrate their accuracy, we follow the following three step procedure. First, we solve the mean deterministic groundwater flow equation without accounting for the small-scale heterogeneity. We then use the computed head to evaluate the mean hydraulic gradient and substitute it into (3.13) 37 Table 3.1. Parameter definitions for three test examples Example 1 Example 2 Example 3 lnITcorrelation structure Exponential / Hole Exponential / HoTe Exponential] Hole an variance 05,10, and 2.0 0.5,0.8,1.0, and 1.5 1.0 and 2.0 an correlation scale A (m) 0.5,1.0, and 2.0 0.5,0.8,1.0, and 1.5 5.0 and 10.0 Mean conductivity K 9 5,25, and 125 5,25,50, and 100 2.0 and kriged K field (m/day) (see Figure 3.2) (see Figure 3.3) (see Figure 3.4) Domain length (m) 80x80 160x160 2000x1500 Grid number[ASM{NSM) 81x81 161x161 251x186 Grid number MCS 256x256 512x512 1024x768 Global recharge (m / day) No recharge No recharge 0.002 West boundary condition Const Head 100m Const Head 100m No flow East boundary condition Const Head 99.2m Const Head 98.4m No flow North boundary condition No flow No flow No flow South boundary condition No flow No flow No flow ASM:approximate spectraF method; NSM:nonstaHonary spectral method; MCS:Monte Carlo simulation. and (3.14) to obtain the local variance values. Finally, we compare the closed-form solutions with the corresponding numerical solutions obtained from the first-order nonstationary spectral method [44, 45, 46, 42, 43] and Monte Carlo simulation (based on 10,000 realizations). 38 (a) 8.0 Mean - ------------- Mean + fluctuations 2.0 °-° 20 ‘ 4’0 6'5 5'0 (b) 8.0 2 z u 5‘ Mean A 1 2 a 5: - ------------- Mean +fluctuations 6-0 ,9. .i .1311 it A’ x i‘ ii f.)- ”‘3 I" 55% g 5 S 4.0 i9 ' #3 'EE 1.; :? ‘3 r. '_ #3 ‘ ,9 z: . 5% _ i . ’ ‘ iii * it 1 . s5 is! i i " ' i! 4 > g 2.01 i g’ 1 if ’f i)°'°u---2tr--ztr—-tir--sii * no i . o C 6.0» ‘ 5' Mean .A 1' . | 3, - ------------- Mean+fluctuations 4.0-i 3%,...- “1'15 5. g g ! __.. ~ 9 I ' l’SgIEEF, 31,57,735 55A i 2.0 :7 3521512: -g fai ”i=3: 'i ., L i égigi — if.“ ii 7541i"; 1! ‘5‘: 0.0; '1‘ ‘ ’ I -2.0 i 560 1600 1560 A 5000 X. (m) Figure 3.1. A realization of nonstationary conductivity distribution plotted along the domain centerline in the horizontal direction for a) example 1, b) example 2, c) example 3 3.7 Results and Discussion Figures 3.5 through 3.13 present the comparative results for all three examples. Fig- ure 3.5 to 3.7 present contour maps showing the complex, nonuniform distributions of the predicted velocity standard deviations. The results are obtained using the approx- imate spectral method (left column), the numerical nonstationary spectral method (middle column), and the Monte Carlo simulation (right column) based on the simple exponential covariance model. Figure 3.8 to 3.10 show the same results in a profile 39 K ac 6‘ Figure 3.2. Spatial distribution of the geometric mean conductivity and the corre- sponding mean head distribution for example 1 40 111 fly :[3 t1 :0 Ewe 5.12%“ 3:). fl U1 (.71 3 Figure 3.3. Spatial distribution of the geometric mean conductivity and the corre- sponding mean head distribution for example 2 41 Mountain block Recharge (0.002111 day) Lake/sea Stage H = 0.0 m [mm 0 500 1 000 1 500 2000 X1 (”1) Figure 3.4. Spatial distribution of the geometric mean conductivity and the corre- sponding mean head distribution for example 3 42 along the domain centerline. The results clearly show that, despite the simplifications and the strong multi-dimensional medium nonstationarities, the simple closed-form solutions reproduce well the corresponding first-order, nonstationary spectral solu- tions and allow capturing both the spatial structure (see Figure 3.5 to 3.7 ) and the magnitude (see Figure 3.8 to 3.10) of the highly nonstationary, complex uncertainty distributions. The closed-form predictions of the velocity uncertainty also match very well with those obtained from the Monte Carlo simulation for three examples. Figure 3.11 to 3.13 presents similar comparisons between the closed-form and numerical solutions based on an alternative an covariance model - the hole expo- nential covariance model. Although the hole-type model looks similar to the simple exponential model in the physical space domain, it is very different in the frequency domain at low wave numbers. Gelhar (1993) suggested that a hole-exponential model may represent the field data better than the simple exponential model. The results show that the closed-form solutions based on this hole model also match very well for all three examples with the corresponding numerical perturbation and Monte Carlo solutions. The surprisingly robust performance of the closed-form solutions for various highly nonstationary situations involving different media zonations, mean conductivity dis- tributions, flow and uncertainty patterns, boundary conditions, and conductivity covariance models suggests that the seemingly crude simplifying assumptions and empirical observations made in deriving the analytical formulas are highly effective. These simple assumptions can indeed capture the dominant factors controlling the spatial uncertainty distributions and make it possible to model velocity uncertainty in complex, nonstationary groundwater systems. The closed-form formulas do introduce errors near media interfaces as they make no explicit use of boundary conditions. The solutions become inaccurate near the zone 43 discontinuities, especially at the prescribed head boundaries (see Figures 3.5 through 3.13). However, these errors are limited to the immediate proximity (within 2 log- conductivity correlation scales) of the discontinuities. We feel these very localized errors are quite acceptable, especially considering what we have gained out of the simplifying assumptions - the enormous flexibility and efficiency in our ability to model complex, composite media systems. It should also be pointed out that our ability to model composite media has the effect of decreasing significantly the variance of deviations from the mean. This is why our first-order solutions match so well with the Monte Carlo simulations even under highly variable conditions. Large-scale changes in log conductivity that increase the variance around a constant mean are now treated as nonstationarities (e.g. trends) since the method does not require that the log conductivity mean be constant. There- fore, the small perturbation assumption becomes much less limiting an assumption as in other stationarity-based perturbation methods. The method’s requirement that all uncertainty must ultimately be related to sta- tionary random fluctuations in each medium component (zone) may still seem to be a significant limitation. In reality, this requirement reflects a very important tacit assumption of stochastic groundwater hydrology. To be specific, suppose we wish to know how predicted velocity uncertainty is affected by spatial variations in hydraulic conductivity. A number of probabilistic methods, including the one described here may be used to infer the (possibly nonstationary) ensemble statistics of velocity from the ensemble statistics of hydraulic conductivity. In order to apply any stochastic methods to practical problems, we need to obtain zone specific estimates of the con- ductivity statistics which form the basis for our stochastic analysis. In practice, these statistics are typically derived from a limited number of field measurements available in a zone. In most situations the only type of parameter nonstationarity that we can 44 hope to infer from field data is a large-scale trend in the local mean (both within the zones or between them). The closed-form formulas can accommodate such trends since they require that the mean removed parameter within a zone be stationary. This suggests that the method’s requirement that uncertainty be related to station- ary random parameters within a zone is really not a practical limitation. We need to do this any way if we want our composite media analysis to rely on data gathered in the field. 45 3.8 Conclusions In this paper we have developed and demonstrated closed-form formulas to predict velocity variances for two-dimensional flow in complex composite media. Three ex- amples were used to illustrate the approximate methodology. The results reveal that, despite the gross simplifications, the analytical expressions are highly effective and robust and reproduce surprisingly well the solutions of corresponding first-order and Monte Carlo predictions. The results also show how the nonuniform log conductivity structure changes significantly the spatial distributions of the velocity variances. In summary, the approximate spectral method makes it possible to model the ve- locity uncertainty in complex composite media. The analysis represents a step closer to our ultimate goal to include a systematic uncertainty analysis as a part of routine groundwater modeling. We are currently in the process of extending the approxi- mate methodology to general 3D, strongly unsteady flow in confined/ unconfined and anisotropic aquifers in the presence of complex sources and sinks. 46 i SE .. E's. _. teaming, .. s. ,1. M’ i 5.: s. t. a El: c... m... . E.~a . Ea: » t 3!. an: 1., ,...., a. m. .1 r. . a. . is . . rm 6 5.5 . . E! 5.5 El: 35 SE Eoi$ 3!. 1.0 .55 El.“ . SE Ea. . 5:22:85 250 2:22 35 E I: E In El: Eggs. .. .. _ a .. . .. . was 55, ER: 5.!“ SEER: 35 SW; .2! 5.3 . hawk: 7“.“ 9 .4? . . a s n. . E I: [was was El: .., El . s E E a ‘37., w. ,. SE u . a . _. is. . it. _. '0... .. I, 3.3 5.3. «is S 5 . . . . .. 8 00522 .98QO Bacozfiwcoz 00522 :26QO 9m§x9aa< Figure 3.5. Spatial distribution of the predicted velocity variances using the closed- form formulas, nonstationary spectral method, and Monte Carlo simulation based on the simple exponential covariance model (example 1) 47 mafia 3.5 5.5 SE; was El: 0 3!. A 25 «no \\ Ely "new: mom: 7“: 7min Si: E la 55 ER: 5.x; EH: El: —l._ _la ER: E I: _I._ ER: man: :5 m5 :5 ins: a: :3 :3 if; 2!. . 25 El: El: AWN... . . E‘s Foul: Eggs“ :3 i$5 .. :3 Es . a: awn; a; . . 5.5 £5 E 5 RE; E35 ER: EB: co=m_:E_m 0.80 2:22 8592 @6on bacozfimcoz 00:55. .96me 2m§x95a< form formulas, nonstationary spectral method, and Monte Carlo simulation based on 48 Figure 3.6. Spatial distribution of the predicted velocity variances using the closed- the simple exponential covariance model (example 2) ooom c or WEE.“ EN: Hitch) as: :5 IE.“ .5: EB.“ co_fi_:E_m 0.80 2:22 £0”. ESE Ex; ER; (5 £55 H.101.) WEB: 5.3 WEE.“ 9 :5 AW 2.5 Em: ...: 3.: 1.x: .8522 :26on Emcozflmcoz o ooom a. £55 . =5 flow-K: Dom X Z . .552 a?! as; .. ...: WEE: Rafi; _ g. “If!” .5: .55 3.x.“ 5g .65 :5 in. 0.. “5522 .9825 Sméxoaaa. form formulas, nonstationary spectral method, and Monte Carlo simulation based on Figure 3.7. Spatial distribution of the predicted velocity variances using the closed- the simple exponential covariance model (example 3) 49 0.4 ' Approximate Spectral Method - ........... . Nonstationary Spectral Method 0.3 Monte Carlo Method Nr=10l. US 021 ....... 0.1”" " 0.0 0.4 _ - — Approximate Spectral Method I ------------- Nonstationary Spectral Method 0'3 .' 0 Monte Carlo Method Nr=10000 I. ‘‘‘‘‘‘‘ W???) 80 X.(m) . 1 1 rrrrrr 50 C Figure 3.8. Centerline profile of the predicted velocity standard deviations using the closed-form formulas, numerical spectral method, and Monte Carlo simulation based on the simple exponential covariance model(example 1) 50 0.6 ’ — Approximate Spectral Method - ------------- Nonstationary Spectral Method Monte Carlo Method Nr=10000 0.5 0.4 p p 0.2 h P . ,,,,,,,,,,, P 0.1, y. OOFLMLJMAALLMALJELLILIAJI‘AAJLJAAJMLLJLALJAAL I 0.4 . — Approximate Spectral Method - ------------- Nonstationary Spectral Method 0 Monte Carlo Method Nr=10000 b§ 0.2 f Figure 3.9. Centerline profile of the predicted velocity standard deviations using the closed-form formulas, numerical spectral method, and Monte Carlo simulation based on the simple exponential covariance model(example 2) 51 . —— Approximate Spectral Method 0.10 :- - ------------- Nonstationary Spectral Method ~ 0 Monte Carlo Simulation Nr=10000 : — Approximate Spectral Method _ - ------------- Nonstationary Spectral Method . 0 Monte Carlo Simulation Nr=10000 .o —L C I Figure 3.10. Centerline profile of the predicted velocity standard deviations using the closed-form formulas, numerical spectral method, and Monte Carlo simulation based on the simple exponential covariance model(example 3) 52 0.4 ; —— Approximate Spectral Method ’ - ----------- - Nonstationary Spectral Method 0.3 L 0 Monte Carlo Metho . Nr=1 0 1-11 _ Approximate Spectral Method _— ----------- - Nonstationary Spectral Method 03 f o Monte Carlo Method Nr=10000 Figure 3.11. Centerline profile of the predicted velocity standard deviations using the closed-form formulas, numerical spectral method, and Monte Carlo simulation based on a hole-type covariance model (example 1) 53 0.6 " 0.5 — Approximate Spectral Method - ------------- Nonstationary Spectral Method 0 Monte Carlo Method Nr=10000 0.4 ’ 0.2 I 0.1 E 0.07m 0.4 Approximate Spectral Method - ............. Nonstationary Spectral Method 0.3 0 Monte Carlo Method Nr=10000 . .% A T'fi—‘fi-VFVVV' '0 20 4 60 80 100 10 140 160 X10") Figure 3.12. Centerline profile of the predicted velocity standard deviations using the closed-form formulas, numerical spectral method, and Monte Carlo simulation based on a hole-type covariance model (example 2) 54 —— Approximate Spectral Method - ------------- Nonstationary Spectral Method 0 Monte Carlo Simulation Nr=10000 .0 . O 00 wits I Approximate Spectral Method 0.10 . - ------------- Nonstationary Spectral Method . 0 Monte Carlo Simulation Nr=10000 0.08 - ‘3" 0.06 - b 100 x. (m) Figure 3.13. Centerline profile of the predicted velocity standard deviations using the closed-form formulas, numerical spectral method, and Monte Carlo simulation based on a hole-type covariance model (example 3) 55 CHAPTER 4 Quantifying Flow Uncertainty in Stochastic Groundwater Models 4.1 Abstract Solving the full version of stochastic perturbation equations is generally considered to be computationally expensive and analytically intractable. This paper presents approximate, closed-form formulas for predicting velocity variances in the presence of nonstationarity caused by hydraulic conductivity trends, nonlinearity of unconfined head distributions, nonstationary transient flows, and deterministic sources and sinks in modeling areas. These formulas rely on a linearization of stochastic nonstationary perturbation equations but do not require the common statistical stationary assump- tions. The formulas are illustrated with a number of two-dimensional flow exam- ples and compared with a corresponding first-order nonstationary numerical analysis (based on solving the full version of perturbation equations) and Monte Carlo simu- lation. The intensive numerical experiments indicate that the closed-form formulas can reproduce well the corresponding solutions of first-order numerical method under 56 presented nonstationary flow situations, except in the proximity of prescribed head boundaries and well locations. These localized regions are limited in 2 to 3 As from the boundaries and in 5 to 10 As from the well locations, where A represents the correlation length of the log hydraulic conductivity. 4.2 Introduction Over the years a number of stochastic techniques had been proposed for analyzing the role of spatial heterogeneity in groundwater flow systems [e.g.,[74, 75, 17, 8, 61, 24, 95, 69, 100, 42, 43] ]. The most important approximations introduced in stochastic groundwater analysis is the assumption of statistical stationarity. The technical definitions of stationarity are concerned with the behavior of the statistical properties of a spatially random function when its time and space origin is shifted [65]. When the concept is justified, stationarity assumptions are advantageous both methodologically and conceptually. More importantly, they enable us to use a variety of analytical techniques to solve stochastic flow problems (see, for example, [18, 8, 95, 69]). The stationarity-based stochastic theories have produced a number of insightful results on field-scale flow and transport and have unified the recent stochastic and classical deterministic modeling framework (see reviews by [16, 7]). However, groundwater flows at most sites are statistically nonuniform and, in some cases, strongly nonuniform in response to one or more of the following factors: 1) nonstationarity (or vertical and / or horizontal trends) in hydraulic conductivity, 2) internal and external sources/ sinks (e.g., pumping and injection wells), 3) distributed recharge, 4) geologic and hydrologic boundaries and the associated nonstationary stresses imposed, and 5) transient mean flow effects. If stochastic modeling is to become a routine tool for real-world groundwater investigations [9, 10], it must be able 57 to accommodate these nonstationarities and complex sources and sinks in addition to just ”random” small-scale heterogeneity since they are parts of almost every realistic application[47, 42, 43, 62,64, 59, 87, 19] . Many studies had derived improved closed-form formulas for predicting flow un- certainty in the presence of nonstationarity such as that caused by linear hydraulic conductivity trends[45, 34], nonlinearity of head in unconfined aquifer flows[15, 17], prescribed head boundaries [71], global recharge based on the known statistical struc- ture of recharge distributions[70, 23], and transient flows [83, 95, 69, 92, 52]. These derived formulas, however, are only applicable in the highly restrictive situations. Previous work in Ni and Li [2005a,2005b] focused on the quantification of the local- ized simple closed form formulas for predicting velocity uncertainty in the presence of complex trending and composite media. Their studies indicated that the formulas can reproduce well the corresponding first-order numerical method and Monte Carlo solutions except in the proximity of prescribed head boundaries. In this study we extend the study of Ni and Li[62, 64] to include more gen- eral nonstationary flow situations such as complex hydraulic conductivity trends, nonstationary transient flows, and deterministic sources and sinks in the uncon- fined aquifers. The approximate closed-form formulas are illustrated with a num- ber of two—dimensional flows and compared with a corresponding numerical spectral method[44, 47] and Monte Carlo simulation. 4.3 Problem Formulation Consider general unsteady flows in heterogeneous multidimensional porous media with multiple source/sinks and log hydraulic conductivity trends. The random log conductivity fluctuations can be approximately related to the piezometric head and 58 velocity fluctuations by the following first-order, mean-removed flow equations: 3271' _ of S 017 21,- ah' 11,7,- 5, , , ago—tr 1a., +1, 6t (.1, i. am, +( i. + km W ‘6’" ‘4'” 6h’ J- u'-(x) = —T (— — 7311' — J-f’) , (4.2) J 9 8x,- h J h’(X) : 0 X 6 PD: Vh'(x) ~ n(x) = 0 x 6 FN. for unconfined aquifers. These equations are written in Cartesian coordinates. In two-dimensional flow x is the position vector for the point ($1,152) in the domain of interest D. A homogeneous condition is defined on the specified head boundary I‘D and the specified flux boundary FN. Note that we consider log hydraulic conductivity f is solely the source of uncertainty applied in the aquifer systems. In equation (4.1), S = S(x) denotes the storage coefficient and T9 = Tg(x) represents the geometric mean of aquifer transmissivity. We assign notations Jj(X) = -sz/62:,- for mean head gradient in aquifers and uJ-(x) = 6f/ij for the spatial variation of mean log conductivity. 6, = 8271/ 611:? shown in equation (4.1) represents the second derivative of mean hydraulic head. The notation U = (S - 871/ at — q) /Tg represents the source and sink terms that combine the effects of transient mean head and deterministic sources and sinks. In two dimensional problems, q = q(x) represents the vertical inflow rate. Equations (4.1) and (4.2) represent the full version of perturbation equations for general two-dimensional flows. It is important to relate the physical properties to the terms shown in equation (4.1). In equation (4.1), the aquifer storage term S/Tg controls the transient processes, it,- term reflects the hydraulic conductivity trend that contributes to the small-scale variability. U term reflects the contribution of 59 sources / sinks and the transient mean flow to the small-scale variability. Three more terms 2Jj/I-l, 11ij / k, and flj/k in equation (4.1) are created to reflect the nonlinear flows in unconfined aquifers to the small-scale variability. The assumption that products of fluctuations can be neglected can only be justified when the fluctuation variances are small [8, 17, 95, 47, 69]. Here the perturbation equations describes the linear, nonstationary transformation from f’ to h’ to 11]. 4.3.1 Approximate Perturbation Equations Equations (4.1) and (4.2) apply to flows in general heterogeneous confined and un- confined aquifers. To obtain explicit results, one must in general evaluate these equa- tions numerically. Most prior research focused on limited special cases for which (4.1) through (4.2) can be reduced to simple equations that allow deriving closed form formulas analytically [48, 17, 70, 23, 45, 83, 95, 69, 92, 52]. In this paper, we evaluate (4.1) and (4.2) approximately under general flow con- ditions. Our approximate solutions are based on the following observations: 1. Although the full versions of perturbation equations contain many terms, their relative importance on the final predicted solutions of variances is quite different. 2. The evaluation of equations (4.1) and (4.2) for general, multi-dimensional flows is difficult. The consideration of the secondary terms in equations (4.1) and (4.2) makes the derivation analytically intractable. 3. All the effects of transient flows, hydraulic conductivity trend, source/sink terms, and unconfined flows influence the variance dynamics through these sec- ondary terms and the mean head gradient J (x, t). (see (4.2)). 60 4. The predominantly J (x, t) terms in equation (4.1) control the nonstationary spatial variance dynamics. Given these observations we propose to perform a major “surgery” on (4.1) and (4.2). In particular, we proposed to drop all secondary terms, retaining only terms that contribute significantly to the velocity variances. This converts the original perturbation equations (4.1) and (4.2) to the following simple equations: 62H 8f’ W — Jj‘a'x—j X E D, (4.3) Bh' u'-(x, t) = -—T (— —- J-f’) . 4.4) J 9 0233' J ( Note that the mean head gradient J,(x, t) in equations (4.3) and (4.4) are obtained from transient mean flow solutions. 4.3.2 Spectral Solutions (Approximate Spectral Method) Spectral methods offer a particularly convenient way to derive velocity statistics from linearized fluctuation equations. Taking advantage of scale disparity between the mean and fluctuation processes and invoking locally the spectral representation, one can solve the approximate perturbation equations (4.3) and (4.4) and obtain the following expression for predicting nonstationary velocity variances in heterogeneous unconfined aquifers [17, 62]: 0,2,..(x, t) = C - 0J2,(x)Tg2(x)J2(x, t), (4.5) 04.1.7: (1—“;°:1)’s..oadw. ....) Note Sff is the dimensionless spectral density function of f’ (x), Sff _-= S f f / a}, a} 61 is the log conductivity variance, w,- is the wave number, 002 = w? +1.12%, 0,2,1 and 032 are respectively the longitudinal and transverse velocity variances. The detailed deriva- tion process of (4.5) is very similar to that used in Gelhar (1993) for homogeneous media and is not repeated here. Equation (4.5) can be easily integrated in the polar coordinate system. The result is the following explicit expressions: 03(x) = 0.3750§(x)T§(x)J2(x), (4.7) 03(x) = 0.1250}(x)Tg2(x)J2(x), (4.8) where 0?, denote the longitudinal velocity variance and 03, represents the transverse velocity variance. These expressions are independent of the specific form of the log conductivity spectrum or covariance function for the isotropic case. In the following section, we use a set of examples to show how these approximate closed form expressions can be used to quantify robustly groundwater velocity un- certainty with surprising accuracy, even in the presence of strongly nonstationary in steady/ transient confined and unconfined aquifers. 4.4 Illustrative Examples and Numerical Consid- eration In our test examples we consider two-dimensional flows in bounded and rectangular areas. The general flow direction for each example is driven by two constant head boundary conditions at the east and west boundaries. The domain sizes for examples are either 80A by 80A or 160A by 160A depending on the complexity of flow patterns. We apply same statistical structure in all test examples to describe the small-scale 62 variability in the aquifers. Here the small-scale fluctuation is modeled stochastically by the exponential spectral density function with a Ink variance of 1.0 and a isotropic correlation scale A of 1.0 m. We divided the test examples into steady and unsteady two major groups (see Table 4.1). In steady state flows, we model a set of examples for both confined and unconfined aquifers. These conditions include uniform flows, trending conductivity fields, sources and sinks, and a more complex flow system involving a strong conductivity trend and sources/sinks in the modleing area. In unsteady flow examples, we apply a periodic head fluctuation at west boundaries to create the transient situation in both confined and unconfined aquifers. We only apply the Monte Carlo simulation in examples 4 and 6, which are considered to be the relatively complex situations in steady state flow and transient state flows. In order to provide appropriate grid space in Monte Carlo simulation to resolve small- scale variability, the grid numbers in Monte Carlo are adjusted to A/3. The number of realizations are 10000 in example 4 and 5000 in example 6. Table 4.1 details the input parameters used for these comparison examples. 4.5 Results and Discussion 4.5.1 Steady state situation Our aim in the first example is to quantify the efl'ect of an unconfined aquifer con- tributing to the estimation errors when we apply approximate spectral method to unconfined steady, and uniform two-dimensional flow problems. Figures 4.1 shows the comparison results of velocity standard deviations along the domain centerline for an unconfined aquifer. The results show that the approximate solutions for veloc- ity STDs inside the domain area have slightly overestimation in longitudinal velocity ST Ds and underestimation in transverse velocity STDs. Note that the nonlinear be- 63 Table 4.1. The definitions of physical and first-order numerical parameters for six examples Example description Size Grids (In) Steady State Example 1 lnk=2 (see Figure 4.1) 80 81 Example 2 Conductivity trend (see Figure 4.2) 80 81 Example 3 Source and sinks 80 81 Global recharge: (see Figure 4.3) R=0, 0.05 and 0.1 m/day Local recharge: (see Figure 4.4) R=0.1, Land 2 m/day Pumping wells (see Figure 4.5) Q=-200 and -1000 m3 / day Example 4 Complex flow (see Figure 4.6): 160 81 Global recharge R=0.001 m/ day Local recharge R=1 m/ day Lake stage=56m, leakance=10 day—l Pumping wells Q=-2000 m3 / day Transient State Example 5 lnk=2 (see Figure 4.7): 160 81 Specific yield:Sy=0.01 and 0.1 Example 6 gomplex flow (see Figure 4.9): 160 81 =0] 1: Global recharge R=0.001 m/ day Lake stage=56m, leakance=2 day"l Pumping wells Q=-2000 m3/day havior of unconfined aquifer(J is not a constant) makes the STDs gradually increase with the increase of 131 direction. When we apply the approximate spectral method to an unconfined aquifer, the head gradients J used for calculating approximate so- lutions are obtained from the head gradients of unconfined aquifers. As the results shown here, the contributions of the extra terms produced by unconfined situation in equation (4.1) are not significant. The unconfined efl'ect has been predominately taken into account by mean head gradients J. We found that there are inaccurate ar- eas near two constant head boundaries. The ranges of the boundary effect are about 3 As from boundaries. 64 In the second example we consider a complex trending situation in steady flow systems. Figure 4.2(a) shows the spatial distribution of large scale mean hydraulic conductivity and the corresponding mean head distributions for confined and uncon- fined aquifers. On top of the trend, the structure of small scale fluctuation is the same as the first example. Our objective in this example is to analyze the trending effects introduced in the confined and unconfined aquifers. Figure 4.2(b) and (c) show the profiles of velocity STDs modeled by approximate spectral method and nonstationary spectral method. In addition to the boundary effects, the results are acceptable. We found that the approximate solutions in the high conductivity area produced more error(in the central area of modeling domain, see Figure 4.2). This is because approx- imate solutions solely rely on the head gradient to predict the propagation of the flow uncertainty, in high conductivity area, the head gradient becomes very small and the trending effect becomes relatively strong. Our third example focuses on quantifying the effects of the sources and sinks in steady confined and unconfined aquifers. Although there are many types of sources/sinks existing in groundwater flow systems, we are especially interested in those that directly contribute to the source/sink terms in equation (4.1). Here we consider individual aquifer system including global recharge, local recharge, or pump- ing wells. Other mixing type sources/ sinks such as lakes or rivers will be discussed in other examples. The locations of local recharge area and the well location are defined in the central area (see Figure 4.4 and 4.5). The comparisons of approxi- mate spectral method and nonstationary spectral method for each model are based on applying different recharge rates or diflerent pumping rates in the modeling area. Figure 4.3 shows the results for different global recharge rates applied in the entire modeling area; while Figure 4.4 and 4.5 show the examples involving local recharge and pumping wells. For the recharge examples, the approximate spectral method and nonstationary spectral method show identical results even though the recharge rate 65 is increased to an extremely high value(36.5m/ year). Figure 4.5 shows the selected results that the recharge situations are replaced by pumping wells located at the do- main center point. Although there are significant differences between approximate and numerical solutions near well locations, we found that the diflerences are lim- ited in a range about 5A from the well location even though the pumping rate was increased to an extremely high value 2000m3 / day (in our example, 2100 m3/day pro- duces dry cells in unconfined aquifer). Away from the well location the approximate spectral method can accurately predict the velocity uncertainties. In example 4, we considered a complex conductivity trend, a local recharge with strong recharge rate at 1 m / day, a lake with extremely high leakance of 10 day—1, and a pumping well at strong pumping rate (2000m3 / day) in the modeling area. Figure 4.6(a) shows the conceptual model and the corresponding mean head distributions for confined and unconfined aquifers. Although the driving head diflerence in this example is fixed, the trend and strong source/sinks introduced in the modeling area yield nonuniform flow patterns. In addition to the approximate spectral method and the first-order numerical method, we apply the Monte Carlo simulation in this example for comparison purpose. Our Monte Carlo simulation includes two major processes: a random field generator based on Fast Fourier Transform algorithm [12, 5] and a efficient flow equation solver using Algebraic Multigrid method (AMG) [73, 77, 13, 58, 76, 31, 78] . Figure 4.6(b) and (c) show the profiles of predicted velocity STDs using approxi- mate spectral method, nonstationary spectral method, and Monte Carlo simulations. Except for the areas near boundaries and well locations, the results indicate that the solutions of approximate spectral method match well with the solutions of cor- responding first-order numerical method and Monte Carlo simulation. Although the approximate spectral method creates solutions based on ignoring secondary terms 66 produced by trend, nonlinearity of unconfined flows, and sources/sinks, the results indicate that these terms are relatively insignificant comparing to the mean head gra- dient J. Except for the special locations such as well and boundaries the mean head gradient are capable to accurately predict the velocity uncertainty. 4.5.2 Transient state situation In example 5 and example 6 we simply apply a periodic head boundary condition on west boundaries to produce a transient state flow pattern. The total modeling time for each example is 4 days. To satisfy the stability requirement for numerical nonstationary spectral method, we use a constant time step 0.1 day to model the transient processes in confined and unconfined aquifers. Figure 4.7 shows the conceptual model for example 5. In general, the storage coef- ficients for confined and unconfined aquifers have few orders of magnitudes difl'erences. In this example these storage coefficients are selected to cover general properties of unconfined aquifers (see Table 4.1 for detail). Figure 4.8 shows the transient processes of predicted velocity STDs at two monitoring points. The results show that the ap- proximate spectral method can predict well the corresponding first-order numerical solutions. Note that the nonstationary spectral method solves the full version of perturbation equation (4.1), not the approximate perturbation equation (4.3). This surprising result emerges the issue of computational efficiency in stochastic modeling. Based on our modeling experience in this example, the approximate spectral method needs few seconds to calculate the results by our Intel Xeon 3.6GHz workstation; however, using nonstationary spectral method takes about 1 hour to calculate the confined solutions and about 1.5 hour to obtain the final unconfined solutions. Figure 4.9 shows the conceptual model that are used for example 6. In this 67 example we include the Monte Carlo simulation to compare the first-order solutions (approximate spectral method and numerical nonstationary spectral method). Figure 4.10 shows transient processes of predicted velocity STDs at two monitoring points. The results clearly show that approximate spectral solutions match well with the cor- responding first order numerical solutions, however, two first-order methods all fail to capture the Monte Carlo solutions especially near high permeability area (moni- toring point (80,40)). Figure 4.11 shows the velocity STDs along the profile defined in Figure 4.9. Although the solutions obtained by approximate spectral method and nonstationary spectral method show identically, the inconsistent results also exhibit between first-order methods and the Monte Carlo simulation. If we compare the dif- ference between example 4 and current example, the difference is only the transient effect introduced in this example. First-order methods match well in example 4, how- ever, first-order results in this example are away from Monte Carlo solutions. For some locations, the Monte Carlo results show totally different patters of predicted velocity STDs. The great difference may come from the high-order terms that are truncated by first-order methods. This issue is beyond the scope of this study and was not discuss here. 4.6 Conclusion This study was designed to extend our ability to predict the uncertainty of ground- water flows in realistic complexity and sizes. We have performed a series of example studies that covered the most possible situations generally exhibited in groundwater flow problems. The approximate spectral method, which was derived locally based on the approximate perturbation equation, can reproduce well the solutions of first- order numerical method except for the locally strong stress areas such as hydrological boundaries and well locations. The inaccurate regions were limited to few correlation 68 80 I ' _ (a Notlow 60 - ._ 5 <2 in ‘1: “1] N '- . g a 8 8 a 5' 8 a E A .5 ° N ' N g, 40 .. II)....A_._ _._..— — .i .—._ _.—.]-.—A. .--g. >< - 8 8 g .1; . Pumping rae 20 _ Q= -1000 m Iday Figure 4.5. The conceptual model and solutions for example 3 (well): (a) the concep- tual model and the selected mean head distribution in an unconfined aquifer, (b) the longitudinal velocity (on) standard deviation along domain centerline ($2 2 40m), and (c) the transverse velocity (0,) standard deviation along domain centerline 73 ' stage = 56 m 1 leakance = 10 da ‘ =50m Head r”. d ;.u - , =- 000 m3/da lo'fl Figure 4.6. The conceptual model and solutions for example 4: (a) the conceptual model and the mean head distribution, (b) the longitudinal velocity (on) standard deviation along profile A—A’, and (c) the transverse velocity (0”) standard deviation along profile A-A’ 74 150 -_ No flow ; e 125 f 6? Z i; . . v Fl A100 :- .E ow direction E E * i” 8 v : ‘0. O 0 II N 75 _- '; 8 x _ a g 50 _:_ " Monitoring wells _ 8 _ Q) _ I 25 :- 0 h 1 1 1 1 I 1 1N01fl01w I 1 1 1 1 I O 50 100 150 x, (m) Figure 4.7. The conceptual model for example 5. The transient flow is driven by the west boundary where the head values are changed with time based on a sinusoidal function. 75 Monitoring well (40,80) Monitoring well (80,80) Sy=0.01 * —— Sy=0.01 - ''''''''' sy=0.1 p - ---------- sy=o_ 1 8i Mean Head (m) i? . 6 01 N A C" v C 01 o .b ! ! ! l ! I 2 (D a l ! ! ! ! I 2 (D Z 2 Time (day) 2 Time (day) Figure 4.8. The modeling results for example 5: (a) the mean head time series at two monitoring wells, (b) the time series of longitudinal velocity (on) standard deviation at two monitoring wells, and (c) the time series of transverse velocity (0”) standard deviation at well locations 76 ' , v, “ Global recharge Lake H: 0.001 m/da . Sta e=56m “ Lea ance =2 da " =50m Head -‘ 3 x. 0' 100 125 Figure 4.9. The conceptual model for example 6. In addition to the transient flow situation, a conductivity trend and multiple source/sinks are applied to the modeling area. 77 Monitoring well (40,80) Monitoring well (80,80) E rT—Vj‘f Mean Head (m) V . - .8- . . 57 55‘ M 54* . . - (b) 0.25, , 0.15 P°°°°°°°°Oooooooo°°°°°°°°°°°°°ooooooo°dl :3 . b : 0.10 : —— ASM : --------- NSM 0.05 - o MCS 0.00 . 2 ‘ * L a (c) 0.25, 0.20} ASM } —— ASM : --------- NSM ’ 0.15} o MCS > » b : 0.10: ........ 0.05 ‘ L 4 . 2 - Time (day) dr- L # l 1‘ ‘, 2 ‘ ‘ Time (day) Figure 4.10. The simulation results for example 6: (a) the mean head time series at two monitoring wells, (b) the time series of longitudinal velocity (on) standard deviation at two monitoring wells, and (c) the time series of transverse velocity (0,) standard deviation at well locations 78 t = 2.0 day Mean Head (m) Figure 4.11. The solutions at t = 2.0 (day) for example 6: (a) the mean head distribution along profile A-A’, (b) the longitudinal velocity (on) standard deviation along profile A-A’, and (c) the transverse velocity (0”) standard deviation along profile A-A’ 79 scales of log hydraulic conductivity. The results in this study also suggested that the secondary terms in confined and unconfined aquifers had relatively weak effects to the final predictions of flow uncertainty. In most flow situations, mean head gradients J are capable to capture the major small-scale dynamics. The capability of closed form formulas turns to computational advantages especially for the examples with realistic nonstationarity and sizes. It is appropriate here to mention the computation time for the presented numerical experiments. When we model the transient state and unconfined aquifer in example 5, the time scale for the approximate spectral method is in the order of minutes (totally 40 time steps), however, the time scale for the nonstationary spectral method was extent to the order of hours. The Monte Carlo simulation was the most expensive one that required few days to obtain 1000 replicates. By adjusting the balance between the accuracy and computational cost, the approximate spectral method represents a step closer to our ultimate goal to include a systematic uncertainty analysis as a part of routine tools in stochastic groundwater modeling. We found that two first-order methods, the approximate spectral method and nonstationary spectral method, all matched well with the solutions of Monte Carlo simulation in steady state flows but failed to capture the solutions of Monte Carlo simulation in transient state flows. The essential differences between first-order and Monte Carlo solutions may require further study to identify the key effects for these differences. This issue is out of the scope of this study and was not discussed here. Due to the complex theories and the highly computational cost, stochastic models are rarely considered in solving realistic flow problems. Most groundwater applica- tions are still using the deterministic model such as MODFLOW. The results in this study provide an opportunity to add one advance feature, the prediction of flow un- certainty, to current deterministic models. Although these simplified descriptions do 80 not capture the true nature of hydrogeologic variability, they may provide reasonable approximations for particular applications. Our current study focuses on developing an adaptive processes to interactively cor- rect the highly nonstationary areas, which are usually not captured by using the ap- proximate spectral method. Taking the advantage of the first-order numerical method that can capture the highly nonstationary dynamics and the advantage of closed form formulas that can efficiently predict moderate nonstationarity, the developed method should has the capability to handle the problems with realistic nonstationarity and sizes. 81 CHAPTER 5 Prediction of Velocity Uncertainty in Complex, N onstationary, and Heterogeneous Porous Media - A Hybrid Spectral Method 5.1 Abstract In our three previous papers, we developed and tested the approximate closed-form formulas to predict the velocity uncertainty under a number of nonstationary flow situations. The intensive numerical experiments revealed that the closed-form formu- las can reproduce well the first-order numerical solutions except in the proximity of prescribed head boundaries and well locations. Based on these findings, this paper presents an efficient hybrid nonstationary spectral method for predicting uncertainties (e.g., velocity variances caused by unmodeled small-scale conductivity variability) in complex groundwater flow systems. This method involves two major computational 82 steps after the deterministic mean flow equation is solved. The first step is to apply a set of closed form formulas (named approximate spectral method in this study) to predict the nonstationary variances for the entire modeling area. Then the first-order numerical spectral method is employed in the second step to correct the ”regional solution” in localized areas where the variance distribution is consider to be highly nonstationary. We illustrate the hybrid method with two synthetic examples. In each example the solutions of hybrid method were compared with the corresponding solutions of approximate spectra method, nonstationary spectral method and Monte Carlo simulation. The results reveal that the hybrid spectral method can efficiently handle large modeling areas and can accurately predict the detail dynamics of high nonstationarity 5.2 Introduction It is now generally agreed that natural subsurface environment is very heterogeneous. Soil heterogeneities, in particular, cause dramatic variation in hydraulic conductivity from point to point within a groundwater formation. Such variation sometimes ap- pear to be random, although many sites also exhibit trends which are related to the sedimentation processes that create stratified deposits of contrasting soils, alluvial fans, deltas, and glacial outwash plains. At first glance it may seem that small-scale variation in hydrogeologic properties should have relatively little effect on the larger- scale flows. In reality, small-scale fluctuations in hydraulic conductivity can have significant large-scale consequences, primarily because of the nonlinear relationship which exists between conductivity and groundwater flows [16, 7]. Over the years a number of stochastic techniques have been proposed for ana- lyzing the role of spatial heterogeneity in groundwater flow systems [e.g.[74, 75, 17, 83 8, 61, 24, 95, 69, 100, 42, 43] ]. Generally, they assume that heterogeneous physical parameters such as hydraulic conductivity are random spatial functions with known statistical properties. It follows that environmental variables which depend on these parameters are also random. For example, uncertainty in hydraulic conductivity in- duces uncertainty in hydraulic head and pore velocities. The flow equations provide a physical basis for relating the moments of random dependent variables (e.g. mean head and velocity, and their associated variances and covariances) to the hydrogeo— logical parameters which are the original source of uncertainty. In practice, however, it is very difficult to derive these moments without making approximations of one kind or another. The art of stochastic groundwater modeling lies in knowing which approximations are most appropriate for a given application. One of the most important approximations introduced in stochastic groundwater analysis is the assumption of statistical stationarity. When the concept is justified, stationarity assumptions are advantageous both methodologically and conceptually. They enable us to use a variety of analytical techniques to solve stochastic flow and transport problems [see, for example, [18, 8, 95, 69]]. However, groundwater flow at most sites are statistically nonuniform and, in some cases, strongly nonuniform in re- sponse to one or more of the following factors: 1) nonstationarity (or vertical and / or horizontal trends) in hydraulic conductivity, 2) internal and external sources/sinks (e.g. pumping and injection wells), 3) distributed recharge, 4) geologic and hydro— logic boundaries and the associated nonstationary stresses imposed, 5) systematically- varying aquifer thickness, and 6) transient mean flow effects. If stochastic modeling is to become a viable tool for real-world groundwater investigations, it must be able to accommodate these nonstationarities and complex sources and sinks in addition to just ”random” small-scale heterogeneity since they are parts of most realistic ap- plications. 84 Numerical analysis can be used to analyze general nonstationary flow problems under complex field conditions. However, those numerical methods are computation- ally demanding when applied to flow problems of realistic size [20, 21, 17, 54, 47, 42, 43, 62], mostly because they all require very fine spatial discretizations in order to resolve (either directly or indirectly) the underlying small-scale heterogeneous dy- namics. The numerical difficulty in calculating the moment, Green’s function, and sensitivity equations was pointed out in [47, 42, 43]. The implementation difficulty in the moment equation methods were stressed by [55, 55, 20, 21, 47, 42, 43]. In order to expect accurate solutions by any of these techniques, the discretization required would have to be significantly smaller than the typical scales of variation of the hydraulic parameter (often order(s) of magnitude smaller than that used in a corresponding deterministic model). As a result, computational effort involved often becomes pro- hibitively expensive, and thus applications of these methods to field situations have been very limited [17, 9, 47, 42, 43, 59]. Due to the restrictive applications of analytical methods and the high cost of nu- merical methods, there must be a way to take the major advantages of analytical and numerical methods. Previous studies in Ni and Li[62, 64, 63] had focused on iden- tifying the accuracy of approximate closed form formulas in complex nonstationary flows. Their systematic analysis suggests that the approximate closed form formulas are actually very efficient and can capture most nonstationary dynamics except for the proximity boundary and strong sources/sinks locations. Motivated by the critical assessments in Ni and Li[62, 64, 63], we present in this paper a hybrid spectral method for predicting velocity uncertainty in highly nonstationary groundwater flow systems. This method, based on solving stochastic perturbation equations, involves two major computational steps after solving the de- terministic mean flow equation. We first apply a set of closed form formulas to predict 85 the nonstationary variances for the entire modeling area. We then employ first-order numerical spectral method to correct the ”regional solution” in localized areas where the variance distribution is highly nonstationary (e.g., around inner boundaries and strong sources/ sinks). The boundary conditions for the local numerical solutions are based on the closed-form formulas. We illustrate the hybrid method with two syn- thetic two-dimensional flow problems and compare that with corresponding numerical spectral method[44, 47] and Monte Carlo simulation. 5.3 Problem Formulation We consider a unsteady flow in a heterogeneous multidimensional porous medium with multiple source/sinks and systematic trends in the log hydraulic conductivity. The random log conductivity fluctuation is approximately related to the piezometric head and velocity fluctuations by the following first-order, mean-removed flow equations: 5 3h’ 62H ah' 3f' , , ah' UJ-(X) = —Tg (a—x' — J bf) (5.2) h'(x) = O x 6 Pp, Vh’(x) - n(x) = 0 x 6 FN, for confined aquifers, and S 012' am 2.],- 6h 1MB,- 0 f’ Eat Hang-ax,” h 01:, ( h E)”- J’s—x,+Uf’ “1963) , Bh’ J , , “J(X) = ‘Tg (5;; " 73h Jaf) a (5 4) 86 h'(x) = 0 x 6 PD, Vh’(x) - n(x) = 0 x 6 FN. for unconfined aquifers. These equations are written in Cartesian coordinates. In two- dimensional flow x is the position vector for the point (2:1, 3:2) in the domain of interest D. A homogeneous condition is defined on the specified head boundary I‘D and the specified flux boundary [‘10. Note that we consider the log hydraulic conductivity f is solely the source of uncertainty applied in the 'aquifer systems. In equation (5.1) and (5.3), S = S(x) denotes the storage coefficient and T9 = Tg(x) represents the geometric mean of aquifer transmissivity. We assign notations Jj(X) = —072/ij for mean head gradient in either confined or unconfined aquifers and pJ-(x) = 0 f / 0:12,- for the spatial variation of mean log conductivity. 6,- : 8271/ 63:; shown in equation (5.3) is the second derivative of mean hydraulic head. The notation U = (S - air/at — q) /Tg includes the combination effects of time-varying mean head and deterministic sources and sinks. In two dimensional problems, q = q(x) represents the vertical inflow rate. Note that the head gradient .1,- in equation (5.1) to (5.4) are different in values because of the nonlinearity of head distribution in unconfined aquifer. The assump— tion that products of fluctuations can be neglected can only be justified when the fluctuation variances are small [8, 17, 95, 47, 69]. Here the perturbation equations describes the linear, nonstationary transformation from f’ to h’ to ui. 5.3.1 The Nonstationary Spectral Method (NSM) Spectral methods offer a particularly convenient way to derive head and velocity statistics from linearized fluctuation equations such as (5.1)-(5.4). The traditional stationarity-based spectral approach [2, 18, 17] is applicable only when the indepen- dent variable fluctuation ( f’) as well as dependent variable fluctuations (h’, u’) are 1 87 all wide-sense stationary. Papoulis [65] shows that the output (e.g., h’, u;) of linear transformations such as (5.1)-(5.4) are stationary only if the input (e.g. f’) is sta- tionary and the transformations are spatially invariant. In the problem of interest here, spatial invariance implies that the fluctuation equations (5.1)-(5.4) should have constant coefficients with the boundaries sufficiently distance to have no effect on head and velocity fluctuations in the region of interest. Such a spatial invariance requirement is clearly not met because of real—world complexities, the heterogeneous trending media, boundary effects, and sources / sinks introduced into aquifer systems. The NSM is a perturbation approach and does not require dependent fluctua- tions to be stationary. This method differs from other classical perturbation meth- ods primarily in the form of the spectral representation of the output variable fluc- tuation. More precisely, the dependent fluctuations are represented as a stochas- tic integral expanded in terms of a set of unknown complex-valued ”transfer func- tions” 212M (x, k, t). The fluctuations have the following Fourier-Stieltjes representation [66, 44, 17, 45, 47, 62, 64, 63]: +00 f’(X) = /_ exp(zkj$j)er(k), (5.5) +00 ”(39 t) = ¢h1(x,k,t) €XP(ij$j)de(k), (5.6) and +00 “3 (X, t) : ijf(xt k: t) (”(1)0195“)de (R), (5'7) where z = (—1)1/2, k, is component j of the wave number k, and de(k) is the random Fourier increment of f’ (x), evaluated at k. The Fourier representation can be viewed as the continuous version of a Fourier series expansion of f (x) The random Fourier increment at a particular wave number is analogous to the random amplitude of one of the terms in the Fourier integral. The symbols 112;, f, I/Juj f are unknown head and velocity transfer functions introduced to account for possible nonstationary flow 88 transformations. These transfer functions must be selected such that h’, u; satisfy the governing perturbation equations. Substituting (5.5)-(5.7) into (5.1)-(5.4) we obtain the following transfer function equations: S 3W] _ 321/01; _ . awhf T, at ‘ azjaz, + (22"? + “1) 6x,- +(zkjpj — k2)r,l)hf — sz-Jj + U, (5.8) 8 wujf : T9 (‘5' Ti? zkj’v/th) (5-9) for a confined aquifer, and S 6:10“ _ 821/)hf . 2Jj 6_z_/1_hf T9 at — 83,-8x]- + (22k) +HJ— — hJ—_—-) 8x,- +(ijflj— ijz—IL'L— k2 — JTM— fi__j —)'t[)hf— ijJj + U, (5.10) _ 3th qu-f _ TQJ(J axj + (£3; h _2k j)whf) (5.11) for an unconfined aquifer. Equations (5.8)-(5.11) are deterministic and complex- valued differential equations. Unlike the classical stationary spectral method, which requires transfer functions to be spatially invariant, the transfer functions introduced here are spatially variant. The transfer functions whf, waif obtained from (5.8)-(5.11) can then be used to derive the head and velocity variances in the same way as the classical stationary spectral method (e.g, [56, 17]): +oo 0,2, (X, t) : whfh‘i ka ”wa (x, 1", t)Sff(k)dka (5'12) +00 0121,- (X, t) : wuif(xa k: t)¢;,f(xa k, t)Sff (k)dka (513) where S f f (k) is the spectral density function of the log hydraulic conductivity [66, 17]. In this study we are especially interested in the velocity variances because the head variance is usually very small and, only some special spectral density functions that exhibit closed form formula for head variance. 89 5.3.2 Stochastic Modeling of Moderately Nonstation- ary Groundwater Systems - Approximate Spectral Method(ASM) Equations (5.1) to (5.4) represent the full version of perturbation equations. In order to introduce further approximations, one needs to clarify the hydrogeologic properties that are corresponding to the terms shown in equations (5.1) - (5.4) . In equation (5.1) and (5.3), p,- reflects the trend effect that contributes to the small-scale vari- ability. u,- / B term reflects the effect caused by varied aquifer thickness in confined aquifers. U term reflects the contribution of sources/ sinks and the transient mean flow to the small-scale variability. For the unconfined aquifer system, three more terms 2Jj / 5,11ij / h, and fii/I-z in equation (5.3) are added to reflect the contribution of unconfined aquifer flow to the small-scale variability. Mathematically, we can treat these terms as the extra stresses that are added on to a very simple flow situation. If we consider a quasi-steady flow situation and follow the observations in Ni and Li [2005c], we can aggressively ignore the secondary terms that are not significant in a moderately nonstationary flow situations. Giving these requirements, equations (5.1) - (5.4) can be reduced to the following simple formation: 82h’ Bf' W — Jib}; X E D, (5.14) Bh’ ....) = -.. (_ .. HI), J 9 8133‘ J ( ) Taking advantage of scale disparity between the mean and fluctuation processes and invoking the spectral representation, one can solve (5.14) and (5.15) to obtain the following expression for predicting velocity variances in a heterogeneous porous 90 media [17, 62, 64, 63]: aux) = oiJ2