UFDC Home  |  Search all Groups  |  UF Institutional Repository  |  UF Institutional Repository  |  UF Theses & Dissertations  |  Vendor Digitized Files |   Help

# Solute transport in heterogeneous porous formations

## Material Information

Title:
Solute transport in heterogeneous porous formations
Physical Description:
xv, 138 leaves : ; 29 cm.
Language:
English
Creator:
Demmy, George Gary, 1966-
Publication Date:

## Subjects

Subjects / Keywords:
Agricultural and Biological Engineering thesis, Ph. D   ( lcsh )
Dissertations, Academic -- Agricultural and Biological Engineering -- UF   ( lcsh )
Genre:
bibliography   ( marcgt )
non-fiction   ( marcgt )

## Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 1999.
Bibliography:
Includes bibliographical references (leaves 134-137).
General Note:
Printout.
General Note:
Vita.
Statement of Responsibility:
by George Gary Demmy, Jr.

## Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Resource Identifier:
aleph - 021555660
oclc - 43687467
System ID:
AA00026382:00001

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 fill-in-the-blank 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

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 LAGRANGIAN-EULERIAN 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.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 Path-Travel 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
area-weighted Eulerian trajectories. ................ 26

2.2 Comparison of target log conductivity values for area-weighted Eulerian
trajectories in space and observed area-weighted 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 zero-slope 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

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) equal-area streamline, C) equal-flow streamline,
D) equal-flow 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 turning-bands 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 turning-bands 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 position-dependent 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 position-dependent 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 time-dependent 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 time-dependent 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 two-dimensional 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 flow-weighted 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 flow-weighted 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 flow-weighted Lagrangian trajectory (FW/LT) statistics correspond to a flux injection, and the area-weighed 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 flow-weighted 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 area-weighted Lagrangian (AL), flow-weighted
Lagrangian (AF) and area-weighted 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 area-weighted Lagrangian (AL), flow-weighted
Lagrangian (AF) and area-weighted 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 area-weighted 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 area-weighted 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 flow-weighted 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 area-weighted
Lagrangian (AL) and flow-weighted 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 flow-weighted 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 flow-weighted Lagrangian (AF). .. 95

4.9 Assumed time-dependent 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 time-dependent 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 velocity-log 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 first-order 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 = x/Vh (1.1)

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 first-order approximation" (CFOA) of the travel time. Widespread adoption of the CFOA through out the stochastic-Lagrangian 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 over-predicts 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 semi-analytical 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 space-stationary velocities along equal flow-weight streamlines and time-stationary velocities along equal area-weight 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 stochastic-Lagrangian 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 stochastic-Lagrangian 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 first-order 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
LAGRANGIAN-EULERIAN FLOW FIELD RELATIONSHIP

Develop an intuitive judgment and understanding for everything.

Miyamoto Musashi1

2.1 Introduction

Consider the steady, well-developed flow of water in an aquifer meeting Bear's [1972] definition of a porous medium. In most natural settings, gentle gradients induce a flow well-quantified 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 two-dimensional 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 [L3T-1]). 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 spatially-averaged 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 trajectory-based, equivalent.

2.2 Flow Field Partitioning Strategies A streamline is a hydrodynamic entity that is at all points tangent to the velocity vector. For two-dimensional 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 two-dimensional flow. For a lucid discussion of streamlines, stream functions, and stream surfaces in two and three-dimensional 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 well-developed 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 "micro-Darcy" 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 sub-domains. 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 sub-domains. 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 sub-domains of equal volume and gives an equal "volume weight" to each sub-domain. The division also creates two sub-IPs of equal area and gives an equal "area weight" to each sub-domain. For a perfectly homogeneous porous medium, the groundwater velocity vectors will be parallel to this line, and the discharges through each sub-domain 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) equal-area streamline, C) equal-flow streamline, D) equal-flow straight line.

the volume of water contained in each sub-domain 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 sub-domains are unequal in this case. Again, the magnitude of the areas through which water enters the sub-domains 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 sub-domain 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 sub-domains 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 finite-difference grid may be thought of as a set of points at equal spacing in space along a collection of area-weighted 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 flow-weight or an equal injection area-weight.

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 area-weighted Eulerian-trajectory oriented field. A Lagrangian field refers to a flow-weighted Lagrangian-trajectory 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 two-dimensional 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

mobile water content 0 = 1/5

From these criteria and Darcy's law, the area-averaged velocity is U = KeJ/O = 1/10. Clearly, the parameter values are artificial and selected for convenience, but values for Ke in md-1 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 commonly-used 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 statistically-independent 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

= exp[nplnk + n22nk/2] (2.6) where n is the order of the moment.

In the absence of a perfectly described physical aquifer, the method of turningbands generated equally-likely realizations of a hypothetical conductivity field with

16

prescribed ensemble statistics [Tompson et al., 1989]. Turning-bands was capable of producing two and three-dimensional 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 two-dimensional 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 turning-bands method compared to the standard normal distribution. Target distribution parameters are 1lnk = 0 and r2nk = 1.

turning-bands 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
- - --t-a
0
0 1 2 3
r/AInk

Figure 2.5: Estimated log conductivity longitudinal covariance estimated from a realization generated by the turning-bands 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 near-field and far-field 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 two-dimensional, 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 two-dimensional 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 finite-element 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 finite-element solution to the flow field was specified, following the results and arguments of Mos6 et al. [1994]. The mixed finite-element 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 "particle-tracking" 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 equally-spaced 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 near-ergodic 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 single-realization, the following method generates a large number of statistically independent realizations

20

generate a large, two-dimensional 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 equally-likely as any other. The trajectories associated with these particles carry an equal area weight, analogous to equal-area 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 area-weighted 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 equally-likely statistically-independent 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 flow-weighted 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 round-off 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 two-pass 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 area-weighted trajectory and FW denotes a flow-weighted trajectory ET denotes an Eulerian trajectory and LT denotes a Lagrangian. A particular weighttrajectory combination is concatenated the with a /. Thus, an area-weighted 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 area-weighted Eulerian trajectory (Figures 2.2 and 2.3 A) AW/LT area-weighted Lagrangian trajectory (Figures 2.2 and 2.3 B) FW/LT flow-weighted Lagrangian trajectory (Figures 2.2 and 2.3 C) FW/ET flow-weighted 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 area-weighted 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 turning-bands 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 zero-slope line through the data using a generalized least-squares 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~ ~ --- f--x
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

0-0

-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
C-4 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.75-0.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 area-weighted 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 --x--x-- ----- -- --- --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 area-weighted Eulerian trajectories in space and observed area-weighted 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 "large-displacement" 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 time-averaged timedependent velocity and the arithmetic averaging implicit with a space-averaged 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 equal-time spacing (or equally-likely observations in time) more "measurements" will come from low velocities, since the particle tracing the trajectory spends more time going slow. For an equal-distance 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 position-dependent 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 position-dependent 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 near-field" and transitions to an asymptotic value "in the far-field." 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 zero-slope 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 far-field 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 low-velocity 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 equal-flow sense. Thus, equal-area weighting preferentially weights low velocity zones. That there are implications for sampling schemes in heterogeneous media, especially in context of multi-level 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 two-dimensional 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 o-nk = 1 set of realizations. Covariances presented as dimensionless lag correlation functions.

that the analytical expression over-estimates 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 zero-lag covariance predicted by placing the input simulation parameters into Equation 2.13. The AW/ET velocity variance appears well-predicted 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 zero-slope 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 equal-area 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 streamline-based 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 time-dependent 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 2--69 -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 time-dependent 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 closed-form predictor behaves similarly to its spatial counterpart: early-time correlation is under-predicted and long-time correlation is over-predicted (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 rank-ordered 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 two-dimensional 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 well-described 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.

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 one-dimensional 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 zero-mean perturbation, e.g., In k = l.ink + f. Let the head gradient likewise be resolved into a mean and zero-mean 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 second-order, 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 second-order, 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 cross-covariance. 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]

= -2n k/d (2.19) where d is the dimension of the domain. This negative correlation between head gradient and log conductivity has a profound effect on the flow. The effective conductivity may be defined as that value when multiplied by the spatial mean gradient results in the observed specific discharge (flow rate divided by discharge area). The a priori effective conductivity is the effective conductivity estimated from Equations 2.18 and 2.19. With d = 2, we have for two-dimensional conductivity fields, the a priori

42

effective conductivity is simply Ke = Kg. Notice that the effective conductivity can be significantly less than the arithmetic mean conductivity Ka = Kgexp[aonk/2].

The definition of the effective conductivity is quite general, and a general effective conductivity expression is denoted Ke = Kgc 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, closed-form expressions are sought for the Eulerian, or spatial, mean velocity U, the Eulerian velocity covariance C,a[r], and its point distribution, given that the underlying conductivity field is large, stationary, lognormally-distributed, exponentially correlated, and characterized by parameters lnk, a12nk, and A2, 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

U cJ= (2.20)
0

where c = 1 for a two-dimensional isotropic aquifer. The c is retained to emphasize that this expression is not the "consistent first-order approximation" of the velocity that would result from dropping the a2nk and < fj > terms in Equation (2.18). For the two-dimensional isotropic aquifer, which is the model case, these two terms happen to be equal in magnitude and opposite in sign, and the consistent first-order

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 ank = 1, their magnitude is one-half 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, first-order expansion of Equation (2.16). Equation (2.17) would contain no f2 term, and Equation (2.18) would contain no aUnk 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 second-order terms and higher, the mean-removed expression, or "velocity perturbation," is S= J (f + j) (2.21)
0

44

If these processes are stationary stochastic processes, they can be represented represented by a Fourier-Stieltjes integral (see e.g., Bakr et al. [1978]) u'=J exp[iM x]dZ,[M] (2.22)

where i = v I, M is the wave number vector, and dZ,, is the complex random amplitude of the velocity perturbation. The Fourier-Stieltjes representation of Equation (2.21) is

KJ
dZ,, = (dZf + dZj) (2.23) Compare this expression to that derived by Shapiro and Cvetkovic [1988] (their Equation (A4))

0
dZ/= K (dZf + dZ) (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/O)2 on the one hand and (O/(KJ))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 first-order series approximation of the inverse velocity is linear in the velocity itself (e.g., 1/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 two-dimensional field is given by (Shapiro and Cvetkovic [1988])

C,[s] = 3ank(KJ/O)2 (exp[-s](s2 + 3 + 3s-4) + s-2/2 3s4) (2.25)

45

where s = x/Alnk. The corresponding expression for a three-dimensional field is

C,[s] = 80ank(KgJ/o)2 (exp[-s](s-2 + 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

approximation full expression -estimated covariance +
0.8
C:
4
0.6

"t 0.4

0.2 ',

0
0 2 4 6 8 10 r/A lnk

Figure 2.17: Comparison of velocity covariance expression Equation (2.25) to the approximation Equation (2.13) for aln k = 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 long-range persistence (see Figure 2.17). This behavior is examined with an expression adapted

46

from Dagan [1984]

Rua[r] = Onk (RIn k k [rdr d[rr2 (2.27) where R1nk is the log conductivity correlation function, RlnkO is the log conductivity and head cross-correlation function, and 'y, is the head variogram. From this expression, glean that the velocity covariance is a process comprised of three dominant subprocesses, namely Rlnk, RInk, and yo. At small displacements, all three of these subcomponents exhibit rapid trend towards zero (see Figure 2.18). The RInk, and 7y

log conductivity
C
O
.4j

0

-0.5
0 2 4 6 8 10 r/Aln k

Figure 2.18: Components of the velocity covariance function.

expressions exhibit effects that persist well beyond that of the log conductivity field. The long-range correlation effects of the head field in the mean direction of flow are cannot be precisely predicted with the simple exponential function Equation (2.13). In order to recover the proper correlation length, the approximation must systematically over-predict the velocity correlation at small separations. The "cost" of the

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 Pink, a2nk, and in k and saturated with mobile water at a constant and homogeneous content 0 is lognormal with approximate parameters

MIn u = In K9Jc (2.28)
0 V1 + b2n k

and

2nu = In[1 + boa k] (2.29) (see Appendix A for details of development). A fundamental assumption in addition to the lognormality hypothesis is that the "real space" moments of u are Pu = U and Ua U2ba,2k. Equation (2.13) is adopted as an approximation of the velocity covarince function.

In lieu of a rigorous proof, the results of a numerical experiment in which statistically independent velocities are drawn uniformly in space from simulated flow fields are presented (see Figure 2.16). The normalization of the velocities in Figure 2.16 than that of previous figures in that the log velocities are normalized by the sample mean and standard deviations to emphasize the lognormality of the velocity distribution. Normalization by the moments predicted by Equations (2.28) and (2.29) would have resulted in a lognormal distribution function that deviates from the standard

48

normal. This indicates that the moments I. = U and a2 = U2bou k 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 equal-area weighted partitioning scheme, and results in N sub-IPs of equal area, but dissimilar discharge. The discharge through sub-IP i is di = OAaui (2.30)

where Aa = A/N. The relative flow-weight of this streamtube i is given by wfi = ui/U (2.31) An equal-flow partitioning scheme results in in N sub-IPs of equal discharge, but of dissimilar area. The discharge through sub-IP i is d = OAaivi
(2.32)
=Q/N

Solving for the area yields

d
Aai = (2.33) Ovi

The relative area-weight of this streamtube i is wai = U/vi (2.34)

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 sub-IPs that water enters the streamtubes and a description of the sub-IP properties is essential to streamtube, and thus flow field, descriptions. If from a collection of N sub-IPs the random selection, or drawing, of one sub-IP 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 sub-IPs are sampled. The average sub-IP discharge is the same for the equal-area and equal-flow partitioning schemes, and equal to d = OAU/N. However, the discharge variance is zero for equal-flow and nonzero for equal-area, 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 equal-flow partitioning schemes. The expected value of v must equal the expected value flow-weighted u. Thus, = V

U
=< -U>
1
U (+ u) (2.35) = U(1 + ba2n k)
K, Jc 2
S Jc (1 + bl Ik)
0

50

The variance of v is given by var[v] =< v2 > < v >2. The expected value of v2 must equal the expected value of flow-weighted u2. That is =< UU2
U
1 9
= 1 exp[3Pini + -1,2nu
U 2 (2.36)
1 exp[3(pinu + c2nu/2) + 312nu] = U2(1 + bak)

Thus, the variance of v is given by var[v] = V2b2, k (2.37)

It is worthy to note that these results do not follow from an assumption that v = u2/U. That is, the velocity along the Lagrangian trajectory is not simply the flow-weighted velocity along an Eulerian trajectory.

Hypothesis: The point velocity distribution along a flow-weighted Lagrangian trajectory in the model stationary Eulerian velocity field is lognormal with approximate parameters

Pin v = In [KJ 1 + b2 nk (2.38) and

O nv= n[1 + b Ink] (2.39)

Moreover, the process is assumed second-order stationary in displacements along the mean direction of flow, that is, at control planes, and is characterized by the covariance function

CV,[r] = V2ba k exp[-br/Alk] (2.40)

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 p, = V = U(1 + bulnk) and a = V2bcr' k. 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 area-weighted Eulerian-trajectories and flowweighted Lagrangian-trajectories are second-order 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 area-weighted Lagrangian and the flow-weighted Eulerian streamtubes.

2.5.4 Control Plane Oriented Velocity Moments

Cvetkovic et al. [1996] postulated a semi-analytical estimator for a nonstationary mean velocity along area-weighted Lagrangian-trajectory observed at control planes. We would like to work directly from the assumed stationary flow-weighted Lagrangian-trajectory velocity field. The general expression for the expected value of the area-weighted Lagrangian-trajectory velocity is

U
'a[x] =< v--v[x]> (2.41) v [0]

Resolve the stationary mean velocities into the mean of the flow-weighted Lagrangiantrajectory field and zero-mean perturbation.

0a[x] =< V(1 + v'[zx) > (2.42) V(l + v'[01

52

Expand the denominator in a series to second order va[X] =
2 Cv[r]
= U(1 + o vf ) (2.43) V2 V2
= U(1 + ba2nk(1 exp[-br/Alnkl)) In a similar fashion, an expression is derived for the nonstationary mean velocity along a flow-weighted Eulerian-trajectory. 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[0]
zil[z] =< Uu[z] >

1< U(1 + u'[O])U(1 + u'[x]) > U C[r] (2.44) = U(1 + )
U2
= U(1 + bulnk exp[-br/Al k])

The nonstationary velocity variances along the respective cross-weighted trajectories are derived in a similar fashion. The general expression for the area-weighted Lagrangian-trajectory velocity variance is

U
var[va,[x]] =< (v[x] 1 a[x])2 > (2.45) v[0]

This equation may be re-written

var[va[x]] =< U V[x]2> U[X]2 (2.46) v[O]

Resolve the first term of the preceding equation into perturbation expressions U V2(1 + v'[x])2
< v[x]2 >= U < V(+[ (2.47) v [0] V(1 + V[0])

53

Expanding the denominator in a series, taking the expected value, and neglecting third-order and higher moments yields

U
< v[]2>= UV(1 + 2bo nk(1 + exp] bx/Alnk])) (2.48) In order to improve the quality of the estimate, the right-hand side (RHS) of the preceding equation is required to match the known extreme values of the left-hand (LHS) side. The LHS at X = 0 is

U
< vov[0]2 > = U v [0] (2.49) = UV

The RHS at X = 0 is
RHS[0] = UV(1 + 2bo,2 k(1 exp[-b0/Al k]))
(2.50)
= UV

The statistical independence of the initial velocity and the velocities as x goes to infinity are exploited to calculate the LHS as x -4+ c LHS[oo] =< > v [0]
= o2 V2 (2.51) = UV(1 + bunk)2 The RHS as x -+ oo is
RHS[oo] = UV(1 + 2ba nk(1 exp[-boo/Alnk]))
(2.52)
= UV(1 + 2ban k) The RHS is deficient a term (ba nk)2. Rewrite the RHS thus

RHS = UV(1 + (2balnk + (blk)2)(1 exp[-br/Al}k])) (2.53)

54

Substituting Equations (2.53) and (2.43) into Equation (2.45) yields

var[va[x]] = UV(1 + (2ba2nk + (bu2nk)2)(1 exp[-br/Alnk]))
(2.54)
-(U(1 + bonk(1 exp[-br/Alk ])))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

var[va[x]] = V2bnk U2bnk (2b n k + (bank)2) exp[--br/Alnlk] (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

55

are not a panacea, however, and the approximations have only been tested for the log conductivity variances tested in this work.

The nonstationary flow-weighted Eulerian-trajectory velocity variance may be derived in a similar fashion. These results are presented directly. The complete expression is

var [uf [x]] = U2 (1 + b n + (2b2nk +3 (b2n2 (bbo ( nk)) (1 exp [-br/Alk))

(U (1 + boInk exp [-br/Alnk]))2
(2.56)

and an approximation of the above obtained by neglecting a few terms is

var [u1 [x]] = U2bo2nk (1 + (2b2nk + (bunk)2) exp [-br/Alnk]) (2.57)

2.5.5 Control Time Oriented Velocity Moments

The displacement-invariant 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 area-weighted Lagrangian-trajectory 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 area-weighted Lagrangian-trajectory point velocity statistics are equal to the flow-weighted Lagrangian-trajectory point statistics. At the injection plane, the area-weighted Lagrangian-trajectory point velocity statistics are identically equal to the area-weighted Eulerian-trajectory 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 Lagrangian-time dependent velocity covariance function.

56

The area-weighted velocity at the IP is U by definition. Since the area-weighted Lagrangian-trajectory 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

T [t]dt= X d)-1 (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

Vh = exp[plnv alnv/2]

= exp[,ilnv] .exp[ainv]1/2
1 (2.59)
= U 1 + ba 1
I1 + ba.2n k
=U

Similarly, the harmonic mean of u[x] is given by

U
Uh = + 2 (2.60)
1 + b ln k

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 time-dependent velocities are identical to the initial space-dependent 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 area-weighted Lagrangian trajectory is C,,[t = 0] U2ba2nk. As mentioned, at large times, the statistical properties of the Lagrangiantrajectories converge. Making a change of variables x = Vht = Ut in Equation (2.40) and substituting the known Ca[t = 0] for the value of Crf[r = 0] in yields an

57

expression for the stationary area-weighted Lagrangian-trajectory velocity covariance

Cva[t] = U2bUank exp[-bUt/Alnk] (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 Lagrangian-trajectory velocity, and that it is quantitatively equal to the area-average velocity of the flow field. This is an important point.

The mean velocity and velocity variance along the flow-weighted 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 area-weighted Lagrangiantrajectories that exhibits a stationary velocity. The mean velocity is v[0]
0pr[t]=< U v[t]> (2.62) Expand the velocities into a mean and a zero-mean perturbation. The mean here is the time average mean, or the harmonic mean of the space average.

Vft [t] = U < U (1 + v' [0]) U (1 + v' [t]) >

= U t1 + (2.63) = U (1 + bInk exp [-bUt/Aln k])

This expression requires no series expansion. For the control plane oriented velocity statistics, the mapping was from flow-weighted to area-weighted trajectories. Equation 2.63 is identical in form to Equation 2.44 with the substitution r = Ut. Recognize

58

the quantitative similarity and write velocity variance equations directly The complete expression is

var [v []]= U2 (1 + bank + (2bonk + 3 (blk)2 + (ba k)) 3 (1 exp [-bUt/Alnk]))

(U (1 + bo, k exp [-br/Alnk]))2
(2.64)

and the approximation is

var [vf [t]] U2bIk (i 1+ (2blnk + (bok )2) exp [-bUt/Alnk]) (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 area-weighted velocities start at U and that the flow-weighted velocities start at V and both progress to U/(1 + ban k). Recognizing the similarity shared by the different mean time estimators, write approximate expressions directly

at t] 1 + b (1 + bor k exp [-bUt/Alnk]) (2.66) and

at [t]= 1 +blk (1 + (2boInk + b 1nk) (exp [-bUt/Alnk])) (2.67) 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

at = U,2bank
U= 2 (2.68)
1 + ban k k

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

59

area-weighted Eulerian-trajectory time-dependent velocity variance approximation is

var[ua[t]]= ( U l )b bnk (1 + (2bo~k + (bo2nk)) exp [-bUt/AXIk]) (2.69) The flux-weighted Eulerian-trajectory time-dependent velocity variance approximation is

var[uf[t]] =

1 Ub2 k) b,2nk (1 + (3bu,2k + 3 (b2nk 2 + (bu2 k )3 exp [-bUt/Alnk]) (2.70)

2.5.6 Log Conductivity Moments

The first two moments of the log conductivity along those trajectories for which the moments are stationary at control planes located along the mean direction of flow or at control times are sought. Working from a stationary Eulerian reference field, the mean and variance along the AW/ET are llnk and a2,k, respectively. Infer that the FW/LT is stationary as well, analogous to the velocity field results. The flow-weighted mean log conductivity perturbation is given by u[0]f Kexp[f]J(1 + j)O
<- > = < f>
U OKJ (2.71) =<(1 + f + f2/2 +...)(1 + j)f> /c Expanding, retaining second-order terms, and taking the expected value yields

< If>= nk(1 1/d)/c (2.72) which, for two-dimensional fields, is lnk/2. Similarly, the variance is

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 trajectory-based cross-covariance 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 = K9J/1 = 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 b = 3/8 and a2nk = {1/2, 1} depending upon the set of simulations. The deviation of the mean of area-weighted Eulerian-trajectory data from 1 are an indication of the quality of the effective conductivity expression.

The area-weighted Eulerian-trajectory mean velocity appears stationary, in accordance with established theory, and its value appears to be well predicted by the effective conductivity relationship. The flux-weighted Lagrangian-trajectory 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 area-weighted

61

Eulerian-trajectory velocity variance estimated from the input parameters to the simulations. That is, by

U2b,2nk = (KgJ/0)2ba0nk (2.74)

As with the mean velocity estimators, the variance estimators employ only the simulation input parameters.

Once again, the area-weighted Eulerian-trajectory variance exhibited its presumed behavior, and supports the well established fact that the Eulerian velocity field is fairly well understood in space. The hypothesized flow-weighted Lagrangiantrajectory velocity variance is somewhat higher that the data for both an k. 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 time-dependent 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 area-weighted Lagrangian-trajectory 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 over-predictive bias.

The estimated harmonic mean of Eulerian trajectories appear to have a stronger under predictive bias. The transition of the flow-weighted Lagrangian-trajectory 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

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 time-dependent velocity variances for the four different collections of streamtubes. The variances are normalized by the area-weighted Lagrangian-trajectory velocity variance estimated from the input parameters to the simulations. The area-weighted Lagrangian-velocity appears to be stationary in the variance. The data seem to lie slightly above the estimated value.

The over-prediction of the initial flux-weighted velocity variance has far less impact upon the time-dependent velocity variances. This combined with what appears to be a good estimate of the long-time Eulerian-trajectory 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 area-weighted Lagrangian-trajectory 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 space-dependent and time-dependent means although there seemed to be a slight over-predictive bias for the flow-weighted trajectory (see Figures 2.6 and2.8). The AW/ET log conductivity variance was "well predicted" in that the turning-bands 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

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 zero-slope line regressed through data. The decimal representation of the fractional value is included for convenience.
a2nk = 1/2
parameter predicted observed AW/ETal, k[z] 1/2 0.497 FW/LTU k[x] 7/16=0.4375 0.484 AW/LTank[t] 1/2 0.498 UInk = 1
parameter predicted observed AW/ET Ink[x] 1 1.002 FW/LTn k[x] 3/4 0.944 AW/LTuk [t] 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

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 pump-and-treat or as extravagant as an surfactant-enhanced 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.

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 well-known results are anticipated, in addition to some that are new. However, the approach is novel in that the work is based upon trajectory-based 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 area-weighted Eulerian-trajectories, 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 co across the entire face of the injection plane (IP) of the model aquifer (see Figure 2.1). Let the mass occupy

65

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 area-weight 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 Mr = coOhLAx (3.1)

The mass weight for any one streamtube i is m 1 (3.2) M N

This value is identically equal to the area-weight 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 mi ai
M L
U (3.3)

vi

67

3.2.2 Injection in Flux

Consider now maintaining the IP at co for some brief interval At. A mass

Mf = coOhLUAt (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

L
mi = coOh-viAt (3.5)
N

enters equal-area streamtube i. This corresponds to a mass weight of mi ui
(3.6)
M U

A mass

mi = co0haiviAt
= co QAt (3.7)
N

enters equal-flow 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) 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

68

is

t[ = (3.9) oiX vi[]

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 area-weighted and flowweighted trajectories is sought. Previous results indicate that the velocity for a collection of area-weighted streamtubes is stationary in time, characterized by mean U and approximate covariance U2ba 2nk exp[-bUt/An k]. The expected value of Equation (3.8) is

=< v[t]dt>

=Uf (3.10) = Ut

Although the velocity of the flow-weighted streamtube collection is not stationary in time, the results of the previous chapter are employed to map the properties of the area-weighted collection to the flow-weighted collection.

=< U fv[t]dt> (3.11) Recognize that this equation may be re-arranged to be the integral of Equation (2.63) in time. Substitute the RHS of Equation (2.63) into the preceding equation and

69

integrate to yield < xf[i] >= Ut + AlInkI2nk(1 exp[-bUt/Alnk]) (3.12) The displacement variance for the area-weighted streamtube collection is given by

var[xa[t]] =< Xz[t]2> 2

=< Ci[t'],t"]dt'dt"> -(Ut)2 = e t C rt at[t t]dt'dt" (U t)2
2 / j Cat [s]dt'ds (Ut)2 (3.13) = 2 (t s)CIat[s]ds (Ut)2 = 2AInk k2nk Ut + lk (exp[-bUt/Aink] 1)

3.3.2 Travel Times

Expressions for the travel time statistics are sought, in analogy to those sought for displacement. The reference collection of velocities from which these expressions will be derived is the stationary velocities of the equal flow streamtubes, or the equal-flow weight Lagrangian trajectories. The mean travel time for these streamtubes is given by

1 dx

1 dx (3.14) Vh

U
The mean travel time for the area-weighted streamtubes is / U dx
a v[] (3.15)
(x + Alnk12 nk(1 exp[-bx/Alnk]))
U

70

The arrival time variance for the flow-weighted collection of streamtubes is

var[tf [x]] =< X> (3.16) ] v [x'] o v [x"](3.16) 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
"] V2(1 + v'[x'])(1 + v'[x"11])

=< (1 v'[x'] + v' [']2 .. .)(1 '[x"] + v'[x"]2 >
V21 (3.17)
= V(1 + 2ao2 + Cvf[x', x"] + HOT)
1 ( Cvif [X', X"])
U2 V2

where HOT is the collection of higher order terms. Inserting this expression into Equation 3.16 yields an expression for the travel time variance for flux-weighted streamlines

var[tf[zx]] = + A (exp[-bx/A2k]- )) (318)

3.4 Solute Transport

Expressions for the statistics of the integrated trajectory properties of displacement and travel time have been derived in the absence of an explicit reference to solute transport to demonstrate that these properties exist simultaneously in the same flow field. However, they are intimately related to solute transport. Consider the transport of an infinitesimal particle whose movement is constrained to follow a streamline. Moreover, at any time or position on the streamline, the particle velocity is identically equal to that of the underlying velocity field. Its displacement from x = 0 in t = T is given by Equation (3.8) and the time required to travel from x = 0 to x = X is given by Equation (3.9). For the same flow field, two collections of streamtubes have been discussed, namely area-weighted and flow-weighted. Placing

71

one particle into each area-weighted 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 flow-weighted 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

72

by a higher velocity. The time-varying mean velocity expression for the flow-weighted Lagrangian-trajectories (Equation (2.63)) is, in fact, an expression for the mean velocity of a plume resulting from a uniform flux injection. The long-time asymptotic slope of the flux-weighted 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 flux-injected 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 area-weighted trajectory displacement variance over-predicts the observed displacement variances, and the long-time 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 Lagrangian-trajectory velocity field in time appears to be well predicted. Thus, the problem may lie with the rather poor description of the short-time 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 in-flux 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

73

10

x4

8 AW/LT + FW/LT x "" Eq. 3.12 ---------Eq 3.10 .................X .44
4 .4 .4

x
2

0 I I
0 2 4 6 8 10
tU/ An k
10 .4

8 AW/LT + FW/LT x Eq. 3.12 -.-------- .X +
E q 3.10 ...................

k .4

XX

( 4
.4
x"

0
0 2 4 6 8 10
tU/ An k

Figure 3.1: Comparison of mean displacements for the two injection modes or streamtube collections. The flow-weighted Lagrangian trajectory (FW/LT) statistics correspond to a flux injection, and the area-weighed Lagrangian trajectory (AW/LT) statistics correspond to a uniform resident injection. The top figure contains the rnk = 1/2 data and the bottom figure contains the cink = 1 data.

74

15

AW/LT + ,* FW/LT x Eq. 3.13 ---------0

0 2 4 6 8 10
tU/Aln k
5

AW/LT + +x FW/LT x L +xt Eq. 3.13 -- -- +x
+X
10 +Xx normalized by +

0 2 4 6 8 10 tU/An k

correspond to a flux injection, and the area-weighed Lagrangian trajectory (AW/LT)

2 2

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 Vh = 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 non-linear 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 mass-weighting of low-velocity 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 near-to-far 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] = AlnkUln2k(1 exp[-bx/Alnk])/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

76

40

30 AW/LT + X FW/LT x Eq. 3.15 ---------- +
Eq. 3.14

0

0 10 20 30 40 X/Alnk
40 .

30 AW/LT + FW/LT x x Eq. 3.15. Eq. 3.14 ...............

20

10

0 10 20 30 40 x/Alk

Figure 3.3: Comparison of mean arrival times for the two injection modes or streamtube collections. The flow-weighted Lagrangian trajectory (FW/LT) statistics correspond to a flux injection, and the area-weighed Lagrangian trajectory (AW/LT) statistics correspond to a uniform resident injection. The top figure contains the aUink -- 1/2 data and the bottom figure contains the 02k =1 data.
In k =I data

77

80
+

tx

60 AW/LT + .x FW/LT x x C14Eq. 3.18 -------- x

S40.*x

20 40

'+Xx.
.4X

S0 0x

+ +
x.
+ X + x..'

20 .*+
++

0 AW/L +

0 10 20 30 40

In k
4-.
+x
60 A"

0 10 20 30 40

78

0.4 1
:n.k = 1/2 -0.3

t 0.2

0.1

-------------------------------------------0 I
0 10 20 30 40
/Al nk

Figure 3.5: Travel time estimate error as function of displacement associated with assuming a uniform resident injection when injection is in flux for Cink = {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 so-called consistent first-order 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" first-order approximation be applied in future works, since this Lagrangian expression is inconsistently based upon the Eulerian mean velocity.

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 closed-forms.

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

80

into collections of noninteracting streamtubes. Dispersion at the streamtube scale is considered to be negligible in its effect on the aggregate-scale 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 (equal-area or equal-flow) and injection mode (uniform resident and in-flux). The reference collection is equal-flow 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 multi-dimensional flow field that is resolved into a collection of noninteracting streamtubes. Transport along an individual streamtube is characterized with a transfer function y[t; x] (see Jury and Roth [1990] and Cvetkovic et al. [1998]). For a conservative and nonreactive solute, the transfer function is simply

y[t; x] = 6 [t T[x]] (4.1) where 6 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 = T[X]. We conceive of this travel time as not only the time required to

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 x P[x]dx (4.2) o V [x]

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 p properties in general are broadly applicable to any distributed reaction parameter, perhaps a first-order 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, p 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 7 and the time sorbed M. The transfer function characterizing this transport and reaction process is (see Cvetkovic et al. [1998])

-y[t; x] = 6 [t 7[x] P[x]] (4.3)

Following the work of Cvetkovic et al. [1998], the adopted reaction parameter model is

P = Pg exp[w] exp[P In k] (4.4)

where Pg is the geometric mean of P, w is a zero-mean, exponentially-correlated normal random variable with variance of and correlation length A, = Alnk 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 A,, = Alnk was specified for convenience.

82

Two specific injection modes are considered: uniform injection in-flux and uniform resident injection. For the uniform injection in-flux, 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

mi[t; x] = mi[0; 0]y[t; x] (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 in-flux, the expected solute breakthrough curve is given by = M <[t; x] > (4.6) For uniform resident injection, the expected solute breakthrough curve is given by = MU (4.7) where U is the arithmetic mean, or spatially averaged, velocity and v[O] is the velocity at the injection point. The temporal moment of this breakthrough curve is sought. The n temporal moment is given by 8n
tn = lim(-1)n (4.8) s-+O aSn

where S is the Laplace transform of S and s the corresponding Laplace space variable. The Laplace transform of Equation (4.3) is given by '[s; x] = exp [-s (T[x] + p[x])] (4.9)

83

Two temporal moments of interest for the in-flux breakthrough curve are the mean tlf =<7> + (4.10) and the variance

t2f = T + 2 + 2a,, (4.11) For uniform resident injection, the moments are tlf = + (4.12) and the variance

t2f =< U(T +)2/v[0] > < U/v[0] >2 < UI/v[0] >2 (4.13) From these equations, it is seen that the mean and variance of the expected reactive travel time distribution is comprised of the first two moments of the nonreactive travel time and the reaction flow path and the cross moment of the travel time and reaction flow path.

4.4 Simulation of Reactive Solute Transport

The results of the simulations described in Chapter 2 are employed to simulate the transport of a reactive solute. As mentioned in that chapter, information was recorded along two particle trajectories: the first was a natural trajectory corresponding to a hydrodynamic streamline, and the second was a straight line parallel to the mean flow direction. Data were recorded at equal increments in time after injection, and at "control planes" perpendicular to the mean flow at equal displacements from the injection plane. Among these data were the velocity, the hydraulic conductivity, and an observation of a lognormal random field with the same statistical and correlation properties as the log hydraulic conductivity field.

84

4.4.1 Reaction Parameter

From these data, trajectory-based 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 = {-1, -1/2, 0,1/2, 1} in log conductivity fields of variances ank = {1/2, 1}. The extreme correlation parameters / = {-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, 0 = {-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, = (1k-,2)a12nk. 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, = 1. In the case of the two-dimensional isotropic log conductivity field, the effective conductivity is equal to the geometric mean. The expected value of the reaction value was fixed < P >= 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 equal-flux streamtubes exhibit an apparently stationary mean reaction parameter in displacement along the mean flow direction, whereas the equal-area streamtubes exhibit distinctly nonstationary behavior that is analogous to that of the log conductivity (see Figures 4.1 and

85

4.2). These observations, of course, are not surprising, since the reaction parameter is a function of the log conductivity. Positive correlation between the log conductivity and the reaction parameter (0 > 0) results in a reaction parameter trajectory mean that is larger than the uniform spatial mean. Negative correlation (3 < 0) results in a reaction parameter trajectory-mean that is smaller than the uniform spatial mean. No correlation (p = 0) results in equal trajectory-based and uniform spatial mean reaction parameter estimates. A decrease in the system variability (e.g., the variance of log conductivity ank) results in a smaller difference between the volume average and the trajectory averages.

4.4.2 Reaction Flow Path

A simple trapezoidal rule was used to calculate the reaction flow path p along a solute particle trajectory. For a collection of n trajectories, there were n p "observations" at each control plane. For an intuitive grasp of what the statistics of the reaction path "means", recall that for linear equilibrium sorption, the reaction path corresponds to the time that the solute is sorbed. The sum of the the nonreactive travel time and the sorbed time equals the total reactive travel time of the solute from the point of injection to the point of observation. See Cvetkovic et al. [1998] for further application of the statistics of the reaction flow path.

The mean and the variance of the reaction flow path are presented for Lagrangian trajectories that carry either an equal flow weight or an equal area weight at a reference plane. Recall that the equal flow collection can be thought to correspond to a uniform injection in-flux and that the equal area collection can be thought to correspond to a uniform injection into the resident fluid. The coupled effect of injection mode and reaction parameter correlation to the log conductivity are illustrated by examining reactive flow fields with five different correlation characteristics. The first case is that in which the reaction parameter is completely determined by the

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

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 fill-in-the-blank 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 LAGRANGIAN-EULERIAN 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 Path-Travel 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 area-weighted Eulerian trajectories in space and observed area-weighted 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 zero-slope 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) equal-area streamline, C) equal-flow streamline, D) equal-flow 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 turning-bands 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 turning-bands 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 position-dependent 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 position-dependent 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 time-dependent 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 time-dependent 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 two-dimensional 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 flow-weighted 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 flow-weighted 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 area-weighted 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 area-weighted 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 time-dependent 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 time-dependent 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 first-order 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 first-order approximation" (CFOA) of the travel time. Widespread adoption of the CFOA through out the stochastic-Lagrangian 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 over-predicts 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. = x/V h (1.1) = x/U (1.2)

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 semi-analytical 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 space-stationary velocities along equal flow-weight streamlines and time-stationary velocities along equal area-weight 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 stochastic-Lagrangian 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 stochastic-Lagrangian 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 first-order 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 LAGRANGIAN-EULERIAN FLOW FIELD RELATIONSHIP Develop an intuitive judgment and understanding for everything. Miyamoto Musashi 1 2.1 Introduction Consider the steady, well-developed flow of water in an aquifer meeting Bear's [1972] definition of a porous medium. In most natural settings, gentle gradients induce a flow well-quantified 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 two-dimensional 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 spatially-averaged 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 trajectory-based, equivalent. 2.2 Flow Field Partitioning Strategies A streamline is a hydrodynamic entity that is at all points tangent to the velocity vector. For two-dimensional 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 two-dimensional flow. For a lucid discussion of streamlines, stream functions, and stream surfaces in two and three-dimensional 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 well-developed 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 "micro-Darcy" 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 sub-domains. 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 sub-domains. 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 sub-domains of equal volume and gives an equal "volume weight" to each sub-domain. The division also creates two sub-IPs of equal area and gives an equal "area weight" to each sub-domain. For a perfectly homogeneous porous medium, the groundwater velocity vectors will be parallel to this line, and the discharges through each sub-domain 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) equal-area streamline, C) equal-flow streamline, D) equal-flow straight line. the volume of water contained in each sub-domain 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 sub-domains are unequal in this case. Again, the magnitude of the areas through which water enters the sub-domains 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 sub-domain 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 sub-domains 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 finite-difference grid may be thought of as a set of points at equal spacing in space along a collection of area-weighted 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 flow-weight or an equal injection area-weight. 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 area-weighted Eulerian-trajectory oriented field. A Lagrangian field refers to a flowweighted Lagrangian-trajectory 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 two-dimensional 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 area-averaged 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 commonly-used 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 statistically-independent 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 equally-likely realizations of a hypothetical conductivity field with C\ nk [r] = af^expf-r/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]. Turning-bands was capable of producing two and three-dimensional 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 two-dimensional 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 -3-2-10 1 2 3 log conductivity Figure 2.4: Log conductivity distribution from a 41400 node realization generated by the turning-bands method compared to the standard normal distribution. Target distribution parameters are fj,\ n k = 0 and af nk = 1. turning-bands 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 turning-bands 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 near-field 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 two-dimensional, 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 two-dimensional 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 finite-element 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 finite-element solution to the flow field was specified, following the results and arguments of Mose et al. [1994]. The mixed finite-element 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 "particle-tracking" 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 equally-spaced 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 near-ergodic 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 single-realization, the following method generates a large number of statistically independent realizations

PAGE 35

20  generate a large, two-dimensional 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 equally-likely as any other. The trajectories associated with these particles carry an equal area weight, analogous to equal-area 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 area-weighted 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 equally-likely statistically-independent 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 flow-weighted 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 round-off 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 two-pass 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 area-weighted 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 flow-weighted 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 turning-bands 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 zero-slope line through the data using a generalized least-squares 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. "fl-|J.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 area-weighted 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 area-weighted 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 "large-displacement" 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 time-averaged timedependent velocity and the arithmetic averaging implicit with a space-averaged 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 equal-time spacing (or equally-likely observations in time) more "measurements" will come from low velocities, since the particle tracing the trajectory spends more time going slow. For an equal-distance 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 position-dependent 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 position-dependent 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 near-field" and transitions to an asymptotic value "in the far-field." 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 zero-slope 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 far-field 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 low-velocity 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 equal-flow sense. Thus, equal-area weighting preferentially weights low velocity zones. That there are implications for sampling schemes in heterogeneous media, especially in context of multi-level 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 two-dimensional 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 over-estimates 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 zero-lag covariance predicted by placing the input simulation parameters into Equation 2.13. The AW/ET velocity variance appears well-predicted 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 zero-slope 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 equal-area 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 streamline-based 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 time-dependent 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 > 1-2 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 time-dependent 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 closed-form predictor behaves similarly to its spatial counterpart: early-time correlation is under-predicted and long-time correlation is over-predicted (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 rank-ordered and subject to the transformation In v  In v \nv' J \n v (2.14)

PAGE 54

39 -3-2-10123 log velocity Figure 2.16: Normalized log initial velocities for two-dimensional 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 well-described 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

PAGE 56

41 may be resolved into a mean, and a zero-mean perturbation, e.g., In A; = fj,\ n k + / Let the head gradient likewise be resolved into a mean and zero-mean 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 second-order, 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 second-order, 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 cross-covariance. 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] =-olJd (2.19) where d is the dimension of the domain. This negative correlation between head gradient and log conductivity has a profound effect on the flow. The effective conductivity may be defined as that value when multiplied by the spatial mean gradient results in the observed specific discharge (flow rate divided by discharge area). The a priori effective conductivity is the effective conductivity estimated from Equations 2.18 and 2.19. With d = 2, we have for two-dimensional conductivity fields, the a priori

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, closed-form 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, lognormally-distributed, 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 two-dimensional isotropic aquifer. The c is retained to emphasize that this expression is not the "consistent first-order approximation" of the velocity that would result from dropping the a? nk and < fj > terms in Equation (2.18). For the two-dimensional isotropic aquifer, which is the model case, these two terms happen to be equal in magnitude and opposite in sign, and the consistent first-order

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 one-half 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, first-order 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 second-order terms and higher, the mean-removed 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 Fourier-Stieltjes integral (see e.g., Bakr et al. [1978]) poo u'= exp[iM-x]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 Fourier-Stieltjes 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 first-order 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 two-dimensional 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 three-dimensional 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 long-range 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] = ^lfc ( R\nk[r\ ^ J l^'J where R\ nk is the log conductivity correlation function, R\ n k is the log conductivity and head cross-correlation function, and 74, is the head variogram. From this expression, glean that the velocity covariance is a process comprised of three dominant subprocesses, namely R\ nk R\nk, and 7^. At small displacements, all three of these subcomponents exhibit rapid trend towards zero (see Figure 2.18). The R\ n k and 7^ o o -0.5 1 1 1 1 1 0 2 4 6 8 10 r/\\nk Figure 2.18: Components of the velocity covariance function. expressions exhibit effects that persist well beyond that of the log conductivity field. The long-range correlation effects of the head field in the mean direction of flow are cannot be precisely predicted with the simple exponential function Equation (2.13). In order to recover the proper correlation length, the approximation must systematically over-predict the velocity correlation at small separations. The "cost" of the

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 equal-area weighted partitioning scheme, and results in N sub-IPs of equal area, but dissimilar discharge. The discharge through sub-IP i is where Aa = A/N. The relative flowweight of this streamtube i is given by An equal-flow partitioning scheme results in in N sub-IPs of equal discharge, but of dissimilar area. The discharge through sub-IP 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 area-weight 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 sub-IPs that water enters the streamtubes and a description of the sub-IP properties is essential to streamtube, and thus flow field, descriptions. If from a collection of N sub-IPs the random selection, or drawing, of one sub-IP 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 sub-IPs are sampled. The average sub-IP discharge is the same for the equal-area and equal-flow partitioning schemes, and equal to d = 8AU/N. However, the discharge variance is zero for equal-flow and nonzero for equal-area, 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 equal-flow partitioning schemes. The expected value of v must equal the expected value flow-weighted u. Thus, = V u =< u> U = jjiA + l*u) (2.35) = U(l + bal k ) = ^(1 + W fc )

PAGE 65

50 The variance of v is given by var[u] =< v 2 > 2 The expected value of v 2 must equal the expected value of flowweighted u 2 That is 2 ^2 =< u l > = ^exp[3/ii nu + ^a 2 nu ] = exp[3(^, nu + a 2 nu /2) + Sa 2 n J (2.36) = U 2 (l + bcl k ) Thus, the variance of v is given by [v} = V 2 bo 2 nk (2.37) It is worthy to note that these results do not follow from an assumption that v  u 2 /U. That is, the velocity along the Lagrangian trajectory is not simply the flow-weighted velocity along an Eulerian trajectory. Hypothesis: The point velocity distribution along a flow-weighted Lagrangian trajectory in the model stationary Eulerian velocity field is lognormal with approximate parameters Mint; = ^ and K g Jc e (2.38) a 2 nv = \n[l + ba 2 nk } (2.39) Moreover, the process is assumed second-order stationary in displacements along the mean direction of flow, that is, at control planes, and is characterized by the covariance function C vf [r] = V 2 ba 2 nk exp[-6r/A lnjfc ] (2.40)

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 area-weighted Eulerian-trajectories and flowweighted Lagrangian-trajectories are second-order 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 flow-weighted Eulerian streamtubes. 2.5.4 Control Plane Oriented Velocity Moments Cvetkovic et al. [1996] postulated a semi-analytical estimator for a nonstationary mean velocity along area-weighted Lagrangian-trajectory observed at control planes. We would like to work directly from the assumed stationary flowweighted Lagrangian-trajectory velocity field. The general expression for the expected value of the areaweighted Lagrangian-trajectory velocity is Resolve the stationary mean velocities into the mean of the flow-weighted Lagrangiantrajectory field and zero-mean 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 + & ^) (2-43) = U(l + ba 2 nk (l-exp[-br/X lnk })) In a similar fashion, an expression is derived for the nonstationary mean velocity along a flow-weighted Eulerian-trajectory. 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 cross-weighted trajectories are derived in a similar fashion. The general expression for the areaweighted Lagrangian-trajectory velocity variance is var[u a [x]] =<-^-(v[x] v a [x]) 2 > (2.45) v[0\ This equation may be re-written 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 third-order 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 right-hand side (RHS) of the preceding equation is required to match the known extreme values of the left-hand (LHS) side. The LHS at X = 0 is <-T7^[0] 2 > = U v[0] 1 J 1 J (2.49) = UV The RHS at X = 0 is RHS[0] = UV{\ + 2bal k (l exp[-W/A lnfc ])) (2.50) = UV The statistical independence of the initial velocity and the velocities as x goes to infinity are exploited to calculate the LHS as x  > oo LHS[oo] =<-^j> = ol + V 2 (2-51) = UV(l + bo? nk ) 2 The RHS as x -> oo is RHS[oo] = UV(1 + 2baf nk (l exp[-6oo/A lnfc ])) = 1^(1 + 2^) The RHS is deficient a term {ba 2 nk ) 2 Rewrite the RHS thus (2.52) RHS = UV(1 + (2ba 2 nk + (ba 2 nk ) 2 )(l exp[-6r/A lnfc ])) (2.53)

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 (l-exp[-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 Eulerian-trajectory 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 displacement-invariant 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 area-weighted Lagrangian-trajectory 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 area-weighted Lagrangian-trajectory point velocity statistics are equal to the flow-weighted Lagrangian-trajectory point statistics. At the injection plane, the area-weighted Lagrangian-trajectory point velocity statistics are identically equal to the area-weighted Eulerian-trajectory 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 Lagrangian-time dependent velocity covariance function.

PAGE 71

56 The areaweighted velocity at the IP is U by definition. Since the areaweighted Lagrangian-trajectory 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 time-dependent velocities are identical to the initial space-dependent 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 area-weighted Lagrangian-trajectory 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 Lagrangian-trajectory velocity, and that it is quantitatively equal to the area-average velocity of the flow field. This is an important point. The mean velocity and velocity variance along the flow-weighted 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 area-weighted Lagrangiantrajectories that exhibits a stationary velocity. The mean velocity is Vft[t]=<^v[t}> (2.62) Expand the velocities into a mean and a zero-mean 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 area-weighted velocities start at U and that the flow-weighted 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*])) (2-67) 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 area-weighted Eulerian-trajectory time-dependent velocity variance approximation is varK[i]] = (y^-J Wnfc (1 + {^< k + {b/c Expanding, retaining second-order terms, and taking the expected value yields < f 1 />=4,(l-l/^)/c (2.72) which, for two-dimensional fields, is crf ak /2. Similarly, the variance is < ^o] /2> < !ffl />2= {
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 trajectory-based cross-covariance 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 Eulerian-trajectory data from 1 are an indication of the quality of the effective conductivity expression. The area-weighted Eulerian-trajectory mean velocity appears stationary, in accordance with established theory, and its value appears to be well predicted by the effective conductivity relationship. The flux-weighted Lagrangian-trajectory 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 area-weighted

PAGE 76

61 Eulerian-trajectory velocity variance estimated from the input parameters to the simulations. That is, by ^ = (V/)V hi (2-74) As with the mean velocity estimators, the variance estimators employ only the simulation input parameters. Once again, the areaweighted Eulerian-trajectory variance exhibited its presumed behavior, and supports the well established fact that the Eulerian velocity field is fairly well understood in space. The hypothesized flow-weighted 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 time-dependent 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 area-weighted Lagrangian-trajectory 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 over-predictive bias. The estimated harmonic mean of Eulerian trajectories appear to have a stronger under predictive bias. The transition of the flow-weighted Lagrangian-trajectory 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 time-dependent velocity variances for the four different collections of streamtubes. The variances are normalized by the area-weighted Lagrangian-trajectory 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 over-prediction of the initial fluxweighted velocity variance has far less impact upon the time-dependent velocity variances. This combined with what appears to be a good estimate of the long-time Eulerian-trajectory 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 area-weighted Lagrangian-trajectory 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 space-dependent and time-dependent means although there seemed to be a slight over-predictive bias for the flowweighted trajectory (see Figures 2.6 and2.8). The AW/ET log conductivity variance was "well predicted" in that the turning-bands 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 zero-slope 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 pump-and-treat or as extravagant as an surfactant-enhanced 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 well-known results are anticipated, in addition to some that are new. However, the approach is novel in that the work is based upon trajectory-based 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 area-weighted Eulerian-trajectories, 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 equal-area streamtube i. This corresponds to a mass weight of A mass rrii = coOhciiViAt (3.6) Q (3-7) = c 0 9h^At enters equal-flow 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 area-weighted and flowweighted trajectories is sought. Previous results indicate that the velocity for a collection of area-weighted 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 flow-weighted streamtube collection is not stationary in time, the results of the previous chapter are employed to map the properties of the area-weighted collection to the flow-weighted collection. Recognize that this equation may be re-arranged to be the integral of Equation (2.63) in time. Substitute the RHS of Equation (2.63) into the preceding equation and =< / v[t]dt> (3.10) = Ut (3.11)

PAGE 84

69 integrate to yield = Ut + X ]nk al k {l-exp[-bUt/X lnk }) (3.12) The displacement variance for the area-weighted streamtube collection is given by var[x a [t]] = 2 3.3.2 Travel Times Expressions for the travel time statistics are sought, in analogy to those sought for displacement. The reference collection of velocities from which these expressions will be derived is the stationary velocities of the equal flow streamtubes, or the equal-flow weight Lagrangian trajectories. The mean travel time for these streamtubes is given (3.13) 2 f \t s)C lat [s}ds (Ut) 2 by (3.14) x U The mean travel time for the area-weighted streamtubes is (3.15) = -(x + X lnk a 2 nk {l exp[-bx/X lnk }))

PAGE 85

70 The arrival time variance for the flow-weighted 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 > = u[x']u[a;"] V 2 {1 + v'[x']){l + v'[x"]) =< 4(1 v'[x'] + v'[xf -...)(!v'[x"\ + v'lx'f -...)> r 1 V2 (3.17) = ^(l + 2a 2 v + C vf [x',x"}+UOT) 1 (, C vf [x',x"Y U2 y L + y 2 where HOT is the collection of higher order terms. Inserting this expression into Equation 3.16 yields an expression for the travel time variance for fluxweighted streamlines var[t/N] = ~f (exp[-WA 2 nfc ] 1) j (3-18) 3.4 Solute Transport Expressions for the statistics of the integrated trajectory properties of displacement and travel time have been derived in the absence of an explicit reference to solute transport to demonstrate that these properties exist simultaneously in the same flow field. However, they are intimately related to solute transport. Consider the transport of an infinitesimal particle whose movement is constrained to follow a streamline. Moreover, at any time or position on the streamline, the particle velocity is identically equal to that of the underlying velocity field. Its displacement from x = 0 in t = T is given by Equation (3.8) and the time required to travel from x = 0 to x = X is given by Equation (3.9). For the same flow field, two collections of streamtubes have been discussed, namely area-weighted and flow-weighted. Placing

PAGE 86

71 one particle into each area-weighted 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 flow-weighted 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 Lagrangian-trajectories (Equation (2.63)) is, in fact, an expression for the mean velocity of a plume resulting from a uniform flux injection. The long-time asymptotic slope of the flux-weighted 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 flux-injected 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 area-weighted trajectory displacement variance over-predicts the observed displacement variances, and the long-time 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 Lagrangian-trajectory velocity field in time appears to be well predicted. Thus, the problem may lie with the rather poor description of the short-time 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 in-flux 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 flow-weighted 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 flow-weighted 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 non-linear 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 mass-weighting of low-velocity 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 near-to-far 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 area-weighed 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 flow-weighted Lagrangian trajectory (FW/LT) statistics correspond to a flux injection, and the area-weighed 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 so-called consistent first-order 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" first-order 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 closed-forms. 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 aggregate-scale 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 (equal-area or equal-flow) and injection mode (uniform resident and in-flux). The reference collection is equal-flow 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 multi-dimensional 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[t-T[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 first-order 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[t-r[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 zero-mean, exponentially-correlated 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 in-flux and uniform resident injection. For the uniform injection in-flux, 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[t-x] (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 in-flux, the expected solute breakthrough curve is given by = M <7M> ( 4 6 ) For uniform resident injection, the expected solute breakthrough curve is given by < S u [t] >= MU < 7[t; x]/v[0] > (4.7) where U is the arithmetic mean, or spatially averaged, velocity and v[0] is the velocity at the injection point. The temporal moment of this breakthrough curve is sought. The n temporal moment is given by t n = lim(-l) n |^<,SM> (4.8) s->0 QS n where S is the Laplace transform of S and s the corresponding Laplace space variable. The Laplace transform of Equation (4.3) is given by j[s; x] = exp [-s (r[x] + fi[x})} (4.9)

PAGE 98

83 Two temporal moments of interest for the in-flux breakthrough curve are the mean t lf = + (4.10) and the variance t 2f = a 2 T + al + 2a T ^ (4.11) For uniform resident injection, the moments are t lf =< Ut/v[0] > + < Un/v[0] > (4.12) and the variance t 2f =< U{t + n) 2 /v[0] > < Ut/v[Q) > 2 < Ufi/v[0] > 2 (4.13) From these equations, it is seen that the mean and variance of the expected reactive travel time distribution is comprised of the first two moments of the nonreactive travel time and the reaction flow path and the cross moment of the travel time and reaction flow path. 4.4 Simulation of Reactive Solute Transport The results of the simulations described in Chapter 2 are employed to simulate the transport of a reactive solute. As mentioned in that chapter, information was recorded along two particle trajectories: the first was a natural trajectory corresponding to a hydrodynamic streamline, and the second was a straight line parallel to the mean flow direction. Data were recorded at equal increments in time after injection, and at "control planes" perpendicular to the mean flow at equal displacements from the injection plane. Among these data were the velocity, the hydraulic conductivity, and an observation of a lognormal random field with the same statistical and correlation properties as the log hydraulic conductivity field.

PAGE 99

84 4.4.1 Reaction Parameter From these data, trajectory-based 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 two-dimensional 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 equal-flux streamtubes exhibit an apparently stationary mean reaction parameter in displacement along the mean flow direction, whereas the equal-area streamtubes exhibit distinctly nonstationary behavior that is analogous to that of the log conductivity (see Figures 4.1 and

PAGE 100

85 4.2). These observations, of course, are not surprising, since the reaction parameter is a function of the log conductivity. Positive correlation between the log conductivity and the reaction parameter (/3 > 0) results in a reaction parameter trajectory mean that is larger than the uniform spatial mean. Negative correlation (/3 < 0) results in a reaction parameter trajectory-mean that is smaller than the uniform spatial mean. No correlation (/? = 0) results in equal trajectory-based and uniform spatial mean reaction parameter estimates. A decrease in the system variability (e.g., the variance of log conductivity of nk ) results in a smaller difference between the volume average and the trajectory averages. 4.4.2 Reaction Flow Path A simple trapezoidal rule was used to calculate the reaction flow path fi along a solute particle trajectory. For a collection of n trajectories, there were n \i "observations" at each control plane. For an intuitive grasp of what the statistics of the reaction path "means" recall that for linear equilibrium sorption, the reaction path corresponds to the time that the solute is sorbed. The sum of the the nonreactive travel time and the sorbed time equals the total reactive travel time of the solute from the point of injection to the point of observation. See Cvetkovic et al. [1998] for further application of the statistics of the reaction flow path. The mean and the variance of the reaction flow path are presented for Lagrangian trajectories that carry either an equal flow weight or an equal area weight at a reference plane. Recall that the equal flow collection can be thought to correspond to a uniform injection in-flux and that the equal area collection can be thought to correspond to a uniform injection into the resident fluid. The coupled effect of injection mode and reaction parameter correlation to the log conductivity are illustrated by examining reactive flow fields with five different correlation characteristics. The first case is that in which the reaction parameter is completely determined by the

PAGE 101

86 1.8 1.6 1.4 1.2 08 0 6 0.4 1.8 1.6 1.4 X +FL +AL +AE -AE o -AL -FL o 10 20 30 40 1.2 ft. 08 06 04 x i S ** x i xx  1 1 iv^^^/^w* +FL +AL x +AE > -AE -AL -FL o 10 20 30 40 10 20 i/A]jt Figure 4.1: Mean reaction parameter values for different values of correlation parameter P and different trajectory collections for the af nk = 1 set of simulations. Top figure is /? = {-1,+1}. Middle figure is f3 = {-1/2, +1/2}. Bottom figure is /3  0 (no correlation). Trajectory collections are area-weighted Lagrangian (AL), flowweighted Lagrangian (AF) and areaweighted Eulerian (AE).

PAGE 102

87 18 1.6  14 12 tu *, **** *** 0.8 0.6 0.4 ( 1.8 16 1.4 1.2 +FL +AL +AE -AE o -AL -FL o 10 20 30 40 ** +FL +AL +AE -AE -AL  -FL Figure 4.2: Mean reaction parameter values for different values of correlation parameter (3 and different trajectory collections for the of nfc = 1/2 set of simulations. Top figure is P = {-1,+1}. Middle figure is 0 = {-1/2, +1/2}. Bottom figure is P  0 (no correlation). Trajectory collections are area-weighted Lagrangian (AL), flowweighted Lagrangian (AF) and areaweighted Eulerian (AE).

PAGE 103

88 conductivity field (ft = {-1,1}). Considered are the two sub-cases, namely that of "perfect" negative correlation (/5 = -1) and "perfect" positive correlation (/5  -1). For /3 =  1, the equal area collection exhibits nonlinear behavior analogous to that of travel time (see Figure 4.3 and 4.4). However, this nonlinearity is enhanced by the correlation of high reaction parameter values in the low conductivity areas. For the equal flow collection, however, the nonlinearity is mitigated, and \i appears to be a linear function of displacement. For /3 = 1, the both trajectory collections exhibit the same linear behavior. Thus positive correlation appears to mitigate the effect of injection mode on the mean behavior. The reaction path variance is strongly affected by the sign of the correlation (see Figures 4.5 and 4.6). For negative correlation, there is a distinct difference between the injection mode. However, positive correlation again mitigates the difference in modes. Moreover, the positive correlation greatly reduces the reaction path variance. For linear equilibrium sorption, this implies a reduction of the reactive travel time variance. The next case considered is that of no correlation between the reaction parameter and the log conductivity (/? = 0). The reaction path mean exhibits the same behavior as does nonreactive travel time. The nonlinearity in the mean for the equal area streamtubes is attributable to the nonstationarity in the associated velocity field. However, the combined effects of heterogeneous velocity and sorption manifest themselves in the reaction flow path variance. The reaction flow path variance is higher than that of the nonreactive travel time, even though the scaled means are the same. The /3  0 reaction flow path variance "lies in between" that of the positive and negative correlation. For the lower "total system variability" case of af nk = 1/2, the magnitude of the effect of correlation appears nearly symmetric (see Figure 4.6). For cinfc  1 5 however, there the effect of negative correlation seems to have a greater impact upon both total variability and the effect of injection mode (see Figure 4.5).

PAGE 104

89 Figure 4.3: Mean reaction path observed values and estimators for different values of correlation parameter ft and different trajectory collections for the crf nk = 1 set of simulations. Top figure is ft = {-1, +1}. Middle figure is /3 = {-1/2, +1/2}. Bottom figure is ft = 0 (no correlation). Trajectory collections are areaweighted Lagrangian (AL) and flowweighted Lagrangian (AF).

PAGE 105

90 0 10 20 30 40 0 10 20 30 40 l/A|n* o J£ 1 1 1 0 10 20 30 40 i/Aijt Figure 4.4: Mean reaction path observed values and estimators for different values of correlation parameter and different trajectory collections for the af nk = 1/2 set of simulations. Top figure is P = {-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).

PAGE 106

91 Figure 4.5: Reaction path variance observed values and estimators for different values of correlation parameter @ and different trajectory collections for the af nk = 1 set of simulations. Top figure is /3 = { 1,+1}. Bottom figure is /3 = {1/2, +1/2}. Trajectory collections are area-weighted Lagrangian (AL) and flow-weighted Lagrangian (AF).

PAGE 107

92 Figure 4.6: Reaction path variance values and estimators for different values of correlation parameter /3 and different trajectory collections for the af nk = 1/2 set of simulations. Top figure is /3 = {-1,+1}. Bottom figure is /3 = {-1/2,-1-1/2}. Trajectory collections are area-weighted Lagrangian (AL) and flow-weighted Lagrangian (AF).

PAGE 108

93 4.4.3 Reaction Flow Path-Travel Time Cross Moment The reaction flow path-nonreactive travel time cross moment a^ T contributes to the variance of the expected solute flux (see Equations (4.11) and (4.13)). For linear equilibrium sorption, this reflects the correlation between the time sorbed and the time traveling. The correlation of the reaction flow path and the travel time is little affected by injection mode (see Figures 4.7 and 4.8). It is, however, sensitive to the correlation parameter @. As with the reaction flow path statistics, the deviation of from the /? = 0 is greater in the case of positive correlation for the same value of /?. The "total system variability" seems to diminish the effect of correlation upon a MT This is most readily apparent from a comparison of the two ft = {1,1} cases for the different values of af nk (compare the top figures in Figures 4.7 and 4.8). The correlations for the correlated case "shift" towards that of the /3  0 in the higher log conductivity variance field, and the f3 = 0 case is slightly lower for the greater af nk In all cases except /3 = l,
PAGE 109

94 -0.2 1 1 1 0 10 20 30 40 Figure 4.7: Correlation between the reaction flow path and travel time for different values of correlation parameter /? and different trajectory collections for the of nJfc = 1 set of simulations. Top figure is j3 = {-1,+1}. Bottom figure is (3 = {-1/2, +1/2}. Trajectory collections are area-weighted Lagrangian (AL) and flow-weighted Lagrangian (AF).

PAGE 110

95 0.8 0.6 0.4 0.2 -0.2 -AL -FL x OAL OFL a +AL +FL o 10 20 30 40 b b 3. b 0.8 06 0.4 0.2 -0.2 ^ '' -AL -FL x OAL OFL +AL +FL 10 20 30 40 Figure 4.8: Correlation between the reaction flow path and travel time for different values of correlation parameter /3 and different trajectory collections for the a \nk  1/2 set of simulations. Top figure is /3 = { 1,+1}. Bottom figure is ft  {1/2, +1/2}. Trajectory collections are area-weighted Lagrangian (AL) and flowweighted Lagrangian (AF).

PAGE 111

96 conductivity in the calculation of the mean, working from Lagrangian properties that are stationary in displacements along the mean flow direction, and examination of the effects of injection mode. 4.5.1 Mean Reaction Flow Path Substitution of Equation (4.4) into Equation (4.2) and taking the expected value yields < The stationary Lagrangian serves as a reference field. The log conductivity field is decomposed into its mean and mean removed perturbation / and the velocity into its mean and mean removed perturbation thus < Define P g = P* exp[/xi nfc ], and expand the nonlinear terms of the previous equation into series, take the expected value retaining second order terms thus = P exp ^ /2] (1 fa, + fo}/2 + al)x (4.16) Recall that aj = af nk The / subscript is used for clarity. This expression captures the linearity of the observed mean \x (see Figure 4.3), but shows a heavy dependence upon the sign of the correlation parameter f3 that is not apparent in the data. The mean reaction flow path for the equal flow streamtubes appears to be dependent only upon the mean value of the reaction parameter, e.g. <>= P >^ Kl l+^>l\ +
PAGE 112

97 For the hypothesized Lagrangian distribution V = U(l + baj) and a\ = baj, thus the Equation reduces to <>= P g ex V [a 2 j2 + P 2 a 2 /2]x/U (4.18) The equal area streamtube mean [i is = p;e Xp[4 /2,,f(2^)^ (4,9) Manipulations similar to that which yielded Equation (4.16) yield < Uv/v[0] >= ^^Hf f (1 + 2a 2 + P 2 o 2 /2 + C v [x] 2(3C vf [x})dx (4.20) 1/(1+00" j) J Q Substituting the simplified relationship C vf [x] = a 2 exp[-a;/(2A)]/2 yields 40^< 1+2 ^ + ^ 2) (4 21) + 2/3AaJ(exp[-x/(2A)] 1) + \o){\ exp[-ox/A]) Notice that the first collection of terms, or the linear portion, of the preceding equation is the same as that in Equation (4.16). Similar arguments are employed to yield the equal area streamtube mean // estimator P g exp[^/2 + /?V;/2] < U,/v[0] >= f x (4 22) + 2^AaJ(exp[-x/(2A)] 1) + Xa){\ exp[-6x/A]) 4.5.2 Reaction Flow Path Variance Adapting an expression for the reaction flow path variance given by Cvetkovic et al. [1998] (their Equation (B4)) to the simplified velocity and velocity-log conductivity functions central to this work, write al\x] = 2P 2 exp[a 2 w }/U 2 / (x s)(C w [s] + p 2 C f [s] 4pC fv [s] + C v [s})ds (4.23) Jo

PAGE 113

98 For the special case of a\ = (4.24) + (2-4P)x/\ f } A simple and intuitive reaction flow path variance estimator for the equal-area streamtube collection is untenable in this framework. This is due in part to the nonlinear reaction flow path mean. Figures 4.3 and 4.4 demonstrate a close match between Equations (4.18) and (4.22) and simulation data. The parameters placed into the estimators are the target values selected for the simulation, and thus represent an a priori estimate of the mean reaction flow path. Equation (4.22) captures not only the nonlinearity of the nonreactive travel time, but the nonlinearities induced by reaction parameter and log conductivity correlation. While the "absolute difference" between Equations (4.18) and (4.22) appear small, consider that the scale of the simulations is fairly large. The solute has traveled over forty correlation scales of the log conductivity. In settings where the the scale of conductivity heterogeneity are large compared to the overall scale, the relative impact of injection mode and reaction parameter correlation to the log conductivity will have a relatively greater effect. The equal-flux streamtube reaction flow path variance estimator (Equation (4.24)) shows relatively poor performance in the more heterogeneous system (Figure 4.5), compared to that in the less heterogeneous system (Figure 4.6). Moreover, this expression does not work well at all for values of positive correlation, yielding physically implausible reaction flow path variances (cr^ < 0) for correlation parameter values greater that 1/2. The lack of "robustness" is perhaps attributable to an inadequate 4.6 Discussion

PAGE 114

99 representation of the myriad interacting processes by the first-order perturbation expansions. Injection of solute in the influent fluid mitigates the effect of correlation on the mean and the nonlinearity of mean nonreactive travel time. For tracer experiments that determine the volume or area of some property based upon mean arrival times of reactive and nonreactive tracers, correlation will have little effect upon the first temporal moment of the solute breakthrough curve if and only if the solute is uniformly injected in flux. However, injections that preferentially favor low-flux areas may exhibit strong nonlinearities, and be particularly sensitive to negative correlations between the local conductivity and reaction parameter. The contribution of the fj, variance to the second moment of the reactive solute travel time is extremely sensitive to correlation between the reaction parameter and the conductivity. Experiments that rely upon higher order moments of tracer experiments are, therefore, questionable if there is no p-ln k correlation information. Whereas injection mode dominates the correlation behavior in the mean, correlation dominates injection mode in the reactive travel time variance. 4.7 Spatial Moments Symmetry would dictate proceeding with a corresponding spatial moment analysis for transport of solute subject to a linear equilibrium sorption in a heterogeneous reaction parameter field. However, this process presents a particular complication that thwarts the development of relatively simple expressions following the same methods those previously employed. In the case of the travel time statistics, the "control structures" are two control planes. A fluid and solute parcel travels between these two control planes, sweeping out a fixed and constant volume, without regard to correlation properties, reaction parameters, or processes. The time required to sweep this volume, of course, is highly dependent upon these factors. For a "clock time"

PAGE 115

100 control structure, a solute parcel will trace different volumes for different processes and correlation structures. The reaction parameter model is comprised of two random processes. The first is that of the log conductivity field and the second is a random component that has the same distributional properties as, but is otherwise uncorrected to, the log conductivity field. In the case of the temporal moments, the weighted covariance functions of the log conductivity and this random component were employed directly. These functions can be parameterized, of course, in terms of the nonreactive travel time. However, these travel time parameterized covariance functions time are not immediately applicable to sorption analysis in the way that covariance parameterized in clock time would be, because the displacement for a given clock time is highly dependent upon the correlation of the reaction parameter to the log conductivity field. This is not to say that sorption-dependent spatial moment analysis is not viable in this framework. In fact, Cvetkovic et al. [1998] analyzed the spatial moments of a solute subject to nonlinear equilibrium sorption in a heterogeneous sorption field in this framework. However, the sorption parameter and the underlying log conductivity field were uncorrelated. Further clock time oriented are developments left for future work. Expressions involving the travel time dependent reaction flow path are developed and applications to and implications for first order decay are discussed. Continuing the theme of this dissertation, expressions based upon the trajectorybased covariance functions are sought. The time-dependent reaction flow path for the selected reaction parameter model is related to the log conductivity covariance function, as shall be demonstrated in the following. As in the previous chapter, transport along Lagrangian trajectories is considered. Assumed are second-order stationary velocities along area-weighted trajectories in time and second-order stationary velocities along flow-weighted trajectories in space.

PAGE 116

101 The log hydraulic conductivity is assumed to be similarly stationary. Moreover, the point random process P g exp[w] is assumed likewise stationary, since it is uncorrelated to the random log conductivity. From these assumptions, the expected value of the time-dependent reaction flow path is < h[t] > = ^ < P g exp[w[t}] exp[p In k[t}} > dt = P g exp[a 2 j2 + /3 2 al k /2}r The variance is given by var[/x[r]] =< f P g exp [w[t]] exp [0 In k[t}} dt Jo (4-26) / P g exp [w[t\] exp [P In k[t]} dt'< h[t] > 2 J o Consider the integrals in the previous equation. The exact expression is denoted by I and is approximated by the truncated series expansions of the exponential terms thus J ~ p g f [ T < {I + w[t}) {I + w[t']) (1+P\n k[t}) (1+ 13 In k[t'])>dtdt' (4.27) Jo Jo Note that a possibly more elegant approach would be to explore the Lagrangian features of the conductivity field, rather than the log conductivity field, and work from a stationary conductivity field covariance function and a similar uncorrelated field. This would obviate the need for the expansions. This suggested inquiry is left to future works. Expanding the above expression, taking the expected value, and dropping products of covariances yields I = P 2 g f f (1 + C w [t,t'} + (3 2 C lnk [t,t'})dtdt' (4.28) 7o Jo A few observations allow the adoption of an approximate expression for the timedependent Lagrangian covariance function C\ n k[t, t']. Relatively simple expressions that capture the behavior of the dominant processes are sought here. These expressions might facilitate future works that deemphasize simplicity for greater rigor and

PAGE 117

102 perhaps a better fit to the data. The Eulerian log conductivity covariance is the foundation upon which the rest of this work has been built. The model field is exponentially correlated in space, with correlation length Ai nfc The Eulerian velocity covariance exhibits a much longer correlation length due to the dissipation of the head across the field (see previous chapters). It is reasonable, therefore, to expect that the Lagrangian log conductivity covariance correlation length will lie somewhere in between that of the Eulerian log conductivity and the Lagrangian velocity. At zero lag, the area-weighted Lagrangian-trajectory log conductivity covariance is equal to the variance of the Eulerian field. This value is af nfc by definition. Thus, an exponential approximation of the Lagrangian log conductivity covariance function is C lDk i at [t] = al k exv[-b'Ut/\ lnk } (4.29) where b' is a function similar to b that relates the correlation length to a correlation time. The approximation b'  1 is made, which is equivalent to assuming that the Lagrangian and Eulerian correlation lengths are equal is space and substituting x = Vht = Ut into the space-dependent Eulerian log conductivity covariance function. Comparison of Equation (4.29) with b' = 1 to covariance functions estimated from log conductivity observations taken along area-weighted Lagrangian-trajectories generated in the numerical experiments described in Chapter 2 shows good qualitative agreement (see Figure 4.9). The data exhibit an apparent dependence of the Lagrangian correlation length upon the variance of the log conductivity field that is not captured by Equation (4. 29). One might expect similar behavior in the Eulerian field a priori, since the appropriate substitution into the Eulerian covariance function would be x = Uht = Ut/(1 + baf nk ). The higher variance of log conductivity results in a lower harmonic mean average velocity, implying it takes longer to to travel a unit length in more heterogeneous fields and a larger Lagrangian correlation time. In

PAGE 118

103 0.8 0.6 c o -£ 0.4 o u 0.2 1/2 1 Lagrangian In k Lagrangian V X -0.2 *++ X + + x xx xx x x xx,. 4 6 Ut/X lnk 10 Figure 4.9: Assumed time-dependent 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). fact, this is the behavior observed in the simulations (see Figure 4.10) The data are compelling enough to suggest further work. However, we shall accept Equation (4.29) with b' = 1 as a suitable approximation. Substituting Equation (4.29) with t = \t t'\ into Equation (4.28), integrating and inserting the result into Equation (4.26) yields var[/i[r]] = U 2 iUT/K k -(l-exp[-Ur/X lnk }) -P?exp[al + p 2 ol k ]T 2 (4.30) Neither this equation nor Equation (4.25) show a dependence upon the sign of (3. For the mean reaction flow path, the simulations show a weak dependence upon correlation parameter (see Figure 4.11). The linear propagation of the mean reaction flow path in time is due to the stationarity of the local reaction parameter along the area-weighted Lagrangian trajectories. Again, it should be noted that "time"

PAGE 119

104 0.8 1/2 1 Lagrangian In k Eulerian In k 1/2 Eulerian In k 1 c o o u 0.6 0.4 0.2 0 +++++ + ++ .0.2 1 1 1 1 0 2 4 6 8 10 Figure 4.10: Hypothesized time-dependent Eulerian log conductivity covariance functions for of nJfc = {1/2, 1} compared to simulation data and Equation (4.29). Figure 4.11: Mean reaction flow path as a function of nonreactive travel time compared to estimator Equation (4.25) for correlation parameter values /5 = {  1,0, 1}.

PAGE 120

105 refers to nonreactive travel time. Reaction flow path variance, however, is strongly affected by correlation, and Equation (4.30) is completely inappropriate (see Figure 4.12). Moreover, this equation can give physically implausible values at large 60 50 40 20 10 0 0 2 4 6 8 10 Ut/\ lnk Figure 4.12: Reaction flow path variance a function of nonreactive travel time for correlation parameter values /3 = {1,0, 1}. values of time. This is a result of the series approximations required for the "raw second moment" term in the variance, and the lack of such approximations in the "first moment squared" term. 4.7.1 Injection Mode The coupled effects of injection mode and correlation on the mean travel time dependent reaction flow path are examined. Using the relationships developed in Chapter 3, the hypothesis that the statistics reaction flow path corresponding to the equal-flow streamtubes are given by flow-weighting the equal-area properties is made.

PAGE 121

106 Thus, the mean reaction flow path for a flux injection is given by < u frtT)>=< u f[pm> (4.3i) In taking the expected value, the cross-correlation of the initial velocity and the that portion of the reaction parameter correlated to the log conductivity field must be considered. Using now-familiar expansions, the expected value of the flowweighted reaction flow path is < ^W] >= ^^^^ V + PC fu [t]/U)dt (4.32) U Jo where C fu is the covariance of the initial velocity perturbation and the log conductivity perturbation. Once again, this covariance shall be approximated with an exponential function of the form (see Appendix B) Ut (4.33) Substituting this term into Equation (4.32) and integrating yields < f ,H >= 5^^! (ur + ^.expI-J^]) (4.34) 4.7.2 Discussion Similar to its space dependent counterpart, the mean time dependent reaction flow path is sensitive to injection mode. However, it is the uniform resident injection that appears to mitigate correlation effects upon the mean reaction flow path. Correlation has a strong effect of [i variance, and it is quite likely that these effects dominate any injection mode effects, as was the case for the space dependent \x.

PAGE 122

107 10 A a. v +1/2 0 -1/2 est +1/2 est 0 est -1/2 ..-'+ sk i .+ -x W* At* ,r A x. >, A J* 7^7 .+ .-x 4 X .+ X .+ x ,+ * *' A J? J* + X X,4*** /Si*--' £XH4 6 10 10 A a. V X/' /+ X' /+J* *' 4 6 Ut/X ]nk ^ -X ,-+ .

PAGE 123

CHAPTER 5 APPLICATION OF CONCEPTS 5.1 Introduction An illustration of a practical application of trajectory based flow and transport modeling is sought. Useful and fairly robust relationships, namely predictors of spatial moments of a continuously injected plume subject to a first order decay, using simpler streamtube concepts (e.g. Jury and Roth [1990]). However, the observed bias in prediction shall be explained using the more complete trajectory based approach. In light of the robustness of the simple relationships developed here, they are used to estimate relative degradation rates for different solute species and applied to field data collected at a field-scale macrodispersion and natural attenuation experiment conducted in Columbus, Mississippi (see Stauffer et al. [1994] and Stauffer et al. [1997]). 5.2 Theory Consider a subsurface source of groundwater contamination that releases solute into the aqueous phase at a constant concentration. An example of such a source is a poorly mobile non-aqueous phase liquid, comprised of relatively insoluble organic compounds, distributed at some near-residual saturation through out a portion of some aquifer. The "constant" assumption applies when the rate of change of a property of interest, such as the aqueous phase concentration at the source, is small in comparison to the temporal scales associated with transport and decay. The actual groundwater flow field is modeld as a steady irrotational flow field that may be resolved into an aggregation of streamtubes. Neglecting the effect of local 108

PAGE 124

109 dispersive processes, an advection-decay equation with streamtube specific parameters is assumed to describe the fate and transport of the solute along the streamtube trajectory. This advection-decay equation is S-*  subject to c[x, t < 0] = 0 and c[x = 0, t => 0] = Cq, where v is an effective advective velocity [LT _1 ], and A: is a first-order decay constant [T -1 ]. The assumption of effective streamtube parameters is the simplification that separates this method from the more general consideration of properties varying along the trajectory. Although much of this language is similar to that used in earlier chapters, the fundamental difference is that the streamtube has constant effective parameters. This is analgous to an assuption of stationarity of properties, but is different in significant ways, as shall be shown. For advective transport in a single streamtube, the resident concentration is equal to the flux-averaged concentration. However, the difference between these detection modes manifests itself when concentrations are averaged over several streamtubes. For the remainder of this discussion, any reference to a concentration will implicitly imply a reference to a resident concentration. The solution to this equation for the given initial and boundary conditions is  = exp [ kx/v] h[vt  x] (5.2) Co where h[x] is the Heaviside step function, with the properties h[x] = 0 for x < 0 and h[x] = 1 for x > 0. 5.2.1 Spatial Moments The defintion of the raw spatial moment of order n of curve c is fj! n = J oox n c[a;](ix (5.3)

PAGE 125

110 A normalized moment of order n is a raw moment of order n divided by the zeroth moment, and is denoted by dropping the prime (i.e., /i n ). The zeroth spatial moment of this concentration curve is given by ^[t] = l[l exp[-fct]] (5.4) The zeroth spatial moment of a contaminant plume related to the total mass in solution. The first normalized spatial moment, or center of mass, is given by rtn 1 exp[-kt](kt + 1)  = 1* J 1 exp[-!fct] ( ] In the limit t  > oo n\  v/k. The concentration given in Equation (5.2) is a function of this asymptotic center of mass v/k. Solving Equation (5.2) with h[-] = 1 for v/k yields (5.6) k In [c/co] Thus, an observation of the concentration at any point of a steady streamtube concentration curve is an observation of its center of mass, given that cq and x are known. Assume that independent samples can be drawn uniformly from the collection of streamtubes that originate from a source Cq located at x  0. From the linearity of first moments, the mean center of mass may be estimated from N samples by the simple estimator 1 N x*"|>mT^ (5 ?) 5.2.2 Relative Degradation Consider two solutes a and b subject to a first-order decay process, characterized by rate constants A; a and k b respectively. The relative degradation rate of solute a to

PAGE 126

Ill b is defined as K = k a /k b (5.8) Let these solutes be released into a flow field at a steady rate, and allow the solutes to reach a "steady state" in which the concentration profiles for each do not vary in time. If within a streamtube i, the streamtube effective velocity Vi is a approximately equal for solutes a and b, the relative degradation rate of a to b in streamtube i is given by kri = (vi/k b )/(vi/k a ) (5.9) That is, the relative degradation rate is given by the ratio of the first normalized spatial moment of the steady-state concentration curve of solute b to that of a. The average relative degradation rate is estimated as the arithmetic mean of k ri Notice that this is not the ratio of the mean first spatial moments. 5.3 Application: Simulation The estimator was tested by simulating the release of a solute into a heterogeneous two-dimensional aquifer using the United States Geological Survey flow and transport code MOC3D [Konikow et al, 1996]. The length units used herein are a dimensionless quantity resulting from the normalization of the fundamental length by that of the correlation length of the exponentially correlated log hydraulic conductivity field. In the interest of clarity, however, this dimensionless quantity will be referred to as a correlation length A. Let the reference Cartesian coordinate system be oriented such that the ^-direction is parallel with the mean direction of flow, and the y-direction is orthogonal to that and lies within the plane formed by the aquifer. Let the origin lie at the lower left corner of the aquifer. The rectangular aquifer extended 15 A x and 25 A in y, and was

PAGE 127

112 subdivided into 10 nodes per A in both direction. Constant heads maintained at the x = OA and x = 15A boundaries imposed a steady mean gradient of J = 0.01. No-flow conditions were maintained at the y = OA and y = 25A boundaries. To simulate the release of a contaminant at constant concentration, the influent water concentration was assumed to be a constant concentration c 0 = 1. This "source" was centered on the x = 0 boundary, and was 10A wide. A turning-bands algorithm generated a synthetic heterogeneous log conductivity field used in the simulations [Tompson et ai, 1989]. The output of the turning-bands generator was a realization of an exponentially-correlated standard normal process (/x = 0,cr = 1). This field was subsequently converted to fields with variances of 0.1, 0.5, and 1.0 via the transformation / = VlnkZ (5.10) where a\ a is the target standard deviation of the log conductivity field / and z is the standard normal process. The aquifer porosity was assumed to be equal to the mobile water content and equal to a constant value n = 0.2. The effective conductivity of a two-dimensional exponentially correlated log conductivity field is, to a good approximation, the geometric mean conductivity K g For \x  0, K g = 1, in units of (dimensionless) correlation lengths per characteristic time. The independence of the effective conductivity from its variance convieniently allows the variance to be changed without substantially changing the bulk flow properties. Thus, the anticipated Darcy flux for the aquifer, regardless of o\ nk was = 0.01A per characteristic time, and the characteristic filtration velocity was vi, = 0.05A per characteristic time. The subscript b refers to a bulk property, to be distinguished from streamtube properties discussed earlier. The solute was assumed to undergo a constant decay process characterized by a first-order rate coefficient k  0.05 in units of an inverse characteristic time defined

PAGE 128

113 by Vf/X. This value yields what might be termed a "bulk" first spatial moment Vb/k = 1. For a homogeneous conductivity field, this is numerically equivalent to the first longitunial spatial moment of the steady-state plume. The steady-state concentration field was "sampled" by taking the reported numerical concentration values as the local resident concentration. The plume centerof-mass was estimated by the estimator presented in Equation (5.7). The "actual" first longitudinal spatial moment was calculated using a trapezoidal rule integration. A fate and transport simulations were carried out in 6 sets of hydraulic conductivity fields corresponding to 6 standard normal process realizations. 5.3.1 Results For the case of a homogeneous conductivity field, any one non-zero concentration measurement used with the estimator given in Equation (5.7) returned a value for the first longitudinal spatial moment that for all practical purposes was equal to Vb/k, and five percent less than that predicted by the numerical integration. As the variance of log conductivity increased, however, both the numerical integration and the estimator tended to systematically predict values greater than Vb/k. Additionally, the estimator systematically underpredicted the value given by the numerical integration, and by an amount which generally increased with variance of the log conductivity. 5.3.2 Discussion This systematic increase in the the first longitudinal spatial moment of the plume with increasing log hydraulic condutivity variance, despite a fairly constant v b is explained by considering the nature of the underlying Lagrianian velocity field. The steady-state concentration plume is, in essence, a map of travel times from the source. As previously stated, Cvetkovic et al. [1996] demonstrated that the statisics of such a field are nonstationary. For small displacements from the source, the mean travel

PAGE 129

114 Figure 5.1: Resident concentration profile averaged along a transect orthogonal to the mean flow direction for plumes for cr, 2 nfc = {0.0, 0.1, 0.5, 1.0}. time propagates dr/dx = l/v h (5.11) where Vh is the harmonic mean velocity. For regions distant from the source region, the mean travel time propagates dr/dx = l/v a (5.12) where v a is the arithmetic mean velocity, and equal in value to Vb, the bulk filtration velocity. The harmonic mean velocity is necessarily less than or equal to the arithmetic mean velocity, which in turn indicates that the mean travel time propagates as a slower rate near the source. Additionally, the correlation of the Lagrangian velocity persists over greater distances than that of the underlying conductivity field. Solute which is introduced in

PAGE 130

115 areas of high local flux tend to move away from the source at relatively high velocities. The net effect of what might be termed a "preferential flow" of solute in these higher velocity streamtubes, is to induce a "tail" on a transversely-averaged resident concentration profile, which increases with variance of log conductivity. These effects are illustrated by comparing concentration profile created by averaging concentrations along planes normal to the mean direction of flow to concentration profiles given by c/c 0 = exp[-kx/v] and c = /c 0 exp[-/cr[:r]], where r is given by t[x] =  (x + \inklk (! exp[-te/Ai nfc ])) (5.13) "a Here Ai n is the correlation length of the log conductivity distribution and b is a shape factor for the conductivity anisotropy. 5.4 Application: Field Data Data from an elaborate experiment conducted at Columbus Air Force Base, Mississippi were analyzed using the results dervied in the previous section. The objective the experiment was to characterize the natural attenuation of certain hydrocarbons in groundwater emanating from subsurface nonaqueous phase liquid (NAPL) sources Stauffer et al. [1997]. A NAPL hydrocarbon source comprised of decane, benzene, toluene, ethyl benzene, para xylene (p-xylene), and naphthalene was emplaced at the site of the MADE [Boggs et a/., 1992] and MADE-2 [Maclntyre et al, 1993] experiments. This source released these poorly-soluble organic solutes into the groundwater, and the resultant plumes were monitored at down-gradient multilevel sampler locations. The center of mass of the steady-state plume was estimated using a trapezoidal rule integration of concentration data in space and by a simple averaging of pointwise observations of concentration via Equation (5.7). Relative degradation and absolute degradation rates of the different constituents are estimated using Equation (5.9).

PAGE 131

116 5.4.1 Results and Discussion: Center of Mass The center of mass of the steady state plumes emanating from the emplaced contaminant source was estimated by two methods. The first might be considered a traditional approach, namely simple trapezoidal rule integrations of the data. A known initial concentration for each species was assumed known. This concentration was calculated from initial source mass fractions and aqueous solubilities given in Stauffer et al. [1997] and assuming Raoult's law (see Table 5.1). The second approach was to use Equation (5.7) (see Table 5.1). The same initial concentration used in the spatial integration was used with the estimator approach. Displacements were taken as the distance from the multilevel sampler to the source zone. benzene toluene ethyl benzene p-xylene naphthalene integration 4l) 21) 3l 2^3 3^8 v/k 5.3 1.7 2.5 1.9 3.1 Table 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). Equation (5.7) gives results similar to those of the trapezoidal rule integration, while offering some compelling computational advantages. While the spatial integration is conceptually simple, it is computationally demanding for unequally spaced data in that the different "trapezoids" must be calculated and summed for each integration. New data require recalculation of these trapezoids. Of course, these calculations are trivial for computer codes, but tedious for "back of the envelope" calculations. Equation (5.7) is a fairly simple calculation for any scientific calculator, and provides an excellent tool for a rapid estimate of the center of mass from sparse and scattered local concentration. New data are easily incorporated into existing sets without the need of a global recalculation.

PAGE 132

117 5.4.2 Results and Discussion: Degradation The relative degradation and absolute degradation of the different contaminants were estimated using Equation (5.9) and degradation rates reported by Maclntyre et al. [1993]. These reported benzene, p-xylene, and naphthalene degradation rates were estimated from another experiment conducted at the site in which a "cocktail" of hydrocarbons was injected into the groundwater and the fate of the resultant plume was monitored (see Table 5.2). This experiment was conducted before the NAPL benzene toluene ethyl benzene p-xylene relative 087 3A L5 2l) estimated 0.0056 0.022 0.0094 0.013 reported 0.0070 NR NR 0.011 Table 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 source was emplaced, and is in that sense independent of the concentration data associated with the natural attenuation experiment. Never-the-less, the rates estimated from the natural attenuation data were similar to those reported by Maclntyre et al. [1993]. The technique used here is quite simple, and does not require regressions or nonlinear fits to transport models. 5.5 General Conclusions These estimators, while fairly simple and apparently robust, are limited in applicability to quasi-steady plumes, and it is not clear the extent to which such plumes exist. Sites where poorly mobile NAPLs have been known to exist for long times are likely candidates, especially if decay rates are known to be large in comparsion to transport time scales.

PAGE 133

118 Local dispersion was neglected, and local dispersion certainly has a significant effect on the local concentration variability. However, its neglect does not appear to diminish the results here to the point of uselessness. The nonstationarity of equal-area Lagrangian trajectories in a steady flow field was neglected. For highly heterogeneous conductivity fields, low velocity areas are disproportionately sampled with respect to the flow average velocity. For a spatially uniform decay process, this would tend to underestimate the plume extent, since the "observed values" will tend to be associated with long travel times. It should be noted that this sampling bias is intrinsic to multilevel samplers. The correlation of velocity and decay coefficient was neglected. A strong negative correlation could result in a significant underprediction of the center of mass in a highly heterogeneous field. Again, it should be noted that this is partially due to the sampling bias assumed by the equal-weight multilevel sampler sampling strategy that will preferentially sample low velocity areas with respect to the flow weighted mean.

PAGE 134

CHAPTER 6 GENERAL CONCLUSIONS Any discipline in which you can publish a paper on boundary conditions is an immature discipline. Sten Berglund 1 Many of the results of this study are quantitatively and qualitatively the same as those of previous studies. For instance, chemical engineers have long known that the "mean residence time" of the "swept volume" of a reactor for a flux injection is L/U, where L is the reactor length and U is the mean filtration velocity. However, hydrologists seek to understand the effects of the heterogeneity of the "reactor" itself and the interplay of the processes as they manifest themselves "between x = 0 and x = L." An old premise served as the conceptual foundation of this work, namely that a system, in this case a heterogeneous flow field, may be observed in two ways: along "straight-line," or Eulerian, trajectories and along "natural" or Lagrangian trajectories along hydrodynamic streamlines. The flow field was decomposed into an aggregation of elements, either as area elements that lie between Eulerian trajectories or hydrodynamic streamtubes that lie in between streamlines. Two criteria defined four different collections of the flow field elements, namely that at some definition plane, each element has either an equal area through which water enters the element, or an equal volume of water enters the element per unit time. The imaginative 1 From a conversation with the author about the acceptance of Demmy et al. [1999] for publication. 119

PAGE 135

120 reader might contrive a plethora of trajectory and trajectory weighting schemes. For compactness, the following notation summarizes the four different trajectories AW/ET equal-area weighted Eulerian, or straight line, trajectory AW/LT equal-area weighted Lagrangian, or streamline, trajectory FW/ET equal-flow weighted Eulerian, or straight line, trajectory FW/LT equal-flow weighted Lagrangian, or streamline, trajectory Myriad methods exist for mapping properties observed along the different trajectories. That is, the value of some property may be observed along the trajectory and recorded at different intervals in space or time. A traditional method is a regular Eulerian gridding that divides the flow field into equally-spaced trajectories and "records" information at regular intervals along the trajectory. Two specific parameterizations are considered: equal intervals in displacement along the mean flow direction and equal intervals in advective travel time. Employing these concepts, the behavior of different flow field properties was examined. These properties, such as log conductivity and local velocity, were "observed" along the different trajectories at "control planes," or planes equally spaced in the mean flow direction, or "control times," or at the displacement that an indivisible fluid parcel would be after traveling for some reference time. For a steady and irrotational flow field resulting from a constant mean gradient applied across a stationary Eulerian log conductivity field of uniform and constant mobile water content, the trajectories for which the statistics certain properties appeared to be second-order stationary were identified. That is, the mean and variance of the property appeared to have a constant value at all observation points along the trajectory. The log conductivity and velocity observed at control planes appeared stationary for the equallyspaced Eulerian trajectories (AW/ET) and Lagrangian trajectories separated by an

PAGE 136

121 equal discharge (FW/LT). For equally-spaced control times, only Lagrangian trajectories equally-spaced at a reference plane (AW/LT) appeared to exhibit stationary log conductivity and velocity statistics. The stationarity of the properties greatly simplifies the development of predictors of solute transport. When taking the expected value of the integral of some property along a trajectory, say the integral of velocity in time, the time-varying velocity function may be replaced with its time average. As implied by the stationarity of the statistics along some trajectory collections, but not others, the numerical statistics of properties varied between the different collections. Adopting simple and accurate models for the statistical properties of some property in a reference collection, the statistics of the properties along other collections, including the observed nonstationarities, were well predicted. This had significant application in the evaluation of the effect of injection mode upon solute transport. Two methods by which solute might be uniformly introduced into a system were considered. The first was a resident injection were a fixed volume at the system inlet was uniformly filled with solute at some reference concentration. The second was to allow influent water at some reference concentration to carry the mass into the system. The displacement and travel time statistics of AW/LT correspond to solute transport associated with a uniform resident injection. The transport of solute injected in flux is described by the displacement and travel time statistics of the FW/LT. The so-called consistent first-order approximations of the inverse velocity were shown to be built upon a flawed premise, but are none-the-less quantitatively correct and robust. The premise is that the Lagrangian flow field may be accurately approximated by the Eulerian flow field. While this may be true for the weakest of heterogeneities, this approximation breaks down rapidly with increasing heterogeneity if it is consistently applied. As inconsistently formulated, the approximation

PAGE 137

122 < 1/v >~ 1/U is quantitatively equal to the inverse of the harmonic mean of the Lagrangian velocity < 1/v >= l/V h Moreover, the time-average Lagrangian velocity is quantitatively equal to the space average Eulerian velocity, and these equalities contribute to the robustness of other approximations. Were the Lagrangian field to be consistently approximated by the Eulerian, travel times would be greatly overpredicted and displacements greatly under-predicted. It is, in fact, the consistent inconsistent approximation of the Lagrangian by the Eulerian that contributes to much of the observed robustness of the Lagrangian transport theory. While not many "new" equations were proposed or so many new phenomena were unearthed, it is hoped that a modicum of understanding was fostered. This is the primary contribution of this work. These concepts were illustrated with an analysis of the transport of an solute experiencing a linear equilibrium sorption in a field characterized by a variable sorption coefficient. While the process has been studied before, it is important to illustrate the new concepts and approaches here with familiar examples. The sorption study was restricted to transport as characterized by mass breakthrough at control planes. An interesting new result is that injection mode strongly affects the propagation of the mean breakthrough time. A uniform resident injection exhibits a nonlinearity in the mean, and this nonlinearity is predictably enhanced by negative correlation between the sorption coefficient and the underlying log conductivity field. Injection in flux, however, results in a linear propagation of mean arrival time in space, and mitigates the effect upon correlation. Arrival time variance, however, is dominated by the correlation between sorption coefficient and the log conductivity field. Injection mode has comparatively little effect. This has profound implications for tracer test analyses based upon the arrival time variance of breakthrough curves. The concepts developed in this work are further illustrated with the development and analysis of an estimator of the center of mass of a plume resulting from

PAGE 138

123 continuously injected solute subject to a first-order decay. The assumption of noninteracting streamtubes with constant "effective" properties allows solution of the advection-decay equation for a single streamtube. An ergodic-like hypothesis allows this solution to a distribution of streamtubes. This is what might be termed a "classic" streamtube approach. While surprisingly robust and perhaps worthy of application in field settings, the estimator shows a bias explained using the trajectory based analysis developed in this work. One measure of quality of a work is the number of questions that it answers. Perhaps another is the number of interesting questions that it begs. Several issues are left for further study. A few specific questions are addressed first, then a few more general. The relationship between an Eulerian "reference" system and its Lagrangian counterpart remains incompletely understood. A primary area that requires attention is the increasing Lagrangian velocity correlation length with increasing log conductivity variance. A great deal of work remains for predicting the spatial moments of a reactive solute. Numerical constraints limited this study to relatively mild heterogeneity. However, a highly heterogeneous medium, might be replaceable by two or more media characterized by lesser heterogeneity, and a systematic "replacement system" would greatly enhance computational efforts. The turning bands method employed in the generation of the random fields was computationally expensive, in comparison to, say, the particle tracking algorithms or data analysis programs, and to the flow solver for weak heterogeneity. The method requires the generation of several independent normal line processes. This problem is well-suited to parallel and distributed computational techniques, and an easy-to-use implementation would certainly find wide use. Several research groups are working on parallel and distributed flow codes, and the research community would benefit greatly from a wider dissemination of or greater access to these codes.

PAGE 139

124 Thinking more broadly, some of the concepts employed here might find application in other fields. The flow of water through a wetland system might be reasonably approximated as a two-dimensional aquifer-like object. Transport of contaminants might be modeled using similar techniques, especially if the transport time scales are relatively short compared to system variability. In general, heterogeneity will imply some sort of "preferential flow," and the identification of preferential flows in heterogeneous systems should be paramount consideration in the design of monitoring systems.

PAGE 140

APPENDIX A EULERIAN VELOCITY DISTRIBUTION PROPERTIES Consider lognormally distributed random variable V with moments < V >= U and var(V) = S. These moments are related to the parameters \i and a 2 via U = exp(At +
PAGE 141

126 The effective conductivity of a porous medium may be denned as the numerical conductivity value that satisfies the Darcy relationship Q = K e J (A.9) where Q is the specific discharge and J is the difference in hydraulic head across the formation. There are several expressions for the effective conductivity of a stationary heterogenous conductivity fields (see Gelhar [1993]). The general form of the effective conductivity is K e = K g f(al k ,g) (A.10) where K g and af nk are the geometric mean and the variance of the log conductivity field, respectively, g is a function of aquifer anisotropy, and / is a function that scales the geometric mean to the proper effective value. For a two-dimensional isotropic aquifer, f{cr 2 nk ) = 1. Thus, for an aquifer of homogeneous mobile water content 6, an estimate of the spatial mean seepage velocity U is given by U= I Y (A.11) The point velocity variance for such a field is S = bU 2 a 2 nk where b = 3/8. Substitution of Equation (A. 11) into this expression yields bK 2 J 2 o 2 S= 9 Q2 (A.12) Substitution of Equations (A. 11) and (A.12) into Equations (A. 7) and (A. 8) 0y/TTbal k a 2 = ln(l + ba 2 nk ) (AAA) I* = ]n \ j£sL= ] (A.13)

PAGE 142

APPENDIX B VELOCITY LOG CONDUCTIVITY CORRELATION A simplified correlation relationship between the velocity field and the log conductivity is sought. Graham and McLaughlin [1989b] and Rubin and Dagan [1992] derived expressions based upon perturbation expansions of Darcy's law and spectral methods. Graham and McLaughlin [1989b] assume a "hole-type" covariance function for the log conductivity field, in order to assure a complementary stationary head gradient field. Rubin and Dagan [1992] assume an exponentially correlated anisotropic log conductivity field. The velocity-log conductivity covariance functions for both analytical cases shows qualitatively similar behavior (see Figure B.l), and an simplified estimator would be applicable for both log conductivity correlation models. Although these rigorouslydeveloped estimators are attractive theoretically, the mundane and oft-neglected task of mapping theoretical results into application often requires simplifying assumptions. Presented for your consideration is the covariance function of Graham and McLaughlin [1989b] (their Equation (B5) where a = 7r/(4A/) and K { is the modified Bessel function of order i. The integral scale of this covariance function is A/. An exponential approximation to this function that reproduces the integral scale, general trend, and point covariance of that given by Graham and McLaughlin [1989b] is sought. A function of the form -(a/2) (C#i[aC]-aC 2 ffo[aC])} (B.l) 29 exp[-x/(2A)] (B.2) 127

PAGE 143

128 Figure B.l: 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 K g Jaf nk /9. The apparent "roughness" of the Rubin and Dagan [1992] function is due to estimation errors associated with inferring their function from graphical data.

PAGE 144

129 meets these criteria (see Figure B.l). The covariance function presented by Rubin and Dagan [1992] is the sum of three two-function products that involve integrations of Bessel functions. For twodimensional porous systems, the simplified function condenses that monumental work into an elegant and concise form. The simplified expression is tested against correlation functions estimated from the velocity and log conductivity data taken from the simulations (see Figure B.2). The simplified estimator gives a reasonable a priori estimate of the velocity and 1 1/2 x Eq. B.2 'xt + <* + x x + I I I 0 5 10 15 20 r/X } Figure B.2: Comparison of the simplified velocity-log conductivity correlation function to simulation data for aj = {1/2, 1}. log conductivity correlation observed along the equal-flux Lagrangian trajectories. The reasonable match to analytical expressions and the observed data enhance our confidence in this simplified estimator. The results of Rubin and Dagan [1992] (their Figure 7) indicate that anisotropy decreases the point covariance and increases the long range correlation effects. Yet 0.8 0.6 b ^ 0.4 b 0.2 _n o

PAGE 145

130 the functional form appears to remain exponential. Therefore, it is suggested that future work focus upon generalizing the simple results found here to three-dimensions and general anisotropy.

PAGE 146

APPENDIX C CONSISTENT FIRST ORDER APPROXIMATION C.l Travel Time Usage of a consistent first-order approximation for travel time permeates the stochastic Lagrangian literature. This approximation is quantitatively correct for collections of equal-flux streamtubes. However, the premise upon which it is built is incorrect. Consider a stationary lognormal process V. The lognormal distribution is selected for convenience, as a random variable with values guaranteed to be greater than zero is required. Since V is stationary by definition, then so is 1/V. Consider the expected value of the integral The order of evaluation may be interchanged, taking the expected value of the argument of the integral thus: where Vh is the point harmonic mean of the random process V. Consider now a perturbation expression for the random process V = V(l + v) where V is the mean of V, and v is a zero-mean random process representing the perturbation of the process about its mean. Inserting this expression into Equation (C.l) yields (C.l) DO (C.2) = x/V h (C.3) 131

PAGE 147

132 Expand (1 + v) in a series, 1 f = T j <{1-v[x] + v[x] 2 )>dx (C.4) v Jo Should this expression be truncated at first order, the leading term < r x/V is left. However, from Equation (C.2), the exact expression is < r >= x/V h < r >= x/V is a well-established result in the literature for the first temporal moment of a solute plume injected in flux Kreft and Zuber [1978]. The answer to this riddle lies in the examination of the properties of the flow field. Even if the deviations are relatively small from the straight-line, or Eulerian, trajectory, the statistics of the properties recorded along the actual, or Lagrangian, trajectory can be quite different. For instance, the mean velocity recorded along equal-flow weight streamlines in an quasi-infinite constant mean gradient flow in a exponentially correlated lognormal log conductivity field is approximately V = (1 + baf nk )U where U is the spatially averaged Eulerian velocity. The variance of the the velocity perturbation along these trajectories is approximately o\ = baf nk Evaluating the expected value in Equation (C.4) and retaining the second order term yields =t(l+ol) (C.5) Substituting our approximate values into the preceding equation yields \ l + h(J Lk) u (C.6) = x/U The so-called small-deviation assumption is a poor one. The approximation of the Lagrangian field by the Eulerian field works because the correlation properties are quite similar and the "coincidence" that the harmonic mean Lagrangian velocity is approximately equal to the arithmetic mean Eulerian velocity.

PAGE 148

133 C.2 Displacements The center of mass of a solute plume uniformly injected in space is known to move = Ut, where U, again, is the spatially averaged Eulerian velocity. If the solute parcels move along Lagrangian trajectories, and the statistics of Lagrangian and Eulerian fields are known to be quite different, why does this work? The answer is that the temporally averaged velocity of equal-area Lagrangian streamtubes is stationary in time and equal to the spatially averaged Eulerian velocity. In fact, the temporally averaged Eulerian velocity is lower than its spatially averaged counterpart. So once again, the small-displacement assumption, if consistently applied, would result in a very poor result. What is remarkable is that the correlation properties of the two fields are so similar, and herein lies the strength of the approximations of Lagrangian fields with Eulerian fields.

PAGE 149

REFERENCES Annable, M. D., J. W. Jawitz, P. S. C. Rao, D. P. Dai, H. K. Kim, and A. L. Wood. Field evaluation of interfacial and partitioning tracers for characterization of effective naplwater contact areas. Ground Water, 36(3):495-502, 1998. Bakr, A. A., L. W. Gelhar, A. L. Gutjahr, and J. R. MacMillan. Stochastic analysis of spatial variability in subsurface flows 1. Comparison of oneand three-dimensional flows. Water Resour. Res., 14(2):263-271, 1978. Bear, J. Dynamics of Fluids in Porous Media. Dover Publications, Inc., New York, 1972. Bellin, A., Y. Rubin, and A. Rinaldo. Eulerian-Lagrangian approach for modeling of flow and transport in heterogeneous geological formations. Water Resour. Res., 30 (ll):2913-2924, 1994. Bellin, A., P. Salandin, and A. Rinaldo. Simulation of dispersion in heterogeneous porous formations: Statistics, first-order theories, convergence of computations. Water Resour. Res., 28(9):2211-2227, 1992. Boggs, J. M., S. C. Young, L. M. Beard, L. W. Gelhar, K. R. Rehfeldt, and E. E. Adams. Field study of dispersion in a heterogeneous aquifer 1. Overview and site description. Water Resour. Res., 28(12):3281-3291, 1992. Cvetkovic, V. D., H. Cheng, and X.-H. Wen. Analysis of nonlinear effects on tracer migration in heterogeneous aquifers using Lagrangian travel time statistics. Water Resour. Res., 32(6):1671-1680, 1996. Cvetkovic, V., G. Dagan, and H. Cheng. Contaminant transport in aquifers with spatially variable hydraulic and sorption properties. Proc. R. Soc. London Ser. A, 454:2173-2207, 1998. Dagan, G. Stochastic modeling of groundwater flow by unconditional and conditional probabilites 1. Conditional simulation and the direct problem. Water Resour. Res., 18:835-848, 1982a. 134

PAGE 150

135 Dagan, G. Stochastic modeling of groundwater flow by unconditional and conditional probabilites 2. The solute transport. Water Resour. Res., 18:813-833, 1982b. Dagan, G. Solute transport in heterogeneous formations. J. Fluid Mech., 145:151177, 1984. Dagan, G. Flow and Transport in Porous Formations. SpringerVerlag, New York, 1989. Dagan, G. and V. Cvetkovic. Spatial moments of a kinetically sorbing solute plume in a heterogeneous aquifer. Water Resour. Res., 29(12):4053-4061, 1993. Dagan, G., V. Cvetkovic, and A. Shapiro. A solute flux approach to transport in heterogeneous formations 1. The general framework. Water Resour. Res., 28(5): 1369-1376, 1992. Demmy, G., S. Berglund, and W. Graham. Injection mode implications for nonreactive solute transport in porous media: Analysis in a stochastic Lagrangian framework. Water Resour. Res., 35:1965-1974, 1999. Domenico, P. A. and F. W. Schwartz. Physical and Chemical Hydrogeology. John Wiley and Sons, New York, 1990. Freeze, R. A. A stochastic-conceptual analysis of one-dimensional groundwater flow in nonuniform homogeneous media. Water Resour. Res., 11(5):725-741, 1975. Gelhar, L. W. Stochastic Subsurface Hydrology. Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1993. Gelhar, L. W. and C. L. Axness. Three-dimensional stochastic analysis of macrodispersion in aquifers. Water Resour. Res., 19(1):161 180, 1983. Graham, W. and D. McLaughlin. Stochastic analysis of nonstationary subsurface solute transport 1. Unconditional moments. Water Resour. Res., 25(2):215-232, 1989a. Graham, W. and D. McLaughlin. Stochastic analysis of nonstationary subsurface solute transport 2. Conditional moments. Water Resour. Res., 25(11):2331 2355, 1989b. James, A. I., W. Graham, K. Hatfield, P. S. C. Rao, and M. D. Annable. Optimal estimation of residual non-aqueous phase liquid saturations using partitioning tracer concentration data. Water Resour. Res., 33(12) :2621-2636, 1997.

PAGE 151

136 Jawitz, J. W. Aquifer contaminant source zone characterization with partitioning tracers and remediation with single-phase microemulsion flushing. PhD thesis, University of Florida, Gainesville, Florida, 1999. Jin, M., M. Delshad, V. Dwarakanath, D. C. McKinney, G. A. Pope, K. Sepehrnoori, and C. E. Tilburg. Partitioning tracer test for detection, estimation, and remediation performance assessment of subsurface nonaqueous phase liquids. Water Resour. Res., 31(5):1201-1211, 1995. Journel, A. and C. Huijbregts. Mining Geo statistics. Academic Press, New York, 1978. Jury, W. A. and K. Roth. Transfer Functions and Solute Movement Through Soil. Birkhauser Verlag, Basel, 1990. Konikow, L. F., D. J. Goode, and G. Z. Hornberger. A three-dimensional method-ofcharacteristics solute-transport model (MOC3D). Water-resources investigations report 96-4267, United States Geological Survey, Restion, Virginia, 1996. Kreft, A. and A. Zuber. On the physical meaning of the dispersion equation and its solutions for different initial and boundary conditions. Chem. Eng. Sci., 33: 1471-1480, 1978. Maclntyre, W. G., M. Boggs, C. Antworth, and T. B. Stauffer. Degradation kinetics of aromatic organic solutes introduced into a heterogeneous aquifer. Water Resour. Res., 29(12):4045-4051, 1993. Mood, A. M., F. A. Graybill, and D. C. Boes. Introduction to the Theory of Statistics. McGraw-Hill Publishing Company, New York, 1974. Mose, R., P. Siegel, P. Ackerer, and G. Chavent. Application of the mixed hybrid finite element approximation in a groundwater flow model: Luxury or necessity? Water Resour. Res., 30(11):3001-3012, 1994. Musashi, M. A Book of Five Rings. The Overlook Press, Woodstock, New York, 1982. Osnes, H. Stochastic analysis of velocity variability in bounded rectangular heterogeneous aquifers. Adv. Water Resour., 21:203-215, 1998. Papoulis, A. Probability, Random Variables, and Stochastic Processes. McGraw-Hill Publishing Company, New York, 1991. Pollock, D. W. Semianalytical computation of path lines for finite-difference models. Ground Water, 26(6):743-750, 1988.

PAGE 152

137 Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling. Numerical Recipes in C. Cambridge University Press, Cambridge, England, 1988. Rubin, Y. Stochastic modeling of macrodispersion in heterogeneous porous media. Water Resour. Res., 26(7):1689-1697, 1990. Rubin, Y. and G. Dagan. Stochastic analysis of boundaries effects on head spatial variability in heterogeneous aquifers 1. Constant head bounary. Water Resour. Res., 24(10):1689-1697, 1988. Rubin, Y. and G. Dagan. Stochastic analysis of boundaries effects on head spatial variability in heterogeneous aquifers 2. Impervious boundary. Water Resour. Res., 25(4):707-712, 1989. Rubin, Y. and G. Dagan. A note on the head and velocity convariances in threedimensional flow through heterogeneous anisotropic porous media. Water Resour. Res., 28(5):1463-1470, 1992. Shapiro, A. M. and V. D. Cvetkovic. Stochastic analysis of solute arrival time in heterogeneous porous media. Water Resour. Res., 24(10): 171 1-1718, 1988. Stauffer, T. B., C. P. Antworth, R. G. Young, W. G. Maclntyre, J. M. Boggs, and L. M. Beard. Degradation of aromatic hydrocarbons in an aquifer during a field experiment demonstrating the feasibility of remediation by natural attenuation. Technical report, Armstrong Laboratory, United States Air Force, Tyndall Air Force Base, Florida, 1994. Stauffer, T. B., J. M. Boggs, and W. G. Maclntyre. Ten years of research in groundwater transport studies at Columbus Air Force Base, Mississippi. In Sayler, G., J. Sanseverino, and K. L. Davis, editors, Biotechnology in the Sustainable Environment, Plenum Press, New York, 1997. Theile, M. R., R. P. Batycky, M. J. Blunt, and F. M. Orr, Jr. Simulating flow in heterogeneous systems using streamtubes and streamlines. SPERE, 11:5-12, 1996. Tompson, A. F. B. Numerical simulation of chemical migration in physically and chemically heterogeneous porous media. Water Resour. Res., 29(ll):3709-3726, 1993. Tompson, A. F. B., R. Ababou, and L. W. Gelhar. Implementation of the threedimensional turning bands random field generator. Water Resour. Res., 25(10): 2227-2243, 1989.

PAGE 153

BIOGRAPHICAL SKETCH George Demmy was born to George and Ellen Demmy 23 November 1966 in Lakeland, Florida. He received a Bachelor of Arts in Physics and German Literature from the Florida State University in 1990, a Master of Engineering in Agricultural Engineering from the University of Florida in 1993, and Doctor of Philosophy in Agricultural and Biological Engineering from the University of Florida in 1999. 138

PAGE 154

I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Wendy Professor of A^ricitltura^and Biological Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. P.S.C Rao/ CoChair Lee A. Rieth Distinguished Professor of Environmental Engineering Purdue University I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Kenneth L. Campbell Professor of Agricultural and Biological Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Kirk Hatfield Associate Professor of Civil Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Michael D. Annable Associate Professor of Environmental Engineering Sciences

PAGE 155

This dissertation was submitted to the Graduate Faculty of the College of Engineering and to the Graduate School and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy. December 1999 M. J. Ohanian Dean, College of Engineering Winfred M. Phillips Dean, Graduate School