A STATISTICAL CONTINUUM APPROACH FOR MASS TRANSPORT IN FRACTURED MEDIA By M A R K DONALD ROBERTSON B . A . S c . U n i v e r s i t y o f B r i t i s h C o l u m b i a 1986 A THESIS S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F THE REQUIREMENTS OFTHE DEGREE OFMASTER OF APPLIED SCIENCE in T H E F A C U L T Y O F G R A D U A T E STUDIES Department o f Geological Sciences We accept this thesis as c o n f o r m i n g to the required standard T H E U N I V E R S I T Y O F BRITISH C O L U M B I A September 1990 Q M a r k D o n a l d Robertson 1990 In presenting degree freely at this the thesis in partial fulfilment University of British Columbia, I agree that the available for copying of department publication this or of reference thesis by this his for and study. scholarly or thesis for her I further purposes financial of DE-6 (2/88) £-BP7BM&£<L OS , I^Ro requirements agree that It gain shall not The University of British Columbia Vancouver, Canada Date the may representatives. permission. Department of be is for an advanced Library shall make permission for granted by understood the that be allowed without it extensive head of my copying or my written ii ABSTRACT The stochastic-continuum model developed by Schwartz and Smith [1988] is a new approach to the traditional continuum methods for solute transport in fractured media. Instead of trying to determine dispersion coefficients and an effective porosity for the hydraulic system, statistics on particle motion (direction, velocity and fracture length) collected from a discretely modeled sub-domain network are used to recreate particle motion in a full-domain continuum model. The discrete sub-domain must be large enough that representative statistics can be collected, yet small enough to be modeled with available resources. Statistics are collected in the discrete sub-domain model as the solute, represented by discrete particles, is moved through the network of fractures. The domain of interest, which is typically too large to be modeled discretely is represented by a continuum distribution of the hydraulic head. A particle tracking method is used to move the solute through the continuum model, sampling from the distributions for direction, velocity and fracture length. This thesis documents extensions and further testing of the stochastic-continuum two-dimensional model and initial work on a three-dimensional stochastic-continuum model. Testing of the model was done by comparing the mass distribution from the stochastic-continuum model to the mass distribution from the same domain modeled discretely. Analysis of the velocity statistics collected in the two-dimensional model suggested changes in the form of the fitted velocity distribution from a gaussian distribution to a gamma distribution, and the addition of a velocity correlation function. By adding these changes to the statistics collected, an improvement in the match of the spatial mass distribution moments between the stochastic-continuum and discrete models was effected. This extended two-dimensional model is then tested under a wide range of network conditions. The differences in the first spatial moments of the discrete and stochastic-continuum models were less than 10%, while the differences in the second spatial moments ranged from 6% to 30%. Initial results from the three-dimensional stochastic-continuum model showed that similar statistics to those used in the twodimensional stochastic-continuum model can be used to recreate the nature of threedimensional discrete particle motion. iv TABLE OF CONTENTS Abstract « List of Tables • List of Figures v »' 1 X Acknowledgement X 1 V 1.0 Introduction 1 1.1 Fracture Characteristics 2 1.1.1 Aperture 1.1.2 Geometry 3 < 1.2 Transport Mechanisms 3 4 1.2.1 Advection 4 1.2.2 Dispersion 6 1.2.3 Diffusion 6 1.3 Continuum Models 7 1.3.1 Homogeneous Media 8 1.3.2 Non-homogeneous Materials 9 1.3.3 Fractured Media 10 1.4 Discrete Method for Network Simulations 11 1.5 The Stochastic-Continuum Method 12 1.5.1 Basic Method 12 1.5.2 Extensions to the Basic Method 13 2.0 The Basic Model 2.1 The Discrete Model 16 16 2.1.1 Aperture 16 2.1.2 Fracture Length 17 2.1.3 Fracture Density 17 2.1.4 Boundary Conditions and Gradient Orientations 18 2.1.5 Computational Method 18 2.1.6 Particle Transport in the Discrete Sub-Domain 19 2.1.7 Movement Statistics in Discrete Sub-Domain 20 2.1.8 Code Verification 21 2.1.9 Test for Representative Statistics 21 2.2 The Continuum Model 23 2.2.1 Solution of the Groundwater Flow Equation 23 2.2.2 Particle Transport in the Continuum Model 23 2.3 Comparsion of Full-Domain Discrete Model to Full-Domain Continuum Model 3.0 Extensions to 2D Model 24 48 3.1 Variable Aperture Model 48 3.2 An Improved Stochastic-Continuum Model 49 3.2.1 Velocity Distributions 49 3.2.2 Correlations in Particle Velocities 51 4.0 Testing the Extended 2D Model 84 4.1 Gradient Orientation 84 4.2 Variation in Aperture Distribution 88 4.3 Variation in Fracture Length 90 4.4 Conclusion 93 vi 5.0 Three-Dimensional Model 5.1 The Discrete Model 122 122 5.1.1 Aperture 123 5.1.2 Geometry 123 5.1.3 Boundary Conditions and Gradient Orientations 124 5.1.4 Computation Methods 125 5.1.5 Particle Motion in the Discrete Network 125 5.1.6 Code Verification 126 5.1.7 Movement Statistics 126 5.2 Continuum Model 127 5.3 Comparsion of the Discrete and Continuum Models 128 6.0 Conclusions and Discussion 135 6.1 Conclusions 135 6.2 Discussion 136 Bibliography 139 vii LIST OF TABLES Chapter 2 Table 2.1 Summary of input parameters for discrete sub-domain model used in Chapter 2 29 Table 2.2 Comparison of directional and velocity statistics for SO, 100 and 150 simulations of the discrete sub-domain model. 500 particles are used in each simulation. 30 Table 2.3 Comparison of directional and velocity statistics for 250, 500 and 1000 particles in the discrete sub-domain. The results are the compilation of 100 realizations 31 Table 2.4 Summary of input parameters for the full-domain models used in Chapter 2 32 Chapter 3 Table 3.1 Summary of input parameters for the discrete sub-domain model with variable fracture aperture 58 Table 3.2 Summary of input parameters for the discrete full-domain model with variable fracture aperture 59 Table 3.3 Summary of the first two correlation coefficients for the short and long correlation models from the uniform and variable fracture aperture discrete sub-domain simulations 60 Table 3.4 Coefficients for the first-order polynomial representation of spatial moment curves for the discrete and gamma(2) simulations 61 Table 3.5 Summary of the far-field percentage differences in the first and second spatial moments for the continuum models described in chapter 3 compared to the average discrete model 62 Chapter 4 Table 4.1 List of simulations carried during testing of the stochastic-continuum model with various input parameters 94 Table 4.2 Summary of the far-field percentage differences in the first and second spatial moments between the discrete and gamma(2) models for 4 orientations of the hydraulic gradient to the network 95 viii Table 4.3 Summary of the first two correlation coefficients for 4 Orientations of the hydraulic gradient to the variable fracture aperture, discrete subdomain simulations 96 Table 4.4 Summary of the far-field percentage differences in the first and second spatial moments between the discrete and gamma(2) models for three values of aperture standard deviation 97 Table 4.5 Summary of the first two correlation coefficients for three values of the standard deviation in fracture aperture 98 Table 4.6 Summary of the far-field percentage differences in the first and second spatial moments between the discrete and gamma(2) models for the three statistically different fracture networks 99 Chapter 5 Table 5.1 Summary of the input parameters for three-dimensional discrete model. .. 130 ix LIST OF FIGURES Chapter 1 Figure 1.1 Fracture trace map of a rock outcrop at Yucca Mountain, Nevada, U.S.A., from Barton and Hsieh [1989] 14 Chapter 2 Figure 2.1a,b One realization of the uncleaned and cleaned discrete sub-domain fracture network 33 Figure 2.2 The discrete network flow domain, showing boundary conditions, gradient orientations and the nomenclature for flow directions, after Schwartz and Smith [1988] 34 Figure 2.3 Fracture length histogram for network described in Table 2.1 (Network A) 35 Figure 2.4 Breakthrough comparisons for the discrete and continuum full-domain models. The discrete curve is the compilation of 10 simulations of the discrete full-domain model 36 Figure 2.5 The discrete full-domain network used for particle transport simulation in Chapter 2 37 Figure 2.6 Particle distribution plots for one realization of the discrete and continuum full-domain models after 70 days. Moment statistics of the plume are listed in the top right corner of the plot. Simulation parameters are listed in Table 2.4 38 Figure 2.7 Histograms of the first and second spatial moments of the discrete fulldomain plumes compiled from 100 Monte Carlo simulations at T = 70 days 39 Figure 2.8a-f Particle distribution plots showing the development of the plume for one realization of the discrete and continuum full-domain models (10-60 days) 40-45 Figure 2.9a First spatial moments ( X Q , ZQ) for the continuum and discrete fulldomain models. The discrete curve is the average of 100 realizations of the discrete network 46 Figure 2.9b Second spatial moments (SYY» ^ZZ^ f ° continuum and discrete full-domain models. The discrete curve is the average of 100 realizations of the discrete network 47 r t n e X Chapter 3 Figure 3.1 Particle distribution plots for variable, aperture discrete and continuum full-domain models, gaussian(O), after 35 days .-. Figure 3.2 Breakthrough curves for the plumes shown in Figure 3.1 65 66 Figure 3.3 Breakthrough curves for the discrete and continuum full-domain, variable aperture simulations. The discrete curve is the average of 10 67 simulations of the discrete full-domain model Figure 3.4a Development of the first moments (Xg, ZQ) of the plume for the variable aperture, discrete and continuum models. The discrete curves are the average of 100 realizations of the discrete network. ... 68 Figure 3.4b Development of the second moments (Svy, S77) f ° variable aperture, discrete and continuum models. The discrete curves are the average of 100 realizations of the discrete network 69 Figure 3.5a Histograms of velocities sampled in the uniform aperture, discrete subdomain model 70 Figure 3.5b Histograms of velocities sampled in the variable aperture, discrete subdomain model 71 Figure 3.6a Plots of the growth in the first horizontal and vertical spatial moments with time for the uniform aperture simulations. Comparisons are made between the curves from the three continuum models (gamma(O), histogram and gaussian(O)) and the discrete model averaged curves. Comparisons are also made between the curves of the two approximate velocity distribution continuum models (gamma(O) and gaussian(O)) and the histogram velocity distribution model curves 72 Figure 3.6b Plots of the growth in the second horizontal and vertical spatial moments with time for the uniform aperture simulations. Comparison are made between the curves from the three continuum models (gamma(O) histogram and gaussian(O)) and the discrete model averaged curves. Comparison are also made between the curves of the two approximate velocity distribution continuum models (gamma(O) and gaussian(O)) and the histogram velocity distribution model curves 73 Figure 3.7a Plots of the growth in the first horizontal and vertical spatial moments with time for the variable aperture (0.2) simulations. Comparisons are made between the curves from the three continuum models (gamma(O), histogram and gaussian(O)) and the discrete model averaged curves. Comparisons are also made between the curves of the two approximate velocity distribution continuum models (gamma(O) and gaussian(O)) and the histogram velocity distribution model curves 74 r i n e xi Figure 3.7b Plots of the growth in the second horizontal and vertical spatial moments with time for the variable aperture (0.2) simulations. Comparisons are made between the curves from the three continuum models (gamma(O), histogram and gaussian(O)) and the discrete model averaged curves. Comparisons are also made between the curves of the two approximate velocity distribution continuum models (gamma(O) and gaussian(O)) and the histogram velocity distribution model curves 75 Figure 3.8 Schematic diagram of paticle motion to illustrate the long and short correlation models. Nodes indicate locations of fracture intersections. 76 Figure 3.9a Plots of the growth in the first horizontial and vertical spatial moments with time for the variable aperture (0.2) simulations. Comparisons are made between the curves from the two correlated continuum models (gamma(long), and gamma(short)) the discrete model (discrete) and the uncorrelated continuum model (gamma(O)) 77 Figure 3.9b Plots of the growth in the second horizontal and vertical spatial moments with time for the variable aperture (0.2) simulations. Comparisons are made between the curves from the two correlated continuum models (gamma(long), and gamma(short)) the discrete model (discrete) and the uncorrelated continuum model (gamma(O)). 78 Figure 3.10a Plots of the growth in the first horizontial and vertical spatial moments with time for the variable aperture (0.2) simulations. Comparisons are made between the curves from the two autoregressive, long correlation, continuum models (gamma(2), and gamma(l)) and the discrete model (discrete). The gamma(l) curve is equivalent to the gamma(long) curve in Figure 3.9a,b. The uncorrelated continuum model (gamma(O)) is shown for reference 79 Figure 3.10b Plots of the growth in the second horizontial and vertical spatial moments with time for the variable aperture (0.2) simulations. Comparisons are made between the curves from the two autoregressive, long correlation, continuum models (gamma(2), and gamma(l)) and the discrete model (discrete). The gamma(l) curve is equivalent to the gamma(long) curve in Figure 3.9a,b. The uncorrelated continuum model (gamma(O)) is shown for reference. ... 80 Figure 3.11a Summary plots of the growth in the first horizontal and vertical spatial moments with time for the variable aperture (0.2) simulations. The progressive improvement of the model, from the gaussian(O) model (equivalent to the Schwartz and Smith model [1988]) through the histogram and gamma(O) models to the gamma(2) model, in the match to the discrete curves is illustrated 81 Figure 3.11b Summary plots of the growth in the second horizontal and vertical spatial moments with time for the variable aperture (0.2) simulations. The progressive improvement of the model, from the gaussian(O) model, (equivalent to the Schwartz and Smith model [1988]), through the histogram and gamma(O) models to the gamma(2) model, in the match to the discrete curves is illustrated 82 xii Figure 3.12 Particle distribution plots for the discrete and gamma(2) full-domain models at T=35 days 83 Figure 3.13 Histograms of the first and second spatial moments of the discrete, variable aperture (0.2), full-domain plume complied from 100 Monte Carlo simulations at T=35 days 84 Figure 3.14 Breakthrough curves for the discrete and continuum, full-domain, variable aperture simulations. Comparisons of the discrete gamma(2) and gaussian(O) models to the average of 10 discrete simulations 85 Chapter 4 Figure 4.1a Plots of the growth in the first horizontal and vertical spatial moments with time for a gradient orientation of 15 degrees 101 Figure 4.1b Plots of the growth in the second horizontal and vertical spatial moments with time for a gradient orientation of 15 degrees 102 Figure 4.2a Plots of the growth in the first horizontal and vertical spatial moments with time for a gradient orientation of 30 degrees 103 Figure 4.2b Plots of the growth in the second horizontal and vertical spatial moments with time for a gradient orientation of 30 degrees 104 Figure 4.3a Plots of the growth in the first horizontal and vertical spatial moments with time for a gradient orientation of 45 degrees 105 Figure 4.3b Plots of the growth in the second horizontal and vertical spatial moments with time for a gradient orientation of 45 degrees 106 Figure 4.4 Particle distribution plots for the discrete and continuum full-domain models with a gradient orientation of 15 degrees 107 Figure 4.5 Particle distribution plots for the discrete and continuum full-domain models with a gradient orientation of 30 degrees 108 Figure 4.6 Particle distribution plots for the discrete and continuum full-domain models with a gradient orientation of 45 degrees 109 Figure 4.7a Plots of the growth in the first horizontal and verticle spatial moments with time for the gamma(2) and discrete models uniform aperture simulations. The gamma(2) and gaussian(O) curves are shown for reference 110 Figure 4.7b Plots of the growth in the second horizontal and verticle spatial moments with time for the gamma(2) and discrete models uniform aperture simulations. The gamma(2) and gaussian(O) curves are shown for reference Ill Figure 4.8a Plots of the growth in the first horizontal and verticle spatial moments with time for the gamma(2) and discrete models variable aperture (0.1) simulations 112 Figure 4.8b Plots of the growth in the second horizontal and verticle spatial moments with time for the gamma(2) and discrete models variable aperture, (0.1), simulations 113 Figure 4.9 One realization of Network B 114 Figure 4.10 Fracture length histograms for Network B 115 Figure 4.11a Plots of the growth in the first horizontal and verticle spatial moments with time for the gamma(2) and discrete models Network B simulations 116 Figure 4.11b Plots of the growth in the second horizontal and verticle spatial moments with time for the gamma(2) and discrete models Network B simulations 117 Figure 4.12a Plots of the growth in the first horizontal and verticle spatial moments with time for the gamma(2) and discrete models Network C simulations 118 Figure 4.12b Plots of the growth in the second horizontal and verticle spatial moments with time for the gamma(2) and discrete models Network C simulations 119 Figure 4.13 Particle distribution plots for the discrete and continuum ful-domain simulations for Network B 120 Figure 4.14 Particle distribution plots for the discrete and continuum full-domain simulations for Network C 121 Chapter 5 Figure 5.1 XZ slice through 3-D network at Y=5.00 m 131 Figure 5.2 Schematic diagram of the three-dimensional fractures to illustrate the division of fractures into quadrants, the local coordinate systems within the quadrants and the boundary conditions of the quadrants. . 132 Figure 5.3 Schematic diagram of particle motion to illustrate the collection of particle motion statistics in 3-D 133 Figure 5.4 Particle plots for the 3-D sub-REV simulations 134 Of making many books there is no end, and much study wearies the body. Eccl 12:12 Praise be to God for getting me through this. I would like to thank my wife and family for their support. And the Natural Science and Engineering Reseach Council for funding this reseach and for providing me with a scholarship. 1 1.0 INTRODUCTION In the field of hydrogeology, we are faced with questions of how to clean up existing subsurface hazardous wastes and how to properly dispose of wastes we are creating today. It is becoming increasingly important to be able to accurately predict the movement of solutes in groundwater. Much work, with considerable success, has gone into the development of models which predict the transport of solutes through porous media. However, prediction of transport through geologic material where fractures are the major fluid conduits is not so well developed. Models of solute transport in fractured media are typically categorized as belonging to either the continuum method or the discrete method. Traditional continuum models, which work well for porous media, rely on assumptions which are not always valid for fractured geologic media. Traditional discrete network models avoid the pitfalls of these assumptions, but are limited by computer resources to modeling relatively small fracture systems. The stochastic-continuum method, introduced by Schwartz and Smith [1988], is a new approach to the continuum method. Instead of trying to determine dispersion coefficients for the medium and an effective porosity, statistics on particle motion (direction, velocity and fracture length) collected from a discretely modeled sub-domain network are used to recreate particle motion in a full-domain continuum model. The discrete sub-domain must be large enough that representative statistics can be collected, yet small enough to be modeled with available resources. Statistics are collected in the discrete sub-domain model as the solute, represented by discrete particles, is moved through the network of fractures. The domain of interest, which is typically too large to be modeled discretely is represented by a continuum distribution of the hydraulic head. A particle tracking method is used to move the solute through the continuum model sampling from the collected 2 distributions for direction, velocity and fracture length. This thesis will describe extensions to the original Schwartz and Smith method and will show that the method is computationally reasonable and accurate in predicting solute transport through fractured media. The remainder of this chapter introduces the concepts used in modeling solute transport in fractured media. An initial discussion of fracture characteristics is followed by descriptions of the traditional continuum and discrete methods. Finally the stochastic-continuum method, as developed by Schwartz and Smith [1988] is described. Chapter 2 describes the basic two-dimensional stochastic-continuum model. Discussion of preliminary results, similar to those of Schwartz and Smith [1988] leads into chapter 3, in which extensions to the basic 2-D model are developed. Chapter 4 analyses the sensitivity of the results, from the improved model, to various input parameters. Chapter S introduces a three-dimensional stochastic-continuum model and reports on some initial results from this model. Chapter 6 brings this thesis to a close with a summary of the major conclusions made in the previous five chapters, and a discussion of the application of the method to real-world situations. 1.1 Fracture Characteristics To understand solute transport in fractured media we first must review the characteristics of naturally-occurring fractures. Fractures are the major fluid conductors in the majority of igneous, metamorphic and highly cemented sedimentary bodies. They are found in nearly all consolidated and many unconsolidated geologic materials. Fractures are the expression of brittle failure due to differential stresses. Sources of such stresses include changes in pore pressure, weight of overburden, temperature, and tectonic loading (Narr and Currie 1982). Fractures are discontinuous, their bounding surfaces are typically rough and irregular, and often there is some degree of infilling between the surfaces. Because of the asymmetry of stress conditions causing failure and 3 the co-planarity of natural weaknesses, fractures often form with similar orientations. Fractures of similar orientation are termed a set. Typically a body of rock will have several fracture sets which interconnect, forming a network. The nature of flow through a fractured rock body is a function of the size and distribution of the apertures and the geometry of the network. 1.1.1 Aperture Aperture is the perpendicular distance between the two bounding surfaces of a fracture. Because fracture surfaces are typically irregular, fracture aperture will vary across the plane of the fracture. The irregularity of a fracture surface will be spatially correlated in that, if we think of the fracture surface as the surface of a planet, elevations at points near each other are likely to be related, ie. both above the average elevation or both below it. Brown and Scholz [1985] described methods of measuring the variability and correlation of fracture surface roughness. Brown et al. [1986] extended this work to fracture apertures between two correlated surfaces. Tsang and Tsang [1987] fitted a gamma distribution to measured apertures within individual fracture planes in granitic rock. Fracture aperture will also vary within a set. Snow [1967] found that the distribution of apertures measured at outcrops could be described by a lognormal distribution. Raven [1986] calculated the hydraulic aperture of over 1400 fractures, from borehole injection test and fracture logs in a monzonitic gneiss. The resulting histogram of apertures, though not exactly lognormal is highly skewed towards higher values. 1.1.2 Geometry The orientation, dimensions and density of the fractures will define the anisotropy and interconnectivity of the network. Due to the preferred orientation of the fracture sets, networks are typically anisotropic with respect to flow. Figure 1.1, from 4 Barton and Hsieh [1989], is a fracture trace map from a rock outcrop at Yucca Mountain, Nevada, U.S.A.. The exposed network of fractures is an example of the nature of anisotropy and interconnectivity commonly found in surface outcrops. In the 1 upper middle of Figure 1.1 the large number of short traces represent an intense swarm of small fractures. Short fractures are typically more abundant than longer ones as is borne out in work by Cruden [1977] and Rouleau and Gale [1987], who fit negative exponential and lognormal distributions to fracture trace length histograms of data measured in the field. The spacing between fractures, or fracture density, has been fit to various probability distributions including uniform, negative exponential, lognormal, and random (Priest and Hudson [1976]; Rouleau and Gale [1987]). Rouleau and Gale [1985] described methods of measuring fracture trace lengths, spacing, and density in the field. 1.2 Transport Mechanisms Solutes are mechanically transported, in the fractures, as they are carried along by the groundwater in response to the local hydraulic gradient. Mechanical transport is typically divided into advective and dispersive components. Advection can be thought of as transport at the mean flow rate and dispersion as a random component, due to the variation in local velocity. Solutes also undergo diffusion in response to concentration gradients. 1.2.1 Advection Integration of the one-dimensional form of the Navier-Stokes equation for incompressible fluids lead to the cubic law for uniform, laminar flow between parallel plates: It should be noted that with increasing depth of overburden the open length of the fractures will decrease, thus decreasing the interconnectivity of the network. 5 Q=(pg/12^)ib w 3 (1.1) where: Q is the volumetric flow rate [L /T] J •5 p is the fluid density [M/L ] J g is the gravitational constant [L/T ] u is the fluid's dynamic viscosity [M/LT] i is the hydraulic gradient [L/L] b is the aperture [L] w is the width of the fracture [L] Darcy, the father of quantitative hydrogeology, experimentally discovered a similar relationship for flow in porous media: Q=Kia (1.2) where: K is the hydraulic conductivity [L/T] a is the cross-sectional area of flow [L ] For porous media the hydraulic conductivity is related to the permeability (k) of the material by: K=kpg/ M (1.3) Comparing equations 2 and 3 with equation 1 we see that the intrinsic permeability of a parallel plate conduit is b^/12. As mentioned earlier, fracture surfaces are not smooth parallel plates but rather irregular surfaces. For the range of fracture closures expected during elastic deformation, Brown [1987] found that the actual flow rate through synthetically created rough walled fractures is 70-90% of that which would be predicted by the cubic-law. 6 1.2.2 Dispersion Dispersion of the solute in fractured rock occurs on three scales. 1) At the microscopic scale, velocity variations occur due to friction between the fluid and the bounding wall of the fracture. Solutes closer to the fracture walls will travel slower than those near the center of the channel. 2) At the scale of the individual fractures, variations in aperture within the plane of the fracture will create a myriad of flow paths with varying velocities and lengths. 3) At the network scale, mixing occurs at fracture intersections creating another level of possible flow paths for the solute, thus creating dispersion at a larger scale. 1.2.3 Diffusion At low groundwater velocities, diffusion due to concentration gradients becomes important. Fick's first law states that the mass of diffusing substance passing through a given cross-section per unit time is proportional to the concentration gradient F = -D (dC/dl) d (1.4) where: F is the mass of solute per unit area per unit time [M/LT] D is the diffusion constant [L /T] 2 d C is the solute concentration [M/L ] 3 dC/dl is the concentration gradient which is negative in the direction of flow. Diffusion in geological materials is reduced, in comparison to diffusion in free water, due to the longer flow paths around the solid particles. Effective diffusion constants for solutes in geologic materials range from 1-50% of their value in free water (Freeze and Cherry 1979). 7 1.3 Continuum Methods Because it is d i f f i c u l t to model a geologic m e d i u m at the pore scale, it is convenient to describe the m e d i u m as one or more homogeneous units, or c o n t i n u u m elements, w i t h volume-averaged properties. T o describe volume averaged properties f o r geologic units it is necessary to understand the concept o f a Representative E l e m e n t a l V o l u m e or R E V . A simple example o f the concept can be shown f r o m the determination o f porosity. A t any i n d i v i d u a l point i n a fractured rock the porosity is either 1 (solid) or 0 (void). A s the volume of investigation increases the porosity fluctuates but approaches a stable value somewhere between these two extremes. The m i n i m u m size v o l u m e to reach this stable value is termed the R E V for porosity. Permeability a n d d i s p e r s i v i t y as described below are both material constants w h i c h are based o n the assumption o f an R E V . In porous m e d i a , R E V ' s are quite small, usually o n the order o f centimeters. In fractured r o c k , the size o f the R E V ' s w i l l depend on the density of the fractures and the d i s t r i b u t i o n o f the fracture apertures. H o w e v e r , they are i n general m u c h larger than those f o r porous m e d i a . T h e v a l i d i t y o f using the c o n t i n u u m method f o r a fractured material w i l l be dependant on whether it can be approximated as an equivalent porous media. T h e c r u c i a l question is whether the R E V ' s f o r the material are smaller than the d o m a i n to be m o d e l e d . In sparsely fractured rock it is possible that the R E V ' s w i l l be larger than the d o m a i n to be modeled. L e a v i n g these questions of v a l i d i t y aside for a moment, I w i l l first r e v i e w the traditional c o n t i n u u m equations used to model transport through a porous m e d i u m . 8 1.3.1 Homogeneous Media In a homogeneous, isotropic porous medium it has been experimentally shown (Bear 1961) that mechanical dispersion can be represented as a Fickian process. In onedimension Fick's law for dispersion can be expressed as: F m =- D ( / d c m D = m d x ) av W (1.6) where: F m is the mass of solute per unit area per unit time due to mechanical dispersion [M/L T] 2 Dm is the dispersion coefficient [L /T] a is the dispersivity, an experimentally determined material constant[L] v is the fluid velocity [L/T] The most general form of the three-dimensions dispersion coefficient for an anisotropic medium is a fourth rank tensor with 36 non-zero components Bear and Verruijt [1987] By combining advection, dispersion and diffusion we can mathematically describe solute transport in a continuum. The conventional form of the one-dimensional advection-dispersion equation in a homogeneous, isotropic, saturated, porous medium under steady-state uniform flow conditions is: Dj(d C/dl ) - vj(dC/dl) = (dC/dt) 2 2 D-D m +D d (1.7) (1.8) where: 1 is the distance variable [L] The plume for an instantaneous point source will become elliptical in shape as it leaves the near source region and can be described by a gaussian distribution. For 9 simple geometries this equation can be used to analytically solve for solute concentration as a function of time and position, given initial and boundary conditions and continuum material constants (permeability, dispersivity, and diffusion constants). For more complex geometries finite difference or finite element methods are used. Large scale tracer tests for relatively homogeneous porous media have indicated that dispersivity is variable at small scales but approaches a constant value at a large scale. Thus dispersivity in such cases can be modeled as a simple function of distance, though the form of the equation is rarely known a priori. A further complication is the difficulty in obtaining values for dispersivity, as field-measured values tend to be orders of magnitude larger than laboratory-measured values for the same material (Freeze and Cherry, 1979). This is due to the fact that lab samples cannot reproduce the large-scale inhomogeneities found in situ. 1.3.2 Non-homogeneous Materials Unfortunately most geologic media, including fractured rock, are nonhomogeneous at some scale and the Fickian model of dispersion breaks down. To model transport in non-homogeneous material hydraulic conductivity is described as a spatiallycorrelated stochastic field. Smith and Schwartz [1980] used a particle tracking technique to model transport through a stochastic velocity field created by applying boundary conditions to the stochastic hydraulic conductivity field. Macro-dispersion (dispersion due to REV size or greater inhomogeneities in the material) was accounted for by the variation in the velocity field. Several authors (Gelhar et al. 1979; Gelhar and Axness 1983; Dagan 1982,1984; Molz et al. 1983; Guven et al. 1984; Guven and Molz, 1986; Neuman et al. 1987) have used the distribution of hydraulic conductivity to determine macro-dispersivity tensors, which can be applied to the conventional form of the advection-dispersion equation. A fundamental assumption of these methods is that the random hydraulic conductivity field is stationary. Stationarity implies that there is no 10 spatial trend in the distribution of the parameter. Thus the hydraulic conductivity distribution obtained from sampling spatially separate points can be used to statistically define the unknown value of hydraulic conductivity at any given point. The effect of this assumption is that as the solute passes through the medium it will experience a representative distribution of hydraulic conductivities. The resulting plume will approach gaussian behaviour as time and distance increase. Wheatcraft and Tyler [1988] proposed a method for transport in heterogeneous materials, based on fractal flow paths, that never reaches Gaussian behaviour. If transport in heterogeneous material never reaches Gaussian behaviour, then dispersivity can never be considered a constant and a REV for dispersivity cannot be defined. 1.3.3 Fractured Media As mentioned earlier, in order to use the continuum method for fractured rock we must assume that the network can be approximated as a porous medium, in that REV parameters can be determined. Long et al. [1982], Long and Witherspoon [1985] and Shimo and Long [1987] analyzed the effect of fracture geometry on the REV for permeability, to determine when flow through fractured rock can be modeled by continuum methods. They found that increases in fracture interconnectivity decreased the size of the REV, increasing the likelihood of equivalent porous media flow, while a broader distribution of aperture and smaller distributions of fracture orientation increased the size of the REV, decreasing the likelihood of equivalent porous media flow. Similar work for mechanical transport was carried out by Endo et al. [1984] for simple fracture networks. They found that a fracture system that behaves as a continuum for fluid flux might not behave as a continuum for mechanical transport. 11 1.4 Discrete Method for Network Simulations T h e discrete method f o r f r a c t u r e d m e d i a avoids the need to determine equivalent porous m e d i a parameters b y m o d e l i n g solute transport through a n u m e r i c a l representation o f a fracture network. A fracture network i n t w o - d i m e n s i o n s is represented as a series o f inter-connecting line segments w i t h assigned apertures. In three-dimensions the line segments are replaced b y t w o - d i m e n s i o n a l f i n i t e difference or finite element g r i d d e d planes. B y a p p l y i n g appropriate boundary conditions the hydraulic head a n d velocity can be solved f o r , at each nodal location. Because the exact location a n d aperture o f each fracture i n a body of rock is impossible to determine, statistical methods are used to create realizations o f the network. S a m p l i n g f r o m k n o w n distributions f o r fracture aperture, length, orientation and density we c a n create possible realizations o f what the actual network may look l i k e . B y a p p l y i n g a r a n d o m w a l k technique we c a n determine what the d i s t r i b u t i o n o f a solute plume w i l l look l i k e after travelling through the realization. F r o m the first and second spatial moments o f the plume d i s t r i b u t i o n , (mean a n d variance), it is possible to back calculate the e f f e c t i v e m a c r o - d i s p e r s i v i t y f o r an equivalent porous media. U s i n g the Monte C a r l o technique o f sampling f r o m numerous realizations, we can calculate the average value a n d variance f o r the moments o f d i s t r i b u t i o n . U n f o r t u n a t e l y the discrete method is l i m i t e d b y computational constraints o f existing c o m p u t e r resources to small numbers o f fractures. We are not only limited b y R A M (random access memory) w i t h respect to the size o f the d o m a i n we can m o d e l , but also by execution speed as to the number o f stochastic realizations that c a n be reasonably solved. 12 1.5 The Stochastic-Continuum Method 1.5.1 Basic Method To overcome the computational limitations, yet retain the parametric simplicity, of the discrete method, Schwartz and Smith (1988) introduced the stochastic continuum method. The idea behind the stochastic continuum method is to collect a set of statistics on particle motion, (velocities, directions of motion away from fracture intersections, and distances between fracture intersections), by tracking particles in a discrete sub-domain of the actual problem area and then, applying rules based on those statistics, move particles through the larger domain, which can be modeled as a continuum. The key assumptions of the model are: 1) As in the traditional continuum method the parameters we use, (particle motion statistics), must be representative of the network. In the same way that we defined REV's for dispersivity etc., we must define REV's for the motion statistics. As we increase the size of the domain of investigation the motion statistics should approach a constant value. For example, the mean velocity along horizontal fractures will approach a constant as the size of the sub-domain in which the statistics are collected increases. Thus our first assumption is that the discrete sub-domain model is larger than the REV for each of the statistics collected. It is also important that this sub-domain can be modeled discretely by available resources. 2) The statistics that we collect must be sufficient to recreate the nature of discrete particle motion in the continuum model. The behaviour of particle motion in the 13 continuum model should in a statistical sense mimic the motion of the particles in the discrete network. As with porosity, as we increase the size of the discrete sub-domain model the movement statistics should approach a constant value. Thus by modeling various sizes of the discrete sub-domain we can determine the size of the REV's and thus determine the validity of assumption 1. To test assumption 2, and thus the validity of the method, we compare break-through curves and moment statistics for equivalent-sized continuum and discrete numerical models which are greater than the size of the REV. The fracture statistics used to create the full-domain discrete models, against which the stochasticcontinuum model is compared, are the same as those used to create the sub-domain discrete models. 1.S.2 Extensions to the basic method The goal of this thesis is to expand on the Stochastic-Continuum method. Areas of investigation include: 1) Comparison of the stochastic-continuum model to discrete models where the size of the domain is greater than the stochastic-continuum REV (assumption 2 test). Due to computer limitations, Schwartz and Smith [1988] were unable to carry out such tests. 2) Sensitivity analysis of the solute distribution to: a) variable orientation of the hydraulic gradient b) variable fracture aperture and c) variable fracture length. 14 Work by Schwartz and Smith [1988] was limited to uniform apertures and one network realization. 3) The development of a three dimensional stochastic-continuum model. In three dimensions the problem becomes more complicated, and more realistic, as the particles are able to move across the plane of the fracture. Decisions must be made on what statistics will embody the nature of particle motion in 3-D. Ultimately it would be desirable to compare the model to an existing field site, as testing the model against discretely modeled simulations is somewhat artificial. Unfortunately tracer tests, in fractured material, have been limited to single fractures, or a small set of fractures on a sub-REV scale, (Abelin, 1986; Cacas et. al. 1990). However, comparison to full-scale discrete network simulations, which act as a stand-in for an actual network, indicate that the model will prove useful for modeling of large scale fracture flow systems. F i g u r e 1.1 Fracture trace map o f a rock outcrop at Y u c c a M o u n t a i n , N e v a d a , U . S . A . , f r o m Barton and H s i e h [1989] 16 2.0 T H E BASIC M O D E L As described in the introduction, the stochastic-continuum method is made up of two parts: 1) the discrete sub-domain model from which particle movement statistics are collected and 2) the full-domain continuum model in which the particles are moved by rules based on the statistics collected in the discrete model. This chapter reviews the design of the basic stochastic-continuum model described by Schwartz and Smith [1988], and extends the testing of the model carried out by Schwartz and Smith [1988] to larger domains. 2.1 The Discrete model The two-dimensional discrete model used in this thesis is very similar to the one developed by Schwartz and Smith [1988]. It models two sets of orthogonal fractures in a rectangular domain, and moves particles based on the random-walk method. Analysis of networks with variable fracture orientations is postponed until Chapter 5 as the threedimensional discrete model described there allows for more freedom in the orientation of individual fractures, as well as in the mean orientation of the sets. The remainder of this section will describe the basic components of the discrete model used in this thesis. 2.1.1 Aperture In the two-dimensional model, fracture apertures are assumed to be constant over the entire length of the fracture, variability in fracture aperture across the fracture plane will be considered in the 3-D model. Previous work on the stochastic-continuum method, (Schwartz and Smith, 1988), considered networks in which fracture aperture was held constant across the network. One of the initial goals of this work is to extend the analysis of the model to allow for variable apertures for the individual fractures. As will be seen in Chapters 3 and 4, the effect of variable aperture on the resulting 17 contaminant plume is quite pronounced. Apertures are sampled from a lognormal distribution. For the uniform aperture models a mean value of 32 u was chosen as a representative size of fractures commonly found in granites and similar intrusive bodies at depth (Gale, 1982; and Neretnieks, 1982). 2.1.2 Fracture Length Fracture lengths are sampled from the tail of a negative exponential distribution, given a minimum length. The minimum length cutoff is a simplifying assumption of the model to reduce the computational burden of solving for the flow system. The smallest fractures do not contribute a great deal to the overall flow regime as they typically are poorly connected to the flow system. The smallest fractures do, however, contribute to the retardation of the contaminant by diffusion into them, as described by Neretnieks [1982]. The addition of appropriate retardation factors due to micro-fracture diffusion, in addition to absorption and adsorption and chemical effects, will be left for future work. 2.1.3 Fracture Density Fractures are positioned randomly in an area larger than the modeling domain to allow for fractures which are centered outside of the domain to extend into the domain (Figure 2.1a). This procedure is followed to avoid producing lower than average fracture densities near the domain boundaries. Rouleau and Gale [1985] measured a total fracture density of 6.4 m~^, for the four fracture sets observed in their research in fractured granites at the Stripa Mine, (Sweden). However, many of these fractures were found to be dry, (Dverstorp and Andersson, [1989]), suggesting a lower fracture density for fluid flow and transport modeling. The number of fractures generated in the stochastic-continuum discrete model is set to create a total fracture density of 1.4 m~', which gives a fairly well connected network. The fracture density measured in the IS model was calculated using the method of Rouleau and Gale [1985]. Once all the fractures have been generated, dead end and isolated fracture segments are removed, or cleaned off, and node numbers are assigned to each intersection (Figure 2.1b) . 2.1.4 Boundary Conditions and Gradient Orientations To simplify the collection of statistics on particle movement, the discrete subdomain model is solved assuming a uniform, steady-state flow field. The hydraulic head at the four corners of the domain are calculated so as to produce an uniform flow field at the desired angle to the fracture network. From the four corner values, constant heads at boundary nodes are interpolated and entered into the model, (Figure 2.2). The statistics characterizing particle motion through the network will vary with the orientation of the hydraulic gradient. Thus to apply the statistics collected in the discrete model to a continuum model with a non-uniform flow field, statistics must be collected for various orientations of the mean hydraulic gradient to the fracture network. Schwartz and Smith [1988] chose fifteen degrees as the interval over which to interpolate the movement statistics. Because of the symmetric nature of the network, seven orientations of the gradient in the discrete model, between 0 and 90 degrees, are deemed sufficient to delineate all possible gradient orientations, (Figure 2.2). The magnitude of the gradient is kept constant for all seven orientations by changing the constant head boundary values. In the bottom left hand corner of Figure 2.2 the nomenclature for flow direction, (vertically upwards being direction 1 etc.), is shown. 2.1.5 Computational Method Nodal relationships (distance and aperture along the segments between fracture intersections) are stored in a finite difference formulation. The resulting system of Due to the fine mesh spacing of the discrete model (0.025 m), it is possible that fractures which appear to be connected in Figure 2.1a may not be truly connected, and thus are removed in the cleaned network, Figure 2.1b. 19 equations, (nodal relationships and boundary conditions), is non-symmetric and typically has a very large band-width. To reduce the amount of storage required for such large band-widths, a sparse matrix direct solver is used, which requires only the non-zero elements to be stored (Duff et al., 1987). Hydraulic heads are determined at the nodal points, from which the gradients between fracture intersections can be calculated. The calculated gradients, along with the fracture apertures and values of fluid density and viscosity are substituted into the cubic-law equation (Equation 1.1) to calculate velocities in the fractures. 2.1.6 Particle Transport in the Discrete Sub-domain Once the velocities are known a random-walk method is used to transport solute through the network. The random-walk method assigns the mass of solute to numerous particles, each particle representing a small percentage of the total mass. The particles move through the fractures, by advection, at the rate of the local flow velocity in the fractures. A t fracture intersections a frequency histogram, based on the percentage flow in each direction, is sampled to determine in which direction the particle will move. This probabilistic splitting of the contaminant plume is equivalent to complete mixing in the fracture intersections. A n instantaneous point source for introducing the particles into the model is chosen, as sensitivity analysis carried out by Schwartz and Smith [1988] showed no appreciable difference in the statistical results between point sources and line sources. The source location is chosen at a fracture intersection which is well connected to the flow system. The connectivity of the intersection is determined by checking that the exiting velocities from the intersection are above a specified minimum value. Particle movement is confined to an interior region of the domain to avoid boundary effects due to truncation of fractures and imposed boundary conditions as described by Long et. al. [1982]. 20 2.1.7 Movement Statistics in Discrete Sub-domain model Each time a particle enters a fracture segment a record is made of the direction, velocity and fracture length (distance between fracture intersections). From a single run of the discrete sub-domain model, the following statistical distributions are compiled: 1) Fracture length histogram (distance between intersections) for each set. 2) A frequency histogram cataloging the number of times particles moved in each of the 4 possible direction as indicated in Figure 2.2 3) Mean and standard deviation of the logarithm of velocities sampled for the 4 directions. From their initial work Schwartz and Smith [1988] deemed that these three sets of statistics were sufficient in capturing the nature of particle motion in the discrete network. As the particles will pass through only a sub-set of all the fractures forming the network, it is necessary to determine whether a histogram based on all the fracture lengths will introduce bias into the results. Comparisons of continuum model simulations, at various gradient angles, using two different fracture length histograms, (generated and sampled), showed no appreciable bias. Thus histograms based on all the fractures generated are used for all orientations of the hydraulic gradient in the continuum simulations, (Figure 2.3). At this point it should be noted that the nomenclature for fracture lengths is revised, from the total length of a fracture which may intersect several other fractures, to the distance between fracture intersections of fracture segments. Because of this the mean length of the fracture histograms will be much smaller than the value entered into the model. 21 2.1.8 Code Verification The code used here was modified from the original code of Schwartz and Smith [1988]. To check that there were no errors in the computer code simple networks of infinite length fractures were modeled and compared to analytical solutions (de Jong, 1972) . 2.1.9 Test for Representative Statistics As discussed in Chapter 1, the first assumption of the model is that the statistics collected are representative, in that the sub-domain model is greater than or equal to the size of the REV (Representative Elemental Volume) for solute transport. Smith and Schwartz [1988], found that a minimum of 900 fractures, (450 in each set), in a 10 meter by 10 meter sub-domain was necessary to collect REV statistics for the network geometry (fracture length, aperture and density) they simulated (Table 2.1). From Table 2.1 it can be seen that the two fracture sets are identical in nature, thus there will be a natural symmetry in the directional and velocity statistics. For example, in Table 2.2 for a gradient at 45° particles should be as likely to move in direction 2 as in direction 3, with equivalent mean velocities. There is also symmetry between gradient orientations of 0° and 90°, 15° and 75°, and 30° and 60°. Smith and Schwartz [1988] used this natural symmetry to obtain at least two values for each statistic collected from their single discrete sub-domain realization. By averaging the values, statistical noise due to the stochastic nature of the model is reduced. In this thesis with increased computing resources it is possible to run multiple realizations of the discrete sub-domain, the results of which can be combined to more effectively remove statistical noise. To determine how many realizations are necessary to remove the majority of the statistical noise, statistics from 50, 100 and 150 simulations of a given discrete network were compared. All input parameters for these realizations are given in Table 2.1. 22 Symmetry comparisons were made for the directional choice and mean velocity statistics. Standard deviations in velocity were included in Table 2.1 to give some idea of the range over which velocities vary for a given direction. From Table 2.2 there appears to be a slight improvement in overall symmetry between 50 and 100 realization, but little or no improvement between 100 and 150 realizations. From these comparisons it was decided that the compilation of statistics from 100 realizations is sufficient to remove the majority of the statistical noise for the network described in Table 2.1. The number of particles in the simulations will also affect how representative the statistics are. To determine how many particles are necessary to produce representative statistics, simulations using 250, 500 and 1000 particles for each of 100 MC (Monte Carlo) realizations were carried out. The results are tabulated in Table 2.3. Table 2.3 indicates a sizeable improvement in symmetry in going from 250 to 500 particles, but little improvement in going from 500 to 1000 particles, suggesting that 500 particles is a valid choice for the given sub-domain. Figure 2.3 graphically shows the symmetry in directional probability statistics for 100 realizations of 500 particles each. To run 100 realizations, of 500 particles each, for the discrete sub-domain described in Table 2.1, requires approximately 3.2 hours on a Sun SparcStation 4/330. The problem requires 5.5 Megabytes of memory. The computer limitations of the earlier work by Schwartz and Smith [1988] was related to charges for allocating large portions of RAM memory on a main-frame computer. A single realization of a 5 Megabyte discrete network simulation cost approximately $300 per run in 1987 (L. Smith, personal communication). At present up to 16 Megabytes of RAM are available on the Sun Workstation. 23 2.2 T h e C o n t i n u u m M o d e l Once the particle motion statistics have been collected from the discrete subdomain simulations, they are passed to the continuum model of the stochastic-continuum method. The stochastic-continuum model utilizes a continuum representation of the flow domain in which particles are moved by rules based on the discrete sub-domain statistics. Before applying the particle tracking method, the groundwater flow equation must be solved to determine the hydraulic head distribution within the domain. 2.2.1 Solution of the Groundwater Flow Equation To keep the coding simple, a regular grid, five point, finite difference model was chosen. Constant flux or constant head conditions are imposed on the boundaries and heads are solved for the interior nodes. Because of the regularity of the grid the bandwidth of the resulting matrix of linear equations is much smaller than for the discrete model, and thus the equations are much easier to solve. To simplify the computer coding and model verification an iterative SOR method was chosen to solve the linear equations (Huyakorn and Pinder, 1983). This solution method is different than that used by Smith and Schwartz [1988], but will produce equivalent results. The particle tracking method requires as input the magnitude and direction of the local hydraulic gradient. To determine the local hydraulic gradient each rectangle of the finite difference mesh is divided into two triangles. At any point within a triangle the magnitude and direction of the hydraulic gradient is constant, being determined from the corner head values. 2.2.2 Particle transport in the continuum model Particles are introduced into the continuum model as an instantaneous point source. The source position is chosen so as to avoid loss of solute out of the sides of the domain during the simulation. Particles are moved in steps along imaginary fractures, until the time period desired is completed. At each imaginary intersection the direction, 24 velocity and length of the next fracture step must be determined. The direction of particle motion must be sampled first as the velocity and fracture length histograms to sample from are delineated on the basis of the direction chosen. The directional probability distribution to sample from is determined by the orientation of the local hydraulic gradient. For gradient directions other than the seven orientations of the discrete sub-domain, probabilities on directional choices must be interpolated. As the particle is not allowed to travel back in the direction from which it last came, the histogram is resampled if this direction was chosen. Next the distance to a downgradient fracture intersection is sampled from the appropriate histogram for the given fracture set. Finally a log velocity is sampled for the given direction assuming a gaussian distribution. The mean and standard deviation for the distribution are interpolated from the values for the mean and standard deviation of the two nearest orientations of the hydraulic gradient for the seven values collected in the discrete subdomain. No memory is kept of the previous fracture steps taken by the particle. 2.3 Comparison of Full-Domain Discrete Model to Full-Domain Continuum Model (Uniform Fracture Aperture Case) To test whether the continuum model can mimic transport in a fracture network, a full-domain discrete model, twice the length of the sub-domain,(10 x 20 meters), was generated to permit direct comparison to a full-domain continuum model. The parameters for both full-domain models are listed in Table 2.4. To simplify comparisons both models were run with a uniform gradient parallel to the horizontal fracture set (0°). The execution times for the simulations are comparable, 15 minutes for the continuum and 12 minutes for the discrete model. These times include both the time for solving the flow equations and particle tracking. The advantage of the continuum model can be seen in the fact that it required only 1.5 megabytes of memory for 50000 particles, while the discrete model required 15.2 megabytes. 25 Schwartz and Smith [1988] found a reasonable match between the breakthrough curves for particles exiting across the down-stream boundaries of their discrete and continuum models. The computer generated test domain they used was equal in size to the discrete sub-domain described in Table 2.1. Because of computer limitations at the time, they were only able to simulate one realization of the smaller discrete network. Figure 2.4 shows breakthrough comparisons for the discrete network and continuum full-domain models. Even in the full-domain simulations, there can be considerable variation in breakthrough from one discrete realization to the next generated from the same statistical distribution. To represent an average breakthrough curve for the discrete network simulations, 10 discrete realizations were combined to form the discrete curve shown in Figure 2.4. Each of these simulations used 50,000 particles. The match is quite good for all but the tail of the plume and would suggest that the continuum model is doing a fairly good job of mimicking the discrete particle motion. The mismatch in the tail can be seen from a different perspective by comparing plots of the plume distribution for the continuum and discrete models. The network used for the discrete simulation described in the following is shown in Figure 2.5. Figure 2.6 is a comparison of one realization of the particle distribution in the discrete model to the continuum model, for a simulation time of 70 days. The arrow in the lower left hand corner indicates the direction of the applied hydraulic gradient. The asterisk in the middle of the left-hand side indicates the point injection location. To limit the number of particles exiting the domain to 1% of the total, a maximum simulation time of 70 days was chosen. To be able to represent the density as well as the location of the particles, the size of the points plotted in Figure 2.6 are made proportional to the number of particle falling within a square of dimension 0.05 meters The computer limitations relate only to the modeling of the discrete network. Computer burden for the actual continuum portion of the stochastic-continuum model is quite small. 26 around the center of the point. The discrete nature of the flow paths in the fracture network can be seen in the discrete plume plot, ss particles are restricted to actual fractures. Each realization of the discrete network will, therefore, produce a unique discrete plume. The shaded square below the injection point in the discrete plume plot indicates a large number of particles which have been travelling together in direction 3 at a very low velocity. These particles have not dispersed as they have not passed through many intersections. The continuum model cannot recover the exact distribution of an individual network plume, rather it produces an "average" plume distribution, similar to what we would get if we combined the plume distributions of multiple realizations of the discrete network. Plume statistics in the upper left hand corner of the plots record the first and second moments of the plume distribution. X is positive horizontally to the right and Z is positive vertically upwards on the plots. The origin of the coordinate axes is located in the middle of the left side boundary. X ^ and ZQ define the centroid of the plume with respect to the X,Z coordinate system, and are referred to as the first horizontal and first vertical moment of the plume, respectively. The second horizontal and second vertical moment S ^ x and S ^ , respectively, define the spread of the plume with respect to the XZ coordinate axes. The continuum model fails to reproduce the tail of the distribution seen in the discrete model. Histograms for the first and second moments calculated from 100 realizations of the discrete network are shown in Figure 2.7. The mean and standard deviation of the moments are noted in the upper right hand corner of each plot. The moments of the continuum plume in Figure 2.6 can be compared to the mean values of the moments in Figure 2.7, to determine how well the continuum model reproduces the average discrete fracture network plume. The first moments, X ^ and ZQ, for the continuum model (13.59, -0.06) show a reasonably good match to the first moments of the average discrete model (13.93, -0.01). The very small values ZQ indicate that there is no bias in the model with respect to motion up or down. This is as expected, as the 27 statistics for particle motion in directions 1 and 3, for a horizontal gradient across the fracture network, should be equivalent. The second horizontal moment, S^x> for the continuum plume (5.51) is below the mean value for the discrete simulations (7.40) by approximately one standard deviation (1.74) while the second vertical moment, S ^ , f ° r the continuum plume (1.17) is above the discrete value for the discrete simulations (0.69) by approximately 3 standard deviations (0.16). Because the magnitude of * s m u c n less than S ^ x the difference in value for the vertical moments is not as obvious from comparison of the particle plots in Figure 2.6. For hydrogeologic problems the difference in the magnitude of the second moments may be a more relevant comparison. l 1 In terms of magnitude the difference in S ^ x is 1.89 m over 13.6 m , while the difference in ^>jJL' 0-485 over 13.6 meters. Thus in terms of magnitude, the difference s in S - £ x IS approximately 4 times larger than the difference in S z z . The different nature of the development of dispersion over time for the two models can be observed in Figures 2.8a-f. Movement in the discrete model is limited to a finite number of flow paths defined by the fractures, creating an irregular distribution of the particles. The splitting of the solute at intersections can be observed by following the leading edge of the solute within individual fractures as time progresses. Dispersion in the continuum model creates a smooth distribution of the particles as there are an infinite number of possible flow paths. Figures 2.9a,b show the corresponding development of the first and second moments of the plumes, the discrete values being the mean of 100 discrete network simulations. The underestimate of S x x °y the continuum model is shown in Figure 2.9b. The growth of S ^ x f ° r t n e discrete network plume can be divided into three sections. The initial section from 0 to 30 days shows non-linear growth in the near-source region. The region through which the solute has moved up to this point is sub-REV. From 30 days to 60 days the growth in S x X 1S linear as would be expected as the region through 28 which the solute has passed is greater than the REV for solute transport. After 60 days a decrease in the growth of S X Y i is caused by the loss of particles out of the right-hand side of the domain. Because of this loss analysis of results should be confined to a maximum of 60 days for these simulations. Comparing the three methods of analysis, (breakthrough curves, particle distribution plots, and moment development plots), it can be seen that each provides complementary information. For the question of when the contaminant will reach a drinking well the start of the breakthrough curves is the most important piece of information. If, however, remediation of the aquifer by pumping out the contaminant is proposed, then the distribution of the plume is important and the particle plots or moment diagrams are needed. In the first case, in determining when the contaminant will reach the well, the model of Smith and Schwartz [1988] does a sufficient job. In the second case, however, the failure to match the long tail of the plume may cause an under-estimation of the time needed to pump out all the contaminant. The next chapter deals with the question of why the continuum model fails to reproduce the long tail of the discrete model. The physical mechanisms behind the long tail are determined, and statistics which capture these mechanisms are added to the basic model. 29 Table 2.1 Summary of input parameters for the discrete sub-domain model used in Chapter 2. Parameter (nr) (nc) (AX) (AZ) number of rows number of columns horizontal grid spacing vertical grid spacing Value Background Grid Fracture System (nhorzf) number of horizontal fractures (nvertf) number of vertical fractures ^hfP horizontal fracture length (/x «) mean vertical fracture length (minhfl) minimum horiz. fracture length (minvfl) minimum vert, fracture length (^ ) log mean horiz. fracture aperture (/i p ) log mean vert, fracture aperture (t7 ) log std. dev. horiz. fracture aperture (a j ) log std. dev. vert, fracture aperture m e a n y aph a V a p h ap 1 General Parameters (mc) number of monte carlo realizations (np) number of particles (At) time step size (Ah/AX) gradient in direction of flow 400 400 0.025 m 0.025 m 450 450 0.4 m 0.4 m 1.20 m 1.20 m -4.5 m -4.5 m 0.0 m 0.0 m 50/100/150 250/500/1000 25 days 0.01 30 Table 2.2 Comparison of directional and velocity statistics for 50, 100 and 150 simulations of the discrete sub-domain model. 500 particles are used in each simulation. 50 SIMULATIONS 0 15 30 45 60 75 90 1 2 3 4 0.12 0.76 0.11 0.00 0.05 0.72 0.23 0.00 0.02 0.62 0.36 0.00 0.01 0.49 0.50 0.01 0.00 0.36 0.62 0.02 0.00 0.23 0.72 0.05 0.00 0.10 0.78 0.12 Mean Log Velocity 1 2 3 4 3.75 3.37 3.76 3.85 3.78 3.38 3.68 4.17 3.83 3.42 3.58 4.01 3.90 3.48 3.48 3.93 3.97 3.58 3.43 3.83 4.05 3.68 3.38 3.79 4.28 3.77 3.36 3.75 Std. Dev. Log Velocity 1 2 3 4 0.34 0.14 0.35 0.06 0.35 0.15 0.29 0.05 0.33 0.17 0.24 0.20 0.27 0.19 0.20 0.30 0.20 0.24 0.17 0.33 0.13 0.29 0.15 0.35 0.00 0.33 0.14 0.33 Directional Choices % 100 SIMULATIONS 0 15 30 45 60 75 90 Directional Choices 1 2 3 4 0.12 0.77 0.12 0.00 0.04 0.72 0.24 0.00 0.02 0.62 0.37 0.00 0.01 0.49 0.40 0.01 0.00 0.36 0.62 0.02 0.00 0.23 0.73 0.05 0.00 0.11 0.78 0.11 Mean Log Velocity 1 2 3 4 3.76 3.36 3.76 4.10 3.77 3.38 3.68 4.17 3.82 3.42 3.57 3.99 3.91 3.48 3.48 3.90 3.99 3.58 3.42 3.81 4.16 3.68 3.38 3.79 4.18 3.77 3.36 3.77 Std. Dev. Log Velocity 1 2 3 4 0.34 0.14 0.34 0.04 0.34 0.15 0.29 0.12 0.33 0.16 0.24 0.19 0.28 0.19 0.20 0.29 0.19 0.23 0.17 0.32 0.12 0.29 0.15 0.34 0.00 0.33 0.14 0.33 45 60 75 90 % ISO SIMULATIONS 0 15 1 2 3 4 0.12 0.77 0.12 0.00 0.04 0.72 0.23 0.00 0.02 0.61 . 0.37 0.00 0.01 0.46 0.49 0.01 0.00 0.37 0.62 0.02 0.00 0.23 0.73 0.05 0.00 0.11 0.77 0.12 Mean Log Velocity 1 2 3 4 3.76 3.37 3.76 4.23 3.78 3.38 3.68 4.18 3.83 3.42 3.67 4.01 3.92 3.50 3.48 3.92 3.99 3.58 3.42 3.82 4.13 3.69 3.38 3.79 4.21 3.77 3.37 3.77 Std. Dev. Log Velocity 1 2 3 4 0.34 0.14 0.34 0.05 0.34 0.15 0.29 0.12 0.33 0.17 0.24 0.20 0.27 0.19 0.19 0.27 0.20 0.24 0.17 0.33 0.12 0.29 0.15 0.34 0.06 0.33 0.14 0.33 Directional Choice* 96 30 31 Table 2.3 Comparison of directional and velocity statistics for 250, 500 and 1000 particles in the discrete sub-domain. The results are the compilation of 100 realizations. 250 PARTICLES 0 15 30 45 60 75 90 1 2 3 4 0.12 0.77 0.11 0.00 0.04 0.73 0.23 0.00 0.02 0.62 0.37 0.00 0.00 0.49 0.50 0.01 0.00 0.36 0.62 0.01 0.00 0.23 0.73 0.04 0.00 0.11 0.78 0.12 Mean Log Velocity 1 2 3 4 3.76 3.37 3.75 4.01 3.81 3.38 3.68 3.98 3.82 3.42 3.57 3.99 3.88 3.49 3.48 3.86 4.02 3.58 3.42 3.84 4.19 3.69 3.38 3.79 4.16 3.78 3.37 3.76 Std. Dev. Log Velocity 1 2 3 4 0.33 0.14 0.33 0.04 0.34 0.15 0.29 0.08 0.32 0.16 0.24 0.19 0.27 0.19 0.19 0.26 0.20 0.24 0.16 0.32 0.09 0.29 0.15 0.34 0.04 0.34 0.14 0.34 Directional Choices % . 500 PARTICLES 0 15 30 45 60 75 90 1 2 3 4 0.12 0.77 0.12 0.00 0.04 0.72 0.24 0.00 0.02 0.62 0.37 0.00 0.01 0.49 0.50 0.01 0.00 0.36 0.62 0.02 0.00 0.23 0.73 0.05 0.00 0.11 0.78 0.11 Mean Log Velocity 1 2 3 4 3.76 3.36 3.76 4.10 3.77 3.38 3.68 4.17 3.82 3.42 3.57 3.99 3.91 3.48 3.48 3.90 3.99 3.58 3.42 3.81 4.16 3.68 3.38 3.79 4.18 3.77 3.36 3.77 Std. Dev. Log Velocity 1 2 3 4 0.34 0.14 0.34 0.04 0.34 0.15 0.29 0.12 0.33 0.16 0.24 0.19 0.28 0.19 0.20 0.28 0.19 0.23 0.17 0.32 0.12 0.29 0.15 0.34 0.01 0.34 0.14 0.33 Directional Choices % 1000 PARTICLES 0 15 30 45 60 75 90 1 2 3 4 0.12 0.77 0.12 0.00 0.05 0.71 0.24 0.00 0.01 0.62 0.37 0.00 0.00 0.49 0.50 0.00 0.00 0.36 0.62 0.01 0.00 0.23 0.73 0.04 0.00 0.12 0.77 0.11 Mean Log Velocity 1 2 3 4 3.76 3.37 3.76 4.27 3.81 3.38 3.68 4.24 3.85 3.42 3.58 4.05 3.96 3.49 3.50 3.94 4.06 3.58 3.44 3.84 4.12 3.68 3.38 3.81 4.27 3.75 3.37 3.78 Std. Dev. Log Velocity 1 2 3 4 0.34 0.14 0.34 0.07 0.34 0.15 0.29 0.12 0.31 0.16 0.23 0.19 0.27 0.19 0.19 0.28 0.17 0.23 0.16 0.33 0.13 0.29 0.15 0.34 0.04 0.33 0.14 0.34 Directional Choices % 32 Table 2.4 Summary of input parameters for full-domain models used in Chapter 2. DISCRETE MODEL Parameter (nr) number of rows (nc) number of columns (AX) horizontal grid spacing (AZ) vertical grid spacing Value Background Grid Fracture System (nhorzf) number of horizontal fractures (nvertf) number of vertical fractures (/ijjfj) mean horizontal fracture length (/i f|) mean vertical fracture length (minhfl) minimum horiz. fracture length (minvfl) minimum vert, fracture length (/i ) log mean horiz. fracture aperture (^ p ) log mean vert, fracture aperture (a ) log std. dev. horiz. fracture aperture (<7 ) log std. dev. vert, fracture aperture y apn a V apn A P J L General Parameters (np) number of particles (Ah/AX) gradient in direction of flow 400 800 0.025 m 0.025 m 904 904 0.4 m 0.4 m 1.20 m 1.20 m -4.5 m -4.5 m 0.0 m 0.0 m 50000 0.01 CONTINUUM MODEL Parameter (nr) number of rows (nc) number of columns (AX) horizontal grid spacing (AZ) vertical grid spacing Value Background Grid General Parameters (np) number of particles (Ah/AX) gradient in direction of flow 100 200 0.1 m 0.1 m 50000 0.01 a) uncleaned network 0 5 meters 10 b) cleaned network Figure 2.1 One realization of the discrete sub-domain fracture network. 34 F i g u r e 2.2 T h e discrete network f l o w d o m a i n , showing boundary c o n d i t i o n s , gradient orientations and the nomenclature f o r f l o w directions, f r o m Schwartz and S m i t h [1988]. 0.25 "1 1 HORIZONTAL SET i ~] i 1 1 1 1 1 1 1 1 1 1 r -i 1 r -i i 0.2 13 o *S 0-15 CO no cd 0.1 o cu ft 0.05 0 L 1 1 1 0 1 1 0.2 i i 1 0.4 0.6 fracture length 0.25 i 1 r n VERTICAL SET r 1 1 r T 0.8 (m) ~i 1 r 0.2 ca -t-> o 0.15 0.4 0.6 fracture length 0.8 i_ (m) Figure 2.3 F r a c t u r e length histograms for the network described i n Table 2.1 (network A). Time (days) Figure 2.4 B r e a k t h r o u g h comparisons for the discrete and continuum f u l l — d o m a i n u n i f o r m aperture simulations. The discrete curve is the average of 10 simulations of the discrete f u l l — d o m a i n model. 0 5 10 meters 15 Figure 2.5 The discrete full-domain fracture network used for the simulations in Chapter 2. 20 DISCRETE MODEL Plume Stats. xc= zc= sxx= 12.74 0.50 11.69 70 days — > 0° CONTINUUM MODEL Plume Stats. xc= 13.59 0 5 10 meters 15 20 Figure 2.6 Particle distribution plots for the discrete and continuum full-domain models at T=70 days. Moment statistics of the plume are listed in the top left-hand corner of the plot. Simulation parameters are listed in Table 2.4. 1st HORIZONTAL MOMENT (XC) 12 13 14 15 16 2nd HORIZONTAL MOMENT (SXX) 5 10 15 1st VERTICAL MOMENT (ZC) -1 -0.5 0 0.5 1 2nd VERTICAL MOMENT (SZZ) 0.4 0.6 0.8 1 1.2 1.4 Figure 2.7 Histograms of the first a n d second spatial moments of the discrete f u l l — d o m a i n plume compiled f r o m 100 Monte Carlo simulations at T=70 days. D I S C R E T E MODEL Plume xc= zc= sxx= szz= Stats. 1.61 0.14 0.57 0.13 10 c\ays 0° — 7 " CONTINUUM MODEL Plume xc= zc= sxx= szz= Stats. 2.06 -0.01 0.58 0.15 10 clays •> I 1 0 I 1 0° 1 1 ! 1 ! 1 I 1 5 ! 1 1 1 1 1 1 1 1 1 1 1 1 1 i 1 1 i i ' I 1 i 1 i 1 1 i i I 10 15 20 meters Figure 2.8a Particle distribution plots for one realization of the discrete and continuum full-domain models at T=10 days. DISCRETE MODEL Plume Stats. xc= 3.41 zc= 0.30 sxx= 1.90 szz= 0.17 20 (Jays — 7 o° CONTINUUM MODEL Plume Stats. xc= 3.99 zc= -0.02 sxx= 1.38 szz= 0.32 20 days — 0 > 0 ° 5 10 meters 15 Figure 2.8b Particle distribution plots for one realization of the discrete and continuum full-domain models at T=20 days. 20 DISCRETE Plume xc= zc= sxx= szz= MODEL Stats. 5.21 0.31 3.46 0.32 30 c\ays — 7 0° C O N T I N U U M MODEL Plume xc= zc= sxx= szz= Stats. 5.92 -0.02 2.21 0.48 30 days — 0 > 0 ° 5 10 meters 15 Figure 2.8c Particle distribution plots for one realization of the discrete and continuum full-domain models at T=30 days. 20 DISCRETE MODEL Plume Stats. xc= 6.99 zc= 0.36 sxx= 5.57 szz= 0.33 40 ctays — 7 " 0° CONTINUUM MODEL Plume Stats. xc= 7.84 zc= -0.03 sxx= 3.05 szz= 0.65 40 (Jays — 7 0 0 ° 5 10 meters 15 Figure 2.8d Particle distribution plots for one realization of the discrete and continuum full-domain models at T=40 days. 20 DISCRETE MODEL Plume Stats. xc= 8.89 zc= 0.45 sxx= 7.79 50 days — 7 " 0° CONTINUUM MODEL Plume Stats. xc= 9.76 0 5 10 meters 15 Figure 2.8e Particle distribution plots for one realization of the discrete and c o n t i n u u m f u l l - d o m a i n models at T=50 days. 20 DISCRETE MODEL Plume, Stats. xc= 10.83 zc= 0.51 sxx= 9.87 60 cteys —7 o° CONTINUUM MODEL Plume Stats. xc= 11.69 0 5 10 meters 15 Figure 2.8f Particle distribution plots for one realization of the discrete and c o n t i n u u m f u l l - d o m a i n models at T=60 days. 20 46 1st Horizontal Moment (Xc) i r i r discrete continuum 10 -t-3 0 0 20 40 60 1st Vertical Moment (Zc) "i CO u 1 CD 0 O 1 r i 1 r i L -1 -2 i 0 20 40 time (days) 60 Figure 2.9a First spatial momenta (Xc, Zc) for the continuum and discrete, fulldomain uniform aperture simulations. The discrete curves are the average of 100 realizations of the discrete network. 2nd Horizontal Moment (Sxx) 10 n —i r 1 1 1 1 1 1 1 1 1 1 r discrete continuum 8 6 CM X X Vi 2 - 1 I 1 I 0 1 1 20 1 i 60 40 2nd Vertical Moment (Szz) 0 - 0 Illllllllllll llllll N i t—-—r i 1 i i i 1 i 1 i i i 1 i 20 l i 1 - i 40 time (days) 1111U1111111 0 60 Figure 2.9b Second spatial moments (Sxx, Szz) for the continuum and discrete, full—domain uniform aperture simulations. The discrete curves are the average of 100 realizations of the discrete network. 48 3.0 E X T E N S I O N S T O 2 - D M O D E L As discussed in Chapter Two, the basic stochastic-continuum model does not reproduce the long tail of the solute plumes found in the discrete full-domain simulations. As shown in the next subsection the differences in the second moments between the discrete and continuum models become even larger when aperture variability is added to the models. 3.1 Variable Aperture Model Table 3.1 lists the input parameters for a variable aperture, discrete sub-domain model. The domain is 10m by 10m. The standard deviation in logarithmic fracture aperture is 0.2. In this case, 95% of the fracture apertures created will fall between 13 and 78 u with a mean aperture of 32 fi. As was done in the uniform aperture model, movement statistics are collected for 7 angles of the hydraulic gradient with respect to • the orientation of the fracture network. By introducing aperture variability in the model a greater range in flow path velocities has been created. Because of this increased variability in velocity the number of simulations necessary to remove the statistical noise in the particle motion statistics will increase. The results of a symmetry analysis suggests that 200 simulations of 500 particles each will produce representative statistics. As in the uniform case a larger scale discrete model, 10m by 20 m, was created to compare with the full-domain continuum model, (Table 3.2). Figure 3.1 shows particle plots at T=35 days for both the discrete and continuum models. As can been seen in the discrete plot, the centroid of the plume moves approximately 1.4 times as far in 35 days compared to the uniform aperture simulations for the same time period (Figures 2.8c and d). The overall higher velocity of the plume is due to the fact that the majority of the solute is channeled along the 49 higher velocity paths, which due to the variability in aperture have even higher velocities. The spread of the plume has also greatly increased. For the same time period the uniform aperture fracture networks gave an average second horizontal moment of 3.08 m compared to an average of 9.2 m' for the variable aperture case. Figure 3.2 shows the breakthrough curves for the two plumes illustrated in Figure 3.1. Figure 3.3 shows the breakthrough curve for an average of 10 discrete realizations compared to the breakthrough curve for the continuum model. From Figure 3.3, it can be seen that the continuum model fails to predict both the initial breakthrough and the long tail of the plume. The temporal development of the first and second spatial moments is shown in Figures 3.4a and b. The difference in the spread of the plume is quite apparent from the second horizontal moments ( S ^ x ) plotted in Figure 3.4b. Figure 3.4a shows that although the continuum model does a poor job in reproducing the second spatial moments, it does capture the centroid (X^, ZQ) of the plume. 3.2 An Improved Stochastic-Continuum Method Two factors, both relating to the model used to select particle statistics, were found to account for the difference between the discrete and continuum plumes. The first has to do with the assumption of a gaussian distribution of particle velocities and the second with correlations between successive velocities as solutes move through a discrete fracture network. 3.2.1 Velocity Distributions Figure 3.5a shows histograms on particle velocity for the uniform aperture case, while Figure 3.5b shows similar plots for the variable aperture case. The results are a compilation of 100 simulations. The number of velocities sampled for all but direction 4, which is rarely traveled in, is quite high as noted at the top of the plots. The histograms are normalized by the total number of samples, making the area under each 50 curve equal to 1. The logarithmic velocities are multiplied by -1 so that the plot will have a positive ordinate, thus velocities are decreasing to the right. Overlain on the histograms are the best fit gaussian and three-parameter gamma (a,A,a) distributions. Gamma densities provide a fairly flexible class for modeling non-negative random variables (Rice, 1988). The a parameter defines the shape of the distribution, while the A parameter determines the scale of the distribution, without changing its' shape. The a parameter determines the offset of the distribution from the frequency axis (abscissa). The two curves can be identified by their characteristic shapes, the gaussian curve is symmetric while the gamma curve is skewed to the right. The parameters for the threeparameter gamma distribution are calculated from the mean, standard deviation and skew of the histogram data: a = Hh) 1 (3-D A = (a)- /<7 (3.2) 5 a = n - a/A (3.3) where: a is the estimated mean a is the estimated standard deviation 7 is the estimated skew Ignoring direction 4 for the moment, it can be seen that the gamma distribution fits the histograms better than the gaussian distribution for all cases except direction 2 in the uniform aperture case, for which a gamma distribution could not be fitted. Sampling from a gamma distribution for velocities is much like sampling from a gaussian distribution (Press et al. 1988). The improvement in the moment statistics by sampling particle velocities from the histogram or from a fitted gamma distribution, as compared to sampling from a gaussian distribution, can be seen in comparing Figures 3.6a,b and 3.7a,b, for the uniform and variable aperture cases, respectively. The curve marked 51 gaussian(O) is the original continuum curve (Figures 2.9a,b and 3.4a,b) created by sampling from a gaussian distribution. The curves labeled histogram and gamma(O) are the results of sampling directly from the histogram and from a fitted gamma distribution, respectively. Comparisons of the gaussian(O), gamma(O) and histogram curve indicate that the gamma distribution represents the histogram data better than the gaussian distribution from the perspective of a transport simulation. The effect of whether the velocity histogram in direction 4 is represented by a gaussian or gamma distribution is negligible as it is sampled infrequently. For the uniform aperture case, the difference in S ^ x between the continuum and discrete simulated plumes is decreased significantly by using a gamma distribution, (Figure 3.6b). However, the difference in S x x f ° r variable aperture simulations is still quite large, (Figure 3.7b), suggesting that there is still something missing from the statistical-continuum model. The good fit of the gaussian distribution to the uniform aperture velocity histogram in direction 2, which is the principal direction of motion for these simulations, shows why Schwartz and Smith [1988] found relatively good matches between their discrete and continuum breakthrough curves as their simulations were done for a uniform aperture network. 3.2.2 Correlations in Particle Velocities Analysis of the velocity statistics from the discrete sub-domain model indicates that the velocities in successive fractures segments along a flow path are correlated. Correlations between consecutive fracture segment flow velocities when the flow directions are different, (ie. going from fracture segment 2 to fracture segment 3 in Figure 3.8), were found to be small, thus calculations of correlation in the discrete subdomain were restricted to velocities in similar direction. This analysis of crosscorrelation was, however, restricted to simulations in which the hydraulic gradient was 52 kept parallel to the horizontal fracture set. As will be discussed in Chapter 4 crosscorrelations may be important for other orientations of the hydraulic gradient. Two correlation models, referred to as the short model and the long model, were investigated. The short model calculates the correlation between velocities in directly connected successive fracture segments, while the long model does not restrict the correlations to velocities in directly connected fracture segments. For the travel path shown in Figure 3.8, the short model will record a lag(l) correlation between the velocities in the first and second fractures segments only, while the long model will record a lag(l) correlation between the first and second and the second and fourth fractures segments. The first four correlation coefficients, (lag(l) - lag(4)), are calculated for both models. For the example above, the long model would record the correlation between the velocity in the first and fourth fracture segments as a lag (2) correlation, while the short model would not record any lag (2) correlations. The fact that the short model will not record velocity correlations over intermediate non-parallel travel segments limits the effect of the correlation to a fairly short distance, as most flow paths sample frequently from both fracture sets. The long model more realistically extends the correlation over these intermediate steps, thus extending the effect of the velocity correlation. The correlation values for the short model are typically higher than the values for the long model, as shown in Table 3.3 for both the variable aperture and uniform aperture cases. Due to the increased channelization of the flow, the particle velocity correlations in the variable aperture simulations are higher than those in the uniform aperture simulations. Correlation values for direction 4 were set to zero as not enough samples were collected for this direction. The results in Table 3.3 are from 100 realizations of the discrete network model as described in Tables 1.1 and 2.1. To apply the correlations in the continuum model, standard auto-regressive (AR) models for gaussian random variables (Box and Jenkins [1976]) were converted by a 53 probability transformation to a gamma distribution (Rubenstein [1981]). The reason for the two step process is to avoid the approximations in auto-regressive gamma models such as the Wilson-Hilferty transformation (Lettenmaier and Burges, 1977), and the method used by Gaver and Lewis [1980]. The process is somewhat slower than the approximate methods but with the speed of modern computers is not prohibitive. For the lag one auto-regressive gaussian model (AR(1)) the velocity at time step t (V ) is t conditioned on the velocity of the previous time step (V^ _jp. If we allow Z to be the t deviation in velocity from the mean velocity (/i): Z =V- u (3.4) then: Z= #l)Z _ t (t #1) a 2 a 1} +a t (3.6) P O ) = (T (l-#1) ) 2 (3.5) 2 2 (3.7) where: a is a Normal (gaussian) random deviate t N(0,<7 ) A p(l) is the lag(l) correlation coefficient To convert to a gamma (a,A,a) distribution the following transformation is applied, (Rubenstein [1981]): X = F (»(Z /a )) + a _1 t z where: X is the gamma (a,A,a) distributed AR(1) particle velocity [L/T] F~* is the gamma (a,A) inverse Cumulative Distribution Function (CDF) $ is the standard normal CDF Examples of numerical models for the inverse gamma CDF can be found in Becker and Chambers [1984]. (3.8) 54 Figures 3.9a and b show the moments of the particle distribution for both the short and long auto-regressive gamma models compared with the discrete and gamma(O) curves for the variable aperture case. The better match in S ^ x °f t n e gamma(long) curve to the discrete curve is due to the ability of the long model to carry correlations over a longer distance. As a result of these comparisons, the long model was considered superior to the short model, and will be used exclusively for the remainder of this thesis. The equation for the AR(2) model is slightly more complicated: Z = #l)Z _ t (t 1} + * 2 ) Z _ + (1 - pOWl) - p(2)tf2))a (t 2) (3.9) t (KD = (p(l)(l-p(2))/(l-p(l)) (3.10) <K2) = (p(2)-p(l) )/((l-p(l) ) (3.11) cr = r (l-u(l)p(l)-c>(2)p<2)) (3.12) 2 2 2 a 2 2 z Figures 3.10a and b shows the first and second moments of the plume for the AR(1) and AR(2) models, with the discrete and gamma(O) curves for reference, for the variable aperture case. The AR(1) and AR(2) curves are designated gamma(l) and gamma(2), respectively. A slight improvement in the match of the gamma(2) S x x curve to the discrete curve over the gamma(l) curve can be seen. As a result of the comparisons documented in this chapter a new twodimensional stochastic-continuum model equivalent to the gamma(2) model is proposed. Figures 3.11a,b summarize the improvements in the fit of the continuum model spatial moments to the spatial moments of the discrete model. The particle plots for the gamma(2) and discrete model are shown in Figure 3.12. Figure 3.12 visually shows an improvement in the distribution of the plume over the gaussian(O) model shown in Figure 3.1. It should be kept in mind that the discrete plume shown is the product of only one realization of the discrete network. Figure 3.13 shows the distribution of the first and second spatial moments for the 100 discrete full-domain simulations, at T=35 days. Comparison of Figures 3.12 and 3.13 shows that the discrete plume simulated in 55 Figure 3.12 has above average first horizontal moments (9.65 m compared to the average discrete value of 9.82 m for X^) and below average second horizontal moments (8.38 m compared to the average discrete value of 8.90 m for S ^ x ) - Thus the discrete plume has travelled further but is less spread out than the average discrete plume, which the continuum model is designed to simulate. Figure 3.14 shows an improvement in the match of the continuum model breakthrough to the discrete breakthrough curve by going to a gamma(2) model. The discrete curve is the average of 10 discrete simulations. The fact that the spatial moment curves shown in Figure 3.11a,b are linear after an initial time allows us to describe the difference between the continuum and discrete curves as a percentage of the discrete value. A slight downwards deviation in the discrete curve after 30 days is due to the loss of a small proportion of the particles out the right side of the domain. Due to their linear nature both the first and second spatial moment curves can be represented by first order polynomials of the form: Y =A T d Y d c = c A T B + + B d c (3.13) <- > 3 14 where: T is the time variable Y is the spatial moment O^Q^C^XX^ZZ^ c,d are subscripts indicating the continuum or discrete curve Table 3.4 lists the values of A and B for the discrete and gamma(2) spatial moment curves. The first vertical moment, (ZQ), is not included as it should be equal to zero for all time. Any variation in ZQ in the data is due to statistical noise. As discussed earlier the magnitude of the difference in the second moments is a more meaningful comparison than the percentage differences. By adding the two second moment curves together a more meaningful percentage comparison can be made. This combination proves quite useful for comparing second moments between simulations with different orientations of the hydraulic gradient in the next chapter. The combination of the two second moment curves is made so that the differences between the continuum and discrete moments add together. By combining Equations 3.13 and 3.14 the continuum curves can be expressed in terms of the discrete curve. Y = A / A + (B - (A /A )B ) c c d c c d (3.15) d By substituting the values of A and B given in Table 3.4 into Equation 3.15 we get X S C c = 0.94 X XXc = 0 8 4S S Zc - ° S + ZZc> = S 0 8 Values of A / A for ( S ^ x ^ZZ^ + c d a + 0.31 + ZZd + 7 1 S 3 d XXd Z < XXc C ( XXd S r e + 0 1 °S (3.16) <- > 2 3 17 <- > 1 2 ZZd) 3 18 + 0 - (- ) 2 2 3 19 calculated to give the maximum combined difference. As time increase the constant terms become negligible and the far-field percentage difference between the moment curves asymptotically approaches a value of (1-A /A ), or (1-C) where C is the first coefficient in Equations 3.16-3.19. Thus the c d far-field difference between the gamma(2) and discrete models for XQ, S X > S X Z Z , A N C * ( S + S ) are 6%, 16%, 29% and 17%, respectively. Values of the far-field percentage xx zz difference for the first and second spatial plume moments for the continuum models tested in this chapter as compared to the discrete model are tabulated in Table 3.5. The improvement in the far-field percentage difference for (Sxx ^ZZ^ ^ + r o mt n e gaussian(O) (57%) to the gamma(2) (17%) model is quite significant. The poorer match for X c using the gamma(2) model (6%), compared to 0% for the gaussian(O) model is within reasonable error. Questions arise as to what is still missing from the gamma(2) model which would further improve the match of the mass distribution to the discretely modeled distribution. Possible causes of the differences, which were not investigated to any great length, include cross-correlation of the velocities in different directions, correlations between directional choice and velocity, and correlations between fracture length and velocity. The next chapter will test how well this new model simulates solute transport fractured media, for various networks aperture distributions, and orientations of the hydraulic gradient. 58 Table 3.1 S u m m a r y o f input parameters f o r discrete s u b - d o m a i n model w i t h variable fracture apertures. Parameter Value Background Grid (nr) number of rows (nc) number of columns (AX) horizontal grid spacing (AZ) vertical grid spacing 0.025 0.025 Fracture System (nhorzf) number of horizontal fractures (nvertf) number of vertical fractures (/i fj) mean horizontal fracture length (/i £j) mean vertical fracture length (minhfl) minimum horiz. fracture length (minvfl) minimum vert, fracture length (/i ) log mean horiz. fracture aperture (/i ) log mean vert, fracture aperture (cr p ) log std. dev. horiz. fracture aperture ( aph) ' * f aperture R y aph apv a a n o gs t d d e V -v e r t > 400 400 r a c t u r e General Parameters (mc) number of monte carlo realizations (np) number of particles (At) time step size (Ah/AX) gradient in direction of flow m m 450 450 0.4 m 0.4 m 1.20 m 1.20 m -4.5 m -4.5 m 0.2 m 0.2 m 200 500 15 days 0.01 59 Table 3.2 Summary of input parameters for full-domain models with variable fracture apertures. DISCRETE MODEL Parameter Value Background Grid (nr) number of rows (nc) number of columns (AX) horizontal grid spacing (AZ) vertical grid spacing 400 800 0.025 m 0.025 m Fracture System (nhorzf) number of horizontal fractures (nvertf) number of vertical fractures ^hfl) horizontal fracture length (/i fj) mean vertical fracture length (minhfl) minimum horiz. fracture length (minvfl) minimum vert, fracture length (/i jj) log mean horiz. fracture aperture (/i ) log mean vert, fracture aperture ( aph) *" ^ " horiz. fracture aperture ^aph) '°*» * ^ ' f aperture m e a n v ap apv a stc st< ev - e v v e r t - r a c t u r e General Parameters (np) number of particles (Ah/AX) gradient in direction of flow 904 904 0.4 m 0.4 m 1.20 m 1.20 m -4.5 m -4.5 m 0.2 m 0.2 m 50000 0.01 CONTINUUM MODEL Parameter (nr) number of rows (nc) number of columns (AX) horizontal grid spacing (AZ) vertical grid spacing Value Background Grid 100 200 0.1 m 0.1 m General Parameters (np) number of particles (Ah/AX) gradient in direction of flow 50000 0.01 Table 3.3 Summary of the first two correlation coefficients for the short and long correlation models from the uniform and variable fracture aperture discrete sub-domain simuations. UNIFORM APERTURE Long Correlation Model Dir 1 Dir 2 Dir 3 Dir 4 r(l) 0.10 0.51 0.12 0.00 r(2) 0.02 0.31 0.01 0.00 Short Correlation Model Dir 1 Dir 2 Dir 3 Dir 4 r(l) 0.34 0.64 0.34 0.00 r(2) 0.12 0.46 0.13 0.00 VARIABLE APERTURE Long Correlation Model Dir 1 Dir 2 Dir 3 Dir 4 r(l) 0.42 0.68 0.39 0.00 r(2) 0.25 0.51 0.22 0.00 Short Correlation Model Dir 1 Dir 2 Dir 3 Dir 4 r(l) 0.60 0.80 0.65 0.00 r(2) 0.48 0.69 0.51 0.00 Note: Correlation values for Direction 4 have been arbitrarily set to zero as insufficient statistics are collected for this direction. 61 Table 3.4 Coefficients for first-order polynomial representation of spatial moment curves for the discrete and gamma(2) simulations. Model X C s xx s zz S XX ZZ +S Discrete A 0.278 0.326 0.042 0.368 Discrete B 0.170 -2.190 -0.100 -2.290 Gamma(2) A 0.262 0.274 0.030 0.304 Gamma(2) B 0.470 -1.720 0.050 -1.670 Table 3.5 Summary of the far-field percentage differences in the first and second spatial moments for the continuum models described in chapter 3 compared to the average discrete model. Continuum Model AX (%) Gamma(2) 6 16 29 17 Gamma(long) Gamma(l) 5 23 29 23 Gamma(short) 6 39 19 36 Gamma(O) 6 44 29 42 Histogram 19 45 43 45 Gaussian(O) 0 61 19 57 r AS V V ASw AS -v+AS Y 63 DISCRETE MODEL Plume Stats. xc= 11.05 ZC= 0.54 • • V " " : - 35 clays — 7 " 0° CONTINUUM MODEL Plume Stats. 0 5 10 meters 15 Figure 3.1 Particle distribution plots for variable aperture discrete and c o n t i n u u m f u l l - d o m a i n , gaussian(O), models at T=35 days. 20 Time (days) Figure 3.2 Breakthrough curve f o r the plumes shown i n Figure 3.1 Time (days) Figure 3.3 Breakthrough curves for the discrete and continuum f u l l domain, variable aperture simulations. The discrete curve is the average of 10 simulations of the discrete f u l l domain model. 66 1st Horizontal Moment (Xc) i i r i r 1 1 1 1 r i r L J I I L J L - — —r 1 discrete continuum gaussian(O) 10 - CD 0 J I I L 0 J I I 10 30 20 1st Vertical Moment (Zc) 2 OT -4-> a i i i r I i i i i i I I I !_ -i i : r j I I L 1 0 - 1 fcJ -2 0 J 10 20 30 time (days) Figure 3.4a Development of the first moments (Xc.Zc) for the variable aperture, discrete and continuum models. The discrete curves are the average of 100 realizations of the discrete network. 10 2nd Horizontal Moment (Sxx) I i i i i | i i i i | i i i i | i time (days) Figure 3 . 4 b Development of the second moments (Sxx.Szz) for the variable aperture, discrete and continuum models. The discrete curves are the average of 1 0 0 realizations of the discrete network. r DIRECTION 2 DIRECTION 1 2 2 3 9 4 8 s a m p l e s T—i—i—i—|—i—i—i—i—|—i—i—i—i—|—i—i—i—i—|—n 0.06 gaussian ^ = 3 7 5 1.50889e+06 s a m p l e s T—I—I—I—|—I—I—I—I—|—I—I—I—I—|—I—I—I—I—|—I yu= - 0.04 0.64 3.37 a= <r= 0.34 7= - 0.15 0.1 0A5' 7= 0.2f gaussian 0.05 0.02 j i i i I i i i 3 -log(velocity) ( m / s ) DIRECTION 3 228616 samples DIRECTION 4 0.25 y.= - 1 i i i 4 i 1 l i i i 5 I i_ 6 —log(velocity) ( m / s ) ~i—i—i—I—1—i—i—i—i—I—i—i—I—I—1—i—I—I—I—I—r 0.06 i 331 samples i—i—i—i—i—i—i—i—i—i—i—i—i—i—|—i—i—i—r 0.2 gamma 0.15 0.04 3.96 fi= 3.75 CT= 0.267= 1.2?: - a= 2.5 0.1 gaussian X= 6.0 gamma—£> 0.02 a= 3.5 0.05 3 4 5 -log(velocity) ( m / s ) 6 ± 3 J 4 I 1_J I I I I 5 -log(velocity) ( m / s ) Figure 3 . 5 a Histograms of velocities sampled i n the uniform aperture, discrete sub-domain model. I L 6 DIRECTION 1 2 5 5 7 3 3 s a m p l e s 0.06 - ' ' ' l _' ' 1 I 1 1 1 1 1 I 1 1 1 I — ' 1 1 DIRECTION 2 0.08 r T—I—I—I—I—i—i—i—i—I—i—i—i—I—I—i—i—i—i—I—r /x= 3.43 gaussian 0.04 a = Q 3 Q 1.2131 l e + 0 6 s a m p l e s gamm 0.06 0.02 h a= 8.2 X= 9.5 . 0.02 a= 2.3 - 3 -log(velocity) ( m / s ) DIRECTION 3 0.06 I I i | I i i i i i_i 4 I i i i 5 i i i 6 -log(velocity) ( m / s ) 236452 samples | 0.30-1 0.04 J—i—I_J—I—i—i—i—i—I—i "i <7= 7= 0.70 \= 8.9 a= 2.1 3.14 gaussian 7= 0.58 «= 12.1 fj,= DIRECTION 4 I—i—i—i—|—I—i—i—i—|—r "i i—i 1737 s a m p l e s |—i—i—i—i—|—i—i—i—i—|—i—i—i—i—|—r- 0.15 M= 3.71" cr= 0.44 0.04 7= 0.65-1 0.1 <x= 9.5 0.02 gamma 0.05 J 3 4 5 -log(velocity) ( m / s ) 6 gaussian i i i I 3 i i i i I 4 i i i i i i i i 5 -log(velocity) ( m / s ) Figure 3.5b Histograms of velocities sampled in the variable aperture, discrete sub-domain model. 6 70 1st Horizontal Moment (Xc) 0 I 0 1 I I I I I I 20 I I I I 40 I L 60 1st Vertical Moment (Zc) 0 T i 1 1 1 i 1 1 1 1 1 1 r J i i I i i i I i i i 1 i 20 40 time (days) : 60 Figure 3.6a Plots of the growth in the first horizontal and vertical spatial moments with time for the uniform aperture simulations. Comparisons are made between the curves from the three continuum models (gamma(O), histogram and gaussian(O)) and the discrete model averaged curves. Comparisons are also made between the curves of the two approximate velocity distribution continuum models (gamma(O) and gaussian(O)) and the histogram velocity distribution model curves. 71 2nd Horizontal Moment (Sxx) 10 discrete gamma(O) 8 — histogram gaussian(O) 6 X vi 0 0 20 40 60 2nd Vertical Moment (Szz) N N VI 40 time (days) Figure 3.6b Plots of the growth in the second horizontal and vertical spatial moments with time for the uniform aperture simulations. Comparisons are made between the curves from the three continuum models (gamma(O),histogram and gaussian(O)) and the discrete model averaged curves. Comparisons are also made between the curves of the two approximate velocity distribution continuum models (gamma(O) and gaussian(O)) and the histogram velocity distribution model curves. i "i i 1st Horizontal Moment (Xc) i i i i i i i i i i l l l r discrete 10 0 — gamma(O) — histogram gaussian(O) J I L J l l 10 0 20 l J L J L 30 1st Vertical Moment (Zc) i 1 r "i 1 1 I L J I I r i 1 1 i i i r 1 0 •1 -2 J 0 I 10 20 time (days) l 30 Figure 3.7a Plots of the growth in the first horizontal and vertical spatial moments with time for the variable aperture (0.2) simulations. Comparisons are made between the curves from the three continuum models (gamma(O), histogram and gaussian(O)) and the discrete model averaged curves. Comparisons are also made between the curves of the two approximate velocity distribution continuum models (gamma(O) and gaussian(O)) and the histogram velocity distribution model curves. 73 0 10 20 30 t i m e (days) Figure 3.7b Plots of the growth in the second horizontal and vertical spatial moments with time for the variable aperture (0.2) simulations. Comparisons are made between the curves from the three continuum models (gamma(O), histogram and gaussian(O)) and the discrete model averaged curves. Comparisons are also made between the curves of the two approximate velocity distribution continuum models (gamma(O) and gaussian(0))and the histogram velocity distribution model curves. Figure 3.8 Schematic diagram of particle motion to illustrate the long and short correlation models. Nodes indicate locations of fracture intersections. 75 lst Horizontal Moment (Xc) i 1 1 r i 1 1 r i 1 1 r ~\ r I J I I.. I J L i r discrete gamma (long) gamma (1) gamma (short|) 10 - gamma(O) u cu -t-i CD 0 1 J L I L I I I 10 0 20 30 1st Vertical Moment (Zc) w 1 a; 0 ~i 1 1 r 1 1 1 r J i i L J I I i_ T i r i i_ 1 0 10 20 time (days) J \ 30 Figure 3.9a Plots of the growth in the first horizontal and vertical spatial moments with time for the variable aperture (0.2) simulations. Comparisons are made between the curves from the two correlated continuum models (gamma(long), and gamma(short)) the discrete model (discrete) and the uncorrelated continuum model (gamma(O)). 76 0 10 20 30 time (days) Figure 3.9b Plots of the growth in the second horizontal and vertical spatial moments -with time for the variable aperture (0.2) simulations. Comparisons are made between the curves from the two correlated continuum models (gamma(long), and gamma(short)) the discrete model (discrete) and the uncorrelated continuum model (gamma(O)). 7? 1st Horizontal Moment (Xc) "i 1 1 "i r 1 1 r I L 1 l 1 1 r ~i r I 1 L J L ~i r discrete gamma(2) gamma(l) gamma(i) 10 gamma(O) co SH -t-> 0> CJ 0 w ti 1 CO 0 CO o i i i J i 0 I 10 "i 1 i r J I I 1_ J 20 30 1st Vertical Moment (Zc) ~i 1 i 1 1 r 1 r I L -1 -2 0 10 i i i i 20 time (days) J 1 30 Figure 3.10a Plots of the growth in the first horizontal and vertical spatial moments with time for the variable aperture (0.2) simulations. Comparisons are made between the curves from the two auto—regressive, long correlation, continuum models (gamma(2), and gamma(l)) and the discrete model, (discrete). The gamma(l) curve is equivalent to the gamma(long) curve in Figure 3.9a,b. The uncorrelated continuum model (gamma(O)) is shown for reference. 78 0 10 20 30 time (days) Figure 3.10b Plots of the growth i n the second horizontal and vertical spatial moments with time for the variable aperture (0.2) simulations. Comparisons are made between the curves from the two auto—regressive, long correlation, continuum models (gamma(2), and gamma(l)) and the discrete model (discrete). The gamma(l) curve is equivalent to the gamma(long) curve i n Figure 3.9a,b. The uncorrelated continuum model (gamma(O)) is shown for reference. 79 1st Horizontal Moment (Xc) i 1 1 r 1 1 1 r -i 1 1 r I L J I I L discrete gamma(2) gamma(O) 10 histogram gaussian(O) Ui %-< CD CD O X 0 J I I L J I 10 0 20 J L J L 30 1st Vertical Moment (Zc) Ui CD CD O N "i 1 1 r J i i L i 1 1 r i 1 r j i i i J I I 1 0 1 0 1 10 20 time (days) L 30 Figure 3.11a Summary plots of the'growth in the first horizontal and vertical spatial moments with time for the variable aperture (0.2) simulations. The progressive improvement of the model, from the gaussian(O) model, (equivalent to the Schwartz and Smith Model [1988]), through the histogram and gamma(O) models to the gamma(2) model, in the match to the discrete curves is illustrated. 80 0 10 20 30 t i m e (days) Figure 3.11b Summary Plots of the growth in the second horizontal and vertical spatial moments with time for the variable aperture (0.2) simulations. The progressive improvement of the model, from the gaussian(O) model, (equivalent to the Schwartz and Smith Model [1988]), through the histogram and gamma(O) models to the gamma(2) model, in the match to the discrete curves is illustrated. DISCRETE MODEL CONTINUUM MODEL Plume Stats. xc= 9.65 0 5 10 meters 15 2 Figure 3.12 Particle distribution plots for the discrete and gamma(2) full-domain models at T=35 days. 1st HORIZONTAL MOMENT (XC) n 1st VERTICAL MOMENT (ZC) 1 — i — | — i — i — i — i — | — i — i — i — i — | — j 1 I I I | I I I I | I I /x= 9.82 10 - I I | I 1 I T" /j,= -0.08- <i= 1.88 _ 10 - c= 0.60 7= -0.21 7= Jl 0.05 " 11 5 - Jl 0 ' J ' i * ' i ' 0 i 10 i i i i i i i I i i i i 15 i i i i JL i i I i i i i i i i I -1 2nd HORIZONTAL MOMENT (SXX) T i i i i _ i I i i i i i i i i i 0 2nd VERTICAL MOMENT (SZZ) r jx= 8.90 8 - a= 2.97 - 10 - 7= 1.78 - 6 - 4 - 5 - \ - 2 u 0 I I I I I I I I I 10 I I I I I 15 I I I I I 20 I I I I 0 25 Figure 3.13 Histograms of the first and second spatial moments of the discrete, variable aperture (0.2) full—domain plume compiled from 100 Monte Carlo simulations at T=35 days. 0 20 40 60 80 100 Time (days) Figure 3.14 Breakthrough curves for the discrete and continuum, full—domain variable aperture simulations. Comparisons of the discrete gamma(2) and gaussian(O) models to the average of 10 discrete simulations. 84 4.0 TESTING THE EXTENDED 2-D MODEL Up to this point, simulations of solute transport by the stochastic-continuum method have been restricted to varying a limited set of input parameters. In this chapter three of these input parameters, (gradient orientation, aperture distribution and fracture length distribution), are varied to observe how well the improved stochastic-continuum model, gamma(2), reproduces plumes created in full-domain discrete simulations under a range of network conditions. The ability of the stochastic-continuum model to simulate the mass distribution observed in the discrete network simulations will be analyzed by comparing the development of the first and second spatial moments of the plumes, with respect to the XZ coordinate system, from full-domain discrete and continuum models. As before the discrete moment curves will be the average of 100 Monte Carlo simulations. Particle plots of the plumes are also given for some simulations to show the nature of the mass distribution. Statistics for the continuum model are collected from 200 realizations of the appropriate discrete sub-domain model. The travel time through the sub-domain network is adjusted to confine statistical sampling to a region within the domain unaffected by boundary effects, as discussed earlier. Table 4.1 lists the simulations carried out in this chapter, giving reference to the corresponding plots. If not specifically stated, the input parameters for the simulations are equivalent to those given in Table 3.1, for the variable aperture (0.2) simulations tested under a 0° hydraulic gradient orientation. 4.1 Gradient Orientation Because the two fracture sets described in Table 3.1 are statistically identical, varying the gradient orientation over a range of 45° will represent all possible orientations of the flow field with respect to the network. Four orientations, (0 ,15°,30 ,45°), over this range were chosen for analysis. The 0° orientation o o 85 simulations have already been described in Chapter 3. Plots of the temporal development of the spatial moments of the plumes for the continuum, (gamma (2) and gaussian(O)), and discrete models at a gradient orientation of 0° were shown in Figures 3.11a and b. Far-field percentage differences in the spatial moments, summarized in Table 3.5, show a difference in of 6% and a difference in S x x +S z z of 17%. The far-field percentage difference, as developed in Chapter 3, is defined as (1-A /A ) c d where A and A^ are the slopes of the continuum and discrete moment curves, c respectively. Similar curves for simulations at 15°,30° and 45° are shown in Figures 4.1a,b, 4.2a,b and 4.3a,b, respectively. Corresponding particle plots for all four angles are given in Figures 3.13, 4.4, 4.5 and 4.6. The network used is the same for all four discrete simulations. The results at 45° are limited to a maximum time of 25 days as a large proportion of the particles exit the domain after this time. Even at 25 days a significant amount of mass has exited the domain for the discrete simulation at 45° (Figure 4.6). There is not as great a loss in the continuum simulation as the particles at the front of the continuum plume do not travel as far. This fact is seen in the underestimation of the second moments by the continuum model in Figure 4.3b. The effect of the loss of particles on the calculated moments for the discrete model can be seen in Figures 4.3a and b. The loss of particles affects the vertical moments more than the horizontal moments, as the particles lost will have travelled the greatest distance vertically but are approximately evenly distributed in terms of their horizontal position. In Figure 4.3b the two sets of curves ( S a n x x d S ) should be symmetric due to the symmetry of the z z network at 45°. At 25 days the values of S x x and S z z for the gamma(2) model are approximately equal (2.33 and 2.32, respectively). However, the values of S x x and S z z for the discrete model at 25 days are noticeably different (2.77 and 2.32, respectively). The reduction in the value of S z z as compared to the more correct value of S x x is due to the loss of particles out of the bottom of the domain. The loss of particles is also 86 noted in the curves of the first moments. Values of (4.91) and ZQ (-4.87) are approximately symmetric for the gamma(2) model at 25 days. For the discrete model the matching values are 4.96 for X ^ and -4.63 for ZQ. Comparing these first moment results to the second moment values, it can be seen that the effect of the loss of particles is greater for the second moments than the first. A similar loss of particles is noted in the particle distribution plots for the 30° discrete network simulation (Figure 4.5) at 35 days. The effect of this loss shows up in a decrease in the slope of the second vertical moment, discrete curve (S^) after 30 days (Figure 4.2b). This loss of particles also affects the values of the other moments to a lesser degree. For example, in Figure 4.2a and b, a slight increase in the slopes of both ZQ and SxX» is noted after 30 days. The increase in the slope of S X x is due to the fact that the particles exiting have average horizontal positions. Removing these particles increases the weight of the particles that have extreme horizontal positions in the calculation of S^x- The increase in the negative slope of ZQ is caused by the removal of the particles which had travelled the greatest distance vertically downwards. Any one individual discrete plume distribution (Figures 3.13, 4.4, 4.5, and 4.6) will not represent the average discrete plume distribution. Thus comparison of the discrete and continuum models are best made by comparing the temporal development of the first and second spatial moment curves, or the far-field percentage differences calculated from these curves. The discrete curves in these plots are the average values from 100 realizations of the discrete model. Table 4.2 summarizes the far-field percentage differences between the spatial moments of the discrete and continuum simulations. Due to the loss of particles in the 45° discrete simulation, values for AZ^-. and AS^z are not reported. However, the values of AZ^ and ASz^, due 1 0 t n e symmetry at 45°, should be equivalent to the values of AX^ and ASXx> respectively. Values for AX<^ and A S x x a r e calculated from the linear portion of the curves prior to 87 20 days, to avoid the effect of the loss in particles. The results, however, are questionable as the mean distance travelled at 20 days (5.6 m) is not much greater than the distance travelled in the discrete sub-domain to collect the statistics (approximately 4.3 m). Also there will be some loss of particles out of the bottom of the domain even at 20 days. For the same reasons as in the 4 5 ° case, values for A Z ^ and A S 3 0 ° simulations are not listed in Table 4.2. Values of A X ^ and A S X X Z Z for the are calculated from the linear portion of the curves prior to 30 days. The mean travel distance of the plumes at 30 days is 8.1 meters for the 3 0 ° simulations. For all orientations both the gamma(2) and gaussian(O) first moment curves match the discrete first moment curves quite well. Far-field percentage differences range from 0% to 9% for the gamma(2) model, as listed in Table 4.2. The difference in the second horizontal moments in general increase for the 1 5 ° , 3 0 ° and 4 5 ° simulations (28%, 29% and 22%, respectively) relative to the difference for the 0 ° simulation of 16% between the gamma(2) and discrete models. The far-field percentage differences between the gaussian(O) and discrete horizontal second moment curves also increase as the hydraulic gradient orientation moves away from the horizontal. The value of 22% for the far-field percentage difference in the S x x for the 4 5 ° simulation is somewhat questionable, as mentioned above, but may show that there is an improvement in the match of the model at symmetric orientations ( 4 5 ° ) compared to the non-symmetric, off-fracture orientations ( 1 5 ° and 3 0 ° ) . match in S z z Due to the lack of information, little can be said about the as the hydraulic gradient changes (however, see discussion on symmetry information below). The two values of far-field percentage differences for the 0 ° and 15° simulations are approximately equal (29% and 30%, respectively). value of 22% for A S Z Z at 4 5 ° (see above for discussion of the validity of this value) may indicate that the far-field percentage differences in S orientation of the A symmetric z z will decrease as the hydraulic gradient shifts away from the horizontal. 88 By symmetry of the network we can assume that the far-field percentage differences of A S X X at 75° and 90° will be equivalent to the far-field percentage differences of A S ^ z at 15° and 0°, respectively. Thus the values of far-field percentage differences of A S X X for 0°, 15°, 30°, 45°, 60°, 75° and 90° will be 16%, 28%, 29%, 22%, ?, 30% and 29%. The value at 60° is unknown but would be expected to lie between 22% and 30%. Values for A S ^ z for the same series of angles would be in the reverse order. From the simulations tested on this network the maximum farfield percentage difference is approximately 30%. The quality of the predictions of the spatial moments by the stochastic-continuum model is still not perfect; 30% may be considered too large of an error for some situations. However, compared to other continuum models for transport in fractured media these results could be considered acceptable. As discussed in the previous chapter there are several possible source of the differences between the discrete and stochastic-continuum mass distributions. Crosscorrelation between velocities of different directions is one possible source of the differences as is described in the next section. Table 4.3 lists the first and second correlation values (r(l) and r(2)) for the four orientations of the hydraulic gradient. As would be expected the correlations in direction 2, (horizontally right), decrease while the correlations in direction 3, (vertically down) increase as the gradient angle approaches 45°. Cross-correlations between successive fracture segment velocities when the direction changes from 2 to 3 or viceversa are expected to also increase. The degradation in the match of the second moment curves as the hydraulic gradient moves away from the horizontal could be partially due to the lack of cross-correlation in the continuum model. 4.2 Variation in Aperture Distribution Two values, (0.0 and 0.2), for the standard deviation in logarithmic aperture have already been tested for a hydraulic gradient orientation of 0°. The temporal 89 development of the spatial moments for the uniform case, (0.0), are shown in Figures 4.7a,b, which include curves from the gamma(2) model, which were not shown earlier. The far-field percentage difference for S x x between the gamma(2) and discrete model is 9%. Similar curves for a standard deviation in aperture of 0.2 have already been illustrated in Figures 3.11a,b. The far-field percentage difference for S x x between the gamma(2) and discrete model, for the 0.2. aperture variability case is 16%. An intermediate value of the standard deviation in log aperture (0.1), was also tested. For this case, 95% of the fracture apertures created will fall between 20 and 50 u, the mean being 32 u. The temporal development of the spatial moments for this case can be seen in Figures 4.8a,b. The far-field percentage difference S x x between the gamma(2) and discrete model for the 0.1 aperture variability is 15%. The far-field percentage differences in the first and second spatial moments between the discrete and gamma(2) models for all three simulations are listed in Table 4.4. Values for the first and second correlation coefficients for the three simulations are listed in Table 4.5. The rate of rise in the first correlation coefficients in direction 2, for the three simulations in order of increasing standard deviation in aperture, (0.51, 0.64, 0.69), approximately matches the rate of rise in the percentage far-field difference in A S X X , (9%, 15%, 16%). This similarity suggests that the difference in the second moments, A S X X and AS2z» between the gamma(2) and discrete simulation curves is related to a correlation effect, possibly the missing cross-correlation in velocities of different directions as discussed earlier. The close match of the far-field percentage differences in the second moments for the 0.1 and 0.2 variable aperture simulations show that the match of the gamma(2) model mass distributions to discrete model distributions will be consistent for small values of standard deviation in log fracture aperture (0.0 to 0.2). 90 4.3 Variation in Fracture Length The ability to test the model for various fracture networks is somewhat limited as the discrete model is constrained to producing two orthogonal fracture sets. However, the length of the fractures distributions for the two sets can be varied. Two variations of the fracture lengths were tested, labelled Network B and Network C. Network A is the original network as described in Table 3.1. For Network B, the mean length of the fractures in set 1 (parallel to the direction of flow) was increased to 0.6 m, while the mean fracture length of the perpendicular set was decreased to 0.2 m. Both sets of fractures in Network A have a mean length of 0.4 m. For Network C, the mean lengths for the two sets are reversed, the set parallel to the direction of flow being shorter on average than the perpendicular set. The uncleaned and cleaned fracture networks for one realization of Network B are illustrated in Figures 4.9a and b. The horizontal elongation of the network in the horizontal plane can be seen when compared to fracture Network A shown in Figures 2.1. Corresponding plots for Network C can be observed by rotating the plots in Figures 4.9 by 90°. Fracture histograms for Network B, Figure 4.10, also show the elongation in the horizontal direction. Equivalent histograms for Network C can be observed by switching the two histogram sets shown in Figure 4.10. Figures 4.11a,b and 4.12a,b show the development of the first and second spatial moments of the discrete and continuum simulation plumes for Networks B and C, respectively. Corresponding curves for Network A are found in Figures 3.11a,b. The gamma(2) model matches the first moments quite well for Network B, the far-field percentage difference in is 8%, close to the value of 6% for Network A. Figure 4.11b shows the first simulation where the continuum model overestimates the first horizontal moment for the entire length of the simulation. The actual far-field percentage difference is the smallest of all the simulations run so far at 4% and the slope of the curves indicate that the continuum model will eventually underestimate the first 91 horizontal moment of the discrete simulation. The large increase in the percentage differences for AS 2 Z m t n e far-field Network B simulation (46% compared to 29% for Network A ) is misleading as the total value of S ^ z for the discrete simulation is quite small (0.55 m at 20 days). The visibly large difference between the gamma(2) and discrete S ^ x curves in Figure 4.11b is developed in the near source region. In the near-source region the discrete model has fewer choices of flow paths to move the solute through than in the far-field region due to the discrete nature of the network and the fact that source is a point. The continuum model, however, has in effect the same number of flow paths throughout the domain. This effect causes an overestimation of 5>XX by the continuum model in the near-source region which can be observed in all of the previous simulations but is not so pronounced as in the Network B simulation. The reason for the large near-source overestimation of S x x m Network B is the decrease in the number, of possible flow paths near the source in Network B caused by the horizontal elongation of the network. The small far-field percentage difference in ^XX ^**ZZ + reflects the general improvement in match of the second moments relative to the far-field percentage difference of 17% for &$XX ^ZZ + Network A simulations. from the Table 4.6 summarizes the far-field percentage differences between the gamma(2) and discrete first and second spatial moment curves. Figure 4.13 and 4.14 show the plume distributions for simulations with Network B and C , respectively. The elongated nature of the discrete network in the horizontal direction produces a very elongated plume as seen in Figure 4.13. Individual discrete simulations are strongly controlled by the available flow paths near the source. The non-linear portion of the S ^ x curve (Figure 4.11b) extends to 10 days at which time the mean position of the plume is slightly more than four meters from the source (Figure 4.11a). Comparisons to Network A show the same non-linear region ending at approximately 10 days at which time the mean position of the plume is less than three meters from the source (Figures 3.11a and b). These results along with the discussion in 92 the previous paragraph, suggest that Network B has a larger near-source region and solutes travel further through the network before the temporal development of the second moments become linear. For the Network C simulations the far-field percentage difference for AX^> shows a large increase (21%) while the far-field percentage differences for A S AS 2 decrease to 0% and 21% respectively. The value of A S Z X X +AS Z Z X X and also shows a decrease to 4% compared to 17% for Network A. The results, however, are inconclusive as the large amount of lateral spreading reduces the distance over which the plume can travel before losing particles out of the top and bottom of the modeling domain. The loss of particles can be seen in the change in slope of S z z after 40 days. The values of far-field percentage difference mentioned above were calculated in the linear portion of the moment curves before 40 days. At 40 days the mean position of the average discrete plume is 3.7 meters which is approximately the same distance over which statistics were collect for this network. The nature of the plume distribution in Figure 4.14 reflects the vertically elongated nature of the discrete network. The flow paths require increased travel perpendicular to the direction of the mean hydraulic gradient, creating the larger spreading of the plume vertically. Further testing in a larger domain is necessary to determine if the far-field percentage differences calculated for the simulation are truly representative of the errors in the stochastic-continuum estimate of the moments to be expected with the gamma(2) model in Network C. The results from Network B and C show that as the length of fractures in the direction of flow increase relative to the length of cross fractures the match in the spatial moments of the solute distribution between the gamma(2) and discrete models generally improve. 93 4.4 Conclusion In conclusion the gamma(2) model is able to reproduce discrete network plumes for a variety of network conditions. Percentage errors in the first moments are less than 10%, while percentage errors in the second moments range from 6% to 30%. The addition of cross-correlation to the gamma(2) model may decrease the error, especially for conditions where the hydraulic gradient is not parallel to one of the fracture set orientations. Further testing of the model under different network conditions will be left for future research to be carried out with a three-dimensional stochastic-continuum model. The design and initial results of a three-dimensional stochastic-continuum model are discussed in the next chapter. Table 4.1 List of simulations carried out during testing of the stochastic-continuum model with various input parameters. Variation in Gradient Orientation Gradient Angle Figure Numbers 1st and 2nd Spatial Moments 0 15 30 45 Particle Plot 3.13 4.4 4.5 4.6 3.11 a,b 4.1 a,b 4.2 a,b 4.3 a,b Variation in Aperture Distribution Std. Dev. in log Aperture (m) Figure Numbers 1st and 2nd Spatial Moments 0.0 0.1 0.2 Particle Plot 4.7 a,b 4.8 a,b 3.11 a,b 3.13 Variation in Fracture Length Distribution Network A B C Set 1 Set 2 Figure Numbers mean (m) cutoff (m) mean (m) cutoff (m) 1st and 2nd Spatial Moments Particle Plot 0.4 0.6 0.2 1.2 1.8 0.6 0.4 0.2 0.6 1.2 0.6 1.8 3.11 a,b 4.11 a,b 4.12 a,b 3.13 4.13 4.14 95 Table 4.2 Summary of the far-field percentage differences in the first and second spatial moments between the discrete and gamma(2) models for 4 orientations of the hydraulic gradient to the network. Gradient Angle AX (%) C AZ C AX^+AZ C AS X X AS +AS AS VV (%) 16 29 17 15^ 28 30 28 30 c 29 45 c Loss of particles may affect results. 22 zz 96 Table 4.3 Summary of the first two correlation coefficients for 4 orientations of the hydraulic gradient to the variable fracture aperture, discrete sub-domain simuations. 0° Dir 1 Dir 2 Dir 3 Dir 4 r(D 0.41 0.69 0.40 0.00 r(2) 0.22 0.51 0.21 0.00 15° Dir 1 Dir 2 Dir 3 Dir 4 r(D 0.38 0.68 0.45 0.31 r(2) 0.20 0.51 0.25 0.16 30° Dir 1 Dir 2 Dir 3 Dir 4 r(D 0.37 0.64 0.51 0.44 r(2) 0.19 0.45 0.31 0.00 45° Dir 1 Dir 2 Dir 3 Dir 4 r(D 0.36 0.57 0.57 0.37 r(2) 0.19 0.37 0.37 0.21 97 Table 4.4 Summary of the far-field percentage differences in the first and second spatial moments between the discrete and gamma(2) models for three values of aperture standard deviation. Standard Deviation in Log Aperture AX (%) C AS^y (%y AS (%7 Z Z AS +AS (%) X X 0.0 5 9 38 13 0.1 4 15 43 17 0.2 6 16 29 17 Z Z 98 Table 4.5 Summary of the first two correlation coefficients for three values of the standard deviation in fracture aperture. 0.0 Dir 1 Dir 2 Dir 3 Dir 4 r(D 0.10 0.51 0.12 0.00 r(2) 0.02 0.31 0.01 0.00 0.1 Dir 1 Dir 2 Dir 3 Dir 4 r(D 0.22 0.64 0.19 0.00 r(2) 0.08 0.46 0.06 0.00 0.2 Dir 1 Dir 2 Dir 3 Dir 4 r(D 0.41 0.69 0.40 0.00 r(2) 0.22 0.51 0.21 0.00 Table 4.6 Summary of the far-field percentage differences in the first and second spatial moments between the discrete and gamma(2) models for three statistically different fracture networks. Network AX/(fe AS v A S ^ V cif % AS Y+Aa 7 7 Y x x (%) zz Network A 6 16 29 17 Network B 8 4 46 6 Network C 21 0 12 4 100 1st Horizontal Moment (Xc) i i i i r 1 1 r i 1 1 r 1 r J I I L J L discrete gamma(2) gaussian(O) 10 0 1 0 1 1 L I 10 I I I 20 30 1st Vertical Moment (Zc) 10 20 time (days) Figure 4.1a Plots of the growth in the first horizontal and vertical spatial moments with time for a gradient orientation of 15 degrees. 10 2nd Horizontal M o m e n t (Sxx) i 1 1 1 1 1 1 1 1 1 1 1 1 [— ' ' I discrete gamma(2) 8 gaussian(O) 6 - 2 - x x Vi 0 0 10 2 1.5 1 =vi 20 L J 30 2nd Vertical M o m e n t (Szz) o 0 10 20 time 30 (days) Figure 4.1b Plots of the growth in the second horizontal and vertical spatial moments with time for a gradient orientation of 15 degrees. l_ 103 ~i 1 10 1 1st Horizontal Moment (Xc) 1 1 1 1 1 1 1 1 1 r r — discrete -- gamma(2) gaussian(O) 8 co CD CD 6 - a, o X 2 0 J u L 10 0 l I l l 20 30 1st Vertical Moment (Zc) 0 CO I "i r i r J L -2 CD CD a -4 o N -6 j 0 i i L 10 i L 20 time (days) j i i L 30 Figure 4.2a Plots of the growth in the first horizontal and vertical spatial moments with time for a gradient orientation of 30 degrees. 104 2nd Horizontal Moment (Sxx) 8 "i 1 1 r i 1 r 1 1 r "i r l I L J L discrete gamma(2) gaussian(O) 6 X X in 0 J I I 10 0 i L 20 30 2nd Vertical Moment (Szz) i r 1 1 1 j i i. - -1 • • 1 1 ' - - 1 1 1 r i 1 1 r i r L J I I L J L 3 2 N N Vi 1 0 0 10 J 20 time (days) 30 F i g u r e 4.2b Plots of the growth i n the s e c o n d h o r i z o n t a l a n d v e r t i c a l s p a t i a l m o m e n t s with t i m e for a gradient o r i e n t a t i o n of 30 degrees. 105 T 0 1 1st Horizontal Moment (Xc) 1 1 1 1 i 1 1 1 1 1 10 20 1 1 1 30 time (days) Figure 4.3a Plots of the growth in the first horizontal and vertical spatial moments with time for a gradient orientation of 45 degrees. r 106 2nd Horizontal Moment (Sxx) 6 "i 1 1 2 0 Vi J I L 0 1 1 r l l L 10 "i 1 1 r "i r J L 20 n r J L 30 2 n d V e r t i c a l Moment (Szz) "i 4 - 2 - 0 1 - 6 N N 1 discrete - gamma(2) - gaussian(O) — X X vi 1 1 1 1 1 1 J 0 1 1 r i I L 10 ~i 20 time 1 1 r 30 (days) F i g u r e 4.3b Plots of the growth i n the s e c o n d h o r i z o n t a l a n d v e r t i c a l s p a t i a l m o m e n t s with t i m e for a g r a d i e n t o r i e n t a t i o n of 45 degrees. CONTINUUM MODEL Plume Stats. 9.30 -2.56 6.82 1.42 35 days ^5 0 15" 10 meters 15 20 Figure 4.4 Particle distribution plots for the discrete and continuum full—domain models with a gradient orientation of 15 degrees. 108 DISCRETE MODEL CONTINUUM 0 MODEL 5 10 meters 15 20 Figure 4.5 Particle distribution plots for the discrete and continuum full-domain models with a gradient orientation of 30 degrees. DISCRETE MODEL Figure 4.6 Particle distribution plots for the discrete and continuum full—domain models with a gradient orientation of 45 degrees. 1st H o r i z o n t a l 0 M o m e n t (Xc) 20 40 60 1st V e r t i c a l M o m e n t ( Z c ) co CD CD O 1 0 -1 J N 0 I I I I I 20 I I 40 time I I I I 60 (days) Figure 4.7a Plots of the growth in the first horizontal and vertical spatial moments with time for the gamma(2) and discrete models uniform aperture simulations. The gamma(2) and gaussian(O) curves are shown for reference. L 2nd Horizontal Moment (Sxx) 10 i i 1 n r r "i 1 J 1_ r discrete gamma(2) 8 gamma(O) gaussian(O) 6 CM X X 0 s 3 J J 0 L 20 2 E- i - 0-5 0 L 40 60 2nd Vertical Moment (Szz) T r I0 20 i r J L 1 r 40 time (days) 60 Figure 4.7b Plots of the growth in the second horizontal and vertical spatial moments with time for the gamma(2) and discrete models uniform aperture simulations. The gamma(2) and gaussian(O) curves are shown for reference. 112 1st H o r i z o n t a l M o m e n t "i 1 r i i (Xc) r i i r discrete gamma(2) 10 0 0 20 40 1st V e r t i c a l M o m e n t 20 F i g u r e 4.8a (Zc) 40 time (days) P l o t s of t h e g r o w t h i n t h e f i r s t h o r i z o n t a l a n d v e r t i c a l s p a t i a l m o m e n t s with time for the gamma(2) and discrete models variable a p e r t u r e (0.1) s i m u l a t i o n s . 60 60 113 2nd Horizontal Moment 10 "i 1 i r (Sxx) r i i r discrete gamma(2) 8 6 0 0 20 40 2nd Vertical Moment 20 (Szz) 40 time (days) Figure 4.8b Plots of the growth in the second horizontal and vertical spatial moments with time for the gamma(2) and discrete models variable aperture (0.1) simulations. 60 60 a) uncleaned network 0 5 meters b) cleaned network Figure 4 . 9 One realization of Network B 10 HORIZONTAL SET 0.25 T ~i 1 r T ~i 1 r 1 1 r 0.2 o 0.15 0.2 0.4 0.6 0.8 fracture length (m) VERTICAL SET 0.25 i 1 r 1 — I i i 1 1 1 1 1 1 1 1 1 0.2 cd o 0.15 CD OZ) CU 0.1 o 0.05 0 _i 0 0.2 i_ J 0.4 0.6 0.8 fracture length (m) Figure 4.10 Fracture length histograms for Network B I L_ 116 1st H o r i z o n t a l M o m e n t "i i i i I i i i i 1 1 (Xc) i i i 1 I L I i i i r J I I L discrete gamma(2) 10 - 5 - co CD CD X J 0 I I L 0 J I I L 5 J 10 15 1st V e r t i c a l M o m e n t co f~t CD 1 CD 0 E- -i-> a o N •1 "i 1 1 J I L r "i r i 20 (Zc) i i r i i i i 1 r i i '- =- 0 J i i i 5 10 time i 15 (days) Figure 4.11a Plots of the growth i n the first horizontal and vertical spatial moments -with time for the gamma(2) and discrete models Network B simulations. i 20 117 2nd Horizontal Moment 10 i i r T i i i r (Sxx) i i i r J I I L ~i 1 1 r J I I L discrete gamma(2) 8 6 CNJ X X vi 0 J 1 I L 0 5 10 15 2nd Vertical Moment 2 1.5 i 1 1 r n i i r 20 (Szz) ~i 1 1 r J i i i T i i r i i i 1 N Vi 0.5 0 0 5 10 time i 15 (days) Figure 4.11b Plots of the growth in the second horizontal and vertical spatial moments with time for the gamma(2) and discrete models Network B simulations. - 20 118 1st H o r i z o n t a l "i 1 1 1 i - Moment i (Xc) r l i r discrete gamma(2) 10 CO U CD CD X 1 0 0 20 40 1st V e r t i c a l M o m e n t 60 (Zc) CD O N 20 40 time (days) Figure 4.12a Plots of the growth in the first horizontal and vertical spatial moments with time f o r the gamma(2) and discrete models Network C simulations. 60 119 2nd H o r i z o n t a l M o m e n t (Sxx) 8 -i 1 r- 1 "i 1 r l 1 r discrete gamma(2) [- 2 0 J 0 L J 20 L 40 60 2nd V e r t i c a l M o m e n t (Szz) n l 0 r i T r i r J L FJ 0 L j L 20 40 time (days) Figure 4.12b Plots of the growth in the second horizontal and vertical spatial moments with time for the gamma(2) and discrete models Network C simulations. 60 DISCRETE MODEL Plume Stats. xc= 11.93 zc= 0.50 sxx= 11.05 szz= 0.65 25 days —> 0° CONTINUUM MODEL Plume Stats. xc= 10.27 zc= 0.01 sxx= 10.91 szz= 0.35 0 5 10 meters 15 20 Figure 4.13 Particle distribution plots for the discrete and continuum full—domain simulations for Network B. DISCRETE MODEL Figure 4.14 Particle distribution plots for the discrete and continuum full—domain simulations for Network C. 122 5.0 THREE-DIMENSIONAL MODEL Although modeling fracture networks in two-dimensions is a useful step towards understanding the mechanisms involved in the transport of solutes, to properly simulate transport in a fracture network the model must be able to recreate the three-dimensional interconnectivity of the network. The original plan of this thesis was to quickly test the two-dimensional stochastic-continuum model for various network conditions and then move on to developing and testing a three-dimensional model. A fair amount of time was spent initially developing the three-dimensional model, however, after only a few runs were carried out in three-dimensions the focus of the research shifted to developing the gamma(2) two-dimensional stochastic-continuum model. The relatively poor match in the second moments with the gaussian(O) two-dimensional model for simulations with variable fracture aperture, as shown in Figure 3.4b, had not been expected. Because of this change in focus, only the outline of a three-dimensional analysis was carried out. This chapter documents the initial research carried out on the three-dimensional stochastic-continuum model. As in the two-dimensional method, the model is made up of two parts: the discrete model in which the particle motion statistics are collected, and the continuum model in which the particles are moved by rules based on the statistics collected. 5.1 The Discrete Model A three-dimensional discrete network transport code, similar to the one used by Smith et al. [1985], was modified to allow for collection of particle motion statistics. The following sub-sections detail the significant features of the discrete model. 123 5.1.1 Aperture Due to the irregularity of fracture surfaces, the value of aperture (distance between the bounding fracture surfaces) varies within the plane of the fracture. As described in Chapter 1 the values of aperture are also spatially correlated, in that adjacent apertures are close in value. In three-dimensions, fractures are approximated by rectangular planes. Each plane is in turn represented by a finite element mesh made up of triangular elements. Aperture values are specified at the node points, (vertices of the triangular elements) and are assumed to vary linearly across the triangular elements. A lognormal distribution of aperture within the fracture plane is chosen. To effect the correlated nature of the apertures, a two-dimensional auto-correlation function is used with the generation technique based on the matrix decompostion approach (Clifton and Neuman, 1982). The aperture at each node is conditioned on the values of aperture at neighbouring nodes. An input correlation length determines over what distance the effect of the correlation extends. Variation in fracture aperture within the plane of a single fracture has been studied in detail by Tsang and Tsang [1987]. 5.1.2 Geometry Fractures can vary in orientation with respect to the X and Z axes but must be parallel to the Y axis, resulting in a constant strike for all fractures in the Y direction. Within a given set, individual fracture orientations (dip angles) are sampled from a normal distribution given a mean and standard deviation in the angle with respect to the X axis. This allows for the fractures in a set to vary in orientation about the mean orientation of the set. Fracture length and widths are sampled from the tail of a negative exponential distribution, given a minimum cutoff length. Unlike the twodimensional model, where dead-end fracture segments were eliminated, fracture segments which only connect to one other fracture are retained as flow may enter and exit across different portions of the same fracture intersection. Due to the increase in 124 the required number of nodes per fracture, (an increase of approximately 2 orders of magnitude), the number of fractures that can be modeled with available resources is greatly reduced. For the initial testing of the model two fracture sets of 200 fractures each were chosen. Fracture set one has a mean dip of 60°, while fracture set two has a mean dip of 0°. The standard deviation in dip angle from both fractures is 0.5°. The input parameters for the discrete simulation presented here are listed in Table 5.1. The density of the fracture network, calculated using the method of Rouleau and Gale [1985], is 0.43 m . This density is approximately one third of the density used in the -1 two-dimensional simulations. An example of the fracture network used here is shown in a two-dimensional XZ slice through the network at Y=5.0 m in Figure 5.1. Only those fractures which intersect the Y=5.0 m plane are shown. 5.1.3 Boundary Conditions and Gradient Orientations The model will allow for either constant head or constant flux boundary conditions for each of the six sides of a box shaped domain. The value of constant head or constant flux must be constant for the entire side of the domain. Because of this limitation uniform gradients across the domain can only be created parallel to the coordinate axes. This is not a major drawback as, unlike the two-dimensional model, the orientation of the fracture network can be changed to simulate various orientations of the hydraulic gradient with respect to the network. For the simulations presented here only one orientation of the uniform hydraulic gradient, parallel to the X axis, was tested. To create this uniform gradient, appropriate constant head values are assigned to all nodes on the X-minimum and X-maximum boundaries. Nodes on the Y and Z minimum and maximum boundaries were set to a constant flux of zero. 125 5.1.4 Computational Method Fractures in three-dimensions are subdivided into rectangular quadrants, similar to the fracture segments in the two-dimensional model. Fracture plane 1 in Figure 5.2 is divided into three quadrants labelled A,B and C. Partial intersections such as that between fracture planes 1 and 2 (boundary 4b) are extended across both fractures to delineate the quadrants (boundary 4a and 4c). A record is kept of where the fracture intersection begins and ends along partial intersecting boundaries so that the code will know which quadrants boundary nodes lie along the actual intersection. The Galerkin method (Zienkiwicz and Morgan, 1982) for linear triangular elements is used to create the finite element system of equations for each quadrant. The resulting sub-matrices for each quadrant are then combined into the global matrix by matching and re-numbering the duplicate nodes at fracture intersections. Once the global matrix has been assembled the hydraulic head at the nodes is calculated using an iterative conjugate gradient method (Press, 1988). 5.1.5 Particle Motion in the Discrete Network Transport of the particles in three-dimensional networks is more complicated than in the two-dimensional model as velocities will vary within the plane of the fracture. Particles are moved in small time steps across the quadrants at the rate determined by the local fluid velocity. Velocities within the triangular elements are calculated from the hydraulic head values at the nodes. The time step is set so that a particle will move a maximum of one tenth of the smallest dimension of the element before the local velocity is recalculated. Thus variations in the local transport velocity over short distances are correctly simulated. When a particle reaches the boundary of a quadrant one of four possible quadrant boundary conditions may occur, as is described with the aid of Figure 5.2. At a zero-flux boundary (boundaries 1 and 3) the particle is reflected back into the quadrant. Quadrant boundaries which intersect the domain 126 boundaries (boundary 5) may be zero-flux or constant-head boundaries. For constant head boundaries the particle is removed from the domain. At fracture intersections boundaries (boundaries 2 and 4b) complete mixing is assumed. A frequency table, based on the percentage flow into each quadrant is sampled to determine which fracture quadrant the particle will move into. At imaginary boundaries such as boundary 4a, particles are transferred directly from quadrant to the next within the same fracture. 5.1.6 Code Verification The code utilized in this thesis was a modification of a code originally developed by Smith et al. [1985] (C. Mase, per. comm.). To check for errors in the code simple networks for which the flow volumes in individual fractures could be analytically determined were tested. 5.1.7 Movement Statistics One of the tasks in developing a three-dimensional stochastic-continuum model is to determine the best set of statistics to be collected to represent particle motion in the discrete network. Due to the variation in the local velocity within the fracture quadrants the travel path of the particles will not be linear. For each fracture quadrant a local XY coordinate system is used to measure the distance along, (dx), and across, (dy), the quadrant between the entrance and exit locations of the particle (Figure 5.3). These two 2 2 parameters, (dx and dy), along with the average velocity (time/\/(dx +dy )) define the movement of particles across the fracture quadrant, similar to the length of the fracture segment and velocity recorded for each fracture travel segment in the two-dimensional model. As the strike of the fractures must be parallel to the Y axis, all fracture intersections occur along the local X-minimum and X-maximum quadrant boundaries. Thus for a given fracture set particles must enter from one of these local X boundaries. 127 As noted before particles may enter and exit across the same boundary, though typically particles will move from one boundary to the other. The statistics for particle motion can be categorized into two entrance directions for each fracture set, or a total of four entrance directions for the two sets used in the simulations presented here. As in the two-dimensional model the choice of direction is recorded as particles cross intersections. The choice of direction for the case where particles are transferred directly across an imaginary boundary (eg. boundary 4a, Figure 2.5) are recorded as if the boundary was an ordinary fracture intersection. For each of the four directions, histograms for dx and dy and approximate gaussian velocity distribution are created. The assumption of a gaussian distribution for velocities was made before the revisions to the two-dimensional model showed that velocity data was skewed for the two-dimensional simulations. The importance of high correlations in velocities for the two-dimensional model had also not been identified at the time of this work on the three-dimensional model, thus velocity correlations were not included in the statistics collected. Future work should investigate the validity of the gaussian assumption for the distribution of velocities and determine the effect of velocity correlations on particle motion. 5.2 The Continuum Model The design of the continuum model for the initial testing of the threedimensional model is somewhat simplified compared to its' two-dimensional counterpart. Determining the distribution of hydraulic head in the continuum is unnecessary as the magnitude and direction of the hydraulic gradient are constant. Particles are transported through the domain by rules based on the statistics collected in the discrete model and the magnitude and direction of the local hydraulic gradient. The rules are very similar to those used in the two-dimensional model. Particles are moved in steps across imaginary fracture quadrants, until the time period desired is completed. At each 128 imaginary intersection the direction, velocity and local travel lengths (dx and dy) for the next step are sampled from the statistical distributions collected in the discrete model. As the orientation of the fractures vary within a set the orientation of the imaginary quadrant is sampled from the same distribution used to create the fractures in the discrete model. 5.3 Comparison of the Discrete and Continuum Models Due to time constraints, testing of the stochastic-continuum model was restricted to comparisons of one realization of the solute distribution from the discrete and continuum models. The statistics used in the continuum model come directly from the single simulation of the discrete model. Although the testing of the three-dimensional model is not as rigorous as that done for the two-dimensional model some information about the choice of statistics collected can be obtained. The three-dimensional domains for the discrete and continuum models are 10.0 by 10.0 by 10.0 meters. The input parameters for the discrete simulation are listed in Table 5.1. Fracture set one dips at 60° towards the minimum X boundary of the domain, while fracture set two is horizontal. Two hundred fractures were generated for both fracture sets. The input mean fracture length for both sets (0.8 m) is three times the length used in the two-dimensional model. The widths of the fractures, in the Y direction, has the same distribution as the fracture lengths. Figure 5.4 shows particle plots for the discrete and continuum models. The input location noted by the asterisk is at (2.58, 2.50, 2.51). The three-dimensional distributions of the solute for both models are shown in four XZ slices through the domain (Figure 5.4). All particles between the limits given for the Y coordinates of the slice are transposed onto the one plane. The size of the points plotted is constant thus occurrences of multiple particles at the same location are not represented. The time of the simulation was set to 25 days. Visual comparison of the plumes in Figure 5.4 shows a reasonable match of the solute distribution in the two models. No actual numerical comparison, such as mean and variance, of the distribution was carried out due to limited research time. The match between the two solute distributions suggests that the statistics chosen here are a good first step to determining the best statistics to collect for representing particle motion in discrete network. The discrete simulation tested here required 27 Megabytes of Ram storage and took approximately 5 hours to run on a Sun4 SparcStation. 130 Table S.l Summary of input parameters for three-dimensional discrete model. Parameter Value Fracture Domain xbl) length in the X-direction [ybl) length in the Y-direction zbl) length in the Z-direction General Fracture Parameters dxinit) grid spacing in local X-direction dyinit) grid spacing in local Y-direction nsets) number of fracture sets Fracture Set 1 ;nf) number of fractures > ) mean dip angle > _g) standard deviation in dip angle mean fracture length minfl) minimum fracture length >f ) mean fracture width minfw) minimum fracture width > ) log mean fracture aperture > ) log std. dev. fracture aperture l ) correlation length in local X \py) correlation lenght in local Y ang a w ap ap a p x 10.0 m 10.0 m 10.0 m 0.25 m 0.25 m 2 200 60.0 0.50 0.8 m 2.40 m 0.8 m 2.40 m -4.5 m 0.0 m 0.0 m 0.0 m Fracture Set 2 nf) number of fractures > _) mean dip angle • ang) °dard deviation in dip angle UQ) mean fracture length minfl) minimum fracture length Pr ) mean fracture width minfw) minimum fracture width > ) log mean fracture aperture > ) log std. dev. fracture aperture l ) correlation length in local X l ) correlation lenght in local Y an a sta w ap ap a p x a p y 200 0.0 0.50 0.8 m 2.40 m 0.8 m 2.40 m -4.5 m 0.0 m 0.0 m 0.0 m General Parameters mc) number of monte carlo realizations np) number of particles At) time step size Ah/AX) gradient in direction of flow 1 5000 25 days 0.01 131 0 meters F i g u r e 5.1 X Z s l i c e t h r o u g h 3 - D 5 - 1 0 network at Y= 5.00 m. 132 Figure 5.2 Schematic diagram of the three—dimensional fractures to illustrate the division of fractures into quadrants, the local coordinate systems within the quadrants and the boundary conditions of the quadrants. 133 Figure 5.3 Schematic diagram of particle motion to illustrate the collection of particle motion statistics in 3-D. DISCRETE MODEL * 25 days —> 0 - 2.5 • 25 days 0° 2.5 - 5 5 - 7.5 CONTINUUM MODEL 125 days ~ > 0° 0 5 10 meters Figure 5.4 Particle plots for the 3-D sub-REV simulations. 135 6.0 CONCLUSIONS AND DISCUSSION 6.1 Conclusions 1. Full-domain testing of the original Schwartz and Smith [1988] model, gaussian(O), bore out the good match that they found between the initial portion of the breakthrough curves for the stochastic-continuum and discrete models. It was, however, noted that the stochastic-continuum model was unable to reproduce the long tail of the plumes simulated in the discrete model. The difference in the plumes was more noticeable in the plots of particle distribution and plots of the temporal development of the spatial moments. 2. When variability in fracture aperture was added to the model, a large increase in the difference between the second moments of the stochastic-continuum and discrete simulated plumes was noticed. Analysis of the velocity statistics collected in the discrete sub-domain showed that the distribution was skewed to lower velocity values. To better fit the velocity data the assumption of a gaussian distribution was replaced with that of a gamma distribution. Further analysis showed that there was a high degree of correlation between solute velocities in successive fracture segments when the successive fracture segments are of the same orientation. The addition of the gamma velocity distribution and velocity correlation structure greatly improve the ability of the stochastic-continuum model to match the second moments of the discrete model plumes. 3. Comparisons were made between the spatial plume moments from the extended stochastic-continuum model, gamma(2), and the discrete network model for various network conditions. The far-field percentage difference in the first moments were less than 10%, while the far-field percentage difference in the second moments ranged from 6% to 30%. 136 4. The difference in the far-field second moments increased for hydraulic gradient orientations which were at an angle to the orientation of both the fracture sets. The increase is possibly due to correlation in successive flow path velocities which are not in the same direction. This cross-correlation is not accounted for in the stochasticcontinuum model as tested. 5. Of the two additional fracture networks tested, the stochastic-continuum model performed well for the network were the hydraulic gradient was parallel to the longer fracture set. For the network where the hydraulic gradient was parallel to the shorter fracture set, results were inconclusive as the domain was too small to obtain representative results. In general the match of the spatial moments between the stochastic-continuum and discrete models improves as the length of the fractures parallel to the hydraulic gradient increase relative to the cross-fractures. 6. Further testing of the two-dimensional model is warranted under different network conditions, (non-orthogonal fracture sets, 3 or more fracture sets, etc.), however, such testing may be incorporated in future work on the 3-D model. 7. Initial results from a three-dimensional model give some indication of what statistics will represent particle motion in three-dimensions. Statistics collected in a single realization of the discrete model were used to approximately recreate the discrete network plume in a stochastic-continuum model. 6.2 Discussion Application of the model has been primarily focussed on predictions of solute transport in fractured media. Other possible applications of the stochastic-continuum method may include: 1) Analysis of the distribution of grout injected into a fractured rock body and 2) Withdrawal of contaminants from a fractured rock body by pumping. 137 The remainder of this discussion will focus on how to apply the stochastic-continuum method to a potential field site for prediction of solute transport. The first step in applying the model to an actual site is to statistically describe the fractures at the site. Necessary fracture parameters are orientation, frequency, length and aperture. Fracture orientations can be measured from surface and underground exposures and from orientated drill core. The poles to the fracture planes are typically plotted on lower-hemisphere equal area nets. By contouring these plots the mean and variance in the orientation of the fracture sets can be determined (Rouleau and Gale, 198S). Some interpretation must be done to decide upon the numberof fracture sets which best represent the data. Several methods of measuring fracture frequency or spacing have been presented in the literature (Snow, 1969; Priest and Hudson, 1976; Cruden, 1977; Hudson and Priest, 1979; and Rouleau and Gale, 1985). The stochastic-continuum models as presented here assume a random fracture spacing , however, this could easily be changed to fit other 4 distributions. Measurements of fracture lengths are best made from fracture trace maps, such as Figure 1.1. Rouleau and Gale [1985] describe a method of measuring fracture trace lengths in underground exposures where censoring of the data is a problem, as fractures disappear into the walls of the underground opening. As discussed earlier, fracture trace length distributions are typically skewed towards higher lengths. The method of inputting fracture lengths in the model can be adjusted to fit the distribution of fracture lengths measured in the field. Measuring fracture apertures is somewhat problematic. In obtaining access to the fracture the in situ stress across the fracture opening is usually decreased. Under Fracture centers are randomly distributed. 138 decreased stress the fracture opening, or aperture, will often increase in size. Several authors have suggested methods of physically measuring fracture apertures (Snow, 1969; and Brown et al. 1986). A different approach is to measure the effective hydraulic aperture from in situ or laboratory flow tests (Snow, 1967; Gale, 1982; and Raven, 1986). From a combination of these measurements an approximate distribution for fracture aperture may be obtained, and input into the model. Once all the pertinent fracture parameters distributions have been determined, they are entered into the discrete sub-domain portion of the model. By sensitivity analysis, the minimum size of the discrete sub-domain to obtain representative particle motion statistics is determined. The dimensions and boundary conditions of the continuum portion of the model are set to simulate the domain of interest. Particles are introduced into the model at the desired location and moved through the domain by rules based on the statistics collected in the discrete sub-domain model. Predictions of the mass distribution from the stochastic-continuum model must be tempered with some idea of the accuracy of the model (eg. the 30% underestimate in the variance for the two-dimensional model). Errors due to the approximations in the choice of parameter distributions used to simulated the fractures must also be taken into account. A third source of error is the spatial variation in the fracture parameters within the domain. A fracture parameter, such as aperture, may not be homogeneously distributed, as is assumed in the continuum method, but may rather trend to larger or smaller values in a given direction. This assumption of homogeneity of the fracture parameters is a basic assumption of all continuum methods. Whether or not it is a reasonable assumption is a question which must be answered for each individual site. 139 BIBLIOGRAPHY Abelin, H., Migration in a Natural Fracture: An In Situ Experiment * Natural Fracture. Ph.D. Thesis Dept. of Chem. Eng., Royal Institure of Technology. S 100, 44, Stockholm, Sweden, 1986. na Barton, C . C , and P.A. Hsieh, Physical and Hvdrologic-Flow Properties of Fractures. Inter. Geol. Congress Field Trip T385, 1989. Bear, J., and A. Verruijt, Modeling Groundwater Flow and Pollution. D. Reidel Publishing Company, Dordrecht, Holland, 1987. Becker, R.A., and J.M. Chambers, S: An Interactive Enviroment for Data Analysis and Graphics. Wadsworth, Belmont, CA., 1984. Box, G., and G. Jenkins, Time Series Analysis Forecasting and Control. Holden-Day, San Francisco, pp. 575, 1976. Brown, S.R., and C H . Scholz, Broad Bandwidth Study of the Topography of Natural Rock Surfaces. J. Geophys. Res.. 90(B14), pp 12,572-12,582, 1985. Brown, S.R., R.L. Kranz, and B.P. Bonner, Correlation Between the Surfaces of Natural Rock Joints. Geophys. Res. Letters, 13(13), pp. 1430-1433, 1986. Brown, S.R., Fluid Flow Through Rock Joints: The Effect of Surface Roughness. J. Geophys. Res.. 92(B2), pp. 1337-1347, 1987. Cacas, M . C , E. Ledouc, G. de Marsily, B. Tillie, A. Barbreau, E. Durand, B. Feuga, and P. Peaudecerf, Modelling Fracture Flow with a Stochastic Discrete Fracture Network: Calibration and Validation 1.. Water Resour. Res., 26(3), pp. 479-489, 1990. Clifton, P.M., and S.P. Neuman, Effects of Kriging and Inverse Modelling on Conditional Simulation of the Avra Vallev Aouifer in Southern Arizona. Water Resour. Res., 18(4), pp. 1215-1234, 1982. Cruden, D.M., Describing the Size of Discontinuities. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr., 14, pp. 133-137, 1977. de Josselin de Jong, G., and S.C. Way, Dispersion in Fissured Rock, unpublished report, 30 pp., N.M. Inst, of Min. and Technol., Socorro, 1972. Duff, I., A. Erisman, and J. Reid, Direct Methods for Sparse Matrices. Oxford University Press, Oxford, 1987. Dverstorp, B., and J. Andersson, Application of the Discrete Fracture Network Concept with Field Data: Possibilities of Model Calibration and Validation. Water Resour. Res., 25(3), pp. 540-550, 1989. Endo, H.K., J.C.S. Long, C R . Wilson, and P.A. Witherspoon, A Model for Investigating Mechanical Transport in Fracture Networks. Water Resour. Res., 20(10), pp. 1390-1400, 1984. 140 Freeze, R . A . , and J . A . Cherry, Groundwater. Prentice-Hall, Inc., Englewood Cliffs, N . J . , U.S.A., 1979. Gale J . E . , Assessing the Permeability Characteristics of Fractured Rock. Geol. Soc. of America Special Paper 189, 1982. Gaver, D.P., and P.A. Lewis, First-Order Autoregressive Gamma Sequences and Point Processes. Advanced Applied Problems, 12, pp. 727-745, 1980. Hudson, J . A . , and S.D. Priest, Discontinuities and Rock Mass Geometry. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr., 16, pp. 339-362, 1979. Hull, L . C . , J.D. Miller, and T . M . Clemo, Laboratory and Simulation Studies of Solute Transport in Fracture Networks. Water Resour. Res., 23(8), pp. 15051513, 1987. Huyakorn, P.S., and G . F . Pinder, Computational Methods in Subsurface Flow. Academic Press, Inc., London, 1983. Lettenmaier, D.P., and S.J. Burges, A n Operational Approach to Preserve Skew in Hvdrologic Models of Long-Term Persistence. Water Resour. Res., 13(2), pp. 281-290, 1977. Long, J.C.S., J.S. Remer, C R . Wilson, and P.A. Witherspoon, Porous Media Equivalents for Networks of Discontinuous Fractures. Water Resour. Res., 18(3), pp. 645-658, 1982. Long, J.C.S., and P.A. Witherspoon, The Relationship of the Degree of Interconnection to Permeability in Fracture Networks. / . Geophys. Res., 90(B4), pp. 30873098, 1985. Narr W., and J.B. Currie, Origin of Fracture Porosity - Examples from Altamont Field. Utah. The Amer. Assoc. of Petrol. Geologists Bull., 66(9), pp. 1231-1247, 1982. Neretnieks, I., The Influence of Microfissures in Crystalline Rock on Radionuclide Migration. Predictive Geology with Emphasis on Nuclear-Waste Disposal, G . de Marsily, D . Merrian ed. Pergamon Press, 1982. Press, W.H., B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, Numerical Recipes in C . The Art of Scientific Computing. Cambridge University Press, Cambridge, 1988. Priest, S.D., and J . A . Hudson, Discontinuity Spacing in Rock. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr., 13, pp.135-148, 1976. Raven, K . G . , Hydraulic Characterization of a Small Ground-Water Flow System in Fractured Monzonitic Gneiss. National Hydrology Research Institute, Cananda, Paper No. 30, IWD Scientific Series No. 149, 1986. Rice, J . A . , Mathematical Statistics and Data Analysis. Wadsworth and Brooks/Cole, Belmont, C A . 1988. 141 Rouleau, A . , and J . E . Gale, Statistical Characterization of the Fracture System in the Stripa Granite. Sweden. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr., 22, pp. 353-367, 1985. Rouleau, A . , and J . E . Gale, Stochastic Discrete Fracture Simulation of a Groundwater Flow in to an Underground Excavation in Granite. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr., 24(2), pp. 99-112, 1987. Rubenstein, R . V . , Simulation and the Monte Carlo Method. Wiley, New York, 1981. Schwartz, F.W., L . Smith, and A.S. Crowe, A Stochastic Analysis of Macroscopic Dispersion in Fractured Media. Water Resour. Res., 19(5), pp. 1253-1265, 1983. Schwartz, F.W., and J . L . Smith, A Continuum Approach for Modeling Mass Transport in Fractured Media. Water Resour. Res., 24(8), pp. 1360-1372, 1988. Shimo, M . , and J.C.S. Long, A Numerical Study of Transport Parameters in Fracture Networks. Flow and Transport Through Unsaturated Fractured Rock, AGU Monograph, 42, pp. 121-131, 1987. Smith, L . , and F.W. Schwartz, A n Analysis of the Influence of Fracture Geometry on Mass Transport in Fractured Media. Water Resour. Res., 20(9), pp. 12411252, 1984. Smith, L . , C W . Mase, F.W. Schwartz, and D . Chorley, A Numerical Model for Transport in Networks of Planar Fracture. Hydrology of Rocks of Low Permeability, International Association of Hydrologists Memoires, Volume X V I I , Part 2 Proceedings, Tucson, Arizona, Congress, 1985. Snow, D . T . , The Frequency and Aperture of Fractures in Rock. Int. J. Rock Mech. Min. Sci., 7, pp. 23-40, 1970. Tsang, Y.W., and C . F . Tsang, Channel Model of Flow Through Fractured Media. Water Resour. Res., 23(3), pp. 467-479, 1987. Wheatcraft, S.W. and S.W. Tyler, A n Explanation of Scale-Dependent Disoersivitv in Heterogeneous Aquifers Using Concents of Fractal Geometry. Water Resour. Res., 24(4), pp. 566-578, 1988. Zienkiewicz, O . C , and K . Morgan, Finite Elements and Approximation. John Wiley and Sons, New York, 1983.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A statistical continuum approach for mass transport...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A statistical continuum approach for mass transport in fractured media Robertson, Mark Donald 1990
pdf
Page Metadata
Item Metadata
Title | A statistical continuum approach for mass transport in fractured media |
Creator |
Robertson, Mark Donald |
Publisher | University of British Columbia |
Date Issued | 1990 |
Description | The stochastic-continuum model developed by Schwartz and Smith [1988] is a new approach to the traditional continuum methods for solute transport in fractured media. Instead of trying to determine dispersion coefficients and an effective porosity for the hydraulic system, statistics on particle motion (direction, velocity and fracture length) collected from a discretely modeled sub-domain network are used to recreate particle motion in a full-domain continuum model. The discrete sub-domain must be large enough that representative statistics can be collected, yet small enough to be modeled with available resources. Statistics are collected in the discrete sub-domain model as the solute, represented by discrete particles, is moved through the network of fractures. The domain of interest, which is typically too large to be modeled discretely is represented by a continuum distribution of the hydraulic head. A particle tracking method is used to move the solute through the continuum model, sampling from the distributions for direction, velocity and fracture length. This thesis documents extensions and further testing of the stochastic-continuum two-dimensional model and initial work on a three-dimensional stochastic-continuum model. Testing of the model was done by comparing the mass distribution from the stochastic-continuum model to the mass distribution from the same domain modeled discretely. Analysis of the velocity statistics collected in the two-dimensional model suggested changes in the form of the fitted velocity distribution from a gaussian distribution to a gamma distribution, and the addition of a velocity correlation function. By adding these changes to the statistics collected, an improvement in the match of the spatial mass distribution moments between the stochastic-continuum and discrete models was effected. This extended two-dimensional model is then tested under a wide range of network conditions. The differences in the first spatial moments of the discrete and stochastic-continuum models were less than 10%, while the differences in the second spatial moments ranged from 6% to 30%. Initial results from the three-dimensional stochastic-continuum model showed that similar statistics to those used in the two-dimensional stochastic-continuum model can be used to recreate the nature of three-dimensional discrete particle motion. |
Subject |
Hydraulic fracturing -- Statistical methods |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-11-02 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0052808 |
URI | http://hdl.handle.net/2429/29740 |
Degree |
Master of Applied Science - MASc |
Program |
Geological Sciences |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1990_A7 R62.pdf [ 6.77MB ]
- Metadata
- JSON: 831-1.0052808.json
- JSON-LD: 831-1.0052808-ld.json
- RDF/XML (Pretty): 831-1.0052808-rdf.xml
- RDF/JSON: 831-1.0052808-rdf.json
- Turtle: 831-1.0052808-turtle.txt
- N-Triples: 831-1.0052808-rdf-ntriples.txt
- Original Record: 831-1.0052808-source.json
- Full Text
- 831-1.0052808-fulltext.txt
- Citation
- 831-1.0052808.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052808/manifest