UNIVERSITY OF MINNESOTA 
ST. ANTHONY FALLS HYDRAULIC LABORATORY 
Project Report No. 314 
HETEROGENEITY IN HYDRAULIC CONDUCTIVITY: 
ITS IMPACT UPON SOLUTE TRANSPORT 
AND CORRELATION STRUCTURES 
by 
Deane T. Morrison 
and 
Roko Andricevic 
LEGISLATIVE COMMISSION ON MINNESOTA RESOURCES 
State of Minnesota 
St. Paul, Minnesota 
June, 1991 
Minneapolls, Minnesota 
     
ACKNOWLEDGEMENT 
This study was completed with the support from the Legislative 
Commission on Minnesota Resources and computations are performed under a 
grant from the Minnesota Supercomputer Institute. We appreciate the 
comments and support from Dr. Efi Foufoula-Georgiou. Carolyn M. 
Wasikowski, from the Minnesota Supercomputer Center, produced a graphical 
visualization. 
The University of Minnesota is committed to the policy that all persons 
shall have equal access to its programs, facilities, and employment without 
regard to race, religion, color, sex, national origin, handicap, age, or 
veteran statuB. 
i 
ABSTRACT 
This study presents numerical simulation of conservative solute transport 
in randomly heterogeneous porous media. Heterogeneity is associated with 
the hydraulic conductivity field, which is described by three parameters: 
mean Iky, variance at, and isotropic correlation length Ay, where Y is the 
natural logarithm of the hydraulic conductivity. Transport simulations are 
performed using a particle tracking random walk (PTRW) method which is 
suitable for treating convection and dispersion processes in a computationally 
efficient manner. The objectives of this study are three-fold. First, the 
correlation structures of the velocity and concentration fields are analyzed as 
a function of the predefined Ay Second, the large-scale, spatial effects of 
the variable flow field on the developing solute plume are analyzed. Third, 
nonergodic effects are investigated. 
The numerical results indicate that in most of the performed simulations 
nonergodic effects occur. Stochastic theory predictions of longitudinal and 
transverse mixing differ from the computational results. These differences 
may be caused by the highly anisotropic velocity correlation structure, which 
exhibits an apparent hole--effect in the transverse direction, 
ii 
TABLE OF CONTENTS 
ACKNOWLEDGEMENTS 
ABSTRACT 
LIST OF FIGURES 
LIST OF TABLES 
1. INTRODUCTION 
2. HETEROGENEITY IN POROUS MEDIA 
2.1 
2.2 
2.3 
2.4 
Introduction 
Hydraulic Conductivity 
2.2.1 Field Findings 
2.2.2 Statistical Representation 
2.2.3 Correlation Structure 
Other Properties 
2.3.1 Porosity 
The Ergodic Hypothesis 
2.4.1 Length Scales in Porous Media Modeling 
2.4.2 Monte Carlo Simulations 
3. GOVERNING EQUATIONS 
3.1 
3.2 
3.3 
Introduction 
Flow Equations 
Transport Equation 
3.3.1 Convection 
3.3.2 Dispersion 
4. NUMERICAL METHODS 
4.1 
4.2 
4.3 
4.4 
Introduction 
Hydraulic Conductivity Field Generation 
4.2.1 Turning Bands Method 
Velocity Field Generation 
4.3.1 Galerkin Finite Element Method 
Concentration Field Generation 
4.4.1 Particle Tracking Random Walk 
iii 
Page No. 
i 
ii 
v 
ix 
1 
4 
4 
4 
5 
5 
10 
12 
12 
12 
13 
14 
15 
15 
15 
16 
16 
16 
18 
18 
18 
18 
19 
20 
20 
21 
~ 
I 
I 
5. 
6. 
7. 
METHODS FOR ANALYZING RESULTS 
5.1 Introduction 
5.2 Correlation Structures 
5.2.1 Hydraulic Conductivity 
5.2.2 Velocity 
5.2.3 Concentration 
5.3 Spatial Moments 
5.3.1 Velocity 
5.3.2 Spreading 
5.3.3 Macrodispersivity 
5.4 Temporal Moments 
5.4.1 Breakthrough Curves 
NUMERICAL SIMULATIONS 
6.1 Introduction 
6.2 Physical Domain 
6.3 
6.2.1 Hydraulic and Transport Parameters 
Correlation Structures 
6.3.1 Hydraulic Conductivity 
6.3.2 Velocity 
6;3.3 Concentration 
6.4 Spatial Moments 
6.4.1 Velocity 
6.4.2 Spreading 
6.4.3 Macrodispersivity 
6.5 Temporal Moments 
6.5.1 Mean Arrival Times 
6.5.2 Standard Deviation of Arrival Times 
CONCLUSION 
REFERENCES 
APPENDIX A 
iv 
Page No. 
24 
24 
24 
24 
25 
26 
27 
28 
29 
30 
31 
32 
33 
33 
33 
35 
38 
38 
41 
49 
52 
52 
54 
65 
65 
70 
75 
80 
83 
89 
LIST OF FIGURES 
Figure No. 
2-1 Measurements of permeability of COreS taken from a borehole 
in the Mt. Simon aquifer in cent,;a1 Illinois [from Bakr, 
1976]. . 
2-2 Frequency histograms for K and inK at the Borden site 
[from Sudicky, 1986]. 
2-3 Method for determining correlation length. 
6-1 Flow diagram for numerical methods used in the numerical 
simulations. 
6-2 Size and shape of flow domain, including injection point, 
boundary conditions, and mean flow direction. 
6-3 Comparison between the longitudinal and transverse 
correlation structures for .inK and the assumed exponential 
correlation function, (Ay=5m). 
6-4 Comparison between the longitudinal and transverse 
correlation structures for .inK and the assumed exponential 
correlation function, (Ay=lOm). 
6-5 Correlation structure of the longitudinal velocity in the 
longitudinal and transverse directions, (Ay=5m, ui=1.0). 
6--6 Correlation structure of the transverse velocity in the 
longitudinal and transverse directions, (Ay=5m, ui=1.0). 
6-7 Correlation structure of the longitudinal velocity in the 
10ngitl1:dinal and transverse directions, (Ay=lOm, ui=1.0). 
6-8 Correlation structure of the transverse velocity in the 
longitudinal and transverse directions, (Ay=10m, ui=1.0). 
v 
I 
..., 
Figure No. 
6-9 Correlation structure of the longitudinal velocity in the 
longitudinal direction, (Ay::::5m, O"Y=0.5,l.O,1.5). 
6-10 Correlation structure of the transverse velocity in the 
longitudinal direction, (Ay =5m, O"Y:=O.5,l.O,1.5). 
6-1i Longitudinal and transverse correlation structures for 
concentration at 180 days, (Ay =5m, O"Y=l.O). 
6-12 Longitudinal correlation structures for Y(x),. Vl(X), and 
c(x,t=180), (Ay =5m, O"y::::l.O). 
6-13 Horizontal trajectory of the plume center of mass, (A y :=5m, 
O"Y:=1.5). 
6-14 Horizontal displacement of plume center of mass as a 
function of time since injection, (Ay =10m, O"Y=0.5,1.0,1.5). 
6-15 Longitudinal displacement variances about plume center, 
(Ay =5m, O"Y=0.5,1.0,1.5). 
6-16 Comparison of numerical and theoretical [Dagan, 1984] 
longitudinal displacement variances, (Ay =5m, O"Y::::::0.5,1.5). 
6-17 Longitudinal displacement variances about plume center, 
(Ay =10m, O"Y=0.5,1.0,1.5). 
6-18 Comparison of numerical and theoretical [Dagan, 1984] 
longitudinal displacement variances, (Ay =10m, O"Y=0.5,1.5). 
6-19 Transverse displacement variances about plume center, 
(Ay =5m, O"Y=O.5,l.O,1.5). 
6-20 Comparison of numerical and theoretical [Dagan, 1984] 
transverse displacement variances, (Ay =5m, O"Y=0.5,1.5). 
vi 
Figure No. 
6-21 Transverse displacement. variances about plume center, 
(Ay=10m, ut==O.5,1.0,1.5). 
6-22 Comparison of numerical and theoretical [Dagan, 1984] 
transverse displacement variances, (Ay=lOm, ut=O.5,1.5). 
6-23 Longitudinal macrodispersivity and theoretical asymptotic 
values [Gelhar and Axness, 1983], (Ay=5m, ut==0.5,1.0,1.5). 
6-24 Longitudinal macrodispersivity and theoretical asymptotic 
values [Gelhar and Axness, 1983], (Ay=lOm, ut=O.5,l.O,1.5). 
6-25 Transverse macrodispersivity and theoretical asymptotic 
values [Ge1har and A:x:ness, 1983], (Ay=5m, ut=0.5,1.0,1.5). 
6-26 Transverse macrodispersivity and theoretical asymptotic 
values [Gelhar and Axness, 1983], (Ay=10m, ut=0.5,1.0,1.5). 
6-27 Mean arrival time through Face 1, (Ay==5m., ui=0.5,1.0,1.5). 
6-28 Mean arrival time through Face 1, (Ay=::10m, 
uf==O.5,1.0,1.5). 
6-29 Mean arrival time through Face 2, (Ay==5m, ut==0.5,1.0,1.5). 
6-30 Mean arrival time through Face 2, (Ay==10m, 
ut=O.5,1.0,1.5). 
6-31 Standard deviation of the mean arrival time through Face 1, 
(Ay ==5m, uf=0.5,1.0,1.5). 
6-32 Standard deviation of the mean arrival time through Face 1, 
(Ay=10m, uf==0.5,1.0,1.5). 
vii 
I 
Figure No. 
6-33 Standard deviation of the mean arrival time through Face 2, 
(Ay=5m, O"f=O.5,1.0,1.5). 
6-34 Standard deviation of the mean arrival time through Face 2, 
(Ay=lOm, O"f=O.5,1.0,1.5). 
A-I Two-dimensional example for calculating K for flow 
simulations. 
viii 
CHAPTER 1. INTRODUCTION 
The emphasis in recent groundwater investigations has shifted away 
from groundwater quantity to groundwater quality, due to an increasing 
assault by chemicals, radioactive waste, and pesticides. Often by the time 
the subsurface contamination is discovered the aquifer water supply is 
damaged beyond repair. Groundwater pollution increases the cost of 
providing clean drinking water. Once a pollutant is discovered it is vitally 
important to be able to predict its movement. This is a difficult problem 
because the dissolved contaminant migrates through a porous material which 
displays a large degree of spatial variability. A description of the subsurface 
material is needed in order to accurately predict the movement of 
contaminants within flow systems. 
Porous media are heterogeneous and nonuniform. Naturally interspersed 
in the ground may be lenses, layers, and fractures composed of materials 
ranging from sands and gravels to clay and rocks. Field samplings reveal 
that soils and geological formations vary in an irregular manner in space. 
Porous materials vary at a scale so small that it is impossible to measure 
aquifer properties accurately over any great distance. Therefore, in the past, 
measurements of hydrogeological variables were made at the small scale and 
these values were ( assumed to represent a much larger area. These 
deterministic values were then inserted into the equations which describe 
groundwater flow and solute transport. . However, this deterministic approach 
did not adequately describe the contaminant plume behavior witnessed in the 
field. Because subsurface material is so variable, a stochastic approach to 
groundwater flow and transport prediction has been investigated in the last 
fifteen years. 
Stochastic modeling assumes that the values of the variable 
hydrogeological parameters cannot be determined at the large-scale. 
Therefore, in this approach, these parameters are assumed to be random 
variables. The value at each location for these variables comes from a 
prespecified (estimated from available measurements) probability density 
function. A number of researchers have developed stochastic theories for 
describing subsurface flow and transport phenomena. 
Freeze [1975] performed computer simulations with randomly variable 
hydrogeological parameters and concluded " ... the most realistic representation 
of a naturally occurring porous medium is a stochastic set of macroscopic 
elements ... " These macroscopic elements are regions which contain uniform, 
but random, values for the parameters. Although this study lacked some 
realistic features (i.e., spatial correlation) it was the first to attack 
groundwater modeling with stochastic methods. 
1 
Matheron and de Marsily [1980] demonstrated for the special case of a 
stratified porous medium that the rate of spreading of a contaminant plume 
may not be Fi cki an. , Because of this non-Fickian behavior the classical 
convection-dispersion transport equation may not be appropriate. However, 
field studies have found that, in some heterogeneous porous materials, plume 
spreading is non-Fickian at early stages of transport but that after a plume 
has traveled a sufficient distance the spreading rate does become constant 
(obeys Fickian law). . 
Gelhar and Axness [1983] developed stochastic solutions for the general 
case of a three-dimensional, anisotropic heterogeneous porous medium using 
small perturbation techniques and the ergodicity assumption. Their results 
predict asymptotic plume spreading rates. The small perturbation technique 
assumes that the variability of a random variable is small so that a 
first-order approximation of the variable is adequate. The ergodic hypothesis 
states that a single heterogeneous aquifer appro:ximates the entire ensemble of 
aquifers with the same statistical characteristics. This allows the acquired 
. resUlts to be generally valid and not a function of a particular porous 
medium structure. 
Dagan [1984] found expressions which relate plume spreading to the 
spatial characteristics of the porous medium and to travel time. This theory 
also assumes small perturbation and ergodicity. Because travel time is 
included) the theory is more general than that of Gelhar and Axness [1983]. 
In order to confirm these theories, field-scale experiments and numericaI 
simulations have been performed. 
Recent field-scale experiments have provided the physical evidence 
needed to test stochastic theories. A natural gradient tracer test at the 
Borden aquifer in Ontario, Canada has provided valuable information and has 
been the subject of many studies, including those by Freyberg [1986], Sudicky 
[1986], Mackay et al. [1986], and Barry et al. [1988]. A tracer test at the 
Twin Lake aquifer was conducted and the data was used to investigate 
fie1d-scale spreading rMoltyaner and Killey, 1988a, b]. Field tests may take 
years to complete and few have been conducted. 
Numerous researchers have performed numerical simulations in recent 
years using the stochastic approach to the groundwater flow and transport 
problem. Smith and Schwartz [1980, 1981] performed two-dimensional 
numerical simulations using the Monte Carlo approach to generate hydraulic 
conductivity fields. They found that transport results were uncertain even 
when the statistical structure of the porous media is known. For example, 
the time it takes a plume to completely pass an area is highly dependent on 
the type of heterogeneity assumed for the porous media. The more 
heterogeneous the material the more uncertain the predictions are. 
2 
, 
1 
1 
Salandin and Rinaldo [1990] performed Monte Carlo numerical 
simulations which were similar to that of Smith and Schwarb [1980, 1981]. 
However, they preserved the prespecified spatial correlation structure of the 
hydraulic conductivity fields. Their results from transport simulations seem 
to indicate that the first..-order approximation of the small perturbation 
assumption has a robust solution for large degree of variability. 
Rubin f1990), following the work of Graham and McLaughlin [1989] and 
Dagan [1984, analytically derived the pore-water velocity covariance tensor 
and expressed it in terms of the parameters of the hydraulic conductivity 
field. He showed that the velocity correlation is anisotropic when the porous 
medium is assumed to be isotropic. He applied the covariance derivations to 
field data from the Borden test site and found agreeable results. 
Thompson and Gelhar [1990] examined large...-scale dispersive plume 
behavior in a variable three-dimensional flow field. They performed 
numerical simulations using single realizations, and three different degrees of 
heterogeneity. The numerical simulations indicated that stochastic theory 
accurately predicted the experimental longitudinal movement and mixing, but 
differed markedly from the experimental transverse mixing. The simulations 
also demonstrated that significant nonergodic effects may occur in individual 
experiments. In the more variable flow fields, asymptotic plume behavior 
often did not develop within the flow domain. They were uncertain whether 
or not this was due to incomplete plume development or inaccuracies in 
stochastic theories for extremely variable porous media. 
In this study, numerical simulations are performed in order to address 
the following issues. First, to investigate what effects the degree of 
heterogeneity in hydraulic conductivity has upon the transport process. 
Second, to investigate what effects the spatial correlation of hydraulic 
conductivity has on the velocity field and the transport process. Third, to 
determine to what extent the spatial correlation of the hydraulic conductivity 
propagates through to the velocity and concentration spatial correlation 
structures. 
3 
1 ___________ _ 
----------------- -------- -----
CHAPTER 2. HETEROGENEITY IN POROUS MEDIA 
2.1 Introduction 
Subsurface porous materials display a large degree of variability. 
Correspondingly, the hydrogeological parameters used to describe a porous 
medium, such as hydraulic conductivity and porosity, are variable and 
heterogeneous. One problem with trying to predict the movement of solutes 
in the ground is how to describe and quantify this heterogeneity. Perhaps 
the most important parameter to quantify is the hydraulic conductivity due 
to the extreme spatial variability and apparent impact on the subsurface flow. 
2.2 Hydraulic Conductivity 
The hydraulic conductivity may be defined as a "combined property of 
the porous medium and the fluid flowing through it" and gives an indication 
of the "ability of a1· uifer material to conduct. water through. it under 
hydraulic gradients II Bear, 1972]. The hydraulic conductivity, K, is 
dependent on the geometry of the porous medium through the permeability, 
K" and dependent on the fluid through the kinematic viscosity, v, following 
K=~ 
v 
(2-1) 
where g is the acceleration of gravity. The dimension of K, is length2, of v is 
length2·time-1, and g is length·time-2, yielding the dimension of K, 
length.time-1• In subsurface flow and transport problems the hydraulic 
conductivity is usually considered the important aquifer parameter instead of 
permeability because K and K, are interchangeable. 
In natural formations the hydraulic conductivity is anisotropic and 
therefore exhibits a tensorial property 
Kll K12 K13 
K = K21 K22 K23 
KS1 KS2 KS3 
(2-2) 
If the Cartesian coordinate axes are aligned with the principal axes of 
hydraulic conductivity, K is reduced to three diagonal components, Kl, K2, 
and Ks. For the sake of simplicity, and without loss of generality, this 
alignment is assumed for the remainder of this work. 
4 
1 
., 
2.2.1 Field Findings 
Numerous investigators have studied the spatial variability of hydraulic 
conductivity. Some. in.vesti~ato. rs, inc.luding .Freeze .. [1975], Delhomme. J1979], 
Hoeksema and Kitanidis l1985], and Gelhar [1986], have collecte and 
analyzed data from a variety of field studies. Other investigators have been 
involved in. large-scale tracer tests which coincided with measurements of 
hydraulic c ..on.d. uctivity, includin.g Sudick.y [1.98.6] .. and Freyber.g. [;1.986] at the 
Borden aquifer in Ontario, Canada, and Moltyaner and Killey [1988a, bl at 
the Twin Lake site, also in Ontario. A common characteristic gathered from 
these studies is that the hydraulic conductivity is highly heterogeneous. For 
example, K can vary three orders of magnitude within a few horizontal or 
vertical meters [Gelhar, 1986; Sudicky, 1986]. Figure 2-1 illustrates the 
variability of permeability within a relatively homogeneous sandstone aquifer. 
The range of measured conductivities is immense, from meters per day in 
coarse sands to fractions of millimeters per day in clay [de Marsily, 1986]. 
Table 2-1 gives an example of potential values of K. In sedimentary media, 
with more or less horizontal stratification, two distinct hydraulic 
conductivities are present, a large horizontal conductivity and a smaller 
vertical conductivity, and the ratio between the two is generally between 1 
and 100 [de Marsily, 1986]. These field findings indicate that the hydraulic 
conductivity is spatially-correlated and random. 
2.2.2 Statistical Representation 
Since natural formations do not possess uniform properties. over large 
areas and sufficient data are usually not available to completely describe the 
small-scale variability, many researchers represent the hydraulic conductivity 
as a stochastic random field. Several researchers have developed large-scale 
transport theories based on the stochastic approach, including Freeze [1975], 
Gelhar and Axness [1983], Dagan [1984], and Gelhar [1986]. Others have 
conducted numerical experiments based on a stochastic representation of the 
hydraulic conductivity including Smith and Schwartz [1980, 19811, Tompson et 
al. [1988], Graham and McLaughlin [1989], Salandin and Rinaldo [1990], and 
Tompson an,d Gelhar [1990]. 
Let x be the Cartesian coordinate vector of a point with Xt,X2 or 
Xt,X2,Xa its components in the two- or thre~imensional space. Let f(x) 
denote a probability density function (pdf) at that point. To completely 
characterize a random field the joint probability density function is required, 
f{xt,x2, ... ,xn}, where Xl,X2,,,.,Xn are the coordinates of points 1,2, ... ,n. If the 
joint pdf is invariant with respect to a translation through space, then the 
stochastic random field is stationary, or statistically homogeneous. If 
statistical homogeneity is restricted to the first two moments, the random 
field is weakly-stationary. With this assumption the mean of the random 
field is constant through space and the covariance depends only on the 
separation distance between two points, not on their positions. 
5 
i 
I t 
0 
,...; 
>-
f-
;;::!o 
OJ· 
<:(N 
W 
~ 
0: 
W P-q 
<..!l 
0 
~ 
0 
0 
r) 
I 
0 
Fig. 2-1 
.•. _. ____ +___ , ___ _ '+- ___ •.• M_ ..... _-+-~~ 
100 200 
DIST A"ICE • FEET 
Measurements of permeability of cores taken from 
a borehole in the Mt. Simon aquifer in central 
Illinois [from Bakr, 1976]. 
6 
TABLE 2-1 
Medium K (em/s) 
Gravels 10° 
Sands 10-3 
Silts 10-7 
Clay 10-10 
7 
-, 
-i 
Many investigators have concentrated on finding the pdf which best 
describes the spatial distribution of the hydraulic conductivity field, K(x). 
Freeze [19751 summarized many years worth of field measurements and a 
large body ot laboratory work and concluded that the log-normal distribution 
seems. to fit hydraulic conductivity field data. Since then, numerous other 
investigators have supported this conclusion. Hoeksema and Kit ani dis [1985] 
analyzed hydrogeological properties from 31 regional aquifers and found the 
natural logarithm of the hydraulic conductivity measurements generally passed 
normality tests. Figure 2-2 indicates that the logarithm of hydraulic 
conductivity data more closely resembles a Gaussian distribution than does 
the original data. At present the assumption of log-normally distributed 
hydraulic conductivity is accepted as a general tenet in natural 
hydrogeological formations [Dagan, 1989]. 
If the distribution of K(:x) is assumed to be log-normal, then the 
variable 
Y(x) ~ lnK(x) (2-3) 
is distributed normally, with mean /ly and variance cry 
(2--4) 
The parameters /ly and cry are defined as 
/ly ~ E(Y(x)) = E(lnK(x)) (2-5) 
cry = E(Y(x) - /ly)2 ~ E[lnK(x) - E(lnK(x))J2 (2-D) 
where E( ) denotes expectation. 
degree of variability of K(x). 
conductivity field is defined as 
The variance of lnK(x) represents the 
The geometric mean of the hydraulic 
(2-7) 
Typical values of /ly range from -5 to -12 for gravels and sands to -21 to 
-28 for silts and clays (units of K are cmjsec) [de Marsily, 1986]. The 
typical values of Kg for these material types range from 100 to 1O~3 and 1O~7 
to 1O~1o (cm/sec). The degree of variability, cry, may range from near zero 
for homogeneous deposits, to five or even higher for extremely variable porous 
media [Dagan, 1989]. The identification of the correlation structure of lnK(x) 
is the final step towards a complete statistical deSCription of the hydraulic 
conductivity field. 
8 
'--~.---.~-.-- -----.--~--..•. -- ... -------.-----~~---.~-.. ~--~~-.------- --------.---.-~--.. --.-... ~----
28.0 
>-- 15,0 
o 
c 
(l) 
:J 10,0 
CT 
(l) 
h-
L!..- 5.0 
o,O~L---~=c=L~~--~~~-L~~~==L-J 
-7.0 -5.0 -5.0 
Y;::ln(K) 
-4.0 -3.0 
Q.o~----------~------------~----------~ 
20.0 
>- 15.0 
o 
c 
(l) 
:J 10.0 
tT 
(l) 
'-
LL 5.0 
O.OL-~--~~~~--~~~~--~~~~--~~ 
0.001 
Fig. 2-2 
0.011 0.021 . 0.031 
HydrBul ic Conductivity (cm/s) 
Frequency histograms for K and lnK at the 
Borden site [from. Sudicky, 1986]. 
9 
2.2.9 Correlation Structure 
Besides using field measurements to describe the randomness of the 
hydraulic conductivity, some investigators use these measurements to find a 
correlation structure for K(x). However, this work is difficult because a large 
number of measurements is needed, particularly at small distances. Hoeksema 
and Kitanidis [1985] assumed an exponential correlation function and found 
that there was good corresponden~e between it and the measurements of 
inK(x). At the Borden site a.n analysis o. f the spatial correlatio. n structure of 
inK (x) revealed that the correlation function could be approximated as 
exponential [Sudicky, 1986J. Because no conclusive agreement exists on the 
most appropriate correlatlOn structure of the hydraulic conductivity, and 
because it is easy to use, the exponential function is often employed in 
subsurface flow and transport modeling. 
If weak stationarity is assumed for the hydraulic conductivity field then 
the exponential correlation structure for lnK(x) is a function of the separation 
distance between two points, Xl and X2, not their positions. Let r be the 
separation distance between Xl and X2 
where r i is the separation distance in the i-th prinCipal directions. The 
spatial correlation of the hydraulic conductivity field is then described by 
py(r) = exp(-[(rl/AY1)2 + (r2/AY2)2 + (ra/AYa)2Ph) (2-9) 
where py(r) is the· correlation between two points separated by a distance of 
rand Ay i is the correlation length of Vex) in the i-th principal direction. 
The correlation length is used to parameterize the spatial correlation and is 
an approximate measure of the distance at which correlation ceases (the 
correlation becomes zero as the distance between two points approaches 
infinity). A larger correlation length implies greater spatial continuity in 
hydraulic conductivity. In this study, the correlation length is defined to be 
equal to the separation distance between two points at which the correlation 
becomes equal to e-i . Figure 2-3 is an example of how the correlation length 
is found. . 
In some natural formations it has been found that the hydraulic 
conductivity correlation structure in the plane parallel to layering is 
approximately isotropic and that the correlation length in the horizontal plane 
is much larger than the correlation length normal to layering rGelhar, 1986]. 
For example, at the Borden site the horizontal correlation length is 2.8m and 
the vertical correlation length is O.lm [Sudicky, 1986]. If the vertical 
correlation length of an isotropic hydraulic conductivity field is assumed to be 
much smaller than the horizontal correlation length then the exponential 
correlation function reduces to 
10 
1.0 
0.9 
0.8 
0.7 
-~ 
---Q.. 0.6 
.. 
z 
0 
- 0.5 f-o 
:s 
rzJ 
0:: 0.4 0:: 
0 e- t 
..... u 
..... 0.3 
0.2 
0.1 
0.0 I 
0 1 2 A. 3 4 5 6 7 
SEPARATION DISTANCE, r 
Fig. 2-3 Method for determining correlation length. 
py{r) :;: exp[-r/ Ay] (2-10) 
where rand Ay are the horizontal separation distance and correlation length, 
respectively. 
Under the assumptions of an exponential covariance function and 
log-normal distribution, the random isotropic hydraulic conductivity field is 
statistically characterized by the three parameters: /Jy, O"f' and Ay. If these 
parameters are known, or if they are estimated from available data, then 
stochastic studies of groundwater flow and transport can be' performed. 
2.3 Other Properties 
Other hydrogeological properties, such as porosity, are also spatially 
variable in porous media. Field measurements have been made for these 
properties as well as for hydraulic conductivity. However, these studies 
indicate that the effects of their variability upon solute transport are 
substantially smaller than those of the hydraulic conductivity. 
2.3.1 Porosity 
Porosity, n, which is dimensionless, may be defined as the ratio between 
the volume of water able to circulate in· a porous medium to the total 
volume of the porous medium [de MarsHy, 1986]. Freeze [1975] suggested 
that the porosity can be approximated with a normal distribution. Gelhar et 
al. [1979] showed analytically that the porosity variability has only a minor 
effect in large-scale transport problems. Porosity measurements at the 
Borden site reveal a mean value of 0.33 and a standard deviation of 0,017, 
which, relative to the variability of the hydraulic conductivity, is very small. 
Dagan [1989] concluded that the variability of the porosity may be four 
orders of magnitude smaller than the variability of .fuK. These results seem 
to indicate that in the study of solute transport in heterogeneous porous 
media, porosity is of little importance compared to the hydraulic conductivity. 
2.4 The Ergodic Hypothesis 
Stochastic studies of flow and transport phenomena should satisfy the 
ergodic hypothesis. The hypothesis requires a random field (for example the 
hydraulic conductivity random field) to include all possible states; or 
variabilities, of the physical property it represents. This requirement insures 
that any conclusions reached from numerical simulations are general and not 
a function of a particular random field. Theoretically, ergodicity is satisfied 
only when a domain size approaches infinity or when an infinite number of 
realizations are produced which have finite domain size. Researchers have 
recommended some guidelines which should be followed in order to satisfy 
12 
ergodicity in numerical simulations. These guidelines include the careful 
selection of certain length scales, and/or performing multiple, or Monte Carlo, 
simulations. These guidelines are presented next. 
2.4.1 Length Scales in Porous Media Modeling 
Three length scales are often used in numerical models of large-scale 
porous formations: a representative elementary length (AX), a correlation 
length (A), and the flow domain size (L). The flow domain is' discretized 
into small zones, or representative elementary volumes (REV), which have 
length Ax per side. Within each REV a non-varying v81ue for every 
hydrogeological parameter is assigned. The second length scale is Ay. The 
third length scale is the length of the domain, L. Two dimensionless ratios, 
Ay/ Ax and L/ Ay, are important if the ergodic hypothesis is to be satisfied. 
The ratio Ay/ Ax must be greater than one so that the statistical 
structure of the correlation is preserved. When this ratio is increased, the 
correlation structure of the hy:draulic conductivity is more accurately 
discretized. Smith and Schwartz [1980] worked with ratios between 1 and 2, 
while Tompson et al. [1988] conducted three-dimensional numerical transport 
simulations with Ay/l:::.x=2. Ababou et al. [1988] investigated ergodic 
requirements and suggested a conservative lower limit of 4 to 5. Graham 
and McLaughlin [1989] used ratios of 4 and 8. In this study the limits 
suggested by AbaDou et al. [1988] will be followed. 
The second ratio, L/ Ay, must be quite large to insure that enough 
uncorrelated samples appear in the flow domain. Furthermore, this ratio 
must be large if large-scale transport and asymptotic plume behavior are of 
interest. Ababou et al. [1988] suggested 40 to 50. Salandin and Rinaldo 
[1990] performed two-dimensional transport simulations with L/ Ay=18 in one 
direction and 36 in the other. Graham and McLaughlin [1989] used ratios 
ranging from 11 to 44 in one direction and 6 to 24 in the other. Tompson 
et al. [1988] produced three-dimensional hydraulic conductivity fields with 
L/ Ay as low as 12.5. In this study the ratio will be much higher than one 
and the relation 
Ax < Ay « L (2-11) 
will be satisfied. 
13 
2.4.2 Monte Oarlo Simulations 
Along with the careful selection of the three . length scales, the Monte 
Carlo method can help insure that ergodicity is satisfied. This technique is 
powerful, conceptually simple, and, with high speed computers, increasingly 
accessible. The technique is also the proper way for verifying ensemble 
stochastic theories [Tompson et al., 1988.;1. The Mo.nte Carlo .. method consists 
of performing repeated simulations of the flow and transport problem with 
different realizations of the spatially-correlated random hydraulic conductivity 
field. Each realization has. the same input characteristics (1Jy! uy! and Ay), 
but are otherwise entirely independent. 
Several researchers have utilized the Monte Carlo method in stochastic 
models of subsurface flow and transport. Freeze [19751 analyzed 
one-dimensional groundwater flow using random, but uncorrelated, hydraulic 
conductivity fields. Smith and Schwartz [1980] investigated two-dimensional 
mass transport behavior with random, spatially-corre1ated hydraulic 
conductivity fields. Because of the development of faster computers, 
researchers are now able to apply the Monte Carlo method to larger domains, 
in order to investigate large-scale . effects caused by a random, 
spatially-correlated hydraulic conductivity field. 
14 
CHAPTER 3. GOVERNING EQUATIONS 
3.1 Introduction 
As was previously described, the hydraulic conductivity is the most 
important hydrogeological parameter. In this study it is represented as an 
isotropic, weakly.:..stationary, spatially.-correlated random field. In this chapter 
the fundamental equations which describe subsurface flow and solute 
movements are introduced. The importance of hydraulic conductivity in these 
equations is also analyzed. 
3.2 Flow Equations 
In investigations of large heterogeneous flow systems the conservation of 
mass equation and DarcY's law describe water movement. Darcy's law relates 
the flux of water through an area to the hydraulic conductivity. If the 
Cartesian coordinates axes are aligned with the principal axes of the hydraulic 
conductivity Darcy's law is 
q(x) =: - K(x) V¢>(x) (3-1) 
where q(x) is the specific discharge (dimension of length·time~l) and ¢>(x) is 
the hydraUlic head tdimension of length). By combining. the conservation of 
mass equation 
V.q(x) =: 0 (3-2) 
with (3-1), the steady-state flow equation, with constant porosity, is 
V· [K(x) V¢>(x)] =: 0 (3-3) 
Given the random distribution of the hydraulic conductivity and suitable 
boundary conditions, (3.,..3) is solved numerically for the hydraulic head. 
With the hydraulic head known, the velocity is determined from Darcy's law 
v(x) =: - -} [K(x) V¢>(x)] (3-4) 
where v(x) is the. pore-water velocity. 
15 
I 
-; 
3.3 Transport Eguation 
The general equation describing transport of a nonreactive, conservative 
solute in a three-dimensional saturated porous medium, with constant 
porosity, is 
8c(x,t) + V. [c(x,t) vex)] - V· [D(x) Vc(x,t)] ::::: 0 (3-5) 
at 
where c(x,t) is the concentration (dimension of mass·length~3) and D(x) is 
dispersion coefficient tensor (dimension of length2• time~l). The components of 
D(x) are 
D(x) ::::: aT ! v(x)! I + (a - a ) v(x) ·v(x) + Dm (3--6) 
L T Iv(x)! 
where aL and aT are the longitudinal and transverse local dispersivities 
(dimension of length), ! vex) I is the magnitude of the velocity, I is the 
identity matrix, and Dm is the coefficient of molecular diffusion (dimension of 
length2.time~1). In very slow moving groundwater systems molecular diffusion 
contributes very little to the dispersion process and is not considered. in this 
study. The transport equation (3-5) is used to solve for the spatial and 
temporal distribution of concentration once the components of D(x) are known 
and the spatial distribution of velocity is found from the methods described 
in Section 3.2. The transport equation includes two mechanisms of solute 
transport: convection and dispersion. 
3.3.1 Convection 
Convection is the process in which a dissolved substance is carried along 
with the fluid velocity. In (3-5) the second set of terms on the left-side 
represent convection. Convective movement is directly proportional to the 
velocity, which in turn. is defend. ent on the distribution of the hydraulic 
conductivity field through (3-3 and (3-4). Although K(x) is not explicitly 
included as a variable, it indirectly influences convection through vex). 
3.3.2 Dispersion 
Dispersion is a mDong phenomenon caused by pore-water velOcity 
variabilities within the porous medium. It is the physical process which 
describes the tendency of a plume to spread out from the path that would be 
expected from convective movement alone and is represented by the third set 
of terms on the left side of (3-5). 
16 
Field investigations show in a consistent manner that the spreading 
experienced in the field is far more than that found in the laboratory. For 
example, Gelhar [1986] compiled measurements from around the world and 
found that, in general, the spreading rates witnessed in the field were orders 
of magnitude larger than those determined in the laboratory. These findings 
indicate that dispersion is scale-dependent and that the convection-dispersion 
equation can not adequately describe natural transport behavior from a 
deterministic approach. Just in the last fifteen years have these problems 
been addressed from a stochastic approach. 
In the stochastic modeling approach dispersion consists of two 
components: local dispersion and macrodispersion. Local, or pore-scale, 
dispersion is caused by velocity fluctuations through the pores of the medium. 
In terms of this study the scale of these fluctuations are much smaller than 
the REV. The local dispersivities, aL and aT' parameterize this local 
dispersion. Macrodispersion is a result of velocity fluctuations caused by 
large-scale heterogeneity of the hydraulic conductivity. These fluctuations are 
on the same scale as Ay . The parameter which describes the macro dispersion 
is the macrodispersivity, A(t). Gelhar and Axness [1983] obtained closed 
form expressions which relate the asymptotic macrodispersivity; A( t-too) , to o-f 
and Ay . Their analysis included two- and three-dimensions, anisotropic 
hydraulic conductivity conditions, and various degrees of variability of lnK(x). 
They compared their theoretical formulations to field tracer tests and found 
that asymptoticmacrodispersivity may only be reached after a plume has 
traveled dozens of correlation lengths of lnK(x). Dagan [1984] developed 
equations which relate the spread of a plume to. o-f' Ay, and the travel 
distance of the plume. His equations are more general than those developed 
by Gelhar and Axness [1983], however the local dispersion is ignored. 
17 
CHAPTER 4. NUMERICAL METHODS 
4.1 Introduction 
A wide variety of physical quantities in hydrologic research have been 
represented by artificial random fields. These synthetic fields are the driving 
force in large-scale simulations of heterogeneous phenomena. Since the 
hydraulic conductivity is assumed to be a spatially--correlated random field, a 
synthetic generation of such a field is needed. A random hydraulic 
conductivity field is used as input for the flow equation and Darcy's law to 
determine the spatial distribution of the velocity field. The velocity field is 
used as input for the convection-dispersion equation to solve for the spatial 
and temporal distribution of concentration. 
4.2 Hydraulic Conductivity Field Generation 
Several numerical techniques are available to produce replicates of the 
spatially--correlated random hydraulic conductivity field, including the Nearest 
Neighbor method [Bartlett, 1975]~ the standard. Four.ier integrat.ion m. ethod 
[Meija and Rodriguez-Iturbe, 1974J, and the turning bands method [Matheron, 
1973J. However, prohibitive computational costs or storage requirements may 
result when the number of points, or nodes N, in the random field are large. 
For example, the Nearest Neighbor method, as applied by Smith and Freeze 
[1979], requires inversion of an N x N banded matrix and does not guarantee 
stationarity. The standard Fourier integration method requires small matrix 
storage but high computational costs for N greater than 10,000 [Tompson et 
al. , 1989]. The turning bands method has been found to be accurate and 
economic for producing large random fields and is employed in this study for 
the generation of the hydraulic conductivity fields. 
4.2.1 Turning Bands Method 
The turning bands method was originally developed by Matheron [1973] 
and its practical use for developing a spatially--correlated random hydraulic 
conductivity field was advanced by Mantoglou and Wilson [1982] and 
Tompson et al. [1987, 1989]. A version developed by Tompson [1990] is used 
in this study. A detailed description of the turning bands technique is 
beyond the scope of this study, but can be found in Tompson et al. [1987, 
1989]. 
The turnin~ bands method generat. es an isotropic, spatia.lly--correlated 
random field, Zlx), from a standard normal distribution and a prespecified 
correlation structure 
18 
I 
1---- .- .... 
Z(x) : N(O,1) (4-1) 
If the correlation structure of the hydraulic conductivity field is assumed to 
be represented by the exponential function then the spatial correlation of Z(x) 
can be parameterized by Ay. Point values of K(x) are determined from 
point values of Z(x) by first calculating a new variable Y(x) 
Y(x):= /1y + Uy Z(x) (4-2) 
where /1y and uy are the mean and variance of inK (Section 2.2.2). Y(x) is 
distributed norm<!lly and has correlation length Ay. Point values of K(x) are 
found by 
K(x) := exp[ Y(x)] = exp[ /1y + Uy Z(x)] (4-3) 
Given a correlation length of inK, the effects on flow or transport 
caused by different uy can be investigated in the following way: 1) A 
realization of Z(x) is created, 2) Different realizations of K(x) are produced 
from one realization of Z(x) by altering u{- in (4-a), and a) The flow and-
transport equations are solved for the velocity field v(x) and the 
concentration field c(x,t) given each distribution of K(x). Likewise, the 
effects on flow and transport caused by different correlation lengths can be 
investigated by producing realizations of· Z(x) with a new Ay. 
4.a Velocity Field Generation 
Given the distribution of the hydraulic conductivity field from (4-3), 
various numerical methods can be used to solve the general flow equation 
(a-3) for the hydraulic head field, ¢(x). Some of the more frequently 
empfoyed methods include finite difference, finite element, boundary element, 
and analytical element. In this study, the Galerkin weighted residuals 
approach to the finite element method is employed [Hromadka II et al., 
19851. A brief description of the method is given; a more detailed description 
of the Galerkin finite element method can be found elsewhere [Pinder and 
Gray, 1977; Huyakorn and Pinder, 1983; Bear and Verruijt, 1987]. 
19 
4.3.1 Galerkin Finite Element Method 
In the Galerkin finite element method the physical system is discretized 
into triangular subregions) or elements. Values of the random hydraulic 
conductivity field are assigned uniformly within each element so that K(x) is 
discontinuous between elements. Because the turning bands method assigns 
its output to evenly spaced nodes) and this finite element method utilizes 
triangles, the hydraulic conductivity values for each triangular element need 
to be averaged, Appendix A describes how this averaging was achieved. 
Once the boundary and initial conditions are selected the Galerkin finite 
element method solves the flow equation (3-3). Since the initial head field is 
only an approximate guess a non.....zero residual results on the right-side of 
(3-3). A steady-state solution for the hydraulic head distribution is reached 
when the average residual for the entire system is acceptably smalL 
A velocity field v(x) is determined from the discrete nodal distributions 
of ¢(x) and K(x)) and Darcy's law (3-4). Some interpolation scheme must 
be adopted. In this study, cell velocities are calculated using linear 
interpol.ation, so that each. cell contains a constant velocit:y. Though this is 
not the most accurate interpolation scheme [Goode, 1990j, it is a suitable 
approximation that offers considerable computational savings, because 
evaluation of each cell velocity would only need to be done once. Tompson 
et aL [1988] recommended an interpolation approach that computes cell 
velocities in terms of mid-node conductivities. Their approach is adopted for 
this study. 
4.4 Concentration Field Generation 
The purpose of a numerical transport model is to accurately solve the 
convection-dispersion equation (3-5) for the spatial and temporal distribution 
of concentration c(x,t). The numerical methods often used include the finite 
element, finite difference, and the particle tracking random walk. 
Numerical dispersion is an inherent problem in the finite element and 
finite difference methods. A measure of the numerical dispersion is the grid 
Peclet number, which is conservatively approximated by 
(4-4) 
if molecular diffusion is ignored. When PeG is greater than ten numerical 
dispersion appears [Huyakorn and Pinder, 1983]. Kinzelbach [1988] gave an 
example of potential inaccuracies when aLI aT»lO. In order to circumvent 
these numerical problems the size of the REV would need to be reduced, 
which would increase the computer cost and storage, to impractical limits. 
20 
The particle tracking random walk (PTRW) method is free from 
numerical dispersion [Kinzelbach, 1988]. Another advantage of the PTRW 
method is that concentration distributions only need to be calculated when it 
is of interest, instead of at every time step. However, this method does not 
predict point concentration values as accurately as the other methods. 
Numerous versions of the PTRW have been developed which solve (3-5), 
including those by Ahlstrom et al. [1977], Prickett et al. [1981], Tompson et 
al. [1988], and Kinzelbach [1988]. Many studies have used the PTRW 
method in stochastic models, including Tompson et al. [1988], Kinzelbach 
[19881, Rubin.. [1990], Salandin a.nd Rin.aldo [19.901, and. Tom. ps.on an. d Gelhar 1990]. Because this study is interested in overall mean plume behavior, and 
not the accurate prediction of point concentrations, the PTRW method has 
been chosen as the transport model. The numerical scheme used is based on 
the techniques described by Prickett et al. [1981]' Kinzelbach [1988] and 
Tompson et al. [1988]. A detailed description of the particle tracking random 
walk method is beyond the scope of this study, but can be found in 
Kinzelbach [1988] and Tompson et al. [1988]. 
4.4.1 Particle Tracking Random Walk 
The particle tracking random walk method represents an injected mass 
as a large collection of particles, each of which is assigned an identical 
unchanging mass value. In each time step every particle is moved by a 
convective displacement and a dispersive displacement. The convective step 
is deterministic and moves particles in the direction of the pore-water 
velocity at that point. The dispersive step is random and is based on the 
physical process of dispersion. At the end of each time step the position of 
each particle is known so that the average plume statistics may be calculated. 
These statistics are collected at certain time intervals so that the time rate 
of change of these statistics may be analyzed. Once a particle exits the 
domain it is no longer transported and the total mass in the domain is 
reduced proportionally. 
A brief mathematical description of the PTRW method follows. The 
velocity field distribution v(x) is defined for every point in the domain and is 
constant within each cell, or REV. Each cell also contains a constant value 
of the dispersion coefficient, D(x), which is determined from (3-6). Let xP(t) 
represent the position of the p-th particle at time t. In each time step .6.t 
the particle is transported a convective displacement, Ll~, and a dispersive 
displacement, LlX~ 
21 
-< 
I 
i 
The deterministic convective displacement is a function of the velocity at 
:xP(t) , v(xP(t)), and the gradient of the dispersion tensor at xP(t) 
Ax~ -= [ v(xP(t)) + V.D(:xP(t))] . At (4..-6) 
The gradient terms in (4..-6) are important near stagnation points. If these 
terms are neglected an unphysical buildup of particles may occur in regions of 
low hydraulic conductivity [Kinzelbach, 1988], The random dispersive 
displacement in each time step is, in two dimensions, 
Axa -= [ 2 aL V1(XP(t)) ·At p/2 • Z1 + 
[ 2 aT V2(:xP(t)) • At p/2 ,. Z2 (4-7) 
where Z1 and Z2 are two random values from a standard normal distribution. 
The selection of At is an important choice in the PTRW simulation 
process. The computational expense of this method is proportional to the 
number of time steps. If a very small time step is used, the computational 
expense becomes prohibitive. If too large a time step is used, overshoot 
errors in (4-5) may occur [Tompson and Gelhar, 1990]. The cell Courant 
Number is the ratio between the average convective aisplacement and the 
grid spacing 
c -= v(xloAt 
c x 
where v ( x) is the mean pore-water velocity. Use of 
force Cc greater that one result in overshoot problems. 
choice of At will be made such that C c < 1. 
(4-8) 
times steps which 
In this study the 
In the PTRW method the total solute mass is represented by a 
hypothetical group of particles. By increasing the number of particles the 
solution to (3-5) becomes more consistent and reliable, but the computational 
expenses will increase proportionally. Though point concentration predictions 
become much more accurate with a larger number of particles, the accuracy 
of the overall plume behavior does not increase to the same degree. Since 
average plume behavior is of interest in this study a minimal number of 
particles will be used. The use of approximately 1000 particles has been 
found to be a substantial number for an accurate representation of overall 
plume behavior [Tompson et al., 1988]. 
22 
, 
..... 
I 
I 
--1 
As with flow simulations, transport simulations require initial and 
boundary conditions. Initial conditions involve describing the distribution of 
particles at the beginning of the simulation, xP(t:::;O). These particles may be 
initially distributed at a point, on a line, or within an area. The injection 
may be either instantaneous or continuous over time. 
The boundary conditions involve describing the cells on the outermost 
perimeter of the domain. Two distinct t;ypes of boundaries ma1 be identified, 
the free flow and the no flow boundary <rAhlstrom et al., 1977J. A free flow 
boundary allows particles to leave the domain, thus eliminatmg them from 
further simulation. When this occurs the total mass within the domain is 
decreased proportionally and mass in the system is no longer conserved. A 
no flow boundary does not allow a particle to pass through it. Such a 
boundary will reflect particles back· into the system, like a billiard ball 
[Tompson et al., 1988]. 
23 
CHAPTER 5. METHODS FOR ANALYZING RESULTS 
5.1 Introduction 
Two types of methods are used to analyze the results of the numerical 
simulations performed. The first involves the correlation structure of the 
random fields, Y(x), v(x), and c(x,t). The hydraulic conductivity correlation 
structure is assumed to be isotropic, weakly--stationary, and represented by 
the exponential function. From the flow and transport simulations the 
velocity and concentration fields are produced. The correlation structure of 
each of these fields may be an~yzed. The second analytical method used is 
moment analysis. This moment analysis includes spatial and temporal 
moments. 
5.2 Correlation Structure 
The random fields generated by the turning bands technique are 
analyzed to insure that they are accurate replicates of Y(x) :;= inK(x). The 
correlation structure of each random field is calculated to insure that it 
approximates the correlation structure assumed. Likewise the correlation 
structures of the velocity and concentration fields are analyzed. The three 
correlation structures are then compared. The methods for determining the 
correlation structures of Y(x), v(x) , and c(x,t) are now given. 
5.2.1 Hydraulic Conductivity 
The correlation structure of the hydraulic conductivity field is 
determined by calculating. the correlation between values of Y(x). Let 
Y(xn)m represent the value of Y{x) at point xn' n=1,2, ... ,N, in the m-th 
Monte Carlo realization, m:::::l,2, ... ,M, where N is the number of nodes in the 
field and M is the number of Monte Carlo realizations performed. The 
sample mean of Y(x) for M realizations is 
M N 
p,(Y) = . 1 E E Y(xn)m 
M·N m=l n=1 
(5-1) 
24 
and the sample variance about that mean is 
M N 
Var(Y):= 1 lJ . lJ [Y(xn)m - J.t(Y)]2 (5~2) 
M·N m:=1 n:=1 
With the sample mean and variance known the correlation structure for Y{x) 
is determined in each of the principal directions. Let r i be the separation 
distance between two points in the ~th direction. The correlation between 
all points separated by ri is 
(5-3) 
Var{Y) 
It follows that 
(5-4) 
Once the correlation is determined for numerous separation distances the 
sample correlation length may be determined. As previously defined, the 
correlation length is equal to the separation distance at which the correlation 
is equal to e~l. The correlation structure of Y{x) should be similar to the 
assumed exponential function, so that the sample correlation length should be 
near the assigned correlation length, Ay. If the correlation structure of Y{x) 
is similar to the assumed correlation structure then the turning bands method 
accurately produces a random, spatially-correlated hydraulic conductivity field. 
5.2.2 Velocity 
The velocity field correlation structure is evaluated in much the same 
way as Y{x), except that at each point the velocity has i components, v i(x). 
The sample mean and variance for each of these components, J.t{ v i) and 
Var(v i)' are determined in the same way as for Y(x). In a three-dimensional 
problem the velocity correlation structure has nine components and can be 
represented in matrix form by 
PV1(r1) p~1(r2) PV1(rp,) 
PV2(r1) pv/r2) PV2(rs) 
PV3(r t) pvs(r2) pvs(ra) 
25 
(5-5) 
The component Pv (r2) is the correlation function of Vl(X) in the 2nd 
1 
principal direction. The correlation structure of v i(x) is 
p (r.) = ~~--.--------------
vi 3' (5-6) 
Var(v i) 
where r j is the separation distance between two points in the j-th direction. 
5.2.9 Concentration 
The determination of the concentration correlation structure is different 
from the other two random fields. First, the concentration is 
time-dependent, so that a specific time must be selected to calculate the 
correlation. Second, the sample sb~e (or number of observations of c(x:,t)) is 
not equal to a constant. Some cells will contain zero concentration and these 
are not included in determining the concentration correlation structure, 
pc(r,t). In order to get a large enough sample of nonzero concentration 
values, the time at which the correlation is calculated is chosen to be large. 
The mean and variance of the concentration field, J.L(c) and Var(c), are 
calculated similarly to the other fields. If the total number of cells which 
contain non--zero concentration in the m-th realization is N cm' then the mean 
of c(x, t) for M realizations is 
1 M 
J.L(c) = - E 
M m=l 
and the sample variance about that mean is 
1 M 
Var(c) = - E 
M m=l 
Ncm 
1 . E [c( xn' t) m - /1-( c)] 2 
Ncm n=l 
(5-7) 
(5-8) 
The correlation structure for concentration, pc(r,t), is determined the same as 
(5-3), except that Ncm replaces N. 
26 
Because of potential anisotropy in the hydraulic conductivity, velocity, 
and concentration correlation structures, an anisotropy correlation ratio, IR, is 
introduced which quantifies the anisotropy. This ratio relates the correlation 
length in the mean direction of flow to the correlation length in one of the 
transverse principle directions. For example, in a two-ilimensional problem 
the anisotropic ratios for Y(x), v(x) , c(x,t) are defined as 
IR(Y) ;::: Ay / Ay 2 
IR(V1) ;::: AV IAv IR(V2) = AV IAv 
11 12 21 22 
lR(c) ;::: AC/AC2 
(5-9) 
(5-10) 
(5-11) 
Thus, in a two.....rumensional problem, the velocity will have two anisotropy 
ratios while the hydraulic conductivity and concentration field will have one. 
5.3 Spatial Moments 
This study compares the results from numerical simulations to the 
stochastic theories which describe ensemble, large-scale plume behavior. . A 
popular method for analyzing ensemble plume behavior is spatial moment 
analysis. This method is easily implemented with the particle tracking 
random walk numerical model and spatial moments are used to quantify 
plume velocity, spreading ability, and rate of spreading. Several investigators 
have used spatial moment analysis, including Smith and Schwartz [1980], 
Tompson et a1. [1988], Rubin [1990], Salandin and Rinaldo [1990], and 
Tompson and Gelhar [1990]. 
The particle tracking random walk model simulates an injected mass 
with a large collection of discrete particles. These particles are transported 
according to the steps outlined in Section 4.~,1. Let xP(t), or 
(x~(t),xf(t),x~(t)), represent the position of the pth particle at time t. In a 
large three-dimensional domain the general moment equation is defined as 
1 Np i j k 
MiJ'k(t) = --- lJ [x~(t)] [xf(t)] [x~(t)] 
Np p;:::l 
(5-12) 
where Np is the number of particles injected. As long as no particles have 
left the domain 'the normalized mass, Mooo(t), is equal to one. 
27 
5.3.1 Velocity 
At time t the position of the plume center of mass in each direction, 
xc(t), or (xCl(t),XC2(t),Xcs(t)), is found by setting i,j, or k equal to one in 
(5-12). For example, the center of mass location in the 2nd principal 
direction and for the m-th realization is found by setting j equal to one and 
i and k equal to zero 
(5-13) 
where x~(t)m is the position in the 2nd principal direction of the p-th 
particle in the m-th realization. If M Monte Carlo simulations are performed 
then the ensemble center of mass is ' 
1 
x (t) := -C2 M (5-14) 
The ensemble center of mass in the Xl and Xa directions are found similarly. 
The plume velocity is calculated from the first spatial moment and is 
time rate of change of the plume displacement 
xc(t)m ~ xc(t=O)m 
t 
(5-15) 
where m is the .Monte Carlo realization. The ensemble mean velocity for M 
realizations is 
M 
v(t) := _1_ ~ v(t)m 
M m:=l 
(5-16) 
In the early stages of transport (small t) each realization will yield 
different early time, or small-scale, velocities because the plume has not had 
time to experience all the heterogeneity in the field. However, after a large 
travel time in a weakly stationary heterogeneous medium, an ensemble of 
measurements could be gathered and safe comparisons to theory and field 
results can be made. The asymptotic velocity of the plume is the time rate 
of change of the displacement after the plume has traveled a sufficient 
distance. The ensemble asymptotic velocity v(t-loo). 
28 
, 
.... 
5.9.2 Spreading 
The three-dimensional spread of solute about the center of mass at time 
t is equal to the variance of c(:x:,t) about :X:c(t). In the m-th realization the 
displacement variance in the i-th direction, Var(:x:i(t))m' about the center of 
mass xci(t)m is 
Np 
Var(x .(t)) ::::: --±- ~ [xP(t) x (t) ]2 Z m ~N £.J i m - ci m p:::::l p 
i:::::l,2,3 
The ensemble displacement variance for M Monte Carlo realizations about the 
ensemble center of mass xci(t) is 
M N 1 P p2 
Var(xi(t)) ::::: -- }J }J [xi(t)m - xci(t)] 
M·N m:::::l p=1 p 
i:::::l,2,3 (5-18) 
Although moments higher than second-order may be calculated they are often 
unreliable [Freyberg, 1986] and only the zeroth-, first-, and second-order 
spatially averaged moments (i,j,k:::::O,I,or 2) are used in this study. 
The time rate of change of the variance is not constant in 
spatially-variable porous media. This has been found in field studies 
[Freyberg, 1986; Gelhar and. Axness., .. 1983] and theoretical computations 
[Gelhar and Axness, 1983; Dagan, 1984]. Using a linearized, first-order 
approximation and an isotropic exponential covariance function for the 
logarithm of hydraulic conductivity, Dagan (1984) found expressions for the 
time dependence of the variance for two- and three-dimensional steady-state 
flow. The time-dependent relationships in two dimensions are 
Var(xl(t)) ::::: 0"-f.'xy.{1.5 + 2·t + 3·[Ei(-t) -
In( t) - 'Y] + 3· t-2• ( e-t . (1 + t) - I)} (5-19) 
Var(x2(t)) - 0"-f.'xy.{-1.5 - [Ei(-t) -
In( t) - 'Y] - 3· t-2• ( e-t . (1 + t) - I)} (5-20) 
where Ei( ) is the negative exponential integral, 'Y is the Euler constant 
(0.5772), and t is a dimensionless characteristic time 
29 
t ';::; t·v d t--+oo) 
Ay 
(5-21) 
where Vl(t--+OO) is the ensemble asymptotic plume velocity in the direction of 
mean flow. 
The theoretical predictions of Dagan [1984] have been substantiated by 
tracer tests and by numerical experiments. Freyberg [1986] compared the 
longitudinal and transverse displacement variances from the Borden site to 
Dagan'S theoretical predictions and found good agreement between the two. 
Numerical experiments performed by Salandin and Rinaldo (1990) show that 
Dagan's linear theory yielded acceptable results for a large range of variability 
in tnK. As has been shown by the theoretical studies, some numerical 
simulations, and especially a couple of field experiments, the rate of spread of 
a solute plume is not necessarily constant, i.e. the transport is non.-Fickian. 
It has been found though that after a long travel time, or distance, transport 
can become Fi cki an. 
5.3.3 Macrodispersivity 
The displacement variances calculated by spatial moments can be used 
to detect large-scale dispersive behavior. The macrodispersivity parameterizes 
the spreading caused by large-scale heterogeneity of the hydraulic 
conductivity. In terms of spatial moments, the macrodispersivity is found 
from the expression [Freyberg, 1986] 
Var(xi(t)) .... Var(xz{t';::;O)) 
A~t) ';::; -------
Z 2.t'Vt(t--+oo) 
where ai the local dispersivity (aL or aT)' 
i ,;::;1,2,3 (5-22) 
Of particular interest ~s the asymptotic macrodispersivity Ai(t--+oo). 
Gelhar and Amess [1983] solved the stochastic flow and transport equations 
using analytical spectral theory, and related Ai(t--+w) to the parameters of the 
hydraulic conductivity field. For example, for two..-dimensional flow parallel 
to layering, and assuming an isotropic exponential correlation function for 
lnK, the expressions for longitudinal and transverse asymptotic 
macrodispersivities are 
(5-23) 
30 
uy aL [ 1 + 3(aT/aL) ] 
A2( t--+ro) :::; ~'---~~--'---'----
8 ¢2 
(5-24) 
where ¢ is the ratio between the effective hydraulic conductivity to the 
geometric mean 
(5-25) 
Ke depends on the number of dimensions in the problem, In two-
dimensional numerical simulations ¢ is equal to one. Note that Al(t--+ro) is 
only dependent on uy and Ay, while A2(t..-+ro) is dependent on the local 
dispersivities but not Ay. This indicates that longitudinal dispersion is 
controlled by convection but the transverse dispersion is dependent on local 
dispersion [Ge1har and Axness, 19831. The second-order spatial moments 
determined from the numerical simulations will be compared to both the 
time-dependent variances predicted by Dagan [1984] (5-17 and 5-18) and the 
asymptotic macrodispersivity values predicted by Gelhar and Axness [1983] 
(5-23 and 5-24). 
5.4 Temporal Moments 
The second method for describing solute movements in heterogeneous 
porous media is to consider the solute arrival at a fixed location. The 
arrival time provides an interpretation of longitudinal solute movement that is 
complementary to spatial moments [Shapiro and Cvetkovic, 19881. These 
solute arrival times are displayed graphically in breakthrough curves. 
Breakthrough curves help characterize uncertainties in transport predictions 
rSmith and Schwartz, 1981]. If Monte Carlo simulations are performed the 
hreakthrough curves can be used to show worse case scenarios of transport 
[Dagan, 1989]. Breakthrough curves will be developed for each numerical 
simulation. 
31 
5.4.1 Breakthrough Curves 
By using Monte Carlo simulations, arrival times at a particular location 
can be analyzed in terms of their means and standard deviations. Let the 
time it takes for a certain percentage, (P::::1%, 25%, 50%) 75%, 99%) of 
particles to arrive at plane x::::xp be denoted Tm(xp'P), where m is the mth 
of M Monte Carlo simulations. Then the average time it takes P percentage 
particles to arrive at plane xp is 
M 
T(xp'P) :::: _1_ E Tm(xp'P) 
M m;::::l 
The standard deviation about this mean time is 
M 
ST(xp'P) ;:::: (--L E {T (x ,P) - T(x ,P)]2 )1/2 M m::::1 m p p 
(5-26) 
(5-27) 
The mean arrival time depends on O"y. For example, as O"y increases the 
expected arrival time of the first particle, T(xp'O%) decreases because the 
solute follows preferred paths of higher conductivities [Smith and Schwartz, 
1980]. Conversely, the breakthrough times for P>20% increase for higher O'y 
because the solute encounters zones of lower conductivity [Shapiro and 
Cvetkovic, 1988]. Thus) for higher O'y, the slower particles move more 
slowly while the faster particles move more quickly. This behavior probably 
reflects the skewed log normal hydraulic conductivity distribution. Further, 
as the spatial continuity of the hydraulic conductivity increases (higher Ay) 
the mean arrival time decreases, probably due to the more extensive low 
hydraulic conductivity zones and log normal distribution [Smith and Schwartz, 
1980]. The variability of arrival times, ST(xp'P), increases by increasing O'y 
and/or Ay , while later positions of the breakthrough curve (higher P) are 
more uncertain than earlier positions [Smith and Schwartz, 1980]. 
32 
-I 
CHAPTER 6. NUMERICAL SIMULATIONS 
6.1 Introduction 
The concepts developed in previous sections are now incorporated into a 
large-scale heterogeneous flow system in order to investigate the effects that 
variability in hydraulic conductivity has upon flow and contaminant transport. 
Some limitations have been placed upon the scope and size of the problem. 
First) the problem is limited to transport of conservative solute through a 
two-dimensional, steady-state saturated flow region. The hydraulic 
conductivity field is distributed log-normally and is weakly-stationary) so that 
the first two moments characterizing lnK(x) are everywhere the same in the 
flow region. The exponential correlation function is assumed for lnK(x). 
The porous media is assumed to be isotropic. The three parameters which 
quantify the heterogeneity of K(x) are ry' O"y) and Ay The mean is held 
constant. The effects of the degree of heterogeneity and the spatial 
continuity of lnK(x) upon flow and transport are investigated. A total of 
thirty Monte Carlo simulations are performed. The contaminant is injected 
at a point source instantaneously. Figure 6-1 provides a flow diagram of the 
step by step procedures for the numerical simulations. 
6.2 Physical Domain 
A domain size has been selected according to the recommendations of 
Ababou et al. [1988] and (2-11). In particular the size of the representative 
elementary volume (REV) is much smaller than the correlation length Ay, 
and Ay is much smaller than the size of the domain in both directions. The 
size of the REV is the same everywhere in the domain and is equal to 
Llx=Llxl=Llx2=1. All length units in the remainder of this study are meters 
and the time units are in days. Two correlation lengths have been 
investigated: Ay =5m and Ay =10m. The size of the domain is held constant, 
Lf==200m and L2=100m, so that there are 20,000 cells. Therefore, when 
Ay =5m, the dimensionless length scale ratios outlined in Section 2.4.1 are 
33 
Fig. 6-1 
ITurning Band \ 
1 
Z(x) ~ ~(O,1) 
! 
I{(x) = exp[ l1-y + Z(x) oy] 
J 
I Galerkin F.E·I 
! 
V(x) 
! 
P article Tracking Random Walk I 
1 
O(x,t) 
1 I Moment Analysis I 
Flow diagram for numerical methods used 
in the numerical simulations. 
34 
Ay/Ax;::: 5 L1/ Ay = 40 (6-1) 
and when Ay=lOm they are 
(6-2) 
The hypothetical domain is shown in Figure 6-2. The geometric mean of the 
hydraulic conductivity, Kg, has been set to 8.64 m/day ((10t4 m/sec), so 
fly= -9.21. For each Ay, three different ay are used: 0.5, 1.0, and 1.5. 
Table 6-1 outlines the six different sets of K(x) fields to be used. 
6.2.1 Hydraulic and Transport Parameters 
In the flow problems, unconfined flow is developed from specifications of 
a mean hydraulic gradient. Constant head boundaries are chosen at the right 
and left boundaries and are 
(6-3) 
The top and bottom boundaries (X2=0.Om and x2=100.0m) are free-flow 
boundaries. The mean hydraulic head gradients in each direction are 
-J 1 = p(xl=200) - p(xt=O) _ 18 - 20 / L1 - 200 = -0.01 m m (6-4) 
J2 ;::: 0.0 m/m 
which result in a mean flow in the Xi direction. The porosity is n=0.144 so 
the mean pore-water velocity is 
V i (x) 1 - 1 = - n J 1 Kg = - 0.144 (-0.01) (8.64) = 0.6 m/day (6-5) 
v 2 ( x) = 0.0 m/ day 
All transport simulations start with an initial injection of 5000 particles 
at a point. At the beginning of each simulation the position of each particle 
p is x~(t=0)=4.5m and x~(t;:::0)=49.5m, which is in the middle of a cell. 
The longitudinal and transverse local dispersivities are aL=0.2m and 
aT=0.02m. 
35 
~ 
0 
<t 
~ 
:r: 
S E-o 
0 Z 0 < 
.-1 E-o 
rn 
Z 
c:,;:, 0 
0') u 
FREE FLOW 
LONGITUDINAL DIRECTION )0 
TRANSVERSE DIRECTION 
r 
SOURCE 
• 
X2 
Xl 
Fig. 6-2 
.. 
DIRECTION OF MEAN FLOW 
FREE FLOW 
200m 
Size and shape of flow domain, including injection point, 
boundary conditions, and mean flow direction. 
0 
<t 
til 
:r: 
E-o 
Z 
< E-o 
rn 
Z 
0 
u 
Problem Set 1 
fly 
Ay 
2 O'y 0.5 
, 
\ 
TABLE 6~1 
2 3 
I 
-9.21 
5 
1.0 1.5 
37 
4 5 6 
10 
0.5 1.0 1.5 
Each transport problem is continued until at least 99% of the particles 
exit the domain through the right boundary. Unfortunately, it is 
economically unfeasible to continue the transport until all particles leave the 
domain because some particles may be caught in vast regions of low 
conductivity. The particles in each Monte Carlo realization are displaced 
according to (4-5), (4-6), and (4-7). The time step is chosen to be 1 day 
so that the Courant Number (4-8) is less than one 
C -_ vl(x).ilt c Ax;:::: 0.6 (6-6) 
6.3 Correlation Structure 
The methods for analyzing the correlation structures of the random 
hydraulic conductivity, velocity, and concentration fields were outlined in 
Section 5.2. These fields are analyzed in two directions; the direction of the 
mean flow (Xl), or the longitudinal direction, and the direction perpendicular 
to the mean flow (X2), or the transverse direction. Because velocity has two 
components, Vl(:X:) and V2(:X:) , the correlation structure of each component in 
both the longitudinal and transverse directions is evaluated. 
6.3.1 Hydraulic Oonductivity 
The isotropic correlation structure of the hydraulic conductivity is 
parameterized by >'Y and follows the exponential distribution (2-9). Figures 
6-3 and 6..-4 show the computed longitudinal and transverse correlation 
structures of lnK(:x:) for >'y=5m and >'y=10m, respectively, along with the 
assumed exponential distribution. The x-axis in these correlation figures and 
all that follow are in terms of a dimensionless lag. This dimensionless lag 
represents separation distance normalized by the correlation length. For 
example, in Figure 6-3, a dimensionless lag of 2 is the same as a separation 
distance of 10m (2. >'y;::::10),' while in Figure 6..-4 this same dimensionless lag 
is the equivalent of a separation distance of 20m, since here >'y;::::10m. These 
figures indicate that the turning bands method accurately reproduces an 
isotropic correlation structure, though the accuracy is lower for >'y:::=10m .. In 
each figure, the longitudinal and transverse correlation lengths are equal,' 
so that the anisotropy correlation ratio for inK '(5-9) is IR{Y):::=1. 
38 
,.-... 
s... 
-
:>-
~ 
-Z 
0 
-j 
til 
c:;., 0:: 
~ 0:: 
0 
u 
~ 
1.0 
0.9 
0.8 
0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 
r~ 
~ 
o py{rt ) 
o py(r2) 
Exponential 
Ay =5m 
0.0 
~ 
- ~ 
-------------- ~ ~ ~ 
-0 1 
• I 
0 1 
Fig. 6-3 
2 3 4 5 
NORMALIZED DISTANCE, r lAy 
Comparison between the longitudinal and transverse 
correlation structures for 1nK and the assumed exponential 
correlation function, ()'y-5m }. 
6 7 
_-1~~~-L-
i.o 
0.9 
0.8 
0.7 
-s.. 0.6 ....., 
ct 
. 
z 0.5 0 
-E-< j 0.4 
I'.::rl 
0:: 
~ 0:: 0.3 0 0 
u 
0.2 
0.1 
0.0 
-0.1 
I 
0 1 
Fig. 6-4 
o py{rt } 
o py(r2) 
Exponential. 
Ay =10m 
=s 0 ~ ~ - -~ - e 0  
2 3 4 5 6 
NORMALIZED DISTANCE, r/')...y 
Comparison between the longitudinal and transverse 
correlation structures for fuK and the assumed exponential 
correlation function, (Ay-lOm). 
7 
6.9.2 Velocity 
The correlation structure of the two components of the random velocity 
field are each evaluated in the longitudinal and transverse directions. Thus, 
four correlation structures are determined and the correlation of the velocity 
can be expressed in matrix form as 
PV1(r1) PV1(r2) 
PV2(r1) PV2(r2) 
Figures 6-5 and 6-{) display the four components of the velocity correlation, 
for Ay:=5m and uf:=1.0.· Figure 6-5 shows that an isotropic hydraulic 
conductivity correlation structure produces an anisotropic longitudinal velocity 
correlation structure. The correlation length of the longitudinal velocity in 
the longitudinal direction, AVU' is approximately 2.4· Ay, or 12m, while the 
correlation length of the longitudinal velocity in the transverse direction, AV ' 12 
is approximately 0.8· Ay, or 4m. The correlation length is determined by 
finding the separation distance at which the correlation is equal to e~l. 
anisotropic ratio for the longitudinal velocity (5-10) is 
lR(v) ~ 12m ~ 3 
1 - 4m -
The 
(6-8) 
The transverse velOCity correlation structure is displayed in Figure 6-6. The 
correlation length appears to be approximately the same in both directions, 
AV =AV £0.8· Ay:=4 m, and the anisotropic ratio is 21 22 
(6-9) 
These results from numerical simulation are in close agreement with the 
study of Rubin [19901, who analytically derived the velocity covariance 
function using the smal perturbation approach. The velocity anisotropy was 
also noticed by Graham and McLaughlin [1989]. 
The . longitudinal and transverse velocity correlations for Ay=10m are 
displayed in Figures 6-7 and 6-8. Figure 6-7 shows that the longitudinal 
velocity correlation length in the longitudinal direction is approximately 
1.8· Ay, or 18m, while in the transverse direction it is approximately 0.8· Ay, 
or 8m. The anisotropy ratio for Ay=10m is 
41 
1.0 
0.9 
0.8 
0.7 
-
s.. 
- 0.6 -> Q. 
z 0.5 
0 
-~ 0.4 
tzl 
~ 
,;0.. 0: 0.3 ~ 0 
u 
0.2 
0.1 
0.0 
-0.1 
I 
0 2 
Fig. 6-5 
----- PVl (r t) 
e-e PVl(r2) 
Ay =5m 
a~= 1.0 
- ---- ------------
4 6 8 10 12 
. NORMAUZED DISTANCE, rj'Ay 
Correlation structure of the longitudinal velocity in . the 
longitudinal and transverse directions~ (Ay =5m, cri-LO). 
14 
1.0 
0.9 
0.8 
0.7 
-s.. 
-
0.6 CIt 
> Q. 
z 0.5 
0 
-e-. j 0.4 
ril 
~ 0:: 
~ O! 0.3 0 
u 
0.2 
0.1 
0.0 
-0.1 
I 
0 2 
Fig. 6-6 
------ PV2(rt ) 
8--e PV2(r2) 
Ay =5m 
a~= 1.0 
-- e=---- .--'" ..--===i 
..... . .---- ---. 
4 6 8 10 12 
NORMALIZED DISTANCE, r/'Ay 
Correlation structure of the transverse velocity in the 
longitudinal and transverse directions , ().y-5~,(Jf-1. 0). 
14 
----------------_._._-_._ .. 
1.0 
0.9 
0.8 
0.7 
,..... 0.6 s.. 
--> Q.. 0.5 
. 
z 
0 0.4 
-E-c j 
0.3 til 0:: 
Jo!::>- ~ 
Jo!::>- 0 0.2 u 
0.1 
0.0 
-0.1 
-0.2 I 
0 
" 
B--II PVl(rt ) 
--- PVl(rZ) 
Ay=10m 
a~= 1.0 
----'\.- --- ------
1 
Fig. 6-7 
2 345 6 7 
NORMALIZED DISTANCE. r lAy 
Correlation structure of the longitudinal velocity in the 
longitudinal and transverse directions, (>'y=10m, O"t= 1. 0). 
8 9 
1.0 
0.9 
0.8 
0.7 
,... 0.6 s.. 
-
Of 
> ~ 0.5 
. 
z 
0 0.4 
-E-o j 
rz:l 0.3 
~ ~ 
.~ 0:: 
0 0.2 u 
0.1 
0.0 
-0.1 
-0.2 
'\ .... 
I 
0 1 
Fig. ·6--8 
2 3 4 5 6 
NORMALIZED DISTANCE. r lAy 
------- PV·ir l) 
.-. PV2(r2 ) 
Ay=10tn 
a~=1.0 
7 8.' 
Correlation structure of the transverse velocity in the 
longitudinal and transverse directions, (Ay-lDm, (Trl.D). 
9 
which seems to be smaller than the anisotropy ratio for the case of Ay=5m 
(IR(Vl) ~ 3). The smaller anisotropy ratio for the case of Ay=10m is a result 
of decreasing longitudinal correlation length (Ay=5m q Avu=2.4. Ayi Ay =lOm 
q A =1.8· Ay), while the transverse correlation length, remains unchanged 
V11 
(AV12=0.8.Ay), both relative to the Ay .. The smaller IR(Vl) for larger Ay is 
the result of the mean flow damping the velocity correlation in the 
longitudinal direction, as Ay gets larger. The correlation structures for the 
transverse .veloc.ity are shown in FiaUIe 6-.8. The structures appear similar, 
so that the anisotropic ratio for V2(X) is 
(6-11) 
just as for Ay=5m. 
The longitudinal and transverse velocity correlation structures in the Fig 
6-7 longitudinal direction (p (rt) and p (r1)) were determined for the cases 
VI V2 
of uf=O.5 and uf=1.5, and it was found that uf has no effect on the 
correlation length (Figures 6-9 and 6-10). This may be because the same 
realization of Z(x) is used to create K(x), by rescaling uf in (4-3). 
Independent hydraulic conductivity fields, with different uf, may show pv(r) 
to be dependent on uf. The analysis of the velocity correlations indicate 
that the velocity covariance functions are anisotropic when hydraulic 
conductivity fields used to generate the velocity fields are considered isotropic. 
This is because the constant head boundaries establish a preferential direction 
of water flow in the longitudinal direction. 
It is important to emphasize that the velocity correlation structure 
exhibits the hol~ffect in three [pv (ri), p (rt), p (r2)] out of the four 1 V2 V2· 
correlation functions of the velocity field (the only one which does not show 
the hole-effect is p (rt)). By definition the hole-effect has occurred when V1 
the spatial correlation between two points becomes negative. A directional 
hol~ffect has been witnessed in some mineral deposits [Journel and 
Huijbregts, 1978]. A velocity field hole-effect signifies a velocity field which 
periodically repeats itself with clusters of large and low values. For example, 
the hole-effect for the longitudinal velocity in the transverse direction 
represents conduits of high velocity alternating with sections of low velocity. 
The impact of the velocity hole-effect upon solute transport is not clear at 
this time. 
46 
/----------------------
1.0 
0.9 
0.8 
0.7 
--$.0 0.6 
--> Q. 
z 0.5 
0 
-~ j 0.4 
ril 
JolS.. 0:: 
~ 0:: 0.3 0 
u 
0.2 
0.1 
0.0 
-0.1 
I 
0 
Fig. 6-9 
2 4 6 8 10 
NORMALIZED DISTANCE. r1/AY 
[3--£J a~= 0.5 
G-€> a~=1.0 
I:r----A a~= 1.5 
Ay =5m 
12 
Correlation structure of the longitudinal velocity in the 
longi~udinal direction, (Ay-5m,CTf-O.5,1.O,1.5). 
14 
1.0 
0.9 
0.8 
0.7 
~ 
-s.. 0.6 
-
Dr 
> Cl. 
z 0.5 
0 
-E-t 
:s 0.4 
~ 
~ 0:: 
00 0:: 0.3 0 
u 
0.2 
0.1 
0.0 
-0.1 
I 
0 2 
Fig. 6-10 
4 6 8 10 
NORMALIZED DISTANCE. rt/AY 
B-EJa~=O_5 
G--E) a~=-1_0 
~ a~=1_5 
Ay=5m 
12 
Correlation structure of the transverse velocity in the 
longitudinal direction, CAy-5m, CTf=0.5,1.0,1.5). 
14 
6.9.9 Concentration 
Because the concentration distribution is tim~ependent a specific time 
is chosen to calculate the correlation structure. The largest possible time was 
selected in order to allow the plume to spread out sufficiently, yet small 
enough to still have all the particles in the domain. It was found that 180 
days allowed this to happen. The particle distribution of all 30 Monte Carlo 
realizations are collected after 180 days and the concentration correlation 
structure is calculated according to the procedures outlined in Section 5.2.3. 
Figure 6-11 displays the concentration correlation structures in both the 
longitudinal and transverse directio~s, for Ay=5m and O"y=1.0. The 
correlation lengths in each direction are approximately '\1=4Ay=20m, and 
AC2=1.3Ay=7m. The anisotropic ratio (5-11) is 
lR(c)!Y 20m t!! 3 
- 7m- (6-12) 
This ratio is the same as IR(Vl), but the concentration correlation lengths 
nearly double, from 12m to 20m in the longitudinal direction, and from 4m 
to 7m in the transverse direction. These longer correlation lengths indicate 
that the distribution of concentration is more spatially continuous than either 
the hydraulic conductivity or the velocity. The correlation of concentration 
in the transverse direction exhibits the hole-dfect. This hole--effect indicates 
that the contaminant plume is made up of alternating longitudinal rows of 
high and low concentration. These areas cannot be located in the same 
regions of the high and low velocity conduits because the correlation lengths 
are different. . 
The differences in the Y(x), Vl(X), and c(x,t=180) correlation structures 
in the longitudinal direction are shown in Figure 6-12; with Ay=5m. All 
three exhibit positive correlation. The correlation lengths increase from 5m 
for Y(x) , to about 12m for V1(X), and to approximately 20m for c{x,t=180). 
The isotropic hydraulic conductivity field produces anisotropic velocity and 
concentration correlation structures. The anisotropic ratio for the 
longitudinal velocity is approximately 3, regardless of the correlation length of 
inK(x). The anisotropic ratio for concentration remains the same as Vl{X), 
about 3. 
Both the longitudinal and transverse velocities display a hole-effect in 
the transverse direction. The concentration correlation structure also displays 
a hole--effect in the transverse direction. The correlation structures for the 
longitudinal velocity, transverse velocity, and concentration field are all 
similar in the transverse direction. 
49 
1-----
.:3 TWYtfE " n-
1.0 
0.9 
0.8 
0.7 
-
0.6 s.. 
-~ 0.5 . 
z 
0 
- 0.4 ~ .. 
rzl 0.3 0::: 
Crt 0::: 
'0 0 
.U 0.2 
0.1 
0.0 
-0.1 
-0.2 
I 
0 2 
Fig. '6-11 
4 6 8 10 
NORMALIZED DISTANCE, r /'Ay 
----- pc(rl,l= 180) 
--- pc(r2 ,l=180) 
Ay =5m 
a~=1.0 
12 14 
Longitudinal and transverse correlation structures for 
concentration at 180' days, (.Ay =5m, O"f=LO). 
""":: 
s... 
........ 
Q" 
. 
z 
0 
-E-t 
:s 
tzl 
Of .~ 
....... a:: 
0 
t,) 
1.0 
0.9 
0.8 
0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 
~ ----.-
py(rt} 
PVl (r1) 
pc(rh l=180) 
Ay=5rn 
0.0-- - - - - - - -
-O.l~--~--~--~~--~--~--~--~--~--~~--~--~---r~ 
o 2 
Fig. 6-12 
4 6 8 10 12 
NORMALIZED DISTANCE. rt/'Ay 
Longitndinalcorrelation structures for Y(x) , Vl(X) , and 
c(x,t=180) , (Ay=5m,uf--l.O). 
14 
6.4 Spatial Moments 
The methodology for spatial moment analysis was outlined in Section 
5.3. Spatial moments are used to investigate large-scale, ensemble plume 
behavior, not point predictions. This method can help evaluate the 
assumptions of small perturbation and ergodicity, on which much of transport 
theory is based. Small perturbation is the foundation of the theoretical 
developments by Gelhar and Axness [1983] and Dagan [1984]. Small 
perturbation assumes that the level of heterogeneity, (Jy, is small. In most 
cases it is assumed to be smaller than one. Although the stochastic theories 
do not apply when (Jy is equal to 1.5, this study will compare the theoretical 
predictions to the numerical results. The purpose of the ergodic hypothesis 
was outlined in Section 2.4. When the ergodic hypothesis. is fulfilled, 
conclusions reached from numerical simulations are general and not dependent 
on a particular realization. By increasing the correlation length of lnK(x), 
and thereby decreasing the number of independent samples in the domain, 
this hypothesis can be investigated. The ensemble plume characteristics most 
often determined from spatial moment analysis are the plume velocity, 
spreading, and macrodispersivity, The first is found from first-order moments 
of the plume and the last· two are determined from the second-order 
moments. 
6.4.1 Velocity 
The mean movement of the plume center of mass is determined from 
(5-14). Figure 6-13 gives an example of the position of the center of mass 
as a function of time. The x-axis (Xl) is in the direction of mean flow and 
the y-axis (X2) is perpendicular to the mean flow. The axes pass through 
the point in the domain where the particles are injected, x~(t=0)=4.5m and 
x~(t=0)=49.5m. Each data point represents the mean position of 30 Monte 
Carlo realizations at a specific time (t=1,5,10,20,40,60, ... ,200 days). In Figure 
6-13 the correlation length of lnK(x) is 5m and (JY=1.5. This figure shows 
that the center of mass does not deviate substantially from its expected path, 
which is in the mean flow direction. The small deviation which does exist is 
a result of a finite number of realizations. If the number of Monte Carlo 
realizations were increased, the movement of the mass center would more 
closely follow the direction of mean flow. By decreasing the degree of 
heterogeneity the center of mass moves slightly less in the transverse 
direction. By increasing Ay to 10m, with the same (Jy, there is no 
difference in the position of the center of mass. Thus, Ay and cry have no 
effect on the position of the center of mass in these numerical experiments. 
52 
-E 
-til 
X 
'at 
~ 
10 
8 
6 
4 
2 
0 
-2 
-4 
-6 
-8 
_. __ .-----
Ay=5m 
(T~ = 1.5 
-10 I ' I 
o 10 20 30 40 50 60 70 80 90 100 110 ·120 
Fig. 6-13 Horizontal trajectory o~Xth~~e center of mass, CAy-Sm, 
0"f-~.5). 
Because there is little displacement in the X2 principal direction, the 
total mean displacement is assumed to equal the displacement in the Xl 
principal direction 
(6-13) 
Figure 6-14 shows the mean displacements for Ay:=10m. The slope of these 
lines is the ensemble plume velocity. This figure indicates that the ensemble 
plume velocity is nearly constant throughout the simulation. An increase in 
uy has very little effect on the average displacement of the mass center or 
the ensemble plume velocity. This result is the same for Ay=5m. 
From the constant head boundary conditions the average pore-water 
velocity imposed was v 1 ( x) ~0.6 m/day and v 2 ( x) :=0.0 m/day. After 
200 days of transport, or 120m, the experimental velocity ranged from 0.56 
m/day for Ay:=5m to 0.54 m/day for Ay:=10m. However, the average 
velocity for all simulations after each plume had traveled a large time (after 
100-120 days) is apprOximately equal to the assumed velocity. Thus, the 
ensemble asymptotic plume velocity is equal to the imposed pore-water 
velocity 
V 1 ( t 4(0) = 0.6 m/day V 2 ( t 4(0) = 0.0 m/ day (6-14) 
6.4.2 Spreading 
The spreading of the plume is described by the displacement variance 
about the center of maSs. These are determined in both the longitudinal 
and transverse directions by (5-18). The longitudinal displacement variances 
Var(xl(t)) are shown in Figure 6-15. This figure is for Ay:=5m and contains 
the longitudinal displacement variances for each of the three degrees of 
heterogeneity. As the variance of inK (x) increases the displacement variance 
increases proportionally. After 120 days all three plumes appear to have 
reached a constant rate of increase, which indicates the rate of spreading of 
the plume in the longitudinal direction has reached a constant. Dagan [19841 
developed an expression for Var(xl(t)), (5-19), which is dependent on Ay and 
uf. Figure 6-16 compares the theoretical predictions from Dagan [1984] 
(5-19) to numerical results for C1f:=O.5 and C1f=1.5 , and Ay:=5m. Both 
numerical results show fairly good agreement at early times but begin to 
depart from theory later. Because Dagan's theory is based on the assumption 
_=ri 
-6 
......... 
E-o 
Z 
tEl 
::s 
tEl 
U 
:s 
~ Il. 
Cfl rn 
-Q 
120 
G-EJ a~=O.5 
G--€) a~= 1. 0 
100 --1 /:r---f). a~= 1.5 
80-l Ay= 10m 
60 
40 
20 
o F 
o 20 40 60 80 100 120 140 160 180 
Fig. 6-14 
TIME (days) 
Horizontal displacement of plume center of mass as a 
function of time since injection, (Ay =10m, O"f=O.5,l.O,1.5). 
200 
-at 
0) 
II---B u~=O.5 
--- u,=l.O 
15001 +-----+ u,= 1. 5 
-
.. 
,S 
-ri:I j / Co) 'A,y=5m Z 
< 
-~ 1000 
> 
...l 
< 
Z 
-C 
::> 
e-
-0 
. :..J.Z 500 9 
°r=oe-=r= 
o 20 
Fig. :6-15 
40 60 80 100 120 140 160 180 
TIME (days) 
Longitudinal displacement variances' about plume center, 
{Ay ,5m, uf-O.5,l.O,1.5}. 
• 
200 
~ .... tmt 
CTt 
-;r 
• a~=O.5 
• a~=1.5 
1500-1 Dagan [1984] 
-l 
-'" e 
-C3:l 
u J Ay=5m z 
< 
-~ 1000 
> /+ 
...l 
< 
Z 
- ~ Q / • ::> 
E-r 
-
1 /" 
.c." 
z 500 0 
...l 
o~~ 
o 20 
Fig. 6-16 
40 60 80 100 120 140 160 
TIME (days) 
Comparison of numerical and theoretical [Dagan, 1984] 
longi~udinal displacement variances, (Ay -5m, cr-Q--=O.S,L5). 
• 
• 
• 
• 
180 200 
of small perturbation, a value of O'Y=1.5 may be too high a perturbation. 
Because the spreading predicted for low perturbation (O'Y=0.5) is also too low 
after large times of transport, the difference may be the result of the 
hole-effect witnessed in the velocity fields (Section 6.3.2). This hole..-effect 
may have an impact on the transport process after large times. 
The effect upon the longitudinal displacement variance by increasing the 
correlation length of lnK(x) to 10m is displayed in Figure 6-17. The 
spreading of the plume, for a given O'Yl increases by increasing Ay. Unlike 
the case for Ay=5m, it does not appear that the slopes of these lines have 
reached a constant asymptotic value. Figure 6-18 shows a comparison 
between Dagan's [1984] theoretical results and the computational results for 
Ay =10m. As with Ay =5m, the graph with O'Y=0.5 agrees well at early 
times but the two lines diverge after large travel distance. For O'Y=1.5 a 
pronounced decrease in longitudinal displacement variance exists throughout 
the simulation. After 100 days the rate of spreading appears to be the same 
as the theoretical rate, for O'Y=1.5. 
The displacement variances are now analyzed in the transverse direction. 
Figure 6-19 shows the transverse displacement variance about the mass center 
for the three degrees of heterogeneity. These variances, and the spreading, 
are about ten times smaller than the longitudinal displacement variances for 
each case. The spreading increases nearly linearly until 160 days and then 
levels off to a near constant value. In Figure 6-20, two cases of O'y are 
compared to Dagan's [1984] solution. While the numerical results are 
comparable for small perturbation the computational results are much higher 
after 100 days for the larger perturbation. However, after a large travel time 
the rate of increase in the transverse displacement variances does appear to 
level off to a low, constant value. 
The effect upon the transverse spreading by increasing Ay is now 
analyzed. In Figure 6-21 the transverse displacement variance is shown for 
each of the three degrees of variability for Ay=10m. Since the variance 
should be proportional to O'f it appears in the case of O"Y=1.5 the numerical 
results are substantially smaller than that predicted. In Figure 6-22, Dagan's 
[1984] solution is compared to two of the three cases, O'Y=0.5 and O'Y=1.5. 
For the lower perturbation there is good agreement but not exact. However, 
as with the spreading in the longitudinal direction for O'Y=1.5, the transverse 
variance is substantially over-predicted by Dagan's theory. 
58 
- -- -- ---~-----~------- - ----~- ---- ~--------~---.------~- ---~--
f .... """--~--~_"'O"'::"""-
• 2000 -l .-. a~= 0.5 / .-. a~= 1.0 
/ 
/ +---+ a~= 1.5 
/ .......... Of 5 1500 
• // ri1 u Z 
.+ .• 
< 
-
Ay =10in ~ -< 
> 
...l 1000 
< z 
-Cit C 
to ::> 
E-< 
-a 
z 
0 500 ...J 
o 1""-=-·7 
o 20 40 60 80 100 120 140 160 180 200 
TIME (days) 
Fig. '6-17 Longitudinal displacement variances about plume center, 
(Ay =10m, O"Y=O.5,1.0,L5). 
2500 ---1 
• a~=O.5 
• a~=1.5 
Dagan [1984] 
.... 2000J ./ • 
e 
--~ l ./ • 
u 
z 
~ 1500 I 
./ • 4( 
> 
....J 
./ i 1000~ Ay=10m • 0') /' • 0 • E-o 
-0 
z 
9 
500 
o r=~ 
o 20 40 60 80 100 120 140 160 180 200 
TIME (days) 
Fig. 6-18 Comparison of numerical and theoretical [Dagan, 1984] 
longitudinal displacement variances, (Ay-10m, a-Q-0.5,1.5). 
Q:) 
l--"'-
-
N 
e 
-
140 -1 -------- a~= a . 5 
----- a~= 1. a 
+-----. a~= 1. 5 
. 120 
r:rl 100 
u 
z 
< 
-~ 80 < > ~ r:rl Ay=5m 
CI) 
0:: 60 r:.::J 
:> 
CI) 
z 
< ~ 40 ~ 
20 
o r"'--
o 20 
Fig .. 6-19 
/ // 
~/ 
/~~ 
.. ~ 
40 60 80 100 120 140 160 
TIME (days) 
Transverse displacement variances about plume center, 
('xy .Sm, CTf=O.S,l.O,l.S). 
180 200· 
r-' -,--
-
.. 
e 
-Cz1 
U 
Z 
< 
-0:: 
< 
> 
Cz1 
en 
0:: 
f;iI 
~ > to.:> en 
Z 
< 0:: 
e--
140 • a'=O.5 
• a'=1.5 • 
Dagan [1984] 
120 
• 
100 
• 
80 
Ay=5m 
60 
40 
~ 
20 
o r-= 
o 20 
Fig. 6-20 
40 60. 80 100 120 
TIME (days) 
140 160 
Comparison of numerical and theoretical [Dagan, 1984] 
transverse displacement variances, (Ay-5m,uf-0.5,1.5). 
• • 
180 200 
-.... 
i
i






IiIIIiiiIirIi _______ . ___ ._._ . ..,--______ _ 
d) 
c:.o 
140-1 --- u~=O.5 
----. a~= 1. 0 
+---+ u~= 1.5 
120 
~ 
e 
-r:.:l 100 
u 
z 
< 
-~ 
~ 
IZ1 
rn 
~ 
r:.:l 
> til 
Z 
na 
e-. 
80 
~ Ay=10m 
60 
40 
20 
°re~ 
o 20 
Fig. 6-21 
40 60 80 100 120 140 160 
TIME (days) 
Transverse displacement variances about plume center, 
(Ar 10m, O"f-O.5,l.O,1.5). 
180 200 
o In r?*::: -
o 20 
Fig. :6-'22 
40 60 80 100 120 
TIME (days) 
140 160 
Comparison of numerical and theoretical IDagan, 1984] 
trans:-erse displacement variances; (Ay -10m, lTf-O.S,L5). 
180 200 
6.4.3 Macrodispersivity 
The slope of the variance in each direction is parameterized by the 
longitudinal and transverse macrodispersivities, Al(t) and A2(t), respectively. 
The macrodispersivity is proportional to the time rate of change of the 
displacement variance and it is determined from (5-22). Gelhar and Axness 
[1983] developed expressions £,or the asym,pto,tic l,ongitudinal and transverse 
macrodis, persivities, A1(t, -too) and A2(t-too), (5-,23) and ~5-24)", respec,tive,ly. In 
Figure 6-23 the numerical results of Al(t) and A2(t) are shown for Ay =5m, 
along with the theoretical asymptotic values from Gelhar and Axness [1983]. 
For l1f=O.5 the macrodispersivity appears to have reached a constant value 
while in the other two cases tlle macrodispersivity is still increasing. The 
longitudinal macrodispersivities for A=10m are shown in Figure 6-24, along 
with the theoretical asymptotic values. In none of the cases have the plumes 
reached asymptotic spreading rates. For the highly variable hydraulic 
conductivity fields (l1f=l.O and l1f=1.5) it appears that the plumes would 
need to travel much larger distance to achieve the predicted 
macrodispersivity. These results indicate that the ergodic hypothesis has not 
been satisfied for Ay =10m. 
The plume macrodispersivity is now analyzed in the transverse direction. 
The numerical transversemacrodispersivities, A2(t), for Ay =5m are shown in 
Figure 6-25. The experimental macrodispersivities are significantly larger 
than the predicted asymptotic values yet are consistent with numerical 
results from Tompson and Gelhar [1990]. By increasing the correlation length 
of .fuK(x) to 10m in Figure 6-26 the differences remain between the 
numerical and theoretical predictions. The numerical differences from, the 
stochastic theory may come from the observed hole-effect in the velocity 
correlation structure. 
6.5 Temporal Moments 
The third method for analyzing the results is by looking at temporal 
moments, or breakthrough curves. These curves help characterize the 
longitudinal solute movement by displaying the, arrival time of the solute at a 
barrier, or face. In this two-dimensional simulation the face is perpendicular 
to the mean flow direction. Two faces are chosen. Face 1 is located at 
xp=xl=80m and Face 2 is at xp=x1=200m (the right boundary). The mean 
first arrival time at Face 1, T(xl=80m,O%), is the average time it takes for 
the first particle to cross xl=80m for each of the thirty Monte Carlo 
realizations. The mean last arrival time at Face 2, T(x2=200m,99%), is the 
average time for 99% of the particles to leave the domain, for each of the 
realizations. Arrival times between these extreme times are also collected, 
65 
C!) 
C!) 
-e 
---~ 
-> 
-tf.) 
~ 
til 
Q.. 
rn 
a 
o 
0:: 
t.> 
< 
::2 
~ 
< z 
-o 
::> 
E-< 
-C,!) 
z 
:3 
11.0~----~----~----~----~--~----~----~----~----~----~ 
• a~=O.5 
10.0 
9.0 
• a~=1.0 
• a~= 1.5 
Gelhar and Axness[1983] 
Ay=5m 
8.0 
7.0 
• - --+ 
• 
6.0 • 
• 
• 
5.0 · ._--------------------~--- -----------. 
• 
4.0 • 
• • • • • 3.0 
• • • • 
----.--- -----------------. --------
2.0 • • • • 
• • 
1.0 
• 
o.o~----~----~----~----~----~~--~----~----~----~----~ 
o 20 
Fig .. 6-23 
40 60 80 100· 120 140 160 
" 
TIME (days) 
Longitudinal macrodispersivity and theoretical asymptotic 
values [Gelharand Axness, 1983], (Ay-5m, 0'f-0.5,1.0,1.5). 
180 200 
1------'-----
..-
E 
-~ 
-> 
-en 
0:: 
rzl 
tl. 
en 
-0 
0 
0:: 
t) 
< ::s 
....l 
< Q') z 
-.;r 
-0 
::> 
E-o 
-t!J 
Z 
0 
....l 
20.0 
18.0 
16.0 
14.0 
12.0 
10.0 
8.0 
6.0 
4.0 
2.0 
• 0"~=O.5 
• u~=1.0 
• a~= 1.5 
-Ay =10m 
Gelhar and Axness [1983] 
--------------------------------------------
u~=1.5 
u~=1.0 
.. .. 
.. 
.. 
.. e 
.. 
e 
e 
.. e 
u,=O.5 .. ____ e ______________________ • 
---- -----. . 
• 
.. 
• 
• 
.. 
• 
• 
• • 
• • 
• • 
O.O-rr----~----~----~----~----~--~----~----~----~----~ 
o 20 
Fig. 6-24 
40 60 80 100 120 140 160 
TIME (days) 
Longitudinal macrodispersivity and theoretical asymptotic 
values [Gelhar and Axness, 1983], (>.y-lOm, CTf=O.5,l.O,1.5). 
180 200 
~ 
cP 
1.0 
• a'=O.5 
0.9-1 • a~= 1.0 
Ay=5m 
..-... • a~= 1.5 5 0.8 Gelhar and Axness [1983] • • ~ • 
- 0.7 • • > • • 
-en 
~ • Prl 0.6 • • 0.. 
rn 
-Q 
0 0.51 • • ~ • • • • u • • < • • ::s 0.4 • 
Pel 
rn 
• ~ Pel 0.3 • > rn 
• • • • z • • • • • < 0.2 • • Ct! • E-
0.1 • <1'=0.5 a'=1.0 a'=1.5 
-
-------------------------------------------~ ~ ~0.0 .- - - - - - - - - - - - - - - -' - - - - '- - - - - - - - - - - - - - - - - - - - - - - - . 
I I 
0 20 
Fig. 6-25 
I I I I I I I I 
40 60 80 100 120 140 160 180 
TIME (days) 
Transver~e macrodispersivityand theoretical asymptotic 
values [Ge1har and Axness, 1983]; CAy-5m, O"f=O.5,1.0,1.5). 
200 
1.3 
1.2 
-E 1.1 
-i: 1.0 
-> 
- 0.9 CIJ 
0:: 
r:il 0.8 0.. 
CIJ 
-Q 0.7 0 
0:: 
u 0.6 < ::s 
Czl 0.5 
en CIJ co O:! 
Czl 0.4 
> CIJ 
z 0.3 < 0:: 
E-t 0.2 
ii' 
ii' 
0.1 'i! iI: 
Ii: 
ii: 0.0 iii 
I': 
Ii! 
0 
iii 
:jl 
Ii' 
!i 
Ii: 
iii 
:i: 
Iii 
I!, 
!I 
II 
,Ii 
I': 
II; 
'i' I, 
'I: 
il: 
• a~=O.5 
• a'=1.0 
• a~=1.5 Ay=10m 
Gelhar and Axness [1983] 
• 
• • 
• 
+ 
+ • + 
• • • 
• • • 
• 
• 
• 
• • • • • 
• •• 
• • • • 
• 
• • 
• 
•• ... __ ',_ _ _ _ _ _ _ _ _ _ _ _ _ a~=O.5 a'=1.0a~=1.5 
-------------------------
I: 
20 
Fig. 6-26 
40 60 80 100 120 140 160 180 
TIME (days) 
Transverse macrodispersivity and theoretical asymptotic 
values [Gelbar and Axness, 1983], (Ay-lOm, 0"t-0.5,l.O,1.5). 
200 
. __ .---L.J 
specifically for P=25%, 50%, and 75%. Beside the mean arrival times the 
standard deviations, ST(xl=80m,P%) and ST(x2=200m,P%), about the mean 
arrival times may be determined. The calculation of the mean and standard 
deviation, for thirty Monte Carlo realizations, is explained in (5-26) and 
(5-27), respectively. 
6.5.1 Mean Arrival Times 
The breakthrough curve through Face 1 for Ay=5m is shown in Figure 
6-27. The y-axis is the normalized concentration where c represents a 
number of particles and Co is the total number of particles injected into the 
domain (5000). The x-axis is the normalized time. All three degrees of 
heterogeneity are included. The first arrival time, T(xl=80m,0%), appears to 
slightly decrease by increasing O"y. Thus the front of the plume proceeds 
more quickly by increasing O"y. This concurs with the results of Smith and 
Schwartz [1980] and shows that the solute follows preferred paths of higher 
hydraulic conductivity. Further, the tails of the plume are substantially 
longer by increasing O"y, which indicates that a more heterogeneous porous 
medium contains the solute longer. Figure 6-28 is the breakthrough curves 
for Ay =lOm. As was found by Smith and Schwartz [1980], all arrival times 
increase by increasing the spatial continuity. The increase in arrival times is 
larger for larger O"y. Interestingly, for the larger correlation length of .fuK(x), 
the first arrival does not depend on O"y. 
Breakthrough curves through Face 2 are shown in Figure 6-29 and 
Figure 6-30. In Figure 6-29, for Ay =5, the differences in first arrival times 
are more pronounced at the right boundary. In a comparison between Figure 
6-27 and Figure 6-29 it appears that (Jy has a larger influence on first 
arrival the further it travels. This indicates that the plume continuously 
follows conduits of higher hydraulic conductivity. The tails of the plumes 
become more different as the transport distance increases. In Figure 6-30, 
for Ay=10m, all arrival times are larger than for Ay=5m. The differences in 
first arrival times begin to show for Face 2, whereas for Face 1 there had 
not been sufficient time for (Jy to influence the first arrival. 
70 
1.00 
---- a~=0.5 
----
a~= 1.0 
+---+ a~= 1.5 
-.. 
0 
0 
" 0.75 CJ 
--Z 
0 
...... 
e-
< n:: 
E-
Z 
~ 0.50 u 
Z 
0 
U 
Q 
~ ~ 
J--L N 
O.25J jJj Face 1 ...... ...:I < Ay =5m ::s 0:: 
0 
Z 
0.00 I I+'1 -I 
5 10 15 20 25 30 35 
NORMALIZED MEAN TIME (T-v lAy) 
Fig. 6-27 Mean arrival time through Face 1, CAy-5m, 0"..y-0.5,l.O,1.5). 
-;r 
~ 
C\ 
1.00 
.......... 
0 
C> 
" 
C) 0.75 
--z 
0 
-E-« 
< ~ 
e-
Z 
PrJ 
u 
z 
0 
u 
0 
PrJ 
N 
-.....l 
< :a 
~ 
0 
Z 
0.50 
O.25J 
0.00 Ii" i 
5 
ill 
Fig. 6-28 
10 
I /' /" I 
/ / ./ 
15 20 25 
NORMALIZED MEAN TIME (T'v lAy) 
Mean arrival time through Face 1, (Ay=10m, 
~~=O.5,1.0,1.5). 
-'--
I 
~ a~=O.5 
--- a~=1.0 
+--+ a~= 1.5 
Face 1 
Ay =10m 
\. 
30 35 
1II,i 
!: 
Ii 
I ~ 
,i: I,: 
Ii' 
" ii, 
(!: 
i" 
iii, 
i!'! 
~ i; 
-:t 
~ 
-
0 
0 
" 0 
--Z 
0 
-e-
< ~ 
E-o 
Z 
r:.J 
t.> 
Z 
0 
U 
Q 
r:.J 
N 
-....l 
< ~ 
0:: 
0 
Z 
1.00 _ LMi i:::>tO' 
0.75 
o.so 
0.25 
~ a~=O.5 
~ a~=1.0 
+-----+ a~= 1.5 
0.00 I ' (,-; ,"', 
10 15 20 25 30 35 40 
NORMALIZED MEAN TIM~ (r'v lAy) 
Face 2 
Ay=5m 
45 so 
Fig. 6-29 Mean arrival time through Face 2, (Ay-5m, 0'.y-O.5,l.O,1.5). 
55 
r , 
iii 
:\; 
Iii 
I! 
Ii 
I! 
II' 
!I 
I!: ,. i' il 
~ 
1.00 
, / '/ Y 
11-11 a'=0.5 
-.-. (12=1 0 Y .• 
-
/ / ./ +-+ 0"=1.5 
0 
CJ 
.......... 0.75 CJ 
........ 
Z 
0 
-~ 
O:! 
E-
Z 
rz:l 0.50 u 
z 
0 
t.) 
0 
rz:l I O.25J jjj Face 2 Ay =10m 
0 
z 
0.00 I \"'r-I 10 1 I 
Fig. 6-30 
NORMALIZED MEAN TIME (T'v lAy) 
Mean arrival time through Face 2, (Ay-10m, 
(T~O.5,1.O,1.5). 
55 
6.5.2 Standard Deviation of Arrival Times 
The standard deviation of arrival times gives an indication of the 
uncertainty in predicting arrival times. Figure 6~31 and Figure 6~32 show 
the standard deviations of arrival times at Face 1 for Ay=5m and Ay=10m, 
res.p.ectively. B.y in .. c.reasing the. degree of heteroge .. neity all t.imes of arrival 
become more variable. The tails of the plume are more uncertain than the 
fronts. These results match those of Smith and Schwartz [1980]. Figure 
6~31 indicates that by increasing the spatial continuity of the hydraulic 
conductivity, greater unpredictability is introduced. Figures 6-33 and 6-34 
are the standard deviations at Face 2. As the plume travels further, the 
uncertainty continues to increase. For example, in the case of Ay=lOm and 
(JY=1.5, the mean time for the entire plume to pass Face 1 is 330 days and 
the standard deviation is 180 days. At Face 2 the mean time for last arrival 
increases to 520 days and the standard deviation increases to 250 days. The 
significant variability in the character of breakthrough curves indicates that 
the prediction of concentration distributions is difficult even when the 
parameters which describe the random hydraulic conductivity field are known 
accurately. 
75 
I 
I, 
I, 
iii 
-.:r 
a:. 
1.00 
~ ------ a~='O.5 
~ a~=1.0 
-
~ a~=1.5 
0 
CJ 
"'-C) 0.75 
--z 
0 
-E-
< 0::: 
e-
z 
t:El 
to) 
z 
0 
u 
Cl 
r£l 
N . 
.-. 
...J 
< ~ 
~ 
0 
z 
0.50 
O.25J 1 j j Face 1 Ay =5m 
0.00 I I'"' "'i K 
a 1 2 
Fig. 6~31 
3 4 5 6 7 8 9 10 11 
NORMALIZED STANDARD DEVIATION (Sr'v lAy) 
Standard deviation of the mean arrival time through Face 1, 
(Ay=5m, CT-f=O.5,LO,1.5). 
12 13 
-
0 () 
" () ......... 
z 
0 
-t-o 
< 0:: 
E-
Z 
t':cJ 
U 
Z 
0 
u 
0 
t':cJ 
-.;r N 
-.;r ::l 
< :s 
0:: 
0 
Z 
·1'i
1
, 
r: 
JI 
Ili 
" 
Ii,!, 
I
,il 
\i 
1.00 _, ;w ::;:oott 
0.75 
0.50 
0.25 
Face 1 
Ay=10m 
+-+ ai=1.5 
---
0.00 -, ... l' 
o 1 2 
Fig. 6-32 
3 4 5 6 7 8 9 10 11 
NORMALIZED STANDARD DEVIATION (s,.." lAy) 
Standard deviation of the mean arrival time through Face 1, 
(Ay =10m, CT..y-0.5,l.O,1.5). 
12 13 
Ii, 
W 
W 
!:i ~ 
I'~ , I )! 
" 
~ 
I' Ii! 
I~ II 
I'i. :i . 
~ ~ I: 
I' i~ 
:~ 
'I ~ II~ ,: 
Iii 
l~ 
..;:( 
(p 
1.00.- .... 
-
o 
CJ 
i> 0.75 
---z 
o 
-~ 
0:: 
E-c 
Z 
Czl 
u 0.50 
z 
o 
U 
0: 
rxl 
N 
:3 
« 
~ 0.25 
o 
z 
Face 2 
Ay=5m 
+-+ a~=1.5 
--- a~=1.0 
0.00 I I 1- 1* N'I 
234 
Fig. 6-33 
5 6 7 8 9 10 11 12 13 14 15 16 17 1 B 19 20 
NORMALIZED STANDARD DEVIATION (s,.·v lAy) 
Standard deviation of the mean arrival time through Face 2, 
(Ay =5m1 u-f-O.5,l.O,1.5). 
~ 
<:0 
1.00 _ » ;;w 
-. 
0 
0 
.......... 
o 0.75 
~ j 
e-
< ~ 
e-
z 
tEl 0.50 u 
z 
0 
u 
Q 
tEl 
N 
-...l 
< ~ 0.25 ~ 
0 
z 
0.00 I -[ 1-
234 
Fig. 6-34 
I 
Face 2 / / Ay=10m 
/ 
+-+ a~=1.5 
--- a~=1.0 
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
NORMALIZ.ED STANDARD DEVI~TION (~'V/Ay) 
Standard deviation of the mean arrival time through Face 2; 
(Ar l0m1 (J.y-0.5 11.0 11.5). 
CHAPTER 7. CONCLUSIONS 
Numerical simulations were performed to investigate the effects that the 
random, spatially.-correlated, hydraulic conductivity field has upon large..-scale 
solute transport. The two-dimensional analysis involves conservative solute 
moving through a steady..-state, saturated flow region. The hydraulic 
conductivity fields are· assumed to be isotropic, weakly..-stationary, 
log-normally distributed, and the exponential distribution is used to describe 
the correlation structure of the natural logarithm of the hydraulic 
conductivity, inK. The random fields are parameterized by the mean Ily, 
variance O"f' and correlation length Ay. The mean of tnK is held constant, 
so that effects of the spatial continuity and degree of variability of tnK can 
be investigated. Thirty Monte Carlo realizations are used for each set of 
parameters (Ay=5,10m, and O"f=0.5,1.0,1.5, a~ shown in Table 6-1). 
The turning bands method reproduces the random, spatially.-correlated 
hydraulic conductivity fields accurately. The correlation length of each of the 
cases are the same as that assumed. From the isotropic K(x) fields, 
anisotropic velocity fields are produced. The longitudinal velocity has a 
correlation length in the longitudinal direction (A ) that is approximately 
vu 
two times Ay. The other three components of the velocity covariance 
function [pv (r2)' P (rl), Pv (r2)] are similar to each other, have correlation 1 V2 . 2 
lengths less than Ay, and display an apparent hole-effect. These results 
match analytical solutions. The correlation length of the concentration field 
in the longitudinal direction (Ac1) is about three times the hydraulic 
conductivity correlation length. In the transverse direction, the hydraulic 
conductivity, longitudinal and transverse velocities, and the concentration 
correlation lengths are all approximately· the same. The ratio between the 
longitudinal and transverse concentration correlation lengths is about three. 
Results from the spatial moment analysis show that the asymptotic 
mean plume velocity reaches the mean pore-water velocity, which is imposed 
by constant head boundary conditions. The plume spreading observed in this 
study departs from theory in several cases. For low Ay, the experimental 
results match theory very well for early times. However, after large travel 
distances. the theoretical predictions . underestimate. the observed sprea.ding, 
both in the longitudinal and transverse directions. This may be due to the 
hole-effect witnessed in the correlation structure of the velocity field. This 
80 
hole--effect behavior may have an impact on the transport phenomenon only 
after the plume has traveled a sufficient distance. Fickian conditions are only 
reached for uf:=O.5. For larger >'y the same pattern exists. Experimental 
results match theory in early stages, but the numerical spreading becomes 
larger at later times. Results from the temporal moment analySis indicate 
that the plume tailing is highly dependent on the correlation length and the 
variance of lnK(x). Even when these parameters are known there is great 
variability. 
Asymptotic plume spreading seems to be achieved only for the least 
heterogeneous case (>'y=5m and uf=O.5). When uf is increased above one, 
the asymptotic plume spreading is not reached. This has usually been 
associated with the fact that stochastic theories assume small perturbation, 
and are not valid for uf greater than one. However, this study suggests 
that in some cases, due to the highly anisotropic velocity field, ergodicity 
may not satisfied. 
Therefore, one of the results of this study is that the guidelines for the 
ergodic hypothesis may need to take into account the velocity correlation 
length instead of the hydraulic conductivity correlation length. In an attempt 
to satisfy ergodicity in this study, Monte Carlo realizations were performed 
with a relatively large domain. But the expected asymptotic behavior was 
not reached in cases with large >'y or large ufo The velocity correlation 
length, instead of the hydraulic conductivity correlation length, may be the 
parameter for which the ergodic hypothesis should be tested, because it is the 
velocity field that is the direct input for transport models. The size of the 
domain is decreased by more than two times when the velocity correlation is 
used and ergodicity may not be satisfied for this size domain. Independent 
realizations of the hydraulic conductivity field should be used to inspect the 
. effect of uf upon the velocity correlation structure. In this study, uf had 
no effect on the velocity correlation structure because the same realizations 
were used to produce the hydraulic conductivity fields. 
81 
     
REFERENCES 
Ababou, R., L.W. Gelhar, and D. McLaughlin. (1988). "Three-dimensional 
flow in random porous media," Technical Report Number 313, Volumes 
1 and 2, R.M. Parsons Laboratory, Department of Civil Engineering, 
M.LT., Cambridge, Massachusetts. 
Ahlstrom, S., H. Foote, R. Arnett, C. Cole, and R. Serne. (1977). "Multiple 
component mass transport model: Theory and numerical 
implementation," Rep. BNWL 2127, Battelle Pacific Northwest 
Laboratory, Richland, Washington. 
Bakr, A.A. (1976). "Effect of spatial variations of hydraulic conductivity on 
groundwater flow," Ph.D. dissertation, New Mexico Institute of Mining 
and Technology, Socorro, New Mexico 
Bakr, A.A., L.W. Gelhar, A.L. Gutjahr, and J.R. MacMillan. (1978). 
Stochastic analysis of spatial variability in subsurface flows. Part 1. 
Comparison of one- and three-dimensional flows," Water Resources 
Research, Vol. 14, pp. 263-275. 
Barry, D.A., J. Coves, and G. Sposito. (1988). "On the Dagan model of 
solute transport in groundwater: Application to the Borden site," Water 
Resources Research, Vol. 24, pp. 1805-1817. 
Bartlett, M.S. (1975). The Statistical Analysis of Spatial Pattern. Chapman 
·and Hall, London. 
Bear, J. (1972). Dynamics of Fluids in Porous Media. American Elsevier, 
New York. 
Bear, J. and A. Verruijt. (1987). Modeling Groundwater Flow and Pollution. 
D. Reidel, Dordrecht, Holland. 
Dagan, G. (1984). "Solute transport in heterogeneous porous formations," 
Journal of Fluid Mechanics, Vol. 145, pp. 151-177. 
Dagan, G. (1986). "Statistical theory of groundwater flow and transport: 
Pore to laboratory, laboratory to formation, and formation to regional 
scale," Water Resources Research, Vol. 22, pp. 120S-134S. 
Dagan, G. (1987). "Theory of solute transport by groundwater," Annual 
Review of Ffuid Mechanics, Vol. 19, pp. 183-215. 
Dagan, G. (1989). :Flow and Transport in 'Porous Formations. Springer-
Verlag, New York. 
83 
Dagan, G. (1990). "Transport in heterogeneous porous formations: Spatial 
moments, ergodicity, and effective dispersion," Water Resources 
Research, Vol. 26, pp. 1281-1290. 
de Marsily, G. (1986). Quantitative Hydrogeology. Academic Press, San 
Diego. 
Delhomme, J.P. (1979). "Spatial variability and uncertainty in groundwater 
flow parameters: A geostatistical approach," Water Resources Research, 
Vol. 15, pp. 269-280. 
Freeze, R.A. (1975). "A stochastic..-conceptual analysis of one..-dimensional 
groundwater flow in nonuniform homogeneous media," Water Resources 
Research, Vol. 11, pp. 725-741. 
Freeze, R.A. and J.A. Cherry. (1979). Groundwater. Prentice-Hall, Englewood 
Cliffs, New Jersey. 
Freyberg, D.L. (1986). "A natural gradient experiment on solute transport in 
a sand aquifer. Part 2. Spatial moments and the advection and 
dispersion of nonreactive tracers," Water Resources Research, Vol. 22, 
pp. 2031-2046. 
Gelhar, L.W., A.L. Gutjahr, and R.L. Naff. (1979). "Stochastic analysis 
of macrodispersion in a stratified aquifer," Water Resources Research, 
Vol. 15, pp. 1387-1397. 
Gelhar, L.W. and C.L. Axness. (1983). "Three-dimensional stochastic analysis 
of macrodispersion in aquifers," Water Resources Research, Vol. 19, pp. 
161-180. 
Gelhar, L.W. (1986). "Stochastic subsurface hydrology from theory to 
applications," Water Resources Research, Vol. 22, pp. 135S-145S. 
Goode, D.J. (1990). "Particle velocity interpolation in block-eentered 
finite difference groundwater flow models," Water Resources Research, 
Vol. 26, pp. 925-940. 
Graham, W. and D. McLaughlin. (1989). "Stochastic analysis of nonstationary 
subsurface solute transport. Part 1. Unconditidnal moments," Water 
Resources Research, Vol. 25, pp. 215 ... 232. 
Hoeksema, R.J. and P.K. Kitani.dis .. (1985). "Analysis of the spatial structure 
of properties of selected aquifers," Water Resources Research, Vol. 21, 
pp. 563-572. 
Hromadka,T.V., T.J. Durbin, and J.j. Devries. (1985). Computer Methods in 
Water Resources. Lighthouse Publications, Mission Viejo, California. 
84 
Huyakorn, P.S. and G.F. Pinder. (1983). ComputAtional Methods in 
Subsurface Flow. Academic Press, San Diego. 
Journel, A.G. and Ch.J. Huijbregts. (1978). Mining. Geostatistics. Academic 
Press, San Diego. . 
Kinzelbach, W. (1988). "The random walk method in pollutant transport 
simulations," Groundwater Flow and Quality Modeling. E. Custodio, A. 
Gurgui, and J.P.t. Ferreira (Eds.). D. Reidel, Hingham, Massachusetts. 
Mackay, D.M., D.t. Freyberg, P.V. Roberts, and J.A. Cherry. (1986). "A 
natural gradient experiment on solute transport in a sand aquifer. Part 
1. Approach and· overview of plume movement," Water Resources 
Research, Vol. 22, pp. 2017-2029. 
Mantoglou, A. and J.t. Wilson. (1982). "The turning bands method for the 
simulation of random fields using line generation by a spectral method," 
Water Resources Research, Vol. 18, pp. 1379-1394. 
Matheron, G. (1973). "The intrinsic random functions and their applications," 
Advanced Applied Probability, Vol. 5, pp. 438--468. 
Matheron, G. and G. de Marsily. (1980)- "Is transport in porous media 
always diffusive? A counterexample,' Water Resources Research, Vol. 
16, pp. 901-917. 
Meija, J., and 1. Rodriguez-Iturbe. (1974). "On the synthesis of random fields 
from the spectrum: An application to the generation of hydrologic 
spatial processes," Water Resources Research, Vol. 10, pp. 705-711. 
Moltyaner, G.L. and RoW.D. Killey. (1988). "Twin Lake tracer tests: 
Transverse dispersion," Water Resources Research, Vol. 24, pp. 
1613-1627. 
Moltyaner, G.L. and RoW.D. Killey. (1988). "Twin Lake tracer tests: 
Longitudinal dispersion," Water Resources Research, Vol. 24, pp. 
1628-1637. 
Neuman, S.P. (1982). "Statistical characterization of aquifer heterogeneities: 
An overview," Special Paper Number 189, Geological SOCiety of 
America, pp. 81-102. 
Neuman, S.P., C.L. Winter, and C.M. Newman. (1987). "Stochastic theory of 
field-scale Ficldan dispersion in anisotropic porous media," Water 
Resources Research, Vol. 23, pp. 453--466. 
Pinder, G.F. (1973). "AGalerkin~finite element simulation of groundwater 
contamination of Long Island, New York," Water Resources Research, 
Vol. 9, pp. 1657-1669. 
85 
Pinder, G.F. and W. Gray. (1977). Finite. Element Simulation in Surface 
and Subsurface Hydrology. Academic Press, Orlando, Florida. 
Prickett, T.A., T.G. Naymik, and C.G. Lonnquist. (1981). "A random walk 
solute transport model for selected groundwater quality evaluations," 
Illinois State Water Survey Bulletin 65, Champaign, Illinois. 
Rubin, Y. (1990). "Stochastic modeling of macrodispersion in heterogeneous 
porous medla," Water Resources Research, Vol. 26, pp. 133~141. 
Salandin, P. and A. Rinaldo. (1990~. "Numerical experiments on dispersion 
in heterogeneous porous media,' Computational Methods in . Subsurface 
Hydrology, pp. 495~501. 
Shapiro, A.M. and V.D. Cvetkovic. (1988). "Stochastic analysis of solute 
arrival time in heterogeneous porous media," Water Resources Research, 
Vol. 24, pp. 1711~1718. 
Smith, L. and RA. Freeze. (1979). "Stochastic analysis of steady state 
groundwater flow in a bounded domain. Part 2. Two~imensional 
simulations," Water Resources. Research, Vol. 15, pp. 1543-1559. 
Smith, L. and F.W. Schwartz. (1980). "Mass transport. Part 1. A stochastic 
analysis of macroscopic dispersion," Water Resources Research, Vol. 16, 
pp. 303-313. 
Smith,L. and F.W. Schwartz. (1981). "Mass transport. Part 2. Analysis of 
uncertainty in predication," WaterResources Research, Vol. 17, pp. 
351-369. 
Sudicky, E.A. (1986). "A natural gradient experiment on solute transport 
in a sand aquifer: Spatial variability of hydraulic conductivity and its 
role in the dispersion process," Water Resources Research, Vol. 22, pp. 
2069-2082. 
Tompson, A.F.B., R Ababou, and L.W. Gelhar. (1987). "Application and 
use of the three~imensional turning bands random field generator: 
Single realization problems," Technical Report Number 313, R.M. 
Parsons Laboratory, Department of Civil Engineering, M.LT., 
Cambridge, Massachusetts. 
Tompson, A.F.R, E.G. Vomvoris, and L.W. Gelhar. (1988). "Numerical 
simulation of solute transport in randomly heterogeneous porous media: 
Motivation, model development, and application," Technical Report 
Number 316, RM. Parsons Laboratory, Department of Civil Engineering, 
M.LT., Cambridge, Massachusetts. 
86 
1 
I 
I 
~ 
1~-~-~-- . -.--.-.-. ___ -.-
Tompson, A.F.B., R. Aba.bou, and L.W .. Gelhar. (1989). "Implementation of 
the three-dimensional turning bands random field generator," Water 
Resources Research, Vol. 25, pp. 2227-2243. 
Tompson, A.F.B. (1990). "Turn2d: Program for generating a two..-dimensional, 
stationary, correlated anisotropic random field," Personal 
Communication. 
Tompson, A.F.B. and L.W. Gelhar. (1990). "Numerical simulation of solute 
transport in three-dimensional, randomly heterogeneous porous media," 
Water Resources Research, Vol. 26, pp. 2541-2562. 
Uffink, G.J.M. (1988). "Modeling of solute transport with the random walk 
method," Groundwater Flow and Quality Modeling. E. Custodio, A. 
Gurgui, and J.P.L. Ferreira (Eds.). D. Reidel, Hingham, Massachusetts. 
87 
     
l 
I 
I 
J 
APPENDIX A. HYDRAULIC CONDUCTIVITY INTERPOLATION 
The hydraulic conductivity field produced from the turning bands 
method assigns values to nodes. Input for the flow simulations requires a 
value within a triangle, A slight variation on the usual interpolation scheme 
is used in this study, A two-dimensional example is given, 
For a triangular element with equal sides (Figure A-1) the element 
mean is usually the geometric mean of the three nodes 
lis 
KE :;:: (K, ,.K'+l ,.Ki '+1) I,J 1 ,J ,J 
However, if perpendicular bisectors, of lines connecting the nodes are drawn it 
is seen that node i,j should have twice as much influence than the other two 
nodes, To account for this KE should be skewed to account for this 
influence 
KE :;:: (K, ,·K, ',.K'+l ,.K,. '+1)1/ 4 I,J I,J 1 ,J I,J (A-2) 
Computations using this adjusted geometric mean appear to smooth out the 
velocity transition from element to element, Overall the effect on flow 
computations is expected to be small. 
89 
J 
J 
Fig. A-I 
~I ~~_~ ~ _~~~"'~o._.~~~~,~~~,~"-- _"_ 
Two-dimensional example for calculating 
K for flow simulations. ' --
90