HIGH RESOLUTION CHARACTERIZATION OF AQUIFERS TO IMPROVE FLOW AND
TRANSPORT MODELS OF HIGHLY HETEROGENEOUS MEDIA
By
Mine Dogan Diker

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Geological Sciences - Doctor of Philosophy
2013

ABSTRACT
HIGH RESOLUTION CHARACTERIZATION OF AQUIFERS TO IMPROVE FLOW AND
TRANSPORT MODELS OF HIGHLY HETEROGENEOUS MEDIA
By
Mine Dogan Diker
Aquifers are the primary sources of clean drinking water. Pollution in aquifers is one of
the most challenging and important environmental problems. It is not only extremely complex to
map but also difficult to remediate. Flow and transport of water and pollutants in porous media
requires detailed characterization of the properties of the media. The main property which
controls the flow and transport is hydraulic conductivity (K), which can be defined as the ability
of the media to let the water flow through. Intensive studies to map the distribution of hydraulic
conductivity are necessary to model the plume migration. Conventional in-situ aquifer
characterization techniques are invasive and lack the necessary high resolution. Therefore, novel
methods are required to improve the methods to monitor and simulate the flow and transport
through aquifers.
This study introduces a combination of novel techniques to provide the necessary
information related to porous media. The proposed method was tested at a highly heterogeneous
site called the Macro Dispersion Experiment (MADE) site in Mississippi. The MADE site is a
very well studied site where several large scale tracer tests were conducted in the 1980s and
1990s. The tracers used for these tests were monitored using more than 300 multi-level sampler
(MLS) wells. Concentration measurements showed that the majority of the mass stayed near the
injection area, whereas minute concentrations were measured further down-gradient. This
behavior is significantly different from the simulations created using models based on the

Advection-Dispersion Equation (ADE). This behavior and the inability to explain this using most
models has led to a major debate in the hydrologic science community.
The hypothesis of this study is that the ADE based models can reproduce simulations of
the measured transport when the models are parameterized with sufficient high-resolution
hydraulic conductivity data. Two novel high resolution characterization methods, the direct-push
high resolution hydraulic conductivity (HRK) tool and 3D full-resolution ground penetrating
radar (GPR) were combined to generate 3D K fields using fractal stochastic methods. This study
demonstrated that the complementary geophysical data can be used to reduce the K variance by
dividing the aquifer into hydrofacies. This approach, in combination with a fractional
differencing filter, simplifies the statistically complex distribution of K. Fractional differencing
was also capable of removing the long range dependence in vertical K profiles to investigate the
underlying K distribution. The 3D K fields were then used to test the ADE based modeling
approach at the site and resulting concentrations were compared to one of the large scale tracer
experiments. The simulations in this study resulted in mass distributions comparable to those
measured during the tracer test experiments. They successfully reproduced the extent of the
plume in both 1D and 2D using K fields based on solely field data. Additional tests emphasized
the importance of high-resolution data to parameterize K models to successfully simulate flow
and transport using the ADE model.

Copyright by
MINE DOGAN DIKER
2013

For my parents, Emine and Sami Dogan, for all the love, trust, and support they provided.
I am proud for being able to make them proud of me.

v

ACKNOWLEDGEMENTS

First and foremost, I would like to state that there are no words to express my gratitude
for my advisor Dr. Remke L. van Dam. He supported me, encouraged me, challenged me and,
most importantly, he inspired me. Together, we collected millions of GPR traces during many
cold and humid nights in Mississippi. He became not only an advisor, but also a colleague and a
most supportive friend throughout this process. His expertise as an earth scientist guided me
towards becoming a better scientist with an optimistic character. Without him, my experience as
a graduate student would not be as bright and satisfying. “Beep, evet!”
Second, I would like to acknowledge my co-advisor, Dr. David W. Hyndman, for being a
great role model for my academic personality and guiding me through the great deal of stress that
every graduate student experiences. His confidence in me made every task seem doable and
every target reachable. It was a great gift and privilege having not only one, but two wonderful
advisors.
I would like to thank Dr. Mark M. Meerschaert for his endless patience for my questions,
explanations of challenging statistical concepts, and belief that I could make this project work.
His creative insight about our work was exciting and gave me clarity about how to proceed. I
also would like acknowledge Drs. Brian Hampton, Kamini Singha, and David Benson for their
contribution to my research. They were the ones who made my research complete and strong.
I would like to thank the Hydrogeology Lab members, especially my friends, Erin King
and Alex Kuhl, who have been very encouraging and supportive throughout this process. Gizem
and Mustafa Ali, for listening to my endless explanations and practice talks , just their existence
made many things easier than they actually were. And, my graduate school sister, Dr. Nicole D.

vi

LaDue, for being the magnificent best friend, guiding me through the cultural differences, and
challenges. Thanks lady!
Many thanks dear Jackie, Jennae, and Heidi of our departmental office for all of their
help handling paperwork and managing deadlines. Their smiling faces always made me thankful
to be in the Department of Geological Sciences.
I would like to acknowledge Dr. Gulcin Ozurlan Agacgozgu, who is a role model for me
as an intelligent, and confident female scientist, and Dr. Argun Kocaoglu since they kept
supporting me all the way from Turkey.
My husband, Kaya, thank you for joining me on this immense journey, and following me
all the way to Michigan. Your existence made East Lansing a real home for me and you have
helped me grow stronger from this experience.
Finally, my family, I love you to the pieces. My father, Sami Dogan, has been always the
greatest role model and my biggest supporter. I could not be who I am without him. Thanks,
Dad! My mother, Emine Dogan, made the biggest sacrifice of accepting my decision to live far
away, always missing her only daughter. I want to thank to my brother, Malik, for leading and
encouraging me to venture out of my safe shell and become a world citizen. Wish you see the
better ones from your sons, my beloved nephews, Saman Eren and Bilge Cinar.

vii

TABLE OF CONTENTS
LIST OF TABLES .......................................................................................................................... x
LIST OF FIGURES ....................................................................................................................... xi
Chapter 1 Prologue ......................................................................................................................... 1
1.1 Introduction ................................................................................................................. 1
1.2 Aquifer Characterization ............................................................................................. 3
1.3 The MADE Site ........................................................................................................... 6
1.4 Methods ..................................................................................................................... 10
1.5 Outline of the thesis/dissertation ............................................................................... 14
REFERENCES ................................................................................................................ 18
Chapter 2 Hydrostratigraphic analysis of the MADE site with full-resolution GPR and directpush hydraulic profiling ................................................................................................................ 21
Abstract............................................................................................................................ 21
2.1 Introduction ............................................................................................................... 22
2.2 Methods ..................................................................................................................... 24
2.3 Results ....................................................................................................................... 27
2.4 Conclusions ............................................................................................................... 31
2.5 Acknowledgments ..................................................................................................... 32
APPENDIX ..................................................................................................................... 34
REFERENCES ................................................................................................................ 41
Chapter 3 Hydraulic conductivity fields: Gaussian or not? .......................................................... 43
Abstract............................................................................................................................ 43
3.1 Introduction ............................................................................................................... 44
3.2 Statistical analysis ..................................................................................................... 48
3.3 Model validation ........................................................................................................ 56
3.4 Discussion.................................................................................................................. 58
3.5 Conclusions ............................................................................................................... 59
3.6 Acknowledgments ..................................................................................................... 60
APPENDICES ................................................................................................................. 61
APPENDIX 3.1: Fractional difference filter ............................................................. 62
APPENDIX 3.2: Flow-Transport Simulations .......................................................... 65
REFERENCES ................................................................................................................ 71
Chapter 4 Novel characterization method provides a major advance for flow and transport
prediction ...................................................................................................................................... 76
Abstract............................................................................................................................ 76
4.1 Introduction ............................................................................................................... 77
4.2 Methods ..................................................................................................................... 78
4.3 Results ....................................................................................................................... 81
4.4 Discussion and Conclusions ...................................................................................... 83
REFERENCES ................................................................................................................ 86

viii

Chapter 5 Quantifying the value of different data sets and modeling schemes for flow and
transport simulations ..................................................................................................................... 89
5.1 Introduction ............................................................................................................... 89
5.2 Complementary GPR data ......................................................................................... 92
5.3 Cost-benefit analysis ................................................................................................. 94
5.3.1 Transient flow ................................................................................................... 95
5.3.2 Flowmeter K ..................................................................................................... 98
5.3.3 Effect of conditioning data density................................................................... 99
5.3.4 Complementary GPR data .............................................................................. 103
5.4 Conclusions ............................................................................................................. 104
REFERENCES .............................................................................................................. 107
Chapter 6 Synthesis..................................................................................................................... 109
REFERENCES .............................................................................................................. 115

ix

LIST OF TABLES
Table 1.1 Inventory of collected GPR data. ................................................................................. 14
Table A2.1 Single-sample K-S tests of four HRK profiles for facies and layers identified in fullresolution 3D GPR data (see Figures 2.1 and 2.3a for locations). Columns, “SEGMENT", name
of segment (facies/sub-facies/layer), "N_SAMPLE", number of HRK measurements used,
"p_FOR_LOG_K", p values for log10 K using single-sample K-S test with 95% CI. ................. 36
Table A2.1 (cont'd)....................................................................................................................... 37
Table A2.2 Statistical comparison of K values for four identified radar facies and adjacent
segments, where p values were calculated using two-sample K-S tests. Columns, "SEGMENT",
names of compared consecutive segments (facies/sub-facies/layers), "p_FOR_LOG_K", p values
for log10 K using two-sample K-S test with 95% CI. ................................................................... 37
Table A2.2 (cont'd)....................................................................................................................... 38
Table 5.1 Quantitative measures for the comparison of simulation scenarios and field experiment
2
(R values were calculated comparative to the MADE-1 experiment t=503 days in the first row).
....................................................................................................................................................... 95

x

LIST OF FIGURES
Figure 1.1 Number of articles over the years which comprise "aquifer characterization" as
keywords. For interpretation of the references to color in this and all other figures, the reader is
referred to the electronic version of this dissertation. ..................................................................... 3
Figure 1.2 Vertically exaggerated (x5) hydraulic conductivity cross-section along A-A' using
1560 flowmeter K measurements [Boggs, 1991] (see Figure 1.4 for location). ............................. 6
Figure 1.3 Multi level sampler and flowmeter well setup at the MADE site (photo courtesy:
stanford.edu). .................................................................................................................................. 8
Figure 1.4 Site map showing the locations for flowmeter K measurements with cross-section line
(A-A’ in Figure 1.1), multilevel sampler wells, and test boundary. .............................................. 9
Figure 1.5 Panoramic view of the site showing the source location (metal pipes in the centerfront) and GPR data collection setup (on the right) (photo courtesy: Kaya Diker). ..................... 12
Figure 1.6 Site map showing HRK and GPR measurement locations, and the original MADE
test boundary. ................................................................................................................................ 12
Figure 1.7 Side view of DP High Resolution K (HRK) probe (photo courtesy: Kaya Diker). .... 13
Figure 2.1 (a) Map of the MADE site on Columbus Air Force Base (AFB) with GPR
measurement lines. The GPR survey coordinates are shown in blue. The blowup of the ICA cube
2
(144 m ) shows DP sites and locations of the 3D GPR cubes in Figures 2.1b and c (yellow
shaded area) and the profile in Figure 2.2 (red dashed line); viewing angles are indicated with
arrows in corresponding colors. Full-resolution 3D GPR data cubes at (b) 100 MHz and (c) 50
MHz are shown with no vertical exaggeration. An envelope was used to render negative
amplitudes transparent; number labels in Figure 2.1b are discussed in the text. .......................... 25
Figure 2.2 Interpretation of GPR and HRK data at line 105E (see Figure 2.1 for location). (a)
GPR profile with red (positive) to blue (negative) amplitude scale using combined 100 MHz and
50 MHz data; black triangles indicate the zone where the two data sets were averaged. In
addition to processing mentioned in the text, these data were plotted with an energy decay gain.
(b) Continuous reflections identified using an automated picking algorithm and interpreted GPR
facies (color shaded). (c) Qualitative interpretation of GPR facies with HRK data; facies
boundaries are marked by horizontal lines. .................................................................................. 28
Figure 2.3 (a) 3D GPR facies boundaries (shaded to visualize topography) with collocated DP
HRK profiles. Descriptive statistics of these K data for (b) all saturated material above the
aquitard, (c) GPR facies, (d) sub-facies, and (e) layers. Box plots show the sample median,
interquartile range, and positions of extreme values. (f) Variance of log10 K and ln K values for
the data in Figures 2.3b–2.3e, respectively. The red line is a power fit through the medians of the
variance values for each group (horizontal axis in log scale). ...................................................... 31

xi

Figure A2.4 (a) Cross plot of log10 K mean and variance for high resolution K data from four
DP profiles (see Figure 2.1 for location and Figure 2.2 for data). (b) Cross plot of log10 K mean
and the slopes of linear trends of the change in log10 K with depth. In these plots, each facies is
shown by a different symbol and color. The number next to each symbol indicates the DP profile
(1, 2, 3, and 4 represent 111108A, 111108B, 111108C, and 121108A, respectively). ................ 39
Figure 3.1 Layout of Macro Dispersion Experiment test site, showing key features of MADE
experiments, as well as the locations of GPR data collected for this project. The inset of the
12×12 m ICA (Intensively Cored Area) cube shows the locations of the four HRK profiles and
the 2D transect discussed in this paper. ........................................................................................ 47
Figure 3.2 Histogram of ln K for HRK profile 121108A (see map in Figure 3.1) before (a) and
after (b) applying the fractional difference filter (1) with d=0.9. The filtered data are organized
into a unimodal distribution with a sharper peak and a heavier tail than the best fitting Gaussian
probability density function (black line). ...................................................................................... 49
Figure 3.3 Autocorrelation function for ln K from HRK profile 121108A (see Figure 3.1) before
(a) and after (b) applying the fractional difference filter (3.1) with d=0.9. Autocorrelations inside
dashed lines are statistically zero. ................................................................................................. 51
Figure 3.4 Simulated ln K field without (a) and with (b) GPR facies (dashed lines), conditioned
on four HRK profiles (vertical black lines). Histogram (c, d) of one column (white line at x = 172
m) from simulated ln K field (a, b, respectively) after applying the fractional difference filter (1)
with d=0.9. The histogram (c) fits a Gaussian model, but the histogram (d) from facies
simulation (b) deviates from Gaussian shape, similar to measured HRK data (Figure 3.2b)....... 53
Figure 3.5 Fractionally differenced ln K data (a) from the shallowest facies at horizontal location
174 m fits a Gaussian distribution. Deepest facies (b) at horizontal location 170 m deviates from
the Gaussian model. These probability plots show the sorted data on the horizontal axis, and the
corresponding model percentiles for the best fitting Gaussian model on the vertical axis. If the
data fits this model, the points will follow the reference line, with some random scatter. ........... 55
Figure A3.6 Simulated plume 500 days after the end of injection, using a constant K value in
each GPR facies (top), the K field from Figure 3.4a without GPR facies (middle), and the K field
from Figure 3.4b with GPR facies (bottom). ................................................................................ 67
Figure A3.7 Normalized concentration profiles (top) 500 days after the end of injection, and
normalized concentration breakthrough curves (bottom), measured at location x=172 m, from the
three simulations illustrated in Figure A3.6. ................................................................................. 68
Figure 4.1 (a) Map of the MADE site with the test boundary (dashed line) and sampling
locations; gray shaded rectangle shows the model domain used for simulations. (b) Model
domain with HRK, injection, and observation borehole locations. .............................................. 80
Figure 4.2 Relative mass distribution profiles 501 days after the injection. The black line shows
the mean of the simulations along with 2σ error range (in gray shaded area); the red line shows
the profile for the MADE-1 observations. .................................................................................... 81

xii

Figure 4.3 Vertically integrated contour maps of relative mass 503 days after the start of
injection for MADE-1 experiment (a), the mean of the 6 simulations that met the head criteria (b)
a sample simulation which matches the head change criteria best (c). All three contour maps
were created following the same interpolation procedure, as discussed in the text. .................... 82
Figure 5.1 A detailed map of the model domain with locations of flowmeter K (orange) and
HRK measurements (blue) and multi-level sampler (grey) wells. The area of 3D GPR data is
given by the rectangle with green dashed line. Injection wells are given with red stars. See Figure
4.1 for the location within the larger MADE site. ........................................................................ 91
Figure 5.2 Vertically integrated contour maps of relative mass; (a) MADE-1 experiment at
t=503 days, (b) mean of 6 simulations, (c) best head change simulation. K fields for the
simulations in (b) and (c) were based on a combination of HRK data and GPR facies (these maps
can be compared with Figure 4.2, which presents contour maps of relative mass for K fields
based on just HRK data). .............................................................................................................. 92
Figure 5.3 Relative mass distribution profiles for the MADE-1 experiment at t=503 days (red
line), and the mean and 2σ range of 6 simulations (black line and gray shaded area, respectively).
Simulations used (a) K fields based on just HRK data and (b) K fields based on a combination of
HRK data and GPR facies............................................................................................................. 93
Figure 5.4 Vertically integrated contour maps of normalized mass; (a) MADE-1 experiment at
t=503 days, (b) base simulation, (c) transient flow simulation. .................................................... 97
Figure 5.5 Normalized mass distribution profiles of the MADE-1 experiment at t=503 days (red
line), base simulation (black line), and transient flow scheme (blue line). .................................. 97
Figure 5.6 Vertically integrated contour maps of normalized mass; (a) MADE-1 experiment at
t=503 days, (b) base simulation, (c) simulation based on just flowmeter K data. ........................ 99
Figure 5.7 Normalized mass distribution profiles of the MADE-1 experiment at t=503 days (red
line), base simulation (black line), and simulation based on just flowmeter K data (blue line). .. 99
Figure 5.8 Vertically integrated contour maps of normalized mass: (a) MADE-1 experiment at
t=503 days, (b) base simulation based on all 25 HRK profiles, and simulations based on (c) 20 ,
(d) 15, (e) 10, and (f) 5 randomly selected HRK profiles. The randomly selected HRK profiles
used to generate the K fields are shown with blue symbols; see Figure 5.1 for the location of all
25................................................................................................................................................. 101
Figure 5.9 Normalized mass distribution profiles of the MADE-1 experiment at t=503 days (red
line), base simulation (black line), and simulations based randomly selected number of HRK
profiles (blue line). (a) 80% of HRK profiles (#’s 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 14, 15, 16, 17, 18,
19, 20, 22, 23, 24), (b) 60% of HRK profiles (#’s 2, 3, 4, 8, 10, 11, 12, 13, 14, 16, 17, 18, 20, 21,
24), (c) 40% of HRK profiles (#’s 2, 5, 7, 11, 12, 14, 17, 18, 19, 22) (d) 20% of HRK profiles
(#’s 8, 11, 14, 15, 24). See Figure 5.1 for the locations of the HRK profiles. ............................ 102

xiii

Figure 5.10 Vertically integrated contour maps of normalized mass: (a) MADE-1 experiment at
t=503 days, (b) base simulation, and (c) simulation based on a combination of HRK data and
GPR facies. ................................................................................................................................. 103
Figure 5.11 Normalized mass distribution profiles of the MADE-1 experiment at t=503 days
(red line), base simulation (black line), and simulation based on a combination of HRK data and
GPR facies. ................................................................................................................................. 104

xiv

Chapter 1
Prologue
1.1 Introduction
It is an extensive challenge to understand all aspects of underground water such as
storage, flow, transport, and interactions with surrounding material, since aquifers constitute the
most important natural medium for storage, purification and filtering of fresh water. Hydraulic
conductivity (K) is the main property that controls water flow in porous media [e.g., Dagan,
1989] K can simply be defined as the ability of any porous medium to allow the fluids to flow
through. It depends on the pore and grain size, their distribution, and sorting. To solve
environmental problems, accurately parameterized models that describe the flow of water and
solutes in aquifers are critical, since models that can consistently and realistically simulate flow
and transport rely on accurate representations of K. Although there have been many studies on
these topics involving in-situ measurements of aquifer properties and pollution monitoring , the
methods utilized in these studies have certain limitations and are mostly invasive. Furthermore,
the lack of an in-situ and non-invasive method to directly measure K is a major obstacle to create
realistic simulations of flow and transport through aquifers.
Even though several characterization methods exist, resolution still remains as an issue.
Characterization methods can be combined with Gaussian geostatistical methods to develop 3D
K fields to overcome the resolution problem, which at sites with low levels of heterogeneity
provides acceptable results [Freyberg, 1986; Mackay et al., 1986, Garabedian et al., 1991;
LeBlanc et al., 1991, Hyndman et al., 2000; Phanikumar et al., 2005]. Parameterization of K
fields for flow and transport modeling in highly heterogeneous aquifers, however, remains

1

challenging. Various stochastic and geostatistical methods that have been proposed to improve
simulations of highly heterogeneous K fields generally rely on conventional in-situ measurement
methods that provide limited spatial sampling, and have large support volumes with unknown
geometries. Therefore, innovative approaches to measure and interpret K are necessary to truly
understand flow and transport through highly heterogeneous media.
The main goal of this study is to develop a novel approach for predicting flow and
transport in highly heterogeneous porous media with the aid of novel measurement and
stochastic methods. If successful, this approach will provide a significant contribution to
ameliorate real-world problems related to ground water and aquifers. The principal aspect of this
research involves high-resolution characterization of hydraulic conductivity and aquifer
structure. Descriptive statistics and distribution tests were used to confirm the consistency of the
data sets used in this study. This comparative approach forms the foundation of this research,
since there is no known direct relationship between geophysical parameters and K.
Consecutively; geophysics-derived information was integrated with in-situ K measurements
using stochastic methods. Fractal methods, which are capable of representing the long-memory
nature and connectivity of K, were used to populate simulated 3D K fields to run flow and
transport models.
In this thesis, detailed K data from a heterogeneous site were used to model fluid flow
and solute transport utilizing the advection-dispersion equation. A combination of novel and
traditional measurements and data processing methods was explored for detailed characterization
of shallow, unconfined aquifers.

2

1.2 Aquifer Characterization
Aquifer characterization is a broad topic addressed by many researchers from different
disciplines, including earth science, hydrology, environmental engineering, and statistics.
Researchers from these fields have been employing many different tools including laboratory
measurements [Illman et al., 2010], in-situ measurements [Kabala, 1993], macro-scale field
experiments [Garabedian et al., 1991; Boggs, 1991], geophysical imaging [Dafflon et al., 2011],
and computer simulations [Gómez-Hernádez and Wen, 1998]. Figure 1.1 shows the numbers of
articles in Google Scholar for the keyword "aquifer characterization”. In recent years, the
number of the articles has increased significantly, and technological advances promote more
interest in this matter. Aquifer characterization methods, involving the measurement and/or

1000

806
567

500
9

13

26
1985-1988

1989-1992

99

63

1981-1984

325

1900-1980

167
2009-2012

2005-2008

2001-2004

1997-2000

0
1993-1996

Number of Articles

estimation of K, will be discussed from here on.

Years
Figure 1.1 Number of articles over the years which comprise "aquifer characterization" as
keywords. For interpretation of the references to color in this and all other figures, the reader is
referred to the electronic version of this dissertation.
There are several different methods to determine hydraulic conductivity. Empirical
approaches involving pedotransfer functions rely on textural analysis for defining the physical

3

properties of the medium such as pore/grains/fracture size, distribution, and organization.
Textural analysis involves collecting core samples from the aquifer and sieving these samples to
get grain size distributions [Masch and Denny, 1966]. Although several methods, such as
freezing the samples, exist to collect nearly undisturbed cores, it is still an extremely challenging
process especially for unconsolidated sediments.
Alternative experimental methods to determine K require either laboratory measurements
(following Darcy's law) on the core samples or more complicated procedures for in-situ
experiments. In both cases, the measured values are affected by the core removal or
experimentation procedure since re-sorting due to the vibration and the loss of pore water are
almost unavoidable. Flowmeter tests are the most commonly used in-situ experimentation
methods since they allow in-situ measurements of bulk parameter K. These are capable of
measuring K averaged over an uncertain volume and have been widely used to map the hydraulic
properties of aquifers. Although several existing studies demonstrate the success of flowmeter
tests in mildly heterogeneous media, they also have many limitations. Beside the invasiveness,
flowmeter tests are relatively slow, and capable to measure the bulk parameter for an unknown
volume. Therefore, these tests can provide K information with a limited resolution and, are not
very sensitive to low K.
Indirect measurement procedures also exist, and include slug tests, infiltration and tracer
experiments. These methods intend to derive K distribution using other measured parameters
rather than attempting to measure the actual K. Infiltration experiments based on the amount of
water infiltrated through a limited surface area and only capable of providing the average K
estimates for limited and generally unknown depth. Slug tests require addition or removal of a
known volume of water from a designated well, with detailed measurements of the associated

4

changes in the water levels. Tracer experiments intend to derive a [or the] K distribution by
measuring concentration distribution of the injected tracer. They are capable of providing the
largest scale K distribution and may require a large amount of sampling.
Several other indirect measurements in the boreholes include well logging methods, that
do not require the material removal. However, they require complementary data and an
interpretation step since a known relationship between K and the measured parameters (i.e.,
resistivity, seismic velocity, gamma ray exposure) does not exist.
Geophysical methods are another set of techniques which can only provide indirect
information related to K distribution. However, they often provide the advantage of noninvasiveness and higher resolution than the borehole based methods. The potential benefits of
geophysical methods in addition to, or as replacements of, current in-situ methods are very
significant. Geophysical methods can be effectively used for measuring bulk properties and
mapping certain structures such as fractures, gravel beds, and clay lenses, which may have strong
effects on the movement of groundwater either as preferential flow paths or as retardant barriers
[Overmeeren, 1998; Streich et al., 2006]. Innovative approaches to collect and interpret
geophysical data to map these structures, can also contribute to the efforts to derive a quantitative
affinity between geophysics-derived parameters and hydraulic conductivity. Besides,
improvements of geophysical methods provide data for the non-invasive monitoring of pollution
and remediation processes [Brewster and Annan, 1994; Halihan et al., 2005].
In recent years, several novel borehole based approaches have emerged to delineate
aquifer characteristics at higher resolution than previously possible. These approaches include insitu measurements using direct-push technology, which is a significantly faster alternative to the
conventional borehole methods, especially in unconsolidated sediments. A recently developed

5

direct-push tool, called High Resolution K (HRK) tool, which couples direct-push permeameter
(DPP) and direct-push injection logger (DPIL), provides cm-scale vertical resolution and partly
overcomes an important obstacle related to resolution issues [Liu et al., 2009]. This relatively
less destructive tool involves a probe with an injection port and two pressure sensors. Water is
injected through the port and the difference in back pressure is measured via pressure sensors.
These backpressure measurements are then converted to relative K values. Surface-based
geophysical methods can be used to provide complementary data to direct-push methods which
give high resolution data in vertical profiles.
1.3 The MADE Site
The method developed in the context of this study was applied at the Macro Dispersion
Experiment (MADE) site on the Columbus Air Force Base of Mississippi. The MADE site is a
highly heterogeneous, unconfined shallow aquifer. Previous research suggests that the sediments
consist of meandering fluvial deposits over braided fluvial deposits over a fine-grained sand
layer inter-bedded with clay and silt [Bowling et al., 2005].

Elevation [m]

0

Distance [m]
100
150

50

200

250

62
60
58
56
54
52
-4

-2

0

2

log10 (K [m/day])

Figure 1.2 Vertically exaggerated (x5) hydraulic conductivity cross-section along A-A' using
1560 flowmeter K measurements [Boggs, 1991] (see Figure 1.4 for location).
The MADE site is one of the most heavily studied uncontaminated sites, and it represents
a high level of heterogeneity in many contaminated sites [Zheng et al., 2010]. It is a secure

6

facility that has been used to explore flow regimes in an aquifer using macro scale tracer tests
and borehole measurements. The research at the site was started with measurements of ground
water flow in 1984; later, the MADE-1 experiment was conducted between October 1986 and
June 1988 [Boggs, 1991]. It was an extensive natural gradient experiment involving the
preliminary flowmeter measurements of K. Figure 1.1 shows a cross section created using these

K measurements along a transect (see Figure 1.4 for location), and the variance of ln K was
calculated as 4.5. The experiment was monitored using 328 multi-level sampling wells (Figure
1.3) [Zheng et al., 2010; Boggs et al., 1992]. A mixture of several different chemicals, including
bromide (in the form of calcium bromide, Cinj=2500 mg/l), pentafluorobenoic acid (PFBA,
Cinj=400 mg/l), trifluoromethylbenzoic acid (TFBA, Cinj=400 mg/l), and orthofluorobenzoic
3

acid (DFDA, Cinj=400 mg/l) was used as a tracer. The total mixture volume was 10.03 m and
groundwater extracted from the aquifer was used to prepare the mixture. The injection setup
consists of five boreholes installed on a linear array with approximately 1 m separation (Figure
1.4) approximately perpendicular to the average down gradient direction. Injection was carried
out through a 0.6 m screened interval between the 57.5 and 58.1 m elevations above sea level
(7.4-8 m below the surface) during 48.5 hours.
Sampling was done on a snapshots basis, and the sampled area was increased through
time for each snapshot. A total of 11446 samples were collected for 8 snapshots, and analyzed
for the concentrations of the injected chemicals. Bromide (total injected mass of 25.0 kg) was the
primary tracer due to its highly conservative nature and consistency of measured concentrations.
Analysis of the measured concentrations suggested rather an interesting outcome, extremely
different than what was expected. More than 20 % of the injected mass stayed around the

7

injection area and did not move any further than 10 meters even 500 days after the injection. In
contrast, very low but detectable concentrations moved extremely fast down gradient. This
interesting behavior resulted in a concentration profile with several peaks and a heavy tail which
did not match the expected Gaussian shape (as predicted by advection-dispersion equation
models). Thereafter, the flow behavior at the MADE site became a notorious problem for
hydrologists. Two more macro-scale (MADE-2 and -3) and several smaller scale tracer tests
were conducted to further investigate the issues at the site.

Figure 1.3 Multi level sampler and flowmeter well setup at the MADE site (photo courtesy:
stanford.edu).
The MADE-2 experiment was conducted between June 1990 and September 1991
[Boggs, 1993]. The same injection procedure was followed but a different mixture was used as
tracer including tritium (Cinj=55,610 pCi/ml),

14

C labeled p-xylene (Cinj=2770 pCi/ml), benzene

(Cinj=68.1 mg/l), p-xylene (Cinj=51.5 mg/l), naphthalene (Cinj=7.23 mg/l), and o3

dichlorobenzene (Cinj=32.8 mg/l). The total injected mixture volume was 9.7 m and observed

8

using 258 multi-level sampling wells. The sampling procedure for the second experiment was
improved with additional multi-level sampling wells installed along two transects known as
fences, in the vicinity of injection area. However, the recovered mass was still not comparable to
the injected mass and the concentration profiles presented the same interesting behavior.

300

A′
250

Distance [m]

200

150

100

N
50

0

MADE
test boundary

Injection wells
MLS wells
Flowmeter K wells
0

A

50
100
Distance [m]

150

Figure 1.4 Site map showing the locations for flowmeter K measurements with cross-section line
(A-A’ in Figure 1.1), multilevel sampler wells, and test boundary.
Another attempt was made with the MADE-3 experiment. This time, the injection
procedure was passive, involving a trench filled with sand and tracer mixture. The natural flow

9

was allowed to dissolve the tracer and transport [Boggs et al., 1995]. The natural condition of the
sediments in the injection area was destroyed and homogenized irreversibly with this
experiment. Several smaller scale tracer tests, including a push-pull test and an extensive core
sampling procedure following blue dye injection [Liu et al., 2010], were also carried on recently.
Several approaches have been explored to simulate the flow regime at the MADE site
following the experiments. Most of the earlier research concluded that the advection-dispersion
equation (ADE) based approaches were not capable of reproducing the observed transport
behavior at the MADE site [Zheng et al., 2010]. Modeling efforts based on ADE resulted in
smooth concentration curves that were not comparable to the experiment outcomes. More
exhaustive approaches followed the classic ADE method, including preferential flow paths
approach, mass transfer approach between mobile and immobile domains, utilizing the fractal
version of the ADE [Benson et al., 2001; Zheng and Gorelick, 2003; Harvey and Gorelick,
2000]. However, these approaches either did not include adequate parameterization due to lack
of necessary data or involved naturally undeterminable features and calibration procedures. In
conclusion, a successful modeling effort based on solely field data and capable of reproducing
the observed transport behavior is still not available.

1.4 Methods
Geophysical methods such as electrical resistivity tomography (ERT) and ground
penetrating radar (GPR) are minimally-invasive and can often provide high resolution data. They
have been used for mapping hydrologically different structures and more recently for tracking
fluid migration during tracer tests [Johnson et al., 2007]. Borehole and surface GPR have been
successfully used to map sedimentary structures [Neal, 2004], estimate porosity from GPR
velocity fields [Klotzsche et al., 2010], track Dense Non-Aqueous Phase Liquid (DNAPL) and

10

saline tracers [Brewster and Annan, 1994; Birken and Versteeg, 2000; Hwang et al., 2008], and
monitor remediation [Halihan et al., 2005]. However, most of these studies lack comparatively
high-resolution in-situ measurements of K, which in the absence of a direct relationship is
necessary to determine the connection between GPR signal response and K.
GPR reflections originate at interfaces between geological materials with distinctive
values of dielectric permittivity, which is directly related to water content, which for saturated
media is governed by porosity. The correlation between hydraulic conductivity (which is related
to porosity) and dielectric properties of media has been discussed in recent publications [Chen et

al., 2001; Lambot et al., 2006; Kowalsky et al., 2005; Klotzsche et al., 2010], but no direct
relationship has been established. Morin [2006] found an inverse relationship, which contrasts
the common belief that K is higher for large porosities. Therefore, in this research, GPR data was
used to delineate main unit boundaries and different sedimentary facies.
GPR and direct-push (DP) hydraulic conductivity data collection was completed during
two field campaigns in 2008 and 2009 (Figure 1.5 and see Figure 1.6 for the locations). Multiple
GPR data collection techniques were applied, with most collection using: 2D profiles, fullresolution 3D and 4D data cubes, and cross-borehole tomography (Figure 1.6). GPR data along
2D profiles were collected for broader scale site characterization and three distinct regions were
selected for full resolution GPR data cubes for detailed characterization (Table 1.1). Regions to
collect full resolution 3D GPR data were selected from previous tracer test sites to take
advantage of existing grain size and tracer test data.

11

Figure 1.5 Panoramic view of the site showing the source location (metal pipes in the centerfront) and GPR data collection setup (on the right) (photo courtesy: Kaya Diker).

300

250

ICA cube

Distance [m]

200

150
MLS cube
DTA cube
100
N
STA cube
MADE
test boundary
50
HRK wells
Injection wells
2D GPR profiles
3D GPR cubes
0
0
50
100
Distance [m]

150

Figure 1.6 Site map showing HRK and GPR measurement locations, and the original MADE
test boundary.

12

58 high resolution vertical K profiles were obtained with the new DP High Resolution K
(HRK) probe (Figure 1.7). This tool, which was developed for rapid characterization of
unconsolidated shallow aquifers [Liu et al., 2009], is advanced into the subsurface while water is
injected out of a small screened port located a short distance behind the tool tip. The injection
rate and injection induced back pressure are recorded every 1.5 cm. The ratio of these quantities
is transformed into K following the approach described by Liu et al. [2009]. Although the
calibration of the transform equation is the subject of ongoing work, the spatial patterns of K,
which are of greatest interest in this study, would not change with different transformation
parameters.

Figure 1.7 Side view of DP High Resolution K (HRK) probe (photo courtesy: Kaya Diker).

13

Table 1.1 Inventory of collected GPR data.
Area
MADE site
from tracer tests
Intensively Cored Area
(ICA)

Size

Data Type

250x150 m 2D profiles
12x12 m

Multi-Level Sampler Area
(MLS)

25x25 m

Dipole Tracer Test Area
(DTA)

5x12 m

Source Trench Area
(STA)

30x10 m

Frequencies
50-100 MHz

3D data cubes
Vertical Radar Profiling
3D data cubes
Zero Offset Profiling
Multi Offset Gathering
4D data cubes
Zero Offset Profiling
Multi Offset Gathering

50-100-200-250 MHz
50-100 MHz
50-100 MHz
50-100 MHz
50-100 MHz
50-100 MHz
50-100 MHz
50-100 MHz

3D data cubes

50-100 MHz

1.5 Outline of the thesis/dissertation
This study aims to provide a novel approach to create flow-transport models solely based
on field data to reproduce the outcomes of MADE-1 and -2 experiments. This introductory
chapter provides the history and the technical details related to the experiments conducted at the
MADE site. Below, I present a brief synopsis of the following chapters in this dissertation.
Chapter 2 presents a qualitative comparison of GPR and K data on which the remaining
chapters are based. This chapter includes a statistical analysis of K data as a proof of concept.
The objective of this second chapter is to investigate the potential of the geophysical data to
delineate the hydrostratigraphically different units and to demonstrate the correlation between
the GPR reflections and K data. Full-resolution 3D GPR derived facies boundaries were used to
define the boundaries of different hydrostratigraphic structures. Descriptive statistics and
distribution tests were then used to examine the power of subdivisions to reduce the variances of

K distributions. This chapter shows that GPR reflections coincide with anomalies in K profiles

14

and separate units with distinct K characteristics. The chapter presents a new approach for
building high-resolution 3D hydrostratigraphic frameworks of heterogeneous aquifers and to
improve the stochastic K field creation procedure to better simulate transport through such highly
heterogeneous aquifers.
The proof of concept is followed by Chapter 3, which introduces the quantitative and
statistical analysis of HRK data along with the fractional differencing approach. The MADE site
sediments possess a very complex and widely varied K distribution. Henceforth, conventional
statistical approaches cannot be successfully utilized to model the K distribution. In this Chapter,
I show that fractional differencing filter serves to remove the long range dependence and reveals
the underlying K distribution. Fractal stochastic methods, which are not commonly used for
aquifer characterization, were then employed to simulate K fields containing the naturally
observed connectivity of the sedimentary structures. This chapter introduces the procedure
followed to simulate fractal stochastic K fields and demonstrates the value of fractal methods as
well as GPR derived hydrostratigraphy. It concludes that the use of fractal methods in
combination with GPR derived facies boundaries, provides improved representations of K
created using exclusively data based mixing of Gaussian fields. Chapter 3 is followed by an
Appendix that presents a comparative analysis of conditional and unconditional K fields with and
the effects of the hydrostratigraphic information obtained from GPR. This comparison includes
flow and transport simulations to explore the best suitable method for stochastic K field
generation.
Chapter 4 presents the outcomes of 3D flow transport simulations for an area (25x45 m)
near the MADE1-3 injection site. This chapter builds on the statistical methods developed in
Chapter 3, but are now applied in 3D. The research in this chapter is unique since it provides the

15

first 2D representations of simulated plumes in comparison to measured concentrations from
MADE-1 experiment. For this chapter both data sets published after the experiments were
revisited and reanalyzed. Resulting simulations were used to create 2D vertically integrated
contour maps and 1D profiles of relative mass distribution. Amount of recovered mass is not
comparable to the experimentation due to the limitations of reproducing the sampling procedure
in the field. However flow and transport simulations, based on the presented approach, are
significantly capable of representing the spatial extent of the observed plume as well as the
multi-peak and heavy tail behavior of the observed relative mass distribution along the down
gradient direction.
Chapter 5 combines the 3D flow simulations near the injection area with the GPR
hydrofacies approach developed in Chapter 2 and investigates the effect of this additional piece
of information on the outcomes of flow-transport simulations. This chapter demonstrates that the
complementary GPR data improves the representation of tail behavior. The second part of this
chapter includes a comparative analysis for simulated plumes for several different scenarios
including: use of flowmeter K data, a transient flow scheme, and changing densities of HRK
profiles. This part of the study showed that: (1) flowmeter K data alone is not sufficient to create

K fields to reproduce the plume behavior successfully, (2) no significant difference exist between
steady state and transient flow schemes in simulated plume shape and characteristics, (3) HRK
profiles have a significant effect on the simulated plume.
Finally, Chapter 6 lays out a synthesis of the previous chapters and tabulates the
conclusions of this research.

16

REFERENCES

17

REFERENCES
Benson, D. A., R. Schumer, M. M. Meerschaert, & S. W Wheatcraft (2001), Fractional
dispersion , Lévy motion , and the MADE tracer tests. Transport in Porous Media, 42(1908),
211-240.
Boggs, J. M. (1993), Database for the second Macrodispersion Experiment (MADE-2). EPRI
Topical Rep. TR-10 (Vol. 2072).
Boggs, J. M., S. C.Young, L. M. Beard, L. W. Gelhar, K. R. Rehfeldt & E. E. Adams (1992),
Field study of dispersion in a heterogeneous aquifer: 1. Overview and site description. Water
Resources Research, 28(12), 3281-3291. doi:10.1029/92WR01756.
Boggs, J. M. (1991), Database for the first Macrodispersion Experiment (MADE-1).
Boggs, J.M., J. A. Schroeder, and S. C.Young (1995), Data to support model development for
natural attenuation study. Report No. WR28-2-520-197. TVA Engineering Laboratory,
Tennessee Valley Authority, Norris, Tennessee.
Bowling, J. C., A. B. Rodriguez, D. L. Harry, & C. Zheng (2005), Delineating alluvial aquifer
heterogeneity using resistivity and GPR data. Ground Water, 43(6), 890-903.
Brewster, M. L., & A. P. Annan (1994), Ground-penetrating radar monitoring of a controlled
DNAPL release: 200 MHz radar. Geophysics, 59, 1211-1221.
Birken, R., & R. Versteeg (2000), Use of four-dimensional ground penetrating radar and
advanced visualization methods to determine subsurface fluid migration. Journal of Applied
Geophysics, 43(2-4), 215-226. doi:10.1016/S0926-9851(99)00060-9.
Chen, J., S. Hubbard, & Y. Rubin (2001), Estimating the hydraulic conductivity at the south
oyster site from geophysical tomographic data using Bayesian Techniques based on the
normal linear regression model. Water Resources Research, 37(6), 1603.
doi:10.1029/2000WR900392.
Dafflon, B., J. Irving, and W. Barrash (2011), Inversion of multiple intersecting high-resolution
crosshole GPR profiles for hydrological characterization at the Boise Hydrogeophysical
Research Site. Journal of Applied Geophysics, 73.4: 305-314.
Dagan, G. (1989), Flow and transport in porous formations. Springer-Verlag GmbH & Co. KG.
Garabedian, S. P., D. R. LeBlanc, L. W. Gelhar, and M. A. Celia (1991), Large-scale natural
gradient tracer test in sand and gravel, Cape Cod, Massachusetts: 2. Analysis of spatial
moments for a nonreactive tracer. Water Resources Research, 27: 911–924.
Gómez-Hernádez, J. J., and X. -H. Wen (1998), To be or not to be multi-Gaussian? A reflection
on stochastic hydrogeology, Adv. Water Resour. 21(1), 47–61.

18

Freyberg, D. L. (1986), A natural gradient experiment on solute transport in a sand aquifer: 2.
Spatial moments and the advection and dispersion of nonreactive tracers. Water Resources
Research, 22: 2031–2046.
Halihan, T., S. Paxton, I. Graham, T. Fenstemaker, & M. Riley (2005), Post-remediation
evaluation of a LNAPL site using electrical resistivity imaging. Journal of Environmental
Monitoring : JEM, 7(4), 283-287. doi:10.1039/b416484a.
Harvey, C., & S. M. Gorelick (2000), Rate-limited mass transfer or macrodispersion: which
dominates plume evolution at the Macrodispersion Experiment (MADE) site? Water
Resources Research, 36(3), 637–650. American Geophysical Union.
Hyndman, D. W., M. J. Dybas, L. Forney, R. Heine, T. Mayottee, M. S. Phanikumar, G. Tatara,
J. Tiedje, T. Voice, R. Wallace, D. Wiggert, X. Zhao, and C. S. Criddle (2000), Hydraulic
characterization and design of a full-scale biocurtain. Ground Water, 38: 462–474.
Hwang, Y. K., A. Endres, S. D. Piggott, & B. L. Parker (2008), Long-term ground penetrating
radar monitoring of a small volume DNAPL release in a natural groundwater flow field.
Journal of Contaminant Hydrology, 97(1-2), 1-12. doi:10.1016/j.jconhyd.2007.11.004.
Illman, W. A., J. Zhu, A. J. Craig, and D. Yin (2010), Comparison of aquifer characterization
approaches through steady state groundwater model validation: A controlled laboratory
sandbox study, Water Resources Research, 46, W04502, doi:10.1029/2009WR007745.
Johnson, T. C., P. S. Routh, W. Barrash, & M. D. Knoll (2007), A field comparison of Fresnel
zone and ray-based GPR attenuation-difference tomography for time-lapse imaging of
electrically anomalous tracer or contaminant plumes. Geophysics, 72(2), G21-G29.
doi:10.1190/1.2431638.
Kabala, Z. J. (1993), The dipole flow test: A new single-borehole test for aquifer
characterization, Water Resources Research, 29(1), 99–107, doi:10.1029/92WR01820.
Klotzsche, A., J. van der Kruk, G. Meles, J. Doetsch, H. Maurer, & N. Linde (2010), Fullwaveform inversion of crosshole ground penetrating radar data to characterize a gravel
aquifer close to the Thur River, Switzerland. Ground Penetrating Radar (GPR), 2010 13th
International Conference on (pp. 1–5). IEEE.
Kowalsky, M. B., S. Finsterle, J. Peterson, S. S. Hubbard, Y. Rubin, E. Majer, A. Ward, et al.
(2005), Estimation of field-scale soil hydraulic and dielectric parameters through joint
inversion of GPR and hydrological data. Water Resources Research, 41(11), 1-19.
doi:10.1029/2005WR004237.
Lambot, S., E. C. Slob, M. Vanclooster, & H. Vereecken (2006), Closed loop GPR data
inversion for soil hydraulic and electric property determination. Geophysical Research
Letters, 33(21), 1-5. doi:10.1029/2006GL027906.
LeBlanc, D. R., S. P.Garabedian, K. M.Hess, L. W.Gelhar, R. D.Quadri, K. G.Stollenwerk, and
W. W. Wood (1991), Large-scale natural gradient tracer test in sand and gravel, Cape Cod,

19

Massachusetts: 1. Experimental design and observed tracer movement. Water Resources
Research, 27: 895–910.
Liu, G., J. J. Butler Jr, G. C. Bohling, E. Reboulet, S. Knobbe, & D. W. Hyndman (2009), A new
method for high-resolution characterization of hydraulic conductivity. Water Resources
Research, 45(8), W08202. American Geophysical Union.
Liu, G., C. Zheng, G. R. Tick, J. J. Butler Jr., and S. M. Gorelick (2010), Relative importance of
dispersion and rate-limited mass transfer in highly heterogeneous porous media: analysis of a
new tracer test at the Macrodispersion Experiment (MADE) site. Water Resources Research
46: W03524. doi:10.1029/2009WR008430.
Mackay, D. M., D. L. Freyberg, P. V. Roberts, and J. A. Cherry (1986), A natural gradient
experiment on solute transport in a sand aquifer: 1. Approach and overview of plume
movement. Water Resources Research, 22: 2017–2029.
Masch, Frank D., and Denny, Kleber J. (1966), Grain size distribution and its effect on the
permeability of unconsolidated sands. Water Resources Research, 2.4: 665-677.
Morin, R. (2006), Negative correlation between porosity and hydraulic conductivity in sand-andgravel aquifers at Cape Cod, Massachusetts, USA. Journal of Hydrology, 316(1-4), 43-52.
doi:10.1016/j.jhydrol.2005.04.013.
Neal, A. (2004), Ground-penetrating radar and its use in sedimentology: principles, problems and
progress. Earth-Science Reviews, 66(3-4), 261-330. doi:10.1016/j.earscirev.2004.01.004.
Van Overmeeren, R. A. (1998), Radar facies of unconsolidated sediments in The Netherlands: A
radar stratigraphy interpretation method for hydrogeology.Journal of Applied Geophysics,
40.1: 1-18.
Phanikumar, M. S., D. W. Hyndman, X. Zhao, and M. Dybas (2005), A three-dimensional model
of microbial transport and biodegradation at the Schoolcraft, Michigan site. Water Resources
Research, 41: W05011.
Streich, R., J. Van Der Kruk, and A. G. Green (2006), Three-dimensional multicomponent
georadar imaging of sedimentary structures. Near Surface Geophysics, 4.1: 39-48.
Zheng, C., M. Bianchi, & S. M. Gorelick (2010), Lessons Learned from 25 Years of Research at
the MADE Site. Ground Water, 1-14. doi:10.1111/j.1745-6584.2010.00753.x.
Zheng, C., & S. M. Gorelick (2003), Analysis of solute transport in flow fields influenced by
preferential flowpaths at the decimeter scale. Ground Water, 41(2), 142-155.

20

Chapter 2
Hydrostratigraphic analysis of the MADE site with full-resolution GPR and direct-push
hydraulic profiling

Abstract
Full-resolution 3D Ground-Penetrating Radar (GPR) data were combined with high-resolution
hydraulic conductivity (K) data from vertical Direct-Push (DP) profiles to characterize a portion
of the highly heterogeneous MAcro Dispersion Experiment (MADE) site. This is an important
first step to better understand the influence of aquifer heterogeneities on observed anomalous
transport. Statistical evaluation of DP data indicates non-normal distributions that have much
higher similarity within each GPR facies than between facies. The analysis of GPR and DP data
provides high-resolution estimates of the 3D geometry of hydrostratigraphic zones, which can
then be populated with stochastic K fields. The lack of such estimates has been a significant
limitation for testing and parameterizing a range of novel transport theories at sites where the
traditional advection-dispersion model has proven inadequate.1

This chapter is based on Dogan, M., R. L. Van Dam, G. C. Bohling, J. J. Butler Jr., and D. W.
Hyndman (2011), Hydrostratigraphic analysis of the MADE site with full-resolution GPR and
direct‐push hydraulic profiling, Geophysical Research Letters, 38, L06405,
doi:10.1029/2010GL046439.

21

2.1 Introduction
The transport of solutes through an aquifer is primarily controlled by medium properties,
in particular the spatial distribution of hydraulic conductivity (K) [e.g., Gelhar, 1993;

Fleckenstein and Fogg, 2008]. Studies of mildly heterogeneous aquifers have demonstrated that
solute transport can be reasonably modeled using the classical advection-dispersion equation
(ADE) with limited K data [e.g., Mackay et al., 1986; Hess et al., 1992], and that
hydrostratigraphic analysis of core material improves transport predictions [e.g., Phanikumar et

al., 2005]. In contrast, studies in highly heterogeneous aquifers have shown that the classic ADEbased approach with K data from conventional field methods does not accurately simulate
transport in such systems [e.g., Eggleston and Rojstaczer, 1998; Whittaker and Teutsch, 1999].
Indeed, three large-scale natural gradient tracer experiments performed at the MAcro Dispersion
Experiment (MADE) site (Figure 2.1) on Columbus Air Force Base, Mississippi, USA, displayed
pronounced non-Gaussian behavior [Boggs et al., 1992; Zheng, 2006]. The MADE aquifer
consists of highly heterogeneous unconsolidated fluvial sediments (ln K variance = 4.5 from
borehole flowmeter data [Rehfeldt et al., 1992]), underlain by a clay aquitard at ~12 m depth.

Bowling et al. [2005] used 2D Ground-Penetrating Radar (GPR) lines and information from a
nearby quarry to identify three main facies above the aquitard: a meandering fluvial system over
a braided fluvial deposit over a fine-grained sand interbedded with clay and silt.
Several approaches have been proposed for simulating the observed tracer transport at
MADE by incorporating preferential flow paths [e.g., Zheng and Gorelick, 2003] or mass
transfer between mobile and immobile domains [e.g., Harvey and Gorelick, 2000] into the ADE,
or by using a fractional form of the ADE [e.g., Benson et al., 2001]. Although these approaches
may provide reasonable representations of the average plume behavior, they do not accurately

22

replicate concentration histories at observation wells nor can they be parameterized using
available data. Novel high-resolution characterization methods, however, may provide the
necessary subsurface data to greatly improve traditional transport simulations in such highly
heterogeneous systems and aid in the assessment of alternative transport theories.
GPR is a common noninvasive method for high-resolution exploration of spatial
variability in the shallow subsurface [e.g., Jol, 2009], but it does not provide direct information
about K [Hubbard and Rubin, 2000]. Previous efforts to use GPR to improve flow and transport
models have combined GPR facies analysis with modeled K fields and stochastic simulations
[Rauber et al., 1998; Moysey et al., 2003; Ezzy et al., 2006; Engdahl et al., 2010]. Such studies
have not directly combined surface 3D GPR data with high-resolution in-situ K estimates to
develop hydrofacies models for heterogeneous aquifers, which is the focus of this paper.
In this paper, we present results of a recent field demonstration at the MADE site where
full-resolution 3D GPR and cm-scale Direct-Push (DP) K data were collected. Following a
description of the approach and methods, we discuss the general reflection patterns in the GPR
data cubes. We then present the results of a GPR facies analysis for a 2D plane where four DP K
profiles were obtained. Following a qualitative comparison of these collocated data sets, we
present the results of statistical tests to evaluate whether GPR facies are also distinct
hydrostratigraphic units. The results of this field demonstration indicate that the combination of
methods presented here is a promising approach for characterizing 3D hydrostratigraphic
structures. These structures, which can then be populated by stochastic simulation of K fields,
can serve as the basis for flow and transport models of highly heterogeneous aquifers, such as at
the MADE site.

23

2.2 Methods
GPR is an excellent method to image shallow sedimentary structures because the signal
response is controlled by textural properties [Neal, 2004]. In saturated low-loss media, such as
sand and gravel, the variable most directly linked to GPR signal propagation and reflection is
porosity, which is governed by sediment characteristics such as grain size, sorting, and packing.
GPR has traditionally been used for 2D and pseudo-3D characterization, but recent studies have
demonstrated the added value of full-resolution data, with less than quarter wavelength (λ)
sampling for in- and cross-line directions [Grasmueck et al., 2005]. Full-resolution GPR
maximizes the potential to characterize 3D subsurface structures. Its vertical resolution depends
on signal wavelengths, which depend on by frequencies of propagating waves and dielectric
permittivities of the medium. For example, vertical resolution is ~0.145 m (1/4 λ) for 100 MHz
signals in saturated sediments with a relative dielectric permittivity of 23 (EM velocity ~0.058
m/ns). The lateral resolution depends on the Fresnel zone, which gets larger with increasing
depth and decreasing frequency.
We used 2D GPR lines to characterize the stratigraphy over the region where three
natural-gradient tracer experiments [Zheng, 2006] were conducted. We then collected fullresolution 3D GPR data around the Intensively Cored Area (ICA – Figure 2.1a) where a singlewell, push-pull tracer test was recently performed [Liu et al., 2010]. A total of 3.8 km of GPR
lines were collected in the ICA cube using 50 and 100 MHz antennae, with step sizes of 0.2 and
0.1 m, respectively, which is less than the 1/4 λ required spacing. Line spacing was equal to step
size, thus forming a regular grid of GPR traces. Data were collected using a sampling interval of
800 ps over 550 and 400 ns time windows, and 16 and 32 stacks for the two frequencies,
respectively. Accurate positioning was achieved using guidance ropes and odometer-wheel

24

Figure 2.1 (a) Map of the MADE site on Columbus Air Force Base (AFB) with GPR
measurement lines. The GPR survey coordinates are shown in blue. The blowup of the ICA cube
2
(144 m ) shows DP sites and locations of the 3D GPR cubes in Figures 2.1b and c (yellow
shaded area) and the profile in Figure 2.2 (red dashed line); viewing angles are indicated with
arrows in corresponding colors. Full-resolution 3D GPR data cubes at (b) 100 MHz and (c) 50
MHz are shown with no vertical exaggeration. An envelope was used to render negative
amplitudes transparent; number labels in Figure 2.1b are discussed in the text.

25

triggering. All GPR data were collected with a Sensors and Software pulseEKKO 100 system
(1000V transmitter) at night and on weekends to avoid flight-time interference with a
communication station adjacent to the site (Figure 2.1a).
GPR data were processed with a background removal filter (dewow, 14 and 8 ns for 50
and 100 MHz, respectively) followed by a band-pass filter to eliminate high-frequency noise.
Static corrections were then applied to flatten the reflection from the top of the saturated zone, as
the measured water table gradient was only ~0.0003 (3.3 cm over 111 m). The reflection times
for the saturated zone were converted to depths based on the average measured velocity of 0.058
m/ns, from CMP and cross-borehole data.
A GPR facies approach [Van Overmeeren, 1998] was used to identify zones with distinct
reflection characteristics. The primary criteria used to define GPR facies were reflection
terminations, dip angle, amplitude, and continuity in 3D. We then compared GPR facies with
high-resolution vertical K profiles that were obtained with the new DP High-Resolution K
(HRK) probe. This tool, which was developed for rapid characterization of unconsolidated
shallow aquifers [Liu et al., 2009], is advanced into the subsurface while water is injected out of
a small screened port located a short distance behind the tool tip. The injection rate and injectioninduced back pressure are recorded every 1.5 cm. The ratio of these quantities is transformed
into K following the approach described in Liu et al. [2009]. Although the calibration of the
transform equation is the subject of ongoing work, the spatial patterns of K, which are of greatest
interest in this study, would not change with different transformation parameters. We used the
Kolmogorov-Smirnov (K-S) test and box plots to evaluate differences among K distributions for
different GPR facies and layers.

26

2.3 Results
Cutouts of the 100 and 50 MHz GPR data, with 10 and 12 m of signal penetration,
respectively, clearly image the details of 3D structures (Figures 2.1b and c). The 100 MHz northsouth oriented cut at 97 m East shows two ~2 m thick packages with northward dipping
reflections between 4 and 8 m depth (

and

in Figure 2.1b). These structures likely

represent large-scale clinoforms associated with channel bar migration. This interpretation is
corroborated by a reflection pattern along the perpendicular cut at 170 m North that resembles
trough cross-stratification (

in Figure 2.1b). The GPR reflections from the deepest portion of

the cube are dominated by sub-horizontal continuous reflectors, but the signal is notably
attenuated for the 100 MHz data. The 50 MHz data, which depict the same dipping clinoforms,
have reasonable signal strength to the top of the clay aquitard (Figure 2.1c).
We conducted facies analysis across a transect at 105 m East (Figure 2.2), where the
general reflection pattern is comparable to the plane at 97 m East. For this analysis, 100 MHz
data were used from the water table to 8 m, and 50 MHz data were used below 9 m. The average
data from both frequencies were used between 8 and 9 m depth (Figure 2.2a), since picks from
both were consistent. To define the spatial distribution of GPR stratigraphy, we developed an
algorithm for automated picking of peak amplitudes and identification of laterally continuous
reflections (Figure 2.2b). Decisions on how reflections connect and terminate were aided by 3D
analysis of the data. Using the procedure outlined earlier, the GPR data were separated into four
GPR facies (Figure 2.2b). Facies A (green) consists of sub-horizontal reflections, and can be
divided into two sub-facies those appear to be separated by an erosional surface. The underlying
Facies B (brown) contains the most notable clinoform sets, and these can again be divided into

27

two sub-facies. The lowermost two facies are characterized by laterally continuous subhorizontal reflections. Facies C (blue) has several internal clinoform structures (based on 3D
analysis of the GPR data) and gently dipping bounding surfaces; Facies D (tan) has primarily
horizontal reflections.

Figure 2.2 Interpretation of GPR and HRK data at line 105E (see Figure 2.1 for location). (a)
GPR profile with red (positive) to blue (negative) amplitude scale using combined 100 MHz and
50 MHz data; black triangles indicate the zone where the two data sets were averaged. In
addition to processing mentioned in the text, these data were plotted with an energy decay gain.
(b) Continuous reflections identified using an automated picking algorithm and interpreted GPR
facies (color shaded). (c) Qualitative interpretation of GPR facies with HRK data; facies
boundaries are marked by horizontal lines.

28

Figure 2.2c presents a qualitative comparison of the GPR data and the HRK profiles
along this transect. The GPR data show good correlation with the HRK profiles, as is evident
from numerous small-scale anomalies that coincide with GPR reflections. In addition, GPR
facies are consistent with the main K zones across the interpreted portion of the aquifer. Facies A
has generally high K values, consistent with coarse-grained sediments. The relatively constant K
values in Facies A in profiles 111108A-C reflect an upper HRK measurement threshold of
roughly 10 m/d; the actual K values are likely higher and more variable than indicated. In the
zone with prominent clinoforms, Facies B shows declining K with depth. Facies C shows
constant to increasing conductivities with depth, whereas Facies D is characterized by generally
high K values. There is a clear transition into the low-K aquitard at the bottom of HRK log
111108A (Figure 2.2c); in other logs, DP probe advancement was halted at the top of the clay.
The observations from this qualitative evaluation suggest that GPR and HRK methods
can be used in tandem for high-resolution hydrostratigraphic analysis (Figure 2.3a). To evaluate
this possibility quantitatively, we statistically analyzed the K data within GPR facies, sub-facies,
and layers. Boxplots in Figures 2.3b-e visualize the descriptive statistics of the K data for GPR
facies, sub-facies, and layers (collectively called 'segments'), respectively. One-sample K-S tests
with 95% confidence intervals (CI) rejected the null-hypothesis that the K data from each
segment have a log-normal distribution (see Table A2.1 of the appendix); therefore, the common
assumption that K distributions are log-normal is not valid for these data.
To evaluate differences in K distributions, we used two-sample K-S tests with 95% CI.
These tests show that the distribution of K data for each of the four GPR facies is distinct (see
Table A2.2 of the Appendix). Similar K-S tests were used to test the difference of K distributions
between adjacent segments. These tests show that all adjacent layers as well as sub-facies/facies

29

have statistically distinct distributions of K (see Table A2.2 of the Appendix). Although this
analysis suggests that the individual GPR facies and layers can be translated into statistically
distinct hydrofacies, it does not indicate that both are necessarily equivalent. When the K data are
separated by individual HRK profiles in a plot of mean versus variance, the between-profile
variation in log K means for each facies is generally significantly smaller than the between-facies
variation (see Figure A2.4a of the Appendix). There is, however, considerable overlap between
the mean log K values for Facies B and C, which is expected since they have opposite trends of

K with depth. Indeed, a plot of mean log K versus the slope of log K values with depth in each
facies clearly separates the facies into clusters (see Figure A2.4b of the Appendix). Figure A2.4a
also shows that the variance of log K is low for Facies A (affected by K truncation discussed
earlier) and D but higher for Facies B and C, which we argue is related to the depositional
environment.
Modeling of flow and transport through heterogeneous aquifers would greatly benefit
from detailed characterization of K. Figure 2.3f demonstrates that as the aquifer is split into
facies, sub-facies, and layers based on our stratigraphic analysis of full-resolution GPR data, the
total variance in K is drastically reduced. Most of this reduction occurs in the first two splits into
facies and sub-facies. Therefore, subdivision into layers may not be required to develop realistic
3D K fields at this site. A wide range of stochastic methods can be used to distribute the K data
through facies/sub-facies shown in Figure 2.3a. Although a stochastic K field could be developed
for all DP data shown in Figure 2.3a, it would clearly not be possible to fully capture the
geometry of the GPR stratigraphy based on the K data alone. The value of separating the cube
into facies, sub-facies or layers (Figure 2.3f) can then be quantified using simulations of tracer
tests through the stochastically derived K fields.

30

Figure 2.3 (a) 3D GPR facies boundaries (shaded to visualize topography) with collocated DP
HRK profiles. Descriptive statistics of these K data for (b) all saturated material above the
aquitard, (c) GPR facies, (d) sub-facies, and (e) layers. Box plots show the sample median,
interquartile range, and positions of extreme values. (f) Variance of log10 K and ln K values for
the data in Figures 2.3b–2.3e, respectively. The red line is a power fit through the medians of the
variance values for each group (horizontal axis in log scale).
2.4 Conclusions
Accurate predictions of transport through highly heterogeneous aquifers would greatly
benefit from a method to characterize the detailed structure of aquifers; this would be an
important first step to populate 3D K fields with high vertical and horizontal resolution. Recently
developed DP methods can provide high-resolution K (HRK) data in vertical profiles, yet they
cannot provide sufficient spatial density to establish lateral connectivity. In this paper, we
present the first comparison of full-resolution 3D GPR and HRK data to develop high-resolution

31

hydrofacies for highly heterogeneous sediments. Four GPR facies that were identified at the
MADE site were determined to be distinct hydrofacies based on statistical analysis of several
collocated HRK profile data. The division of these facies into smaller segments (sub-facies,
layers) results in zones with lower variance in K. These zones can then be used to generate
stochastic fields with less uncertainty than previously possible.
We found good agreement between full-resolution GPR stratigraphy and HRK profiles,
thus forming a solid foundation for hydrostratigraphic characterization of this site. Our approach
provides an opportunity to reconstruct 3D subsurface structures with their correct geometries and
hydrologic attributes. The strong connection between the HRK data and GPR facies indicates
that at many sites, a 3D K field could be generated using GPR data tied to a few HRK profiles.
Clearly, as facies change character laterally, additional HRK profiles are needed to intercept
these facies and to capture a representative distribution of K values.
It should be noted that the approach presented in this paper is not without limitations. The
vertical and lateral resolution of GPR data is finite, and reduces with increasing depth. Limiting
the variance of K within hydrofacies units and layers nevertheless provides the basis for betterconstrained stochastic fields for solute transport simulations. In future research, the hydrofacies
model described here will be used to test a range of emerging transport theories.

2.5 Acknowledgments
This work was supported by grants from NSF (EAR-0738938 and EAR-0738955); any
opinions, findings, and conclusions or recommendations expressed are those of the authors and
do not necessarily reflect the views of the NSF. Acknowledgment is made to the Donors of the
American Chemical Society Petroleum Research Fund for partial support of this research (PRF48515-G8). We thank G. Tick, K. Diker, E. Reboulet, S. Knobbe, U. Schneidewind, T.

32

Viencken, and K. Singha for field support, C. Zheng, P. Dietrich, G. Liu, and M. Meerschaert for
insightful discussions, K. Hubbard and Geoprobe Systems for equipment loans, and S. Reed and
Columbus AFB personnel for logistical support. Anonymous reviewers are acknowledged for
their constructive comments and suggestions.

33

APPENDIX

34

Stratigraphic analysis of full-resolution ground-penetrating radar (GPR) data at the
MAcro Dispersion Experiment (MADE) aquifer in Mississippi identified 4 radar facies, 2 of
which can be subdivided in sub-facies, and 30 individual layers (collectively called 'segments').
In two tables, we present statistical analyses of hydraulic conductivity (K) data obtained with the
HRK tool for these facies and layers. For each facies- and layer boundary in these tables, data
from a transition zone with 1/4 wavelength of the GPR frequency (1/8 wavelength on top and
bottom) was excluded from the analysis.
A one-sample Kolmogorov-Smirnov (K-S) test was used to check if K data in any of the
segments is log-normally distributed. Probability values larger than 0.05 indicate log normal
distribution for a 95% confidence interval. Combining K data from four HRK profiles (111108A,
111108B, 111108C, 121108A), Table A2.1 provides the number of data points and p values of
the entire aquifer, the four facies (A: layers 1-9, B: layers 10-22, C: layers 23-26, D: layers 2730), the sub-facies, and individual layers (see Figure 2.2 for location relative to GPR fullresolution data).
Table A2.2 provides the p values calculated using two-sample K-S tests with 95%
confidence intervals for the statistical comparison of facies and consecutive sub-facies/layers,
respectively. The null hypothesis of this test is that data sets of 2 samples come from the same
continuous distribution.
Separated by facies and DP profile, Figure 2.4 presents cross plots of log10 K mean
versus (a) variance and (b) the slopes of linear trends of the change in log10 K with depth.

35

Table A2.1 Single-sample K-S tests of four HRK profiles for facies and layers identified in fullresolution 3D GPR data (see Figures 2.1 and 2.3a for locations). Columns, “SEGMENT", name
of segment (facies/sub-facies/layer), "N_SAMPLE", number of HRK measurements used,
"p_FOR_LOG_K", p values for log10 K using single-sample K-S test with 95% CI.
SEGMENT N_SAMPLE p_FOR_LOG_K
ALL_DATA
2163
3.24E-28
Facies_A
503
6.47E-134
Facies_B
533
3.76E-33
Facies_C
302
1.06E-136
Facies_D
655
1.06E-31
Facies_A1
337
9.70E-109
Facies_A2
128
8.73E-26
Facies_B1
265
1.11E-03
Facies_B2
232
1.04E-41
Facies_C
302
1.06E-136
Facies_D
655
1.06E-31
Layer_01
55
9.08E-30
Layer_02
24
1.21E-07
Layer_03
26
1.84E-18
Layer_04
22
5.39E-11
Layer_05
65
8.40E-21
Layer_06
43
1.36E-29
Layer_07
13
3.29E-09
Layer_08
38
1.40E-11
Layer_09
48
3.26E-08
Layer_10
22
1.40E-04
Layer_11
33
2.68E-02
Layer_12
12
1.44E-02
Layer_13
38
1.85E-04
Layer_14
14
1.97E-10
Layer_15
21
5.78E-08
Layer_16
33
2.86E-08
Layer_17
15
7.84E-14
Layer_18
12
1.41E-09
Layer_19
18
5.94E-08
Layer_20
36
2.06E-13
Layer_21
12
1.13E-02
Layer_22
60
8.76E-20
Layer_23
10
1.13E-07
Layer_24
22
1.12E-16
Layer_25
92
6.56E-60

36

Table A2.1 (cont'd)
SEGMENT N_SAMPLE p_FOR_LOG_K
Layer_26
64
1.16E-17
Layer_27
79
6.41E-11
Layer_28
70
1.50E-08
Layer_29
123
4.47E-12
Layer_30
155
1.94E-18

Table A2.2 Statistical comparison of K values for four identified radar facies and adjacent
segments, where p values were calculated using two-sample K-S tests. Columns, "SEGMENT",
names of compared consecutive segments (facies/sub-facies/layers), "p_FOR_LOG_K", p values
for log10 K using two-sample K-S test with 95% CI.
SEGMENT
Facies_A-B
Facies_A-C
Facies_A-D
Facies_B-C
Facies_B-D
Facies_C-D
Facies_A1-A2
Facies_A2-B1
Facies_B1-B2
Facies_B2-C
Facies_C-D
Layer_01-02
Layer_02-03
Layer_03-04
Layer_04-05
Layer_05-06
Layer_06-07
Layer_07-08
Layer_08-09
Layer_09-10
Layer_10-11
Layer_11-12
Layer_12-13
Layer_13-14

p_FOR_LOG_K
8.06E-110
3.11E-158
5.09E-94
9.50E-44
8.80E-49
2.91E-149
8.67E-08
1.27E-21
4.20E-17
1.54E-11
2.91E-149
5.32E-16
3.26E-12
5.40E-05
5.02E-02
7.84E-05
2.26E-01
1.64E-01
1.97E-06
6.80E-05
2.39E-03
1.13E-02
1.23E-03
6.34E-07

37

Table A2.2 (cont'd)
SEGMENT
Layer_14-15
Layer_15-16
Layer_16-17
Layer_17-18
Layer_18-19
Layer_19-20
Layer_20-21
Layer_21-22
Layer_22-23
Layer_23-24
Layer_24-25
Layer_25-26
Layer_26-27
Layer_27-28
Layer_28-29
Layer_29-30

p_FOR_LOG_K
1.54E-08
7.16E-05
2.90E-10
5.73E-07
1.88E-07
1.09E-01
4.84E-06
2.73E-06
6.31E-03
1.24E-02
3.10E-02
1.44E-07
2.22E-10
2.05E-04
7.31E-12
2.83E-03

38

Figure A2.4 (a) Cross plot of log10 K mean and variance for high resolution K data from four
DP profiles (see Figure 2.1 for location and Figure 2.2 for data). (b) Cross plot of log10 K mean
and the slopes of linear trends of the change in log10 K with depth. In these plots, each facies is
shown by a different symbol and color. The number next to each symbol indicates the DP profile
(1, 2, 3, and 4 represent 111108A, 111108B, 111108C, and 121108A, respectively).

39

REFERENCES

40

REFERENCES
Benson, D. A., R. Schumer, M. M. Meerschaert, and S. W. Wheatcraft (2001), Fractional
dispersion, Lévy motion and the MADE tracer tests, Transp. Porous Media, 42(1-2), 211240.
Boggs, J. M., and E. E. Adams (1992), Field study of dispersion in a heterogeneous aquifer: 4.
Investigation of adsorption and sampling bias, Water Resources Research, 28(12), 33253336.
Bowling, J. C., A. B. Rodriguez, D. L. Harry, and C. Zheng (2005), Delineating alluvial aquifer
heterogeneity using resistivity and GPR data, Ground Water, 43(6), 890-903.
Eggleston, J., and S. Rojstaczer (1998), Identification of large-scale hydraulic conductivity
trends and the influence of trends on contaminant transport, Water Resources Research,
34(9), 2155-2168.
Engdahl, N. B., G. S. Weissmann, and N. D. Bonal (2010), An integrated approach to shallow
aquifer characterization: combining geophysics and geostatistics, Comput. Geosci., 14, 217229.
Ezzy, T. R., M. E. Cox, A. J. O'Rourke, and G. J. Huftile (2006), Groundwater flow modelling
within a coastal alluvial plain setting using a high-resolution hydrofacies approach; Bells
Creek plain, Australia, Hydrogeol. J., 14(5), 675-688.
Fleckenstein, J. H., and G. E. Fogg (2008), Efficient upscaling of hydraulic conductivity in
heterogeneous alluvial aquifers, Hydrogeol. J., 16(7), 1239-1250.
Gelhar, L. W. (1993), Stochastic Subsurface Hydrology, 390 pp., Prentice-Hall, New Jersey.
Grasmueck M., R. Weger, and H. Horstmeyer (2005), Full-resolution 3D GPR imaging.
Geophysics, 70(1), K12-K19.
Harvey, C. F., and S. M. Gorelick (2000), Rate-limited mass transfer or macrodispersion: Which
dominates plume evolution at the Macrodispersion Experiment (MADE) site?, Water
Resources Research, 36(3), 637-650.
Hess, K. M., S. H. Wolf, and M. A. Celia (1992), Large-scale natural gradient tracer test in sand
and gravel, Cape Cod, Massachusetts 3. Hydraulic conductivity variability and calculated
macrodispersivities, Water Resources Research, 28(8), 2011-2027.
Hubbard, S. S., and Rubin, Y. (2000), Hydrogeological parameter estimation using geophysical
data: A review of selected techniques, J. Contam. Hydrology, 45, 3-34.
Jol, H. M. (2009), Ground penetrating radar: Theory and Applications, 524 pp., Elsevier,
Amsterdam.

41

Julian, H. E., J. M. Boggs, C. Zheng, C. E. Feehley (2001), Numerical simulation of a natural
gradient tracer experiment for the natural attenuation study: flow and physical transport,
Ground Water, 39(4), 534-545.
Liu, G., J. J. Butler Jr., G. C. Bohling, E. Reboulet, S. Knobbe, and D. W. Hyndman (2009), A
new method for high-resolution characterization of hydraulic conductivity, Water Resources
Research, 45, W08202.
Liu, G., C. Zheng, G. R. Tick, J. J. Butler, and S. M. Gorelick (2010), Relative importance of
dispersion and rate-limited mass transfer in highly heterogeneous porous media: Analysis of a
new tracer test at the Macrodispersion Experiment (MADE) site, Water Resources Research,
46, W03524.
Mackay, D. M., D. L. Freyberg, P. V. Roberts, and J. A. Cherry (1986), A natural gradient
experiment on solute transport in a sand aquifer 1. Approach and overview of plume
movement, Water Resources Research, 22(13), 2017-2029.
Moysey S., J. Caers, R. Knight, and R. M. Allen-King (2003), Stochastic estimation of facies
using ground penetrating radar data, Stochastic Environmental Research and Risk
Assessment, 17(5), 306-318.
Neal, A. (2004), Ground-penetrating radar and its use in sedimentology: principles, problems and
progress, Earth-Sci. Rev., 66(3-4), 261-330.
Phanikumar, M. S., D. W. Hyndman, X. Zhao, and Dybas M. J. (2005), A three-dimensional
model of microbial transport and biodegradation at the Schoolcraft, Michigan, site, Water
Resources Research, 41, W05011.
Rauber, M., F. Stauffer, P. Huggenberger, and T. Dracos (1998), A numerical three-dimensional
conditioned/unconditioned stochastic facies type model applied to a remediation well system,
Water Resources Research, 34, 2225-2233.
Rehfeldt, K. R., J. M. Boggs, and L. W. Gelhar (1992), Field study of dispersion in a
heterogeneous aquifer 3. Geostatistical analysis of hydraulic conductivity, Water Resources
Research, 28(12), 3309-3324.
Van Overmeeren, R. A. (1998), Radar facies of unconsolidated sediments in The Netherlands: A
radar stratigraphy interpretation method for hydrogeology. J. App. Geophysics, 40, 1-18.
Whittaker, J., and G. Teutsch (1999), Numerical simulation of subsurface characterization
methods: application to a natural aquifer analogue, Adv. Water Resour., 22(8), 819-829.
Zheng, C. (2006), Accounting for aquifer heterogeneity in solute transport modeling: a case
study from the macrodispersion experiment (MADE) site in Columbus, Mississippi, in
Handbook of Groundwater Engineering, edited by J. W. Delleur, 18 pp., CRC Press.
Zheng, C., and S. M. Gorelick (2003), Analysis of solute transport in flow fields influenced by
preferential flow paths at the decimeter scale, Ground Water, 41(2), 142-155.

42

Chapter 3
Hydraulic conductivity fields: Gaussian or not?

Abstract
Hydraulic conductivity (K) fields are used to parameterize groundwater flow and transport
models. Numerical simulations require a detailed representation of the K field, synthesized to
interpolate between available data. Several recent studies introduced high resolution K data
(HRK) at the MAcro Dispersion Experiment (MADE) site, and used ground-penetrating radar
(GPR) to delineate the main structural features of the aquifer. This paper describes a statistical
analysis of these data, and the implications for K field modeling alluvial aquifers. Two striking
observations have emerged from this analysis. The first is that a simple fractional difference
filter can have a profound effect on data histograms, organizing non-Gaussian ln K data into a
coherent distribution. The second is that using GPR facies allows us to reproduce the
significantly non-Gaussian shape seen in real HRK data profiles, using a simulated Gaussian ln

K field in each facies. This illuminates a current controversy in the literature, between those who
favor Gaussian ln K models, and those who observe non-Gaussian ln K fields. Essentially, both
camps are correct, but at different scales.2

This chapter is based on Meerschaert, M.M., M. Dogan, R. L. Van Dam, D. W. Hyndman, and
D. A. Benson (2013), Hydraulic conductivity fields: Gaussian or not? Water Resources
Research, doi:10.1002/wrcr.20376.

43

3.1 Introduction
Groundwater flow and transport simulations require a densely defined hydraulic
conductivity (K) field to populate the model grid. Because it is not practical to collect 2-D or 3-D
data at this resolution, stochastic simulation methods are commonly used to interpolate between
measured data values. Stochastic K field simulation requires a statistical analysis of the available

K data, to ensure that the synthesized K field resembles the data in terms of its distribution and
correlation structure. The two main simulation steps are: (1) generate an uncorrelated noise field;
and (2) apply an appropriate filter to impose a correlation structure. Since random number
generators produce only uncorrelated noise, both steps are necessary. To parameterize the
simulation model, the process is reversed: (1) apply an appropriate inverse filter to the raw data
to remove the correlation; and (2) examine the filtered, uncorrelated data to determine its true
underlying distribution. Unless the data is filtered properly to remove correlations, the data
histogram can significantly misrepresent the underlying distribution, since a histogram of
correlated data need not reflect the true underlying distribution. In this chapter, we will see a
remarkable example of this simple and well-known fact.
Hydraulic conductivity data from the Macro Dispersion Experiment (MADE) site, at the
Columbus Air Force Base in Mississippi, clearly show a high level of heterogeneity in hydraulic
properties [Rehfeldt et al., 1992; Zinn and Harvey, 2003; Llopis-Albert and Capilla, 2009]. The
site was recently revisited to obtain K measurements with much higher spatial resolution than
previous measurements [Bohling et al., 2012; Liu et al., 2009]. Vertical columns (profiles) of
hydraulic conductivity data were measured at approximately 1.5 cm depth increments, using a
new direct-push profiling method that couples the direct-push injection logger (DPIL) and the
direct-push permeameter (DPP) [Butler Jr et al., 2007; Liu et al., 2009, 2012]. This novel high-

44

resolution K (HRK) tool was advanced into the subsurface, while water was injected out of a
small screened port located a short distance behind the tool tip. The injection rate, and injectioninduced back pressure, were recorded every 1.5 cm, and the ratio of these quantities was then
transformed into K estimates [Liu et al., 2009]. The cm-scale spatial resolution of the resulting K
data is orders of magnitude finer than the data considered in previous studies [Rehfeldt et al.,
1992; Meerschaert et al., 2004; Bohling et al. 2012] analyzed the resulting K data, and compared
those measurements to previous flowmeter-based K estimates collected at lower resolution
across the same site.
A parallel data collection effort used ground- penetrating radar (GPR) to image the
related sedimentary structures in the aquifer, called facies, by identifying distinct reflection
characteristics, such as reflection terminations, dip angles, amplitudes, and continuity. Such GPR
facies have been shown to correlate with hydrogeological units [e.g., see Van Overmeeren, 1998;

Heinz and Aigner, 2003; Schmelzbach et al., 2011]. Full-resolution 3D GPR data using 50 and
100 MHz antennae were obtained with step sizes (and line-spacings) of 0.2 and 0.1 m,
respectively, using a Sensors and Software pulseEKKO 100 system. Data processing and
analysis to extract facies boundaries was detailed in Dogan et al. [2011]. The map in Figure 3.1
outlines the GPR data collection site, and the location of the four HRK profiles that form the
basis for our study. The ICA (Intensively Cored Area) cube was the site of a push-pull tracer test
described in Liu et al. [2010], see also Zheng et al. [2011]. The MLS (Multi-Level Sampler) cube
was the site of the MADE-5 tracer test reported in Bianchi et al. [2011].
The modeling of hydraulic conductivity fields at the MADE site has been the focus of
intensive study and modeling for over twenty years. The geostatistical analysis of Rehfeldt et al.
[1992] documented a high level of heterogeneity, indicated by the variance of 4.5 for ln K in

45

their multi-Gaussian model, as well as anisotropy, indicated by horizontal and vertical
correlation scales of 12.8 m, and 1.6 m, respectively. Silliman and Wright [1988] and Rubin and

Journel [1991] argued that a Gaussian model with a single covariance function cannot reproduce
the preferential pathways (connected regions with the highest ln K values) observed in real
aquifers. Gómez-Hernádez and Wen [1998] continued this argument against the multi-Gaussian
model, and cautioned against drawing broad conclusions on the basis of one dimensional data
distributions. Renard and Allard [2011] survey several methods for characterizing connectivity,
and note that the multi-Gaussian model alone is often insufficient to reproduce the connectivity
observed in real aquifers. Significant deviations from a Gaussian profile were noted by Painter
[1996] and Meerschaert et al. [2004], and some alternative non-Gaussian models were proposed.

Zinn and Harvey [2003] point out that even in a model with Gaussian ln K profiles, deviation
from the usual multi-Gaussian model can lead to connected features. Salamon et al. [2007]
discuss the non-monotone variograms in MADE ln K data, and recommend a sequential
Gaussian simulation methodology with a non-monotone covariance structure, to reproduce this
“hole effect.” Llopis-Albert and Capilla [2009] use a gradual conditioning algorithm to produce
non-Gaussian ln K fields based on flowmeter, head, and concentration data from MADE-2. This
controversy between Gaussian and non-Gaussian ln K fields has profound implications for flow
and transport modeling. Heavy tailed ln K distributions support novel approaches including the
CTRW [Berkowitz et al., 2006], fractional ADE [Benson et al., 2013], and some related
stochastic hydrology models [Cushman and Ginn, 2000; Neuman and Tartakovsky, 2009], while
Gaussian ln K models are more consistent with the traditional ADE, mobile-immobile, and dualdomain models.

46

300

LEGEND
MADE-2 test boundary
MADE-3 source trench
2D GPR lines
3D GPR cubes
HRK profiles

N

ICA cube

Northing [m]

200

100
MLS cube

121108A
111108C
111108B
111108A

0
0

100

200
Easting [m]

300

Figure 3.1 Layout of Macro Dispersion Experiment test site, showing key features of MADE
experiments, as well as the locations of GPR data collected for this project. The inset of the
12×12 m ICA (Intensively Cored Area) cube shows the locations of the four HRK profiles and
the 2D transect discussed in this paper.
The two main findings of this study are that: (1) a fractional difference filter can be useful
to reveal the true underlying distribution of highly correlated vertical columns of HRK data; and
(2) using GPR facies, a multi-Gaussian simulation method with an appropriate operator scaling
correlation structure applied to each facies can reproduce the significantly non-Gaussian profiles
seen in columns of filtered HRK data. There remains a significant debate in the literature
between those who favor Gaussian models, and others who believe that a non-Gaussian approach
is needed. In our view, both groups are correct, albeit at different scales. Within a single facies,
an appropriate multi-Gaussian model can be effective, and when different facies are combined, a
non-Gaussian profile with a sharper peak and a heavier tail will emerge.

47

3.2 Statistical analysis
Many studies have analyzed the statistical properties of low-resolution K data profiles;
see Meerschaert et al. [2004] for a brief review. A typical field experiment collects K data at a
vertical resolution of 1 − 3 m. Since the vertical resolution of the new HRK data is orders of
magnitude finer, it is useful to reconsider the results of past analysis. For relatively homogeneous
aquifers, it has been common to employ a log-normal distribution for K: the distribution of ln K
is assumed to be normal, and aquifer heterogeneity is inferred from the variance of ln K [Rehfeldt

et al., 1992]. A more detailed analysis suggests a departure from normality, with a sharper peak
and heavier tails [e.g., Lu et al., 2002; Meerschaert et al., 2004]. This deviation becomes more
significant for aquifers that display a higher degree of heterogeneity.
Typical values for ln K are highly correlated, leading many researchers to employ models
such as a fractional Brownian motion. The MADE site is highly heterogeneous, with ln K
variance greater than 4.5. Several novel models have been proposed to try and capture this
combination of non-Gaussian distributions and strong correlations [e.g., Painter, 1996; Herrick

et al., 2002; Molz and Boman, 1993; Kohlbecker et al., 2006].
Figure 3.2a shows a histogram of ln K data from HRK profile 121108A (see map in
Figure 3.1). The histogram suggests a complex underlying distribution, widely varying with
several peaks, and no simple discernible shape. A fractional difference filter was applied to
remove the correlations, resulting in the histogram in Figure 3.2b. The filtered data show no
significant serial correlations, indicating that fractionally integrated noise is a reasonable model
for this vertical column of ln K data. Fractional models have been applied in hydrology since
pioneering work of Hurst [1951] on flood levels of the Nile river. In these models, observations

48

Xn are related to a sequence of independent and identically distributed random variables Zn
(white noise) by the fractional difference relation
∞

Zn = ෍ wj Xn-j
j=0

(3.1)

where the fractional binomial coefficients w୨ can be computed recursively using wn =1 and

wj =w ൫j-1-d൯/j for j≥1. If the underlying noise sequence Zn is Gaussian, then the sequence Xn
j-1
is a fractional Brownian motion with Hurst index H=d-1/2. See the Appendix 3.1 for more
details.

Figure 3.2 Histogram of ln K for HRK profile 121108A (see map in Figure 3.1) before (a) and
after (b) applying the fractional difference filter (1) with d=0.9. The filtered data are organized
into a unimodal distribution with a sharper peak and a heavier tail than the best fitting Gaussian
probability density function (black line).

49

Figure 3.2b was obtained using a fractional difference filter with d=0.9. The parameter d
was gradually increased until the autocovariance plot showed no significant correlations, see
Figure 3.3. The same value of d was effective in removing correlation in all four vertical HRK
columns 111108A, 111108B, 111108C, and 121108A (see Figure 3.1) that formed the basis for
our study. A fractional difference filter of order d=0.89 was used in Lu et al. [2002] to remove
correlations in laboratory ln K data from a vertical sandstone core; a value of d=0.9 was found
suitable for a sandstone slab in Major et al. [2011]; Meerschaert et al. [2004] used d=0.74 for
lower resolution ln K values from three horizontal profiles in a sandstone facies at a site in Utah.
The effect of fractional differencing on the histogram is striking. The filtered HRK data in
Figure 3.2b form a coherent shape, with a sharper peak and a heavier tail than a Gaussian (the
data fail the Anderson-Darling test for normality, p<0.0005. It is known that correlation can
distort a histogram, but we have never seen such a clear example in real data. The fractional
difference filter transforms a highly complex histogram into a form amenable to statistical
modeling, by removing the correlation. This is the first major finding of our statistical analysis:
A simple fractional difference filter is sufficient to capture and remove the correlation structure
of a vertical ln K profile. This filter reveals the underlying noise distribution needed to design a
faithful ln K field simulation.
The dramatic transformation between Figures 3.2a and 3.2b has not been observed
previously, perhaps because the available data were either more homogeneous (e.g., laboratory
studies of a sandstone slab) or more widely spaced (e.g., flowmeter data from field studies) than
the data considered in this study. Since our data are closely spaced, many similar K values tend
to clump together due to high correlations, creating histogram peaks. These high correlations are
evident in Figure 3.3a. Mathematically, this strong correlation is a fractional integration. Since

50

each vertical section spans several different facies, with significantly different material
properties, multiple peaks can occur in a single HRK profile. The fractional difference filters out
the correlations by reversing the fractional integration.

Figure 3.3 Autocorrelation function for ln K from HRK profile 121108A (see Figure 3.1) before
(a) and after (b) applying the fractional difference filter (3.1) with d=0.9. Autocorrelations inside
dashed lines are statistically zero.
Next, we discuss our simulation scheme. Since the ln K data exhibit long range
dependence, with a shorter correlation length in the vertical direction, we applied the anisotropic
random field generator of Benson et al. [2006]: Fourier transformed Gaussian white noise on a
1.5 cm grid was multiplied by a Fourier filter ψሺkሻ=

ሾ ∑i Ci |k·θi |2 ሿ-

H+1
2

with Hurst index

H=0.4, θ1 horizontal, θ2 vertical, and correlation length parameters C1 =10 and C2 =1 to produce
anisotropic ln K fields with a longer correlation length in the horizontal direction. In this

51

simulation, any horizontal row or vertical column of simulated data represents a fractional
Brownian motion with Hurst index H=0.4. The horizontal autocorrelation parameter C1 was
chosen to match measured autocorrelations between the four vertical HRK profiles 111108A,
111108B, 111108C, and 121108A used in this study, but due to the 4 m horizontal spacing, this
represents only a very rough fit. The simulated ln K field was then adjusted to impose the same
mean and standard deviation as the log-transformed HRK data. The conditioning algorithm of

Benson et al. [2013] was then applied, to make the simulated ln K field agree with available
HRK data. Figure 3.4a shows the results of this simulation procedure for the combined HRK
data, without subdividing into GPR facies.
Next the HRK data were segregated by facies, using the GPR method discussed in
Section 3.1. Both the mean and the standard deviation of ln K were found to vary significantly
between facies. Separate ln K fields with the same mean and standard deviation as the logtransformed HRK data were generated over the entire model domain for each facies, using the
same method of Benson et al. [2006, 2013] with the same filter, and the same white noise
sequence as in Figure 3.4a, with dip angle θ1 matched to the orientation of GPR reflections for
each facies. Then, GPR facies boundaries were used to cut out the relevant portions of the
simulated ln K field for each facies, resulting in the ln K field shown in Figure 3.4b. The
multiscaling fractal filter used in this simulation methodology produces enhanced connectivity,
as compared to a traditional multi-Gaussian model. Connectivity is further enhanced by our
facies approach, since ln K statistics vary by facies.

52

Probability Density

Probability Density

3

(c)

2
1
0
-1.5

-1

-0.5
0
0.5
1
Fractionally Differenced ln K

3

1.5

(d)

2
1
0
-1.5

-1

-0.5
0
0.5
1
Fractionally Differenced ln K

1.5

Figure 3.4 Simulated ln K field without (a) and with (b) GPR facies (dashed lines), conditioned
on four HRK profiles (vertical black lines). Histogram (c, d) of one column (white line at x = 172
m) from simulated ln K field (a, b, respectively) after applying the fractional difference filter (1)
with d=0.9. The histogram (c) fits a Gaussian model, but the histogram (d) from facies
simulation (b) deviates from Gaussian shape, similar to measured HRK data (Figure 3.2b).

53

Figure 3.5 shows the Gaussian fit to fractionally differenced ln K data in a single facies,
using the facies boundaries shown in Figure 3.4b. The data from facies A (shallowest) at
horizontal location 174 m had the smallest standard deviation (σ=0.0110). The probability plot in
Figure 3.5a shows that these data fit a Gaussian distribution reasonably well, except for a single
outlier (0.7244, removed). The histogram (not shown) is similar to Figure 3.4c. Figure 3.5b
shows the corresponding plot for facies D (deepest) at horizontal location 170 m, which had the
largest standard deviation ( σ=0.3637 ). Since the points on the probability plot in Figure 3.5b
show a significant and systematic deviation from the reference line, a lack of fit to the Gaussian
model is indicated. The histogram (not shown) is similar to Figure 3.4d. We attribute this
deviation from the Gaussian model in our data to the existence of sub-facies and smaller
sedimentary variations with significantly different material properties [Dogan et al., 2011]. In
this study, we employ only a few of the most definitive and connected GPR reflection
boundaries, to subdivide the model domain into four distinct facies. However, the full
geostatistical analysis reported in Dogan et al. [2011] did uncover additional substructures.

Zhang et al. [2013], and others referenced in Section 3.5 of that paper, find that sub-facies
heterogeneity has only a secondary influence on transport, hence the importance of accurately
modeling sub-facies is unclear.
As noted by Silliman and Wright [1988] and further discussed in Gómez-Hernádez and

Wen [1998], a multi-Gaussian simulation with a single covariance function will not produce
continuous regions where the highest or lowest ln K values occur. However, in our model, the
facies with the highest or lowest mean ln K value produce just such features. This is no
contradiction, because our model employs a different multi-Gaussian mean and covariance
structure in each facies. In our opinion, the “hole effect” in the variograms of Rehfeldt et al.

54

[1992] and Salamon et al. [2007] can be the result of combining data from distinct facies, which
will naturally cause a deviation from a single multi-Gaussian model with a fixed mean and
covariance structure. Furthermore, combining the simulated multi-Gaussian ln K values from
different facies does produce the kind of non-Gaussian histogram, with a sharper peak and a
heavier tail, frequently seen in column data.

Figure 3.5 Fractionally differenced ln K data (a) from the shallowest facies at horizontal location
174 m fits a Gaussian distribution. Deepest facies (b) at horizontal location 170 m deviates from
the Gaussian model. These probability plots show the sorted data on the horizontal axis, and the
corresponding model percentiles for the best fitting Gaussian model on the vertical axis. If the
data fits this model, the points will follow the reference line, with some random scatter.

55

This simulation methodology used to produce the ln K field in Figure 3.4b produces
results similar to the indicator geostatistics method of Fogg et al. [1998], which has been
successfully applied in both groundwater hydrology [Weissmann et al., 1999] and surface water
hydrology [Rubin et al., 2006]. The idea of combining fractal simulation methods with a facies
model is already present in Lu et al. [2002]. The difference in our approach is that we use GPR to
determine the facies boundaries. Since the actual facies boundaries are known, there is no need
to resort to an indicator simulation method to synthesize the facies boundaries. Ritzi [2000] notes
that lithofacies data can also be used to determine facies boundaries, but if aquifer lithology is
not available at sufficient resolution to parameterize a flow and transport model, then a
combination of GPR facies and HRK profiles can provide a useful modeling approach for highly
detailed K field synthesis.
3.3 Model validation
If a simulated ln K profile exhibits the statistical features of a measured ln K data profile,
then this validates the simulation methodology. The histogram shown in Figure 3.4c represents a
single column (81st column, at 172 m) of values from the simulated ln K field without GPR
facies in Figure 3.4a, fractionally differenced with d=0.9 as in Figure 3.2b. Without facies, the
fractionally differenced simulated HRK profile fits a Gaussian distribution, and hence does not
resemble the measured HRK data. Figure 3.4d shows the corresponding histogram from a single
column of the simulated ln K values with GPR facies in Figure 3.4b, fractionally differenced
with d=0.9. Before fractional differencing, the simulated profile histogram (not shown) appeared
similar to Figure 3.2a. After fractional differencing, the histogram of simulated values in Figure
3.4d appears quite similar to the corresponding histogram in Figure 3.2b, with a sharper peak and
a heavier tail than the best-fitting Gaussian. There are also some significant differences between

56

Figure 3.4d and Figure 3.2b, including a higher peak and some asymmetry in Figure 3.2b, but the
overall shape seems to support our conclusion that combining simulated Gaussian ln K values
from different facies can reproduce a significantly non-Gaussian shape, similar to what is seen in
real HRK data profiles. Even though the simulated noise is normal in each facies, the histogram
in Figure 3.4d does not fit a normal probability density (Anderson- Darling test p<0.0005). This
is due to the well-known fact that a mixture of Gaussian random variates with different mean
and/or standard deviation cannot be normally distributed. Indeed, many non-Gaussian
distributions that have been used to model ln K data, including the Laplace and symmetric stable,
are Gaussian mixtures [Kotz et al., 2001; Guadagnini et al., 2012; Riva et al., 2013a]. We
conclude that GPR facies are useful in this simulation, as they provide a data-based procedure for
delineating statistically distinct regions of K values, leading to the more sharply peaked and nonGaussian profile evident in Figure 3.4d. The facies approach also allows us to preserve observed
correlation structures and angles.
In order to gain a practical appreciation for the accuracy of d estimates, we then
simulated a number of statistically identical ln K fields, and applied automatic d estimation to the
resulting ln K profiles. Using a standard maximum likelihood estimation routine for fractional
ARIMA models, we found typical estimates of the d parameter to vary from the true (input)
value of d=0.9 by േ0.2 in those simulations. Hence we cannot rule out other values of d
(including d=1.0, a simple difference), and the estimated d value from any single profile should
only be taken as a rough indicator of the true value. However, since the value d=0.9 resulted in
no significant serial correlation in any of the four HRK profiles in this study, this value was
deemed adequate for our purposes. It is certainly possible that more significant variations in d
could emerge on a larger scale, or at a different site.

57

We believe that the departure from a Gaussian distribution, commonly observed in many
ln K data profiles from alluvial aquifers, can be attributed to mixing. Although our simulated ln

K field is based on Gaussian noise, the distribution of any single column exhibits a significant
non-Gaussian shape, because different facies are mixed. This leads to the second major finding
of our statistical analysis: A simulation that uses GPR facies, with a fractional Brownian motion
within each facies, generates ln K fields whose fractionally differenced vertical profiles have a
strongly non-Gaussian distribution, with a sharper peak and a broader tail, consistent with nonGaussian ln K models applied in previous studies [Meerschaert et al., 2004; Painter, 1996]. The
GPR data are valuable in this simulation method, since they delineate facies boundaries that
allow the Gaussian simulation to reproduce non-Gaussian ln K profiles.
3.4 Discussion
Modeling and simulation of K fields is challenging, especially in highly heterogeneous
aquifers including the MADE site, where the ln K fields exhibit anisotropy [Boggs et al., 1990;

Riva et al., 2008], long-range correlations [Neuman, 2003; Ritzi et al., 2000], non-monotone
variograms [Ritzi et al., 2004; Salamon et al., 2007], and a significantly non-Gaussian shape
[Ritzi et al., 2004]. The standard model of ln K is based on a normal distribution, but many
studies have found significant deviations from the Gaussian shape in increments of low
resolution ln K field data, with a sharper peak and/or heavier tails [Meerschaert et al., 2004;

Painter, 1996, 2001]. Some researchers have suggested that accurate representation of the K data
at the smallest scale may be a critical component of solute transport simulation, particularly
regarding the distribution and long-range dependence [Zheng et al., 2003; Ritzi et al., 2004; Dai

et al., 2004; Ramanathan et al., 2008]. Based on the statistical analysis reported in this paper, we
find that the observed non-Gaussian distribution of fractionally differenced ln K data from

58

alluvial aquifers can be reproduced using a Gaussian ln K field in each facies. The combination
of ln K data from different facies at different depths combines into a data profile with a nonGaussian shape (e.g., Painter, 2001). The non-Gaussian ln K distributions used in previous
studies are also Gaussian mixtures of this type. It is well known that mixing of data from
different populations changes the histogram shape, but it is usually impossible to reconstruct the
Gaussian components. This has led to popular indicator geostatistics methods that synthesize
facies boundaries (e.g., Weissmann et al., 1999). Using GPR facies, it does seem possible to
delineate the actual facies boundaries without resorting to simulations, and thereby reduce
measured ln K data to a reasonably Gaussian form. This allows a simple method for interpolating
highly variable and non-stationary ln K fields. Another advantage of facies modeling is laid out
in [Winter et al., 2002, 2003]: It facilitates efficient perturbation-based stochastic methods based
on locally homogeneous ln K fields. Riva et al. [2013b] has reported a significantly nonGaussian distribution of log permeability for the two faces parallel to bedding in a relatively
homogeneous sandstone slab, while the distribution on the other four faces was close to
Gaussian. Hence the Gaussian facies model promoted here may not be universally applicable.
3.5 Conclusions
In this paper, ground-penetrating radar (GPR) reflections were used to delineate facies
boundaries, and a high resolution fractal ln K field was simulated within each facies to
interpolate between available K data. There were two main findings of this study. First, a
fractional difference filter can be useful to capture the correlation structure of ln K profiles. The
unfiltered data histogram from one profile is severely distorted, but the filter uncovers a coherent
noise distribution, required for simulation design. Second, GPR data can be used to delineate
facies boundaries for the K field model. While the overall distribution of ln K profiles in a typical

59

alluvial aquifer deviates significantly from Gaussian, it is reasonable to model the ln K field
within each GPR facies as Gaussian. The deviation from Gaussian in the combined profile is the
result of mixing, since the combination of data from different Gaussian distributions will no
longer fit a Gaussian model. In past research, many investigators have assumed a Gaussian
model for ln K, while many others have presented strong evidence for non-Gaussian alternatives.
Our analysis indicates that both groups are correct, albeit at different scales, consistent with the
findings of Lu et al. [2002]. A Gaussian model with an appropriate correlation structure can be
adequate for a single facies. For a highly heterogeneous aquifer, comprised of significantly
different facies, the combination of ln K values with a different mean and variance in each facies
will produce significantly non-Gaussian profiles.
3.6 Acknowledgments
Support for this research was provided by NSF grants DMS-1025486, DMS-0803360,
EAR-0738938, EAR-0738955, NIH grant R01-EB012079, and PRF grant 48515-G8. Any
opinions, findings, and conclusions or recommendations expressed are those of the authors and
do not necessarily reflect the views of the funding agencies. We thank Geoff Bohling and Jim
Butler Jr. at the Kansas Geological Survey for insightful discussions. Thanks also to Hans-Peter
Scheffler at the University of Siegen for providing the computer simulation code for operator
scaling random fields.

60

APPENDICES

61

APPENDIX 3.1: Fractional difference filter
The fractional difference filter was pioneered by Hurst [1951] to remove correlation in
river flood level data. It has now become a standard tool in one dimensional time series analysis
[Brockwell and Davis, 1991] and multidimensional spatial statistics [Beran, 1994]. Given a
correlated time series Xn (or a spatial series collected at equally spaced intervals along a one
dimensional line), the backward shift operator BXn ൌX
facilitates a simple notation for the
n-1
fractional difference
∞
d X =ሺI-Bሻd X = ෍ w X ,
∆ n
n
j n-j
j=0

(A3.1)

where the fractional binomial coefficients

d ሺ-1ሻj Γሺd+1ሻ
wj =ሺ-1ሻj ൬ ൰ =
j
j!Γሺd-j+1ሻ

(A3.2)

using the natural extension of the integer order binomial coefficients. Using the well-known
property Γሺx+1ሻ=xΓሺxሻ of the gamma function, one can also write

wj =

െ݀ሺ1 െ ݀ሻ ··· ሺj െ 1 െ ݀ሻ
j!

(A3.3)

from which the recursive formula wj ൌw ൫j‐1‐d൯/j follows. Hence in the special case where d
j-1
is a positive integer, the sum (A3.1) is finite, since wj =0 when j ൐ ݀ . Integer order derivatives
are defined as the limit of difference quotients using these operators. In the same way, fractional
derivatives can be defined as the limit of fractional difference quotients [Meerschaert and

Sikorskii, 2012].

62

In time series and spatial statistics, an integer order difference is also useful to remove
trends, since for example the first order difference of a linear trend is a constant, and the second
order difference of a quadratic trend is also a constant. Since the goal is to filter out the
correlation (and possibly a trend), an effective fractional difference filter will output an
uncorrelated white noise Zn =∆d Xn . This is tested in practice by computing the sample
autocorrelation defined by ρሺhሻ =γሺhሻ/γሺ0ሻ where the sample autocovariance is defined by
ො
ො
ො
N
1
ഥ
ഥ൯ሺZn -Zሻ with Z= ෍ Zn
ഥ
γሺhሻ= ෍൫Zn+h -Z
ො
N
n=1
n=1
N-h

(A3.3)

for a data set of length N. Standard statistical theory [Brockwell and Davis, 1991] shows that, for
large N, the sample autocorrelation of an uncorrelated white noise at any lag h is approximately
normally distributed with mean zero and variance 1/N. Since this random quantity lies between

±1.96/√N approximately 95% of the time, the autocorrelation plots in Figure 3.3 show dashed
lines at ±1.96/√N. Then the correlation in the data is judged to be statistically insignificant
(statistically zero) if 95% of the sample autocorrelations ρሺhሻ lie within these bounds, and the
ො
remaining sample correlations do not lie very far outside these bounds. In this case, there is no
compelling evidence to contradict the (null) hypothesis that the correlation is zero at any lag

h ് 0. The data Zn in Figure 3.2b was obtained from equation (3.1) using the data Xn from
Figure 3.2a, and the optimal value of d=0.9 was determined by increasing d gradually until the
autocorrelation (Figure 3.3) was reduced to be statistically zero. It is also possible to obtain an
estimate of d for a single time series using standard maximum likelihood estimation routines for
fractional ARIMA models [Brockwell and Davis, 1991].

63

For spatial data in 2 or 3 dimensions collected at equally spaced grid points, a fractional
difference filter can be applied in each coordinate. The order d of the fractional difference filter
can vary with the coordinate to remove spatial correlations. The entire data set can be used to
estimate the order(s) of the fractional difference [Beran, 1994; Guo et al., 2009], which
facilitates a more accurate estimate of the d parameter(s).

64

APPENDIX 3.2: Flow-Transport Simulations3
The work presented in this chapter demonstrates the advantages and the details of
fractional differencing and fractal stochastic simulation methods for creating K fields. It also
emphasizes the contribution of GPR derived facies boundaries to define hydrologically different
parts of the aquifer. Although the K fields generated using a combination of GPR and HRK data
sets using fractal statistical methods exhibit a better statistical representation of K data, it is
important to also see the effects of these different K realizations on flow and transport
simulations. Therefore, this section involves flow transport simulations through the simulated K
fields to investigate the distinct effects of each data set.
Macro-scale tracer experiments at the MADE site resulted in the well-known nonGaussian behavior of tracer concentrations with multiple peaks and a heavy tail, as explained in
detail in Chapter 1.3. However, most of the earlier modeling efforts were not successful to
reproduce both components of the behavior simultaneously. This additional study also aims to
test whether it is possible to reproduce these components using new data sets (HRK and GPR)
and K fields.
Based on K fields discussed in this chapter, MODFLOW and MT3D software were used
to model hydraulic heads and conservative tracer transport. The model domain had 160x570
(width x height) cells representing a 16 by 9.8 m (width x height) vertical cross section of the
aquifer. The entire model was saturated during the simulation, with specified head boundaries on
the two sides of the domain, and no-flow boundaries above and below. The tracer was injected at
x=166 m with 100 mg/l concentration over the entire depth for two days to minimize potential
boundary related artifacts. Figure A3.6 (top) shows the simulated plume through the simple K
This appendix is not included in the manuscript.

65

model, assuming a constant K values in each facies, set to the median of the measured K values
in that facies, 500 days after the injection. Figure A3.6 (middle) shows the same plume simulated
through the K field (shown in Figure 3.4a) populated without using facies boundaries. Figure
A3.6 (bottom) shows the plume simulated through the K field populated using both HRK data
and facies boundaries. Flow and transport through the top and bottom facies is significantly
faster, with more dispersion than the middle two facies.

66

Figure A3.6 Simulated plume 500 days after the end of injection, using a constant K value in
each GPR facies (top), the K field from Figure 3.4a without GPR facies (middle), and the K field
from Figure 3.4b with GPR facies (bottom).

67

C/Cinj

0.03

↑ 0.041347
↑ 0.045132

0.02

↑ 0.065026

0.01
0
164
0.03

C/Cinj

(a)
simple K model
without facies
with facies

166

168

170

172
174
Distance [m]

↑ 0.041704
↑ 0.023139

176

178

180

(b)
simple K model
without facies
with facies

0.02
0.01
0

0

500

1000
Time [days]

1500

2000

Figure A3.7 Normalized concentration profiles (top) 500 days after the end of injection, and
normalized concentration breakthrough curves (bottom), measured at location x=172 m, from the
three simulations illustrated in Figure A3.6.
In order to achieve a detailed comparison of the three simulations illustrated in Figure
A3.6, concentration histories and profiles were created. Figure A3.7 (top) shows the normalized
concentration profiles for the simulated plumes. Figure A3.7 (bottom) shows the normalized
concentration histories at location x=172 m, black dashed lines in Figure A3.6. These curves
demonstrate the impact of each data set (GPR and HRK) on plume behavior. The simple K
model represents what can be accomplished using only GPR facies. The resulting curves show
smooth, symmetrical variations in concentration that do not compare favorably to the field
measurements at the MADE site. This demonstrates the value of detailed HRK conditioning data
for capturing natural plume roughness. The simulated curves without facies seem more realistic,
but significantly smoother than the model with facies. The breakthrough curve has multiple

68

peaks and a heavy tail, consistent with tracer tests at the MADE site. I conclude that HRK data is
important for characterizing micro-scale variations, and GPR facies are useful for delineating
meso-scale variation. Together, they allow a simple Gaussian ln K model in each facies to
reproduce a realistic K field, consistent with the measured HRK data, and resulting in plausible
simulated plume behavior.

69

REFERENCES

70

REFERENCES
Benson, D. A., M. M. Meerschaert, B. Baeumer, and H. P. Scheffler (2006), Aquifer operatorscaling and the effect on solute mixing and dispersion, Water Resources Research, 42(1),
doi:10.1029/2004WR003755.
Benson, D. A., M. M. Meerschaert, and J. Revielle (2013), Fractional calculus in hydrologic
modeling: A numerical perspective, Advances in Water Resources, 51, 479–497,
doi:10.1016/j.advwatres.2012.04.005.
Beran, J. (1994), Statistics for long-memory processes. Chapman & Hall/CRC Press, Boca
Raton, Florida.
Berkowitz, B., A. Cortis, M. Dentz, and H. Scher (2006), Modeling non-Fickian transport in
geological formations as a continuous time random walk, Rev. Geophys., 44, RG2003.
Bianchi, M., C. Zheng, G. R. Tick, and S. M. Gorelick (2011), Investigation of small-scale
preferential flow with a forced-gradient tracer test, Ground Water, 49(4), 503–514.
Boggs, J. M., S. C. Young, D. J. Benton, and Y. C. Chung (1990), Hydrogeologic
characterization of the MADE site, epri en-6915, research project 2485-5, interim report,
Electric Power Research Institute, Palo Alto, CA.
Bohling, G. C., G. Liu, S. J. Knobbe, E. C. Reboulet, D. W. Hyndman, P. Dietrich, and J. J.
Butler Jr (2012), Geostatistical analysis of centimeter-scale hydraulic conductivity variations
at the MADE site, Water Resources Research, 48(2), doi:10.1029/2011WR010791.
Brockwell, P. J. and R. A. Davis (1991), Time Series: Theory and Methods, 2nd ed. New York:
Springer-Verlag.
Butler Jr, J. J., P. Dietrich, V. Wittig, and T. Christy (2007), Characterizing hydraulic
conductivity with the direct-push permeameter, Ground Water, 45(4), 409–419.
Cushman, J. H., and T. R. Ginn (2000), Fractional advection-dispersion equation: A classical
mass balance with convolution-Fickian flux, Water Resources Research, 36(12), 3763– 3766.
Dai, Z., R. W. Ritzi Jr., C. Huang, Y. N. Rubin and D. F. Dominic (2004), Transport in
heterogeneous sediments with multimodal conductivity and hierarchical organization across
scales, J. Hydrology, 294(1–3), 68–86.
Dogan, M., R. L. Van Dam, G. C. Bohling, J. J. Butler Jr, and D. W. Hyndman (2011),
Hydrostratigraphic analysis of the MADE site with full-resolution GPR and direct-push
hydraulic profiling, Geophys. Res. Lett., 38(6), doi:10.1029/2010GL046439.
Fogg, G. E., C. D. Noyes, and S. F. Carle (1998), Geologically based model of heterogeneous
hydraulic conductivity in an alluvial setting, Hydrogeology J., 6(1), 131–143.

71

Gómez-Hernádez, J. J., and X.-H. Wen (1998), To be or not to be multi-Gaussian? A reflection
on stochastic hydrogeology, Adv. Water Resour. 21(1), 47–61.
Guadagnini, A., M. Riva, and S. P. Neuman (2012), Extended power-law scaling of heavy-tailed
random air-permeability fields in fractured and sedimentary rocks. Hydrol. Earth Syst. Sci.
16, 3249–3260, doi:10.5194/hess-16-3249-2012.
Guo, H., C. Y. Lim, and M. M. Meerschaert (2009), Local Whittle estimator for anisotropic
random fields. J. Multivariate Anal. 100(5), 993–1028.
Heinz, J., and T. Aigner (2003), Hierarchical dynamic stratigraphy in various Quaternary gravel
deposits, Rhine glacier area (SW Germany): implications for hydrostratigraphy, International
Journal of Earth Sciences, 92(6), 923–938.
Herrick, M. G., D. A. Benson, M. M. Meerschaert, and K. R. McCall (2002), Hydraulic
conductivity, velocity, and the order of the fractional dispersion derivative in a highly
heterogeneous system, Water Resources Research, 38(11), 1227.
Hurst, H. E. (1951), Long-term storage capacity of reservoirs, Trans. Amer. Soc. Civil Eng., 116,
770–808.
Kohlbecker, M. V., S. W. Wheatcraft, and M. M. Meerschaert (2006), Heavy-tailed log
hydraulic conductivity distributions imply heavy-tailed log velocity distributions, Water
Resources Research, 42(4), doi:10.1029/2004WR003815.
Kotz, S., T. J. Kozubowski, and K. Podgorski (2001), The Laplace Distribution and
Generalizations: A Revisit with Applications to Communications, Economics, Engineering,
and Finance, Birkhäuser, Boston.
Liu, G., J. J. Butler Jr, G. C. Bohling, E. Reboulet, S. Knobbe, and D. W. Hyndman (2009), A
new method for high-resolution characterization of hydraulic conductivity, Water Resources
Research, 45(8), doi:10.1029/2009WR008319.
Liu, G., C. Zheng, G. R. Tick, J. J. Butler Jr, and S. Gorelick (2010), Relative importance of
dispersion and rate-limited mass transfer in highly heterogeneous porous media: Analysis of
a new tracer test at the Macrodispersion Experiment (MADE) site, Water Resources
Research, 46(3), W03524.
Liu, G., J. J. Butler Jr, E. Reboulet, and S. Knobbe (2012), Hydraulic conductivity profiling with
direct-push methods, Grundwasser, 17(1), 19–29.
Llopis-Albert, C., and J. E. Capilla (2009), Gradual conditioning of non-Gaussian transmissivity
fields to flow and mass transport data. 3 Application to the Macrodispersion Experiment
(MADE-2) site, on Columbus Air Force Base in Mississippi (USA) J. Hydro. 371, 75–84.
Lu, S., F. J. Molz, G. E. Fogg, and J. W.Castle (2002), Combining stochastic facies and fractal
models for representing natural heterogeneity, Hydrogeology J., 10(4), 475–482.

72

Major, E., D. A. Benson, J. Revielle, H. Ibrahim, A. Dean, R. M. Maxwell, E. Poeter, and M.
Dogan (2011), Comparison of Fickian and temporally nonlocal transport theories over many
scales in an exhaustively sampled sandstone slab, Water Resources Research, 47(10),
W10519.
Meerschaert, M. M., T. J. Kozubowski, F. J. Molz, and S. Lu (2004), Fractional Laplace model
for hydraulic conductivity, Geophys. Res. Lett., 31, 1–4.
Meerschaert, M. M. and A. Sikorskii (2012), Stochastic Models for Fractional Calculus. De
Gruyter Studies in Mathematics 43, De Gruyter, Berlin.
Molz, F. J., and G. K. Boman (1993), A fractal-based stochastic interpolation scheme in
subsurface hydrology, Water Resources Research, 29(11), 3769–3774.
Neuman, S. P. (2003), Relationship between juxtaposed, overlapping, and fractal representations
of multimodal spatial variability, Water Resources Research, 39(8), 1205.
Neuman, S. P., and D. M. Tartakovsky (2009), Perspective on theories of non-Fickian transport
in heterogeneous media, Adv. Water Resour., 32, 670-680.
Painter, S. (1996), Evidence for non-Gaussian scaling behavior in heterogeneous sedimentary
formations, Water Resources Research, 32(5), 1183–1195.
Painter, S. (2001), Flexible scaling model for use in random field simulation of hydraulic
conductivity, Water Resources Research, 37(5), 1155–1163.
Ramanathan, R., R. W. Ritzi,and C. Huang (2008), Linking hierarchical stratal architecture to
plume spreading in a Lagrangian-based transport model, Water Resources Research, 44(4),
W04503.
Rehfeldt, K. R., J. M. Boggs, and L. W. Gelhar (1992), Field study of dispersion in a
heterogeneous aquifer: 3. Geostatistical analysis of hydraulic conductivity, Water Resources
Research, 28(12), 3309–3324.
Renard, P., and D. Allard (2011), Connectivity metrics for subsurface flow and transport, Adv.
Water Resour., 51, 168–196.
Ritzi, R. W., D. F. Dominic, A. J. Slesers, C. B. Greer, E. C. Reboulet, J. A. Telford, R. W.
Masters, C. A. Klohe, J. L. Bogle, and B. P. Means (2000), Comparing statistical models of
physical heterogeneity in buried-valley aquifers, Water Resources Research, 36(11), 3179–
3192.
Ritzi, R. W. (2000), Behavior of indicator variograms and transition probabilities in relation to
the variance in lengths of hydrofacies, Water Resources Research, 36(11), 3375–3381.
Ritzi, R. W., Z. Dai, D. F. Dominic, Y. N. Rubin (2004), Spatial correlation of permeability in
cross-stratified sediment with hierarchical architecture, Water Resources Research, 40(3),
W03513.

73

Riva, M., A. Guadagnini, D. Fernàndez-Garcia, X. Sanchez-Vila, and T. Ptak (2008), Relative
importance of geostatistical and transport models in describing heavily tailed breakthrough
curves at the Lauswiesen site, J. Contam. Hydrol., 101(1-4), 1–13.
Riva, M., S.P. Neuman, and A. Guadagnini (2013), Sub-Gaussian model of processes with heavy
tailed distributions applied to permeabilities of fractured tuff. Stoch. Environ. Res. Risk
Assess. 27, 195–207, doi:10.1007/s00477-012-0576-y.
Riva M., S.P. Neuman, A. Guadagnini (2013), Anisotropic scaling of Berea sandstone log air
permeabilities statistics, Vadose Zone Journal, in press, doi: 10.2136/vzj2012.0153.
Rubin, Y. and A. G. Journel (1991), Simulation of non-Gaussian space random functions for
modeling transport in groundwater. Water Resources Research, 27, 1711–1721.
Rubin, Y., I. A. and Lunt, and J. S. Bridge (2006), Spatial variability in river sediments and its
link with river channel geometry, Water Resources Research, 42, W06D16.
Salamon, P., D. Fernàndez-Garcia, and J. J. Gómez-Hernández (2007), Modeling tracer transport
at the MADE site: the importance of heterogeneity. Water Resources Research, 43, W08404.
Schmelzbach, C., J. Tronicke, P. Dietrich (2011), Three-dimensional hydrostratigraphic models
from ground-penetrating radar and direct-push data. Journal of Hydrology, 398, 235-245.
Silliman, S. E., and A. L. Wright (1998), Stochastic analysis of paths of high hydraulic
conductivity in porous media. Water Resources Research, 24, 1901–1910.
Tennekoon, L., M. C. Boufadel, D. Lavalee, and J. Weaver (2003), Multifractal anisotropic
scaling of the hydraulic conductivity, Water Resources Research, 39(7), 1193,
doi:10.1029/2002WR001645.
Van Overmeeren, R. A. (1998), Radar facies of unconsolidated sediments in the Netherlands: A
radar stratigraphy interpretation method for hydrogeology, Journal of Applied Geophysics,
40(1), 1–18.
Weissmann, G. S. and G. E. Fogg (1999), Multi-scale alluvial fan heterogeneity modeled with
transition probability geostatistics in a sequence stratigraphic framework, J. Hydrology,
226(1), 48–65.
Winter, C.L., D.M. Tartakovsky, and A. Guadagnini (2002), Numerical solutions of moment
equations for flow in heterogeneous composite aquifers. Water Resources Research, 38(5),
10.1029/2001WR000222.
Winter, C.L., D.M. Tartakovsky, and A. Guadagnini (2003), Moment differential equations for
flow in highly heterogeneous porous media, Surveys in Geophysics, 24(1), 81–106.
Zhang, Y., C. T. Green, and G. E. Fogg (2013), The impact of medium architecture of alluvial
settings on non-Fickian transport, Adv. Water Resour., 54, 78-99.

74

Zheng, C., and S. M. Gorelick (2003), Analysis of Solute Transport in Flow Fields Influenced by
Preferential Flowpaths at the Decimeter Scale, Ground Water, 41(2), 142–155.
Zheng, C., M. Bianchi, and S. M. Gorelick (2011), Lessons learned from 25 years of research at
the MADE site, Ground Water, 49(5), 649–662.
Zinn, B., and C. F. Harvey (2003), When good statistical models of aquifer heterogeneity go bad:
A comparison of flow, dispersion, and mass transfer in connected and multivariate Gaussian
hydraulic conductivity fields, Water Resources Research, 42(3), 1051.

75

Chapter 4
Novel characterization method provides a major advance for
flow and transport prediction

Abstract
Predicting the fate and transport of solutes in aquifers is a societal grand challenge [Cantor,
1997] that requires characterization of aquifer properties. Although many studies have
reasonably explained transport through mildly heterogeneous aquifers, transport through highly
heterogeneous aquifers has yet to be accurately predicted. The complex transport behavior in
highly heterogeneous aquifers has fueled an ongoing debate in the hydrology community for
more than two decades [e.g., Berkowitz et al., 2002; Boggs et al., 1992; Dagan et al., 1992;

Dagan and Neuman, 1997; Zimmerman et al., 1998]. A wide range of modeling approaches have
been used to simulate tracer and contaminant transport at such sites, ranging from the classical
advection-dispersion equation (ADE) to dual-domain mass transfer and methods that impart
preferential flow paths to describe the observed complex behavior. A striking outcome of this
research to date has been the inability of approaches to reasonably predict transport in highly
heterogeneous systems based solely on field data. Here we demonstrate that when supported by
direct push high-resolution characterization data, advection dispersion equation can accurately
predict flow and transport without the need for calibration or addition of features that have not
been observed in the field. The ramifications for practical issues, such as the design of effective
remediation schemes and reliable risk assessments, are profound.

76

4.1 Introduction
Hydraulic conductivity (K) is the main property that controls solute transport in
subsurface flow systems [e.g., Dagan, 1989]. Thus, accurate predictions of flow and transport
through porous media require good representations of the K distribution [e.g., Molz et al., 1986].
Examples of successful applications of ADE models based on statistical distributions of K data
from mildly heterogeneous aquifers include those from Borden Ontario [e.g., Freyberg, 1986;

Mackay et al., 1986], Cape Cod, MA [e.g., Garabedian et al., 1991; LeBlanc et al., 1991], and
Schoolcraft, MI [Hyndman et al., 2000, Phanikumar et al., 2005].
However, flow and transport through highly heterogeneous sites has proven to be much
more challenging [e.g., Eggleston and Rojstaczer, 1998; Whittaker and Teutsch, 1999]. The
MAcro-Dispersion Experiment (MADE) site in Columbus, Mississippi is one of the best-studied
highly heterogeneous sites (variance of ln K = 4.5) [Boggs et al., 1992]. Much of the debate
about the applicability of the ADE at highly heterogeneous sites started with the experiments at
the MADE site. The outcome of the natural gradient tracer tests at this site highlighted the
limitations of existing and new modeling approaches.
Two large-scale tracer tests were performed at the site; in this paper we focus on the most
studied experiment called MADE-1, conducted from October 1986 to June 1988 [Boggs, 1991].
This test involved injection of bromide solution through an array of five boreholes,
approximately perpendicular to natural gradient of groundwater flow. The injection was followed
by collection of thousands of samples from 258 monitoring wells that were later analyzed for
tracer concentrations. The resulting data showed non-Gaussian tracer migration that was very
different than predicted using the ADE model based on K data from borehole flowmeters. Low
concentrations were detected far down gradient starting soon after the injection, yet more than

77

20% of the mass stayed within 6-7 meters of the injection wells for the duration of the
experiments [Adams and Gelhar, 1992].
Several approaches have been proposed to explain the observed tracer transport at the site
including: incorporating preferential flow paths [Zheng and Gorelick, 2003], multi-indicator
models [Fiori et al., 2013], mass transfer between mobile and immobile domains [Harvey and

Gorelick, 2000], and using a fractional form of the ADE [Benson et al., 2001]. Although some of
these approaches reasonably represented the average plume behavior in one dimension, they did
not accurately replicate several important characteristics of the plume including the spatial
extent. In addition, these methods would be difficult to use in a predictive sense as they are not
easily parameterized based solely on field data. We postulate that the difficulty in predicting
solute transport at heterogeneous sites like MADE is at least partly due to the absence of
sufficient high-resolution hydraulic property data to parameterize the models.
4.2 Methods
Here we present an approach to reproduce spatial extent of the tracer plume measured
during the MADE experiments. To this effect, we collected a novel suite of high-resolution K
data using a recently developed in-situ measurement method called the High-Resolution
Hydraulic Conductivity (HRK) tool [Liu et al., 2009]. This direct push tool can collect a 10 m
long profile in two to four hours, providing hydraulic conductivity estimates with an
unprecedented 1.5 cm vertical resolution. This allows important characterization of much smaller
scale transport features than previously possible [Dogan et al., 2011; Bohling et al., 2012].
In this study, we exploit the power of this novel tool along with fractal stochastic
methods to populate 3D hydraulic conductivity fields. Several studies have demonstrated that

78

fractional differencing can remove long range dependence, which allows the investigation of the
underlying distribution of K data [e.g., Benson et al., 2001]. Moreover, fractals are well-suited to
represent the connectivity of natural phenomena. Parameters for the fractal K field generation
were calculated based on autocorrelation and variogram analysis of the HRK data [Meerschaert et

al., 2013], which were also used as hard conditioning values for the stochastic fields.
The flow and transport model domain for this study includes a region of high density
HRK data near the tracer injection area of the MADE experiments (Figure 4.1). The model grid
was defined with over 3.2 million 0.25 x 0.25 x 0.05 m (length x width x height) cells oriented
with the long axis parallel to the average downgradient direction of the observed plumes. The
east and west edges of the model, which are approximately along flow paths, were assumed to be
no flow boundaries. The north and south edges of the model were assumed constant head
boundaries, and were assigned head values based on the average measured heads during the first
503 days of the MADE-1 tracer experiment. The maximum measured water level increase (0.64
m; Boggs et al. [1992]) across the five injection wells was used as a test of the reasonableness of
simulated K fields. Of the 20 generated stochastic K fields, the 6 that provided maximum
simulated head increases within 25% of the measured value (0.64 ±0.16 m) were selected for
transport simulations.
Transport simulations, based on the classical ADE, were performed in MT3D [Zheng and

Wang, 1999] without any calibration or parameter estimation. Parameters for the transport
simulations were defined based on literature as follows: longitudinal dispersivity of 0.05 m
[Gelhar, 1993], transverse horizontal and vertical dispersivity ratios of 0.1 and 0.01; effective
-6

molecular diffusion coefficient of 10 ; vertical and horizontal anisotropy in hydraulic
conductivity of 1 since we assume that the anisotropy will be introduced by our highly

79

heterogeneous and high-resolution K fields; and constant porosity of 0.35 [Adams and Gelhar,
1992]. Transport simulations were based on the 3rd order Total Variation Diminishing scheme,
which is mass conservative with minimal numerical dispersion. Total injected mass in the six
simulations averaged 25.1 kg with a standard deviation of 0.143 kg, which compares favorably
with the injected mass of 25.0 kg during the experiment.
300
(a)

Injection wells
MLS wells
Flowmeter K wells
HRK profiles

250
45

(b)

40

200

35
30

150

Distance [m]

Distance [m]

GPR cube

100

25
20
15

50

0

N

0

10

MADE
test boundary

50
100
Distance [m]

5
150

0

0

10
20
Distance [m]

Figure 4.1 (a) Map of the MADE site with the test boundary (dashed line) and sampling
locations; gray shaded rectangle shows the model domain used for simulations. (b) Model
domain with HRK, injection, and observation borehole locations.

80

4.3 Results
We first compare the simulated versus observed transport using 1D mass distribution
profiles, which has been the community standard for this site [e.g., Adams and Gelhar, 1992;

Fiori et al. 2013]. We then provide the first ever 2D comparison. Only concentration
measurements and simulated values that were larger than specified detection limits (0.01 mg/l for
the presented snapshot) were used in this analysis. To replicate the experimental procedure we
sampled the simulations at the location and depth of every multi-level sampler included in the
field experiment. The sampled concentrations were then interpolated in 3D using a 1 x 1 x 0.5 m
(dx, dy, dz) grid, and integrated vertically (along z-axis) and horizontally (along y-axis)

0.03

0.06
0.04

0.02

0.02

0.01

0

0

5

10

15

20
25
Distance [m]

30

35

40

Relative mass

Relative mass

following the procedure in Adams and Gelhar [1992].

0
45

Figure 4.2 Relative mass distribution profiles 501 days after the injection. The black line shows
the mean of the simulations along with 2σ error range (in gray shaded area); the red line shows
the profile for the MADE-1 observations.
The t=503 day snapshot was chosen for comparison as it is the most commonly presented
in publications with MADE-1 simulations (Figure 4.2) in comparison to the measured relative
mass distribution profile. The shape of the mean simulated curve is strikingly similar to the
experiment. The location, extent of the peaks and the tailing behavior of the mean simulated
curve reasonably represent the measured mass distribution, although the average of the

81

simulations shows somewhat heavier tails than was observed. Both could be improved by
optimization or calibration of porosity, but this would reduce the predictive power of the
approach.
Distance [m]
0
10
20
(a)

Distance [m]
0
10
20
(b)

Distance [m]
0
10
20
(c)

40

40

30

30

30

20

20

20

10

Distance [m]

40

10

10

0
Relative mass (x10-3 )
0

1

2

3

0
Relative mass (x10 -3)
5

10

15

0
Relative mass (x10 -3)
2

4

6

8

Figure 4.3 Vertically integrated contour maps of relative mass 503 days after the start of
injection for MADE-1 experiment (a), the mean of the 6 simulations that met the head criteria (b)
a sample simulation which matches the head change criteria best (c). All three contour maps
were created following the same interpolation procedure, as discussed in the text.
Despite the large quantity of research that has focused on the MADE site tracer test data,
literature contains no examples comparing simulation results with experimental data in two
dimensions. Figure 4.3a shows a map of vertically integrated relative mass for the experiment
data, which can be compared to the mean of 6 simulations that match the head change criteria
(Figure 4.3b) and to the simulation that best matches the observed head increase during the tracer
injection (∆h=0.637 m; Figure 4.3c). The measured extent of the plume, defined by the relative
mass parallel and transverse to flow as 22 and 17 m, respectively, is similar to the simulated

82

values of 28 and 15 m (Figure 4.3a, c). Also, the measured plume center of mass is very similar
to those of the simulated plumes.
4.4 Discussion and Conclusions
In this paper we demonstrate that, contrary to common belief, the classic ADE can
reproduce the measured distribution of MADE tracer concentrations. Our approach is different
from previous research to address this problem in that it uses a novel set of high-resolution K
data in combination with a fractal method to generate stochastic K fields. No calibration or
parameter estimation was used to improve the fit to measured data.
As with any transport modeling approach, there are differences between simulated and
observed concentrations. The total mass of the simulations is approximately twice the measured
mass. Other studies have shown similar observations, which may be due to a range of factors
including entrapment of tracer in the vadose zone [Adams and Gelhar, 1992]. We tested the
influence of transient flow behavior using almost two years of temporal water level
measurements, and found that this did not explain the difference. Other possible reasons include
differences between the simulated and field sampling procedure, which involves pumping
sample out from different levels of boreholes where small scale mixing of water with higher and
lower tracer concentrations will occur. Another possible explanation for the difference is due to
the absence of HRK profiles in close proximity to the injection wells. However, these data could
not be collected as the injection area was dug out and homogenized for the MADE 3 tracer test in
the 1990s. Inclusion of K data from flowmeter measurements from this area did not significantly
improve the recovered mass discrepancy.

83

Overall, the approach presented here successfully reproduced nearly all aspects of the
observed tracer plume without calibration. This includes the heavy tails that have previously only
been reproduced using transport theories such as dual domain mass transport or imposed high K
pathways. Predictive solute transport based solely on field data, such as done here, has immense
value for improving the ability to design more effective remediation schemes.

84

REFERENCES

85

REFERENCES
Adams, E. E., and L. W. Gelhar (1992), Field study of dispersion in a heterogeneous aquifer: 2.
Spatial moments analysis, Water Resources Research, 28(12), 3293–3307,
doi:10.1029/92WR01757.
Benson, D. A., R. Schumer, M. M. Meerschaert, and S. W. Wheatcraft (2001), Fractional
dispersion, Levy motion, and the MADE tracer tests. Transport in Porous Media 42: 211.
Berkowitz, B., J. Klafter, R. Metzler, H. Scher (2002), Physical pictures of transport in
heterogeneous media: Advection-dispersion, random-walk, and fractional derivative
formulations. Water Resources Research, 38: 1191.
Bohling, G. C., G. Liu, S. J. Knobbe, E. C. Reboulet, D. W. Hyndman, P. Dietrich, and J. J.
Butler Jr. (2012), Geostatistical analysis of centimeter-scale hydraulic conductivity variations
at the MADE site, Water Resources Research, 48: W02525.
Boggs, J. M. (1991), Database for the First Macrodispersion Experiment (MADE-1). Electric
Power Research Institute.
Boggs, J. M., S. C.Young, L. M. Beard, L. W. Gelhar, K. R. Rehfeldt, and E. E. Adams (1992),
Field study of dispersion in a heterogeneous aquifer: 1. Overview and site description. Water
Resources Research, 28, 3281-3291.
Cantor, K. P. (1997), Drinking water and cancer. Cancer Causes and Control 8: 292-308.
Dagan, G. (1989), Flow and Transport in Porous Formations. Springer-Verlag GmbH & Co.KG.
Dagan, G., and S. P. Neuman (1997), Subsurface flow and transport: A stochastic approach.
International Hydrology Series.
Dagan, G., V. Cvetkovic, and A. Shapiro (1992), A solute flux approach to transport in
heterogeneous formations. 1. The general framework. Water Resources Research, 28: 1376.
Dogan, M., R. L. Van Dam, G. C. Bohling, J. J. Butler Jr., D. W. Hyndman (2011),
Hydrostratigraphic analysis of the MADE site with full-resolution GPR and direct-push
hydraulic profiling. Geophysical Research Letters 38: L06405.
Eggleston, J., and S. Rojstaczer (1998), Identification of large-scale hydraulic conductivity
trends and the influence of trends on contaminant transport, Water Resources Research, 34:
2155–2168.
Fiori, A., G. Dagan, I. Jankovic, and A. Zarlenga (2013), The plume spreading in the MADE
transport experiment: Could it be predicted by stochastic models? Water Resources Research,
49.

86

Freyberg, D. L. (1986), A natural gradient experiment on solute transport in a sand aquifer: 2.
Spatial moments and the advection and dispersion of nonreactive tracers. Water Resources
Research, 22: 2031–2046.
Garabedian, S. P., D. R. LeBlanc, L. W. Gelhar, and M. A. Celia (1991), Large-scale natural
gradient tracer test in sand and gravel, Cape Cod, Massachusetts: 2. Analysis of spatial
moments for a nonreactive tracer. Water Resources Research, 27: 911–924.
Gelhar, L. W. (1993), Stochastic subsurface hydrology. Prentice-Hall Inc., Englewood, NJ.
Harvey, C., and S. M. Gorelick (2000), Rate-limited mass transfer or macrodispersion: Which
dominates plume evolution at the Macrodispersion Experiment (MADE) site? Water
Resources Research, 36: 637–650.
Hess, K. M., S. H. Wolf, and M. A. Celia (1992), Large-scale natural gradient tracer test in sand
and gravel, Cape Cod, Massachusetts: 3. Hydraulic conductivity variability and calculated
macrodispersivities. Water Resources Research, 28: 2011–2027.
Hyndman, D. W., M. J. Dybas, L. Forney, R. Heine, T. Mayottee, M. S. Phanikumar, G. Tatara,
J. Tiedje, T. Voice, R. Wallace, D. Wiggert, X. Zhao, and C. S. Criddle (2000), Hydraulic
characterization and design of a full-scale biocurtain. Ground Water 38: 462–474.
Liu, G., C. Zheng, G. R. Tick, and J. J. Butler Jr., and S. M. Gorelick (2010), Relative
importance of dispersion and rate-limited mass transfer in highly heterogeneous porous
media: Analysis of a new tracer test at the macrodispersion experiment (MADE) site. Water
Resources Research, 46: W03524.
Liu, G., J. J. Butler Jr., G. C. Bohling, E. Reboulet, S. Knobbe, and D. W. Hyndman (2009), A
new method for high-resolution characterization of hydraulic conductivity. Water Resources
Research, 45: W08202.
Koltermann, C. E., and S. M. Gorelick (1996), Heterogeneity in sedimentary deposits: A review
of structure-imitating, process-imitating, and descriptive approaches. Water Resources
Research, 32: 2617–2658.
LeBlanc, D. R., S. P. Garabedian, K. M. Hess, L. W. Gelhar, R. D. Quadri, K. G. Stollenwerk,
and W. W. Wood (1991), Large-scale natural gradient tracer test in sand and gravel, Cape
Cod, Massachusetts: 1. Experimental design and observed tracer movement. Water
Resources Research, 27: 895–910.
Mackay, D. M., D. L. Freyberg, P. V. Roberts, and J. A. Cherry (1986), A natural gradient
experiment on solute transport in a sand aquifer: 1. Approach and overview of plume
movement. Water Resources Research, 22: 2017–2029.
Meerschaert, M. M., M. Dogan, R. L. Van Dam, D. W. Hyndman, and D. A. Benson (2013),
Hydraulic conductivity fields: Gaussian or not? Water Resources Research, 49,
doi:10.1002/wrcr.20376.

87

Molz, F. J.,O. Güven, J. G. Melville, R. D. Crocker, and K. T. Matteson (1986), Performance,
analysis, and simulation of a two-well tracer test at the mobile site. Water Resources
Research, 22: 1031–1037.
Phanikumar, M. S., D. W. Hyndman, X. Zhao, and M. Dybas (2005), A three-dimensional model
of microbial transport and biodegradation at the Schoolcraft, Michigan site. Water Resources
Research, 41: W05011.
Sudicky, E. A. (1986), A natural gradient experiment on solute transport in a sand aquifer:
Spatial variability of hydraulic conductivity and its role in the dispersion process. Water
Resources Research, 22: 2069–2082.
Whittaker, J., and G. Teutsch (1999), Numerical simulation of subsurface characterization
methods: application to a natural aquifer analogue. Advances in Water Resources, 22: 819829.
Zheng, C., and S. M. Gorelick (2003), Analysis of solute transport in flow fields influenced by
preferential flowpaths at the decimeter scale. Ground Water 41: 142–155.
Zheng, C., and P. Wang (1999), MT3DMS: A modular three-dimensional multispecies transport
model for simulation of advection, dispersion, and chemical reactions of contaminants in
groundwater systems; documentation and user's guide.
Zimmerman, D. A., G. de Marsily, C. A. Gotway, M. G. Marietta, C. L. Axness, R. L.
Beauheim, R. L. Bras, J. Carrera, G. Dagan, P. B. Davies, D. P. Gallegos, A. Galli, J.
Gómez-Hernández, P. Grindrod, A. L. Gutjahr, P. K. Kitanidis, A. M. Lavenue, D.
McLaughlin, S. P. Neuman, B. S. RamaRao, C. Ravenne, and Y. Rubin (1998), A
comparison of seven geostatistically based inverse approaches to estimate transmissivities for
modeling advective transport by groundwater flow. Water Resources Research, 34: 13731413.

88

Chapter 5
Quantifying the value of different data sets and modeling schemes for flow and transport
simulations
5.1 Introduction
Depending on the level of heterogeneity, flow and transport modeling through aquifers
can be an extremely challenging task. Previous research has shown that improvement of
modeling results can be achieved in multiple ways, including: (1) different modeling approaches
to account for some of the behavior observed in tracer experiments, (2) a better knowledge of the
(statistical) distribution of hydraulic conductivity (K) by collecting more or better field data, and
(3) addition of complementary data that give information on the structural characteristics of the
aquifer.
Regarding the different modeling approaches, when the classic ADE approach does not
produce acceptable results, alternatives have been proposed. These approaches include, amongst
others, preferential flow paths [Zheng and Gorelick, 2003] and mass transfer between mobile and
immobile domains [Harvey and Gorelick, 2000]. The difficulty with these approaches is that
they are difficult to parameterize based on field data and are thus non-predictive.
Regarding the statistical distribution of subsurface properties, direct measurements of K
have traditionally been done using flowmeters. This is a common, easily accessible, and
relatively cheap method that requires installation of a groundwater well. This method is not
sensitive to low K, but is very capable of measuring high K values [Bomana et al., 1997].
Disadvantages of this method include an unknown support volume, sensitivity to well bore
disturbance, and low vertical resolution. In recent years a new method, the High Resolution

89

Hydraulic Conductivity (HRK) tool has been developed, which allows measurements of K for
smaller support volumes and higher vertical resolution than flowmeters. This method, which is
sensitive to low K, but cannot yet measure K above 60 m/d, has been discussed in considerable
detail in the previous chapters and in Liu et al. [2009, 2012]. A disadvantage of this method is
that it is not yet not widely accessible and relatively expensive.
Many researchers have shown that complimentary data can greatly improve modeling
results. These data can be derived from geophysical measurements, such as electrical resistivity
[e.g., Atekwana et al., 2000; Cassiani et al., 2006], nuclear magnetic resonance (NMR)
measurements [Legchenko et al., 2002], seismic [Hyndman, et al., 1994; Hyndman and Gorelick,
1996], and surface and cross borehole ground-penetrating radar (GPR) [e.g., Tronicke et al.,
2002]. GPR is very sensitive to textural changes, porosity and water content, and can be used to
obtain structural information of an aquifer; GPR can thus provide complementary data to direct
measurements of K.
Of the different methods presented above, each has different levels of detail, cost, and
benefit. Previous chapters in this dissertation have discussed the relative benefits of various
alternative approaches and added data sets for modeling flow and transport at the MADE site.
The objective of this chapter is to compare the various additions or improvements to modeling
the flow and transport behavior at the MADE site, as observed during the MADE-1 tracer test
[Boggs, 1991]. The comparisons presented here are only a subset of possibilities, and are based
on the available data in the injection area (Figure 5.1).
The first comparison of modeling results in this chapter is between simulations through K
fields parameterized with HRK data and through fields parameterized with a combination of
HRK data and GPR facies. This section is a natural follow up on Chapter 4. The next section of

90

this chapter is a cost-benefit analysis of different options, which focuses on (1) the effects of
transient simulation (of the water table), (2) a K field parameterized with only flowmeter data
measured during the MADE studies [Boggs, 1991, 1992], (3) the amount of HRK information,
and 4) GPR facies. All results are compared with the MADE-1 measurements and a single base
simulation that used steady-state flow and parameterization based on 25 HRK profiles (this base
simulation was discussed in Chapter 4 – Figure 4.3c).

Injection wells
MLS wells
Flowmeter K wells
HRK profiles
45
40

GPR cube
24
25

35

22

Distance [m]

18

25
17

15

16

11

20

10

10

5

5
0

19

20

21

30

(b)
23

15
12
9
6
4

14

13

8

7

3

2

1
0

5

10
15
20
Distance [m]

25

Figure 5.1 A detailed map of the model domain with locations of flowmeter K (orange) and
HRK measurements (blue) and multi-level sampler (grey) wells. The area of 3D GPR data is
given by the rectangle with green dashed line. Injection wells are given with red stars. See Figure
4.1 for the location within the larger MADE site.

91

5.2 Complementary GPR data
In this research, 3D full-resolution GPR data were used as complementary data to
measurements of K using the HRK tool. GPR data were collected with 50 and 100 MHz antennas
using a trigger wheel and ropes for positional guidance (details in Chapter 2). The profile
separation was selected as equal to step size and smaller than a 1/4 of the GPR signal wavelength
[Grasmueck et al., 2005]. Data were processed using a de-wow function and a low pass filter.
Time-depth conversion was made using the velocities derived from common-mid-point and cross
borehole measurements. An automated search algorithm was used to pick the reflectors.
Reflectors were then grouped using another automated algorithm to create the facies boundaries.
Distance [m]
0
10
20
(a)

Distance [m]
0
10
20
(b)

Distance [m]
0
10
20
(c)

40

40

30

30

30

20

20

20

10

Distance [m]

40

10

10

0
Relative mass (x10-3 )
1

2

3

0
Relative mass (x10 -3)
2 4

6 8 10

0

Relative mass
0.01

0.02

Figure 5.2 Vertically integrated contour maps of relative mass; (a) MADE-1 experiment at
t=503 days, (b) mean of 6 simulations, (c) best head change simulation. K fields for the
simulations in (b) and (c) were based on a combination of HRK data and GPR facies (these maps
can be compared with Figure 4.2, which presents contour maps of relative mass for K fields
based on just HRK data).
The K field simulation procedure was similar to the one for the scenario with HRK data
alone (see Chapter 4), but now separated by facies (see Chapter 3). For each facies, K fields were

92

generated using appropriate depth interval from the same noise field. K fields for each facies
were then stitched together to create a 3D model for the flow and transport simulations in
MODFLOW and MT3D [Zheng and Wang, 1999]. From the simulation results, 2D vertically
integrated contour maps and 1D concentration profiles were created following the same
procedure as in Chapter 4.

(a)

0.08

0.02

0.06

0.015

0.04

0.01

0.02

0.005

0

0

5

10

15

20
25
Distance [m]

30

35

0
45

40

Relative mass

0.1

0.025

(b)

0.08

0.02

0.06

0.015

0.04

0.01

0.02

0.005

0

Relative mass

0.025

0

5

10

15

20
25
Distance [m]

30

35

40

Relative mass

Relative mass

0.1

0
45

Figure 5.3 Relative mass distribution profiles for the MADE-1 experiment at t=503 days (red
line), and the mean and 2σ range of 6 simulations (black line and gray shaded area, respectively).
Simulations used (a) K fields based on just HRK data and (b) K fields based on a combination of
HRK data and GPR facies.
Similar as with the K fields based on just the HRK data (Chapter 4) in these simulations
the recovered mass is significantly larger than the mass recovered during the experiment.
However, the extent of the plume was very comparable to the measured data (Figure 5.2) and the
location of maximum relative mass was similar to what was observed (Figure 5.3). Most

93

importantly, the facies-based approach resulted in a much better characterization of the plume
tail (Figure 5.3b). This improvement is because the facies boundaries allow one to determine the
high-K zones of the aquifer in more detail. In summary, this comparison shows that the
simulations that incorporate structural information from GPR-derived hydrofacies better
reproduce the part of the tracer that rapidly moved down-gradient.

5.3 Cost-benefit analysis
The objective of this section is to compare the relative benefits and disadvantages of
different modeling approaches, use of flowmeter K data, different densities of HRK information,
and addition of a complementary data set (GPR facies). The results are compared to the MADE1 experiment data and to the base simulation, which had the best match to the head-change
criterion (Figure 4.3b). This base simulation used steady state flow, all 25 HRK data for
conditioning, and no flowmeter data or GPR facies. For all compared K fields, the same noise
field was used.
The results are presented qualitatively using 2D vertically integrated mass distribution
maps (contour maps) and 1D mass distribution profiles. Additionally, for quantitative
comparison of the results, Table 5.1 lists the details of each simulation (down-gradient distance
to maximum mass peak, peak and tail shapes, injected mass, and mass recovery for each
simulation), and compares the RMS error between simulation and field data. All contour maps
and 1D concentration profiles have been normalized by total mass so that the results can be more
easily compared.

94

Table 5.1 Quantitative measures for the comparison of simulation scenarios and field experiment
2
(R values were calculated comparative to the MADE-1 experiment t=503 days in the first row).
2

R
Peak
Peak
Tail
loc. [m] Total
(0-45m) (0-25m) (25-45m)
Data
(MADE-1, t=503
days)

11.5

1

1

1

Base simulation

11.5

0.899

0.880

0.101

Transient flow

11.5

0.879

0.843

0.000

Flowmeter K

19.5

0.002

0.001

0.002

GPR data

11.5

0.830

0.783

0.760

80% HRK

11.5

0.924

0.895

0.445

60% HRK

11.5

0.893

0.874

0.102

40% HRK

8.5

0.397

0.236

0.267

20% HRK

8.5

0.712

0.624

0.049

Mass [grams]
Injected

Simulated

Recovered

25000.00
(100 %)

N/A

5326.5
(23.7 %)

25183.49
(100 %)
24917.69
(100 %)
25173.43
(100 %)
25181.85
(100 %)
25161.81
(100 %)
25177.01
(100 %)
25161.48
(100 %)
25174.02
(100 %)

23722.07
(94.2 %)
14342.35
(57.5 %)
12193.20
(48.4 %)
24634.42
(97.8 %)
19296.71
(76.7 %)
23045.80
(91.5 %)
21670.72
(86.1 %)
25111.05
(99.7 %)

21969
(87.2 %)
21899.35
(87.9 %)
14469
(57.5 %)
24363
(96.7 %)
18583.3
(76.8 %)
21632.7
(85.9 %)
25013.7
(99.4 %)
27940.6
(110.9 %)

5.3.1 Transient flow
Since the flow-transport simulations are computational experiments, simplification is
required. Simplification is also beneficial due to lower computational costs. In hydrologic
simulations it is common to use steady state solutions. However, under natural conditions, flow
and transport is more complicated than is typically achieved by simplified models. Unconfined
aquifers are especially challenging to model with precision, since such aquifers are usually very
shallow, highly heterogeneous and open to recharge through the surface via infiltration and
discharge. Thus, flow and transport simulations in unconfined aquifers may require transient
solutions to reproduce the effect of changing water table depth. Transient solutions require

95

periodically collected piezometric measurements to specify the changes in water table depth
along boundaries of models. Computational power and time for transient simulations are also
more intense than for steady state simulations.
Comparative analysis of flow-transport simulations involving transient and steady state
flow schemes are presented in this section. The objective of this comparison is to test the effect
of possible tracer entrapment in vadose zone due to water level changes, as was hypothesized by

Adams and Gelhar [1992]. The transient flow simulation was based on the piezometric
measurements collected during the MADE-1 experiment. This simulation used the same K field
as the base simulation; only the head conditions were changed.
The results shown in the contour maps of normalized mass (Figure 5.4) and 1D
concentration profiles (Figure 5.5) suggest that the transient flow had very little effect. The
transient flow simulations were almost identical to the ones using steady state flow. It is
hypothesized that this small difference is in part due to the fact that the injection screen was
located at a depth of 7.4-8 meters (around 3.8 meters below the depth of the water table at the
start of the experiment and simulation). Considering that the transient simulation took
significantly more computer time and required the additional piezometric data, this modification
of the approach presents very little additional value.
Quantitative measures presented in Table 5.1 show that the transient flow simulation is
very similar to the base simulation, except the percentage of mass that stayed in the model
domain after 503 days of simulation (57.5 % versus 94.2% for the base simulation). Higher
gradient periods during the flow simulation may be the cause of this difference since the seasonal
changes in the water table promote faster flow in the model domain and more mass leaves
through the northern boundary of the domain. Although mass entrapment in vadose zone could

96

be another reason, it was not captured by the simulation and numerical sampling procedure of the
simulated plume.

Distance [m]
0
10
20
(a)

Distance [m]
0
10
20
(b)

Distance [m]
0
10
20
(c)

40

40

30

30

30

20

20

20

10

Distance [m]

40

10

10

0
0
0
Normalized mass (m/mrec)Normalized mass (m/mrec ) Normalized mass (m/mrec )
0.01

0.02

0.01

0.02

0.01

0.02

Normalized mass (m/m )
rec

Figure 5.4 Vertically integrated contour maps of normalized mass; (a) MADE-1 experiment at
t=503 days, (b) base simulation, (c) transient flow simulation.

0.1

data
steady state flow
transient flow

0.08
0.06
0.04
0.02
0

0

5

10

15

20
25
Distance [m]

30

35

40

45

Figure 5.5 Normalized mass distribution profiles of the MADE-1 experiment at t=503 days (red
line), base simulation (black line), and transient flow scheme (blue line).

97

5.3.2 Flowmeter K
Direct measurements of K have traditionally been done using flowmeters. At the MADE
site many flowmeter measurements were conducted during the MADE field experiments in the
1980s and 1990s (see Figure 5.1 for the locations of flowmeter wells in the study domain).
Flowmeter measurements in heterogeneous media are not very sensitive to low K, but very
sensitive to high K zones [Bomana et al., 1997]. Disadvantages of this method include the
unknown support volume, low vertical resolution, and averaging over a bulk volume of different

K zones.
Here, a direct comparison is given between simulations through K fields based on (1)
flowmeter data, and (2) the more recently introduced HRK tool. The stochastic K fields were
created using the same fractal methods discussed in Chapter 4 with the flowmeter K data as hard
conditioning points.
The results show that for the simulation based on the flowmeter data (Figure 5.6c and
5.7), the mass moves much faster down-gradient than for the simulations based on the HRK data
(Figure 5.6b and 5.7). In fact, a significant portion of the tracer mass moved outside the model
domain (Table 5.1). Flowmeter measurements are not sensitive to the low K, which explains why
the average K is higher than HRK.
Quantitative measures provided in Table 5.1 suggest that the flowmeter K based
simulation cannot reproduce the observed plume shape at the MADE site. This simulation
resulted in the most incomparable plume shape, peak location, and tail behavior of all scenarios
tested.

98

Distance [m]
0
10
20
(a)

Distance [m]
0
10
20
(b)

Distance [m]
0
10
20
(c)

40

40

30

30

30

20

20

20

10

Distance [m]

40

10

10

0
0
0
Normalized mass (m/mrec)Normalized mass (m/mrec ) Normalized mass (m/mrec )
0.01

0.02

0.01

0.02

0.01

0.02

Normalized mass (m/m )
rec

Figure 5.6 Vertically integrated contour maps of normalized mass; (a) MADE-1 experiment at
t=503 days, (b) base simulation, (c) simulation based on just flowmeter K data.

0.1

data
HRK
flowmeter K

0.08
0.06
0.04
0.02
0

0

5

10

15

20
25
Distance [m]

30

35

40

45

Figure 5.7 Normalized mass distribution profiles of the MADE-1 experiment at t=503 days (red
line), base simulation (black line), and simulation based on just flowmeter K data (blue line).

5.3.3 Effect of conditioning data density
A random number generator was used to create subsets of HRK measurement profiles in
the modeling domain. The number of HRK profiles used for the K field parameterization was
decreased from 25 (the base simulation) to 20 (80%), 15 (60%), 10 (40%), and 5 (20%). HRK

99

profile #1, which is located very near to the injection location (see Figure 5.1) was removed from
each of parameterized fields. The same noise field was used to generate each K field.
A visual comparison of 2D vertically integrated mass distribution maps (Figure 5.8c-f)
shows that decreasing the number of HRK profiles has a significant effect on both the plume
extent and shape. Although the locations of randomly selected HRK profiles have an effect on
the results, it is clear that the accuracy of modeled plumes generally decreases with number of
HRK profiles. As expected, the plume simulated using 80 % of the HRK profiles was most
comparable to the base simulation.
The effects of reducing the HRK profile density on the peak and tail shapes of simulated
plumes is also clear from the 1D mass distribution profiles (Figure 5.9). HRK profiles located
near the injection location control the quality of peak representation. The simulation with 40% of
HRK profiles has no nearby HRK profiles (Figure 5.8e), which resulted in a highly distorted
peak shape (Figure 5.9c). It is also notable that although 1D mass distribution profiles (in this
case: 60% HRK profiles; Figure 5.9b) may look nearly unchanged from the base simulation and
the situation with more HRK profiles, the corresponding integrated mass distribution map
(Figure 5.8d) shows a strong deviation from the base simulation.
The simulations with 20 and 40 % of the HRK profiles provided poor correlation with the
base simulation results for both the peak and the tail shape of the mass distribution profiles
(Table 5.1). The mass recoveries for these simulations also significantly deviated from the base
simulation.

100

Distance [m]

40

Distance [m]
0
10
20
(a)

40

Distance [m]
0
10
20
(b)

40

30

30

30

20

20

20

10

10

Distance [m]
0
10
20
(c)

10

0
0
0
Normalized mass (m/mrec )Normalized mass (m/mrec) Normalized mass (m/mrec )

Distance [m]

40

0.01
0.02
Distance [m]
0
10
20
(d)

40

0.01
0.02
Distance [m]
0
10
20
(e)

40

30

30

30

20

20

20

10

10

0.01
0.02
Distance [m]
0
10
20
(f)

10

0
0
0
Normalized mass (m/mrec )Normalized mass (m/mrec) Normalized mass (m/mrec )
0.01

0.02

0.01

0.02

0.01

0.02

Figure 5.8 Vertically integrated contour maps of normalized mass: (a) MADE-1 experiment at
t=503 days, (b) base simulation based on all 25 HRK profiles, and simulations based on (c) 20 ,
(d) 15, (e) 10, and (f) 5 randomly selected HRK profiles. The randomly selected HRK profiles
used to generate the K fields are shown with blue symbols; see Figure 5.1 for the location of all
25.

101

Norm. mass (m/mrec) Norm. mass (m/mrec) Norm. mass (m/mrec) Norm. mass (m/mrec )

0.1 (a)

data
HRK 100%
HRK 80%

0.05

0
0.1 (b)

data
HRK 100%
HRK 60%

0.05

0
0.1 (c)

data
HRK 100%
HRK 40%

0.05

0
0.1 (d)

data
HRK 100%
HRK 20%

0.05

0

0

5

10

15

20
25
Distance [m]

30

35

40

45

Figure 5.9 Normalized mass distribution profiles of the MADE-1 experiment at t=503 days (red
line), base simulation (black line), and simulations based randomly selected number of HRK
profiles (blue line). (a) 80% of HRK profiles (#’s 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 14, 15, 16, 17, 18,
19, 20, 22, 23, 24), (b) 60% of HRK profiles (#’s 2, 3, 4, 8, 10, 11, 12, 13, 14, 16, 17, 18, 20, 21,
24), (c) 40% of HRK profiles (#’s 2, 5, 7, 11, 12, 14, 17, 18, 19, 22) (d) 20% of HRK profiles
(#’s 8, 11, 14, 15, 24). See Figure 5.1 for the locations of the HRK profiles.

102

5.3.4 Complementary GPR data
To allow comparison with the previous sections, one simulation for a K field partially
based on GPR facies is presented. Rather than using the average of 6 simulations (Section 5.2),
only the simulation for the noise field that best matched the head-change criterion is presented
here.
As already discussed in Section 5.2, when GPR facies are used to help construct the K
field, the mass contour maps do not differ much from those based on just HRK data (Figure
5.10). However, the tail of the tracer plume is much better characterized (Figure 5.11).
Distance [m]
0
10
20
(a)

Distance [m]
0
10
20
(b)

Distance [m]
0
10
20
(c)

40

40

30

30

30

20

20

20

10

Distance [m]

40

10

10

0
0
0
Normalized mass (m/mrec )Normalized mass (m/mrec) Normalized mass (m/mrec )
0.01

0.02

0.01

0.02

0.01

0.02

Figure 5.10 Vertically integrated contour maps of normalized mass: (a) MADE-1 experiment at
t=503 days, (b) base simulation, and (c) simulation based on a combination of HRK data and
GPR facies.

Correlation coefficients provided in Table 5.1 shows that the simulation with GPR data
exhibits the best representation of heavy tail behavior observed in the mass distribution curves.
However, the peak representation is not quite as good as the base simulation. This can be

103

explained by the lack of GPR data around the injection area (Figure 5.1). As mentioned earlier in
Chapter 1, the source area was dug out during the MADE-3 macro scale experiment in the 1990s.
Therefore, 3D GPR data were not collected in the vicinity of the source wells, which limits the

Normalized mass (m/m )
rec

effective use of GPR-derived facies boundaries around the injection area.

0.1

data
w/o GPR
w/ GPR

0.08
0.06
0.04
0.02
0

0

5

10

15

20
25
Distance [m]

30

35

40

45

Figure 5.11 Normalized mass distribution profiles of the MADE-1 experiment at t=503 days
(red line), base simulation (black line), and simulation based on a combination of HRK data and
GPR facies.
5.4 Conclusions

The first section of this chapter introduces the use of GPR data in combination with HRK
data and presents the results of simulations for both types of K fields. The mean of the flow and
transport simulations involving the GPR data exhibits a better representation of the tailing
behavior. However, for the simulations based on the facies boundary information the standard
deviation of the simulation means was somewhat higher than for the simulations without facies.
This might be explained by the fact that the criterion used to select the noise fields was based on
a maximum head-change (Chapter 4). Here, the same noise fields were used without taking into
account any effects on the head changes (due to a low K facies around the injection depth, they
increased from a mean of 0.67 m to 1.03 m with a larger standard deviation).

104

In the second section of this chapter, the effects of complementary data, transient flow,
use of conventional flowmeter K data, and different densities of HRK data were examined.
Based on these comparisons, the following conclusions can be drawn:

•

Good representations of peak behavior were obtained in the base simulation, in transient
flow simulations, and in the simulation with the aid of GPR data.

•

Flowmeters are the most common tool to collect K information in field investigations, but
the results presented here show that they are incapable of representing the K field in
sufficient detail, resulting in wrong plume behavior in the transport simulations.

•

The best representation of tail behavior was obtained using GPR facies as
complementary data. This outcome shows that in cases where heavy tail behavior exists,
GPR facies boundaries can help identify the preferential flow paths through an aquifer.

•

Reducing the density of HRK profiles resulted in less comparable simulation results. In
this case, 60% of HRK profiles were required to obtain reasonable 1D mass profiles, but
80% were required for a good 2D mass distribution. However, it is important to point out
that the locations of the selected HRK profiles had a critical effect on the simulation
results. Omission of the HRK profiles close to the source area had a strong negative
effect on representation of the peak behavior.

105

REFERENCES

106

REFERENCES

Adams, E. E., and L. W. Gelhar (1992), Field study of dispersion in a heterogeneous aquifer: 2.
Spatial moments analysis. Water Resources Research, 28(12), 3293–3307,
doi:10.1029/92WR01757.
Atekwana, E. A., W. A. Sauck, and D. D. Werkema (2000), Investigations of geoelectrical
signatures at a hydrocarbon contaminated site. Journal of Applied Geophysics, 44: 167-180.
Boggs, J. M., S. C. Young, L. M. Beard, L. W. Gelhar, K. R. Rehfeldt, and E. E. Adams (1992),
Field study of dispersion in a heterogeneous aquifer: 1. Overview and site description. Water
Resources Research, 28, 3281-3291.
Boggs, J. M. (1991), Database for the first macrodispersion experiment (MADE-1).
Bomana, G. K., F. J. Molz, and K. D. Boonec (1997), Borehole flowmeter application in fluvial
sediments: methodology, results, and assessment. Ground Water, 35: 443–450. doi:
10.1111/j.1745-6584.1997.tb00104.x.
Cassiani, G., V. Bruno, A. Villa, N. Fusi, A. M. Binley (2006), A saline trace test monitored via
time-lapse surface electrical resistivity tomography. Journal of Applied Geophysics, 59: 244259.
Grasmueck, M., R. Weger, and H. Horstmeyer (2005), Full-resolution 3D GPR imaging.
Geophysics, 70: K12-K19.
Harvey, C., and S. M. Gorelick (2000), Rate-limited mass transfer or macrodispersion: Which
dominates plume evolution at the macrodispersion experiment (MADE) site? Water
Resources Research, 36: 637–650.
Hyndman, D. W., J. M. Harris, and S. M. Gorelick (1994), Coupled seismic and tracer test
inversion for aquifer property characterization. Water Resources Research, 30.7: 1965-1977.
Hyndman, D. W., and S. M. Gorelick (1996), Estimating lithologic and transport properties in
three dimensions using seismic and tracer data: The Kesterson aquifer. Water Resources
Research, 32.9: 2659-2670.
Legchenko, A., J. M. Baltassat, A. Beauce, and J. Bernard (2002), Nuclear magnetic resonance
as a geophysical tool for hydrogeologists. Journal of Applied Geophysics, 50: 21-46.
Liu, G., J. J. Butler Jr., G. C. Bohling, E. Reboulet, S. Knobbe, and D. W. Hyndman (2009), A
new method for high-resolution characterization of hydraulic conductivity, Water Resources
Research, 45(8), doi:10.1029/2009WR008319.
Liu, G., J. J. Butler Jr., E. Reboulet, and S. Knobbe (2012), Hydraulic conductivity profiling with
direct-push methods, Grundwasser, 17(1), 19–29.

107

Tronicke, J., P. Dietrich, U. Wahlig, E. Appel (2002), Integrating surface georadar and crosshole
radar tomography: A validation experiment in braided stream deposits. Geophysics, 67:
1516-1523.
Zheng, C., and S. M. Gorelick (2003), Analysis of solute transport in flow fields influenced by
preferential flowpaths at the decimeter scale. Ground Water 41: 142–155.
Zheng, C., and P. Wang (1999), MT3DMS: A modular three-dimensional multispecies transport
model for simulation of advection, dispersion, and chemical reactions of contaminants in
groundwater systems; documentation and user's guide.

108

Chapter 6
Synthesis

This dissertation presents a set of techniques to improve the characterization of and flow
and transport simulations through highly heterogeneous aquifers. The Macro Dispersion
Experiment site in Mississippi was selected to test the developed methods due to existing
extensive data sets and previous research. The biggest challenge for this site is the highly
heterogeneous nature of the porous media and the observed non-Gaussian transport during the
large scale tracer tests in 1980s and 1990s [Boggs, 1991; Boggs et al., 1992]. Flow meter K data,
collected prior to those tracer experiments, suggested that the K distribution was highly
heterogeneous. In previous studies, this behavior could not be successfully modeled using
advection and dispersion as the primary mechanisms. Indeed, this extraordinary behavior led
most of the hydrological research community to believe that the advection-dispersion equation
(ADE) is not capable of modeling the flow and transport in aquifers above a certain
heterogeneity level.
This study was aimed to test the hypothesis that ADE is capable of modeling flow and
transport even in highly heterogeneous aquifers as long as the model is supported by sufficient
high resolution K data. The required high resolution K fields were created using fractal stochastic
methods based on two high-resolution data sets: high-resolution direct-push K profiles and 3D
full-resolution ground-penetrating radar (GPR). Below, I present a brief summary of the main
outcomes of each Chapter, followed by a discussion of some new questions resulting from this
work.
GPR, a high resolution non-invasive geophysical method, is capable of defining textural
properties such as porosity, packing, and sorting in a spatial context. These textural properties

109

are related to K in an indirect way. Chapter 2 [based on Dogan et al., 2011] presents an
integrative study to investigate the consistency between GPR derived structures (or facies) and K
data. This work was performed in an area where a recent push pull test was performed.
Descriptive statistics and distribution tests performed in Chapter 2 proved that the different
facies separated statistically different K distributions, with a reduced level of K variation within
each facies. Therefore, this study demonstrated that GPR-derived structural information can be
successfully used to define hydrofacies.
High-resolution spatial information is necessary to create better flow and transport
models, as hypothesized above. However, no invasive or non-invasive method exists to measure

K fields at a spatial resolution that does not require some form of interpolation. Therefore,
stochastic methods are generally required to generate parameter fields to fill in the gaps between
measurement locations with similar levels of heterogeneity. To avoid drawbacks of common
stochastic methods, including the inability to create parameter fields that display the connectivity
of real sedimentary deposits, fractal stochastic methods were used in this dissertation.
Connectivity is an important property of aquifers that can define connected high K paths as
preferential flow paths. Chapter 3 [based on Meerschaert et al., 2013] presents the details of the
fractal K field generation technique that was developed and shows that the fractional differencing
filter can remove the long-range dependence observed in K data. This finding shows that
fractional differencing filters can be successfully used to further investigate the underlying
distribution of K data. Parameter fields are presented that are more natural looking and have a K
distribution that better matches observations from HRK profiles. Flow simulations through these
fields resulted in improved representation of the transport characteristics typical for this site,
such as a heavy tail and multiple concentration peaks.

110

The focus of the next two Chapters was on the area around the original injection wells of
the Macro Dispersion Experiment tracer tests. Also, the approach demonstrated in Chapter 2 was
applied to a 3D domain. Flow and transport simulations were performed following injection
procedure and sampling strategy of the MADE-1 experiments. Several different K fields were
generated using the method explained in Chapter 3. The concentrations derived from the
simulations, excluding GPR facies to derive the parameter fields, were presented in Chapter 4.
Vertically integrated contour maps of relative mass distribution created in this study were an
innovation compared to earlier modeling efforts at the MADE site, as no other study has
compared simulation results with experiment data in more than 1D. The extent of the modeled
plumes and maximum mass zones were very comparable to the ones measured during the
experiment. The simulations also did a good job of characterizing the fast moving heavy tail of
the plume. These findings show that the proposed method using K fields based on just measured
field data is able to define the parameters necessary for designing effective remediation solutions
for real-world problems.
Chapter 5 presents a comparative study to define the value of different data sets for
improving modeling results and provides insight into the required data sets for successful flow
and transport modeling in heterogeneous media. This Chapter emphasizes the effect of
complementary GPR data to better reproduce the heavy tail behavior of MADE plume. It also
highlights the inability of flow meter K data to produce simulation results comparable to the
observed data. It presents a comparison of the steady state and transient flow schemes, which
concludes that the tracer entrapment in the vadose zone is not a strong argument to explain the
poor mass recovery of the modeling efforts. Finally, an analysis of the effects of HRK

111

conditioning density shows that 60% to 80% of HRK profiles are needed to produce acceptable
simulation results.
Despite the new developments and improvements of modeling of flow and transport in
heterogeneous media presented in this dissertation, there are a few unanswered questions.
Several of these points can be summarized as follows, and may be answered with future work:

•

As was shown in this dissertation, GPR facies that are derived from surface
measurements are very useful to build hydrostratigraphic frameworks that can improve
flow and transport modeling (Chapters 2 and 5). However, facies boundaries cannot be
directly used to obtain K information. Future work could make use of inversion
approaches and interpretation of borehole measurements to develop a more direct link
between geophysical data and hydraulic properties of interest.

•

The 3D simulations had excellent correlation with the transport observations (Chapter 4
and 5) and fractal methods provided a better tool to generate K fields than sequential
Gaussian methods (Chapter 3). One of the two main components of fractal stochastic
parameter field generation is the parameterization (fractal dimension, correlation
structures and statistical measures) of the fractal method, which can be obtained using
HRK data. However, the second component, the noise field, is randomly generated and of
Gaussian nature. This Gaussian noise field may not be the best solution to this step, but
further work would be needed to investigate whether this simplification has important
consequences.

•

One major disagreement between the simulations and the field experiments involves the
mass recovery. All simulations had significantly larger mass recovery than the
experiments, despite identical sampling and interpolation procedures. One possibility for

112

this observation has been considered in Chapter 5 (transient flow), but this did not
address the discrepancy. Other possible explanations include the lower mean K
throughout the model domain due to the upper of limit of the HRK tool, and an inability
to reproduce the mixing in sampling wells as would have occurred during the sampling
procedure in the field.
Finally, this study demonstrated that the ADE is capable of modeling flow and transport even in
highly heterogeneous aquifers such as the MADE site without parameter calibration or
optimization. However, high resolution parameter fields are required to enable these modeling
efforts. The novel high resolution characterization tool HRK can be used in combination with
high resolution geophysical methods to provide the necessary data. Even though some steps of
the presented research need improvement, it presents a straightforward method to create
uncalibrated, reproducible flow and transport models based on solely field data. Moreover, these
models can reproduce the non-Gaussian transport behavior at the infamous MADE site.

113

REFERENCES

114

REFERENCES

Dogan, M., R. L. Van Dam, G. C. Bohling, J. J. Butler Jr, and D. W. Hyndman (2011),
Hydrostratigraphic analysis of the MADE site with full-resolution GPR and direct-push
hydraulic profiling, Geophys. Res. Lett., 38(6), doi:10.1029/2010GL046439.
Meerschaert, M. M., M. Dogan, R. L. Van Dam, D. W. Hyndman, and D. A. Benson (2013),
Hydraulic conductivity fields: Gaussian or not? Water Resources Research, 49,
doi:10.1002/wrcr.20376.
Boggs, J. M., S. C.Young, L. M. Beard, L. W. Gelhar, K. R. Rehfeldt & E. E. Adams (1992),
Field study of dispersion in a heterogeneous aquifer: 1. Overview and site description. Water
Resources Research, 28(12), 3281-3291. doi:10.1029/92WR01756.
Boggs, J. M. (1991), Database for the first Macrodispersion Experiment (MADE-1).

115