PB87-141479
STOCHASTIC PREDICTION OF DISPERSIVE CONTAMINANT
TRANSPORT
Department of Civil Engineering
Cambridge, MA
Dec 86
EJBD ^^=^^^E^^^^=^==E=^^^^^^^=
ARCHIVE
EPA
600-
2-
114 U.S. DEPARTMENT OF COMMERCE
National Technical Information Service
-------
P887-141479
EPA/600/2-86/114
December 1986
STOCHASTIC PREDICTION OF DISPERSIVE CONTAMINANT TRANSPORT
by
Efstratios G. Vomvoris and Lynn W. Gelhar
Massachusetts Institute of Technology
Cambridge, Massachusetts 02139
f*
3"
— Cooperative Agreement
4 CR-81 1135-01-2
Project Officer
Joseph F. Keely
Applications and Assistance Branch
Robert S. Kerr Environmental Research Laboratory
Ada, Oklahoma 74820
ROBERT S. KERR ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
ADA, OK 74820
REPRODUCED BY
U S. DEPARTMENT OF COMMERCE
NATIONAL TECHNICAL
INFORMATION SERVICE
SPRINGFIELD. VA 22161
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing/
REPORT NO.
EPA/600/2-86/114
3 RECIPIENT'S ACCESSION NO, 1-
-' * -- *• '
TITLE AND SUBTITLE
STOCHASTIC PREDICTION OF DISPERSIVE CONTAMINANT
TRANSPORT
5 REPORT DATE
December 1986
6. PERFORMING ORGANIZATION CODE
AUTHORIS)
Efstratios G. Vomvoris and Lynn W. Gelhar
8 PERFORMING ORGANIZATION REPORT NO
PERFORMING ORGANIZATION NAME AND ADDRESS
Massachusetts Institute of Technology
Department of Civil Engineering
Cambridge, Massachusetts 02139
10 PROGRAM ELEMENT NO
CBPC1A
11 CONTRACT/GRANT NO
CR-811135
2 SPONSORING AGENCY NAME AND ADDRESS
Robert S. Kerr Environmental Research Lab. - Ada, OK
Office of Research and Development
U.S. Environmental Protection Agency
Ada, Oklahoma 74820
13. TYPE OF REPORT AND PERIOD COVERED
Final (09/01/83 - 04/30/86)
14 SPONSORING AGENCY CODE
EPA/600/15
5 SUPPLEMENTARY NOTES
6 ABSTRACT
The objective of this research agreement was to develop mathematical models
to quantify the concentration variability observed in field measurements of
concentration plumes. The concentration variability is attributed mainly to
spatial heterogeneity of the hydraulic conductivity field. Since there is limited
information about actual distributions of hydraulic conductivity in a given site,
the log-hydraulic conductivity is modeled as a three-dimensional anisotropic
stationary random process. It is shown that the concentration variance, which
measures the intensity of the variations, is proportional to the mean concentra-
tion gradient and to the variance and correlation scales of the log-hydraulic
conductivity; it is inversely proportional to the local dispersivity values. The
main contribution to the concentration variability results from the coefficients
corresponding to the longitudinal mean concentration gradients. Analysis of a
portion of the data from the Canadian Forces Base, Borden, Ontario, shows the
applicability of the developed results in field situations. Simple analytical
examples demonstrate the way to use the results in a predictive context.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.lDENTIFIERS/OPEN ENDED TERMS
COSATi field/Croup
18. DISTRIBUTION STATEMENT
RELEASE TO PUBLIC.
19 SECURITY CLASS iTIut Reporti
UNCLASSIFIED
21 NO 01
104
20 SECURITY CLASS (Tinspage/
22 PRICE
UNCLASSIFIED
EPA Form 2220-1 (R.». 4-77) PREVIOUS EDITION is OBSOLETE
-------
DISCLAIMER
"The information in this document has been funded wholly
or in part by the United States Environmental Protection
Agency under Cooperative Agreement CR-811135-01-2 to the
Massachusetts Institute of Technology. It has been subjected
- to the Agency's peer and administrative review, and it has
been approved for publication as an EPA document."
Mention of trade names or commercial products does not constitute
endorsement or recommendation for use.
ii
-------
FOREWORD
EPA is charged by Congress to protect the Nation's land, air and water
systems. Under a mandate of national environmental laws focused on air and
water quality, solid waste management and the control of toxic substances,
pesticides, noise and radiation, the Agency strives to formulate and imple-
ment actions which lead to a compatible balance between human activities and
the ability of natural systems to support and nurture life.
The Robert S. Kerr Environmental Research Laboratory is the Agency's
center of expertise for investigation of the soil and subsurface environment.
Personnel at the Laboratory are responsible for management of research pro-
grams to: (a) determine the fate, transport and transformation rates of
pollutants'in the soil, the unsaturated and the saturated zones of the
subsurface environment; (b) define the processes to be used in character-
izing the soil and subsurface environment as a receptor of pollutants; (c)
develop techniques for predicting the effect of pollutants on ground water,
soil, and indigenous organisms; and (d) define and demonstrate the applica-
bility and limitations of using natural processes, indigenous to the soil
and subsurface environment, for the protection of this resource.
This report contributes to that knowledge which is essential in order
for EPA to establish and enforce pollution control standards, predict the
extent of existing contamination, and to develop ways to restore the
environment, when contaminated, in a reasonable and cost-effective manner.
Clinton W. Hall
Director
Robert S. Kerr Environmental
Research Laboratory
ill
-------
CONTENTS
Abstract vi
Figures vii
Tables ix
Symbols x
Section
1. Introduction 1
2. Conclusions 3
3. Recommendations 4
4. Natural Variability and the Stochastic Approach 5
5. Analysis of Concentration Variability 15
6. Application Examples 38
7. Field Data Analysis: The CFB Borden Site 47
8. Discussion of Results 67
References °^
Appendices
A. The Mathematics of the Approximate Analytical Evaluations ... 74
B. The Integration Over Wave-Number 84
C. Transformation to Polar Coordinates 88
D. Mean Concentrations Fitted at the Borden Site Data 91
Preceding page blank
-------
ABSTRACT
The objective of this research agreement was to develop mathematical
models in order to quantify the concentration variability observed in field
measurements of concentration plumes.
The concentration variability is attributed mainly to the spatial
heterogeneity of the hydraulic conductivity field. Since there is limited
information about the actual distribution of the hydraulic conductivity in a
given site, the log-hydraulic conductivity is modeled as a three-dimensional
anisotropic stationary random process. It is characterized by its mean and
its covariance function or spectrum. Implementing Darcy's equation and the
convective-dispersive equation at the local scale, relationships among the
concentration variations and characteristic parameters of the porous medium
are found. Approximate analytical expressions are also developed.
It is shown that the concentration variance, which measures the intensity
of the variations, is proportional to the mean concentration gradient and to
the variance and'correlation-scales of the log-hydraulic conductivity; it. is
inversely proportional., to the local dispersivity values. The main
contribution to the c6ncentration variability results from the coefficients
corresponding to the longitudinal mean concentration gradients.
There is a very high correlation in the concentration values along the
direction of the flow, which is still under investigation. The orientation
of the hydraulic conductivity principal axes does not have a significant
effect on the concentration variance values.
Analysis of a portion of the data from the CFB Borden site shows the
applicability of the developed results in field situations. Simple analytical
examples demonstrate the way to use the results in a predictive context.
This report was submitted in fulfillment of Cooperative Agreement No.
CR-811135-01-2 by the Massachusetts Institute of Technology under the
sponsorship of the U.S. Environmental Protection Agency. This report covers
a period from September 1, 1983 to April 30, 1986, and work was completed as
of April 8, 1986.
vi
-------
FIGURES
Page
(a) Hydraulic conductivity data from laboratory cores
from a borehole in the Mt. Simon aquifer Illinois
(Bakr, 1976) 6
(b) Hydraulic conductivity profile from a core in the
CFB Borden site, Canada (Sudicky, 1985b) 6
4.2 Concentration measurements from the Borden site.
(a) Depth averaged data (b) Vertical profiles
(Freyberg et al. 1983 and personal communication) .... 8
4.3 Schematic graph of a stationary one dimensional
covariance function of fcnK showing the correlation
scale \ 10
4.4 Schematic graph and comparison of the measured,
predicted with field scale dispersion and predicted
with laboratory scale dispersion concentration profiles . 10
5.1 Coordinate systems for the case X2 = X2*; the
dashed elipse represents schematically the e"1
level of the correlation function 13
5.2 Comparison of sample autocorrelation function (crosses)
for hydraulic conductivity (from grain-size analysis)
with theoretical autocorrelation function (solid line)
(From the CRNL, Ontario experiment, Moltyaner, 1985) ... 24
5.3 Comparison of Che exponential correlation function and
the "hole-type" correlation function 25
5.4 Comparison of the variance for an anisotropic aquifer
(correlation scales l^, 12* ^ and an equivalent
isotropic (correlation scale i. = (41^2^3) ) 29
5.5 Comparison of the approximate analytical expressions
and the numerical results for the TJJ coefficients .... 30
(a) TH(.) along the x direction for different e
(b) T22(«) along the x direction
(c) T22(«) along the y direction
5.6 Test of the validity of the approximate result for
the variance TH (1, 1, e, 1;0) with respect to the
e values 32
6.1 Normalized mean concentration and standard deviation at
x = 50 ra, z = 0 ra 39
6.2 Isoconcentration lines and standard deviations at z = 0 m 41
vii
-------
Page
Correlation function for a point at y = z = 0 m 42
(a) Correlation along x
(b) Correlation along lines in the x-y plane at
various angles with x
(c) Lag-correlation contours in the x-y plane
6.4 Normalized concentration and one standard deviation
interval along a longitudinal section (line coordinates
y = z = 0) 44
6.5 Normalized concentration and one standard deviation
interval along a vertical section (line coordinates
x = xc - center of mass of plume, y = 0) 45
6.6 • Normalized concentration and two standard deviation
interval along another vertical section (line
coordinates x = xc - 5 m, z = 1 ra) 46
7.1 Instrumentation, location of cores and hydraulic head
distribution on August 10, 1984 (solid line) and November
15, 1984 ( ) at the tracer-test site (after Sudicky
1985b). The circled samplers show the locations of the
wells used in this analysis. V denotes the well used for
the vertical section analysis 48
7.2 Isoconcentration contours at depth z = -5.31 ra from the
Borden site data (original data from 0. Freyberg, personal
communication) 49
7.3 Schematic presentation of the methods used for the
estimation of the mean concentrations 58
7.4 The measured, fitted and plus-minus one standard deviation
concentration profiles along the longitudinal section at
z = -5.31 m 59
7.5 The mean concentration curves fitted at the wells on the
longitudinal section (c in the vertical vs. z) 60
7.6 The measured, fitted and plus-minus one standard deviation
concentration profiles along the vertical section at (x,y)
= (38.3, 16.08) 64
7.7 The values plotted in Fig. 7.6 shown in a semi-logarthmic
plot 65
C.I Schematic graph of the coordinate transformation geometry . . 90
viii
-------
TABLES
Page
Approximate analytical results for the
concentration covariance 27
5.2 Effect of the stratification orientation on the
concentration variance 34
7.1 Existing Information on the Borden site parameters .... 48
7.2 Effect of the input parameters on effective
hydraulic conductivity 53
7.3 Effect of the input parameters on the
macrodispersiviti.es 55
7.4 Parameters used in the variance estimation 56
7.5 Concentration gradients in the vertical section
at (x,y) = (38.3m., 16.08m.) 56
7.6 Estimation of the range of variations in the
longitudinal section » 62
7.7 Estimation of the range of variations in the vertical
section " 63
D.I Parameters used to fit mean concentrations at the wells . 91
ix
-------
LIST OF SYMBOLS
A auxiliary variable defined in Eq. (5.27)
All» A22» A33 macrodispersivities along axes x, y, and z
B auxiliary variable defined in Eq. (5.27)
c concentration of the transported solute (volumetric)
dZfOO Fourier-Stieltjes wave amplitude of the random process f(x_)
E[ • ] expectation operator
E.. local bulk dispersion coefficient tensor
E coefficient defined in Eq. (A.I)
F mean of log-hydraulic conductivity
f(x) perturbation of log-hydraulic conductivity field
and generic symbol for a random field in Section 4.
G mean concentration gradient along x.; G. = .
j J J J
K.. coefficients used in the evaluation of K. . ;
1J defined in Eq. (7.11) 1J
J hydraulic gradient in direction x.
K hydraulic conductivity
1C defined from exp[E(lnK)j
1C effective hydraulic conductivities
K(t) auxiliary variable defined in Eq. (5.16)
k^ the three-dimensional wave number vector; (kj , k2 , k3)
i characteristic correlation scales of the hole type
i exponential covariance function
-------
m(x) mean of a random field
m rate of mass injection
M total mass injected
n porosity
q specific discharge in the x direction
R_,(x,, x_) covariance function of the random field (f(x_)
Sf,(2c. , £9) spectrum of the random field (f(:c)
T " nondimensional coefficients used in the evaluation of
JP the concentration covariance; defined in Eq. (5.33)
i
T. dimensional coefficients used in the evaluation of
JP the concentration covariance; defined in Eq. (5.33)
V mean pore velocity
£ (xj, x2, x3) = (x, y, z) coordinate vector and
coordinate axes in the three-dimensional space
x location of center of mass of a plume
cu , cu, longitudinal and transverse local dispersivities
6 coefficient defined in Eq. (A.I)
Y flow factor = q/K J.
f^il' 2i2^ semi-variogram of a random process
ij Kronecker delta, equals 1 if i = j
e defined from a,/i^
z defined in Eq. (A.I)
6 auxiliary variable defined in Eq. (5.16)
K ratio of transverse to longitudinal local dispersivity
X correlation scale of the negative exponential
covariance function
xi
-------
PI , p2 aspect ratios of the covariance of the £nK
defined in Eq. (5.33c)
p (xj, Xj) correlation function of the random process g(x)
SB
2
o the variance of a random field g(x)
g ~
1(1 rotation angle of the concentration plume with
respect to the mean flow direction
(overbar) mean of a random field
' (prime) perturbation of a random field
xil
-------
SECTION 1
INTRODUCTION
Field observations of solute transport in groundwater indicate that
dispersivities are scale dependent and can range over many orders of
magnitude. This discrepancy between the laboratory and field experience has
been attributed mainly to the spatial variation of the flow properties of the
natural porous earth materials. Therefore, over the last decade a substantial
portion of the research in the area of subsurface hydrology has focused on
understanding the effects of the spatial heterogeneity of natural formations
on flow and -solute transport.
Continuously expanding computer capabilities have facilitated the use of
deterministic numerical models with parameters which may be spatially variable
over a grid and have, in principle, made it possible to obtain solutions that
previously could be described only qualitatively. However, such solutions are
limited to the specific application and are unable to generate relationships
valid over a broad range of problems. More importantly, in order for a
numerical model to reproduce the actual detail of a given field situation, it
must have sufficient grid resolution to represent the velocity field producing
the mass transport and information about the controlling parameters, such as
hydraulic conductivity and porosity, at each grid point. While the former may
be practical for small idealized flow systems (say over tens of meters), the
number of'nodes required to resolve the actual three-dimensional heterogeneity
at a more practical field scale (say 1 km) can easily exceed the capacity of
even the largest supercomputers. Furthermore, the actual measurement of the
required parameters with that degree of resolution is a totally impractical
task.
It is the lack of information about the spatial distribution of the
controlling parameters that has led to modeling them as realizations of
random fields. A random field can be viewed as a generalization of the random
variable concept. Before measuring the value of a random variable the only
information known is a range for its possible values (defined by its mean and
variance). In a random field it is further assumed that there is a
correlation in space among the different values, so that knowledge of the
value at one location carries information about the possible values for the
hydraulic conductivity at different locations in such a way that they fall
into the prespecified ranges (defined by the mean and the variance) and have
the appropriate correlation in space (defined by the correlation function).
This approach is coined the stochastic approach. The fundamental assumption
is that for each one realization the classical laws of flow and mass transport
hold at some representative scale.
-------
In Che stochastic approach, the mean or expected value behavior of the
concentration field is found to be equivalent to the classical dispersive
transport equation. The theory also predicts the resulting field-scale
dispersion coefficients in terms of the statistical characteristics of the
hydraulic conductivity variability. However, the mean concentration
represents the average behavior of the ensemble rather than the actual
realization of the aquifer. Therefore, a classical transport model will
simulate a smoothed approximation of an actual field simulation; the actual
concentration will vary around the smooth mean. In order to realistically
present the results of model simulations, it is necessary to understand and
quantify this error in classical models.
The stochastic approach can be used to quantify the variation around the
mean concentration. The main focus of this research is the development of the
theory to describe such variability. The degree of variability will be
characterized by the concentration variance. It will be shown that the
concentration variance depends on measurable or determinable quantities
characterizing the aquifer material, such as the correlation scales and the
variance of the log-hydraulic conductivity, the inclination of the stratifi-
cation, and the local dispersivity values. It also depends on the mean
concentration gradient which, through the macrodispersivities, represents the
overall effect of the aquifer heterogeneity on the mass transport. Section 6
shows how the mean concentration gradient can be estimated in a predictive
context with numerical or analytical models. Section 7 demonstrates how it
can be estimated through a simple curve fitting with existing plume
measurements, by the analysis of the field data from the Borden site in
Canada.
The stochastic approach is complementary to classical deterministic
modeling. The predicted concentration variance can be used as a realistic
and physically-based calibration target for the large-scale numerical models.
Furthermore, one can use any classical numerical or analytical model that
reproduces the general characteristics of the contamination event and, in
addition, obtain a measure of the expected variability around the model
results. Finally, knowledge of the possible deviations of the actual plume
from the predicted one can improve substantially the design of monitoring
schemes.
The present study is focusing on quantifying the variation of the
concentration measurements or model predictions, in heterogeneous porous
formations. Examples demonstrate the practical application of these results.
-------
SECTION 2
CONCLUSIONS
The concentration variability resulting from the heterogeneity of the
natural aquifers can be quantified using stochastic analysis in terms of
measurable or determinable properties of the hydraulic conductivity field. It
can be expressed mathematically as the sum of the products of the concentration
gradients weighted by appropriate coefficients, which depend on the variance and
the correlation scales of the log-hydraulic conductivity, the local dispersivity
values and the orientation of the stratification. The dominant coefficient
corresponds to the gradient along the mean flow and is inversely proportional to
the local dispersivity and proportional to the correlation scales.
The effects of the orientation of the stratification in a preliminary
analysis appear to be less important, but further investigation is required for
definitive conclusions. The structure of the concentration covariance for the
particular spectrum used results in a singular effect, namely, high correlations
along the mean flow direction.
Approximate analytical results have been developed that establish the role
of the key parameters.
The preliminary analysis of field data from the Borden site shows that the
results of the proposed theory can be applied in practical situations and
produce consistent confidence intervals for the anticipated concentration
variability.
-------
SECTION 3
RECOMMENDATIONS
A local relationship that quantifies the concentration variability
anticipated in mass transport in heterogeneous aquifers has been found. The
relationship involves the gradient of the mean concentration, the correlation
scales of the log hydraulic conductivity field and the local dispersivity
values. As shown in Section 7, it can be readily used and produce quite tight
confidence intervals for the concentration values, even with crude estimation
of the mean concentration field.
The obtained concentration standard deviation can be implemented for the
calibration of numerical models. The discretization inherent in the numerical
models induces a smoothing or averaging of the true concentration field and
makes it impossible to reproduce the measured values. The estimated standard
deviation of the concentration can be used as a measure of the acceptable
deviation of the numerical model results from the measured values.
The preliminary application of the results to a field site was very
encouraging. Further analysis is needed to explore the full three-dimensional
nature of the controlled plume and to follow its movement and its response to
the orientation of the stratification. Important issues such as spatial
averaging of the concentration field introduced by the sampling technique and
its associated filtering properties need to be addressed.
The developed theory established the basic relationships among the aquifer
characteristics, the flow parameters and the concentration variance. However,
the effects of unsteady concentration input, initial conditions, flow reversal,
and unknown history of injection need to be investigated. The singular behavior
of the concentration covariance along the mean flow direction requires
additional investigation before accepted as a valid descriptor of the physical
system.
Finally the issue of the mean concentration field, central to the theme of
all stochastic theories, should be investigated. Methods to enhance the ability
to identify or estimate the mean concentration need to be developed.
-------
SECTION 4
NATURAL VARIABILITY AND THE STOCHASTIC APPROACH
THE THEORETICAL FRAMEWORK
The conventional approach In groundwater modeling is Co treat an aquifer
as an idealized homogeneous formation with a constant value for each physical
parameter so that the idealized aquifer behaves like the actual formation. In
cases where there is predictable heterogeneity one might divide the aquifer
into a number of homogeneous zones, each with a different equivalent value for
any used parameter. The simple case of known anisotropy can be treated
incorporating suitable coordinate transformations.
However,- the actual heterogeneity in a formation which would
conventionally be treated as homogeneous can be very pronounced. Figure 4.1
illustrates this variability for the hydraulic conductivity which is
considered the most important flow property of an earth porous material.
Note that the quantity plotted is the logarithm of the hydraulic
conductivity. The observed spatial variation is well documented in the
literature; "... the jigs and jags of every electric log, the spread of
measurements from every sampling program, and common sense geological
reasoning attest to this fact." (Freeze, 1975). The intensity of the
variability of hydraulic conductivity can be characterized by the
estimated variance of the natural logarithm of the hydraulic conductivity
and which may be as high as 13 for some rocks (Freeze, 1975).
In addition, due to the underlying depositional process, the hydraulic
conductivity field will typically exhibit some spatial structure. Oftentimes
this structure manifests itself as layers or lenses and reflects directional
dependence or anisotropy.
In a practical sense variability of the porous formation results in
variations of the flow field and the head distributions around the ones
predicted for homogeneous systems. It is further propagated to the
distribution of the concentration since the latter is directly affected by the
flow field. Figure 4.2 presents actual concentration data from the field
scale experiment at the Borden site. Notice that intense variation of the
measured values is observed for the vertical sections and an even substantial
one for the depth averaged data, Figure 4.2a.
The fundamental question that emerges is how can these effects be
modeled. If it is assumed that a complete knowledge of the hydraulic
conductivity process is available, the existing and continuously expanding
computer capabilities that permit the implementation of numerical models with
parameters spatially variable across the grid, can produce the required
solutions. However, at the current level of computing capabilities the
allowed resolution of the input field is not fine enough to capture the
anticipated effects. Furthermore, such solutions will be case specific and
unable to establish general relationships among the key parameters valid over
a broad range of problems. But the "Achille's heel" of the sophisticated
numerical modeling of the spatial heterogeneity is that for all practical
-------
100
DISTANCE, FEET
200
213.238
213.0-18 —
217.048
IU
IU
10 IU
Hydraulic Conductivity (cm/a)
Figure 4.1 (a) Hydraulic conductivity data from laboratory cores from
a borehole in the Mt. Simon aquifer, Illinois (Bakr,
1976)
(b) Hydraulic conductivity profile from a core in the CFB
Borden site, Canada (Sudicky, 1985b)
-------
purposes it is impossible to acquire all the detailed data required to adequately
represent the actual variation of the parameters.
The lack of information about the controlling parameters, in addition to the
need for analytical and general relationships, has led to modeling the
appropriate parameters of the porous medium as realizations of random fields.
The fundamental assumption is that for each one realization, and usually we have
only one in groundwater related problems, the classical laws of flow and mass
transport hold at some representative scale.
The concept of a random field or spatial stochastic process is a
generalization of a random variable. If £ represents a point in the three-
dimensional space and f(x) a random variable corresponding to that point jt, the
set of random variables ff(Xi)» f (.£2 ) i ••• f(.2k)» •••} where 2^,3^, ... Xfc»
are the coordinate vectors of points 1,2, ... k, is called a random field and is
simply written as f(jc).
As a random variable is completely characterized by its probability density
function, for the complete probabilistic characterization of a random field f(j>0
one would need to determine the joint probability density function of any set
fCxj), f(x2>, ... f(xji) for anv n and for i » 22 • ••• 2Sn arbitrary points in
the spatial domain of interest. In most applications the entire spatial law can
neither be identified from finite data nor is required; the first two moments of
the distribution suffice to characterize the random fields in a number of
problems. The first two moments are,
- the mean function or first-order moment of the random field
E[f(x)] = m (x) . (4.1)
where E[«] denotes mathematical expectation
- the covariance function or second-order moment,
The zero lag covariance, when x- x~, is the variance, denoted as af2-
Another representation of a second-moment is the variogram (Journel and
Huijbregts, 1978) introduced in the field of geostatistics, and defined as,
Note that Eq. (4.1) and Eq. (4.2) is a generalization of the concepts of mean
and variance of a random variable. Eq. (4.2) incorporates the information about
the spatial structure of the process. It defines the correlation between two
points of the field located at x and x , if it is normalized by the
1 2
variance.
The mean and the covariance can in general depend on the location x, > as can
be seen in Eqs. (4.1) and (4.2). When the mean and covariance are independent
of jc, the process is said to be (second order) stationary or statistically
homogeneous:
-------
X
(m)
I !
1C
9
6
•7
6
5
A
S
i
1
-6
DAY 85
DAY 0
9 IOO 200 "JOG 444 600
C (mg/l)
0 100 ISO 330 iOO ICO
C (mg/l)
y (m)
-ll—
r
o 103
C (mg/l)
Figure 4.2 Concentration measurements from the Borden site, (a)
Depth averaged data, (b) Vertical profiles (Freyberg et
al, 1983 and personal communication)
8
-------
E{f(x)} = m
Rff(V V = Rff(xl" V (4"4)
If the field is said Co be isotropic, the same correlation exists for points
that are separated by the same distance regardless of their orientation in
space, i.e.,
RffQi^x^ = Rff(|-i~-2r (4"5)
Figure 4.3 illustrates graphically a one-dimensional stationary covariance
function. One characteristic length scale is the separation distance at which
the correlation decreases to the e level, designated as the correlation
scale. In the general stationary case the covariance can always be expressed in
terms of the spectral density function or spectrum as follows,
where the spectral density function can be expressed as,
Sff 00 = I 3 /// e ~^'L R
(2ir) -«
and k_ is the three-dimensional wave number vector. The Fourier transform pair
in Eqs. (4.6) and (4.7) implies that knowledge of only one of the two, spectrum
or covariance, is sufficient. They both carry the same information, but it is
displayed differently. The covariance function has a physical meaning that is
easier to grasp: how two actual values of the process under consideration, a
certain distance apart "are expected to be correlated. The meaning of the
spectrum is more abstract: if one could decompose a whole set of statistically
similar processes in trigonometric waves and make a histogram of the encountered
wave-length, that would be the spectrum. It gives the distribution of the
variance over the wave-number domain.
If a process is statistically homogeneous, it will always have a spectrum,
and then it is possible to take advantage of the representation theorem
[Luraley and Panofsky, (1962), Gelhar (1984)],
CO
f(x) = /// e1^ dZf(k)
E{dZf(k) dZ* (k')} = Sff(k)dk , if k = k1
= 0 , if k_# k_' (4.8)
where dZ,(k) is the Fourler-Stieltjes wave amplitude of the process f(:c).
This theorem essentially states that the process f(x) is composed of a sum of
uncorrelated or orthogonal increments in the wave number domain. Furthermore,
the representation in Eq. (4.8) is unique and provides a convenient formalism
for treating differential equations involving stationary variables.
9
-------
s x,-x2
Figure 4.3 Schematic graph of a stationary one dimensional
covariance function of £nK showing the correlation
scale X
fitted (laboratory scale
** dispersion)
fitted (field scale
dispersion)
Figure 4.4 Schematic graph and comparison of the measured,
predicted with field scale dispersion and predicted
with laboratory scale dispersion concentration profiles
10
-------
The spectral methodology will be adopted In the current work. The_
characteristic random fields are decomposed in a known mean function, g, and a
zero-mean stationary perturbation, g1,
g(x) = g(x) + g'U) (4.9)
Utilizing the representation theorem in the differential equations
relating the various fields, one can relate the spectrum of the dependent
variable to the spectrum of the input processes.
Following the traditional deterministic approach, the aquifer is modeled
as a continuum and the accepted expressions for the mass conservation, the
solute conservation and Darcy's equation are assumed to hold at the local
scale.
REVIEW OF RELATED WORK
Groundwater Flow
Historically the problem treated initially in a probabilistic framework
was that of groundwater flow. The approaches developed can be briefly
summarized into the following categories:
Monte Carlo Simulations—
Realizations of sets of aquifer properties are generated through a model,
and the direct problem is deterministically solved for each realization. From
the statistical analysis' of these solutions, stochastic characteristics of the
head and flow fields can be inferred. Freeze (1975) in his pioneering paper
treated the problem of one-dimensional saturated flow in finite domain where
the head at the boundaries was specified. Spatial variation but not spatial
structure was the main limitation. The latter was introduced later in the
work of Smith and Freeze (1979a, 1979b), where the nearest neighbor model
generator was used to generate one- and two-dimensional fields.
Perturbation or Linearization Models—
This was the most widely used theoretical approach. The stochastic
differential equations are solved using "small" pertubations theory or Taylor
expansions in the space or frequency domain. Gelhar (1984) presents a good
review of the spectral methodology which has been used in a series of papers,
namely, Gelhar (1974, 1977), Bakr et al. (1978), Gutjar et al. (1978), and
Mizell et al. (1980). In the aforementioned works head variances and
covariances are obtained for the cases of one-, two- and three-dimensional
flow in confined and unconflned aquifers. Comparable results are obtained in
Gutjahr and Gelhar (1981), Dagan (1982a), Vomvoris (1982), Kitanidis and
Vomvoris (1983), Hoeksema and ICitanidis (1984), and McLaughlin (1985) using
first order analysis in the space domain. The advantage of the latter
approach is that it can be applied in bounded domains, while the spectral
method is tailored for infinite stationary domains. Note that all the
results, regardless of the method used, show that the one- and two-dimensional
cases provide erroneously large variances, due to the unrealistically limited
number of paths that are offered to the flow system.
11
-------
Composite Media Analysis—
This approach was followed by Dagan (1979) and follows closely the
self-consistent method described in Beran (1968). The hydraulic conductivity
field was represented as a collection of homogeneous blocks with random radii
arranged at random locations. A similar method, that of the embedding matrix
was presented in Chirlin and Dagan (1980).
Solute Transport in Groundwater
In the late seventies, the probabilistic framework was extended to include
the solute transport problem where the heterogeneity of the natural formations
is believed to be more influential. An apparent influence is on the
dispersivity values corresponding to field scale plumes where an increase of
two or three orders of magnitude with respect to the anticipated laboratory
measured values was observed. A critical analysis of existing field
observations of solute transport has revealed that dispersivities are scale
dependent, Gelhar et al. (1985). Furthermore, concentration measurements
exhibit significant fluctuations around the values predicted from analytical
or numerical solutions. The review by Anderson (1979) suggests that the
dispersivities used in numerical models are little more than fitting
parameters and as such lack a sound physical basis.
The central questions that the related research has focused on are:
- What are the governing equations for large-scale mass transport in a
heterogeneous aquifer?
- How do the different scales affect Che characteristic parameters of what
seems to be the same phenomenon?
A typical situation is depicted in Figure 4.A. The two smooth curves
correspond to numerical or analytical solutions predicted with laboratory
scale dispersivities and. fitted field scale dispersivities. The actual
measurements exhibit significant fluctuations around the values predicted by
the models. The current work therefore attempts to answer the following
question:
- How can the concentration variability around the predicted solutions be
quantified?
The first attempt to show the effects of spatial heterogeneity was through
Monte Carlo simulations [Schwartz (1977), and Smith and Schwartz (1980)].
From a numerical analysis point of view, the solution of the mass transport
equation is by far more difficult than the solution of the flow equation.
Instabilities in the algorithm are very frequent and numerical dispersion of
the same order of magnitude as the physical dispersion may occur (Ababou
et al. 1985). No quantitative conclusions were drawn from the Monte Carlo
simulations.
Gelhar et al (1979) quantified for the first time the relationships
among field scale or macro-dispersivities and characteristic measures of the
hydraulic conductivity field. The physical situation examined was that of
macrodispersion in a stratified aquifer, and the method used was spectral
analysis and Fourier-Stieltjes representation. The authors developed the
governing equations near and far away from the source, demonstrated the
transition in the dispersivity until it reached its asymptotic value and
showed the skewness of the concentration distribution near the injection vs.
12
-------
the symmetric Gaussian concentration distribution at points far downstream.
Some further aspects of the stratified case were discussed in Matheron and
de Marsily (1980). In a subsequent paper, Gelhar and Axness (1983) treated
the general three-dimensional case of raacrodispersion in aquifers. That work
focused on the asymptotic dispersivity coefficient from the simple case of
isotropic media to the most general one of statistically anisotropic media
with arbitrary orientation of the stratification. The results produced
numerical values of the expected order of magnitude and compared favorably
with the scarce existing field data.
Dagan (1984) produced similar results for the raacrodispersivity following
a Lagrangian approach which considers displacement of particles in random
velocity field generated from a random hydraulic conductivity field. This
analysis for the case of a statistically isotropic medium yields results for
the asymptotic dispersion coefficient which are equivalent to those of Gelhar
and Axness (1983).
A detailed quantitative testing of the macrodispersivity prediction has
not been completed yet [Gelhar, 1985], but a number of carefully designed
field experiments are now underway [Freyberg et al. (1983), LeBlanc (1985),
and Betson et al. (1985)]. Hufschmied (1985) has shown that
macrodispersivities calculated from measurements of the heterogeneity of
hydraulic conductivity field are consistent with those estimated for a
cold-water plume in an aquifer in Switzerland. Encouraging comparisons are
also shown in Sudicky (1985b).
The issue of concentration variation has not been yet widely recognized.
Gelhar et al. (1981) developed the concentration variance for the stratified
case in answering a comment on their 1979 paper. They found that the
concentration variance is proportional to the variance of the log hydraulic
conductivity, the square of the mean concentration gradient, the fourth power
of the correlation scale of the hydraulic conductivity and inversely
proportional to the square of the local transverse dispersivity. The
assumption of perfect layering is very restrictive and thus these results are
primarily of qualitative significance in demonstrating the factors which
influence the concentration variability. These results do indicate that
significant variations of concentration are to be expected across the layers
even after very large mean transport distances.
Gelhar and Gutjahr (1982) have developed the equation for the
concentration variance for one-dimensional transport of a solute undergoing
radioactive decay. They have considered transport in a stream tube with
varying cross-section, thus incorporating the influence of a three-dimensional
flow variation. The results relate the concentration variance to statistical
parameters of the longitudinal velocity field. Dependence of the
concentration variance on the square of the mean concentration gradient and
the square of the correlation scale of the velocity is observed. Due to the
one-dimensional structure of the problem, analytical solutions for the
concentration perturbations can be obtained which could be used in a finite
domain. The radioactive decay results in a finite concentration variance even
when local mixing is eliminated by setting the local dispersivity value equal
to zero.
13
-------
Gutjahr et al. (1985) developed similar results following a slightly
different methodology. They treated the same problem with emphasis on the
development of equivalent macrodispersivity values. Their result for the
concentration variance is still under investigation.
Dagan (1982b, 1984) is the only other available work focusing on the
calculation of the concentration variance. His approach follows closely the
theory of diffusion by continuous movements [Taylor (1921)] referred to above
as the Lagrangian method. In a qualitative sense, Dagan's results are similar
to the ones found in this research work. One of the assumptions made was that
the effects of local disperslvity can be neglected when calculating the mean
motion of the plume. It is thus believed that Dagan's results are applicable
in a region near the source of pollution where the convective component of the
hydrodynamic dispersion is predominant. In contrast, the results in this
study correspond to the far field downstream region. From the analyses of the
mean concentration, it is indicated that the local disperslvity does not have
an important effect on the macrodispersivity values. The opposite effect is
observed when examining the concentration variance. A final point is that the
approach of Dagan is performed for statistically homogeneous, isotropic media,
while the current work is based on the more realistic case of anisotropic
hydraulic conductivity fields.
14
-------
SECTION 5
ANALYSIS OF CONCENTRATION VARIABILITY
THE GOVERNING EQUATIONS
The general equation describing transport of an ideal nonreactive
conservative solute in a saturated, rigid porous medium is, [Gelhar and Axness
(1983)]:
3° j. „ 3c - 1 v l£ (5.1)
n IT+ qi^ Tx. Eij 3x. <5'1J
where, n = porosity
q. = specific discharge in the *^ direction
c = concentration of the transported solute
E .= local bulk dispersion coefficient tensor.
The dispersivity values in Eij will be assumed constant at the scale at
which Eq. (5.1) is valid. Eq. (5.1) incorporates the conservation of mass
equation for a homogeneous fluid,
1 q = 0 (5.2)
3xt Mi
The fields of log hydraulic conductivity, specific discharge and
concentration will be decomposed into a mean and a small perturbation about
the mean. Therefore,
F + f ' (log hydraulic conductivity)
q = q.+ q' i = 1, 2, 3 (specific discharge) (5.3)
c = "c + c1 (concentration)
The log hydraulic conductivity perturbation field, f, is assumed to be a
statistically homogeneous random field with zero mean and known covariance
function. Note that through the Darcy equation, the specific discharge and
the hydraulic conductivity are directly related, [Gelhar and Axness (1983)].
The concentration perturbation, c', becomes a random field also and
substituting Eq^_ (5.3) into Eq. (5.1) the equations for the mean
concentration, c, aiid perturbation, c1, can be found. Without loss of
generality, we can align the Xj axis along the mean flow, which renders
q - q\ * 0 , ?2 = ?3 = 0 (5.4)
The local dispersion tensor may then be approximated in the form
"aLq 0 0
[E
°Tq
15
-------
where cu and OT are the local longitudinal and transverse
disperstvitiesT Eq. (5.4) and Eq. (5.5) can simplify the equations
for the mean concentration and perturbation to:
n |I + q |I + q. i£i = q OL a!| + q ttT (i2! + i!|j (5.6)
1 23
— 222
3c' . , 3c . 3c' 3 c' ,3 c' 3 c' ,, ,s n
n IT + qi 3^ + q JZ ~ q •L —2 ' q aT (—2 + ~~2 J (5<7)
, 3c' , 3c' -
= q! -2— - q' -7— =0
The zero mean term on the right hand side of Eq. (5.7) is of second order in
perturbations and is assumed to be negligible compared to the first order
terras on the left side.
In Gelhar and Axness (1983) Eq. (5.6) and Eq. (5.7) are used to calculate
an expression for the macrodispersivity tensor A . It is shown that far
away from the source, the mean concentration satisfies the classical transport
equation with constant macrodispersivities. In the current work, the
attention focuses on Eq. (5.7). It will be assumed that the mean
concentration field is known.
The effects of the velocity field variations on the local dispersion were
addressed in Naff (1978)'. He retained the actual velocity field in Eq. (5.5)
and found that the effects of the velocity variations on the macrodisper-
sivities were of second order. An extension of his analysis to study these
effects on the concentration variability showed that they were negligible in
this case also. Consequently, Eq. (5.5) is considered a valid representation
of the local dispersion effects, at the current level of approximation.
SPECTRAL SOLUTIONS
The equation for the concentration perturbation, Eq. (5.7), can be
easily solved if it is assumed that the concentration perturbation field is
statistically homogeneous. Invoking the Fourier-Stieltjes representation
theorem [Lumley and Panofsky (1962)] one can find relations between the
spectra of the concentration and log hydraulic conductivity. The simple
case of a steady concentration field will be investigated first.
Steady State Concentration Field
The perturbed quantities are represented as,
00 00
c' = / exp(ik.x)dZc(k); q^ = / exp(ik.x)dZq (k) (5.8)
—o» —oo i
where k=(k ,k ,k.) is the wave number vector, xfU^x^x.^) is the position
16
-------
vector, and the integration is over three-dimensional wave number space.
Combining eq. (5.8) and eq. (5.7) and invoking the uniqueness of the spectral
representation
n o o
=G,dZ_ (5.9)
where G.= - —=—, the mean concentration gradient, is taken to be locally
constant. j
*
The complex conjugate Fourier amplitude of dZ ft, namely dZ , can be easily
found from Eq. (5.9). Taking expectation of dZ d2 results inCthe following
spectral relationship:
S (k) = i -= s-J-i , 0 , (5.10)
where S (_k) = the three-dimensional spectrum of the concentration
cc perturbation,
= the three-dimensional spectrum of the specific discharges
along the axes x.and x (j, i = 1,2,3).
The final step is to establish the relationship between Sq.q (k) and
S,.(jOi the spectrum of the log-hydraulic conductivity field. JIn the
geneTal case, one expects aquifers to exhibit three-dimensional anisotropic
heterogeneity with the mean flow arbitrarily oriented with respect to the
stratification. The coordinate systems might look then as in Figure (5.1) where
it can be observed that 'the mean hydraulic gradient vector does not coincide
with the direction of the mean flow. The covariance function has its principal
axes at the primed system. Using the Darcy equation for a locally isotropic
medium, Gelhar and Axness (1983) showed that the specific discharge spectra can
be expressed in terms of the log-conductivity spectrum as
(5.1 la)
where j,£,m,n = 1,2,3. (Note that summation over m and n is implied.)
A simplified case corresponds to statistically isotropic £nK field, for
which the specific discharge spectrum becomes,
(k) - K J J (6 -
where y = 4/K Ji» and Ji ls the hydraulic gradient in the x direction.
The covariance function of the concentration field can be computed combining
Eq. (5.10) and Eq. (5.11),
17
-------
MEAN FLOW
DIRECTION
Figure 5.1 Coordinate systems for Che case *2 = X2'> cne
dashed elipse represents schematically the e"1
level of the correlation function
18
-------
where R (£) is the covariance function of the concentration field, at lag
£. The variance is found by taking e = 0,
00
Rcc(0) = / Scc(k)dk (5.13)
Unsteady Concentration Field
The concentration perturbation field, c'Ot.t), is, in general, time
varying. Again it will be assumed that the concentration field is locally
stationary (statistically homogeneous) since the time-invariant hydraulic
conductivity field is spatially homogeneous. This assumption is expected to
be particularly valid away from the source. Because the unsteady mean
concent rat io.n solution will be expressed in terms of a moving coordinate
system, 5 ( g = x - qt/n, ?2 = X2' 53 = X3^' 1C is aPPr°Priate Co express
the perturbation transport equation (Eq. (5.7)) in that system as
Jf 'I + "ill = v1 + «^— +-) (5-7a>
3t 'l 1 3Ci L 8«J
where c'(£,t), i.e., _£ and c and are the independent variables. The flow
field is steady in the fixed coordinate system so that q' in (5.7a) is
represented by
q!(x) = f e1-*- dZ (k) (5.14a)
1 ~ i» qi ~
but the time varying c'(^,t) will be represented as
CO
c'(x.t) = / ei-'-dZc (k;t,x.)
= / exp [i(k1?1 + qt/n) + k^ + k^\ dZc(k;t,ci)
— CO
= c'(?i,t) (5.14b)
Note that dZc is considered to vary slowly in space so that the dependence
on gi is only parametric. It should be mentioned parenthetically that the
introduction of the moving coordinate system is the more natural system to
use, but, of course, from a mathematical point of view the two systems are
equivalent. The moving system is introduced in order to take advantage of the
fact that in the moving system, the concentration field will vary slowly in
time.
19
-------
Using Eqs. (5.11a) and (5.14b) in (5.7a) along with the uniqueness of the
representation theorem,
n | I dZ (_k; t)- G.dZ QO + iqk dZ (k_;t)
dt|_£C J °.j 1C
2 22
= ~q lar'ci "*" *]•' 2 "*" 3 ^J "^c(_»^' (5.15)
Denoting,
dZ (k; t) = p
C q , ., A . 2 _,_ ,.2 _,_ . 2 ., ,, .,.
Q ^ ™" I T 1C "T" /v iC TrtCk.TiC.il I T L n J
n I 1 L, 1 ^ 2 3 '
K(t) = - G.dZq.(k)
Eq. (5.15) can be simply written as
|£ I + 0P = K (5.17)
With initial condition c'Ot.O) = 0, or in the wave number domain dZ (k,0) = 0
= p(0), Eq. (5.17) can be solved. The solution is,
-flt C Or
p = e 0C / K(T) e9T dT (5.18)
o
p=e"0CidZ (k) / e9T G.(£,T) dT (5.19)
" qj -- o J
The integral in this last equation can not be evaluated analytically (G.'s
are complicated functions which are unknown in the sense that they depend on the
solution of the concentration perturbation). Because of the exponential term in
Eq. (5.19), the main contribution will, for large 9, come from the region near
the upper limit, T = t. Therefore, G.(£.;T) is expanded around the value at 7
= t. The expansion leads to,
8Gi ,
G.(£,T) = G.(C,T) + r-1- _r (x-t) + . . . (5.20)
J J o T I T L
where it is implied that the time derivative is in the moving coordinate system
— I . This expansion also has the important advantage that it allows the
mean^transport equation to be approximated in the form of the classical
advection-dispersion equation; otherwise the governing equation will be known in
the form of an integro-differential equation wherein the macrodispersive flux is
dependent on the mean concentration gradient at earlier times. This feature is
discussed by Gelhar (1986). Here because the goal is to develop results which
can make use of mean solutions based on the classical transport equation, this
analytical approximation will be explored in detail. Keeping only second-order
derivatives and combining Eqs. (5.19) and (5.20):
20
-------
q. o =c
A change of variables, u=t—r, and the limit for t •*•<», gives.
p = IdZq / e-0u [G.Or.t) - ul |(5.0] du (5.22)
Consequently, dZc(k^,t) has an approximate analytical expression. From this
point and on the procedure of the previous section can be followed. Integrating
Eq. (5.22) the complex amplitude is given from,
dZc(k,t) - dZq [G.^.t) -j- 2 ] (5.23)
J 59
with complex conjugate,
dZ^ (k.t) - I dZ*_ [Cj(i.t) \ - „ ] (5.24)
The spectrum of c'(x,t) is found by multiplying the complex conjugate amplitudes
and taking expectations. Recall,
E [dZ (k;t) dZ* (k;t)l = Srr(k,t) dk;
c - c - cc - - (5>25)
E [dZ (k.) dZ* (k)] = S (k)d]c
qj qz qjqa
This operation leads to the following expression
i i i "*Gi i l 3Go l l
•«<*.<> - i ^,,1 v, ? ? - i^ e. ? ? - n1 BJ e- 1
ac ac Q2 Q*2
where summation over all the values of j and i is implied. To simplify notation
further denote
21
-------
.
n 1
B = J [aLk2 + aT(k2 + k2)] - (5.27)
0 = B + iA
9* = B - iA
08* = B2 + A2 (5.28)
222
0* = B - A + 2iAB
*2 2 2
0 i = B - A - 21AB
Thus, Eq. (5.26) can be written as,
i2Sq,q.Trh2 1GJG*+S
3G. , 3G
- " <5-29'
'<" -' " " n2 V* A2 * B2 '"1".
The relationship between Sq.q and Sff is presented in the previous
section (Eq. 5.1 la)). Notejtnat the inclusion of the time variation has
introduced three additional terras in the concentration spectrum, as can be seen
by the comparison oE Eq. (5.24) with Eq. (5.10). This implies that the "zero"
order behavior of the unsteady case solution will be determined from the steady
state. The results will be different only in the region where the concentration
gradients change rapidly, i.e., near the peak concentration.
THE LOG-HYDRAULIC CONDUCTIVITY SPECTRUM
The intuitive notion for a covariance function structure dictates that it
should have its maximum at zero lag and decay, in some way, as the lag increases
dropping to zero at some distance. A simple model that captures such an image
is a negative exponential. It has been used extensively and appears to be
consistent with the available data [(Bakr (1976) Sudicky (1985a)]. However
within the statistical limitation of field data, other forms of the covariance
can be equally acceptable. Occasionally, the mathematical structure of the
models used to describe the physical behavior require a stronger assumption
about the covariance behavior. Bakr et al. (1978), Gutjahr and Gelhar (1981),
Vomvoris (1982) have used different covariance models, which still show a decay
with distance but, in one sense, reinforce the periodic or stationary nature of
22
-------
the modeled process. The simplest of these latter models is the so-called "hole
type" covariance function (Bakr et al, 1978). Its main characteristic, besides
the negative exponential decay, is the negative correlation after a certain
lag. In the wave number domain this characteristic corresponds to a zero value
of the spectrum at jt = 0, i.e., zero integral scale. Moltyaner (1985, p.
363-22, 23) has observed such negative correlations in field data (Figure 5.2).
In turbulence, where similar mathematical analysis has been used, a popular
spectrum is the following [Hinze (1975), Aldama (1985)].
S(k) - k4exp (- ja2k2) (5.30)
which also has zero integral scale.
For the concentration covariance a similar zero integral scale log hydraulic
conductivity covariance has to be used. It can be seen from Eq. (5.10) that for
k2 = k3 =0, the kj in the denominator will cause the integral to blow up
at kj =0, when one integrates from -» to +». This behavior can be avoided if
Sq.q (0,0,0) = 0, or, from Eq. (5.11), Sff(0_) = 0, which is the
zero-integral requirement postulated.
A simple hole-type, three-dimensional covariance is the following:
2
2 ,, I N -s
= ff (1- s)e ,
2 1 (5.31)
s = —
2 2
42 V2*3 *!*!
a£fV!y - a z 3 (5.32)
•n (1 + ki^i^
The £i's are some characteristic scales of the hole exponential covariance
that can be uniquely related to the e"1 correlation scale. The notation \
will be reserved for the negative exponential covariance function correlation
scale.
2
Figure (5.3) compares the hole type with the exponential, R(^) = a* exp
(-£/X), when they have the same e~^ value for the same distance 5. The two
covariance functions are practically indistinguishable.
23
-------
RUTOCQBRELRTICN FUNCTION
~
!\J
..L
...L
| I J I I I I I I ! I I I I I I
LflC OISTSNCi (CM1
Figure 5.2 Comparison of sample autocorrelation function (crosses)
for hydraulic conductivity (from grain-size analysis)
with theoretical autocorrelation function (solid line).
(From the CRNL, Ontario experiment, Moltyaner, 1985)
24
-------
= a-(s/3s0))exp(-(s/se»
3.0
4.0
5.0
lag (g/x)
8.0
7.0
9.0
Figure 5.3 Comparison of the exponential correlation function and
the hole exponential correlation function
25
-------
APPROXIMATE ANALYTICAL ESTIMATIONS
The covariance function of the concentration field involves the evaluation
of the three-dimensional integral in Eq. (5.12), which, for the general case,
can not be performed in an analytical fashion. However, for some simple
situations, approximate analytical expressions can be found. These expressions,
although not directly applicable in all practical situations, reveal the
structure of the concentration covariance and the influence of the various
parameters.
The covariance function in Eq. (5.12) can be written as:
Rcc (£ } = Tjp,e,l > GjGp (5.33a)
where, the Einstein's notation for repeated indices is used and
K? f. • S
Tjp< WVW*'9'^ = ~2 / exp '-^ ,2
J* q -« k
' 1 Vn < > <* (5'33b)
where the Sff (jc) has the form of Eq. (5.32) and j, p, m, n can take the
values 1,2,3.
One can further convert the integral in Eq. (5.33b) to a function of non-
dimensional variables (Appendix B) so that Eq. (5.33a) becomes
Tjp(p1.P2'e''e'
'•<**.
(i) The concentration variance is inversely proportional to the local
dispersivity value. This is a physically intuitive behavior since the
smaller the local dispersivity the less the capacity of the medium to
attenuate developing tongues of fast moving solute. Gelhar et al.
(1979) and Matheron and de Marsily (1980) have described similar
phenomena .
26
-------
Table 5.L
Case
APPROXIMATE ANALYTICAL SOLUTIONS
T (E=0)
(e=0)
1.
isocropLC i
2
°C
inis. .
2
Y
2
3
J_
KE
1
2
Y
1
9
22
2. Anlsotroplc AnK
1 2 1 \_ \_
e PI P2
2 3
Y
1 I r 1 , _L_
2 6 L 2 + _2
Y 2Pt R
22
R - 2R
3. Isotroplc £nK
a. Rcc(e,0,0)
JL_ 2_ I
2 3 KE
l-B-B2eBEU-B)]
2313
b.
1112
2 . KG 3 ?
-16]
(0,0,5)
»,£,0)
T22(e,0,0)
T2.,(0,?,0)
< = — , R - 1 - p. ,
°L l
I '
El(-x) = Exponential Integral, K2(O = Modified Bessel function of 2nd order, *p2 = 1.
-------
(ii) The concentration variance is proportional to the product of the
correlation scales of the hydraulic conductivities. Higher correlation
scales imply that solute trapped in a more permeable layer tends to stay
there longer, since the formation persists for larger scales• Consequently
a randomly located measuring probe would have higher chances for a "hit or
miss" output.
Vomvoris and Gelhar (1985) presented a comparison of the variances of
two different media, an anisotropic and an isotropic one, having the same
product of correlation scales (Figure 5.4). The favorable comparison of
the two results indicated the option, for variance calculations, of
replacing an anisotropic medium by an equivalent statistically isotropic
one with correlation scale equal to the geometric mean of the anisotropic
correlation scales.
(iii) The concentration covariance along the direction of the mean flow has very
high correlation scales (Figure 5.5). This singular behavior is still
under investigation.
THE NUMERICAL CALCULATION
The equation for the concentration covariance, Eq. (5.12), cannot be
evaluated analytically for a general combination of the input parameters. A
direct three-dimensional numerical evaluation, unless very carefully designed,
can produce unreliable results due to the irregular behavior of the integrand.
A reduction to a two-dimensional integral can be accomplished if one converts
Eq. (5.12) to spherical coordinates and integrates over wave number u. Appendix
B presents the details of the above integration. The result is, Eq. (B.7) and
Eq. (B.U),
2* ir
= / /
w=o g=o
du dg sine F (f5,u)
t
{[A
n
A
13 3
(3+3
)] exp[-
J I
PIB!
exp -
where the different quantities:
M2>
(5.34)
A2, A,E,R,B are functions of
and are defined in Appendix B. The integral in Eq. (5.34) is evaluated first
over angle 0 using the Romberg integration technique and then with a simple
trapezoidal rule over u>«
28
-------
8.
0)
c
o
.2 S.
§- *°
£> SJ
o -:
»-_-«
-e
^
\
v
s — -e
o.i
0.8 1.6 2.4 3.2 4.0
f\j
p
O).
10
ct.'o 0.8
1.5
Figure 5.4 Comparison of Che variance for an anisotropic aquifer
(correlation scales £1, &2, £3) and an equivalent
isotropic (correlation scale i = («,,, £, , Jl
•nW-. ± £,
--._
29
-------
z
o
UJ
o
u
CM
-------
-------
so
c
sj
CD
g_
s.
OJ
d-4.0
2 -2.4
log (a/
approximatQ
analytical
numGricai
-1.6
-0.3
-n
0.0
Figure 5.6 Test of Che validity of the approximate result for the
variance Tu(l,l,e ,1;0) with respect to the e values
32
-------
The approximate analytical expressions in Table 5.1 are evaluated by
comparison with the numerical integration. A comparison of the numerical
results and the approximate analytical is shown in Figure 5.5. The agreement is
very good for small values of eic = of/i- Figure (5.6) compares the expression
for the variance. The analytical approximation gives identical results with the
numerical for values of eic = aj/i UP to 0.10.
Table 5.2 shows the effects of the angles of orientation of the
stratification on the variance. For the chosen correlation scales the effect of
the vertical inclination seems to be stronger. In most cases the values remain
of the same order of magnitude. A notable case is Tj^, the coefficient that
multiplies (-3c/3x1)(-3c/8x3). Its magnitude increases up to an order of
magnitude less than the value of TH« Since the longitudinal gradients are
very much milder than the vertical ones, such an increase could equalize the
effects of the longitudinal and vertical gradients in the estimation of the
variance. Table 5.2 demonstrates the complicated interaction of the correlation
scales, the-local dispersion coefficient and the orientation of the
stratification. Further investigation of these relationships is required.
Figure 6.3, in the next section, shows the results of the numerical
integration for the covariance of an isotropic hydraulic conductivity field.
The presented curves are directly proportional to T^. Very high correlation
scales along the flow direction are observed.
33
-------
Table 5.2
EFFECT OF THE STRATIFICATION ORIENTATION ON
THE CONCENTRATION VARIANCE
*
0
0
0
0
5
5
5
5
10
10
10
10
15
15
15
15
e
= =5 = = =S=
0
• 5
10
15
0
5
10
15
0
5
10
15
0
5
10
15
(a) Pl
Tll
518.5
592.96
653.14
716.22
547.74
598.41
656.95
723.73
583.84
612.45
658.47
725.5
607.27
625.75
668.91
726.47
= 0.02 ,
T22
0.72
1.69
4.57
8.83
0.77
1.7
4.48
8.73
0.89
1.77
4.32
8.34
1.09
1.88
4.22
7.87
P2 = 0.2 ,
T33
0.76
3.51
11.59
24.08
0.81
3.62
11.82
24.57
0.92
3.77
11.87
24.44
1.07
3.79
11.75
23.68
e = 0.1
T12
0
0
0
0
-2.27
-1.92
-1.81
-1.41
-5.83
-5.24
-4.19
-2.5
-10.07
-9.35
-7.65
-4.
, K = 1
0
30.05
64.58
102.0
0
30.87
65.05
102.06
0
32.01
65.79
101.97
0
32.47
65.99
101.31
T23
======
0
0
0
0
0
-O.Ob
-0.07
-0.01
0
-0.27
-0.26
-0.05
0
-0.47
-0.62
-0.23
34
-------
Table 5.2 (continued)
*
0
0
0
0
5
5
5
5
10
10
10
10
15
15
15
15
e
0
5
- 10
15
0
5
10
15
0
5
10
15
0
5
10
15
(b) p
TI
228
241
264
291
230
242
264
291
236
246
266
292
242
250
268
294
!~- 0.
1
.45
.45
.17
.39
.72
.84
.75
.74
.25
.19
.09
.59
.9
.97
.94
.08
05 ,
T22
0.
0.
2.
3.
0.
0.
2.
3.
0.
0.
1.
3.
0.
0.
1.
3.
P2 =
52
9
01
73
54
9
7
57
92
97
6
62
95
93
46
0.2
T33
0.
1.
4.
9.
0.
1.
4.
9.
0.
1.
4.
9.
0.
1.
4.
8.
i £
61
66
68
34
62
67
67
3
66
68
64
18
71
71
57
94
= 0.
Tl
0
0
0
0
-1
-0
-0
-0
-2
-2
-I
-0
-3
-3
-2
-1
1 , K
2
.07
.97
.74
.34
.29
.04
.57
.73
.49
.23
.50
.21
= 1
T13
0
12.4
25.62 •
i
39.89
0
12.44
25.63
39. 8J
0
12.51
25.62
39.62
0
12.53
25.51
39.2
T23
0
0
0
0
0
-0.05
-0.06
-0.01
0
-0.11
-0.14
-0.04
0
-0.17
-0.22
-0.09
35
-------
Table 5.2 (continued)
*
0
0
0
0
5
5
5
5
10
10
10
10
15
15
15
15
e
0
5
10
15
0
5
10
15
0
5
10
15
0
5
10
15
(c) Pl - 0
Tll
2510.15
2587.02
2740.71
2987.12
2519.89
2592.32
2745.17
2989.86
2545.33
2609.93
2758.45
2998.03
2581.12
2637.65
2780.19
3011.38
.05 , P2 =
T22
0.57
4.23
14.75
30.86
0.72
4.32
14.71
30.64
1.14
4.61
14.61
29.99
1.8
5.04
14.42
28.96
0.2 , e
T33
0.75
11.54
42.66
90.41
0.81
11.52
42.39
39.77
0.98
11.45
41.6
87.8
1.24
11.29
40.26
84.69
= 0.01 ,
T12
0
0
0
0
-15.39
-14.74
-12.14
- 7.37
-30.74
-29.29
-24.08
-14.6
-45.66
-43.34
-35.54
-21.52
K = 1
T13
0
135.
.i 275.08
421.07
0
13^.76
274.3
419.7
0
133.95
272.14
415.73
0
132.39
268.3
408.9
T23
0
0
0
0
0
-0.76
-1.18
-0.97
0
-1.5
-2.33
-1.93
0
-2.17
-3.38
-2.81
36
-------
Table 5.2 (continued)
(d) PI = p2 = 0.25 , e = 0.02 , K = 1
T22 T33 T12 T13 T23
0
5
10
15
259.92
263.06
272.34
287.58
0.28
0.39
0.79
1.41
0.37
0.78
1.97
3.79
0
0
0
0
0
8.43
16.92
25.42
0
0
0
0
37
-------
SECTION 6
APPLICATION EXAMPLES
Two hypothetical examples will show the applicability of the results
developed in Section 5. The examples focus on the estimation of the
concentration variances.
The Continuous Concentration Input
If solute mass is injected at a constant rate at the origin (x=y=z=0) in a
three-dimensional aquifer with constant mean groundwater flow velocity, the
resulting steady mean concentration distribution is approximated by the
following form when x» y,z:
r2 Z i / ^ i \
c(x) = m exp-t
where A__, A_. = lateral macrodispersivities along y, z
n = porosity
m = rate of mass injection
V = mean pore velocity (assumed along x) = q/n
x_ = (x,y,z) location vector
The components of the gradient of the concentration are given from,
2 2 .
f -—ill
iA33 x
(6.2)
The concentration variance is then given from the form in eq. (5.33c),
2 2 2 T $1 3i (6.3)
ffc " °f*3 ij 3x£ 3x
where the T . 's will be computed from the approximate analytical ex-
pressions in Table 5.1. The following parameters are used [x^ =x^ t*2 =v >X3=Z^
A = 1.0 m.
A22 = A33 = °'01 ™'
a = 0.01 m.
1 m.
af = 1
The concentration is normalized by the quantity c =m/[ nA A__Vj .
Figure 6.1 shows a cross-section of the plume at x=50 m. and z=0 m.
plotting the mean concentration plus and minus a standard deviation. It is a
38
-------
actual ccnc.?
c
0
c
a
u
o
a
o
a
= CD •
10
3.0
Figure 6.1 Normalized mean concentration and standard deviation at
x - 50 m, z = 0 m
39
-------
logarithmic plot which enables one to observe the increase of the coefficient
of variation at the tails of the plume. The fluctuating line is a sketch to
illustrate a possible concentration realization. At locations y=z=0 ra. due
to symmetry the gradients (3c/3y), (3c/3z) = 0; however, (3c/3x) # 0,
and it is responsible for the approximately 0.25 coefficient of variation
observed at this location. Figure 6.2 shows the isoconcentration lines and
the standard deviations in a horizontal plane located at z=0 m. This form of
plot allows one to visualize how much the contoured configuration of a plume
may vary. Figure 6.3 shows the correlation function of the concentration for
points located at the longitudinal axis of symmetry of
the plume. Along the flow directions extremely high correlation scales
are observed.
Along the center line of the plume, the coefficient of variation can be
written as
! 0'5 <6-4)
vx <• Jo'
c '
2
and for the parameters chosen with the isotropic relationship y=exp(af/6)
ac 7
— = — , x in m. (5.5)
™ X
c
Thus for distances higher than 100 correlation scales the coefficient of
variation will be a few percent. Figure 6.1 shows that this should be viewed
as the minimum concentration variability at a given longitudinal distance from
the source. The implied near source behavior with very large coefficient of
variation must be viewed as only a qualitative feature since the perturbation
approach used to develop the coefficients T.. (Section 5) was based on the
assumption that the concentration variations'1 were relatively small.
The Pulse Injection
As in the previous case, the form of the mean equation is assumed known,
with effective parameters resulting from the aquifer heterogeneity. The mean
concentration distribution for an instantaneous point injection at the origin
is
zUtt). L -rrtsaai.^.^J^w-M
f I ••\J*^/'» A A \ -^ / * ^" 11 1 1 1 "J
n(4irtV) (AnA22 33'
.Y
where M = mass injected
A , A,,, A = macrodispersivities along x, y, z axes respectively
V = mean pore velocity - -'„-
£= (x.y.z) = (Xj.X2.X3)
Denote Vt = x , the location of the center of mass of the plume at any
time. The concentration gradients are given from
40
-------
c- .00001
CD
t-«
X
a
x -axis
Figure 6.2 Isoconcentration lines and standard deviations at
z = 0 m
-------
5
5 2
a i
d
3.
X-AXIS
iso-corrQlation contours
i — o- ®- — ®- — — — — Q --- .
0.3
(D- O- -4D— —• — -©— —
0.5
> o -o .
0. 7 ~
(9 o
0.9
rtl.c) 4.0 0.0 12.0 111. 0
X-AXIS
20. U
6.3 Correlation function for a point at y - z = 0
(a) Correlation along x
(b) Correlation along lines in the (x,y) plane at
various angles with x
(c) Lag-correlation contours In the (x,y) plane
-------
x-x
3c _ - c
3x ~ c 2A
The concentration variances were calculated using the numerical evaluation
for an anisotropic heterogeneous aquifer with parameters,
2
£L= 3m , i2= 1m , H3= 0.15ra , af = 0.5 , o = 0.01m
The macrodispersivities are
A = 0.3ra , A = A__= 0.01m
and the values for T.., T_2, and T~_ are taken from Table 7.4.
The concentration is normalized by
V M/[(4ir)3nAuA22A33]
Figure 6.4 and 6.5 show the mean plus ana minus, a standard deviation for a
longitudinal and vertical section respectively along the axes of symmetry of
the plume. Note that at the maximum concentration a =0. This is
because the concentration gradient is zero at this location. The analysis was
developed assuming that the concentration gradient is constant. This is
clearly not true around the concentration peak, so it is not expected that the
results hold at this region.
Figure (6.6) shows the concentration variance in a vertical section 5 m.
away from the concentration peak in the x direction and 1 m. in the lateral.
The minimum coefficient of variation (at the peak) is not zero anymore.
-------
Longitudinal Section
o
a
42 43 44 45 46 47 48 49 50 51 52 S3 54
-3
36 37 38 39 40
Figure 6.4 Normalized concentration and one standard deviation
interval along a longitudinal section (line coordinates
y = z = 0)
-------
Vertical Section ( x=xc )
a>
o
i l i i l i i l i
-2-1.8-1 6-1 4-1.2-1-08-0.6-0.4-0.2 0
i i
0.2 0.4 0.6 0.8
1 1 2 1 * 1 6 1.8 2
Figure 6.5 Normalized concentration and one standard deviation
interval along a vertical section (line coordinates
x = xc - center of mass of plume, y = 0)
45
-------
Vertical Section (x=40,z=1)
1.8-1 fi-1 4-1.2-1-08-06-04-0.2 0 0.20.40.608 1 1.21.41.618 2
Figure 6.6 Normalized concentration and two standard deviation
interval along another vertical section (line
coordinates x = xc - 5 m, z = 1 m)
46
-------
SECTION 7
FIELD DATA ANALYSIS: THE CFB BORDEN SITE
The results from the stochastic theory will be applied in the controlled
field scale experiment at Canadian Forces Base, Borden, Ontario. This is a
preliminary application since the available data represent two sampling dates
only, namely, October 25, 1982 and November 28, 1983 (Freyberg, 1984, personal
communication). The tracer site is located in a sand quarry down-gradient
from an abandoned landfill (Sudicky, 1985). The tracer injection occured on
August 23, 1982, and comprised of 10.7 kg of chloride solution injected at an
average concentration of 892 tng/1 from a zone between 2.0 m and 3.6 m below
ground surface (Figure 7.1). Table 7.1 summarizes part of the existing
information about the site.
The stochastic theory can be applied to produce estimates of:
(i) effective hydraulic conductivities;
(ii) macrodispersivities; and
(iii) variance and covariance of the concentration deviations.
The input data set required must include:
(a) the mean log-hydraulic conductivity, F, and its standard deviation,
af.
(b) correlation scales of the log-hydraulic conductivity field, £nK;
(c) direction of its principal axes; and
(d) local dispersivicles.
Prior to analyzing the data a preliminary investigation of their information
content is required.
Table 7.1
EXISTING INFORMATION ON THE BORDEN SITE PARAMETERS
a2£ = 0.3 - 0.5
fc. 3 2 - 5 m; i2 s 2 - 5 m; fc_ = 0.10 m
AU s 0.30 m; AZZ =» 0.034 m; A33 - « 1* m
= F = KG = 7.15 x 10~3 cm/s
= -5°
Approximately zero was assumed by Sudicky (1985) without actual analysis
of the vertical spreading; quantitative analysis of data from an earlier
test at the Borden site (Gelhar et al. 1985) shows that vertical spreading
is greater than molecular diffusion (A^ = 1 to 2 mm).
47
-------
60
40 r
g
x
20 L
1 20 .'A'
»—'-4
\
dwa
+ Multilevel Somolers
^ Injection Wells
a Piezometer
• Core Location
-20
20
40
Y(m)
Figure 7.1 Instrumentacion, locacion of cores and hydraulic head
distribution on August 10V 1984 (solid line) and
November 15, 1984 (dashed line) at the tracer-test site
(after Sudicky, 1985b). The circled samplers show the
locations of the wells used in this analysis. V
denotes the well used for the vertical section
analysis.
48
-------
MAX
MIN
31 S
O
CONTOURS
1
2
3
4
5
B
2. 5
ID
6O
100
2OO
260
vo
Figure 7.2 Isoconcentration contours ac depth z = -5.31m, from the
Borden site data (original data from D. Freyberg,
personal communication)
-------
Averaging Effects in Sampling
The theory developed assumes point values of all the processes, where as
it is traditional in groundwater studies, a "point" is defined in a
Representative Elementary Volume (REV) sense, or, at the same scale that the
convective-dispersive equation for the mass transport holds. However, the
data collected are directly affected by the resolution of the sampling
method. For the hydraulic conductivity data at the Borden site, the values
were determined from cores 0.05 m. in height and 0.05 m. in diameter. Even at
this scale the sampling method introduces an averaging of the true process
which results in reduced variances and increased correlation scales (Vanmarke,
1983). A similar analysis was performed at the hydrothermal test site
Aefligen, in Switzerland (Hufschmied, 1985). The so-called "variance
function", x(a)» describes the variance reduction in terms of the correlation
scale. For a one-dimensional case the measured process is given by,
. z + Az/2
f(z) = i / f(z) dz (7.1)
AZ
z - Az/2
and its variance is:
2 2
J_ = Of
f
(7.2)
With the exponential model proposed by Sudicky (1985a) for the vertical
correlation and with i = O.I m, Az = 0.05 in (the height of the sample),
x(0.05) s 0.85. Therefore, the variance of the point values of the
log-hydraulic conductivity should be increased by 15%.
The Spectrum of InK
The estimation of the concentration variance requires a spectrum that
enforces stationarity on the underlying process. Such a spectrum has a zero
integral scale in three dimensions (Section 5). The exponential covariance
function does not possess this property, and this is the reason that a
"hole-type" covariance function was used in the theoretical development. The
difference between these two covariance functions (see Figure 5.3) is so small
as to be practically indistinguishable in terms of statistical inference and,
as will be shown in the following, the shape of the spectrum has minimal
effects on the estimation of the effective hydraulic conductivities and the
macrodispersivities.
50
-------
The important parameter is the correlation scale which is defined by the
lag corresponding to the e~l level of the correlation function. In the
isotropic case the correlation function used is,
p(s) = (1 - y s) e"S , s = | (7.3)
and the e-1 level is given by s = 0.723,
p(0.723) = e"1 , Cg = 11(0.723) (7.4)
Notice that the correlation scale of the above function is different than
the length scale dividing the lag. The opposite holds for the exponential
correlation function. It is easy to verify, since the exponential correlation
function is given from
PE
that
e"S' , s' =f , (7.5)
PE(1) = e"1. 5e = X. (7.6)
Consequently, the correlation scales determined by Sudicky (1985b)
correspond to the X's (in the notation used here), and the correlation scale
for the hole type correlation function, which is used in this analysis, is
given from,
0.723
(7.7)
The above relationship can be expanded to the anisotropic case also, where
the following will hold:
1. = Xi , i = 1,2,3 (7.8)
0.723
Therefore, the values given in Sudicky (1985) Xi=X2=2 m> X3=0.1 m,
would correspond to
A1 = 12 = 2.76 m, Z3 = 0.14 m (7.9)
With these preliminary remarks, required because of the sampling technique
used and the covariance model fitted, the stochastic theory can now be applied
to estimate the quantities mentioned in the introduction.
Effective Hydraulic Conductivities
Gelhar and Axness (1983) have developed the general equations for the
effective hydraulic conductivities in terms of the spectrum of the
log-hydraulic conductivity. The analysis evolved around a linearized version
of the flow equation and was heuristically generalized for the case of af > 1.
The effective hydraulic conductivities along the principal axes of flow are
given from
51
-------
Ku-Kt exptoj ({-g.t)], i = 1,2,3 (7.10)
(no summation over repeated indices), where K£ = exp [E(lnK)] , and the gi;j's
can be calculated from
Si-= / ~SV~ sff (^f) d- (7*1L)
^ k k' a
The primed system of axes is aligned with the principal axes of the
hydraulic conductivity. Freyberg (1985) has observed an approximate 5°
clockwise horizontal inclination of the principal axes of the plume with
respect to the mean flow. This is an indication of horizontal anisotropy in
the hydraulic conductivities, in contrast with the reported values of the
correlation scale, and a counterclockwise alignment of the principal axes of
the plume with respect to the mean flow, or, plume trajectory. Table 7.2
summarizes a series of calculations where the input parameters were the first
seven columns, i.e. of, 4 , A,, *., a = local dispersivity, * = horizontal
angle of x1 and x, 6 = vertical angle of xj and x (where x[ = major principal
axis of hydraulic conductivities, x_= mean flow direction). In columns 8,9,
10 the estimated values of K ,/K., K /'< and K3,/X do not show significant
variations and compare quite well witn tne values of 1.18 and 0.85 reported by
Sudicky (1985). The high values correspond to variances not considered by
Sudicky (1985), i.e., of = 1.5, and a" = 1.0. It can also be observed
that the anisotropy in the effective hydraulic conductivities both in the
horizontal as well as the vertical plane is very mild even for aspect ratios
of the original field of ^/^ = 5/0.5 = 10.
Note that the coefficients gtj , Eq. (7.11), needed for the estimation of
the effective hydraulic conductivities in Table 7.2 were evaluated numerically
for the hole type covariances function, Eq. (5.32). No significant difference
from the values shown in Gelhar and Axness (1983) was observed.
Asymptotic Macroscopic Dispersion
Gelhar and Axness (1983) have developed expressions for the asymptotic
macrodispersivities for the general case of arbitrary orientation of
stratification. These expressions are quite complicated to be evaluated
analytically except for the case of:
hydraulic conductivity isotropy in the horizontal or vertical plane,
negligible local dispersivity,
exponential log-hydraulic conductivity covariance, where approximate
solutions can be obtained.
52
-------
in
U>
Table 7.2
EFFECT OF THE INPUT PARAMETERS ON EFFECTIVE CONDUCTIVITIES
i.
If
I
0.4
1
1
1
0.5
0.5
0.5
0.5
1.5
1.5
1.5
1.5
0.5
0.5
0.5
0.5
0.5
1
1
1
1
Aj (m.)
5
5
2
2
2
2
2
2
2.5
2.5
5
5
2
2
2
2
5
2.5
5
2
2
*2 (m.)
0.5
0.5
0.5
2
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
1
1
*3 (m.)
• 0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
a (m.)
0.01
0.01
0.01
0.01 '
0.01
0.01
0.01
0.01
0.01
0.001
0.01
0.001
0.01
0.001
0.01
0.001
0.001
0.001
0.001
0.001
0.001
4>
0
5
10
5
0
5
0
0
20
20
20
20
20
20
45
30
25
25
20
45
20
Q
33
5
5
0
10
0
0
5
10
0
0
0
0
0
0
0
0
0
0
0
0
0
K/ If
1 1 ' A
11 I
1.22
1.64
1.61
1.59
1.27
1.27
1.27
1.27
2.07
2.07
2.1
2.1
1.27
1.27
1.27
1.27
1.28
1.62
1.64
1.6
1.6
K / K
^11' ^«
22 4
1.14
1.4
1.41
1.59
1.18
1.18
1.18
1.18
1.66
1.66
1.65
1.65
1.18
1.18
1.18
1.18
1.18
1.4
1.39
1.52
1.52
K /K
33 I
0.87
0.72
0.72
0.65
0.85
0.85
0.85
0.85
0.61
0.61
0.61
0.61
0.85
0.85
0.85
0.85
0.85
0.72
0.72
0.67
0.67
-------
Using the same set of input parameters as in Table 7.2, a numerical
evaluation of the asymptotic macrodispersivities was performed assuming the
covariance function of the inK is the hole type, and results are shown on
Table 7.3. Recall from Table 7.1 that the estimated longitudinal dispersivity
for the side is 30 cm. , and that there is a clockwise rotation of the plume
with respect to the mean flow. These two characteristics will be the
calibration targets for the selection of the correlation scales.
Table 7.3 shows the different combinations of input parameters that
produce a longitudinal macrodispersivity of the same order of magnitude as the
observed ones. It is more difficult to match the 5° clockwise rotation; Note
that most of the input combinations produce angles between 1° and 3°. The
predicted asymptotic lateral macrodispersivities are smaller -by 2 to 3 orders
of magnitude. In that case the local dispersivity, which is additive, will
be the prevailing value. Gelhar (1986) in an analysis of the developing
macrodispersivities has shown that the lateral dispersivities attain a maximum
value and then drop to their asymptotic magnitude. Dagan (1984) derived a
similar result following another approach. A plausible speculation,
therefore, is that the measured field value for A22 of about 3 cm, 38 m
downstream from the injection point, is still influenced by the developing
transverse dispersion process and, hence, is much larger than its asymptotic
value, which is the one that the spectral theory estimates.
With the information presented in Tables 7.1, 7.2 and 7.3, it is
considered that a representative set of values for the parameters of the CFB
Borden is as in Table 7.4. The correlation scales have been adjusted for the
hole type covariance. The variance is slightly increased to account for
sampling averaging, and anisotropic behavior is included t'o produce the
longitudinal macrodispersivity.
Variance of the Concentration Field
The theoretical analysis of Section 5 has developed the following form for
the concentration covariance,
RCC - °f 43 T1J(p1,P2,c,K,+ ,8;£) G^. (7.12)
The analysis has shown that the controlling local dispersivity is a^.
(KE = oT/A3) and that the effect of , 6 is minimal. Therefore as a first
approximation the variance of the concentration field will be calculated from:
° " Gf *3 Tiii'V6'''5- ii
GG
c
-------
Ul
Ul
Table 7.3
EFFECT OF THE INPUT PARAMETERS ON MACRODISPERSIV1TIES
2
°f
0.4
1
1
1
0.5
0.5
0.5
0.5
1.5
1.5
1.5
1.5
0.5
0.5
0.5
0.5
0.5
1
1
I
I
tj.)
5
5
2
2
2
2
2
2
2.5
2.5
5
5
2
2
2
2
5
2.5
5
2
2
££.)
0.5
0.5
0.5
2
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
1
1
t£.)
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
a(m)
~~~"™
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.001
0.01
0.001
0.01
0.001
0.01
0.001
0.001
0.001
0.001
0.001
0.001
*
0
5
10
5
0
5
0
0
20
20
20
20
20
20
45
30
25
25
20
45
20
6
5
5
0
10
0
0
5
10
0
0
0
0
0
0
0
0
0
0
0
0
0
M->
0.211
0.276
0.209
0.187
0.187
0.182
0.167
0.126
0.175
0.28
0.185
0.31
0.14
0.224
0.099
0.181
0.236
0.264
0.316
0.306
0.383
Aj,
0.27
0.88
0.38
0.72
0.24
0.25
0.32
0.48
0.91
1.3
l.l
1.9
0.31
0.12
0.27
0.16
0.24
0.79
0.84
0.19
0.11
A33
0.62
1.5
0.79
0.27
0.54
0.56
0.77
0.12
1.1
0.34
1.1
• 0.33
0.76
0.19
1
0.21
0.2
0.31
0.3
0.27
0.24
A1
Al
0
-3
-2
0
0
-0
0
0
-8
-18
-1
—2
-0
-4
-2
-4
-6
-13
-15
-7
-5
2
.1
.4
.13
.56
.8
.1
.4
.94
.7
.7
.6
A1
A13
4.3
14
0
6.6
0
0
2.9
5.9
0
0
0
0
0
0
0
0
0
0
0
0
0
Aft.)
0.211
0.276
0.209
0.187
0.187
0.182
0.167
0.126
0.175
0.28
0.186
0.312
0.141
0.224
0.099
0.181
0.236
0.265
0.317
0.307
0.384
A*1
A22
0.27
0.63
3.6
0.72
0.24
0.25
0.32
0.48
0.46
0.13
0.49
0.12
0.3
0.049
0.23
0.045
0.047
0.092
0.09
0.028
0.032
A*1
A33
0.53
0.88
0.79
2.5
0.54
0.56
0.72
0.92
1.1
3.4
1.1
0.33
0.76
0.19
1.1
0.21
0.2
0.31
0.3
0.28
0.24
*
0
-0.66
-0.66
0.04
0
0
0
0
-2.9
-3.8
-3.4
-4.4
-3.8
-1
-1.2
-1.5
-1.6
-2.9
-2.8
-1.3
-0.85
n
1.2
3.1
0
2.1
0
0
1
2.7
0
0
0
0
0
0
0
0
0
0
0
0
0
1 Units are in
mm.
-------
Table 7.4
PARAMETERS USED IN THE VARIANCE ESTIMATION
a2 = 0.5
f
1. = 3.0 m; 4, = 1.0 m; l = 0.15 m
x » cu =a=0.01m.
Tu = 776.09; T22 = 1.03; T33 = 0.88
Table 7.5
CONCENTRATION GRADIENTS IN THE VERTICAL SECTION
AT (x,y) = (38.3 m., 16.08 m.)
z (m.)
-3.51
-3.81
-4.11
-4.41
-4.71
-5.01
-5.31
-5.61
-5.91
-6.21
-6.51
3c
3x
0.185
0.981
4.824
1.25
6.838
6.64
18.21
1.7
2.16
-13.79
-2.83
ii
ay
-.707
-4.14
-5.61
-7.77
3.64
-4.92
-1.75
1.934
9.731
20.78
20.48
*
R
0.3788
0.579
0.627
0.641
0.411
0.205
0.208
0.087
0.508
0.735
0.825
* R = coefficient of multiple determination, which is equivalent to
the correlation coefficient in univarlate regression analysis
(Neter and Wasserman, 1974).
56
-------
The values for the coefficients T are shown in Table 7.4. They were
calculated using the numerical integration programs developed under the
general name TIJ_ND (where I, J = 1,2,3). The next task in the calculation
of the variance of concentration is the estimation of the mean gradient. Two
sections will be shown, an approximate longitudinal one and a vertical one
(Figure 7.1). The approximate longitudinal section, follows a line at
z = -5.31 m and is rotated at an angle of about 15° clockwise from the
principal axes of the computed macrodispersivities. Although at this depth
the choice of this direction appears justified (Figure 7.2), it was mainly
necessitated by the location of the sampling wells, the objective being to
use the original raw concentration data. Consequently, this direction will
be considered as the x direction for the calculation of the gradients. The
vertical gradient was estimated through analytical curves, of a Gaussian type,
that were fitted to the data at each well along x. The gradient along y was
estimated roughly from the plotted contours. For the vertical section, at
(x,y) = (38.3 m, 16.08 m), the same Gaussian curve fitting as before was used
to estimate the vertical gradient, while the gradients along x and y (Table
7.5) were estimated by fitting a regression plane through that point and its
eight surrounding ones (Neter and Wasserman, 1974). These procedures are
summarized in Figure 7.3.
The longitudinal section with the fitted mean is shown on Figure 7.4.
The series of Figures 7.5(a) through (i) show the fitted curves in the
vertical direction along the wells. Appendix D describes the parameters for
each curve. Finally, Figure 7.2 is used to estimate the gradients along the
y direction. All this information is collected in Table 7.6 and with the
T..'s presented at Table 7.4 the variance of the concentration can easily
be^estimated. Figure 7.4 shows the (c + a , c - a) interval which will
include approximately 67% of all measurements ff concentration is normally
distributed.
Table 7.7 shows the data used for the vertical section at (38.3, 16.08).
The gradients 3c/3x, 3c/3y are calculated from the linear regression as
described above. The c +_ a interval is plotted in Figure 7.6. The effects
of the vertical gradients on the variance are quite strong in this case. A
three-order of magnitude difference in the coefficients T.. and T,_ is
balanced by a four-orders of magnitude difference in the squared gradients
(Table 7.7). The lateral horizontal gradient has a minimal contribution. The
logarithmic version of the previous plot shows the increase of the coefficient
of variation at the tails (Figure 7.7).
57
-------
— : gradient along this direction was estimated by a continuous curve fit
: gradient along this direction was estimated from contour plotting
-R : gradients on this plane were estimated from linear regression
Figure 7.3 Schematic presentation of the methods used for the
estimation of the mean concentration
58
-------
L.ONG I TUD I NAL_ l_ I Nl
Ul
VO
I ' ' ' I ' ' ' I
5 3O 35
X
Figure 7.4 The measured, fitted and plus-minus one standard
deviation concentration profiles along the longitudinal
section at z = -5.31rn
-------
-334 -3.54 -3 64 -4 14 -4.44 -4 74 -504 -S.34 -5 fr» -£• ft ~« 24 -« &4 -6,64
-3.21-3 SI -381-4.11 -4 41 -4 71 -5 01 -S31-SC1-SV1 -«.21-8.51-e.Bt-7.11
-3 IB -34» -37B —4 OB -4 30 -4 «B -4 Co -£. X -5 S» -ffc. -« IB -6*
-31! -342 -372 -4 CC -132 -462 -462 -E 22 -S 52 -£. 82 -« 12 -642 -e 72
Figure 7.5 The mean concentration curves fitted at the wells on
the longitudinal section (c in the vertical vs. z)
+=fitted,O =measured
-------
-312 -342 -3.73 -4 CO -t 33 -4 «3 -4 S3 -122 -513 -S B2 -« 12 -e 42 -« 72
-343 -373 -4O5 -4.33 -4 «3 -4 B3 -123 -5.S3 -6.B3 -«.13 -« 43 —« 73
Figure 7.5 (continued) The mean concentration
curves fitted at the wells on the
longitudinal section
-3 14-3 44-3 74-4 O4-4 34-464—« 64-5 34-5 54-S b< -fl 14-844-674
-------
Table 7.6
ESTIMATION OF THE RANGE OF VARIATIONS
IN THE LONGITUIMNAL SECTION
x (in.)
27.25
31.3
35.3
38.3
40.3
42.3
44.3
46.22
48.35
50
52
y (m.)
6
9
12
16.08
17.94
19.95
22
24
25.95
28
30
z (m.)
-5.39
-5.31
-5.34
-5.31
-5.29
-5.22
-5.22
-5.23
-5.24
-5.3
-5.3
c (mg/z)
(meas. )
1.6
2.3
31.1
169
319
243
208
123
19.3
3
2
c (mg/t)
(fitted) '
2.82E-4
0.17
14.76
183.18
312.97
296.61
149.67
42.20
5.98
0.57
0.06
3c/ 3c/
0.0004
0.18
10.23
55.29
28.34
-38.11
-52.36
-23.81
-4.71
-0.57
-0.07
0
0
0
0
7
0
0
0
9.3
0
0
3c
3z
0
0
15.57
86.96
207.06
233.59
129.9
23.43
35.65
0
0
«c
0.001
0.55
30.29
163.62
86.26
115.00
155.26
70.39
14.40
1.68
0.23
0.001
0.73
45.06
346.80
399.23
411.62
304.93
112.59
20.38
0.29
0.29
«=
-0.0094
-0.38
-15.5
19.55
226.70
181.61
-5.59
-28.19
-8.42
-1.11
-0.16
-------
Table 7.7
ESTIMATION OF THE RANGE OF VARIATIONS
IN THE VERTICAL SECTION
z(m.)
-3.21
-3.51
-3.81
-4.11
-4.41
-4.71
-5.01
-5.31
-5.61
-5.91
-6.21
-6.51
-6.81
-7.11
c(mg/fc) c(mg/A) 3c/3z
(measured) (fitted)
0.4
4.2
9.2
16.4
63.6
112
254
169
306
200
132
13.1
5.6
1 ' 7.8
0.12
0.88
4.84
19.85
60.6
137.9
234.06
296.06
279.17
196.24
102.83
40.17
11.69
2.54
0.85
5.44
25.1
83.57
195.82
310.64
297.92
86.96
-191.33
-326.63
-271.85
-145.53
-53.86
-14.17
3C/3X
0
0.18
0.98
4.82
1.25
6.83
6.64
18.21
1.7
2.16
-13.79
-2.83
0
0
3c/3z
0
-0.70
-4.14
-5.61
-7.77
3.64
-4.92
-1.75
1.93
9.73
20.78
20.48
0
0
°c
0.08
0.77
3.86
. 16.52
19.90
-37.00
35.62
54.50
19.74
33.23
49.00
16.90
5.37
1.41
**«.
0.20
1.65
8.70
36.37
80.51
174.95
269.68
350.56
298.91
229.47
151.84
57.076
17.071
3.95
c -a
i
0.03
0.10
0.98
3.32
40.69
100.93
198.43
241.55
259.42
163.00
53.82
23.26
6.32
1.12
-------
VERTICAL. l_ I Nl
I I
Figure 7.6 xiie measured, fitted and plus-minus one standard
deviation concentration profiles along the vertical
section at (x,y) = (38.3, 16.08)
-------
ooo
i oo
o 10
QL
Z 1
LJ
O
Z
O
o i
O 1
VERTICAL LINE
1 ' ' I ' '
T T
J L
l I I
X
Figure 7.7 The values plotted in Fig. 5.6 in a serai-logarithmic
plot
-------
Conclusions
The analysis of the field data demonstrates the validity of the stochastic
theory results for the concentration variance. The predicted concentration
intervals defined by +_ a are consistent with the observed fluctuations of
concentration around a smooth mean and capture most of the measured values.
The effect of the longitudinal concentration gradient in the calculation of
the variance is dominant except when focusing on vertical sections where the
vertical concentration gradient is higher by 2 or 3 orders of magnitude.
The estimated effective hydraulic conductivities are quite close to the
observed geometric mean of the measured hydraulic conductivities. The
estimated longitudinal macrodispersivity is also comparable to the one
estimated through the method of moments.
The analysis reveals the need for the development of a formal procedure
for the estimation of the mean. Analysis of data from subsequent sampling
dates is needed to show the robustness of the results.
66
-------
SECTION 8
DISCUSSION OF RESULTS
A simple theoretical result was developed to quantify the concentration
variability. When measurements of the concentration exist, it can be thought
of as giving a measure of the range of expected variation of the concentration
around a smooth mean concentration distribution. The theory is also useful as
a predictive tool when the solute distribution is extrapolated for distances
far downstream using a classical dispersive transport model.
The functional form of the concentration variance is easy to implement.
It is a local relationship that involves the mean concentration gradient, the
variance and correlation scales of the log-hydraulic conductivity field, and
the local dispersivity values. The mean concentration gradient can be
estimated either in a predictive context with numerical or analytical models,
or through a simple curve fitting with existing plume measurements. The
former approach was demonstrated in Section 6 and the later by the analysis of
the Borden site data in Section 7.
The variance of the log-hydraulic conductivity can be measured in situ or,
as more information for the behavior of different types of aquifer materials
is accumulated, it can be estimated from information about geologically
similar sites.
The local dispersivities can be also found from a laboratory test of soils
from the site or estimated based on the general properties of the soil (Klotz,
et al. 1980). It should, be noted that the effect of the local dispersivity on
the variance is significant and cannot be neglected. This behavior is in
contrast to that of the macrodispersivity which is not sensitive to the value
of the local dispersivity.
An important factor needed for the variance estimation is the correlation
scales. Analyses such as Sudicky's (1985a,b) show ways to estimate their
values in situ. It is very encouraging that the estimated correlation scales
give a consistent picture for the mean plume when used to estimate the
macrodispersivities and the concentration variance as a measure of the
variation around the mean plume. Gelhar (1985) presents a collection of the
variances and correlation scales for log-hydraulic conductivity or
log-transmissivity estimated either from direct analysis of measurements or
through inverse modeling procedures (Hoeksema and Kitanidis, 1985).
The predicted concentration variance is an appropriate calibration target
for the numerical models. Other sources of error such as measurement errors
or discretization errors associated with the numerical scheme should also be
considered, but typically those effects are not expected to be the predominant
ones. Therefore, if the root mean squared difference between the observed and
simulated concentration fields is used as an index of the model error, the
calibration of the numerical model could be considered adequate when the model
error and the error (standard deviation) predicted from the stochastic theory
are of the same magnitude.
67
-------
Although the results of the comparison with field data are encouraging and
indicate the workability of the stochastic theory to predict concentration
variance, there are a number of approximations and/or refinements which need
to be evaluated before the full range of applicability of the theory can be
established. These features are discussed in the following paragraphs.
The development of the unsteady concentration case in Section 5, invokes
an expansion of the concentration gradient around the current time t. Although
only the first two terras of the expansion have been retained, additional terms
could easily be included and evaluated numerically. Since these terms are
multiplied by the time derivatives of the mean concentration, it is
expected that they will have ajninimal effect on the variance except near
the concentration peak where ac/ax^ i = 1,2,3 changes rapidly. The effects
of these unsteady terms near the peak need to be explicitly evaluated in order
to understand their significance.
The t > o> assumption for the integration of Eq. (5.22) is not required.
Eq. (5.21) could be integrated for any finite time, which would result in an
additional factor [l-exp(-et)] multiplying the complex amplitude in Eq.
(5.23). From the form of Eq. (5.21), it is evident that the asymptotic result
will overestimate the variance for finite time and consequently will usually
be conservative. This finite time effect should be evaluated in future work.
In a heterogeneous aquifer, the local dispersivities are expected to be
variable in space. The effect of local variation of dispersivities can be
analyzed using an approach similar to that of Naff (1978). This effect may be
of significance in determining an effective local dispersivity to be used in
the calculation of the concentration variance.
The behavior of the concentration covariance along the mean flow direction
is somewhat puzzling from a physical point of view. As is shown in Figure
5.5a, the longitudinal correlation scale of the concentration variation
becomes very large. This implies that regions of high or low deviation from
the mean concentration will persist over very large longitudinal distances,
say up to 10 or 100 times the correlation scale of ink. Note that the degree
of longitudinal persistence increases as e=a/fc, decreases. Although
somewhat similar behavior seems to be indicated by other methods (Dagan 1984),
we feel that this behavior is questionable in terms of our intuition about the
physical system. Intuitively we expect a' smaller degree of longitudinal
persistence. This feature requires additional detailed investigation because
the concentration covariance function is a very important factor in the design
of plume monitoring networks. Preliminary results indicate that the high
persistence is caused by the distribution of the spectrum energy at high wave
numbers.
68
-------
REFERENCES
1. Ababou, R., D. McLaughlin, L.W. Gelhar and A.F.B. Torapson. Numerical
Simulation of Saturaced/Unsaturated Flow Fields in Randomly Heterogeneous
Porous Media. Proceedings of International Association of Hydraulic
Research Symposium on the Stochastic Approach to Subsurface Flow,
Montvillargenne, France, June 1985.
2. Aldama, A.A. Rodriguez. Theory and applications of two-and three-scale
filtering approaches for turbulent flow simulation. Ph.D. dissertation,
Massachusetts Inst. of Tech., 397 pp., November 1985.
3. Anderson, M.P. Using Models to Simulate the Movement of Contaminants
Through Groundwater Flow Systems. Critical Reviews in Environmental
Controls, 9(2), pp. 97-156, 1979.
4. Bakr, A.A. Stochastic Analysis of the effect of spatial variations of
hydraulic conductivity of groundwater flow. Ph.D. dissertation, New
Mexico Institute of Mining and Technology, 244 p, 1976.
5. Bakr, A.A., L.W. Gelhar, A.L. Gutjahr and J.R. MacMillan. Stochastic
Analysis of Spatial Variability in Subsurface Flows 1. Comparison of One-
and Three-Dimensional Flows, Vol. 14, No. 2, pp 263-271, April 1978.
6. Beran, M.J. Statistical Continuum Theories, Interscience, New York, 1963.
7. Betson, R.P., J.M. Boggs, S.C. Young, and L.W. Gelhar. Design of field
experiment to investigate ground water saturated zone transport processes,
Electric Power Research Institute, Inc., Interim Report RP 2485-05, 104
pp, 1985.
8. Chirlin, G.R. and G. Dagan. Theoretical head vartograms for steady flow
in statistically homogeneous aquifers. Water Resources Res. 16(6), pp
1001-1015, 1980.
9. Dagan, G. Models of groundwater flow in statistically homogeneous porous
formations, Water Resources Res. 15(1), 47-63, 1979.
10. Dagan, G. Stochastic Modeling of Groundwater Flow by Unconditional and
Conditional Probabilities, part 1, Conditional Simulation and the direct
problem. Water Resour. Res., 18, pp. 813-833, August 1982a.
11. Dagan, G. Stochastic Modeling of Groundwater Flow by Unconditional and
Conditional Probabilities. 2. The Solute Transport. Water Resour. Res.,
Vol. 18, No. 4., pp. 835-848, August 1982b.
12. Dagan G. Solute Transport in heterogeneous porous formations. J. Fluid*
Mechanics, 145, pp. 151-177, 1984. -*•
69
-------
13. Freeze, R.A. A Stochastic-Conceptual Analysis of One-Dimensional
Groundwater Flow in Nonunifonn Homogeneous Media. Water Resour. Res.,
Vol. 11, No. 5, pp. 725-741, October 1975.
14. Freyberg, D.L., D.M. Mackay and J.A. Cherry. Advection and Dispersion in
an Experimental Groundwater Plume. Proceedings of the Hydraulics Division
Specialty Conference in Frontiers in Hydraulics Engineering, A.S.C.E.,
pp. 36-41, July 1983.
15. Freyberg, D.L. A natural gradient experiment on solute transport in a
sand aquifer II. Spatial moments and the advection and dispersion of
non-reactive tracers. Water Resourc. Res., (in submittal) 1985.
16. Gelhar, L.W. Stochastic Analysis of Phreatic Aquifers. Water Resour.
Res., Vol. 10, No. .3, pp. 539-545, 1974.
. i
17. Gelhar, L.W. Effects of Hydraulic Conductivity Variations on Groundwater
Flows. Second International Symposium on Stochastic Hydraulics, Lund,
Sweden, in Hydraulic Problems solved by Stochastic Methods, Water Res.
Publications, Fort Collins, Colorado, pp. 409-431, 1977.
18. Gelhar, L.W. Stochastic Analysis of Solute Transport in Saturated and
Unsaturated Porous Media, (in press 1986).
19. Gelhar, L.W. Stochastic Subsurface Hydrology: from Theory to
Applications. Proceedings of International Association of Hydraulic
Research Symposium on the Stochastic Approach to Subsurface Flow,
Montvillargenne, France, June 1985.
20. Celhar, L.W. Stochastic Analysis of Flow in Heterogeneous Porous Media,
in Fundamentals of Transport Phenomena in Porous Media, J. Bear and M.Y.,
Corapcioglou, eds., Martinus Nijhoff Publishers, Dordrecht, The
Netherlands, 1984.
21. Gelhar, L.W., A.L. Gutjahr and R.L. Naff. Stochastic Analysis of
Macrodispersion in a Stratified Aquifer. Water Resour. Res., Vol. 15,
No. 6, pp. 1387-1397, Dec. 1979.
22. Gelhar, L.W., A.L. Gutjahr and R.L. Naff. Reply to Comments on
Stochastic Analysis of Macrodispersion in a Stratified Aquifer. Water
Resour. Res., 17, pp. 1739-1740, 1981.
>.
23. Gelhar, L.W. and A.L. Gutjahr. Stochastic Solution of the One-Dimensional
Convective-Dispersive Equation. Report No. H-ll, New Mexico Institute of
Mining and Technology, pp. 20, 1982.
24. Gelhar, L.W. and C.L. Axness. Three-dimensional stochastic analysis of
macrodispersion in aquifers. Water Resour. Res., 19, 161-180, 1983.
70
-------
25. Gelhar, L.W., A. Mantoglou, C. Welty and K.R. Rehfeldt. A Review of Field
Scale Physical Solute Transport Processes in Saturated and Unsaturated
Porous Media, Research Project 2485-5, EPRI EA-4190, Electric Power
Research Institute, Palo Alto, California, 116 pp., August 1985.
26. Gutjahr, A.L., L.W. Gelhar, A.A. Baler and J.R. MacMillan. Stochastic
Analysis of Spatial Variability in Subsurface Flows 2. Evaluation and
Application, Water Resour. Res., Vol. 14, No. 5, pp. 953-959, October
1978.
27. Gutjahr, A.L. and L.W. Gelhar. Stochastic Models of Subsurface Flow:
Infinite Versus Finite Domain and Stationarity. Water Resour. Res., Vol.
17, No. 2, pp. 337-350, April 1981.
28. Gutjahr, A.L.., E.J. Bonano and R.M. Cranwell. Treatment of Parameter
Uncertainties in Modeling Contaminant Transport in Geologic Media:
Stochastic Models vs. Deterministic Models with Statistical Parameter
Sampling. Proceedings of International Association of Hydraulic Research
Symposium on the Stochastic Approach to Subsurface Flow, Montvillargenne,
France, June 1985.
29. Gradshteyn, I.S. and I.M. Ryzhik. Table of Integrals, Series and
Products, Academic Press, 1980.
30. Hinze, J.O. Turbulence, McGraw-Hill, 1975.
31. Hoeksema, R.J. and P.K. Kitanidis. An Application of the GeostatistxcaL
Approach to the Inverse Problem in Two-Dimensional Groundwater Modeling.
Water Resour. Res., Vol. 20 (7), pp. 1003-1020,'1984.
32. Hoeksema, R.J. and P.K. Kitanidis. Analysis of the spatial structure of
properties of selected aquifiers, Water Resour. Res., 21, 563-572, 1985.
33. Hufschaiied, P. Estimation of three-dimensional statistically anisotropic
hyraulic conductivity field by means of single well pumping tests combined
with flowmeter measurements. Proceedings of International Association of
Hydraulic Research Symposium on the Stochastic Approach to Subsurface
Flow, Montvillargenne, France, June 1985.
34. Kitanidis, P.K., and E.G. Vomvoris. A Geostatistical Approach to the
Inverse Problem in Groundwater Modeling (Steady State) and One-Dimensional
Simulations. Water Resour. Res., Vol. 14(3), pp. 677-690, 1983.
35. Klotz, D., K.-P. Seiler, H. Moser and F. Neumaier. Dispersivity and
velocity relationship from laboratory and field experiments. Journal of
Hydrology, 45, 169-184, 1980.
71
-------
36. LeBlanc, D.R., Ed., Movement and fate of solutes in a plume of sewage
contaminated groundwater, Cape Cod, Massachusetts: U.S. Geological Survey
Toxic Waste Groundwater Contamination Program Papers presented at the
Toxic Waste Technical Meeting, Tucson, Arizona, March 20-22, 1984,
U.S.G.S. Open File Report 34-475, 1984.
37. Lumley, J.L. and H.A. Panofsky. The Structure of Atmospheric Turbulence,
John Wiley and Sons, Inc., New York, 1964.
38. Matheron, G. and G. de Marsily. Is transport in porous media always
diffusive? A counterexample. Water Resour. Res., 16, pp. 901-917, 1980.
39. McLaughlin, D.B. A Distributed State Space Approach for Evaluating the
Accuracy of Groundwater Predictions, Ph.D. Thesis, Dept. of Civil
Engineering, Princeton University, Princeton, 1985.
40. Mizell, S.A., A.L. Gutjahr and L.W. Gelhar. Stochastic Analysis of
Spatial Variability in Two-Diraensional Steady Groundwater Flow Assuming
Stationary and Nonstationary Heads. Water Resour. Res., Vol. 18, No. 4,
pp. 1053-1067, August 1982.
41. Moltyaner, G.L. Stochastic versus deterministic: A case study.
Proceedings of Intern. Assoc. of Hydr. Res. Symposium on the Stochastic
Approach to Subsurface Flow, Montvillargenne, France, June 1985.
42. Naff, R.L. A Continuum Approach to the Study and Determination of Field
Longitudinal Dispersion Coefficients. Ph.D. thesis, New Mexico Institute
of Mining and Technology, Socorro, New Mexico, pp. 176, 1973.
43. Netter, J. and W. Wasserman. Applied linear statistical models, R. D.
Irwin, Inc., Hornewood, Illinois, 842 p., 1974.
44. Schwartz, F.W. >4acroscopic Dispersion in Porous Media: The Controlling
Factors. Water Resour. Res., Vol. 13, No. 4, pp. 743-752, August 1977.
45. Smith, L., and R.A. Freeze. Stochastic Analysis of Steady State
Groundwater Flow in a Bounded Domain. 1. One-Dimensional Simulations.
Water Resour. Res., Vol. 15, No. 3, pp. 521-528, June 1979a.
46. Smith, L., and R.A. Freeze. Stochastic Analysis of Steady State
Groundwater Flow in a Bounded Domain. 2. Two-Dimensional Simulations.
Water Resour. Res., Vol. 15, No. 6, pp. 1543-1559, Dec. 1979b.
47. Smith, L. and F.W. Schwartz. Mass Transport 3. Role of Hydraulic
Conductivity Data in Prediction. Water Resour. Res., Vol. 17, No. 5, pp.
1463-1479, October 1981.
72
-------
48. Sudicky, E.A. Spatial Variability of Hydraulic Conductivity at the Borden
Tracer Test Site. Proceedings of Intern. Assoc. of Hydr. Res. Symposium
on the Stochastic Approach to Subsurface Flow, Montvillargenne, France,
June, 1985a.
49. Sudicky, E.A. A Natural-Gradient Experiment on Solute Transport in a Sand
Aquifer: Spatial Variability of Hydraulic Conductivity and its Role in the
Dispersion Process. Water Resour. Res., (Submitted, 1985b).
50. Taylor, G.I. Diffusion by Continuous Movements. Proc. of London Math.
Soc., Ser. 2, Vol. XX, pp. 196-212, 1921.
51. Vanmarke, E. Random Fields: Analysis and Synthesis. The MIT Press,
Cambridge, 1983.
52. Vomvoris, E. Groundwater Parameter Estimation: A Geostatistical
Approach. M.S. thesis, Univ. of Iowa, p. 188, July 1982.
53. Vomvoris, E.G. and L.W. Gelhar. Stochastic Analysis of Concentration
Variability in Transport in Aquifers. Proceedings of the Intern. Assoc.
of Hydr. Res. Symposium on the Stochastic Approach to Subsurface Flow,
Montvillargenne, France, June 1985.
73
-------
APPENDIX A
THE MATHEMATICS OF THE APPROXIMATE ANALYTICAL EVALUATIONS
This appendix describes the mathematical solutions used for the
approximate analytical results presented in Table 5.1. The following symbols
will be used repeatedly to provide non-dimensional quantities:
2 2 2^2
ul = klV U2 = V2' U3 = *3*3' U = Ul + "2 + U3
aL aT Ul ^3 ^3
.2 2 2 . -2
u = p. u. + U
,,2 2 2 ,,2 ,-,2 2 2 2 . -2 "2 222. f 2
V =eu +U , V =p1eu +U , V =p1eu +
-------
Changing Che variables with Che aid of Eq. (A.I)
2
III - /" --^O 22* 2*4 (1 - P? I?'2 duldu2du3 (A'3)
-»(l+u) p.u.+eu u
Introducing the variable u = u,/e, Eq. (A. 3) can be written as:
2 1 2 2 2
n +vV 2 2 2 + 2-4 Pl
(1 + V ) P u e + e V
Taking the limit e + 0,
..2 2^2 ..2
V ->• u? + u, = U
(. - Pl* >2 , i
22-4 22. 2-4
PL U +V +P1U +KU
and hence, 1^ takes the form
"
(A.4)
-oo (1 + U) P i U + 1C U
Integrating over u, after changing variables x = Plu,
/ ~~n w du = Z"2 / 2 . 2-4dx = p~ 2 i ~n = p~ ~f^ (A>b)
-co p,U+, u3 = U sin ij>,
I,, = — — f f 5—z- —=—z j ^— U dU d$ (A.8)
e pl 0 0 (I + U ) U (p2 cos % + sin",),)
75
-------
Evaluating separately, for each variable, [Gradshteyn & Ryzhik, 1980,
3.642(3)]
o U+u2)3 4 ' (A.9)
2ir .
o
0 p» cos 4> + sin § P2
which result after substitution in Eq. (A. 8) and Eq. (A. 2),
2
.2
3
°f -2 2 1 J_ JL (A.LO)
Variance; vertical isotropy - !„„.
T22 = E I22
I°2 = I22 («L= a, - 0) - /" —4-3 % duldu2du3 (A. 11)
-=> C 1 T u ; u
Changing to spherical coordinates (u, p, 9), with p2 = 1 ,
U» = U COScj)
u2 = u
u = u
2 2 2
+ sin 0 (1 - R x
(A. 12
sn « do (A. 13)
22 0 (1 + uV 0 0 (p^cos^d, + sin
The three integrals can be evaluated separately,
;" - ?L— du = Ig- (A- I4a)
0 (1 + u2)3 16
/ 7rCos2e d9 - « (A.Ub)
0
76
-------
= -4- + -4 + £ *n JLJLj£%
2Pl2 R2 2 R2 R R - 2R2
The combination of Eqs. (A.14a), (A.14b) and (A.14c) gives the result
presented in Table 5.1. Note that when A1=a2=A3=4» eq. (A.14c) becomes
f sin3* d* = ^ (A.15)
0 J
and, then,
Y
For T33 the only change occurs at Eq. (14b), which becomes
» T
/ sin 9 de = ir (A. 16b)
0
Thus, T22 = T33 as claimed in Table 5.1.
Covariance; isotropic case - T. .(5, 0,0)
From eq. (5.12),
(A. 17)
The spectrum Scc(lc) is even with respect to k. = (kj.kj.kj) for the case
considered here. The exponential can be expressed as,
ei-'-S-= cos(k.£)
and, since the sin(.) is odd it will cancel out. Therefore, the Tn(£) is
given from,
u2
Tu(l> - f- /"cosCu-i') 1 2.4 (1 - ^-)2Sff(u) du (A.19)
E-oa U + £ U U
Introducing u = VE* and takin8 the limit e * °'
C1 €2 £3 "424
77
-------
2 ,
(1 -^ e > * 1 sff^> - S£f(0,u2,u3> (A.20)
u
The uȣ argument in the cosine does not allow for analytical integration
unless only one of the Si's, or, equivalent ly, uj/s, is retained. Hence,
for £= (5,0,0), Eqs. (A.19), (A.20) and (5.32) for SfE(k),
= E / cos(euc) 224 - ~2~3 du du2du3 U<21)
-B> u + K U (1+U)
Integrating over u first,
'KU ^ (A.22)
/ cos(euc) -5 o-T du = -^T e
u + K V KIT
The remaining integral becomes,
CO 2
eK -co (I + U )
Converting to polar coordinates, u2 = U cos$, u3 = U sin
T,,(5,0,0) -E-J /*d4 /V8U " dU (A.24)
11 EK 0 0 (1 + U )
The second integral can be evaluated analytically in terms of the Exponential
integral [Gradshteyn and Ryzhik, 3.353(4)]
2
TU (5,0,0) =E — £[ 1-S-B eBEi(-B)] (A.25)
Note that the term in the brackets equals 1 when 3=0, which verifies the
variance relationship.
For S « 1,
B) = C + InB - B + ••• (A.26)
eB = 1 + g + ...
where C = 0.577215, Euler's constant [Gradshteyn & Ryzhik, 8.214 (1)] and
the term in the brackets becomes
[ 1 - B - B2 (C + Infi) + ...] (A.27)
For B » 1.
eBEi(-B) = - + -^ - \ + .... (A-28)
B B 6
78
-------
[Gradshteyn & Ryzhik, 8.215] and the term in the brackets becomes
i! i ('
Covariance; isotropic case - T..(0,5,0)
In this case Eqs. (A.17) through (A.20) are the same. With £ = (0,5,0),
Eq. (A.21) becomes
09 2
Tu(0,5,0) = E - / cos(u2?) ——^y£ —1~3 du du2 du3 (A.30)
From Eq. (A.22) the integral over u, for 3 = 0, ts known. Changing to polar
coordinates,* . i
«> 2ir
T ,(0,E,0) = E i- f f cos(Ucos^) 0 du d (A-31)
11 CK 0J 0 (L + U2)3
The integral over ,)> is given in Gradshceyn & Ryzhik, [3.715(18)],
2n
J cos(cos,bUc) d^, = 2 TT J (U;;) (A. 32)
0
and, finally,
j"j (U?) !L_ du = I ? ZK (?) (A.33)
from [Gradshteyn & Ryzhik 6.565(4), 8.486(16)] where Jo(.) and K2(.) are a
Bessel function of first kind, zero order and a modified Bessel function of
second kind, second order, respectively. The coefficient T22(0,5,0) is given
from,
2
(?) (A.34)
Y
Note that e2IC(i-) = 2, at 5 = 0, verifying the variance result.
t *•
For the case £ = (0,0,g) the only change is in eq. (A.32) which becomes,'"-
2*
/ cos(U? sin,),) 44 = 2n J (^U) (A.35)
0
[Gradshteyn & Ryzhik, 3.715(9)], therefore,
T11(0,g,0) = Tu(0,0,?) (A.36)
79
-------
Covariance; isotropic case - !„({;,0,0)
As in the variance case, the coefficient T22 can be evaluated when e = 0.
The same arguments for the exponential, Eq. (A.18), hold. Therefore, for
e = 0 and a combination of Eqs. (A.ll), (A.17) and (A.18);
2
« 1 u2
T.-Cj) = E / cosU £') Y~3 ~2~ di^du^u-j (A.37)
-w (1 + u ) u
Changing to spherical coordinates and for g_ = (g,0,0)
2
00 " ! u 32
T00(£) = E f / / cos (u^cosij)) =—=• sin ^cos 9 de d du (A.38)
22 0 0 0 U+uV
The integral over u equals,
u2 -lal 2
f cos(u£COSrji) 5—y du = yV- e I I [ 1 + a - a ] (A. 39)
o (i + u r ' '
where a = ^cose. In addition, the integral over 9 equals IT, from
Eq. (A.lAb). Hence
T.,(6,0,0) = E yj- |"V|a|[ 1 + lal - a2] sin3* d,j, (A.40)
Li 10 Q II
Changing variables to COS(J» = x, and using the following equalities,
/ bx.,, 2N . b,2 2 , ,1 2 >
f e (l-x)dx =e (-^ y) - (r r)
0 b b b
t1 bx,. 2N . b,2 6 . 6 -. , ,1 f> * ,. AM
J e (1 - x )x dx = e (-=• r- + -r) + (-J r> QA.41)
0 b b b b b
,l bx,. 2. 2. b,2 10 . 24 24. f-2 24
J e (l-x)xdx = e (—r r- + —7- r) + V—T + —?)
0 bbbbbb^^
the expression for T22(5,0,0,) becomes,
2
Y 5 SCC
2 3
Expanding, e ^ = 1 - 5 + -jjy - -iy + ..., shows that the term in the brackets
becomes for 5=0,
80
-------
[ TI 1 (A'43)
verifyiag Che variance result.
Covar lance; isotropic case - T22(0,c,0)
The general Eq. (A. 37) is the same. The conversion to spherical coordinates
is introduced in the following way,
U2 = UCOSij)
Uj = usin 0, approaches the value
1/3 ^, verifying thus the asymptotic result for the variance.
81
-------
Covariance; isotropic case - T33(0,0,£)
The general Eq. (A.37) is the same. If the following conversion to spherical
coordinates is introduced,
u3 = u cos
ut = u sinifcose (A.49)
u2 = u sinsin9
the expression for T22 (0,0,5) becomes,
« TI 2ir U2 32
T (0,0,£) = E f f f cos (u^coscjj) ~ « sin ijisin 0 d9d(((du (A.50)
22 000 (1+u )
Since
2ir - 2ir 2
f sin e de = / cos e de = IT (A.51)
0 0
It can be easily seen that
T22(0,0,0 = T22(g,0,0) (A.52)
Covariance; isotropic case - T33(5,0,0)
The general equation is:'
2
oa U-
T^^(t.O.O) =4 f cos (u,c) T-7—2 du1du2du3 (A.53)
— o»
(1 + U'
u2 = u sincose (A. 54)
Changing to spherical coordinates,
u1 = u cose
u2 = u sin^cose
u3 = u sinij)Sin0
Eq. (A. 53) becomes identical to Eq. (A. 50). Hence,
T33U,0,0) = T22(5,0,0) (A.55)
Covariance; isotropic case - T33(0,{;,0)
The general equation is,
2
U_
T (0,5,0) - / cos(u C) - —r^ -
82
du du du (A.56).
-------
Changing to spherical coordinates as in Eq. (A. 44), it can be seen that
T33(0,c,0) = T33(?J0,0) (A.57)
Covariance; isotropic case - T33(0,0,£)
The general equation is,
2
T33(0,0>5) = / cos(u3c) - T^~ 2 du!du2du3 (A<58)
Changing the spherical coordinates as in Eq. (A. 49),
r. =• TT 27T 2
cos •{•sin^ded^du (A. 59)
83
-------
APPENDIX B
THE INTEGRATION OVER WAVE-NUMBER
The spectrum of the concentration can be written as follows,
Eq. (5.10) and (5.11),
k?4J S (k) kk kk
•«<» - * 7^3 -^ I W-VV ^> V ^'l »•»
1 4 123
where E = -=-r- -r- - 5 — , and summation over j,£,ra,n in the branches is assumed.
T *
Note that introduction of spherical coordinates makes the terms inside
the brackets independent of k. Hence, when one evaluates the covariance,
00
R <£> - T exp(lk.e) S (k) dk (B.2)
cc •*• ' ~~ •*• cc — —
— oo
transforming to spherical coordinates and evaluating the integral
over wave number k, results in the following:
n 2- » 22
R <£> - / / l(B,u) /exp(ikcosfj) - ** •> o 9 •> dkd8d" (B'3)
cc 00 0 (l+k-A^)J B'+a-k-
where A, B depend on the angular coordinates from the used transformation
and 1(9 ,u) contains the remaining terms in Eq. (B.2).
Note that due to symmetry of Scc(k_) , i.e.,
Scc(k1 ,k2,k3) = Scc(-k1,-k2 ,-k3)
expanding the exponent, exp(iD) = cosD + isinD, in the integral
Eq. (B.2) cancels the imaginary part.
The spherical coordinate transformation is chosen so that k«£_ = k 5 cosg ,
i.e., depends only in two parameters of integration k,g. Appendix C describes
the details of the used transformation which is similar to the one used by
Bakr et al. (1978) to obtain analytical results for the head field covariance.
"f*—
Denote,
B =
-------
The integrals can be nondiraensionalized using the parameters
"1= ^ ' "2- TJ ' (p3= 1J« e = "I; ' K = ^ ' Ut = Mi (B-5)
Then converting to spherical coordinates and denoting
E2= e2[P2B2-hc(p2C2+D2)]2 (B.6a)
A*= -T (aij"juj]2+ -T l'yPjuj]2+ l-sj'j"/ (B'6b)
P! P2
[a ]. = Transformation matrix between the hydraulic conductivity
-1 principal axes and the axes aligned to the mean flow direction
C, = 5/i (B.6c)
(B.6d)
232
Y IT
the integrals become
2 2 v 2lF
R (£1X1111) = oc£0 f f PF. .(B>(i);xi'l')sinBM djjdu G.G. (B.7a)
cc f 3 0 0 1J 1 J
22
M = fcos(Ru) A " . , / , . du (B.7b)
0 (l+AZu )J pfB^VE-
Eq. (B.7) is the same as Eq. (5.33c) if one denotes the three-dimensional
integral as TIJ(.). The coefficient Fij(.) in the integrand is given
from:
22222 (B'8)
L = B^+C +0
The integral over u in Eq. (B.7b) can be estimated analytically.
From partial fraction expansion, for A2*0, E2#pjA2B2,
85
-------
k2A2
223 p2222
(l+kA) B+Ek
k2-
(k2-Pl)3
(B.9a)
(B'9b)
- P
- P
11 3
(prp2)3
' A
I3
» A
with the aid of the following integration formulas [Gradshteyn and Ryzhik,
1980, 3.737]
dx =
16
3 e
-labl
dx .
4 b
abe
-ab
0 (bN-x')J 2|b|
Eq. (7b) becomes:
Case 1: A2* 0, E2* A2B2Pl2 ,8*0.
E A
A ( + 3
A
exp •
(B-U)
86
-------
If B = 0, Eq. (B.6) becomes,
M- fe < + 3 +3) SXP I" ] (B>12)
E A
If E2 = p 2A2B2 , Eq. (B.6) becomes
» 2
M = /cos (kR) — J-r- k 4 dk
0 E-A* (-L. * kV
E
Finally, if A = 0,
M = 0 (B.14)
The remaining eq. (B.7a) Ls solved numerically using a Romberg integration
scheme for the integration over 3 and the trapezoidal rule for the resulting
one-dimensional integral.
87
-------
APPENDIX C
TRANSFORMATION TO POLAR COORDINATES
The coordinate system (x°, y°, z°) in Figure C-l represents the
original coordinate system to which (x1,x2,x3) or (k1,k2,k3) are aligned. The
objective of the transformation to the polar coordinates is to express the dot
product jt .§£ in the exponential in Eq. (5.12) in terms of k and 0, i.e., two
parameters only. The transformation will be performed in three steps.
Step 1 Create a Cartesian system of axes, (gi,e.l)i where e belongs to the
plane (xo,yo) and |j is on vector |.
The vectors fps.Jt can be expressed in terras of the angles Xi'|» (Figure C-l).
T
|j = [ cosy, sinxcosiji, sinxsm^] , where the superscript "T" stands for
transpose and gj has unit length.
E = [e , ev, 0], since it belongs in (xQ,yo). Invoking its orthogonality with
respect to £, E • 5 = 0, or,
e
n
, cosy, Oj
1
(L - sin
f
Finally, for j, from orthogonality
t • e = 0 and t • 1 = 0
and the following expression can be obtained:
I/0
sin
2 2
1 - sin xcos
(1 - sin
•*)
Note, as a verification, that when 5 || x , i.e., x = °>
1 0 0
0 1 0
0 0 1
[I.S.I] •
For notation purposes, denote,
x y z
" 8 8
•y v V
'x 'y 'z
88
-------
Step 2 Express £, in terms of polar coordinates (k,g,u) in the system (|,e,J)
ip i i i ip
it = [k cosg, k singcosu, k singsinu] = [k.,k2>k_]
where, for the integration purposes, B:0-nr, u:0->-2Ti, k:0+M.
Step 3 Express [kj1,k2',k3']T in terms of [kj.k2.k3]. This operation is
given from:
kl
k2
k3
=
ax ay az
Sx ey Sz
V V V
'x 'y 'z
V
k2'
k3'
It can be easily verified that jt • £ = kgcosg.
T T
JL I - *' I
0
0
kgcosf]
89
-------
Figure C.I Schematic graph of the coordinate transformation
geometry
90
-------
APPENDIX D
MEAN CONCENTRATIONS FITTED AT THE BORDEN SITE DATA
For Che longitudinal section the fitted curve had the following form,
c(x) - K • exp [ - (X " X) ] (D.I)
AV
Through a manual trial and error adjustment of the parameters, the values
selected for Figure 7.5 are:
K = 319 ag/ji
AL = 0.34 m
x = 17.83 m
For the vertical sections the fitted curves had the following form:
(z - Z)2
c(z) = X exp [ -
(38.3)
where 38.3 m = x, the longitudinal coordinate of the center of mass of the
plume, as given in Sudicky (1985). Through a similar manual trial and error
adjustment of the parameters, the values selected for Figures 7.6(a) through
7.6(1) are shown in Table D.I.
TABLE D.I
PARAMETERS USED TO FIT MEAN CONCENTRATION AT THE WELLS
x (m)
35
38.3
40.3
42.3
44.3
46.22
48.3
y (m)
12
16
17.94
20
22
24
25.95
Kv
------- |