COMPUTER SIMULATION OF MASS TRANSPORT IN GROUNDWATER:
AFFECT OF MACROSCOPIC HETEROGENEITIES IN HYDRAULIC CONDUCTIVITY
COMPUTER SIMULATION OF MASS TRANSPORT IN GROUNDWATER: AFFECT
OF MACROSCOPIC HETEROGENEITIES IN HYDRAULIC CONDUCTIVITY
BY
PETER B. MCMAHON, B.S.
THESIS
Presented to the Faculty of the Graduate School of
The University of Texas at Austin
in Partial Fulfillment
of the Requirements
for of
the Degree
MASTER OF ARTS
THE UNIVERSITY OF TEXAS AT AUSTIN
May, 1984
ACKNOWLEDGMENTS
I would like to thank my advisor Dr. John M. Sharp, Jr. and
committee members Dr. Randy Charbeneau and Dr. Charlie Kreitler for their
technical and grammatical review of this project. To Dr. Sharp for his
moral and monetary support, and Dr. Charbeneau for the use of the
Trescott model.
To all of the \onderful people associated with the University of
Texas Geology Department, I say thank you. There a lack of
was never
friendship, stimulating conversation, or armwaving on the fourth floor.
Special recognition belongs to Paul Blanchard, Pamela Tiezzi, and Jay
Vogt.
Finally, I would like to express special gratitude to my parents,
Rachel
and Jerry, for their constant love.
Submitted to Committee: 23 March 1984
III
ABSTRACT
In this study a computer model was used to simulate dissolved
chloride movement through alluvial sediments which border the Canadian
River in Hutchinson County, Texas. Hydraulic conductivity values of the
sediments were required in order to calculate groundwater velocities in
the system. The most realistic representation of conductivity variations
in porous media is expressed by frequency distributions rather than by
averaged values of conductivity. Numerous sedimentological environments
used
exhibit lognormal conductivity distributions; therefore, one was
in this investigation.
A number of conclusions can be based on the results of this study
First, certain conductivity distributions account for the observed spread
of chloride in the aquifer. The best match of'observed chloride disper
sion was obtained with autocorrelated lognormal conductivity distribu
tions. Secondly, the degree of spatial dependence between adjacent con
ductivity values affected numerous results. These include the amount of
chloride dispersion and the extent of uncertainty in calculated hydraulic
head and chloride distributions.
For comparative purposes the chloride distribution was also
modeled using an average conductivity value. Under this condition the
chloride plume moved at an average rate of 10 meters/year. Another
result was that longitudinal and transverse dispersivities of 46 meters
and 9 meters, respectively, were required to obtain a match between
observed and modeled chloride distributions.
IV
TABLE OF CONTENTS
ACKNOWLEDGMENTS iii
ABSTRACT iv
TABLE OF CONTENTS
v
LIST OF TABLES vi
LIST OF FIGURES vii
I. INTRODUCTION 1
11. OBJECTIVES 2
111. THEORY AND MECHANISM OF MASS TRANSPORT 2
IV. GEOLOGY 13
V. ASSUMPTIONS 23
VI. NUMERICAL METHODS 31
VII. RESULTS 41
A. Deterministic Model 41
B. Stochastic (Random) Model 49
C. Stochastic (Autocorrelated) Model 59
VIII. CONCLUSIONS 79
APPENDIX A 84
APPENDIX B 86
BIBLIOGRAPHY 107
VITA 11l
V
LIST OF TABLES
table page
I. Measured (Cl~) in mg/I for the period 198182. 46
11. The summations of velocities in the xand ydirections for the five hydraulic conductivity distributions used in this study. 80
111. Summary of the dispersivity versus alpha trends the five hydraulic conductivity distributions used in this study. for 83
VI
LIST OF FIGURES
figure page
1.
Cartoon illustrating scaledependent dispersion in a heterogeneous medium 10
2.
Map showing location of study site in Texas. 14
3.
of the Canadian River and its
Photograph floodplain.
Picture taken from highway 2277 bridge looking west. 15
4.
Map showing location of study site in Hutchinson County, Texas. 16
5.
Map showing structural features of Texas Panhandle (after Gustavson, 1980a). 17
6.
Stratigraphic columns of the major basins of the Texas 18
Panhandle (after Gustavson, 1980a).
7. Water table surface and locations of monitor wells at 22
study site (after Geogulf, 1981).
8. Cross section from the pond to the Canadian River (after 24
Geogulf, 1981).
9.
Photograph showing pond closure operations. 25
10.
Threedimensional pump test modeling results, delx and
dely equal 1000 feet. 29
11. Threedimensional pump test results, delx and dely equal 100 feet. 30
model. 33
12. Finite Difference grid used in the computer
Photograph of the Peimian redbeds south of highway 2277. 35
13.
14.
Photograph of alluvial sand deposit along the Canadian
an
River just east of the highway 2277 bridge. 36
was
15. Sketch showing how the leakage rate from the pond
determined. 37
VII
figure page
16. Flow chart illustrating sequence in 'which computer
programs were executed in order to produce a complete stochastic simulation. 42
17.
Head distribution, homogeneous K distribution. 43
18.
Chloride distribution, homogeneous K distribution 45
.
19. Photograph of the Canadian River from the highway 2277
bridge looking west. The photograph represents approximately a 30 foot cross section of the river. 51
20. Mean hydraulic head distribution for the random log
normal hydraulic conductivity distribution. 52
21. Standard deviation in hydraulic head for the random lognormal hydraulic conductivity
distribution. 54
22. Mean chloride distribution for the random logdistribution. 55
normal hydraulic conductivity
23.
Example of a random lognormal hydraulic conductivity distribution produced by the program MK. 57
24.
Standard deviation in chloride produced by the random lognormal hydraulic conductivity distribution.
58
25. Mean head distribution produced by the autocorrelated
lognormal hydraulic conductivity distribution, alpha =0.7. 60
26. Standard deviation in head produced by the autocor
related lognormal hydraulic conductivity distribution,
alpha=o.7. 62
27. Mean chloride distribution produced by the autocorre
lated lognormal hydraulic conductivity distribution,
alpha=o.7. 64
28.
Example of an autocorrelated lognormal hydraulic conductivity distribution produced by the program MK. 67
29.
Standard deviation in chloride produced by the auto
correlated lognormal hydraulic conductivity distribution, alpha=o.7. 69
figure page
30.
Mean hydraulic head distribution produced by the autocorrelated lognormal hydraulic conductivity distribution, alpha=0.3.
31.
Standard deviation is hydraulic head produced by the autocorrelated lognormal hydraulic conductivity distribution, alpha=0.3. 71
32.
Mean chloride distribution produced by the autocor
70
distribution, alpha=0.3. 72
related lognormal hydraulic conductivity
Standard deviation in chloride produced by the auto
33.
correlated lognormal hydraulic conductivity distri
bution, alpha=0.3. 73
34. Mean hydraulic head distribution produced by the
autocorrelated lognormal hydraulic conductivity distribution, alpha=0.5. 74
35.
Standard deviation in hydraulic head produced by the autocorrelated lognormal hydraulic conductivity distribution, alpha=0.5. 75
36.
Mean chloride distribution produced by the autocorrelated lognormal hydraulic conductivity
dis
tribution, alpha=0.5, 76
37. Standard deviation in chloride produced by the auto
correlated lognormal hydraulic conductivity distribution, alpha=0.5. 77
I. INTRODUCTION
The potential exists for large scale contamination of ground
water resources in the United States. In fact, nearly 1.7 trillion
gallons of waste are placed in the ground in the U.S. each year (U.S.
EPA, 1977). It would be desirable, therefore, to understand how mass is
transported in groundwater. One approach is computer models to
to use
simulate mass transport in groundwater. Unfortunately, currently applied
modeling techniques do not describe the actual processes of mass trans
port in porous media.
Most computer models are deterministic. A deterministic model
is one in which none of the input parameters are variable (Freeze, 1975).
For example, the hydraulic conductivity range of a porous medium would
be represented by an averaged hydraulic conductivity value. The purpose
of this study was to simulate the movement of dissolved chloride at a
site of groundwater contamination in North Texas by using a stochastic
computer model of mass transport in porous media.
A stochastic model is one in which the input parameters are
described statistically in terms of the mean and standard deviation of
the parameter's distribution in the system being modeled. In this
study hydraulic conductivity was described statistically and used as
input in the stochastic computer model. Furthermore, the hydraulic
conductivity was assumed to follow a lognormal distribution (smith
and Freeze, 1979).
1
II. OBJECTIVES
The objectives of this study were both theoretical and field
oriented.
The theoretical objectives were 1) to quantify the effect of
different lognormal hydraulic conductivity distributions on mean hy
draulic head and dissolved chloride distributions; 2) to examine the un
certainty in head and chloride distributions produced by different log
to
normal conductivity distributions; 3) study the effect of changing the
conductivity distributions on groundwater velocities.
field objectives were 1) to model the observed chloride
My distribution deterministically. In doing this the longitudinal and trans
verse dispersivities required to obtain a match between observed and
calculated chloride distributions would be determined. I also wanted to
model the system stochastically. If the stochastic model could produce
the same amount of mixing as the deterministic model, but with smaller
dispersivities, the mixing must be created by the largescale hydraulic conductivity variations created by the stochastic model.
Largescale conductivity variations in porous media are responsible for fieldscale
of mass
mixing in groundwater (Schwartz, 1977).
III. THEORY AND MECHANISMS OF MASS TRANSPORT
of
The success a computer model in predicting contaminant dis
tributions in groundwater systems is dependent upon the model’s ability
to accurately describe the dispersive character of the natural system.
Mass transport in saturated porous media occurs through the processes
of advection and hydrodynamic dispersion. Advection moves particles
directly with the average velocity of the flowing groundwater. Hydro
dynamic dispersion moves particles by mechanical mixing and molecular
diffusion. Molecular diffusion is the mechanism whereby molecules move
in response to a concentration gradient. Mechanical mixing is created by
microscopic variations in fluid velocity within the pore spaces of the
medium.
models of in media based
Computer mass transport porous are upon
partial differential equation describing
a mass The deter
transport.
ministic model requires using dispersivities in the differential equation
is of the medium and is of
The dispersivity a property porous a measure
the medium’s ability to transport and spread dissolved chemical species
in groundwater. Dispersivities can also be thought of as indicators of
and
the uncertainty in groundwater velocity in the system (Domenico
Palciauskas, 1979). Gillham and Cherry (1982) stated that the observed
degree of spreading of a contaminant in groundwater is generally much
greater than would be predicted by using laboratoryscale dispersivities.
Dispersivities determined in laboratory experiments are on the order of
centimeters while those used to describe fieldscale mixing are on the
order of meters (hang and Anderson, 1982). Fieldscale dispersivities
are determined from trialanderror approximations of field data.
The dispersion portion of the classical governing equation of
mass in which these dispersivities are used, is based on the
transport,
assumption that macroscopic dispersion is analogous to Pick's law of
diffusion. The problem lies in the fact that modeling largescale
dispersion as a diffusion phenomenon is commonly an incorrect approach
to a
different and considerably more complex process (Smith and Schwartz,
Research into the nature of
1980). mass transport in macroscopically
heterogeneous porous media is just beginning to examine this problem in
detail.
The velocity field which provides the advective motion of mass
transport is described mathematically by Darcy’s law and the groundwater
flow equation (Mercer and Faust, 1981).
=k /4 •vh (i)
v
vKvb + R (2)
=
Where, v is the average linear groundwater velocity, K is the hydraulic
conductivity tensor, h is the hydraulic head, $ is porosity, R is a
source/sink term (e.g. 'water wells, rivers, evaporation ponds, evapo
transpiration), and Ss is the specific storage of the medium (see Appendix
A for a complete listing of all symbols and their units used in this
paper). The derivation of these equations follows continuum approach
a
in that the actual ensemble of grains is replaced by a representative
volume
(Freeze and Cherry, 1979, p. 69).
In other words, the parameters and boundary conditions required
to solve these equations cannot be defined at every point in the system
so they must be averaged over a representative volume of the medium.
For example, a constant porosity cannot be defined for a volume of
medium that is smaller than its representative volume and
porous (Freeze
Cherry, 1979, p. 69). This introduces error in the calculations.
Hubbert (1940, p. 59) addressed the conditions under which these equations
provide reasonable results. In general, a good representation of the
bulk fluid motion be obtained for
can many systems of saturated porous
media. This is not the case for many types of fractured media.
Molecular diffusion is described by Pick’s law of diffusion.
F=~D vC (3)
Where, F is the mass flux, D is the coefficient of molecular diffusion,
and C is the concentration of dissolved chemical species.
Ideally, fluid velocities within the individual pores should be
determined in order for mechanical mixing to be accurately described.
As stated above, however, only average velocities can be determined
(Schwartz, 1977). Bear (1979, p. 232) showed that the spatial averaging
of the local, microscopic mass flux (Cv) due to velocity variations with
in pore spaces is:
Cy=CV Cv (4)
+
The first term on the right represents the advective flux, where
average
the average velocity is calculated with equation (1). The second term
fluctuation from Cv, the dispersive
represents the spatially averaged or,
flux. Mechanical mixing is represented by the second term in equation
(4). It is assumed that the dispersive flux is analogous to Pick’s law
of diffusion (Bear, 1979, p, 232). This is represented by equation (5).
(!v = "DvC (5)
The coefficient D, a second order tensor, is the coefficient of mechani
cal dispersion. Scheidegger (1961) proposed the following relationship
for the coefficient of mechanical dispersion.

D,j =**\vi 3 (18)
where, BELT is the time step used in the simulation and ANORM(O) is
a
number between 6 and +6, drawn from normal distribution of numbers
a
having a standard deviation of one and a mean of zero.
Equation (18)
states that dispersion can transport mass up to six standard deviations
from the center of mass of the chloride distribution.
away
The particleinacell method moves each particle a distance v*
BELT during every time step. The velocity in this equation is the aver
age linear velocity determined from equation (1). Throughout this study
porosity was assumed to be 30%. The total displacement for a given time
step for each particle is determined by:
v
+ o+al Mcnt = VtC bt iT +* 2dj*• (19)
c)isf»l ace
brer’) Co)
v
where, all terms have already been defined.
The method used to produce the autocorrelated lognormal conduc
tivity distributions was described by Smith and Freeze (1979). An
autocorrelated lognormal conductivity distribution is one in which the
hydraulic conductivity value at any point in a porous medium is
dependent upon adjacent conductivity values. In a system that exhibits
there is for
a large degree of autocorrelation, a tendency zones of equal
to be For with a small
conductivity grouped together. a system degree of
autocorrelation, zones of high conductivity are as likely to be situated
as
near low conductivity they are to other high conductivity zones. If
Kis lognormally distributed, Y=log(K) is normally distributed. Equation
(20) defines the autocorrelation technique used by Smith and Freeze
(1979) and this study.
=
Xj °