UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
SOLUTE TRANSPORT IN HETEROGENEOUS POROUS FORMATIONS By GEORGE GARY DEMMY, JR. A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1999 For my parents. ACKNOWLEDGMENTS Dissertations have but one name in the author slot, but anyone who has written one will readily testify that it represents the work and support of myriad people. Mary Hall freed my mind from many of the administrative hassles students encounter. Dawn Mendoza appointed, unappointed, and reappointed me numerous times throughout my tenure, and made the business end of things run as painlessly as possible. Ralph Hoffmann kept me in keys and desks and office space, and helped me get things wired when they needed wiring. Wally Hunter, system administrator and Unix hacker, deserves as much praise as I can heap upon him for keeping the network up and running and the data contained with in it safe and available. The quality of an institution and its faculty is reflected in the students they attract. I have enjoyed interacting with the students of the Hydrologic Sciences Academic Cluster, and this interaction has enhanced my efforts. Liyong Li, Yan Zhang, and Xavier Foussereau shared the same computer lab trenches as I, and my work cannot be separated from their influence. Many rewarding research experiences have been associated with collaborations made with Jim Jawitz. Jim is a colleague whose outlook and expertise is somewhat different than my own, but it is precisely these differences that have made the collaborations so rewarding. I shall seek to continue working with Jim through out my career. Dr. Andrew I. James provided access to computer codes he developed in the course of his own research and technical assistance. ill Dr. Kenneth Campbell deserves special recognition not only for his exemplary service on my committee, but also for the example he has set as an engineer, a scientist, and a professor. I have appreciated Dr. Mike Annable's wit and example through the years and have enjoyed interacting with him in a wide variety of settings, whether in the classroom, on the football field, or elsewhere. Dr. Rao has always challenged me with questions that require looking beyond a simple formula and that evade a fillintheblank answer. Should I become scientist of any repute, Dr. Rao would be the man to blame. Dr. Kirk Hatfield is my source. He has been a source of wisdom, sympathy, advice, and funding. I have greatly enjoyed working with Dr. Hatfield as both teaching and research assistant. Dr. Gia Destouni taught a few stochastic subsurface hydrology lectures in the spring of 1994 that would forever change my thinking. She introduced the Lagrangian approach to transport with which I have been obsessed ever since. Dr. Destouni pursues fruitful collaboration vigorously, and this aggressive, yet collaborative, style is one that I now seek to emulate. Dr. Vladimir Cvetkovic has generously hosted three visits to the Kungl Tekniska H6gskolan in Stockholm, Sweden. His considerable generosity flows from his intense interest in stochastic Lagrangian transport. Dr. Sten Berglund has been a friend and mentor through this arduous journey. It is no exaggeration to say that without his help and encouragement this document would not exist. Dr. Wendy Graham deserves special recognition for her tireless faith and support. She has allowed me a free reign to pursue whatever interests I chose. In hindsight, this method coupled, with my own personality quirks and spectrum of interests, may not iv be the most efficient for producing degrees and papers, but it has produced several interesting adventures. George and Ellen Demmy, my parents, have provided love and all manner of support throughout my life, academic and otherwise. They both instilled in me a love of learning and an insatiable curiosity about all things that serve me well, and will remain among my best traits. Finally, I'd like to extend my thanks and love to Celine Bufkin, for her love, encouragement, support, and vision. v TABLE OF CONTENTS page ACKNOWLEDGMENTS ........................... .... iii LIST OF TABLES ......................... ..... ... .... ix LIST OF FIGURES ............... .... .......... x ABSTRACT ................. ..................... xv CHAPTERS 1 INTRODUCTION ........................... 1 1.1 Primary Contribution ............ ..... ..... 1 1.2 Historical Context ........................ 1 1.3 Research Goals ......................... 3 1.3.1 Lagrangian Velocity Mean and Covariance ....... 3 1.3.2 Injection Mode Analysis .................. 4 1.3.3 Heterogeneous Reaction Parameter .......... 4 1.3.4 Application to "Reality" ................. 4 2 LAGRANGIANEULERIAN FLOW FIELD RELATIONSHIP .. 5 2.1 Introduction .................. ....... 5 2.2 Flow Field Partitioning Strategies .............. 7 2.2.1 Simple Subdivision .................. 8 2.2.2 General Subdivision .................. 10 2.3 Descriptions of Fields ................... ... 12 2.4 Numerical Experiments ................... .. 13 2.4.1 The Log Conductivity Field .............. 14 2.4.2 Velocity Field Generation ............... 18 2.4.3 Trajectory Generation ................. 19 2.4.4 Monte Carlo Simulation ................ 19 2.4.5 Sampling Strategy ................... 20 2.4.6 Simulation Results ................... 22 2.4.7 Log Conductivity Statistics .............. 23 2.4.8 Velocity Statistics ................... 29 2.4.9 Head and Head Gradient ............... 40 2.5 Theory .............................. 40 2.5.1 Effective Conductivity ................. 40 vi 2.5.2 Eulerian Velocity Field ................ 42 2.5.3 Streamtube Weights .................. 48 2.5.4 Control Plane Oriented Velocity Moments .... ... 51 2.5.5 Control Time Oriented Velocity Moments ....... 55 2.5.6 Log Conductivity Moments .............. 59 2.6 Evaluation of Simulation and Theoretical Results ....... 60 2.6.1 Mean Velocities in Space ............... 60 2.6.2 Velocity Variances in Space .............. 60 2.6.3 Mean Velocities in Time ................ 61 2.6.4 Velocity Variances in Time .............. 62 2.6.5 Log Conductivity Statistics .............. 62 2.7 General Conclusions .................. 63 3 NONREACTIVE SOLUTE TRANSPORT ............. 65 3.1 Introduction .................. ....... 65 3.2 Injection Mode and Streamtube Selection .......... 65 3.2.1 Uniform Resident Injection .............. 65 3.2.2 Injection in Flux ................... .. 67 3.3 Integrated Streamtube Properties .............. 67 3.3.1 Displacements ................... ... 68 3.3.2 Travel Times ................... .... 69 3.4 Solute Transport ................... ...... 70 3.5 Evaluation of Expressions ................... 71 3.5.1 Displacements ................... ... 71 3.5.2 Travel Times ................... .... 72 3.6 General Conclusions ................... .... 75 4 REACTIVE SOLUTE TRANSPORT ................ 79 4.1 Introduction ...... ................... .. 79 4.2 Temporal Moments ................... .... 79 4.3 Conceptual Overview of Solute Transport .......... 80 4.4 Simulation of Reactive Solute Transport ........... 83 4.4.1 Reaction Parameter .................. 84 4.4.2 Reaction Flow Path .................. 85 4.4.3 Reaction Flow PathTravel Time Cross Moment 93 4.5 Theory of Reactive Solute Transport ............. 93 4.5.1 Mean Reaction Flow Path ............... 96 4.5.2 Reaction Flow Path Variance ............. 97 4.6 Discussion ................... ....... 98 4.7 Spatial Moments. ................... ..... 99 4.7.1 Injection Mode ................... ... 105 4.7.2 Discussion ................... ..... 106 vii 5 APPLICATION OF CONCEPTS ............ ........ 108 5.1 Introduction .......... ................. 108 5.2 Theory .............. .... ............ 108 5.2.1 Spatial Moments ................ . 109 5.2.2 Relative Degradation ............ ... 110 5.3 Application: Simulation .............. . 111 5.3.1 Results ....... ........ .. .. 113 5.3.2 Discussion ........ . . . .... 113 5.4 Application: Field Data .................. 115 5.4.1 Results and Discussion: Center of Mass ....... 116 5.4.2 Results and Discussion: Degradation ......... 117 5.5 General Conclusions ............... . 117 6 GENERAL CONCLUSIONS ................... ... 119 APPENDICES A EULERIAN VELOCITY DISTRIBUTION PROPERTIES .... 125 B VELOCITY LOG CONDUCTIVITY CORRELATION ...... 127 C CONSISTENT FIRST ORDER APPROXIMATION ....... 131 C.1 Travel Time ................... ........ 131 C.2 Displacements .......................... 133 REFERENCES ................................... 134 BIOGRAPHICAL SKETCH ................... ........ 138 V111iii LIST OF TABLES Table page 2.1 Comparison of target log conductivity values and observed values along areaweighted Eulerian trajectories. ................ 26 2.2 Comparison of target log conductivity values for areaweighted Eulerian trajectories in space and observed areaweighted Lagrangian trajectories in time. ......... ... .... ........ .. 29 2.3 Comparison of a priori Eulerian velocity variance to observed values. The decimal representation of the exact fractional values of the predicted variance are provided for convenience. ........... 35 2.4 Comparison of observed log hydraulic conductivity variances to estimators. Observed value is zeroslope line regressed through data. The decimal representation of the fractional value is included for convenience ................ .. ..... .......... 63 5.1 The center of mass in meters of steady plumes estimated by traditional trapezoidal integration of spatial concentration data and by the mean v/k estimator Equation (5.7). .................. .. 116 5.2 Degradation rates relative to naphthalene and absolute degradation rates in reciprocal days as estimated by Equation (5.9) and as reported by MacIntyre et al. [1993]. NR denotes results not reported in MacIntyre et al. [1993]. Naphthalene degradation reported to be 0.0064 d1 .................. ...... ....... 117 ix LIST OF FIGURES Figure page 2.1 Schematic representation of an aquifer. ................. 6 2.2 Simple subdivision of the flow field by one of four strategies: A) equalarea straight line, B) equalarea streamline, C) equalflow streamline, D) equalflow straight line. .................. .... 9 2.3 Four different partitioning strategies for the same flow field: A) area weight Eulerian trajectory, B) area weight Lagrangian trajectory, C) flow weight Lagrangian trajectory, D) flow weight Eulerian trajectory .. .. . ... .. .. .. 11 2.4 Log conductivity distribution from a 41400 node realization generated by the turningbands method compared to the standard normal distribution. Target distribution parameters are ln k = 0 and U~k = 1. 16 2.5 Estimated log conductivity longitudinal covariance estimated from a realization generated by the turningbands method compared to an exponential covariance function. Variance of log conductivity an k = 1 and correlation length Aln k = 1. Spacing for log conductivity values was Alnk/5 ....... ...... ... .... .. ...... 17 2.6 Mean log conductivity as function of displacement in mean flow direction for the four different trajectory collections and Ulnk = 1/2 (a) and a k = 1 (b). ........................... 24 2.7 Log conductivity variance as function of displacement in mean flow direction for the four different trajectory collections and Ulnk = 1/2 (a) and Ink = 1 (b). Variances are normalized by target variance of Eulerian log conductivity field an.k. .................. 25 2.8 Mean log conductivity as function of travel time to control planes along mean flow direction for the four different trajectory collections and nk = 1/2 (a) and l nk= 1 (b) ..................... .... 27 2.9 Log conductivity variance as function of travel time to control planes along mean flow direction for the four different trajectory collections and a1 k = 1/2 (a) and Ink = 1 (b). Variances are normalized by target variance of Eulerian log conductivity field n k. . ... 28 U~YbVYIYI~IIVV VI UIUI~IL ~b VI~l~~UIILI 1~;1l 'Ck 2.10 Comparison of positiondependent mean velocity estimators to data for the four collections of trajectories for ank = 1/2 (a) and r2nk = 1 (b). Velocities are normalized by U .. ............... 30 2.11 Comparison of the positiondependent velocity variance estimators to data for the four collections of trajectories for aInk = 1/2 (a) and aOnk = 1 (b). The variances are normalized by U2bunk. ........ 31 2.12 Comparison of spatial Eulerian (AW/ET) and Lagrangian trajectory (FW/LT) velocity covariance data to covariance model. Data are taken from the a nk = 1 set of realizations. Covariances presented as dimensionless lag correlation functions. ......... ..34 2.13 Comparison of timedependent mean velocity estimators to data for the four collections of trajectories for ank = 1/2 (a) and Ink = 1 (b). Velocities are normalized by U........... . . 36 2.14 Comparison of timedependent velocity variance estimators to data for the four collections of trajectories for aonk = 1/2 (a) and ank = 1 (b). The variances are normalized by U2bok. ............. 37 2.15 Comparison of temporal Lagrangian trajectory (AW/LT) velocity covariance data to covariance model. Data are taken from the aUnk 1 set of realizations. Covariances presented as dimensionless lag correlation functions. ....... ............ ....... 38 2.16 Normalized log initial velocities for twodimensional model aquifer compared to the standard normal distribution. The random variable was normalized by the sample mean and standard deviation to illustrate the lognormality. ................ . . 39 2.17 Comparison of velocity covariance expression Equation (2.25) to the approximation Equation (2.13) for aInk = 1. Covariances are presented as dimensionless correlation functions. ............ 45 2.18 Components of the velocity covariance function. ............ 46 3.1 Comparison of mean displacements for the two injection modes or streamtube collections. The flowweighted Lagrangian trajectory (FW/LT) statistics correspond to a flux injection, and the areaweighed Lagrangian trajectory (AW/LT) statistics correspond to a uniform resident injection. The top figure contains the aink = 1/2 data and the bottom figure contains the an k = 1 data. ....... 73 xi 3.2 Comparison of displacement variances for the two injection modes or streamtube collections. The flowweighted Lagrangian trajectory (FW/LT) statistics correspond to a flux injection, and the areaweighed Lagrangian trajectory (AW/LT) statistics correspond to a uniform resident injection. The top figure contains the oan4k = 1/2 data and the bottom figure contains the Unk = 1 data. Variances are normalized by AnkUnk. ....... ... .......... 74 3.3 Comparison of mean arrival times for the two injection modes or streamtube collections. The flowweighted Lagrangian trajectory (FW/LT) statistics correspond to a flux injection, and the areaweighed Lagrangian trajectory (AW/LT) statistics correspond to a uniform resident injection. The top figure contains the a k = 1/2 data and the bottom figure contains the a k = 1 data. .............. 76 3.4 Comparison of arrival time variances for the two injection modes or streamtube collections. The flowweighted Lagrangian trajectory (FW/LT) statistics correspond to a flux injection, and the areaweighed Lagrangian trajectory (AW/LT) statistics correspond to a uniform resident injection. The top figure contains the an k = 1/2 data and the bottom figure contains the ank = 1 data. Variances are normalized by (U/(Al\kaClnk))2 .. .................... .....77 3.5 Travel time estimate error as function of displacement associated with assuming a uniform resident injection when injection is in flux for aUnk = {1/2,1} ................ ...... .. ...... 78 4.1 Mean reaction parameter values for different values of correlation parameter / and different trajectory collections for the ank = 1 set of simulations. Top figure is P = {1, +1}. Middle figure is / = {1/2,+1/2}. Bottom figure is / = 0 (no correlation). Trajectory collections are areaweighted Lagrangian (AL), flowweighted Lagrangian (AF) and areaweighted Eulerian (AE). ......... 86 4.2 Mean reaction parameter values for different values of correlation parameter / and different trajectory collections for the ank = 1/2 set of simulations. Top figure is / = {1,+1}. Middle figure is = { 1/2, +1/2}. Bottom figure is / = 0 (no correlation). Trajectory collections are areaweighted Lagrangian (AL), flowweighted Lagrangian (AF) and areaweighted Eulerian (AE). ......... 87 4.3 Mean reaction path observed values and estimators for different values of correlation parameter / and different trajectory collections for the nk = 1 set of simulations. Top figure is / = {1, +1). Middle figure is / = {1/2, +1/2}. Bottom figure is 3 = 0 (no correlation). Trajectory collections are areaweighted Lagrangian (AL) and flowweighted Lagrangian (AF). .................. ... 89 xii 4.4 Mean reaction path observed values and estimators for different values of correlation parameter p and different trajectory collections for the k = 1/2 set of simulations. Top figure is / = {1, +1}. Middle figure is P = {1/2, +1/2}. Bottom figure is 3 = 0 (no correlation). Trajectory collections are areaweighted Lagrangian (AL) and flowweighted Lagrangian (AF). ..... ... ................. 90 4.5 Reaction path variance observed values and estimators for different values of correlation parameter / and different trajectory collections for the CUn k = 1 set of simulations. Top figure is 3 = {1, +1}. Bottom figure is 3 = {1/2, +1/2}. Trajectory collections are areaweighted Lagrangian (AL) and flowweighted Lagrangian (AF). .. 91 4.6 Reaction path variance values and estimators for different values of correlation parameter p and different trajectory collections for the crnk = 1/2 set of simulations. Top figure is 3= {1, +1}. Bottom figure is / = {1/2, +1/2}. Trajectory collections are areaweighted Lagrangian (AL) and flowweighted Lagrangian (AF) ......... 92 4.7 Correlation between the reaction flow path and travel time for different values of correlation parameter / and different trajectory collections for the crank = 1 set of simulations. Top figure is / = {1, +1}. Bottom figure is /3 = {1/2, +1/2}. Trajectory collections are areaweighted Lagrangian (AL) and flowweighted Lagrangian (AF). 94 4.8 Correlation between the reaction flow path and travel time for different values of correlation parameter /3 and different trajectory collections for the cank = 1/2 set of simulations. Top figure is / = {1, +1}. Bottom figure is / = {1/2, +1/2}. Trajectory collections are areaweighted Lagrangian (AL) and flowweighted Lagrangian (AF). .. 95 4.9 Assumed timedependent log conductivity covariance function along Lagrangian trajectories (Equation (4.29)) compared to simulation data for variances of log conductivity aink = {1/2, 1}. The Lagrangian velocity covariance function is Equation (2.61). ...... 103 4.10 Hypothesized timedependent Eulerian log conductivity covariance functions for Un2k = {1/2, 1} compared to simulation data and Equation (4.29). ................. ............. 104 4.11 Mean reaction flow path as a function of nonreactive travel time compared to estimator Equation (4.25) for correlation parameter values /3= {1,0,1}. .................. .... ...... 104 4.12 Reaction flow path variance a function of nonreactive travel time for correlation parameter values 3 = {1, 0, 1} .............. 105 xiii 4.13 Mean reaction flow path as a function of nonreactive travel time compared to estimator Equation (4.34) for correlation parameter values = {1/2, 0, 1/2} (top figure) and 0 = {1, 0, 1} (bottom figure). 2nk =1................. .............. 107 5.1 Resident concentration profile averaged along a transect orthogonal to the mean flow direction for plumes for ornk = {0.0, 0.1,0.5, 1.0}. .. 114 B.1 Comparison of velocity and log conductivity correlation functions from Graham and McLaughlin [1989b] and Rubin and Dagan [1992] and a simplified exponential model. All functions are normalized by KJaIk/O. The apparent "roughness" of the Rubin and Dagan [1992] function is due to estimation errors associated with inferring their function from graphical data. .................. 128 B.2 Comparison of the simplified velocitylog conductivity correlation function to simulation data for or = {1/2, 1} ................ 129 xiv Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy SOLUTE TRANSPORT IN HETEROGENEOUS POROUS FORMATIONS By George Gary Demmy, Jr. December 1999 Chair: Wendy Graham Major Department: Agricultural and Biological Engineering This work quantifies relationships between the spatial, or Eulerian, distribution of the properties of a chemically and physically heterogeneous porous medium and those as observed along the natural, or Lagrangian, trajectories that a fluid particle traces in a steady and irrotational flow. From these relationships, expressions that relate the transport of solutes through the porous medium along the natural trajectories to the aforementioned Eulerian distributions are developed. The effects of injection mode upon global measures of transport as reflected by the temporal moments of breakthrough curves and spatial moments of a solute plume are developed. The coupled effects of correlation of a linear equilibrium sorption to the underlying log hydraulic conductivity field and injection mode on the evolving temporal moments of mass breakthrough curve and the coupled effects of correlation of a firstorder decay coefficient and injection mode upon the spatial moments of a solute plume are examined. xv CHAPTER 1 INTRODUCTION 1.1 Primary Contribution The primary contribution of this work is documenting a consistent approach to the development of estimators of solute transport in the Lagrangian framework. This consistent approach is to identify the collection of Lagrangian trajectories associated with a stationary Eulerian velocity field for which the Lagrangian velocities are also stationary. This is a modest extension of work pioneered by myriad hydrologists. The benefit derived from this contribution is enhanced understanding of transport of solutes in heterogeneous velocity fields. 1.2 Historical Context A grant to study the impact of heterogeneous source morphology upon subsequent solute transport initiated this research. In seeking a suitable framework and "toolset" with which to attack this problem, there appeared what seemed to be fundamental inconsistencies in the literature. Given the peculiar tendencies of the author, the "larger issues" could not be addressed until these apparent inconsistencies were resolved. Gedeon Dagan developed a robust theory of solute transport using a stochasticLagrangian framework [Dagan, 1982a, b]. In fact, an entire "school" of stochastic subsurface hydrology sprung from his work. Allen Shapiro and Vladimir Cvetkovic wrote Stochastic analysis of solute arrival time in heterogeneous porous media in 1988 [Shapiro and Cvetkovic, 1988]. The authors made an explicit assumption often implicitly made in the Lagrangian transport literature, namely that a fluid parcel deviates little from its mean trajectory in weakly 1 2 heterogeneous conductivity fields, and that the Lagrangian and Eulerian velocities are approximately equal. They consistently applied this concept, developing expressions for arrival time means and variances for nonreactive solutes in the "Standard Model Aquifer," a hypothetical aquifer characterized by a heterogeneous conductivity field with a prescribed correlation structure and subject to certain prescribed boundary conditions (see Section 2.1). These equations predicted that the mean arrival time for solute subject to a uniform resident injection is given by where angle brackets <> denote an ensemble average operator and Vh is the harmonic mean Lagrangian velocity. From the small deviation assumption, the Lagrangian and Eulerian velocities are approximately the same, thus the estimated harmonic mean Lagrangian velocity is equal to the estimated harmonic mean Eulerian velocity. However, the mean arrival time for uniform resident injection of solute after travelling several integral scales is given by <7>= x/U (1.2) where U is the arithmetic mean Eulerian velocity. This expression is also the "consistent firstorder approximation" (CFOA) of the travel time. Widespread adoption of the CFOA through out the stochasticLagrangian transport literature probably stems Feynman's requirement: we must reproduce what we already know. The harmonic mean of a positive definite process is always less than the arithmetic mean, so Equation (1.1) systematically overpredicts travel times for large displacements. Disturbing was that the estimators derived by Shapiro and Cvetkovic [1988] using clear, consistent, intuitive and logical methods seemed to violate Feynman's rule, especially in light of the success enjoyed by the Dagan school. 3 This harmonic/arithmetic mean discrepancy was explicitly noted by Dagan et al. [1992], who relegated the Shapiro and Cvetkovic [1988] result to applicability in areas "close enough to the input zone" [Dagan et al., 1992 p. 1374]. Cvetkovic et al. [1996] developed a semianalytical expression that described the transition from the "near field" harmonic mean velocity to the "far field" expression based upon the these known "endpoints" and an estimate of the Lagrangian velocity integral scale. Their detailed analysis of the Lagrangian velocity field revealed that the Lagrangian velocity is nonstationary in displacements in space, and resulted in a nonlinear propagation of the mean arrival time with distance. This was strange. Why was the Lagrangian velocity field nonstationary when the Eulerian was stationary? It was commonly assumed that the stationarity of one implied the stationarity of the other. What had started as a simple preliminary literature review resulted in a quandry. The paths followed in pursuit of this quandry led to the "research goals" of the dissertation. 1.3 Research Goals There were four broad objectives specified for this work. All centered around the flow of water and the transport of solutes in aquifers characterized by heterogeneous hydraulic conductivity fields. 1.3.1 Lagrangian Velocity Mean and Covariance A primary goal of this research was to develop a "Lagrangian covariance function" for the "Standard Model Aquifer". This was extended slightly to the mean and covariance of the spacestationary velocities along equal flowweight streamlines and timestationary velocities along equal areaweight streamlines. Quantification of this covariance function greatly simplifies the development of equations, or estimators, that describe the movement of water and solutes in the heterogeneous velocity fields associated with the "Standard Model Aquifer." 4 1.3.2 Injection Mode Analysis The stochasticLagrangian framework is convenient for the analysis of the effects of the method by which solute is introduced into the flow field, or injection mode, upon the subeqent transport. These effects of injection mode upon nonreactive solute transport are characterized by the mean and variance of a breakthrough curve, or mass arrival distribution, and by the mean and variance of mass displacement. A goal of this research was to develop estimators for these quantities. 1.3.3 Heterogeneous Reaction Parameter Set as a goal were estimators for the mean and variance of a heterogeneous "reaction parameter" integrated along trajectories of stationary velocty. These estimators were used to evaluate the coupled effects of injection mode and a heterogeneous linear equilbirium sorption process in a stochasticLagrangian framework. 1.3.4 Application to "Reality" A further goal was to develop a "traditional" estimator for the center of mass of a continuously injected solute plume subject to firstorder decay and to evaluate data from a natural attenuation experiment and numerical experiments. Estimator performance was evaluated using the results of the prevous developments. CHAPTER 2 LAGRANGIANEULERIAN FLOW FIELD RELATIONSHIP Develop an intuitive judgment and understanding for everything. Miyamoto Musashi1 2.1 Introduction Consider the steady, welldeveloped flow of water in an aquifer meeting Bear's [1972] definition of a porous medium. In most natural settings, gentle gradients induce a flow wellquantified by Darcy's law. Knowledge of the water fluid properties, relative permeability of porous medium, and potential gradients in space and time would assure a reasonable estimate corresponding water fluxes.2 One way to record this knowledge would be to create maps of these different properties in space at different times. This framework is commonly called the Eulerian framework, and is the most common approach to hydrologic problems. Another way to record this knowledge is to map properties along moving coordinates in time. This approach is commonly called the Lagrangian framework, and its application may without controversy be described as less widespread than Eulerian approaches. Dagan [1989] and Gelhar [1993] review development of flow of fluids and transport of contaminants in porous media in both of these frameworks. 1 Musashi [1982 p. 49] 2 The concepts associated with the following discussion are quite general, and are applicable to a much larger variety of environmental transport problems. The emphasis upon groundwater follows from an interest in simplicity. 5 6 Consider now flow of water in a more specific model aquifer: the steady, irrotational, and Darcian flow of water of homogeneous fluid properties in a finite rectangular twodimensional aquifer of heterogeneous hydraulic conductivity (see Figure 2.1). Mobile water 0 completely saturates the homogeneous pore space. An orthogonal y no flow boundary \rando conductivity field Cartesian coordinate system lies parallel to the aquifer boundaries, and x and y de)  ,. 0. .raConstant potential boundaries located at x = 0 and =Lmaintainaconstantaverage potential gradient J across the domain, such that water flows into the aquifer at x = 0 and out of the aquifer at i = L at a volumetric discharge rate Q (dimensions [L3T1]). For an intuitive dimensional consistency, allow the aquifer some small and 7 uniform thickness in the vertical direction, say h. The total area of the vertical plane at x = 0 through which water enters the aquifer is A. Define the spatiallyaveraged groundwater seepage velocity at this inlet plane (IP) as U fA q[a]da fA da (2.1) OA For a large class of hydrologic applications, water quantity is the sole concern. For these problems, the "route" taken by the water is of little interest, as is, say the distribution of ages or residence times. Spatial and temporal distributions of hydraulic properties hold some interest, but only in how they relate to effective bulk properties. However, a large class of practical hydrologic problems exist for which the trajectories of water parcels across the flow domain is of great importance. Many of these problems involve the fate and transport of contaminants through the aquifer system. The fundamental problem is this: to relate an Eulerian flow field to its Lagrangian, or trajectorybased, equivalent. 2.2 Flow Field Partitioning Strategies A streamline is a hydrodynamic entity that is at all points tangent to the velocity vector. For twodimensional flow, a streamtube is an object defined by two streamlines. Stream surfaces, the intersection of which form streamlines, define threedimensional streamtubes. Infinite collections of streamlines exist that subdivide the flow field into collections of streamtubes. Again, the interest of clarity suggests a focus upon twodimensional flow. For a lucid discussion of streamlines, stream functions, and stream surfaces in two and threedimensional flows, see Bear [1972]. Impervious boundaries are streamlines, thus the entire model domain is a streamtube by definition. Viewed as a "black box," Eulerian and Lagrangian distinctions have little meaning. This black box approach is the purview of reactor engineering, an important and welldeveloped discipline, from which a great deal of insight and useful mathematics may be drawn. For an excellent review, especially with respect to concepts related to transfer functions, see Jury and Roth [1990]. For a large domain, the scale at which chemical reactions takes place is generally much smaller than that of the domain. That is, processes such as sorption, microbial and chemical degradation, and the like, are dependent upon chemical potentials across distances on the order of the pore scale to perhaps that of a "microDarcy" scale, defined as the smallest volume for which Darcy's law is generally applicable. However, the scale of interest tends to be much larger, say that of a well, or its capture zone or radius of influence, or perhaps the interface across which groundwater discharges into a surface water body, or, perhaps the Floridan aquifer. In order to consider these reactions in detail, it is necessary to resolve the larger domain into smaller subdomains. For modeling real systems, the scale and degree to which the domain is subdivided should be based upon the available data. 2.2.1 Simple Subdivision Consider the water flowing into the model aquifer across the inlet, or injection, plane (IP). Assume it is of interest to divide the domain into two subdomains. There are many ways to do this. One would be to divide the domain "down the middle," that is, at a line parallel to the x axis located y = Ly/2 (curve A Figure 2.2). This division creates two subdomains of equal volume and gives an equal "volume weight" to each subdomain. The division also creates two subIPs of equal area and gives an equal "area weight" to each subdomain. For a perfectly homogeneous porous medium, the groundwater velocity vectors will be parallel to this line, and the discharges through each subdomain will be the same as well. Moreover, there will be no advective transfer of water across the boundary. However, the introduction of conductivity heterogeneity in the transverse direction will almost certainly lead to deviations of the streamlines, and to streamlines crossing the boundary. This implies that while Ly L /2 A B yo C Yo D 0 Figure 2.2: Simple subdivision of the flow field by one of four strategies: A) equalarea straight line, B) equalarea streamline, C) equalflow streamline, D) equalflow straight line. the volume of water contained in each subdomain is constant, the discharge through a plane perpendicular to the partition is a function of plane position in x. Consider now a streamline originating at (0, Ly/2) (curve B on Figure 2.2). For a homogeneous medium, this corresponds to an equal volume partitioning. That is, the domain centerline is a streamline. However, transverse heterogeneity would likely induce changes in the velocity field such that the streamline would not remain on the centerline. For general transverse heterogeneity, the flow rate entering the two domains will not be equal, in general. Thus, the flow weight and the volume weight of the two subdomains are unequal in this case. Again, the magnitude of the areas through which water enters the subdomains are equal. Consider now a streamline originating at some coordinate (0, yo) such that the volumetric flow rate is equal for both domains (curve C on Figure 2.2). This streamline corresponds to the stream function F[ Yo T [0, Ly] T [O, 0] [0, ya] = (2.2) 10 This is an equal partitioning with respect to the flow field, or the stream function. This sort of partitioning arises quite naturally in fluid mechanics and in numerical modeling of heterogeneous flow systems using streamtube approaches (see, e.g., Theile et al. [1996]). The flow weight of each subdomain is equal, although the volume weights are different. A fourth scheme would be to extend a line from (0, yo) to (Li, yo) (curve D Figure 2.2). In this case, the subdomains share only equal flow rates at x = 0. In subsequent discussions, the phrase Eulerian trajectory will refer to such a straight line trajectory, and the phrase Lagrangian trajectory will refer to trajectory of a streamline. The phrase area weight will refer to an equal area partitioning scheme at a definition plane, and the phrase flow weight will refer to an equal discharge partitioning scheme. At this point, it is important to note that no mention has been made of boundary conditions as they relate to transport equations. These partitioning schemes are simply ways of subdividing a flow field. An area weight Eulerian trajectory scheme corresponds to traditional spectral approaches that assume small deviations about a mean trajectory. The use of trajectory is quite general, and does not imply a "natural" path or streamline. 2.2.2 General Subdivision Consider dividing the flow field into N partitions according to the strategies described (see Figure 2.3). Theile et al. [1996] demonstrated that streamlines may represent streamtubes. Therefore, the remainder of this discussion will focus upon streamline properties, with the implication that these properties may be extended to streamtubes as well. A primary purpose of subdividing the domain is to represent the continuum of domain properties with a discrete set of data. A regular finitedifference grid may be thought of as a set of points at equal spacing in space along a collection of areaweighted Eulerian trajectories. An alternative method would be to create spacings based upon travel times along the trajectory, where travel time is defined as .........~ ~ ~ ~ ~ ~ ~~. ......   ..... .A D B C Figure 2.3: Four different partitioning strategies for the same flow field: A) area weight Eulerian trajectory, B) area weight Lagrangian trajectory, C) flow weight  Lagrangian trajectory, D) flow weight Eulerian trajectory. 12 the path integral of the inverse velocity over the trajectory r ol] = ds (2.3) where is a unit vector that is at all points the parallel to the trajectory. It is clear that this sort of Lagrangian description has little meaning in a static system, and is in this sense less general than Eulerian descriptions. Thus, Lagrangian description is a tool reserved for dynamic systems and is a supplement to Eulerian descriptions. It is worthy to note the distinction between a description, a general approach, and a trajectory, a specific entity used in a description. To summarize, a primary objective is to approximate a steady and irrotational flow field continuum with a discrete set of trajectories. The trajectories follow either a straight line (Eulerian trajectory) or a streamline (Lagrangian trajectory). The trajectories carry either an equal flowweight or an equal injection areaweight. 2.3 Descriptions of Fields Consider once again the model aquifer depicted in Figure 2.1. As mentioned in the introduction of this chapter, the aquifer has certain features in which a hydrologist might be interested, e.g., the intrinsic permeability or the porosity. These features may be thought of as a field in which parameters vary continuously. This is a fundamental concept of continuum mechanics and is discussed in a hydrologic context by Bear [1972], Dagan [1989], and Gelhar [1993]. The ultimate description of any such field would be a map of parameter values valid at any point, time, or scale. In lieu of this "ultimate description," stochastic hydrologists characterize spatial heterogeneity with conditional probability density functions that reproduce the dominant characteristics of the field and observed values at given locations. The entire subject of geostatistics is devoted to this end (see Journel and Huijbregts [1978]). An Eulerian field refers to areaweighted Euleriantrajectory oriented field. A Lagrangian field refers to a flowweighted Lagrangiantrajectory oriented field. 13 It is instructive to illustrate these concepts with an example and to interpret the results of this illustration. 2.4 Numerical Experiments The illustration of these concepts requires a heterogeneous flow field and collections of trajectories, both Eulerian and Lagrangian, that correspond to this flow field. The requirement that the flow field be based upon "reference fields" of known properties enhances the possibility of theoretical interpretation drawn from the large body of literature concerning flow and solute transport in idealized porous media. A further requirement that our reference fields bear some resemblance to what might be expected in nature enhances the possibility of useful applicability. A numerical simulation code modeled steady groundwater flow in a hypothetical heterogeneous twodimensional aquifer. A particle tracking code traced Eulerian and Lagrangian trajectories across this aquifer, recording data that included the local groundwater velocity and local hydraulic conductivity at intervals in space and time along the different trajectories. Statistics, specifically mean, variance, and covariance, summarize the results of these "measurements." The conceptual view of the aquifer is that given by Figure 2.1. To generalize the system, correlation length, or integral scale, of the log conductivity field Aln k serves as a characteristic length with which to normalize length scales. The characteristic time Alnk/U serves to normalize time scales, where U is the arithmetic mean velocity defined in Equation (2.1). The definition of an integral scale is the integral from zero to infinity of the correlation function R, or, A = R[s]ds (2.4) 0O 14 The design criteria for the model aquifer are as follows In k correlation length Aln k = 1 effective conductivity Ke = 1 average gradient J = 1/50 mobile water content 0 = 1/5 From these criteria and Darcy's law, the areaaveraged velocity is U = KeJ/O = 1/10. Clearly, the parameter values are artificial and selected for convenience, but values for Ke in md1 and 0 fall in ranges typical for sand [Domenico and Schwartz, 1990]. Two 5000 replicate Monte Carlo simulations of aquifers with log conductivity variances of 1/2 and 1, respectively, comprised the suite of aquifer simulations. The following sections describe the experimental design and the results. 2.4.1 The Log Conductivity Field The foundation upon which to build the illustration is a heterogenous hydraulic conductivity field. A commonlyused stationary lognormally distributed and exponentially correlated isotropic model serves as the Eulerian hydraulic conductivity reference field (see e.g., Freeze [1975] or Gelhar and Axness [1983]). Working from a hypothesized conductivity field is preferable to generating a hypothetical velocity field with specified characteristics (e.g., Bellin et al. [1994]) because an objective of this work is to describe the underlying conductivity field along Lagrangian trajectories and relate its Lagrangian description to its Eulerian properties. For completeness, it is pertinent to briefly review description of a heterogeneous field by treating the process as being "random" Measured values of hydraulic conductivity of virtually any natural porous formation will vary from place to place. In the absence of measurement error, this variability in the measured values reflects the heterogeneity of myriad physical factors and processes that control the hydraulic conductivity. As mentioned, measurements of hydraulic conductivity in the field often follow a lognormal distribution. In these cases, it is useful to conceive of the 15 conductivity field as a realization of a lognormal random process. The process is said to be stationary if the ensemble mean and variance do not vary in space. The process is nonstationary if some trend exists in the mean or variance. In practice, only single realization is available, and the "ensemble" properties are ultimately unknowable. However, if the field is large compared to the scale of variability, one may assume ergodicity, or, that the ensemble parameters may be estimated from a sample drawn from the single realization since there exists a the possibility of a large number of statisticallyindependent observation points. For a more detailed development of these concepts from the hydrologist's perspective, see Gelhar [1993]. For a more formal development of these concepts, see Papoulis [1991]. The model of hydraulic conductivity variability is a stationary, exponentiallycorrelated isotropic lognormal distribution. A lognormal random variable is one for which its natural logarithm is normally distributed. As such, three parameters summarize the statistical properties of the adopted conductivity model, namely the mean of log k Iln k, the variance of log k U2n k, and the correlation length of log k Aln k. The covariance model is given by Clnk[r] = nkexp[r/Alnk] (2.5) where r is the separation between two points in the log conductivity field. Use of this simple covariance model requires stationarity of variance, and usually implies a stationarity of the mean as well. The moments of the "real space" conductivity field are given by In the absence of a perfectly described physical aquifer, the method of turningbands generated equallylikely realizations of a hypothetical conductivity field with 16 prescribed ensemble statistics [Tompson et al., 1989]. Turningbands was capable of producing two and threedimensional fields with the desired statistical properties and was available as a portable and easily modified code. Figure 2.4 illustrates the lognormality of the a typical twodimensional log conductivity field realization. Figure 2.5 compares an estimate of the longitudinal covariance to the hypothesized model. In order to capture the features of the log conductivity spatial variability, the log conductivity standard normal 0.8 S0.6 7. 5 0.4 E 0.2 3 2 1 0 1 2 3 log conductivity Figure 2.4: Log conductivity distribution from a 41400 node realization generated by the turningbands method compared to the standard normal distribution. Target distribution parameters are 1lnk = 0 and r2nk = 1. turningbands algorithm generated five log conductivity observations per correlation length. This level of discretization falls in the range of other similar sets of simulations. Cvetkovic et al. [1998] used a discretization of 4 nodes per Alnk for simulations of fields with oank as large as 1. Bellin et al. [1992] used a discretization of 8 nodes per Alnk for simulations of fields with a1,nk as large as 1.6. 17 estimated In k covariance + exponential covariance 0.8 a 0.6 C 0 o 0.4 0.2 ' 4   ta 0 0 1 2 3 r/AInk Figure 2.5: Estimated log conductivity longitudinal covariance estimated from a realization generated by the turningbands method compared to an exponential covari2 I ance function. Variance of log conductivity Calnk = 1 and correlation length AInk = 1. Spacing for log conductivity values was Alnk/5. Trajectories covered several characteristic lengths Al, k and several characteristic times Aln k/U to capture both nearfield and farfield behavior in space and Lagrangian time. The first requirement was easy to enforce by making the domain size equal to or greater than the required number of correlation lengths. The second requirement was somewhat more problematic, since the distance traveled in some number of Alnk/U varies from trajectory to trajectory. Preliminary simulations indicated that a domain 40Alnk long was adequate to capture travel up to 10 characteristic time scales Alnk/U in 2000 trajectory realizations generated in the most heterogeneous field considered (a2nk = 1). Thus, the selected domain extent was 40Ank The conceptual model aquifer is a large, bounded, primarily twodimensional, steady flow field. As such, it should exhibit a variety of head, head gradient, and velocity nonstationarities due to the boundary conditions (see Figure 2.1). While 18 researchers such as Rubin and Dagan [1988],[1989], and Osnes [1998] have investigated such nonstationarities, the additional complication of a rigorous consideration of these sorts of nonstationarities would unduly obfuscate the very simple concepts illustrated here. To minimize the effect of these nonstationarities, "solute transport" is constrained to take place some distance from the boundaries of the numerical domain. Following the example of Bellin et al. [1992], a 3Alnk "buffer" was established on both sides of the required 40 for a total longitudinal domain length of 46Alnk. Preliminary simulations indicated that the maximum transverse displacement of a streamline in a an k = 1 field was around one third the longitudinal displacement. For a streamline originating on the aquifer centerline parallel to the x axis, an aquifer width of 36Aln k was deemed adequate to contain large transverse displacements over a 40Alnk longitudinal displacement and provide a similar "buffer" to mitigate boundary effects. Combining the "physical" domain requirements with the discretization required to capture the log conductivity variability resulted in a numerical grid 230 by 180 for a twodimensional aquifer and a numerical domain with 41400 nodes. A threedimensional analysis under the same assumptions would have required a numerical domain of more than 7 million nodes. Simulations of this magnitude were untenable in light of the available computational resources. 2.4.2 Velocity Field Generation A mixed finiteelement scheme solved coupled velocity and pressure equations derived from Darcy's law and the continuity equation. The code took input parameters that matched the design criteria and the log conductivity field generated by turningbands. Andrew I. James wrote and tested the core of the code. James et al. [1997] used a derivative version of this code to generate flow fields for synthetic tracer tests. A mixed finiteelement solution to the flow field was specified, following the results and arguments of Mos6 et al. [1994]. The mixed finiteelement scheme is known 19 to provide an accurate velocity field solution for groundwater flow in heterogeneous conductivity fields, one that is more accurate than other, more traditional, schemes [Mos6 et al., 1994]. 2.4.3 Trajectory Generation An objective of this research was to record data, such as local velocity and and conductivity, along different collections of trajectories, both Eulerian and Lagrangian. These trajectories were traced using a "particletracking" scheme based upon a semianalytical method presented by Pollock [1988]. This method assumes that velocity varies linearly over an element, and that the velocity in some direction, say i, is independent of displacements in orthogonal directions. From these assumptions, it is a relatively simple exercise to determine the displacement of a "particle" in a linear velocity field for a given time, or the time required to make a certain displacement. A computer code recorded particle position, time, total path length along trajectory, velocity, local conductivity, and a value from another random field with the same statistical properties as, but uncorrelated to, the log conductivity field at several control planes equallyspaced along the mean flow direction and at several "control times" after injection along both trajectories. This uncorrelated random variable was recorded for reactive transport simulations, to be discussed in subseqent chapters. 2.4.4 Monte Carlo Simulation The possibility of theoretical interpretation of the results is enhanced by requiring nearergodic samples from the population. That is, the samples must reproduce the population statistics. To do this, the number of statistically independent replicates drawn from the velocity and conductivity fields must be quite large. Rather than attempt a massive singlerealization, the following method generates a large number of statistically independent realizations 20 generate a large, twodimensional conductivity field with the desired statistical properties generate a flow field by imposing a constant mean gradient across the saturated conductivity field trace the trajectories of two particles injected at a common "control location", one trajectory is a straight line or Eulerian trajectory and the other a streamline or Lagrangian trajectory, and record properties at equal lags in space and in time along trajectories that pass through "inner domains" that are relatively undisturbed by nonstationarities induced by the boundaries At this juncture, it is perhaps worth mentioning a minor "philosophical" point. These simulations do not represent a real aquifer, per se, but rather a hypothetical aquifer with known properties. The goal is understand the convolution of the physics of the groundwater flow and a known heterogeneity that has analogs in nature, and it is hoped that the reader finds some merit to this modest goal. The 41400 conductivity values needed for any single realization do not represent 41400 measurements taken from some real aquifer. The Monte Carlo aspect of this experiment reenforces this point. 2.4.5 Sampling Strategy Consider "dropping particles" into a flow field, with no particular regard to location. These particles move along streamlines without interaction with the porous medium. Any one location into which one is dropped is as equallylikely as any other. The trajectories associated with these particles carry an equal area weight, analogous to equalarea subdivision of an injection plane. Consider now, the same disposal of particles, but this time an area is assigned to each particle at the injection point such that an equal discharge is associated with each particle, and its associated trajectory. For a homogeneous mobile water content, this area is inversely proportional to the local velocity normal to the area. Thus, each particle may be thought to carry an equal 21 flow weight. In this way, it is possible to "map" from one collection of trajectories to another. At a given control plane or time, the areaweighted mean and variance of some property over N realizations, say v was calculated using the familiar estimators3 S N i=1 N var[va] = N N ) (2.8) i= 1 While these particular estimators seem obvious, it is important to realize that "dropping the particle" into the flow field assigns an equal area weight to each trajectory. For a single realization, this is equivalent to starting each particle with equidistant spacing. These estimators are based upon an assumption of equallylikely statisticallyindependent samples (see, e.g., Mood et al. [1974]). In order to draw even a few hundred independent samples from a single realization, the required flow domain would be enormous. This was a primary consideration for adopting Monte Carlo simulation using a minimal domain. The flowweighted mean property v was calculated via v N= p=u[O] (2.9) Ei=1 Ui[O] 3 Press et al. [1988] have criticized the variance estimator in Equation (2.8) as being prone to accumulation of roundoff error for large data sets. The rather substantial size of the data sets are naturally a concern. Comparision of the variance calculated by a corrected twopass method presented by Press et al. [1988] and the "traditional" estimator yielded no difference between the two methods within the first five significant digits. 22 where ui[0] is the initial velocity recorded for trajectory i. The variance by var[vf = i=1 ui[O](v Vf)2 (2.10) z,1 U[0] An alternative method for generating a collection of streamlines would be to identify the streamline corresponding to the "center" stream function expressed in Equation (2.2). In this case, the estimators of the equal flow statistics would be Equations (2.7) and (2.8). The equal area statistics would be given by a = ] (2.11) 1i= 1/ui[] and var[va] Ei=(Vi a)2/ (2.12) =~= 1/ui[O] The "reference" collection of the trajectories selected uniformly in space was specified in the simulations for convenience. 2.4.6 Simulation Results The following sections present the results of two 5000 replicate Monte Carlo simulations. Detailed discussion of these results is deferred until after the development of some theoretical models with which to analyze these results. However, some heuristics are provided to familiarize the reader with some of the concepts developed in subsequent sections and chapters. A summary of the notation used with the following results provides clarity. AW denotes an areaweighted trajectory and FW denotes a flowweighted trajectory ET denotes an Eulerian trajectory and LT denotes a Lagrangian. A particular weighttrajectory combination is concatenated the with a /. Thus, an areaweighted Eulerian trajectory is denoted as AW/LT. For a given simulation, that is, for one set of 5000 replicates for a log conductivity field of some given arlk,, the statistics of the log conductivity and velocity along the 23 four different trajectories are calculated at intervals in time and in space. Thus, for any one property, say log conductivity, and any one statistic, say the mean, there will be two sets of results, namely that for a2nk = 1/2 and that for o2nk = 1. Associated with each set will be two subsets, namely statistics along the four trajectories recorded at space or time intervals. For convenience, the notation and the trajectories to which each notation corresponds as depicted in Figures 2.2 and 2.3 is summarized here: AW/ET areaweighted Eulerian trajectory (Figures 2.2 and 2.3 A) AW/LT areaweighted Lagrangian trajectory (Figures 2.2 and 2.3 B) FW/LT flowweighted Lagrangian trajectory (Figures 2.2 and 2.3 C) FW/ET flowweighted Eulerian trajectory (Figures 2.2 and 2.3 D) 2.4.7 Log Conductivity Statistics The log conductivity field serves as the "foundation" upon which the simulations are built, as well as the theory. The data for the areaweighted Eulerian trajectories correspond to those design criteria specified. The log conductivity mean and variance for AW/ET appear stationary in space (see Figures 2.6 and 2.7). That is, they appear to "bounce around" their "target" values with no apparent trend. Moreover, the statistics of FW/LT appear likewise stationary, but the mean is somewhat larger and the variance slightly smaller than AW/ET. The magnitude of the differences appear proportional to aOnk. A comparison of the "target" values, e.g., the specified inputs into the turningbands algorithm, and the observed ensemble values along the trajectories is presented in Table 2.1. The "observed" entry in Table 2.1 represents the statistic average along the trajectory, as estimated by regressing a zeroslope line through the data using a generalized leastsquares method. There is close correspondence between the target and observed values, supporting an ergodic hypothesis for the log conductivity statistics. That is, sample statistics are an accurately estimation 24 0.6 0.5 FW/LT + AW/LT x FW/ET 0.4 AW/ET n Eq. 2.72 reference  0.3 *  5 r~ ~  fx XE 0.1 0.1 a 0.2 0 10 20 30 40 x/Al k 0.6 0.5     FW/LT + + + AW/LT x + + N < FW/ET 0.4 + + + AW/ET o xx x Eq. 2.72 xx reference 0.3 x S0.2 x 0.1 00 0.1 x b 0.2 0 10 20 30 40 X/Aln k Figure 2.6: Mean log conductivity as function of displacement in mean flow direction for the four different trajectory collections and a1k = 1/2 (a) and oa = 1 (b). 25 1.1 1.05 ae iFW/LT + o 5 m AW/LT x 0.95 +0 + + o x FW/ET + + +x + X 4 X X +0 C4 0.9 + + x+ !  ........................................4.......... 0.85 a 0.8 0.75 0.7 0 10 20 30 40 X/AIn k 1.1 1.05 FW/LT + AW/LT x S ii IN WM[ 0 N FW/ET 1 .......  AW/ET Sx9 0 reference ... 95 xx + + x Eq. 2.73 0.95 : xx x x x + + + X++ x + + xCN 0.9 + b + 0.85 b 0.8 0.750.7 0 10 20 30 40 X/AXlnk Figure 2.7: Log conductivity variance as function of displacement in mean flow direction for the four different trajectory collections and onk = 1/2 (a) and a2 = 1 (b). Variances are normalized by target variance of Eulerian log conductivity field an k" 26 Table 2.1: Comparison of target log conductivity values and observed values along areaweighted Eulerian trajectories. parameter target observed IAlnk('12nk = 1/2) 0 0.003 nk( 1nfk = 1) 0 0.005 12nk 1/2 0.498 an k 1 1.002 of the population statistics. Again, this is important for a quantitative theoretical analysis of these results. The AW/LT and the FW/ET exhibit distinctly nonstationary behavior in the mean and less distinct behavior in variance as a function of displacement along the mean flow direction (see Figures 2.6 and 2.7). At the injection point, the AW/ET and AW/LT statistics are identical and the FW/ET and FW/LT statistics are identical, because these trajectories carry the same weight and the points in consideration are identical (see e.g., Figure 2.2). At very large displacements, one would expect the statistics of like trajectories to assume similar characteristics, regardless of weight. As the trajectory traverses several integral scales that characterize the heterogeneity in question, "information" about its starting location, is diminished. In the case of the FW/ET, the statistics "start" at the same place as FW/LT, and eventually "end" at the same place as the AW/ET. Similarly, the AW/LT statistics start where do the AW/ET and end where do the FW/LT. In fact, the data exhibit this behavior (see Figures 2.6 and 2.7). The log conductivity statistics in time show somewhat different behavior. Only the AW/LT appear stationary in time (see Figures 2.8 and 2.9). There is an important observation to be gleaned from these results. The AW/LT statistics, in addition to appearing stationary, also appear to have a similar value to the input AW/ET statistics (see Table 2.2). Recall, the AW/LT statistics are clearly nonstationary in space, and are clearly distinct from those of the AW/ET. This observation greatly 27 0.4 FW/LT + AW/LT x FW/ET x AW/ET o est AW/LT 0.2 w++ A x et A +o 0 xx 0.2 0 DD E 0.2 0.4 a 0 2 4 6 8 10 Ut/I 0.4 + FW/LT + + AW/LT x ++ FW/ET w + AW/ET o + + est AW/LT an 0.2 + '5 ++++ E+ 0 Cx 0 m E 0.2 0.4 b 0 00 0 2 4 6 8 10 Ut/I Figure 2.8: Mean log conductivity as function of travel time to control planes along mean flow direction for the four different trajectory collections and ank = 1/2 (a) and ~nk : 1 (b). 28 1.1 I 0 1.05 FW/LT + + x + x O AW/LT x x Oxx+ x ox FW/ET K xx xa+a + X~ AW/ET o .......... +.. .... .s ........ reference  X X )K A TX X XOn a + a x xx a K X X + ++X +00 0 + X X X + S4E 0.95 + 0 ++ xx N++ x xx+ N0 0 + + + N NN NN N ++ 0.9 N + 0.85 a 0.8 0 2 4 6 8 10 tU/Aln k 1.1 1.05 x O FW/LT + x n AW/LT x x x 0 x o o FW/ET X X X 0 XX X+ XN xxx x x x Nxf~o' AW/ET o 1 x o xx    xx reference 0 xx x + N + x~ +xX X Ox ++NX X x + N ++ + + xx95 ++ x+ + x+ xxx + 0.9 + + +' 00+ 0 + 0 + N +0 + ++ N NN 0.85 K NI b 0.8 0 2 4 6 8 10 tU/Alnk Figure 2.9: Log conductivity variance as function of travel time to control planes along mean flow direction for the four different trajectory collections and arnk = 1/2 (a) and an k = 1 (b). Variances are normalized by target variance of Eulerian log conductivity field cr .lnk 29 Table 2.2: Comparison of target log conductivity values for areaweighted Eulerian trajectories in space and observed areaweighted Lagrangian trajectories in time. parameter target observed I.nk( 1nk = 1/2) 0 0.000 Inlnk(12nk = 1) 0 0.017 a1nk 1/2 0.497 o2 k 1 1.000 simplifies relating the Eulerian and Lagrangian fields, as shall be shown. The "longtime" asymptotic mean log conductivity values for both Eulerian and Lagrangian trajectories are somewhat less than their "largedisplacement" counterparts, and the magnitude of the difference is proportional to oalk (compare Figures 2.6 and 2.8). This is consistent with the harmonic averaging implicit with a timeaveraged timedependent velocity and the arithmetic averaging implicit with a spaceaveraged spacedependent velocity. Recall, we map the quantities onto the trajectories by following velocities, whether along "natural" streamlines (Lagrangian trajectories) or not (Eulerian trajectories). A simple heuristic is as follows. Low conductivities correspond to low velocities. For an equaltime spacing (or equallylikely observations in time) more "measurements" will come from low velocities, since the particle tracing the trajectory spends more time going slow. For an equaldistance spacing more measurements will come from high velocities, since the particle "covers more ground" while going fast. This is perhaps more intuitive when discussed in terms of the velocity field itself. 2.4.8 Velocity Statistics The mean velocities along the different trajectories exhibit a similar behavior as do the mean log conductivities along the same trajectories (see Figure 2.10 and compare Figure 2.6). The AW/ET exhibits apparently stationary behavior in the mean and variance (see Figure 2.11) in accordance to established theory (e.g., Dagan [1989]). The AW/LT exhibits nonstationary behavior in the mean as observed by Cvetkovic et al. [1996]. Again, the AW/ET and FW/LT trajectories appear to be stationary 30 1.45 1.4 1.35 a FW/LT + AW/LT x FW/ET 1.3 A AW/ET o Eq. 2.35 .1.25 Eq. 2.43 Eq. 2.44 1.2++ + 1 + reference . 1. + x + + xx 1.15 xxX x x 1.1 1.05 0.95 0 10 20 30 40 /Aln k 1.45 + 1 1.4 1.35 xx FW/LT + i x,, AW/LT x 1.3 FW/ET x AW/ET o x, Eq. 2.35  1.25 Eq. 2.43  Eq. 2.44 1.2 b reference 1.15 : 1.1 1.05 0 0.95 0 10 20 30 40 z/AI k Figure 2.10: Comparison of positiondependent mean velocity estimators to data for the four collections of trajectories for Uank = 1/2 (a) and Uank = 1 (b). Velocities are normalized by U. 31 2 1.8 a FW/LT + AW/LT x FW/ET x AW/ET a 1.6 Eq. 2.37 + + Eq. 2.55 + + Eq. 2.57 + xx+x + Eq. 2.13 L 1.4 ++ .....x..x + x x + 2 S+ + + + x S ." x x x +* S1.2 xx+ x S" + * ... X. x X x TL 1.8 + x* x + FW/LT + + ++ x x Cx +x x AW/LT x x xx + x + x+ FW/ET x x x +x AW/ET S1.6 x xx + Eq. 2.37 Eq. 2.55 S Eq. 2.57 L 1 b Eq. 2.13 ........... 1.4 b S1.24 X 1o du 0.8 0 10 20 30 40 /A, lnk Figure 2.11: Comparison of the positiondependent velocity variance estimators to data for the four collections of trajectories for onk = 1/2 (a) and U nk = 1 (b). The variances are normalized by U2bn k Ink' + 32 in both the mean and variance. The velocity variances exhibit behavior to similar to that of the means for their corresponding trajectories and quite different than the behavior exhibited by the log conductivity variances (compare Figures 2.11 and 2.10 and compare Figures 2.11 and 2.7). There is a qualitative similarity between the nonstationary velocity means and the nonstationary velocity variances. A heuristic is that the velocity statistic "starts" at one value "in the nearfield" and transitions to an asymptotic value "in the farfield." The following sections demonstrate that these endpoints and the transition are predictable. The qualitative dissimilarity between the log conductivity variance and velocity variance is an artifact of examining the log conductivity values on the one hand and the "real space" velocity values on the other. Consider the FW/LT log conductivity mean and variance for the ank = 1 simulation. The FW/LT mean log conductivity is larger than that of the AW/ET, but the variances are roughly the same value. In fact, regressing a zeroslope line through the data results in a "mean" FW/LT variance that is slightly less than that of the AW/ET. However, the "real space" conductivity mean and variance are both functions of both the mean and variance of the log conductivity process (recall Equation (2.6)). A cursory examination of the data in conjunction with Equation (2.6) reveal that the mean and variance of the "real space" conductivity along FW/LT are both higher than their AW/ET counterparts, and the magnitude of this difference is proportional to Unk The mean velocities in Figure 2.10 are normalized by the arithmetic mean velocity U predicted from the input parameters for the flow simulator. That the mean velocity along AW/ET closely follows 1 is a good indication that the effective conductivity expression is accurate. Notice that the FW/LT and farfield AW/LT values are higher than AW/ET, and the magnitude is proportional to the variance of the log conductivity (see Figure 2.10). This behavior is intuitively correct. Consider a 33 particle injected into a flow field. Constraint of the trajectory to an Eulerian trajectory forces the trajectory through lowvelocity areas that a natural streamline, or Lagrangian trajectory, would have otherwise bypassed. Of course, some streamlines pass through even the lowest velocity areas, but these are few in an equalflow sense. Thus, equalarea weighting preferentially weights low velocity zones. That there are implications for sampling schemes in heterogeneous media, especially in context of multilevel samplers, should be intuitively obvious to even the most casual observer. A discussion of some of these implications is made in a subsequent section, after these concepts are developed further in a more quantitative and theoretical framework. It is appropriate to evaluate the adopted expression for the AW/ET longitudinal covariance function suggested by Dagan and Cvetkovic [1993], namely Cua[r] = O nkb[e]U2 exp [b[e]r/Alnk] (2.13) where r is the separation between points along the ET, and b[e] is a function of the anisotropy ratio e. The anisotropy ratio itself is a function of the correlation lengths of the log conductivity in the different directions. The primarily consideration is that b may be derived from theoretical considerations and is not some ad hoc fitting parameter. The value of this function for a twodimensional isotropic aquifer is 3/8. Again, for further background on this function, see Dagan and Cvetkovic [1993]. In the limit r = 0, the covariance reduces to the point variance, which for the model case is o~ = 3U2In k/8. For the remainder of this discussion, a reference to b contains its dependence upon e implicitly. The subscript ua may be decoded as follows: for a subscript tw, the t is the trajectory and property (u for Eulerian velocity, and v for Lagrangian velocity) and the w is the weight (a for area and f for flow). Thus, C,, is the covariance function for velocities along Eulerian trajectories. We present this covariance function and that covariance function estimated from the assumedstationary collection of AW/ETs in Figure 2.12 The perceptive observer will note 34 AW/ET + FW/LT x Eq. 2.13 . 0.8 c 0.6 0 0 8 0.4 X 0.2 + X ' + . 0 5 10 15 r/A lnk Figure 2.12: Comparison of spatial Eulerian (AW/ET) and Lagrangian trajectory (FW/LT) velocity covariance data to covariance model. Data are taken from the onk = 1 set of realizations. Covariances presented as dimensionless lag correlation functions. that the analytical expression overestimates the covariance at short lags and underestimates the covariance at longer lags. This discrepancy and its implications are discussed in some detail in a subsequent section. The velocity variances in Figure 2.11 are normalized by the zerolag covariance predicted by placing the input simulation parameters into Equation 2.13. The AW/ET velocity variance appears wellpredicted by this expression, as evidenced by a relatively good fit to 1 for both values of onk in Figure 2.11. A regression of a zeroslope line through these data indicates that the a priori estimate under predicts the observed values by about five percent for the larger a2nk (see Table 2.3) The apparent stationarity of the velocities along the FW/LTs suggest estimation of the covariance function for these velocities as for the AW/ET (see Figure 2.12). The correlation length of the FW/LT is very close to that of the AW/ET, although the actual value is probably slightly higher. Recall the definition of the correlation length as given by 35 Table 2.3: Comparison of a priori Eulerian velocity variance to observed values. The decimal representation of the exact fractional values of the predicted variance are provided for convenience. parameter predicted observed (G nk = 1/2) 3/1600=0.001875 0.001871 O'u(O nk = 1) 3/800=0.00375 0.00394 Equation (2.4). The correlation function at each considered lag along the estimated functions is systematically larger for FW/LT, even though the magnitude of the difference is small. Integrating these functions would result in a slightly larger value for FW/LT, and thus a larger correlation length. In the case there is no variability in the log conductivity, Eulerian and Lagrangian trajectories coincide. Moreover, an equalarea weight corresponds to an equal flow weight. Thus, the estimated correlation functions from "either set" would be identical. From this conclude that there is some dependence of the correlation length of the FW/LT trajectories upon Unk' Cvetkovic et al. [1996] observed similar behavior in log velocity correlation functions estimated from AW/LTs with a stationarity assumption for the nonstationary AW/LT velocities. Thus, this is an slight extention of this work by identifying the stationary streamlinebased velocity field from which to work, and work from that field. Moreover, this supports their findings of an increasing correlation length with aInk Understanding the relationship between the variability of the log conductivity field and the correlation length of the associated velocity field is of fundamental importance to completely understanding the nature of flow in heterogeneous media. Sadly, these "secrets" remain uncovered and must be left this to future research efforts. The velocity means as a function of time exhibit qualitatively similar behavior as do the log conductivity means in time (see Figure 2.13 and compare Figure 2.8). Again, only the AW/LT exhibits stationary behavior in time for both the mean and variance (see Figure 2.14) as did the AW/LT log conductivity values. The mean and variance of the AW/LT in time are approximately those of the AW/ET in space 36 1.4 'FW/LT + AW/LT x FW/ET = 1.3 AW/ET o Eq. 2.63 Eq. 2.59 1.2 Eq. 2.67 Eq. 2.66 tU/Alnk 1.1 ++FW/LT ++, AW/ET x Eq. 2.63 Eq. 2.590.9 Aa: ... ... 0.9 a 0.8 0.7 anmamam 0 2 4 6 8 10 tU/Aln k 1.4 FW/LT + AW/LT x 1.3 FW/ET AW/ET o S,,Eq. 2.63 'K + Eq. 2.59 1.2 Eq. 2.67 .......... Eq. 2.66 ) + ++ 1.1 X +++;+XX X xxxxxxxX XxxXx 0.9 o " 0 8o ... .... " 0.7 0 2 4 6 8 10 tU/Al k Figure 2.13: Comparison of timedependent mean velocity estimators to data for the four collections of trajectories for u~nk= 1/2 (a) and a"nk = 1 (b). Velocities are normalized by U. 37 FW/LT + AW/LT x 1.8 FW/ET AW/ET o Eq. 2.65 .1.6 Eq. 2.61 Eq. 2.70 N Eq. 2.69 Eq 269 4 1.4 1.2 + ++ ) k+ +. + + x X Xxx+ +.x + x C 4 P 1 ..   ..  , 0.8 .. .. 0.6 a 0.4 0 2 4 6 8 10 tU/An k FW/LT + AW/LT x 1.8 1 FW/ET w AW/ET o Eq. 2.65 1.6 Eq. 2.61 . Eq. 2.70 "14 1 .4 Eq. 2.69 . a 1 + ++ 1.2 + ++ + x x x + xx N 1  0.8 0.8 a b 0.4 0 2 4 6 8 10 tU/Aln k Figure 2.14: Comparison of timedependent velocity variance estimators to data for the four collections of trajectories for onk = 1/2 (a) and oink = 1 (b). The variances are normalized by U2b n k 38 for both simulations. This result is important, because it implies the possibility of mapping from AW/ET to AW/LT. We estimated the covariance function for the stationary AW/LT velocities (see Figure 2.15). As an estimate of the covariance AW/LT Eq. 2.61 0.8 r 0.6 0D+ ++ 8 0.4 0.2+ ++++ 0I 0 2 4 6 8 10 tU/Aln k Figure 2.15: Comparison of temporal Lagrangian trajectory (AW/LT) velocity covariance data to covariance model. Data are taken from the an k = 1 set of realizations. Covariances presented as dimensionless lag correlation functions. function, the substitution r = Ut was made in Equation 2.13 and the result compared to the estimated covariance function (see Figure 2.15). The closedform predictor behaves similarly to its spatial counterpart: earlytime correlation is underpredicted and longtime correlation is overpredicted (compare Figures 2.15 and 2.12). For completeness, one other aspect of the velocity field was examined, namely the distribution of initial velocities (see Figure 2.16). The logarithm of the initial velocity was rankordered and subject to the transformation In v In v In v' lnv (2.14) In v 39 0.8 In V standard normal . 0.6 0.4 0.2 0 " 3 2 1 0 1 2 3 log velocity Figure 2.16: Normalized log initial velocities for twodimensional model aquifer compared to the standard normal distribution. The random variable was normalized by the sample mean and standard deviation to illustrate the lognormality. where Inv is the sample mean and in,, is the sample standard deviation. Working here with the logarithm, rather than the velocity itself, better demonstrates the apparent lognormality of the velocity fields. Recall, the logarithm of a lognormal random variable is itself normally distributed. Figure 2.16 indicates that the transformed log initial velocities are welldescribed by a standard normal distribution. The lognormality of the velocities seems logical since the conductivity field is itself lognormally distributed. However, the velocity is in general a product of random variables: the conductivity, the head gradient, and the mobile water content. The logarithm of the velocity is the sum of random variables. Sums of random variables tend toward normality by the central limit theorem (see Mood et al. [1974]), and thus, products of random variables tend towards lognormality. It would be interesting to test this 40 hypothesis by examining the velocity characteristics of flow fields arising from different log conductivity distributions (e.g., fractal, uniform, white noise, etc.), and this is left to future work. 2.4.9 Head and Head Gradient Head and head gradient information along the trajectories was not recorded. This is a regrettable oversight, since much of the behavior of the flow water and the transport of solutes in heterogeneous media can be summarized as the interplay of head, head gradient and conductivity covariance functions. The literature concerning Eulerian head and head gradient fields is extensive and is summarized in large part by Gelhar [1993]. In fact, much is drawn from this body of literature for some of the theoretical developments in the following sections. That this information was not recorded detracts mainly from the completeness of this work, and further investigations are left to future works. 2.5 Theory 2.5.1 Effective Conductivity As with the simulations, the theoretical developments are founded on the bedrock of the log conductivity field. Section 2.4.7 demonstrates that the simulation results well match the "target" ensemble quantities. From this field is had directly a derived parameter, namely the a priori effective conductivity Ke. An expression for the effective conductivity, similar to that given by Gelhar and Axzness [1983], is derived here to illustrate the mechanics of different techniques employed through out these derivations. Consider a onedimensional Darcy's law of the form q = k d(2.15) dx where q is the specific discharge, k is the hydraulic conductivity, and is the hydraulic head. Consider a large stationary conductivity field such that the log conductivity 41 may be resolved into a mean, and a zeromean perturbation, e.g., In k = l.ink + f. Let the head gradient likewise be resolved into a mean and zeromean perturbation, e.g., d /dx = J(1 + j). Inserting these expressions into Equation (2.15) yields q = exp[Allnk + f]J(1 + j) (2.16) = KgJ exp[f](1 + j) Expand the exponential term in this equation as a Taylor series about zero, truncating the series at secondorder, and expand the product of the perturbation terms q KgJ(1 + f + f2/2)(1 + j) (2.17) = KgJ(1 + f + f2/2 +j + fj + f2j/2) Take the expected value of this expression, dropping terms higher than secondorder, resulting in an expression = KJ(1+ ank/2+ < fj >) (2.18) Consider this equation and the three terms in the parenthesis. The first is a "mean" term, the second is a "contribution" due to the variability of the conductivity field and the third is the log conductivity and head gradient crosscovariance. Gelhar [1993] evaluates this expression for isotropic log conductivity fields using spectral methods. The result, which is independent of the shape of the correlation function, is surprisingly simple [Gelhar, 1993] 
Full Text 
PAGE 1 SOLUTE TRANSPORT IN HETEROGENEOUS POROUS FORMATIONS By GEORGE GARY DEMMY, JR. A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1999 PAGE 2 For my parents. PAGE 3 ACKNOWLEDGMENTS Dissertations have but one name in the author slot, but anyone who has written one will readily testify that it represents the work and support of myriad people. Mary Hall freed my mind from many of the administrative hassles students encounter. Dawn Mendoza appointed, unappointed, and reappointed me numerous times throughout my tenure, and made the business end of things run as painlessly as possible. Ralph Hoffmann kept me in keys and desks and office space, and helped me get things wired when they needed wiring. Wally Hunter, system administrator and Unix hacker, deserves as much praise as I can heap upon him for keeping the network up and running and the data contained with in it safe and available. The quality of an institution and its faculty is reflected in the students they attract. I have enjoyed interacting with the students of the Hydrologic Sciences Academic Cluster, and this interaction has enhanced my efforts. Liyong Li, Yan Zhang, and Xavier Foussereau shared the same computer lab trenches as I, and my work cannot be separated from their influence. Many rewarding research experiences have been associated with collaborations made with Jim Jawitz. Jim is a colleague whose outlook and expertise is somewhat different than my own, but it is precisely these differences that have made the collaborations so rewarding. I shall seek to continue working with Jim through out my career. Dr. Andrew I. James provided access to computer codes he developed in the course of his own research and technical assistance. iii PAGE 4 Dr. Kenneth Campbell deserves special recognition not only for his exemplary service on my committee, but also for the example he has set as an engineer, a scientist, and a professor. I have appreciated Dr. Mike Annable's wit and example through the years and have enjoyed interacting with him in a wide variety of settings, whether in the classroom, on the football field, or elsewhere. Dr. Rao has always challenged me with questions that require looking beyond a simple formula and that evade a fillintheblank answer. Should I become scientist of any repute, Dr. Rao would be the man to blame. Dr. Kirk Hatfield is my source. He has been a source of wisdom, sympathy, advice, and funding. I have greatly enjoyed working with Dr. Hatfield as both teaching and research assistant. Dr. Gia Destouni taught a few stochastic subsurface hydrology lectures in the spring of 1994 that would forever change my thinking. She introduced the Lagrangian approach to transport with which I have been obsessed ever since. Dr. Destouni pursues fruitful collaboration vigorously, and this aggressive, yet collaborative, style is one that I now seek to emulate. Dr. Vladimir Cvetkovic has generously hosted three visits to the Kungl Tekniska Hogskolan in Stockholm, Sweden. His considerable generosity flows from his intense interest in stochastic Lagrangian transport. Dr. Sten Berglund has been a friend and mentor through this arduous journey. It is no exaggeration to say that without his help and encouragement this document would not exist. Dr. Wendy Graham deserves special recognition for her tireless faith and support. She has allowed me a free reign to pursue whatever interests I chose. In hindsight, this method coupled, with my own personality quirks and spectrum of interests, may not iv PAGE 5 be the most efficient for producing degrees and papers, but it has produced several interesting adventures. George and Ellen Demmy, my parents, have provided love and all manner of support throughout my life, academic and otherwise. They both instilled in me a love of learning and an insatiable curiosity about all things that serve me well, and will remain among my best traits. Finally, I'd like to extend my thanks and love to Celine Bufkin, for her love, encouragement, support, and vision. PAGE 6 TABLE OF CONTENTS page ACKNOWLEDGMENTS iii LIST OF TABLES ix LIST OF FIGURES x ABSTRACT xv CHAPTERS 1 INTRODUCTION 1 1.1 Primary Contribution 1 1.2 Historical Context 1 1.3 Research Goals 3 1.3.1 Lagrangian Velocity Mean and Covariance 3 1.3.2 Injection Mode Analysis 4 1.3.3 Heterogeneous Reaction Parameter 4 1.3.4 Application to "Reality" 4 2 LAGRANGIANEULERIAN FLOW FIELD RELATIONSHIP ... 5 2.1 Introduction 5 2.2 Flow Field Partitioning Strategies 7 2.2.1 Simple Subdivision 8 2.2.2 General Subdivision 10 2.3 Descriptions of Fields 12 2.4 Numerical Experiments 13 2.4.1 The Log Conductivity Field 14 2.4.2 Velocity Field Generation 18 2.4.3 Trajectory Generation 19 2.4.4 Monte Carlo Simulation 19 2.4.5 Sampling Strategy 20 2.4.6 Simulation Results 22 2.4.7 Log Conductivity Statistics 23 2.4.8 Velocity Statistics 29 2.4.9 Head and Head Gradient 40 2.5 Theory 40 2.5.1 Effective Conductivity 40 vi PAGE 7 2.5.2 Eulerian Velocity Field 42 2.5.3 Streamtube Weights 48 2.5.4 Control Plane Oriented Velocity Moments 51 2.5.5 Control Time Oriented Velocity Moments 55 2.5.6 Log Conductivity Moments 59 2.6 Evaluation of Simulation and Theoretical Results 60 2.6.1 Mean Velocities in Space 60 2.6.2 Velocity Variances in Space 60 2.6.3 Mean Velocities in Time 61 2.6.4 Velocity Variances in Time 62 2.6.5 Log Conductivity Statistics 62 2.7 General Conclusions 63 3 NONREACTIVE SOLUTE TRANSPORT 65 3.1 Introduction 65 3.2 Injection Mode and Streamtube Selection 65 3.2.1 Uniform Resident Injection 65 3.2.2 Injection in Flux 67 3.3 Integrated Streamtube Properties 67 3.3.1 Displacements 68 3.3.2 Travel Times 69 3.4 Solute Transport 70 3.5 Evaluation of Expressions 71 3.5.1 Displacements 71 3.5.2 Travel Times 72 3.6 General Conclusions 75 4 REACTIVE SOLUTE TRANSPORT 79 4.1 Introduction 79 4.2 Temporal Moments 79 4.3 Conceptual Overview of Solute Transport 80 4.4 Simulation of Reactive Solute Transport 83 4.4.1 Reaction Parameter 84 4.4.2 Reaction Flow Path 85 4.4.3 Reaction Flow PathTravel Time Cross Moment ... 93 4.5 Theory of Reactive Solute Transport 93 4.5.1 Mean Reaction Flow Path 96 4.5.2 Reaction Flow Path Variance 97 4.6 Discussion 98 4.7 Spatial Moments 99 4.7.1 Injection Mode 105 4.7.2 Discussion 106 vii PAGE 8 5 APPLICATION OF CONCEPTS 108 5.1 Introduction 108 5.2 Theory 108 5.2.1 Spatial Moments 109 5.2.2 Relative Degradation 110 5.3 Application: Simulation Ill 5.3.1 Results 113 5.3.2 Discussion 113 5.4 Application: Field Data 115 5.4.1 Results and Discussion: Center of Mass 116 5.4.2 Results and Discussion: Degradation 117 5.5 General Conclusions 117 6 GENERAL CONCLUSIONS 119 APPENDICES A EULERIAN VELOCITY DISTRIBUTION PROPERTIES 125 B VELOCITY LOG CONDUCTIVITY CORRELATION 127 C CONSISTENT FIRST ORDER APPROXIMATION 131 C.l Travel Time 131 C.2 Displacements 133 REFERENCES 134 BIOGRAPHICAL SKETCH 138 viii PAGE 9 LIST OF TABLES Table page 2.1 Comparison of target log conductivity values and observed values along areaweighted Eulerian trajectories 26 2.2 Comparison of target log conductivity values for areaweighted Eulerian trajectories in space and observed areaweighted Lagrangian trajectories in time 29 2.3 Comparison of a priori Eulerian velocity variance to observed values. The decimal representation of the exact fractional values of the predicted variance are provided for convenience 35 2.4 Comparison of observed log hydraulic conductivity variances to estimators. Observed value is zeroslope line regressed through data. The decimal representation of the fractional value is included for convenience 63 5.1 The center of mass in meters of steady plumes estimated by traditional trapezoidal integration of spatial concentration data and by the mean v/k estimator Equation (5.7) 116 5.2 Degradation rates relative to naphthalene and absolute degradation rates in reciprocal days as estimated by Equation (5.9) and as reported by Maclntyre et al. [1993]. NR denotes results not reported in Maclntyre et al. [1993]. Naphthalene degradation reported to be 0.0064 d1 117 ix PAGE 10 LIST OF FIGURES Figure page 2.1 Schematic representation of an aquifer 6 2.2 Simple subdivision of the flow field by one of four strategies: A) equalarea straight line, B) equalarea streamline, C) equalflow streamline, D) equalflow straight line 9 2.3 Four different partitioning strategies for the same flow field: A) area weight Eulerian trajectory, B) area weight Lagrangian trajectory, C) flow weight Lagrangian trajectory, D) flow weight Eulerian trajectory 11 2.4 Log conductivity distribution from a 41400 node realization generated by the turningbands method compared to the standard normal distribution. Target distribution parameters are //i n fc = 0 and of nJb = 1. 16 2.5 Estimated log conductivity longitudinal covariance estimated from a realization generated by the turningbands method compared to an exponential covariance function. Variance of log conductivity 0f nfc = 1 and correlation length Ai n fc = 1. Spacing for log conductivity values was Ai nfc /5 17 2.6 Mean log conductivity as function of displacement in mean flow direction for the four different trajectory collections and af nk 1/2 (a) an PAGE 11 2.10 Comparison of positiondependent mean velocity estimators to data for the four collections of trajectories for of nfc = 1/2 (a) and af nk = 1 (b). Velocities are normalized by U 2.11 Comparison of the positiondependent velocity variance estimators to data for the four collections of trajectories for af nk = 1/2 (a) and af nk = 1 (b). The variances are normalized by U 2 baf nk 2.12 Comparison of spatial Eulerian (AW/ET) and Lagrangian trajectory (FW/LT) velocity covariance data to covariance model. Data are taken from the cr, 2 nfc = 1 set of realizations. Covariances presented as dimensionless lag correlation functions 2.13 Comparison of timedependent mean velocity estimators to data for the four collections of trajectories for af nk = 1/2 (a) and a 2 nk = 1 (b). Velocities are normalized by U 2.14 Comparison of timedependent velocity variance estimators to data for the four collections of trajectories for of nfc = 1/2 (a) and a 2 nk = 1 (b). The variances are normalized by U 2 ba 2 nk 2.15 Comparison of temporal Lagrangian trajectory (AW/LT) velocity covariance data to covariance model. Data are taken from the of n k = 1 set of realizations. Covariances presented as dimensionless lag correlation functions 2.16 Normalized log initial velocities for twodimensional model aquifer compared to the standard normal distribution. The random variable was normalized by the sample mean and standard deviation to illustrate the lognormality 2.17 Comparison of velocity covariance expression Equation (2.25) to the approximation Equation (2.13) for o 2 nk = 1. Covariances are presented as dimensionless correlation functions 2.18 Components of the velocity covariance function 3.1 Comparison of mean displacements for the two injection modes or streamtube collections. The flowweighted Lagrangian trajectory (FW/LT) statistics correspond to a flux injection, and the areaweighed Lagrangian trajectory (AW/LT) statistics correspond to a uniform resident injection. The top figure contains the af nk 1/2 data and the bottom figure contains the of nfc = 1 data xi PAGE 12 Comparison of displacement variances for the two injection modes or streamtube collections. The flowweighted Lagrangian trajectory (FW/LT) statistics correspond to a flux injection, and the areaweighed Lagrangian trajectory (AW/LT) statistics correspond to a uniform resident injection. The top figure contains the af nk = 1/2 data and the bottom figure contains the of nA = 1 data. Variances are normalized by \ 2 nk PAGE 13 4.4 Mean reaction path observed values and estimators for different values of correlation parameter fi and different trajectory collections for the af nk = 1/2 set of simulations. Top figure is fi = {1,+1}. Middle figure is fi = {1/2, +1/2}. Bottom figure is fi = 0 (no correlation). Trajectory collections are areaweighted Lagrangian (AL) and flowweighted Lagrangian (AF) 4.5 Reaction path variance observed values and estimators for different values of correlation parameter fi and different trajectory collections for the of nfc = 1 set of simulations. Top figure is fi = {1,+1}. Bottom figure is fi = {1/2, +1/2}. Trajectory collections are areaweighted Lagrangian (AL) and flowweighted Lagrangian (AF). 91 4.6 Reaction path variance values and estimators for different values of correlation parameter fi and different trajectory collections for the fnfc = V 2 set of simulations. Top figure is fi = {1,+1}. Bottom figure is fi = {1/2, +1/2}. Trajectory collections are areaweighted Lagrangian (AL) and flowweighted Lagrangian (AF) 92 4.7 Correlation between the reaction flow path and travel time for different values of correlation parameter fi and different trajectory collections for the of nfc = 1 set of simulations. Top figure is fi = {1,+1}. Bottom figure is fi = {1/2, +1/2}. Trajectory collections are areaweighted Lagrangian (AL) and flowweighted Lagrangian (AF). 94 4.8 Correlation between the reaction flow path and travel time for different values of correlation parameter fi and different trajectory collections for the af nk = 1/2 set of simulations. Top figure is fi = {1,+1}. Bottom figure is fi {1/2, +1/2}. Trajectory collections are areaweighted Lagrangian (AL) and flowweighted Lagrangian (AF). 95 4.9 Assumed timedependent log conductivity covariance function along Lagrangian trajectories (Equation (4.29)) compared to simulation data for variances of log conductivity af nk {1/2,1}. The Lagrangian velocity covariance function is Equation (2.61) 103 4.10 Hypothesized timedependent Eulerian log conductivity covariance functions for 0\ nk = {1/2,1} compared to simulation data and Equation (4.29) 104 4.11 Mean reaction flow path as a function of nonreactive travel time compared to estimator Equation (4.25) for correlation parameter values fi = {1,0,1} 104 4.12 Reaction flow path variance a function of nonreactive travel time for correlation parameter values fi = { 1,0, 1} 105 xiii PAGE 14 4.13 Mean reaction flow path as a function of nonreactive travel time compared to estimator Equation (4.34) for correlation parameter values P = {1/2,0,1/2} (top figure) and f3 = {1,0,1} (bottom figure). PAGE 15 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy SOLUTE TRANSPORT IN HETEROGENEOUS POROUS FORMATIONS By George Gary Demmy, Jr. December 1999 Chair: Wendy Graham Major Department: Agricultural and Biological Engineering This work quantifies relationships between the spatial, or Eulerian, distribution of the properties of a chemically and physically heterogeneous porous medium and those as observed along the natural, or Lagrangian, trajectories that a fluid particle traces in a steady and irrotational flow. From these relationships, expressions that relate the transport of solutes through the porous medium along the natural trajectories to the aforementioned Eulerian distributions are developed. The effects of injection mode upon global measures of transport as reflected by the temporal moments of breakthrough curves and spatial moments of a solute plume are developed. The coupled effects of correlation of a linear equilibrium sorption to the underlying log hydraulic conductivity field and injection mode on the evolving temporal moments of mass breakthrough curve and the coupled effects of correlation of a firstorder decay coefficient and injection mode upon the spatial moments of a solute plume are examined. xv PAGE 16 CHAPTER 1 INTRODUCTION 1.1 Primary Contribution The primary contribution of this work is documenting a consistent approach to the development of estimators of solute transport in the Lagrangian framework. This consistent approach is to identify the collection of Lagrangian trajectories associated with a stationary Eulerian velocity field for which the Lagrangian velocities are also stationary. This is a modest extension of work pioneered by myriad hydrologists. The benefit derived from this contribution is enhanced understanding of transport of solutes in heterogeneous velocity fields. 1.2 Historical Context A grant to study the impact of heterogeneous source morphology upon subsequent solute transport initiated this research. In seeking a suitable framework and "toolset" with which to attack this problem, there appeared what seemed to be fundamental inconsistencies in the literature. Given the peculiar tendencies of the author, the "larger issues" could not be addressed until these apparent inconsistencies were resolved. Gedeon Dagan developed a robust theory of solute transport using a stochasticLagrangian framework [Dagan, 1982a, b]. In fact, an entire "school" of stochastic subsurface hydrology sprung from his work. Allen Shapiro and Vladimir Cvetkovic wrote Stochastic analysis of solute arrival time in heterogeneous porous media in 1988 [Shapiro and Cvetkovic, 1988]. The authors made an explicit assumption often implicitly made in the Lagrangian transport literature, namely that a fluid parcel deviates little from its mean trajectory in weakly 1 PAGE 17 2 heterogeneous conductivity fields, and that the Lagrangian and Eulerian velocities are approximately equal. They consistently applied this concept, developing expressions for arrival time means and variances for nonreactive solutes in the "Standard Model Aquifer," a hypothetical aquifer characterized by a heterogeneous conductivity field with a prescribed correlation structure and subject to certain prescribed boundary conditions (see Section 2.1). These equations predicted that the mean arrival time for solute subject to a uniform resident injection is given by where angle brackets <> denote an ensemble average operator and V h is the harmonic mean Lagrangian velocity. From the small deviation assumption, the Lagrangian and Eulerian velocities are approximately the same, thus the estimated harmonic mean Lagrangian velocity is equal to the estimated harmonic mean Eulerian velocity. However, the mean arrival time for uniform resident injection of solute after travelling several integral scales is given by where U is the arithmetic mean Eulerian velocity. This expression is also the "consistent firstorder approximation" (CFOA) of the travel time. Widespread adoption of the CFOA through out the stochasticLagrangian transport literature probably stems Feynman's requirement: we must reproduce what we already know. The harmonic mean of a positive definite process is always less than the arithmetic mean, so Equation (1.1) systematically overpredicts travel times for large displacements. Disturbing was that the estimators derived by Shapiro and Cvetkovic [1988] using clear, consistent, intuitive and logical methods seemed to violate Feynman's rule, especially in light of the success enjoyed by the Dagan school. PAGE 18 3 This harmonic/arithmetic mean discrepancy was explicitly noted by Dagan et al. [1992], who relegated the Shapiro and Cvetkovic [1988] result to applicability in areas "close enough to the input zone" [Dagan et al, 1992 p. 1374]. Cvetkovic et al. [1996] developed a semianalytical expression that described the transition from the "near field" harmonic mean velocity to the "far field" expression based upon the these known "endpoints" and an estimate of the Lagrangian velocity integral scale. Their detailed analysis of the Lagrangian velocity field revealed that the Lagrangian velocity is nonstationary in displacements in space, and resulted in a nonlinear propagation of the mean arrival time with distance. This was strange. Why was the Lagrangian velocity field nonstationary when the Eulerian was stationary? It was commonly assumed that the stationarity of one implied the stationarity of the other. What had started as a simple preliminary literature review resulted in a quandry. The paths followed in pursuit of this quandry led to the "research goals" of the dissertation. 1.3 Research Goals There were four broad objectives specified for this work. All centered around the flow of water and the transport of solutes in aquifers characterized by heterogeneous hydraulic conductivity fields. 1.3.1 Lagrangian Velocity Mean and Covariance A primary goal of this research was to develop a "Lagrangian covariance function" for the "Standard Model Aquifer". This was extended slightly to the mean and covariance of the spacestationary velocities along equal flowweight streamlines and timestationary velocities along equal areaweight streamlines. Quantification of this covariance function greatly simplifies the development of equations, or estimators, that describe the movement of water and solutes in the heterogeneous velocity fields associated with the "Standard Model Aquifer." PAGE 19 4 1.3.2 Injection Mode Analysis The stochasticLagrangian framework is convenient for the analysis of the effects of the method by which solute is introduced into the flow field, or injection mode, upon the subeqent transport. These effects of injection mode upon nonreactive solute transport are characterized by the mean and variance of a breakthrough curve, or mass arrival distribution, and by the mean and variance of mass displacement. A goal of this research was to develop estimators for these quantities. 1.3.3 Heterogeneous Reaction Parameter Set as a goal were estimators for the mean and variance of a heterogeneous "reaction parameter" integrated along trajectories of stationary velocty. These estimators were used to evaluate the coupled effects of injection mode and a heterogeneous linear equilbirium sorption process in a stochasticLagrangian framework. 1.3.4 Application to "Reality" A further goal was to develop a "traditional" estimator for the center of mass of a continuously injected solute plume subject to firstorder decay and to evaluate data from a natural attenuation experiment and numerical experiments. Estimator performance was evaluated using the results of the prevous developments. PAGE 20 CHAPTER 2 LAGRANGIANEULERIAN FLOW FIELD RELATIONSHIP Develop an intuitive judgment and understanding for everything. Miyamoto Musashi 1 2.1 Introduction Consider the steady, welldeveloped flow of water in an aquifer meeting Bear's [1972] definition of a porous medium. In most natural settings, gentle gradients induce a flow wellquantified by Darcy's law. Knowledge of the water fluid properties, relative permeability of porous medium, and potential gradients in space and time would assure a reasonable estimate corresponding water fluxes. 2 One way to record this knowledge would be to create maps of these different properties in space at different times. This framework is commonly called the Eulerian framework, and is the most common approach to hydrologic problems. Another way to record this knowledge is to map properties along moving coordinates in time. This approach is commonly called the Lagrangian framework, and its application may without controversy be described as less widespread than Eulerian approaches. Dagan [1989] and Gelhar [1993] review development of flow of fluids and transport of contaminants in porous media in both of these frameworks. 1 Musashi [1982 p. 49] 2 The concepts associated with the following discussion are quite general, and are applicable to a much larger variety of environmental transport problems. The emphasis upon groundwater follows from an interest in simplicity. 5 PAGE 21 6 Consider now flow of water in a more specific model aquifer: the steady, irrotational, and Darcian flow of water of homogeneous fluid properties in a finite rectangular twodimensional aquifer of heterogeneous hydraulic conductivity (see Figure 2.1). Mobile water 9 completely saturates the homogeneous pore space. An orthogonal \\\\\\\\\\\\\\\\\\\\\\\W & no flow boundary random conductivity field A\\\\\\\\\\\\\\\\\\\\\\\ c B no flow boundary average head gradient Figure 2.1: Schematic representation of an aquifer. h2 Cartesian coordinate system lies parallel to the aquifer boundaries, and x and y denote its directions. Parallel and impervious boundaries lie at y = 0 and y L y Constant potential boundaries located at x = 0 and x = L x maintain a constant average potential gradient J across the domain, such that water flows into the aquifer at x 0 and out of the aquifer at x = L x at a volumetric discharge rate Q (dimensions [L 3 T 1 ]). For an intuitive dimensional consistency, allow the aquifer some small and PAGE 22 7 uniform thickness in the vertical direction, say h. The total area of the vertical plane at x = 0 through which water enters the aquifer is A. Define the spatiallyaveraged groundwater seepage velocity at this inlet plane (IP) as jj I A ^ a \ da L da (2.1) = Q_ OA For a large class of hydrologic applications, water quantity is the sole concern. For these problems, the "route" taken by the water is of little interest, as is, say the distribution of ages or residence times. Spatial and temporal distributions of hydraulic properties hold some interest, but only in how they relate to effective bulk properties. However, a large class of practical hydrologic problems exist for which the trajectories of water parcels across the flow domain is of great importance. Many of these problems involve the fate and transport of contaminants through the aquifer system. The fundamental problem is this: to relate an Eulerian flow field to its Lagrangian, or trajectorybased, equivalent. 2.2 Flow Field Partitioning Strategies A streamline is a hydrodynamic entity that is at all points tangent to the velocity vector. For twodimensional flow, a streamtube is an object defined by two streamlines. Stream surfaces, the intersection of which form streamlines, define threedimensional streamtubes. Infinite collections of streamlines exist that subdivide the flow field into collections of streamtubes. Again, the interest of clarity suggests a focus upon twodimensional flow. For a lucid discussion of streamlines, stream functions, and stream surfaces in two and threedimensional flows, see Bear [1972]. Impervious boundaries are streamlines, thus the entire model domain is a streamtube by definition. Viewed as a "black box," Eulerian and Lagrangian distinctions have little meaning. This black box approach is the purview of reactor engineering, PAGE 23 8 an important and welldeveloped discipline, from which a great deal of insight and useful mathematics may be drawn. For an excellent review, especially with respect to concepts related to transfer functions, see Jury and Roth [1990]. For a large domain, the scale at which chemical reactions takes place is generally much smaller than that of the domain. That is, processes such as sorption, microbial and chemical degradation, and the like, are dependent upon chemical potentials across distances on the order of the pore scale to perhaps that of a "microDarcy" scale, defined as the smallest volume for which Darcy's law is generally applicable. However, the scale of interest tends to be much larger, say that of a well, or its capture zone or radius of influence, or perhaps the interface across which groundwater discharges into a surface water body, or, perhaps the Floridan aquifer. In order to consider these reactions in detail, it is necessary to resolve the larger domain into smaller subdomains. For modeling real systems, the scale and degree to which the domain is subdivided should be based upon the available data. 2.2.1 Simple Subdivision Consider the water flowing into the model aquifer across the inlet, or injection, plane (IP). Assume it is of interest to divide the domain into two subdomains. There are many ways to do this. One would be to divide the domain "down the middle," that is, at a line parallel to the x axis located y = L y /2 (curve A Figure 2.2). This division creates two subdomains of equal volume and gives an equal "volume weight" to each subdomain. The division also creates two subIPs of equal area and gives an equal "area weight" to each subdomain. For a perfectly homogeneous porous medium, the groundwater velocity vectors will be parallel to this line, and the discharges through each subdomain will be the same as well. Moreover, there will be no advective transfer of water across the boundary. However, the introduction of conductivity heterogeneity in the transverse direction will almost certainly lead to deviations of the streamlines, and to streamlines crossing the boundary. This implies that while PAGE 24 0 Figure 2.2: Simple subdivision of the flow field by one of four strategies: A) equalarea straight line, B) equalarea streamline, C) equalflow streamline, D) equalflow straight line. the volume of water contained in each subdomain is constant, the discharge through a plane perpendicular to the partition is a function of plane position in x. Consider now a streamline originating at (0,L y /2) (curve B on Figure 2.2). For a homogeneous medium, this corresponds to an equal volume partitioning. That is, the domain centerline is a streamline. However, transverse heterogeneity would likely induce changes in the velocity field such that the streamline would not remain on the centerline. For general transverse heterogeneity, the flow rate entering the two domains will not be equal, in general. Thus, the flow weight and the volume weight of the two subdomains are unequal in this case. Again, the magnitude of the areas through which water enters the subdomains are equal. Consider now a streamline originating at some coordinate (0, y 0 ) such that the volumetric flow rate is equal for both domains (curve C on Figure 2.2). This streamline corresponds to the stream function mv]= mL > ] : mo] (2.2) PAGE 25 10 This is an equal partitioning with respect to the flow field, or the stream function. This sort of partitioning arises quite naturally in fluid mechanics and in numerical modeling of heterogeneous flow systems using streamtube approaches (see, e.g., Theile et al. [1996]). The flow weight of each subdomain is equal, although the volume weights are different. A fourth scheme would be to extend a line from (0,y 0 ) to (L x ,y 0 ) (curve D Figure 2.2). In this case, the subdomains share only equal flow rates at x 0. In subsequent discussions, the phrase Eulerian trajectory will refer to such a straight line trajectory, and the phrase Lagrangian trajectory will refer to trajectory of a streamline. The phrase area weight will refer to an equal area partitioning scheme at a definition plane, and the phrase flow weight will refer to an equal discharge partitioning scheme. At this point, it is important to note that no mention has been made of boundary conditions as they relate to transport equations. These partitioning schemes are simply ways of subdividing a flow field. An area weight Eulerian trajectory scheme corresponds to traditional spectral approaches that assume small deviations about a mean trajectory. The use of trajectory is quite general, and does not imply a "natural" path or streamline. 2.2.2 General Subdivision Consider dividing the flow field into N partitions according to the strategies described (see Figure 2.3). Theile et al. [1996] demonstrated that streamlines may represent streamtubes. Therefore, the remainder of this discussion will focus upon streamline properties, with the implication that these properties may be extended to streamtubes as well. A primary purpose of subdividing the domain is to represent the continuum of domain properties with a discrete set of data. A regular finitedifference grid may be thought of as a set of points at equal spacing in space along a collection of areaweighted Eulerian trajectories. An alternative method would be to create spacings based upon travel times along the trajectory, where travel time is defined as PAGE 26 11 B Figure 2.3: Four different partitioning strategies for the same flow field: A) area weight Eulerian trajectory, B) area weight Lagrangian trajectory, C) flow weight Lagrangian trajectory, D) flow weight Eulerian trajectory. PAGE 27 12 the path integral of the inverse velocity over the trajectory r[l] = / Jo ds (2.3) V[C] a where s is a unit vector that is at all points the parallel to the trajectory. It is clear that this sort of Lagrangian description has little meaning in a static system, and is in this sense less general than Eulerian descriptions. Thus, Lagrangian description is a tool reserved for dynamic systems and is a supplement to Eulerian descriptions. It is worthy to note the distinction between a description, a general approach, and a trajectory, a specific entity used in a description. To summarize, a primary objective is to approximate a steady and irrotational flow field continuum with a discrete set of trajectories. The trajectories follow either a straight line (Eulerian trajectory) or a streamline (Lagrangian trajectory). The trajectories carry either an equal flowweight or an equal injection areaweight. Consider once again the model aquifer depicted in Figure 2.1. As mentioned in the introduction of this chapter, the aquifer has certain features in which a hydrologist might be interested, e.g., the intrinsic permeability or the porosity. These features may be thought of as a field in which parameters vary continuously. This is a fundamental concept of continuum mechanics and is discussed in a hydrologic context by Bear [1972], Dagan [1989], and Gelhar [1993]. The ultimate description of any such field would be a map of parameter values valid at any point, time, or scale. In lieu of this "ultimate description," stochastic hydrologists characterize spatial heterogeneity with conditional probability density functions that reproduce the dominant characteristics of the field and observed values at given locations. The entire subject of geostatistics is devoted to this end (see Journel and Huijbregts [1978]). An Eulerian field refers to areaweighted Euleriantrajectory oriented field. A Lagrangian field refers to a flowweighted Lagrangiantrajectory oriented field. 2.3 Descriptions of Fields PAGE 28 13 It is instructive to illustrate these concepts with an example and to interpret the results of this illustration. The illustration of these concepts requires a heterogeneous flow field and collections of trajectories, both Eulerian and Lagrangian, that correspond to this flow field. The requirement that the flow field be based upon "reference fields" of known properties enhances the possibility of theoretical interpretation drawn from the large body of literature concerning flow and solute transport in idealized porous media. A further requirement that our reference fields bear some resemblance to what might be expected in nature enhances the possibility of useful applicability. A numerical simulation code modeled steady groundwater flow in a hypothetical heterogeneous twodimensional aquifer. A particle tracking code traced Eulerian and Lagrangian trajectories across this aquifer, recording data that included the local groundwater velocity and local hydraulic conductivity at intervals in space and time along the different trajectories. Statistics, specifically mean, variance, and covariance, summarize the results of these "measurements." The conceptual view of the aquifer is that given by Figure 2.1. To generalize the system, correlation length, or integral scale, of the log conductivity field Ai n fc serves as a characteristic length with which to normalize length scales. The characteristic time \\nk/U serves to normalize time scales, where U is the arithmetic mean velocity defined in Equation (2.1). The definition of an integral scale is the integral from zero to infinity of the correlation function R, or, 2.4 Numerical Experiments (2.4) PAGE 29 14 The design criteria for the model aquifer are as follows In A: correlation length X\ n k = 1 effective conductivity K e = 1 average gradient J = 1 /50 mobile water content 6 = 1/5 From these criteria and Darcy's law, the areaaveraged velocity is U = K e J/6 = 1/10. Clearly, the parameter values are artificial and selected for convenience, but values for K e in md~ l and 6 fall in ranges typical for sand [Domenico and Schwartz, 1990]. Two 5000 replicate Monte Carlo simulations of aquifers with log conductivity variances of 1/2 and 1, respectively, comprised the suite of aquifer simulations. The following sections describe the experimental design and the results. 2.4.1 The Log Conductivity Field The foundation upon which to build the illustration is a heterogenous hydraulic conductivity field. A commonlyused stationary lognormally distributed and exponentially correlated isotropic model serves as the Eulerian hydraulic conductivity reference field (see e.g., Freeze [1975] or Gelhar and Axness [1983]). Working from a hypothesized conductivity field is preferable to generating a hypothetical velocity field with specified characteristics (e.g., Bellin et al. [1994]) because an objective of this work is to describe the underlying conductivity field along Lagrangian trajectories and relate its Lagrangian description to its Eulerian properties. For completeness, it is pertinent to briefly review description of a heterogeneous field by treating the process as being "random" Measured values of hydraulic conductivity of virtually any natural porous formation will vary from place to place. In the absence of measurement error, this variability in the measured values reflects the heterogeneity of myriad physical factors and processes that control the hydraulic conductivity. As mentioned, measurements of hydraulic conductivity in the field often follow a lognormal distribution. In these cases, it is useful to conceive of the PAGE 30 15 conductivity field as a realization of a lognormal random process. The process is said to be stationary if the ensemble mean and variance do not vary in space. The process is nonstationary if some trend exists in the mean or variance. In practice, only single realization is available, and the "ensemble" properties are ultimately unknowable. However, if the field is large compared to the scale of variability, one may assume ergodicity, or, that the ensemble parameters may be estimated from a sample drawn from the single realization since there exists a the possibility of a large number of statisticallyindependent observation points. For a more detailed development of these concepts from the hydrologist's perspective, see Gelhar [1993]. For a more formal development of these concepts, see Papoulis [1991]. The model of hydraulic conductivity variability is a stationary, exponentiallycorrelated isotropic lognormal distribution. A lognormal random variable is one for which its natural logarithm is normally distributed. As such, three parameters summarize the statistical properties of the adopted conductivity model, namely the mean of log k /link, the variance of log k af nfc and the correlation length of log k \\ n kThe covariance model is given by where r is the separation between two points in the log conductivity field. Use of this simple covariance model requires stationarity of variance, and usually implies a stationarity of the mean as well. The moments of the "real space" conductivity field are given by where n is the order of the moment. In the absence of a perfectly described physical aquifer, the method of turningbands generated equallylikely realizations of a hypothetical conductivity field with C\ nk [r] = af^expfr/Ainfc] (2.5) <& n >= exp[n//i nJt + n 2 af nk /2] (2.6) PAGE 31 16 prescribed ensemble statistics [Tompson et al, 1989]. Turningbands was capable of producing two and threedimensional fields with the desired statistical properties and was available as a portable and easily modified code. Figure 2.4 illustrates the lognormality of the a typical twodimensional log conductivity field realization. Figure 2.5 compares an estimate of the longitudinal covariance to the hypothesized model. In order to capture the features of the log conductivity spatial variability, the 3210 1 2 3 log conductivity Figure 2.4: Log conductivity distribution from a 41400 node realization generated by the turningbands method compared to the standard normal distribution. Target distribution parameters are fj,\ n k = 0 and af nk = 1. turningbands algorithm generated five log conductivity observations per correlation length. This level of discretization falls in the range of other similar sets of simulations. Cvetkovic et al. [1998] used a discretization of 4 nodes per Ai n fc for simulations of fields with of nfc as large as 1. Bellin et al. [1992] used a discretization of 8 nodes per Ai n fc for simulations of fields with a^ nk as large as 1.6. PAGE 32 17 0.4 0.2 1 \ 1 1 estimated In k covariance + exponential covariance \ t i i 0 12 3 r/\\nk Figure 2.5: Estimated log conductivity longitudinal covariance estimated from a realization generated by the turningbands method compared to an exponential covariance function. Variance of log conductivity af nk 1 and correlation length \\^ = 1. Spacing for log conductivity values was Ai nfc /5. Trajectories covered several characteristic lengths X\ n k and several characteristic times \\ n k/U to capture both nearfield and farfield behavior in space and Lagrangian time. The first requirement was easy to enforce by making the domain size equal to or greater than the required number of correlation lengths. The second requirement was somewhat more problematic, since the distance traveled in some number of \\ n k/U varies from trajectory to trajectory. Preliminary simulations indicated that a domain 40Ai n fc long was adequate to capture travel up to 10 characteristic time scales \\ n k/U in 2000 trajectory realizations generated in the most heterogeneous field considered (of njt 1). Thus, the selected domain extent was 40Ai n jt. The conceptual model aquifer is a large, bounded, primarily twodimensional, steady flow field. As such, it should exhibit a variety of head, head gradient, and velocity nonstationarities due to the boundary conditions (see Figure 2.1). While PAGE 33 18 researchers such as Rubin and Dagan [1988], [1989], and Osnes [1998] have investigated such nonstationarities, the additional complication of a rigorous consideration of these sorts of nonstationarities would unduly obfuscate the very simple concepts illustrated here. To minimize the effect of these nonstationarities, "solute transport" is constrained to take place some distance from the boundaries of the numerical domain. Following the example of Bellin et al. [1992], a 3Ai nfe "buffer" was established on both sides of the required 40 for a total longitudinal domain length of 46Ai nfc Preliminary simulations indicated that the maximum transverse displacement of a streamline in a af nk = 1 field was around one third the longitudinal displacement. For a streamline originating on the aquifer centerline parallel to the x axis, an aquifer width of 36Ai n fc was deemed adequate to contain large transverse displacements over a 40Ai nfc longitudinal displacement and provide a similar "buffer" to mitigate boundary effects. Combining the "physical" domain requirements with the discretization required to capture the log conductivity variability resulted in a numerical grid 230 by 180 for a twodimensional aquifer and a numerical domain with 41400 nodes. A threedimensional analysis under the same assumptions would have required a numerical domain of more than 7 million nodes. Simulations of this magnitude were untenable in light of the available computational resources. 2.4.2 Velocity Field Generation A mixed finiteelement scheme solved coupled velocity and pressure equations derived from Darcy's law and the continuity equation. The code took input parameters that matched the design criteria and the log conductivity field generated by turningbands. Andrew I. James wrote and tested the core of the code. James et al. [1997] used a derivative version of this code to generate flow fields for synthetic tracer tests. A mixed finiteelement solution to the flow field was specified, following the results and arguments of Mose et al. [1994]. The mixed finiteelement scheme is known PAGE 34 19 to provide an accurate velocity field solution for groundwater flow in heterogeneous conductivity fields, one that is more accurate than other, more traditional, schemes [Mose et ai, 1994]. 2.4.3 Trajectory Generation An objective of this research was to record data, such as local velocity and and conductivity, along different collections of trajectories, both Eulerian and Lagrangian. These trajectories were traced using a "particletracking" scheme based upon a semianalytical method presented by Pollock [1988]. This method assumes that velocity varies linearly over an element, and that the velocity in some direction, say i, is independent of displacements in orthogonal directions. From these assumptions, it is a relatively simple exercise to determine the displacement of a "particle" in a linear velocity field for a given time, or the time required to make a certain displacement. A computer code recorded particle position, time, total path length along trajectory, velocity, local conductivity, and a value from another random field with the same statistical properties as, but uncorrelated to, the log conductivity field at several control planes equallyspaced along the mean flow direction and at several "control times" after injection along both trajectories. This uncorrelated random variable was recorded for reactive transport simulations, to be discussed in subseqent chapters. 2.4.4 Monte Carlo Simulation The possibility of theoretical interpretation of the results is enhanced by requiring nearergodic samples from the population. That is, the samples must reproduce the population statistics. To do this, the number of statistically independent replicates drawn from the velocity and conductivity fields must be quite large. Rather than attempt a massive singlerealization, the following method generates a large number of statistically independent realizations PAGE 35 20 generate a large, twodimensional conductivity field with the desired statistical properties generate a flow field by imposing a constant mean gradient across the saturated conductivity field trace the trajectories of two particles injected at a common "control location", one trajectory is a straight line or Eulerian trajectory and the other a streamline or Lagrangian trajectory, and record properties at equal lags in space and in time along trajectories that pass through "inner domains" that are relatively undisturbed by nonstationarities induced by the boundaries At this juncture, it is perhaps worth mentioning a minor "philosophical" point. These simulations do not represent a real aquifer, per se, but rather a hypothetical aquifer with known properties. The goal is understand the convolution of the physics of the groundwater flow and a known heterogeneity that has analogs in nature, and it is hoped that the reader finds some merit to this modest goal. The 41400 conductivity values needed for any single realization do not represent 41400 measurements taken from some real aquifer. The Monte Carlo aspect of this experiment reenforces this point. 2.4.5 Sampling Strategy Consider "dropping particles" into a flow field, with no particular regard to location. These particles move along streamlines without interaction with the porous medium. Any one location into which one is dropped is as equallylikely as any other. The trajectories associated with these particles carry an equal area weight, analogous to equalarea subdivision of an injection plane. Consider now, the same disposal of particles, but this time an area is assigned to each particle at the injection point such that an equal discharge is associated with each particle, and its associated trajectory. For a homogeneous mobile water content, this area is inversely proportional to the local velocity normal to the area. Thus, each particle may be thought to carry an equal PAGE 36 21 flow weight. In this way, it is possible to "map" from one collection of trajectories to another. At a given control plane or time, the areaweighted mean and variance of some property over N realizations, say v was calculated using the familiar estimators 3 1 N i=i 1 N vark] = J^(Y1 v i ~ N ^ ( 2 8 ) i=i While these particular estimators seem obvious, it is important to realize that "dropping the particle" into the flow field assigns an equal area weight to each trajectory. For a single realization, this is equivalent to starting each particle with equidistant spacing. These estimators are based upon an assumption of equallylikely statisticallyindependent samples (see, e.g., Mood et al. [1974]). In order to draw even a few hundred independent samples from a single realization, the required flow domain would be enormous. This was a primary consideration for adopting Monte Carlo simulation using a minimal domain. The flowweighted mean property v was calculated via Ei=i<[o] 3 Press et al. [1988] have criticized the variance estimator in Equation (2.8) as being prone to accumulation of roundoff error for large data sets. The rather substantial size of the data sets are naturally a concern. Comparision of the variance calculated by a corrected twopass method presented by Press et al. [1988] and the "traditional" estimator yielded no difference between the two methods within the first five significant digits. PAGE 37 22 where u;[0] is the initial velocity recorded for trajectory i. The variance by Ei=i w An alternative method for generating a collection of streamlines would be to identify the streamline corresponding to the "center" stream function expressed in Equation (2.2). In this case, the estimators of the equal flow statistics would be Equations (2.7) and (2.8). The equal area statistics would be given by = E£Wm[o] (211) Et=i i/i[o] and varM = ^tM^lum E£i.l/*[0] The "reference" collection of the trajectories selected uniformly in space was specified in the simulations for convenience. 2.4.6 Simulation Results The following sections present the results of two 5000 replicate Monte Carlo simulations. Detailed discussion of these results is deferred until after the development of some theoretical models with which to analyze these results. However, some heuristics are provided to familiarize the reader with some of the concepts developed in subsequent sections and chapters. A summary of the notation used with the following results provides clarity. AW denotes an areaweighted trajectory and FW denotes a flowweighted trajectory ET denotes an Eulerian trajectory and LT denotes a Lagrangian. A particular weighttrajectory combination is concatenated the with a /. Thus, an areaweighted Eulerian trajectory is denoted as AW/LT. For a given simulation, that is, for one set of 5000 replicates for a log conductivity field of some given a? nk the statistics of the log conductivity and velocity along the PAGE 38 23 four different trajectories are calculated at intervals in time and in space. Thus, for any one property, say log conductivity, and any one statistic, say the mean, there will be two sets of results, namely that for af nk = 1/2 and that for af nk = 1. Associated with each set will be two subsets, namely statistics along the four trajectories recorded at space or time intervals. For convenience, the notation and the trajectories to which each notation corresponds as depicted in Figures 2.2 and 2.3 is summarized here: AW/ET areaweighted Eulerian trajectory (Figures 2.2 and 2.3 A) AW/LT areaweighted Lagrangian trajectory (Figures 2.2 and 2.3 B) FW/LT flowweighted Lagrangian trajectory (Figures 2.2 and 2.3 C) FW/ET flowweighted Eulerian trajectory (Figures 2.2 and 2.3 D) 2.4.7 Log Conductivity Statistics The log conductivity field serves as the "foundation" upon which the simulations are built, as well as the theory. The data for the areaweighted Eulerian trajectories correspond to those design criteria specified. The log conductivity mean and variance for AW/ET appear stationary in space (see Figures 2.6 and 2.7). That is, they appear to "bounce around" their "target" values with no apparent trend. Moreover, the statistics of FW/LT appear likewise stationary, but the mean is somewhat larger and the variance slightly smaller than AW/ET. The magnitude of the differences appear proportional to of nfc A comparison of the "target" values, e.g., the specified inputs into the turningbands algorithm, and the observed ensemble values along the trajectories is presented in Table 2.1. The "observed" entry in Table 2.1 represents the statistic average along the trajectory, as estimated by regressing a zeroslope line through the data using a generalized leastsquares method. There is close correspondence between the target and observed values, supporting an ergodic hypothesis for the log conductivity statistics. That is, sample statistics are an accurately estimation PAGE 39 24 0.6 c =1 0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 0.6 0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 +__ + FW/LT + AW/LT x FW/ET AW/ET Eq. 2.72 reference + X ix, S + + + + ** 5* ** + x x X x I L. .I*.... JW.! A. ..S. "flJ.l X FW/LT AW/LT FW/ET AW/ET Eq. 2.72 reference 10 20 30 40 Figure 2.6: Mean log conductivity as function of displacement in mean flow direction for the four different trajectory collections and af nk = 1/2 (a) and af ak = 1 (b). PAGE 40 25 FW/LT AW/LT FW/ET AW/ET reference Eq. 2.73 1.1 1.05 1 0.95 m 0.9 b 0.85 0.8 0.75 0.7 *x V< + + ** X* X $ + 10 20 30 FW/LT AW/LT FW/ET AW/ET reference Eq. 2.73 40 Figure 2.7: Log conductivity variance as function of displacement in mean flow direction for the four different trajectory collections and af nk = 1/2 (a) and a? ak = 1 (b). Variances are normalized by target variance of Eulerian log conductivity field of nk PAGE 41 26 Table 2.1: Comparison of target log conductivity values and observed values along areaweighted Eulerian trajectories. of the population statistics. Again, this is important for a quantitative theoretical analysis of these results. The AW/LT and the FW/ET exhibit distinctly nonstationary behavior in the mean and less distinct behavior in variance as a function of displacement along the mean flow direction (see Figures 2.6 and 2.7). At the injection point, the AW/ET and AW/LT statistics are identical and the FW/ET and FW/LT statistics are identical, because these trajectories carry the same weight and the points in consideration are identical (see e.g., Figure 2.2). At very large displacements, one would expect the statistics of like trajectories to assume similar characteristics, regardless of weight. As the trajectory traverses several integral scales that characterize the heterogeneity in question, "information" about its starting location, is diminished. In the case of the FW/ET, the statistics "start" at the same place as FW/LT, and eventually "end" at the same place as the AW/ET. Similarly, the AW/LT statistics start where do the AW/ET and end where do the FW/LT. In fact, the data exhibit this behavior (see Figures 2.6 and 2.7). The log conductivity statistics in time show somewhat different behavior. Only the AW/LT appear stationary in time (see Figures 2.8 and 2.9). There is an important observation to be gleaned from these results. The AW/LT statistics, in addition to appearing stationary, also appear to have a similar value to the input AW/ET statistics (see Table 2.2). Recall, the AW/LT statistics are clearly nonstationary in space, and are clearly distinct from those of the AW/ET. This observation greatly parameter target observed 0 0.003 0 0.005 1/2 0.498 1 1.002 PAGE 42 27 FW/LT + AW/LT x FW/ET AW/ET o est AW/LT 0.4 0.2 0.2 0.4 m x x xx xx X x X x* ++ ++ + ++ ++V + ++ + + tt ++ + o x 0! rx~ i^XXvx n *" x *xxxV "n xxx 13 Xxx xxxxx xx IK D *x B X xxxxx x xx x x 5! x x xx&x* *, FW/LT + AW/LT x FW/ET AW/ET est AW/LT 4 6 Ut/I 10 Figure 2.8: Mean log conductivity as function of travel time to control planes along mean flow direction for the four different trajectory collections and af nk = 1/2 (a) and a lk = 1 (b). PAGE 43 28 1.1 1.05 PAGE 44 29 Table 2.2: Comparison of target log conductivity values for areaweighted Eulerian trajectories in space and observed areaweighted Lagrangian trajectories in time. parameter target observed MinfcK^ = 1/2) 0 0.000 HinkKk = 1) 0 0.017 1/2 0.497 1 1.000 simplifies relating the Eulerian and Lagrangian fields, as shall be shown. The "longtime" asymptotic mean log conductivity values for both Eulerian and Lagrangian trajectories are somewhat less than their "largedisplacement" counterparts, and the magnitude of the difference is proportional to o^ nk (compare Figures 2.6 and 2.8). This is consistent with the harmonic averaging implicit with a timeaveraged timedependent velocity and the arithmetic averaging implicit with a spaceaveraged spacedependent velocity. Recall, we map the quantities onto the trajectories by following velocities, whether along "natural" streamlines (Lagrangian trajectories) or not (Eulerian trajectories). A simple heuristic is as follows. Low conductivities correspond to low velocities. For an equaltime spacing (or equallylikely observations in time) more "measurements" will come from low velocities, since the particle tracing the trajectory spends more time going slow. For an equaldistance spacing more measurements will come from high velocities, since the particle "covers more ground" while going fast. This is perhaps more intuitive when discussed in terms of the velocity field itself. 2.4.8 Velocity Statistics The mean velocities along the different trajectories exhibit a similar behavior as do the mean log conductivities along the same trajectories (see Figure 2.10 and compare Figure 2.6). The AW/ET exhibits apparently stationary behavior in the mean and variance (see Figure 2.11) in accordance to established theory (e.g., Dagan [1989]). The AW/LT exhibits nonstationary behavior in the mean as observed by Cvetkovic et al. [1996]. Again, the AW/ET and FW/LT trajectories appear to be stationary PAGE 45 30 1.45 1.4 1.35 1.3 1.25 1.2 1.15 1.1 1.05 1 0.95 1.45 1.4 1.35 1.3 1.25 1.2 1.15 1.1 1.05 1 0.95 ,x *** x ** x x x x FW/LT + AW/LT x FW/ET AW/ET Eq. 2.35 Eq. 2.43 Eq. 2.44 reference 10 20 30 40 ++ ++ +\ *J .+:. i.^.$a: x **, *x % x FW/LT + AW/LT x FW/ET AW/ET o Eq. 2.35 Eq. 2.43 Eq. 2.44 reference 10 20 30 40 Figure 2.10: Comparison of positiondependent mean velocity estimators to data for the four collections of trajectories for af nk = 1/2 (a) and af nk 1 (b). Velocities are normalized by U. PAGE 46 31 2 b 1.8 1.6 1.4 b 1.2 0.8  *++*%++ ++X x x+ x *+ + + + x *+ X x + + X x x + *** 3 .:as.vV £3* !_^ FW/LT + AW/LT x FW/ET AW/ET Eq. 2.37 Eq. 2.55 Eq. 2.57 Eq. 2.13 10 20 z/Alnfc 30 40 1.8 1.6 CM s b 1.4 1.2 .' x .' x Xx + x + x x x x *x +t 5 x x xx X *x + x x + X FW/LT + AW/LT x FW/ET AW/ET Eq. 2.37 Eq. 2.55 Eq. 2.57 Eq. 2.13 0.8 10 20 30 40 Figure 2.11: Comparison of the positiondependent velocity variance estimators to data for the four collections of trajectories for cr, 2 nfc = 1/2 (a) and af nk 1 (b). The variances are normalized by U 2 baf nk PAGE 47 32 in both the mean and variance. The velocity variances exhibit behavior to similar to that of the means for their corresponding trajectories and quite different than the behavior exhibited by the log conductivity variances (compare Figures 2.11 and 2.10 and compare Figures 2.11 and 2.7). There is a qualitative similarity between the nonstationary velocity means and the nonstationary velocity variances. A heuristic is that the velocity statistic "starts" at one value "in the nearfield" and transitions to an asymptotic value "in the farfield." The following sections demonstrate that these endpoints and the transition are predictable. The qualitative dissimilarity between the log conductivity variance and velocity variance is an artifact of examining the log conductivity values on the one hand and the "real space" velocity values on the other. Consider the F W/LT log conductivity mean and variance for the af nk = 1 simulation. The FW/LT mean log conductivity is larger than that of the AW/ET, but the variances are roughly the same value. In fact, regressing a zeroslope line through the data results in a "mean" FW/LT variance that is slightly less than that of the AW/ET. However, the "real space" conductivity mean and variance are both functions of both the mean and variance of the log conductivity process (recall Equation (2.6)). A cursory examination of the data in conjunction with Equation (2.6) reveal that the mean and variance of the "real space" conductivity along FW/LT are both higher than their AW/ET counterparts, and the magnitude of this difference is proportional to af nk The mean velocities in Figure 2.10 are normalized by the arithmetic mean velocity U predicted from the input parameters for the flow simulator. That the mean velocity along AW/ET closely follows 1 is a good indication that the effective conductivity expression is accurate. Notice that the FW/LT and farfield AW/LT values are higher than AW/ET, and the magnitude is proportional to the variance of the log conductivity (see Figure 2.10). This behavior is intuitively correct. Consider a PAGE 48 33 particle injected into a flow field. Constraint of the trajectory to an Eulerian trajectory forces the trajectory through lowvelocity areas that a natural streamline, or Lagrangian trajectory, would have otherwise bypassed. Of course, some streamlines pass through even the lowest velocity areas, but these are few in an equalflow sense. Thus, equalarea weighting preferentially weights low velocity zones. That there are implications for sampling schemes in heterogeneous media, especially in context of multilevel samplers, should be intuitively obvious to even the most casual observer. A discussion of some of these implications is made in a subsequent section, after these concepts are developed further in a more quantitative and theoretical framework. It is appropriate to evaluate the adopted expression for the AW/ET longitudinal covariance function suggested by Dagan and Cvetkovic [1993], namely Cua[r) = al k b[e]U 2 exp [b[e)r / \ Xnk ] (2.13) where r is the separation between points along the ET, and b[e] is a function of the anisotropy ratio e. The anisotropy ratio itself is a function of the correlation lengths of the log conductivity in the different directions. The primarily consideration is that b may be derived from theoretical considerations and is not some ad hoc fitting parameter. The value of this function for a twodimensional isotropic aquifer is 3/8. Again, for further background on this function, see Dagan and Cvetkovic [1993]. In the limit r = 0, the covariance reduces to the point variance, which for the model case is a\ = 3[/ 2 of nfc /8. For the remainder of this discussion, a reference to b contains its dependence upon e implicitly The subscript ua may be decoded as follows: for a subscript tw, the t is the trajectory and property (u for Eulerian velocity, and v for Lagrangian velocity) and the w is the weight (a for area and / for flow). Thus, C ua is the covariance function for velocities along Eulerian trajectories. We present this covariance function and that covariance function estimated from the assumedstationary collection of AW/ETs in Figure 2.12 The perceptive observer will note PAGE 49 34 AW/ET + FW/LT x Eq. 2.13 X \ + '. X \ + \ X '. + i\ + *. + 'x.. 0 5 10 15 r/Ai n fe Figure 2.12: Comparison of spatial Eulerian (AW/ET) and Lagrangian trajectory (FW/LT) velocity covariance data to covariance model. Data are taken from the of nfc = 1 set of realizations. Covariances presented as dimensionless lag correlation functions. that the analytical expression overestimates the covariance at short lags and underestimates the covariance at longer lags. This discrepancy and its implications are discussed in some detail in a subsequent section. The velocity variances in Figure 2.11 are normalized by the zerolag covariance predicted by placing the input simulation parameters into Equation 2.13. The AW/ET velocity variance appears wellpredicted by this expression, as evidenced by a relatively good fit to 1 for both values of o? nk in Figure 2.11. A regression of a zeroslope line through these data indicates that the a priori estimate under predicts the observed values by about five percent for the larger o^ nk (see Table 2.3) The apparent stationarity of the velocities along the FW/LTs suggest estimation of the covariance function for these velocities as for the AW/ET (see Figure 2.12). The correlation length of the FW/LT is very close to that of the AW/ET, although the actual value is probably slightly higher. Recall the definition of the correlation length as given by 0.8  06 'M n OJ 8 0.4 0.2 PAGE 50 35 Table 2.3: Comparison of a priori Eulerian velocity variance to observed values. The decimal representation of the exact fractional values of the predicted variance are provided for convenience. para meter predicted observed o 2 u {o'i k = 1/2) 3/1600=0.001875 0.001871 ^(^i 2 nfe = 1 ) 3/800=0.00375 0.00394 Equation (2.4). The correlation function at each considered lag along the estimated functions is systematically larger for F W/LT, even though the magnitude of the difference is small. Integrating these functions would result in a slightly larger value for FW/LT, and thus a larger correlation length. In the case there is no variability in the log conductivity, Eulerian and Lagrangian trajectories coincide. Moreover, an equalarea weight corresponds to an equal flow weight. Thus, the estimated correlation functions from "either set" would be identical. From this conclude that there is some dependence of the correlation length of the FW/LT trajectories upon af nk Cvetkovic et al. [1996] observed similar behavior in log velocity correlation functions estimated from AW/LTs with a stationarity assumption for the nonstationary AW/LT velocities. Thus, this is an slight extention of this work by identifying the stationary streamlinebased velocity field from which to work, and work from that field. Moreover, this supports their findings of an increasing correlation length with af nk Understanding the relationship between the variability of the log conductivity field and the correlation length of the associated velocity field is of fundamental importance to completely understanding the nature of flow in heterogeneous media. Sadly, these "secrets" remain uncovered and must be left this to future research efforts. The velocity means as a function of time exhibit qualitatively similar behavior as do the log conductivity means in time (see Figure 2.13 and compare Figure 2.8). Again, only the AW/LT exhibits stationary behavior in time for both the mean and variance (see Figure 2.14) as did the AW/LT log conductivity values. The mean and variance of the AW/LT in time are approximately those of the AW/ET in space PAGE 51 36 1.4 1.3 1.2 1.1 0.9 0.8 0.7 1.4 1.3 1.2 1.1 FW/LT + AW/LT x FW/ET AW/ET o Eq. 2.63 Eq. 2.59 Eq. 2.67 Eq. 2.66 \ ^effect, ** 4 6 10 0.9 0.8 0.7 FW/LT AW/LT FW/ET AW/ET Eq. 2.63 Eq. 2.59 Eq. 2.67 Eq. 2.66 *x x xxxxx x 1 ; 3 D **w C =G :: 4 6 tu/x Xuk 8 10 Figure 2.13: Comparison of timedependent mean velocity estimators to data for the four collections of trajectories for af nk = 1/2 (a) and a^ nk = 1 (b). Velocities are normalized by U. PAGE 52 37 2 1.8 1.6 Q 1.4 J Oi p 1 b 0.8 0.6 0.4 2 1.8 1.6 1? 1.4 c > 12 FW/LT + AW/LT x FW/ET AW/ET n Eq. 2.65 Eq. 2.61 Eq. 2.70 Eq. 2.69 x x x .. xx x++x + ^* ? x. x xx* b X X X £ x$*£* CM s b 4 6 tU/Kk i 1 10 FW/LT AW/LT J. FW/ET AW/ET i Eq. 2.65 \ Eq. 2.61 Eq. 2.70 Eq. 2.69 A, y V + V ++ \ \ ++ + X + \ x xx v v + + + + + x, x x + + *v X ****x x X jx"Sf *.; 5J 0.8 h b ''....^B^ggn^ s 08 B 58 aB 8 B a SBfigg 8 8 B a b 0.4 4 6 tu/x lnk 10 Figure 2.14: Comparison of timedependent velocity variance estimators to data for the four collections of trajectories for af nk 1/2 (a) and af nk = 1 (b). The variances are normalized by U 2 baf nk PAGE 53 38 for both simulations. This result is important, because it implies the possibility of mapping from AW/ET to AW/LT. We estimated the covariance function for the stationary AW/LT velocities (see Figure 2.15). As an estimate of the covariance o V o u 4 6 tU/K k 10 Figure 2.15: Comparison of temporal Lagrangian trajectory (AW/LT) velocity covariance data to covariance model. Data are taken from the af nfc = 1 set of realizations. Covariances presented as dimensionless lag correlation functions. function, the substitution r = Ut was made in Equation 2.13 and the result compared to the estimated covariance function (see Figure 2.15). The closedform predictor behaves similarly to its spatial counterpart: earlytime correlation is underpredicted and longtime correlation is overpredicted (compare Figures 2.15 and 2.12). For completeness, one other aspect of the velocity field was examined, namely the distribution of initial velocities (see Figure 2.16). The logarithm of the initial velocity was rankordered and subject to the transformation In v In v \nv' J \n v (2.14) PAGE 54 39 3210123 log velocity Figure 2.16: Normalized log initial velocities for twodimensional model aquifer compared to the standard normal distribution. The random variable was normalized by the sample mean and standard deviation to illustrate the lognormality. where lni; is the sample mean and Ei n is the sample standard deviation. Working here with the logarithm, rather than the velocity itself, better demonstrates the apparent lognormality of the velocity fields. Recall, the logarithm of a lognormal random variable is itself normally distributed. Figure 2.16 indicates that the transformed log initial velocities are welldescribed by a standard normal distribution. The lognormality of the velocities seems logical since the conductivity field is itself lognormally distributed. However, the velocity is in general a product of random variables: the conductivity, the head gradient, and the mobile water content. The logarithm of the velocity is the sum of random variables. Sums of random variables tend toward normality by the central limit theorem (see Mood et al. [1974]), and thus, products of random variables tend towards lognormality. It would be interesting to test this PAGE 55 40 hypothesis by examining the velocity characteristics of flow fields arising from different log conductivity distributions (e.g., fractal, uniform, white noise, etc.), and this is left to future work. 2.4.9 Head and Head Gradient Head and head gradient information along the trajectories was not recorded. This is a regrettable oversight, since much of the behavior of the flow water and the transport of solutes in heterogeneous media can be summarized as the interplay of head, head gradient and conductivity covariance functions. The literature concerning Eulerian head and head gradient fields is extensive and is summarized in large part by Gelhar [1993]. In fact, much is drawn from this body of literature for some of the theoretical developments in the following sections. That this information was not recorded detracts mainly from the completeness of this work, and further investigations are left to future works. 2.5 Theory 2.5.1 Effective Conductivity As with the simulations, the theoretical developments are founded on the bedrock of the log conductivity field. Section 2.4.7 demonstrates that the simulation results well match the "target" ensemble quantities. From this field is had directly a derived parameter, namely the a priori effective conductivity K e An expression for the effective conductivity, similar to that given by Gelhar and Axness [1983], is derived here to illustrate the mechanics of different techniques employed through out these derivations. Consider a onedimensional Darcy's law of the form (2.15) where q is the specific discharge, k is the hydraulic conductivity, and PAGE 56 41 may be resolved into a mean, and a zeromean perturbation, e.g., In A; = fj,\ n k + / Let the head gradient likewise be resolved into a mean and zeromean perturbation, e.g., d(p/dx = J(l + ;'). Inserting these expressions into Equation (2.15) yields q = exp[/i[ n k + f]J{l + j) (2.16) = K g Jexp[f](l+j) Expand the exponential term in this equation as a Taylor series about zero, truncating the series at secondorder, and expand the product of the perturbation terms q~K g J(l + f + f/2)(l+j) (2.17) = K g J(l + f + f 2 /2 + j+fj + f 2 j/2) Take the expected value of this expression, dropping terms higher than secondorder, resulting in an expression < q >= K g J(l + al k /2+ < fj >) (2.18) Consider this equation and the three terms in the parenthesis. The first is a "mean" term, the second is a "contribution" due to the variability of the conductivity field and the third is the log conductivity and head gradient crosscovariance. Gelhar [1993] evaluates this expression for isotropic log conductivity fields using spectral methods. The result, which is independent of the shape of the correlation function, is surprisingly simple [Gelhar, 1993] PAGE 57 42 effective conductivity is simply K e = K g Notice that the effective conductivity can be significantly less than the arithmetic mean conductivity K a = K g exp[af nk /2]. The definition of the effective conductivity is quite general, and a general effective conductivity expression is denoted K e = K g c where c is the "correction" due to the aforementioned processes. Under this definition, there is always some effective conductivity value, regardless of nonstationarities in boundaries or flow domains, but its magnitude may depend upon several factors. There is no general a priori effective conductivity estimator. The estimator given here is limited to stationary head and conductivity fields. 2.5.2 Eulerian Velocity Field The relationship between the Eulerian conductivity field and its associated Eulerian velocity field has been one of the more thoroughly examined aspects of subsurface hydrology. For the example problem, closedform expressions are sought for the Eulerian, or spatial, mean velocity U, the Eulerian velocity covariance C ua [r], and its point distribution, given that the underlying conductivity field is large, stationary, lognormallydistributed, exponentially correlated, and characterized by parameters H\nk, cfnfc' and Af nfc and subject to the conditions prescribed for the model aquifer (see Section 2.1 and Figure 2.1). From the development of the effective conductivity relationship in the preceding section, an expression for the mean Eulerian velocity is K Q cJ U = ~T( 2 2 ) where c = 1 for a twodimensional isotropic aquifer. The c is retained to emphasize that this expression is not the "consistent firstorder approximation" of the velocity that would result from dropping the a? nk and < fj > terms in Equation (2.18). For the twodimensional isotropic aquifer, which is the model case, these two terms happen to be equal in magnitude and opposite in sign, and the consistent firstorder PAGE 58 43 expression and effective expression are the same. This damping effect of the negative head gradient and log hydraulic conductivity is largely responsible for the robustness of these perturbation expansions. The truncated terms are not negligible in the sense their magnitude is small. For instance, for a 2 nk = 1, their magnitude is onehalf that of the mean. This brings us to an important philosophical point. A primary usefulness of these expansions is to identify the dominant "components" of the process. A purely "consistent" approach could lead to neglect of some important process. For instance, consider a strict, firstorder expansion of Equation (2.16). Equation (2.17) would contain no f 2 term, and Equation (2.18) would contain no af nk term. The resulting effective expression would seriously underestimate the specific discharge, since the contribution of the conductivity variability would "carry no weight" with respect to the other contributing processes. Analysis of the neglected higher order terms is necessary for a more complete understanding of the process. The general goal of the approach, however, is to identify the dominant subprocesses that contribute to the larger process, and the relative importance, or weight, of each subprocess. Expressions for the Eulerian velocity covariance for the model system have been given by various researchers, including Graham and McLaughlin [1989a] and Rubin [1990]. Shapiro and Cvetkovic [1988] develop an inverse velocity covariance function using a perturbation expansion of the inverse flow velocity and spectral methods. Consider once again the perturbed Darcy's law expression given by Equation (2.16). Dividing this by the homogeneous mobile water content, neglecting secondorder terms and higher, the meanremoved expression, or "velocity perturbation," is u' = ^(f + j) (2.21) PAGE 59 44 If these processes are stationary stochastic processes, they can be represented represented by a FourierStieltjes integral (see e.g., Bakr et al. [1978]) poo u'= exp[iMx]dZ u [M] (2.22) J OO where i \/T, M is the wave number vector, and dZ u < is the complex random amplitude of the velocity perturbation. The FourierStieltjes representation of Equation (2.21) is dZ u = ^(dZ f + dZj) (2.23) u Compare this expression to that derived by Shapiro and Cvetkovic [1988] (their Equation (A4)) dZ l/v = Yj{dZ f + dZ j ) (2.24) The spectra that result from these complex random amplitudes are structurally the same, differing by the square of the mean values by which they are multiplied (e.g., by (KgJ/9) 2 on the one hand and (6/(K g J)) 2 on the other). In general, the correlation structure of the inverse of a random function is not the same as that of the random function itself. The firstorder series approximation of the inverse velocity is linear in the velocity itself (e.g., l/v ~ 1 v ), and thus the velocity correlation characteristics are recovered as an approximation of the inverse velocity correlation characteristics. For the case of the exponentially correlated isotropic conductivity field, the longitudinal velocity covariance of a twodimensional field is given by (Shapiro and Cvetkovic [1988]) C u [s] = 3al k (K g J/6) 2 (exp[s](s2 + 3*3 + 3s~ 4 ) + s" 2 /2 3s~ 4 ) (2.25) PAGE 60 45 where s = x/\\ n kThe corresponding expression for a threedimensional field is C u [s] = Sal k (K g J/6) 2 (exp[s](s2 + 5s" 3 + 12s" 4 + 12s" 5 ) + s" 3 12s" 5 ) (2.26) Equation (2.25) provides a superior match to the data in comparison to Equation (2.13) (see Figure 2.17). However, the expression for the Eulerian covariance 1 0.8 c o u 0.4 _o > 0.2 0 0 2 4 6 8 10 r/Xlnk Figure 2.17: Comparison of velocity covariance expression Equation (2.25) to the approximation Equation (2.13) for of nfc 1. Covariances are presented as dimensionless correlation functions. given in Equation (2.13) captures the dominant features of the velocity covariance structure, including the correlation length and the point velocity variance of Equation (2.25) in a simple and familiar form. As noted, however, Equation (2.13) exhibits some systematic differences from the observed behavior. The data show a rapid initial decay, followed by a longrange persistence (see Figure 2.17). This behavior is examined with an expression adapted PAGE 61 46 from Dagan [1984] o r i 2 (r, M dR ]nk4> [r] d^r}\ Rua[r] = ^lfc ( R\nk[r\ ^ J l^'J where R\ nk is the log conductivity correlation function, R\ n k PAGE 62 47 approximation, however, will realize significant returns in the form of simplicity, as shown in subsequent sections. Previous works have focused upon the moments of the Eulerian velocity field that result from mean gradients imposed upon heterogeneous conductivity fields. These works shall be extended slightly with a hypothesis about a specific case that allows a variety of insights to be drawn. Hypothesis: The point velocity distribution of an Eulerian velocity field induced by an stationary mean gradient J across a lognormally distributed, exponentially correlated log hydraulic conductivity field characterized by parameters n\ n k, PAGE 63 48 normal. This indicates that the moments n u = U and a\ U 2 baf nk may require some further investigation. The adequacy of these estimates are discussed in subsequent sections. 2.5.3 Streamtube Weights The reference velocity field is now to be subdivided into streamtubes. Each streamtube may be considered to carry some "weight" in relation to other streamtubes for a given collection. Consider the partitioning of the flow field into N streamlines separated by an equal distance along the IP. This is an equalarea weighted partitioning scheme, and results in N subIPs of equal area, but dissimilar discharge. The discharge through subIP i is where Aa = A/N. The relative flowweight of this streamtube i is given by An equalflow partitioning scheme results in in N subIPs of equal discharge, but of dissimilar area. The discharge through subIP i is di = 9AaUi (2.30) Wfi = Ui/U (2.31) d = 9AciiVi (2.32) = Q/N Solving for the area yields (2.33) The relative areaweight of this streamtube i is Wai = U/Vi (2.34) PAGE 64 49 The perceptive reader will notice a slight notational change between the two partitioning schemes. A lowercase u represents velocities drawn uniformly in space and a lowercase v represents velocities drawn uniformly in flow. The distinction is important. It is through these subIPs that water enters the streamtubes and a description of the subIP properties is essential to streamtube, and thus flow field, descriptions. If from a collection of N subIPs the random selection, or drawing, of one subIP is equally likely as any other, the weight assigned to any one observation from M samples would be 1/M. Assume that all of the subIPs are sampled. The average subIP discharge is the same for the equalarea and equalflow partitioning schemes, and equal to d = 8AU/N. However, the discharge variance is zero for equalflow and nonzero for equalarea, when the velocity field is heterogeneous. The weights expressed in Equations (2.34) and (2.31) allow the properties of one collection of streamtubes to be mapped to the other collection of streamtubes. These weights combined with the assumed values of the Eulerian velocity field permit derivation of properties of the initial velocities of the equalflow partitioning schemes. The expected value of v must equal the expected value flowweighted u. Thus, PAGE 65 50 The variance of v is given by var[u] =< v 2 > PAGE 66 51 where r is the lag along the mean flow direction. Again this follows from the requirement that the distribution match the assumed "real space" v moments fi v = V = U(l + ba? nk ) and o 2 v = V 2 ba? nk The correlation structure of the FW/LT velocities is similar to that of AW/ET (see Figure 2.12). Employing these velocity distribution hypotheses, it is now possible to derive some relationships between the different flow field partitioning schemes. From observation, it is inferred that the velocities along areaweighted Euleriantrajectories and flowweighted Lagrangiantrajectories are secondorder stationary in displacements along the mean flow direction (recall Figures 2.10 and 2.11). Having adopted or postulated expressions for these first two moments, it is now possible to derive expressions for the nonstationary velocity means and variances of the areaweighted Lagrangian and the flowweighted Eulerian streamtubes. 2.5.4 Control Plane Oriented Velocity Moments Cvetkovic et al. [1996] postulated a semianalytical estimator for a nonstationary mean velocity along areaweighted Lagrangiantrajectory observed at control planes. We would like to work directly from the assumed stationary flowweighted Lagrangiantrajectory velocity field. The general expression for the expected value of the areaweighted Lagrangiantrajectory velocity is Resolve the stationary mean velocities into the mean of the flowweighted Lagrangiantrajectory field and zeromean perturbation. (2.41) v a [x] =< U V{l + v'[x})> (2.42) V(l+v'[0]) PAGE 67 52 Expand the denominator in a series to second order v a [x] =< U{1 v'[0] + v'[0] 2 ){l + v'[x}) > = U(l + & ^) (243) = U(l + ba 2 nk (lexp[br/X lnk })) In a similar fashion, an expression is derived for the nonstationary mean velocity along a flowweighted Euleriantrajectory. In this case, it is the stationary Euleriantrajectory velocity field that serves as a basis. The mean of the perturbation expansion is the Eulerian mean velocity. It is interesting to note that this result requires no series expansion. u f [x] =<^u[x}> = = 17(1 + £=M) = U(l + ba 2 nk exp[br/\ lnk ]) The nonstationary velocity variances along the respective crossweighted trajectories are derived in a similar fashion. The general expression for the areaweighted Lagrangiantrajectory velocity variance is var[u a [x]] =<^(v[x] v a [x]) 2 > (2.45) v[0\ This equation may be rewritten v&r[v a [x}] =< ^v[x] 2 > v a [x} 2 (2.46) v[0\ Resolve the first term of the preceding equation into perturbation expressions PAGE 68 53 Expanding the denominator in a series, taking the expected value, and neglecting thirdorder and higher moments yields < ^Axf >= UV(1 + 2ba 2 nk {\ + exp] bx/X lak })) (2.48) v[0\ In order to improve the quality of the estimate, the righthand side (RHS) of the preceding equation is required to match the known extreme values of the lefthand (LHS) side. The LHS at X = 0 is <T7^[0] 2 > = U PAGE 69 54 Substituting Equations (2.53) and (2.43) into Equation (2.45) yields vax[v a [x]] = UV(l + (2bal k + (kr 2 nfc ) 2 )(l exp[6r/A lnfc ])) (2.54) (U(l + ba? nk (lexp[br/\ lnk }))) 2 Further expansion of this expression offers no intuitively obvious simplifications or insights in its entirety. However, this expression may be approximated by dropping a few terms thus varK[x]] = V Vn* U 2 bal k {2bal k + (6a 2 nfc ) 2 ) exp[6r/A lnfc ] (2.55) Two unconventional techniques were employed in this derivation that shall now be defended. The first is requiring the series expansion to match the known endpoints. Simple expressions that condense the system behavior into an compact form that allows the relative contribution of the component processes to be easily grasped are sought. The series expansion and subsequent truncation is itself an approximation. The systematic error introduced by "missing" a known endpoint is distracting. In the case of the velocity covariance, the shape of the correlation function was "sacrificed" for simplicity, but the correlation length was retained. Had the simplified expression neither matched the shape of the data nor reproduced the correlation length, the expression would have been somewhat more questionable. In this case, the expression is required meet the observed endpoint values at the "cost" of the rigor of the series expansion and truncation. The second simplification is the ad hoc dropping of terms to simplify the expression. Again, simple expressions that allow the relative contribution of the different processes be easily grasped are sought. It would be difficult to hypothesize Equation (2.55) a priori as an empirical expression. A primary function of the series expansion is to indicate the dominant processes and give an estimate of their relative contributions. This hybrid approach is useful for the purposes here. These techniques PAGE 70 55 are not a panacea, however, and the approximations have only been tested for the log conductivity variances tested in this work. The nonstationary flowweighted Euleriantrajectory velocity variance may be derived in a similar fashion. These results are presented directly. The complete expression is var [u, [*]] = t/ 2 (l + ba 2 Xnk + [2baf nk + 3 {ba 2 nk ) 2 + (ftojU) 3 ) (1 exp [6r/A lnfc ])) (L/(l + 6a 2 n ,exp[6r/A ln ,])) 2 (2.56) and an approximation of the above obtained by neglecting a few terms is var [u f [x\] = U 2 bo? nk (l + (2bo 2 nk + {bal k ) 2 ) exp [6r/A lnfc ]) (2.57) 2.5.5 Control Time Oriented Velocity Moments The displacementinvariant velocity covariances for u and v permitted derivation of expressions for the first two moments of the respective velocities for the four different trajectory collections. However, only the areaweighted Lagrangiantrajectory collection appears to be stationary in time. This is perhaps not intuitively obvious, as the velocity is clearly not stationary in displacements along the mean direction of flow (e.g., Equations (2.43) and (2.54)) because there is an implicit distribution of times associated with the establishment of a reference injection plane. After long travel distances and times, the areaweighted Lagrangiantrajectory point velocity statistics are equal to the flowweighted Lagrangiantrajectory point statistics. At the injection plane, the areaweighted Lagrangiantrajectory point velocity statistics are identically equal to the areaweighted Euleriantrajectory point velocity statistics. These relationships can be exploited to derive equations that provide a few insights about how properties change in Lagrangian time, and to derive an expression for a Lagrangiantime dependent velocity covariance function. PAGE 71 56 The areaweighted velocity at the IP is U by definition. Since the areaweighted Lagrangiantrajectory velocity is stationary in time, its mean in time is U regardless of position. At long travel times, the point velocity is stationary in displacements as well. For constant means in space and time, the following relationship holds ?/>(*jr iif (2 58) This expression states that the harmonic mean velocity in space equals the arithmetic mean velocity in time. Returning now to the model aquifer, invoke the hypothesis of lognormality of the velocities along either trajectory. The harmonic mean of v[x] is given by V h = exp[/ii n (Tint,/ 2] = exp[/ii n ] exp[<7i n ] 1/2 1 (2.59) = U Similarly, the harmonic mean of u[x] is given by Uh = 1^2 ( 2 6 ) Although derived by different means, this result is equal to the result given by Shapiro and Cvetkovic [1988]. These harmonic means are the asymptotic mean velocities of the two trajectories in time. The initial timedependent velocities are identical to the initial spacedependent velocities, an artifact of starting all of the trajectories at the injection plane at an implied time t = 0. This indicates that the velocity covariance in time along the areaweighted Lagrangian trajectory is C va [t = 0] = U 2 bo~\ nk As mentioned, at large times, the statistical properties of the Lagrangiantrajectories converge. Making a change of variables x = Vht Ut'm Equation (2.40) and substituting the known C va [t = 0] for the value of C v f[r = 0] in yields an PAGE 72 57 expression for the stationary areaweighted Lagrangiantrajectory velocity covariance C va [t] = U 2 bal k exp[bUt/X lnk ] (2.61) Similar expressions have been derived before by other researchers (see e.g., Dagan [1989]), but the conceptual development is quite different. Rather than seek to approximate the covariance function, an approach to quantify it directly from properties of the system is demonstrated. Moreover, the quantity U is perceived to be the harmonic mean of the Lagrangiantrajectory velocity, and that it is quantitatively equal to the areaaverage velocity of the flow field. This is an important point. The mean velocity and velocity variance along the flowweighted Lagrangiantrajectory is found in a manner similar to those in space, but there are some conceptual issues that should be addressed. It is the collection of areaweighted Lagrangiantrajectories that exhibits a stationary velocity. The mean velocity is Vft[t]=<^v[t}> (2.62) Expand the velocities into a mean and a zeromean perturbation. The mean here is the time average mean, or the harmonic mean of the space average. /<[<] = ^ = U (1 + bal k exp[bUt/\ lnk }) This expression requires no series expansion. For the control plane oriented velocity statistics, the mapping was from flowweighted to areaweighted trajectories. Equation 2.63 is identical in form to Equation 2.44 with the substitution r = Ut. Recognize PAGE 73 58 the quantitative similarity and write velocity variance equations directly The complete expression is var [v f [t]} = U 2 (l + ba 2 nk + (26a, 2 n + 3 (ba 2 nk ) 2 + {ba 2 nk ) 3 ) (1 exp [bUt/ X lnk })) (U (l + ba 2 nk ex V [br/\ lnk \)) 2 (2.64) and the approximation is var[v/[t]] = U 2 baf nk (l + [2ba 2 nk + {ba 2 ak ) 2 ) exp[bUt/X lnk }) (2.65) Calculation of the mean velocities along the Eulerian trajectories is problematic in the absence of a stationary velocity covariance function. However, note that the areaweighted velocities start at U and that the flowweighted velocities start at V and both progress to U/(l + baf nk ). Recognizing the similarity shared by the different mean time estimators, write approximate expressions directly uat [t] = ,Ji 2 (1 + ba 2 nk exp [bUt/X {nk ]) (2.66) and W = T^TTI 1 + (2^infc + *lk) (^P [Mtt/Ain*])) (267) 1 + 0<7 ln fc Again, devoid of a stationary covariance function for the Eulerian velocity field in time, intuitive arguments lead us to anticipate an asymptotic u variance in time to be 2 t r2 l 2 a ut = U h ba \nk ( U \\ 2 ( 2 68 ) The initial velocity variances are known from the covariance functions in space. It is possible to anticipate forms of the nonstationary variance equations similar to those found previously, and include them for the sake of symmetry and completeness. The PAGE 74 59 areaweighted Euleriantrajectory timedependent velocity variance approximation is varK[i]] = (y^J Wnfc (1 + {^< k + {b PAGE 75 60 Finally, note that the AW/LT log conductivity appears stationary in time as does the velocity. Moreover, it appears to have values equal in magnitude to its AW/ET counterpart, as does the velocity. Unlike the velocity, expressions for the trajectorybased crosscovariance of the velocity and log hydraulic conductivity that are necessary for developing the nonstationary log hydraulic conductivity moments were not developed. This is left for future work. 2.6 Evaluation of Simulation and Theoretical Results 2.6.1 Mean Velocities in Space Figure 2.10 presents the control plane oriented mean velocities for the four different collections of streamtubes. The mean velocities are normalized by U = K g J/6 = 1/10 and this value was used in the mean velocity expressions as well. That is, the parameter values are the input parameters to the simulations. Other input parameter values are 6 = 3/8 and of nfc = {1/2, 1} depending upon the set of simulations. The deviation of the mean of areaweighted Euleriantrajectory data from 1 are an indication of the quality of the effective conductivity expression. The areaweighted Euleriantrajectory mean velocity appears stationary, in accordance with established theory, and its value appears to be well predicted by the effective conductivity relationship. The fluxweighted Lagrangiantrajectory mean velocity appears stationary, supporting the stationarity hypothesis. Moreover, the value of the mean velocity appears to be well predicted for both variances of log conductivity. Additionally, the nonstationary mean velocity estimators appear to accurately reproduce the transition between the asymptotic velocities. 2.6.2 Velocity Variances in Space Figure 2.11 presents the control plane oriented velocity variances for the four different collections of streamtubes. The variances are normalized by the areaweighted PAGE 76 61 Euleriantrajectory velocity variance estimated from the input parameters to the simulations. That is, by ^ = (V/)V hi (274) As with the mean velocity estimators, the variance estimators employ only the simulation input parameters. Once again, the areaweighted Euleriantrajectory variance exhibited its presumed behavior, and supports the well established fact that the Eulerian velocity field is fairly well understood in space. The hypothesized flowweighted Lagrangiantrajectory velocity variance is somewhat higher that the data for both af nk A possible reason is the negative correlation between log conductivity and head gradient is not weighted properly in the collection of trajectories skewed towards higher log conductivity values. The other two estimators also suffer from this bias in the asymptotic Lagrangiantrajectory velocity variance, although the transition between initial and final values for both seems to be well characterized. 2.6.3 Mean Velocities in Time Figure 2.13 presents the timedependent mean velocities for the four different collections of streamtubes. Again, the mean velocities are normalized by U as estimated from the simulation input parameters. Only the areaweighted Lagrangiantrajectory mean velocity exhibited a constant mean behavior in time. Its estimated harmonic mean is U, and the estimate seems good, although one might argue that there is a slight overpredictive bias. The estimated harmonic mean of Eulerian trajectories appear to have a stronger under predictive bias. The transition of the flowweighted Lagrangiantrajectory appears to be captured better than those of the Eulerian trajectories. It should be noted that the Lagrangian trajectory transition was based upon an expression for PAGE 77 62 a stationary velocity field, whereas the Eulerian trajectories are ad hoc estimators, based upon a transition between known end points. 2.6.4 Velocity Variances in Time Figure 2.14 presents the timedependent velocity variances for the four different collections of streamtubes. The variances are normalized by the areaweighted Lagrangiantrajectory velocity variance estimated from the input parameters to the simulations. The areaweighted Lagrangianvelocity appears to be stationary in the variance. The data seem to lie slightly above the estimated value. The overprediction of the initial fluxweighted velocity variance has far less impact upon the timedependent velocity variances. This combined with what appears to be a good estimate of the longtime Euleriantrajectory velocity variance yields a better match to the simulations than one might expect from the means. However, there is a strong initial decline in the variance not captured by the simple exponential expressions. The areaweighted Lagrangiantrajectory covariance function exhibits a similar rapid decay. The early time behavior is likely dominated by the conductivity covariance. As the trajectory experiences a wide range of conductivity variability, this effect diminishes, and the longer range correlations related induced by the head field begin to dominate. 2.6.5 Log Conductivity Statistics The predictors of the stationary log conductivity statistics worked well for both the spacedependent and timedependent means although there seemed to be a slight overpredictive bias for the flowweighted trajectory (see Figures 2.6 and2.8). The AW/ET log conductivity variance was "well predicted" in that the turningbands method performed well, as discussed in Sections 2.4.1 and 2.4.7. The log conductivity variance of the FW/LT appeared to be lower than that of the AW/ET, as predicted by PAGE 78 63 Equation (2.73), but the magnitude of under prediction was larger than that observed (see Table 2.4). Table 2.4: Comparison of observed log hydraulic conductivity variances to estimators. Observed value is zeroslope line regressed through data. The decimal representation of the fractional value is included for convenience. < k = 1/2 parameter predicted observed AW/ETai 2 nJt [x] 1/2 0497" FW/LTcrj^x] 7/16=0.4375 0.484 AW/LTajUM 1/2 0.498 parameter predicted observed M/ETal k [x] 1 L002~ FW/LTaf nk [x] 3/4 0.944 AW/LTagJi] 1 1.000 2.7 General Conclusions In an steady irrotational flow field, a fluid parcel trajectory follows a hydrodynamic streamline. The flow field properties as observed along streamlines can differ radically than those observed "uniformly in space." For the model system, a stationary Eulerian log conductivity field, this work identified the trajectory collections for which certain properties appear statistically stationary to second order. For displacements along the mean flow direction, point velocity and log conductivity appear stationary along streamlines separated by equal discharges. For displacements in travel time, these same properties appear stationary along streamlines that are selected uniformly in space. The stationarity of these properties facilitates the development of simple expressions for the statistics of displacement and travel times, as we shall show in following chapters. An immediate result of the preceding analysis is that the bulk of the flow field in the sense of total discharge travels at a velocity that is locally greater than the spatial PAGE 79 64 average, and this disparity increases with heterogeneity of the underlying conductivity field. This might be described as "preferential flow." A practical implication for the transport of kinetically sorbing or interacting solutes is that the local velocity for the majority of the flow may occur in a regime significantly higher than that estimated as the flow field average. An understanding of the medium heterogeneity is essential if laboratory experiments are to be "scaled up." It is intuitively obvious, and has certainly been noted before, that the bulk of the flow passes through an increasingly small portion of the total "swept volume" as the medium heterogeneity increases. Again, the practical implications are likewise obvious. Consider some aquifer remediation effort that involves flushing a "contaminated zone" with some fluid. This might be as simple as the upgradient water in the case of pumpandtreat or as extravagant as an surfactantenhanced microemulsion. In heterogeneous flow systems, the bulk of the flushing fluid contacts a disproportionately small volume of the total swept volume. For contaminants that tend to reside in areas of high conductivity, the serves to enhance the efficacy of remediation efforts. Contaminants that tend to reside in areas of low permeability, conversely, would be more difficult to remove. While these conclusions may be drawn from "common sense," application of the concepts developed here may help to better understand these issues quantitatively. PAGE 80 CHAPTER 3 NONREACTIVE SOLUTE TRANSPORT 3.1 Introduction The previous chapter saw the derivation of the relationship between the velocities that lie along Lagrangian trajectories and the parameters that describe a simple model Eulerian flow field. These Lagrangian relationships are employed to derive a few functions that characterize global measures of solute transport: namely the statistics of mass displacement in time and the statistics of mass arrival times at control planes. The frameworks are are that of Dagan [1982b] and that of Shapiro and Cvetkovic [1988], respectively. The recovery of a few wellknown results are anticipated, in addition to some that are new. However, the approach is novel in that the work is based upon trajectorybased statistical properties of the flow field. 3.2 Injection Mode and Streamtube Selection Demmy et al. [1999] noted the relationship between boundary conditions and sampling of streamlines. A clarification of this relationship and a generalization the results of that work are sought. Unlike the previous chapter, the Eulerian trajectories are not considered in the context of an actual transport trajectory, as they are of little practical importance. However, the results of Shapiro and Cvetkovic [1988] correspond to areaweighted Euleriantrajectories, as will be shown. Therefore, trajectory, Lagrangian trajectory, streamtube, and streamline are used interchangeably. 3.2.1 Uniform Resident Injection Consider the application of mass at a constant density c 0 across the entire face of the injection plane (IP) of the model aquifer (see Figure 2.1). Let the mass occupy 65 PAGE 81 66 some small Ax that is invariant with displacements along the IP. Following the terminology of Kreft and Zuber [1978], this injection mode is defined to be a uniform resident injection. Uniform refers to the homogeneous mass density in space and resident refers to the volume of resident fluid into which the solute is introduced. Subsequent references to a resident injection will refer to a uniform resident injection unless otherwise noted. Now, this injection mode has nothing to do with the flow field per se. That is, the mass "magically appears" and occupies a fixed volume at a fixed concentration without being influenced by local water fluxes. Consider the flow field to be divided into some number of streamtubes, say N. Moreover, let these streamtubes carry an equal areaweight at the IP. For a small depth of injection Ax, this injection mode has assigned an equal mass weight to each streamtube with respect to the total solute mass. The total mass of solute is M r = c Q 6hLAx (3.1) The mass weight for any one streamtube i is ^ = 1 (3.2) M N v This value is identically equal to the areaweight of the streamtube, and constant for all in the collection. Consider now dividing the streamtube into N streamtubes that all carry the same flow weight. The mass in some streamtube i is given by = b (3.3) PAGE 82 67 3.2.2 Injection in Flux Consider now maintaining the IP at c 0 for some brief interval At. A mass Mj = c 0 8hLUAt (3.4) enters the domain and the different streamtube collections. This injection mode is a uniform injection in flux. Uniform, again, refers to the uniform mass density and "in flux" refers to the influent water that carries the solute into the flow domain. A mass rrii = c o 0h^ViAt (3.5) enters equalarea streamtube i. This corresponds to a mass weight of A mass rrii = coOhciiViAt (3.6) Q (37) = c 0 9h^At enters equalflow streamtube i. This corresponds to a mass weight of 1/N. 3.3 Integrated Streamtube Properties Let us return to the notion of the streamtubes as pure hydrodynamic entities. It is possible to discuss properties of these streamtubes with no mention of injection modes. Two properties hold some hydrodynamic interest. The first is the velocity integrated over time. For for some streamtube i this is Xi[t]= [ Vi[t]dt (3.8) Jo A velocity integrated over time is a "displacement" with a dimension of length. The second the inverse velocity integrated over some distance. For some streamtube i this PAGE 83 68 is (3.9) This integral has a dimension of time. These integrals are given in terms of the x velocity component as a function of longitudinal displacement or time. It is suggested that future works examine the details of a more general approach, e.g., displacements along the trajectory. 3.3.1 Displacements Consider a model flow field such as the one given in the previous chapter. Quantification of the displacement and travel time statistics for areaweighted and flowweighted trajectories is sought. Previous results indicate that the velocity for a collection of areaweighted streamtubes is stationary in time, characterized by mean U and approximate covariance U 2 baf nk exp[ bUt/X\ n k\The expected value of Equation (3.8) is Although the velocity of the flowweighted streamtube collection is not stationary in time, the results of the previous chapter are employed to map the properties of the areaweighted collection to the flowweighted collection. Recognize that this equation may be rearranged to be the integral of Equation (2.63) in time. Substitute the RHS of Equation (2.63) into the preceding equation and PAGE 84 69 integrate to yield PAGE 85 70 The arrival time variance for the flowweighted collection of streamtubes is Jo v[x'\ Jo v[x"\ Approximate the inverse velocity covariance contained in the first term of the RHS of Equation 3.16 with the a series expansion of the perturbation equation 1 1 > = PAGE 86 71 one particle into each areaweighted streamtube at x = 0 and t = 0 is analogous to the a uniform resident mass injection. Equations (3.10) and (3.13) express the mean particle displacement and particle displacement variance as a function of time. Equation (3.15) expresses the mean particle travel time at different control planes. Similarly, placing one particle into each flowweighted streamtube at x = 0 and t = 0 is analogous to the a uniform injection in flux. Equation (3.12) expresses the mean particle displacement as a function of time. Equations (3.14) and (3.18) express the mean particle travel time and particle travel time variance at different control planes. In addition to these two injection modes, there are hybrid injection modes. For example, consider if the introduction were somewhat random, or particle mass were some deterministic function of the injection velocity, or some combination. In fact, it is left to the imaginative reader to conceive of different injection mode possibilities and the physical systems to which these correspond. 3.5 Evaluation of Expressions In the previous chapter, a numerical experiment used to test certain hypothesis about the relationship between Lagrangian and Eulerian descriptions of flow fields was described. Those data are employed here. In the comparison of the derived expressions to the data, the functions are based upon the simulation input parameters, and all normalization is carried out using these same parameters. Therefore, the expressions represent a priori estimates of the simulation results, based, of course, on known input parameters. 3.5.1 Displacements Figure 3.1 presents the model estimates of mean mass displacements for the two injection modes. The flux injection preferentially selects, or mass weights, high velocity areas of the flow domain, and the propagation of the center of mass is characterized PAGE 87 72 by a higher velocity. The timevarying mean velocity expression for the flowweighted Lagrangiantrajectories (Equation (2.63)) is, in fact, an expression for the mean velocity of a plume resulting from a uniform flux injection. The longtime asymptotic slope of the fluxweighted mean displacement is U, as is easily verified by taking limit of Equation (2.63) as t > oo. The mean displacement functions predict the observed behavior, and also provide an excellent fit to the data. The displacement variance data are interesting (see Figure 3.2). Injection mode seems to have little effect upon the displacement variance in time. This is not to say that injection mode has little effect upon the spatial distribution of the plume in time. At a given time t, a plume injected in flux will have traveled a greater distance, on average. After traveling the same distance, the displacement variance of the fluxinjected plume will be less than that of the resident injection. Thus, the "spreading as function of mean displacement" is less in the case of the fluxinjected solute. The areaweighted trajectory displacement variance overpredicts the observed displacement variances, and the longtime slopes appear to have different values. After 10 characteristic times, the exponential term in Equation (3.13) retains a little more that two percent of its original value. This difference does not seem large enough to account for the observed discrepancy. The variance of the stationary areaweighted Lagrangiantrajectory velocity field in time appears to be well predicted. Thus, the problem may lie with the rather poor description of the shorttime and longtime correlation exhibited by the observed velocity covariance given by the assumed model. 3.5.2 Travel Times Figure 3.3 shows good agreement between the estimated and observed mean travel times. The injection influx shows earlier arrival than that of the uniform resident injection. Uniform resident injection weights all areas of the IP equally. Large areas that contribute little to the overall flow receive and equal mass weighting as those PAGE 88 73 10 AW/LT FW/LT Eq. 3.12 Eq. 3.10 X* ..+ .X .+ .x ..+ .X .+ .X 4* X .+ x + .* .+ .X .x .X ,+ 4 6 tu/x lnk 10 Figure 3.1: Comparison of mean displacements for the two injection modes or streamtube collections. The flowweighted Lagrangian trajectory (FW/LT) statistics correspond to a flux injection, and the areaweighed Lagrangian trajectory (AW/LT) statistics correspond to a uniform resident injection. The top figure contains the a ink = 1/2 data and the bottom figure contains the af nk = 1 data. PAGE 89 74 15 AW/LT + FW/LT x Eq. 3.13 6 b 10 0* 4 6 tU/\ lnk 10 15 AW/LT + FW/LT x Eq. 3.13 10 +x +x +x +x +x +x +x +x +x t CS H b / ** 4 6 tu/x lnk 10 Figure 3.2: Comparison of displacement variances for the two injection modes or streamtube collections. The flowweighted Lagrangian trajectory (FW/LT) statistics correspond to a flux injection, and the areaweighed Lagrangian trajectory (AW/LT) statistics correspond to a uniform resident injection. The top figure contains the fnfc = V2 data and the bottom figure contains the of nA = 1 data. Variances are normalized by Xi k af nk PAGE 90 75 areas that contribute more significantly to the overall flow. However, as the solute moves through the domain and experiences the full range of velocity variability along its trajectory, the "slow start" is damped out and the instantaneous average point velocity for all for the particles is V and the harmonic mean V h = U. The travel time variances appear to be better matched than the displacement variances (see 3.4). This is reasonable in light of the better apparent match of the assumed velocity covariance function in space than that in time. The uniform injection appears to exhibit less nonlinear behavior than the injection in flux. The uniform injection variance is reasonably expected to be larger than that of the flux injection, due to the skewed massweighting of lowvelocity areas by the uniform resident injection. 3.6 General Conclusions It is clear that the injection mode effects are most prevalent in the "near field" where the initial velocities and travel times or displacements are well correlated. After the plume travels several correlation scales, these injection mode effects diminish. The near field and the neartofar transitional scales are quite interesting, as these intermediate scales are on the order of sites associated with certain agricultural settings such as some feed lots and waste holding facilities, urban land uses such as gas stations and dry cleaners, and their associated remediation efforts, and a plethora of "field scale" experiments. These scales are examined in the perspective of the "error" in travel time estimate associated with assuming a resident injection when the "real" injection is in flux e[x] = Xi nk af nk (l exp[bx/\ lnk ])/x (3.19) given by subtracting Equation (3.15) by Equation (3.14) and dividing the difference by Equation (3.14) (see Figure 3.5) The "near field" might be characterized as the first five to ten correlation lengths, where the injection mode effects are quite prevalent, even for modest log conductivity variability. These effects diminish less rapidly over PAGE 91 76 40 30 J? AW/LT + FW/LT x Eq. 3.15 Eq. 3.14 c 20 10 10 40 20 30 40 30 20 10 AW/LT FW/LT Eq. 3.15 Eq. 3.14 J210 20 30 40 Figure 3.3: Comparison of mean arrival times for the two injection modes or streamtube collections. The flowweighted Lagrangian trajectory (FW/LT) statistics correspond to a flux injection, and the areaweighed Lagrangian trajectory (AW/LT) statistics correspond to a uniform resident injection. The top figure contains the a \nk = 1/2 data and the bottom figure contains the o^ nk = 1 data. PAGE 92 77 80 60 Ji 6 £ 40 AW/LT + FW/LT x Eq. 3.18 tx .+>< fx .fx ,4x + x ...fx ., + x b 20 + x 10 80 20 Z/Alnfc 30 40 + x .. + X .+ xy 60 a 6 li 5 40 Hi b AW/LT FW/LT x Eq. 3.18 + x + x . + x/ + x.^ *' + x.+ xy + x20 + x 5** 10 20 X/Xlnk 30 40 Figure 3.4: Comparison of arrival time variances for the two injection modes or streamtube collections. The flowweighted Lagrangian trajectory (FW/LT) statistics correspond to a flux injection, and the areaweighed Lagrangian trajectory (AW/LT) statistics correspond to a uniform resident injection. The top figure contains the a ink = 1/2 data and the bottom figure contains the PAGE 93 78 0 i 1 1 J 0 10 20 30 40 Figure 3.5: Travel time estimate error as function of displacement associated with assuming a uniform resident injection when injection is in flux for 0f nfc = {1/2, 1}. a "transition" region between 10 and 30 correlation lengths. Beyond 30 correlation lengths, the difference between the two modes is less than a few percent. Much of the Lagrangian transport literature is predicated upon a uniform resident injection, but the socalled consistent firstorder travel time estimate that is regularly employed is quantitatively that for a flux injection. For large displacements, this error is small, and the simplification afforded by this approximation is considerable (compare Equation (3.14) and Equation (3.15)). For many practical applications, an error of ten, or even more, percent might be considered negligible. Therefore, it is not suggested that this approximation be abandoned, but rather phrasing such as "inconsistent" firstorder approximation be applied in future works, since this Lagrangian expression is inconsistently based upon the Eulerian mean velocity. PAGE 94 CHAPTER 4 REACTIVE SOLUTE TRANSPORT 4.1 Introduction The concepts developed in the previous chapters are employed to analyze reactive solute transport. The focus is upon global measures of transport, rather than predictions of local quantities. Particular processes of interest are injection mode effects, the effects of correlation of the local "reaction parameter" to the log conductivity field, and the interplay of these injection mode and correlation effects. A suite of numerical experiments illustrates these effects. Additionally, expressions are developed for temporal moments of mass breakthrough curves for a solute subject to linear equilibrium sorption for the different injection modes that reflect the observed behavior in relatively simple closedforms. Spatial moments of sorbing solutes are not directly treated in this work. The experimental design was targeted directly towards a temporal moment analysis, and this design decision has some implications for spatial moment analysis. These decision and implications are discussed, in hopes that future works might be guided by this analysis. Some of the properties of the travel time dependent reaction parameter are discussed, as are some of the implications for solute transport. 4.2 Temporal Moments The concepts developed in the preceding chapters are applied to examine the effect of injection mode upon the temporal moments of the breakthrough curve of a solute subject to linear equilibrium sorption. A uniform mean gradient induces a steady irrotational flow in a uniformly saturated porous medium characterized by an exponentially correlated lognormal conductivity distribution. This flow field is divided 79 PAGE 95 80 into collections of noninteracting streamtubes. Dispersion at the streamtube scale is considered to be negligible in its effect on the aggregatescale temporal moments (e.g., Dagan [1989]). The sorption process is related to a generic reaction parameter that is heterogeneous and may be correlated to the log conductivity field. The previous chapters demonstrated the relationship between streamtube selection strategies (equalarea or equalflow) and injection mode (uniform resident and influx). The reference collection is equalflow streamtubes, since the associated velocity and log conductivity statistics are stationary in displacements along the mean flow direction. The primary objective of this chapter is to examine the effect of injection mode upon the temporal moments of the mass breakthrough curve of a solute subject to linear equilibrium sorption in a heterogeneous flow and sorption field. A brief review the conceptual model of solute transport is given. General expressions are presented for the first two temporal moments of the breakthrough curves associated with the uniform injection of a solute both in the influent fluid flux and in the resident fluid. The results of the reactive solute transport simulations are presented, followed by the development of some analytical expressions that help explain the observed behavior. 4.3 Conceptual Overview of Solute Transport Consider a general multidimensional flow field that is resolved into a collection of noninteracting streamtubes. Transport along an individual streamtube is characterized with a transfer function 7[£; x] (see Jury and Roth [1990] and Cvetkovic et al. [1998]). For a conservative and nonreactive solute, the transfer function is simply j[t;x] = S[tT[x}] (4.1) where 8 is the Dirac delta function and r is the travel time. The physical interpretation of this transfer function is this: a particle located at x = 0 "arrives" at x X at time t = r[X}. We conceive of this travel time as not only the time required to PAGE 96 81 travel to some point x, but as an intrinsic integrated streamtube property, namely the inverse velocity integrated along the trajectory (e.g., Equation (3.9)). Recall reaction flow path concept of Cvetkovic et al. [1998]. The reaction flow path, parameterized in space, is defined as where P[x] is some reaction parameter of interest. For linear equilibrium sorption, P is a linear partitioning coefficient. The reaction flow path is also an integrated streamtube property. As such, any reference to /x properties in general are broadly applicable to any distributed reaction parameter, perhaps a firstorder decay coefficient or a Langmuir sorption parameter. However, the focus here is upon linear equilibrium sorption to promote a ready and intuitive grasp of the pertinent concepts. For this case, /i may be thought of as the time that the solute is sorbed. Thus the arrival time of a particle undergoing a linear equilibrium sorption process is the sum of the nonreactive travel time r and the time sorbed [i. The transfer function characterizing this transport and reaction process is (see Cvetkovic et al. [1998]) j[t ] x} = S[tr[x}fi[x}} (4.3) Following the work of Cvetkovic et al. [1998], the adopted reaction parameter model is P = P g exp[w] exp[/? In k] (4.4) where P g is the geometric mean of P, w is a zeromean, exponentiallycorrelated normal random variable with variance cr^ and correlation length X w = X\ n k that is uncorrelated to the log conductivity field, and )3 is a strength of correlation parameter relating P to the underlying log conductivity field. The equality X w = Ai nfc was specified for convenience. PAGE 97 82 Two specific injection modes are considered: uniform injection influx and uniform resident injection. For the uniform injection influx, each streamtube receives an equal amount of solute mass (Equation (3.5)), and thus mass weight (Equation (3.6)). Uniform resident injection distributes the mass uniformly in space, but each streamtube receives a different mass and mass weight (Equation (3.3)). The mass breakthrough for some streamtube i is given by m i [t;x]=m i [0;0)j[tx] (4.5) An ergodic condition is assumed such that aggregate properties may be thought of as ensemble properties in the statistical sense. The aggregate mass breakthrough is given by taking the expected value of Equation (4.5). For injection influx, the expected solute breakthrough curve is given by PAGE 98 83 Two temporal moments of interest for the influx breakthrough curve are the mean t lf = PAGE 99 84 4.4.1 Reaction Parameter From these data, trajectorybased P fields were constructed using Equation (4.4). In this manner, any variety of correlation relationships may be analyzed using the same set of simulation data. There is a dearth of information about the structure of sorption fields in natural systems (see, e.g. Tompson [1993] or Jawitz [1999]). Thus, numerically convenient parameter values were selected. This suite of experiments examined correlation parameter values of /3 = {1, 1/2, 0, 1/2, 1} in log conductivity fields of variances of nJt = {1/2, 1}. The extreme correlation parameters {3 = {1,1} are arguably physically implausible, yet they are included to demonstrate the effect of correlation without the "noise" associated with an uncorrelated component. The remaining three values, f) = {1/2, 0, 1/2} might be considered a range of more plausible values. For consistency, the variance of point log reaction parameter was held constant and equal to that of the log conductivity field by requiring a 2 w = (1 /3 2 )of nA .. The expected value of P was 1. This is a slightly different approach than that taken with the log hydraulic conductivity, where the geometric mean was fixed at K g = 1. In the case of the twodimensional isotropic log conductivity field, the effective conductivity is equal to the geometric mean. The expected value of the reaction value was fixed = 1, because in some practical applications P is assumed to be linearly related to some property of interest, and the "amount" of that property is to be estimated based upon the temporal moments of linearly sorbing solutes, or tracers (see e.g., Jin et al. [1995], Annable et al. [1998], and Jawitz [1999]). The point average of reaction parameter P as observed along a Lagrangian trajectory is a function of injection mode, correlation to the log conductivity and the variance of the log conductivity. As with the log conductivity, the equalflux streamtubes exhibit an apparently stationary mean reaction parameter in displacement along the mean flow direction, whereas the equalarea streamtubes exhibit distinctly nonstationary behavior that is analogous to that of the log conductivity (see Figures 4.1 and
