4 ‘ . . , . . ‘ . ‘ .vH-Y‘lin &

1.3..“

3 t
3.... 1.1.x?
‘1“... K4323

i231 - . .
.. $fi¢t§io

sxii...

.. in".
.7 1:. :
:2 1:!

    

. r:
a: n ‘ aw...
‘ $64.1”? . v.31“ , . . .:n .. 1
V . . . . u: ‘ Z . v . .
. .1 k)? .ty-INI . .
:. m3: ‘
I [y '1“ II

 

* LIBRARY
Michigan State
University

 

 

This is to certify that the
thesis entitled

Interactive Ground Water (IGW): Theory and
Illustrative Applications

presented by

Rebekah J. Stephenson

has been accepted towards fulfillment
of the requirements for the

MS. degree In Civil Engineering

 

g7 QR;-

\_,Major Profess r's Signature

 

Mg 5‘ Lab-1'15”-

Date

MSU is an Attir'mative Action/Equal Opportunity Institution

.. .
- ---.-o--n-u-o-o-a-a-u--~---o-o-c-c-c-o-n-o-n-.-o-.-.-.-.--.--._.-n-o—o-.----o-o--n-o-s-¢-a-u-o-a-o--o--4-o--

PLACE IN RETURN sox to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2/05 p:/ClRC/DateDue.indd-p.1

 

 

INTERACTIVE GROUND WATER (IGW): THEORY AND
ILLUSTRATIVE APPLICATIONS

By

Rebekah J. Stephenson

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTERS OF SCIENCE
Department of Civil Engineering

2005

Abstract

INTERACTIVE GROUND WATER (IGW): THEORY AND
ILLUSTRATIVE APPLICATIONS

By

Rebekah J. Stephenson

The software package Interactive Ground Water (IGW) forms a “new paradigm for
groundwater modeling”, in that; both visually and computationally the dynamics of the
real world are replicated. In hopes of bridging the bottleneck of the computational
requirements for a 3-D model, we exploit the 2-D version of the software package IGW
using the “new paradigm”. The driving algorithms and mathematical theory used in IGW
are described and coupled with visual results. In particular, the iterative flow equation
and theory used in IGW are discussed in part to validate the software capabilities and use
of the "new paradigm." We show results for both deterministic and stochastic
hypothetical flow and contaminant transport through examples that are considered novel
in groundwater theory. In addition, we expand on the novel examples and show dynamic
examples in which the internal boundaries of a lake and seepage zone change over time.
We also study the Monte Carlo method and give an example in which the concentration
plume is non-Gaussian due to the internal boundary of a lake that is connected to the

aquifer.

Copyright by
Rebekah Jane Stephenson
2005

‘lll' Ill I‘tl

"l|f\ i

To my loving parents Robert and Mary Jo Stephenson

iv

ACKNOWLEDGENIENTS

Groundwater has fascinated me ever since I became aware of this field of study. I like
predicting and analyzing groundwater because it is a vital source that for the most part is
not visible. Groundwater is a very special sensitive resource located beneath the surface.
In a similar fashion, I have many family and friends who have supported me while
working on this thesis that I wouldn't want to go unnoticed. I would first like to thank my
advisor Dr. Shu-Guang Li for his tremendous support and dedication to this thesis as well
as guiding me in a new direction of scientific thought. Dr. Li first captured my attention
with his software Interactive Ground Water (IGW) with an example of contaminants
entering a well. I have been intrigued ever since and would like to continue to learn more
about groundwater. In addition, I would like to thank my committee members Dr.
Phanikumar Mantha and Dr. Roger Wallace for their support and for their guidance at
critical moments in this process. They have both been great instructors and have
encouraged me in my studies.

Second, I would like to thank my mother and father, Mary Jo and Robert Stephenson,
for being patient and supportive of my dreams. I also greatly appreciate the
conversations that we have about my studies in school and I commend them for always
listening. I would also like to thank Sonya and Eric Snyder, Jody Stephenson, and
James Stephenson for their continual family support and understanding.

I would also like to thank my grandparents, Muriel and James Zielske. Years after my
grandfather created a golf course that included man-made ponds made of the natural rich

blue clay in Ohio, I learned about sources/sinks and leakance. I think back to life on the

farm and golf course in Sandusky, Ohio and I am so thankful that I got to grow up with
my grandparents next door and how they have supported me in all my endeavors. They
still continue to support me and I am blessed to have such accomplished people in my
life.

Next, I would like to thank my fellow colleagues that have been supportive. It has
been a pleasure to work with Andreanne Simard, Soheil Afshari, and Chuen-Fa Ni
(David). Before I started working on this thesis, I attended mini-presentations by these
students and learned about the software package IGW. I am thankfirl for their support
throughout this work. I would also like to thank Cara Brooks and Elif Seckin in the
Mathematics department for giving me support and listening to me as I would discuss
some of the results in this work.

Lastly, I would like to thank my Industrial Mathematics advisor,

Dr. Charles MacCluer for encouraging me to explore environmental problems. And by
his encouragement, I found groundwater modeling and I was able to work with a much

respected professor and my advisor, Dr. Shu-Guang Li.

vi

TABLE OF CONTENTS

 

 

 

Introduction
1 A New Paradigm for Groundwater Modeling- _ 7
1.1 Beyond the “Black Box” .................................................................................. 7
1.2 IGW In Action ................................................................................................. 8
2 Introduction to Ground Water Flow - - -- -- ........... - - - - 9
2.1 IGW Horizontal Model .................................................................................. 10
2.1.1 Storage Coefficient ............................................................................. 11
2. l .2 Transmissivity .................................................................................... 13
2. 1.2. 1 Conductivity ....................................................................................... 13
2.1.2.2 Aquifer Thickness .............................................................................. 14
Outer Iteration or Water Table Iteration ......................................................... 15
2.1.3 Sources and Sinks ............................................................................... 17
2.1.3.1 Explicit Flux: Non- Head-dependent .................................................. 18
Injection and Pumping Wells ......................................................................... 18
Recharge ........................................................................................................ 19
2.1.3.2 Implicit Fluxes: Head-dependent Flux ............................................... 21
Rivers and Lakes ............................................................................................ 22
Drains ............................................................................................................ 26
Evapotranspiration ......................................................................................... 30
2.1.4 Boundary Conditions .......................................................................... 31
2.1.4.1 Mathematical Descriptions for IGW Boundary Conditions ................. 32
2.2 IGW Vertical Model ...................................................................................... 37
2.2.1 Boundary Conditions .......................................................................... 40
2.2.2 Sources/Sinks ..................................................................................... 43
2.2.2.1 Implicit Flux: Head-dependent ........................................................... 43
2.2.2.2 Explicit Flux: Non- Head-dependent .................................................. 43
2.2.3 Water Table Generation and Boundary Conditions ............................. 44
2.3 IGW Profile Model-2V2 D Model ................................................................... 44
2.3.1 Governing Equation ........................................................................... 45
2.3.2 Sources/Sinks ..................................................................................... 46
2.3.3 Boundary Conditions .......................................................................... 48
No-Flow and Flux .......................................................................................... 49
Constant-Head ............................................................................................... 49
Mixed or Head-dependent .............................................................................. 50
2.3.4 Advantages and Limitations ............................................................... 51
2.4 Radially Symmetric F low ............................................................................... 52
2.4.1 Mathematical Set-up ........................................................................... 52
2.4.1.1 Boundary Conditions .......................................................................... 53
2.4.1.2 Derivation of the 2D Radial Symmetric Equation ............................... 54
3 Illustrative Applications of Groundwater Theory- - - 57

 

vii

4

5

6

 

 

 

3.1 Groundwater Flow beneath a Dam ................................................................. 57
3.2 Toth Solution and Local, Intermediate, and Regional Flow Patterns. .............. 63
3.2.1 Toth’s Solution ................................................................................... 64
3.2.2 Flow Patterns with Heterogeneity-Layering ........................................ 67

3 .3 Moving-Boundary Time-Dependent Lake Dynamics ..................................... 71
3.3.1 Source/Sink: Transient Dynamics of the Lake .................................... 71
3.3.2 Moving-boundary-IGW delineation .................................................... 72
3.4 Moving-Boundary Time-Dependent Seepage ................................................. 78
3.4.1 Source/Sink: Transient Dynamics of Seepage ..................................... 79
3.4.2 Moving-boundary-IGW delineation .................................................... 79
3.5 Pumping Well ................................................................................................ 85
Introduction to Contaminant Transport 87
4.1 Governing Equation for Contaminant Transport ............................................. 87
4.2 Discussion of the Transport Equation Terms .................................................. 88
4.2.1 Advection ........................................................................................... 89
4.2.2 Diffusion ............................................................................................ 89
4.2.3 Mechanical Dispersion ....................................................................... 90
4.2.4 Macrodispersion ................................................................................. 91
4.2.5 Radioactive Decay or Biodegradation ................................................. 93
4.2.6 Sorption .............................................................................................. 93
4.2.7 Sources and Sinks ............................................................................... 94
4.2.7.1 Initial conditions ................................................................................. 95
4.2.7.2 Boundary conditions used as sources/sinks ......................................... 95
4.2.7.3 3D IGW Sources/Sinks ....................................................................... 95

4.3 Numerical Methods for the Contaminant Transport Equation ......................... 96
4.3.1 Finite Difference (First Order Upwinding scheme) ............................. 97
4.3.2 Stability .............................................................................................. 97
4.3.3 Modified Method of Characteristics (MMOC) .................................... 99
4.3.4 Particle Tracking .............................................................................. 100
4.3.5 Random Walk .................................................................................. 101
4.3.5.1 Comparison of RW and MMOC in IGW .......................................... 102
4.3.5.2 IGW Improved Methods for 3D Modeling Domains ......................... 104
4.3.6 IGW Hierarchical Modeling ............................................................. 104
Illustrative Applications of Contaminant Transport 106
5.1 Effects of decay, partitioning, and dispersion ............................................... 106
5.1.1 Advection only ................................................................................. 106

5. 1.2 Dispersion ........................................................................................ 107
5.1.3 Decay ............................................................................................... 108

5. 1.4 Partitioning ....................................................................................... 109
5.2 Particle Tracking and Wellhead Protection Areas ......................................... 112
5.3 Numerical Dispersion .................................................................................. 113
5.4 Recovery of Head Contours about a Well Using Hierarchical Modeling ....... 116
Stochastic Groundwater Modeling Using IGW 118
6.1 Deterministic Representation of Heterogeneity ............................................ 119

viii

 

 

 

 

6.1.1 Layers, Zones, Polylines (preferential channels and fiactures) .......... 120
6.2 Stochastic Approach to Modeling Heterogeneity .......................................... 126
6.2.1 Statistical Characterization ............................................................... 127
6.2.1.1 Sample Statistics: Mean, Variance, Covariance ................................ 128
6.2.2 Random Field Representation of Heterogeneity ................................ 128
6.2.2.1 Unconditional Random Field ............................................................ 129
6.2.2.2 Conditional Random Field ................................................................ 129
6.2.2.3 Example of Different Realization Models created with the Same
Statistics ......................................................................................................... 130
6.2.3 Characterizing Aquifer Heterogeneity using Data ............................. 131
6.2.3.1 Exploratory Data Analysis ................................................................ 132
6.2.3.2 Variogram analysis ........................................................................... 134
6.2.3.3 IGW covariance models ................................................................... 138
6.2.4 Generation of Different Realizations ................................................ 140
6.2.4.1 Basic Decomposition ........................................................................ 140
6.2.4.2 Turning-Bands ................................................................................. 143
6.2.4.3 Spectral Algorithm ........................................................................... 151
6.2.4.4 Sequential Gaussian Simulation ........................................................ 154
6.2.5 Recursive Statistical Post-Processing ................................................ 154
6.2.6 Monte Carlo Flow and Transport Modeling ...................................... 155
6.2.6.1 Basic Idea ......................................................................................... 155
6.2.6.2 New Information from Monte Carlo ................................................. 155
6.2.6.3 Monte Carlo Simulation Procedure ................................................... 155
6.2.7 Recursive F orrnulas for Ensemble Statistics ..................................... 156
7 Illustrative Applications of Stochastic Groundwater Theory - 158
7.1 Probability Model ........................................................................................ 158
7.1.1.1 Different realizations and static profiles ............................................ 161
7.1.1.2 Means and Variances ........................................................................ 167
7.1.1.3 Head distributions ............................................................................ 169
7.1.1.4 Concentration distributions ............................................................... 170
7.1.1.5 Log conductivity distribution ............................................................ 172
7.1.1.6 Concentration breakthrough curves ................................................... 174
8 Conclusion - 176
8.1 Research Issues ............................................................................................ 177
8.2 Practitioners Issues ....................................................................................... 178
8.3 Instructional Issues ....................................................................................... 179
9 References 181
10 APPENDICES - 184
10.1 Appendix A ................................................................................................. 185
10.1.1 Seepage beneath a Dam. ................................................................... 186
10.1.2 Toth Solution: Local Model and Regional Model. ............................ 188
10.1.3 Freeze and Witherspoon. .................................................................. 190
10.1.4 Transient Lake .................................................................................. 193

ix

10.1.5 Surface Seepage ............................................................................... 195
10.1.6 Pumping Well .................................................................................. 197
10.1.7 Well Head Protection Area ............................................................... 199
10.1.8 Random Walk .................................................................................. 202
10.1.9 Numerical Dispersion and Hierarchical Modeling ............................. 204
10.1.10 Hierarchical Modeling ...................................................................... 207
10.1.11 Advection only, decay, partitioning, and dispersion of a migrating
plume .................................................................................................... 210
10.1.12 Deterministic heterogeneity using random walk and multiple modelle 3
10.1.13 Log hydraulic conductivity represented as a variable random field. .. 215
10.1.14 Macrodispersion versus local dispersion using efi‘ective parameters .218
10.1.15 Waste pond model with individual fractures ..................................... 221
10.1.16 Probability Model ............................................................................. 224
10.2 Appendix B .................................................................................................. 228

LIST OF TABLES

Table 2.1. Table of boundary conditions and IGW representations ................................. 36
Table 3.1 Input scatterpoint values to form topography. ................................................. 73
Table 3.2 Scatterpoints used in the model, row-wise, to create the topography .............. 82
Table 10.1 Scatter points ............................................................................................. 194
Table 10.2 Scatter points for topography ..................................................................... 196

xi

LIST OF FIGURES
Figure 2.1 We observe a change of flow direction with anisotropy. ............................... 17
Figure 2.2 Pumping well located at (494.45, 361.43) with pumping rate 700 GPM. ....... 19

Figure 2.3 Recharge zone assigned constant value of 2 in/yr (artificial recharge value). (a)
Shows conceptually the recharge basin. (b) Shows the horizontal model set-up,

K=6O m/day, East and West River with constant head values of -5 m. ................... 21
Figure 2.4 Surface water body gaining water in a confined aquifer. .............................. 24
Figure 2.5 Surface water body gaining water in an unconfined aquifer. ........................ 24
Figure 2.6 Surface water body losing water in a confined aquifer. ................................. 24
Figure 2.7 Surface water body losing water in an unconfined aquifer ............................. 24
Figure 2.8 Confined variable thickness losing aquifer. ................................................... 24
Figure 2.9 Disconnected surface water body. ................................................................. 25
Figure 2.10 Disconnected surface water body. ............................................................... 25

Figure 2.11 Illustration of drain elevation placed above the water table. This figure
shows that no water is leaving the aquifer through the drain. There is no impact on
the flow system. ..................................................................................................... 27

Figure 2.12 Profile model of French drain (drain elevation = -2 m). This figure shows
conceptually the water leaving the aquifer and the lowering of the water table. ...... 28

Figure 2.13 Shows the horizontal model set-up with the regional flow pattern towards the
drain. ..................................................................................................................... 28

Figure 2.14 Profile model of wetland region delineated by using the drain feature in IGW.
We set the drain elevation equal to the surface elevation. This figure shows

conceptually the water leaving the aquifer through the wetland. ............................. 30
Figure 2.15 Shows the horizontal model set-up with the regional flow pattern towards

the wetland. ........................................................................................................... 30
Figure 2.16 Illustration of evapotranspiration zones ....................................................... 31
Figure 2.17 No-flow default boundary set for modeling domain in IGW ........................ 33

Figure 2.18 Boundary conditions by sources and sinks and d represents the width of the
aquifer (i.e. for this case 40 meters (the default aquifer thickness». ....................... 36

xii

Figure 2.19 Confined aquifer in the horizontal model corresponds with the vertical model
set-up ..................................................................................................................... 38

Figure 2.20 The vertical model set-up shows the unconfined aquifer from the horizontal
model “flipped” vertically to create a uniform panel. We then can model both
unconfined and confined aquifers in the vertical model explicitly. The vertical
dynamics will only be in the xz-plane. ................................................................... 40

Figure 2.21 Confined aquifer set-up for the vertical model showing different boundary
conditions and sources/sinks. ................................................................................. 42

Figure 2.22 Unconfined aquifer set-up for the vertical model showing different boundary
conditions and sources/sinks. ................................................................................. 42

Figure 2.23 Profile with IGW built-in bounds. This figure shows conceptually the mass
balance obtained from the cross-flow in the IGW profile model. ............................ 47

Figure 2.24 Two-dimensional horizontal model with profile model of the cross-flow in
the IGW Profile Model-2.5D Model. We observe the coupling between the IGW
horizontal model and vertical model in this figure. ................................................. 48

Figure 2.25 Boundary conditions for the profile model corresponding to Figure 2.24. .. 50
Figure 2.26 Plan view of well section with boundary conditions. ................................... 53

Figure 2.27 Depiction of the radial well problem set-up using IGW 2D model. We define
the thickness of the aquifer b as the arc length for a given radius location. We also
give sample boundary conditions ........................................................................... 56

Figure 3.1 Conceptual drawing of the triangular darn set-up without sheetpile. This
figure shows the assumed boundary conditions and sources/sinks for the model. 58

Figure 3.2 Conceptual drawing of the triangular dam set-up with a sheetpile. This figure
shows the assumed boundary conditions and sources/ sinks for the model. ............. 59

Figure 3.3 The numerical solution (flownet) to the seepage beneath the dam in IGW. We
use the IGW feature of particle tracking to view the flow lines beneath the dam. ...61

Figure 3. 4 IGW calculates 3the water balance. The constant-head reservoir zone has a
negative flux (10. 2 m 3/day), which represents the seepage leaving the zone per the
default aquifer thickness (40 m). The constant-head stream zone has a positive flux
(10. 2 m 3/day), which represents the seepage entering the zone. .............................. 62

Figure 3.5 The numerical solution (flownet) to the seepage beneath the dam with an
added sheetpile in IGW. We use the IGW feature of particle tracking to view the
flow lines beneath the dam. The sheetpile causes a longer flow path reducing the

seepage rate. .......................................................................................................... 62

xiii

Figure 3.6 IGW calculates the water balance with the added sheetpile. The constant-head
reservoir zone has a negative flux (8.6 m3/day), which represents the seepage
leaving the zone per the default aquifer thickness (40 m). The constant-head stream
zone has a positive flux (8.6 m3/day), which represents the seepage entering the
zone. ...................................................................................................................... 63

Figure 3.7 Local flow pattern in a shallow aquifer obtained fi'om Toth's solution. .......... 65

Figure 3.8 Regional flow pattern in a deep aquifer obtained fiom Toth's solution. We use
the IGW particle tracking feature to observe the regional flow lines. The deep
aquifer contains the local, intermediate, and regional flow patterns. ....................... 66

Figure 3.9 Toth's solution with extreme undulations. We see a more box-like flow
pattern with greater undulations in a shallow aquifer. The flow line at the 10,000
feet location is the best representation of the box-like flow line in this example. 67

Figure 3.10 IGW numerical solution to the isotropic system considered by Freeze and
Witherspoon. We see how a slight elevation difference causes the flow lines to
extend vertically downward. .................................................................................. 67

Figure 3.11 IGW numerical solution to the anisotropic system considered by Freeze and
Witherspoon. The flow lines are not orthogonal to the head contours. ................... 68

Figure 3.12 IGW numerical solution to the heterogeneous system considered by Freeze
and Witherspoon. The upper zone has a conductivity of 1 m/day and the lower zone
has a conductivity of 10 m/day. .............................................................................. 69

Figure 3.13 IGW numerical solution to the heterogenous system considered by Freeze
and Witherspoon. The upper zone has a conductivity of 100 m/day and the lower
zone has a conductivity of 1 m/day. ....................................................................... 69

Figure 3.14 IGW numerical solution to the heterogenous system considered by Freeze
and Witherspoon. The upper zone has a conductivity of 10 m/day and the lower
zone has a conductivity of 100 m/day. ................................................................... 70

Figure 3.15 IGW numerical solution to the heterogenous system considered by Freeze
and Witherspoon. The left zone and right zone have a high conductivity of 100
m/day, and the middle zone has a low conductivity of 10 m/day. ........................... 71

Figure 3.16 Transient lake profile view at 10 days. We have two lakes in the figure
because of the low stage. This figure corresponds to the horizontal plan view in
Figure 3.17. ........................................................................................................... 74

Figure 3.17 Transient lake plan view at 10 days. We observe two lakes using the IGW
feature of viewing the river/lake depth. .................................................................. 75

Figure 3.18 Transient lake profile view at 50 days. The stage has increased to form one
lake. This figure corresponds to the horizontal plan view in Figure 3.19. ............... 76

xiv

Figure 3.19 Transient lake plan view at 50 days. We observe one lake using the IGW
feature of viewing the river/lake depth. .................................................................. 76

Figure 3.20 Transient lake profile view at 190 days. We start to see the lake level
decreasing almost forming two lakes. This figure corresponds to the horizontal plan
view in Figure 3.21. ............................................................................................... 77

Figure 3.21 Transient lake plan view at 190 days. The stage has decreased almost
forming to separate lakes. ...................................................................................... 78

Figure 3.22 Seepage problem set-up, horizontal model of high water table causing
groundwater seepage. We use the IGW visualization feature of viewing the drain
locations. ............................................................................................................... 80

Figure 3.23 NE to SW Profile model of high water table causing groundwater seepage. 81
Figure 3.24 NW to SE Profile model of high water table causing groundwater seepage. 81

Figure 3.25 Reduced water table at 60 days later. We see that the number of seepage
areas have decreased with the reduction in the water table elevation. ..................... 83

Figure 3.26 NE to SW Profile model of lower water table causing reduced groundwater

seepage. ................................................................................................................. 84
Figure 3.27 NW to SE Profile model of lower water table showing reduced groundwater

seepage. ................................................................................................................. 84
Figure 3.28 Pumping well (radial symmetry) in a vertical model. .................................. 86

Figure 4.1 Finite difi‘erence grid that shows the location of a concentration particle that
would provide a Courant number in the range 0 s C, 51 . This causes the first order

upwinding scheme to interpolate to predict the concentration at the next time step
n+1. ....................................................................................................................... 98

Figure 4.2 Finite difference grid that shows the location of a concentration particle that
would yield a Courant number in the range Cr > 1. This causes the first order
upwinding scheme to extrapolate to predict the concentration at the next time step

n+1 ........................................................................................................................ 99
Figure 4.3 Plume migration using Random Walk in a homogeneous medium versus the

plume migrating using MMOC. ........................................................................... 103
Figure 4.4 Hierarchical model with boundary conditions. ........................................... 105
Figure 5.1 Plume migration with no dispersivity, decay, or partitioning ....................... 107

Figure 5.2 Longitudinal dispersivity = 100 m and Transversal dispersivity = 5 m. ...... 108

XV

Figure 5.3 Decay with lambda = .005 1/day at 3 years ................................................. 109

Figure 5.4 Partitioning coefficient Kd = 1 ml/g ........................................................... 110
Figure 5.5 Static profiles of different plume migration properties. ............................... 111
Figure 5.6 Wellhead protection area using backward particle tracking. ........................ 113

Figure 5.7 Numerical dispersion directly related to grid resolution. We discretize
throughout the domain using a coarse grid. In hot spots of the domain (is.
contaminant transport areas), we discretize the submodel more finely. ................. 114

Figure 5.8 Numerical dispersion with very fine grid. We use the large model to capture
the boundary conditions of the child model. We then use the submodel to model
regions of detailed interest. .................................................................................. 115

Figure 5.9 In the parent model the drawdown curves are not showing at all. By using
hierarchical modeling and forming a child model about the well, IGW is able to
visually recover the drawdown curves .................................................................. 117

Figure 6.1 Clockwise fi'om upper left corner, we have mean conductivity = 20 m/day with
anvariance 1, 2, 3, and 4. ................................................................................. 118

Figure 6.2 Heterogeneity was formed deterministically by the fingering formation and the
lenses. We use the random walk transport solver to model the plume migration. We

see that the fingering and the lenses do tend to spread the contamination. ............ 121
Figure 6.3 Waste pond set-up with boundary conditions. ............................................. 123
Figure 6.4 Hydraulic conductivity and illustration of correlation scales in the media for

the waste pond. .................................................................................................... 124
Figure 6.5 Waste pond with leakage at 150 years. ........................................................ 124

Figure 6.6 Each of the submodels have a mean conductivity value of 20 m/day, ln K
variance of 4, and a correlation scale of 30 m. We observe each of the submodels
are different due to the "seed number” or random generator. Each of the submodels
are statistically equivalent. ................................................................................... 130

Figure 6.7 Variogram used to show the difference between a stationary and nonstationary
process. The figure also shows the sill, range, and nugget. .................................. 136

Figure 6.8 Schematic of field and turning-band lines (from Mantoglou and Wilson,
1982) .................................................................................................................. 145

Figure 6.9 Projection of vector v onto a turning-bands line (from Mantoglou and Wilson,
1981) .................................................................................................................. 147

xvi

Figure 7 .1 Flow field with representative scatterpoints used for a conditional Monte Carlo

simulation. ........................................................................................................... 160
Figure 7.2 Variogram of scatterpoints .......................................................................... 161
Figure 7.3 Realization 100 plume migration at the end of 4800 days ............................ 162
Figure 7.4 Static profile from original location to Monitoring Well 1 ........................... 163
Figure 7.5 Realization 104 plume migration at the end of 4800 days ............................ 164
Figure 7.6 Static profile from original location to Monitoring Well 1 ........................... 165
Figure 7.7 Realization 104 plume migration at the end of 4800 days ............................ 165
Figure 7.8 Static profile from original location to Monitoring Well 1 ........................... 166
Figure 7.9 Realization 107 plume migration at the end of 4800 days ............................ 166
Figure 7.10 Static profile from original location to Monitoring Well 1 ......................... 167
Figure 7.11 Means and Variances after 100 realizations. ............................................. 168
Figure 7.12 Head at Monitoring Well 2, away from internal boundary. ....................... 169
Figure 7.13 Head at Monitoring Well 1, near the internal boundary. ............................ 170
Figure 7.14 Concentration in a non- Gaussian distribution at lake boundary. ............... 171
Figure 7.15 Concentration in a non-Gaussian distribution away from the lake. ............ 171
Figure 7.16 Ln K distribution at Monitoring Well 2. .................................................... 173
Figure 7.17 Ln K distribution at Monitoring Well 1 (near the lake) ............................. 174
Figure 7.18 Concentration breakthrough curve at Monitoring Well 2. .......................... 175
Figure 10.1 Triangular dam with sheetpile ................................................................... 186
Figure 10.2 Porous medium ......................................................................................... 188
Figure 10.3 Freeze and Witherspoon Anisotropic Case 2 ............................................. 190
Figure 10.4 Freeze and Witherspoon Anisotropic Case 3 ............................................. 190
Figure 10.5 Freeze and Witherspoon Anisotropic Cases 4 and 5 .................................. 191
Figure 10.6 Freeze and Witherspoon Anisotropic Case 6 ............................................. 191
Figure 10.7 Lake in cross sectional view ..................................................................... 193

Figure 10.8 Scatter points for topography and seepage ................................................ 195

Figure 10.9 Pumping well ............................................................................................ 197
Figure 10.10 Well protection area ................................................................................ 200
Figure 10.11 Random walk .......................................................................................... 202
Figure 10.12 Numerical dispersion .............................................................................. 205
Figure 10.13 Hierarchical modeling ............................................................................. 208
Figure 10.14 Advection only, decay, partitioning, and dispersion of a migrating plume211
Figure 10.15 Deterministic heterogeneity using random walk and multiple models ..... 213
Figure 10.16 Log hydraulic conductivity ..................................................................... 216
Figure 10.17 Macrodispersion versus local dispersion using effective parameters ........ 219
Figure 10.18 Waste pond model with individual fractures ........................................... 222
Figure 10.19 Probability model ................................................................................... 225
Figure 10.20 Variogram .............................................................................................. 226
Figure 10.21 Undulations for Toth examples ............................................................... 229

xviii

Introduction

For the past couple of decades, groundwater researchers have tried to bridge the gap
between theory and real-world applications. A new paradigm for groundwater modeling
has been established in the software package Interactive Ground Water (IGW). This
software package eliminates some of the inefficiencies present in the existing "old
paradigm" for groundwater modeling that is used in most commercial software programs
available. Under the new paradigm, a groundwater modeler is able to change a
parameter, and is not forced to go through all the “pre-steps” (i.e. building, interpolating
the data, etc). Also, the results are obtained using a "parallel computing" scheme, which
can be seen at the modeler’s choice of time step for the duration of the simulation. We
refer to "parallel computing" not in its usual context, but as a way of "structuring
computation -one that allows seamless data routing and dynamic integration of
groundwater flow modeling, solute transport modeling, data processing, analyses,
mapping, and visualization [Li and Liu, 2004]." This allows the groundwater modeler to
establish whether or not the set of parameter values is reasonable. Whereas, in other
software programs, the visualization is seen offline after the model has completed the
desired simulation length. “Modeling under the new paradigm continually provides and
displays results that have been intelligently processed, organized, overlaid, and displayed.
It seamlessly and dynamically merges heterogeneous geospatial data into graphical
images-integrating related data to provide a more complete view of complex
interrelationships. It provides a quick connection between modeling
concepts/assumptions and their significance/implications [Li and Liu, 2004]." We use

the software package Interactive Groundwater (IGW) as a visualization and

computational tool to bridge the two worlds of theory and real-world applications. We
also show by innovative examples the advantages of this efficient visual software
package in relation to representing the porous medium, contaminant transport, and
associated parameter values. In this volume, we use IGW 3 for 2D and 2.5D modeling.
In the next volume, we will explore the capabilities of the 3D modeling version IGW 4
[Liu, 2003]. Images in this thesis are presented in color.

To summarize, we start with the basics of subsurface flow and the underlying
algorithms used to develop the physical processes. We consider the three different
modeling domains that IGW uses for each of the following different models: horizontal
model (xy-domain), vertical model (xz-domain), and the profile model or 2.5 D model
(coupling of the Jay-domain and the xz-domain). We then solve historically famous
problems in fluid dynamics and groundwater, namely, flow beneath a triangular dam and
Toth’s solution [Toth 1962, 1963]. We address issues such as groundwater seepage and
local, intermediate, and regional flow patterns. We then extend Toth’s results by using
the results from Freeze and Witherspoon, that is, the presence of anisotropy and
heterogeneity in flow patterns [Freeze and Witherspoon, 1966].

As the computational cost has decreased with technological improvements, numerical
modeling has become indomitable. Numerical software programs such as MODFLOW
and MODFLOW-based software packages have become readily available. Many great
works have been produced in educating both researchers and practitioners about
groundwater numerical methods. But the efficiency was and is still lagging, particularly,
with the sequential solution of first modeling the flow, and then using the flow results to

obtain the transport. Moreover, defining the set-up of a real-world physical application is

done “cell by cell”. For example, a river/lake is assigned as “river cells”. The river
cannot expand or decrease without extensive external programming. Even though the
initial numerical set-up may be representative of a region at a particular time, the
transient physical dynamics are heavily restricted.

The elastic property of the medium is an ever-present dynamic that is easily
represented in IGW. The physical processes of a surface water body contracting and
expanding in an aquifer are important when modeling regions, especially in places where
the watershed basin does not have a defined discharge point (i.e. glacial trough or closed
depression). Much of the elasticity of the medium is time dependent or transient.
Examples of such processes of great importance are recharge, baseflow/streamflow
interactions, and evapotranspiration. This dynamical moving-boundary transient
behavior is captured in an expanding/contracting lake region and seepage area.

As we now have an extensive history of environmental disasters with pollution and
harmful waste practices, we introduce the theory of contaminant transport with hopes of
using numerical techniques to remedy the historical past. We first need to detect and
protect. The detection of plume migration is very complex due to the heterogeneity of
the medium, internal boundary conditions, and the chemical properties of the plume. We
present the mathematical set-up of the full contaminant transport differential equation,
namely the advection-diffusion equation. We then look at different properties of the
media and the chemical properties of the plume by studying the effects of the decay
coefficient, partitioning coefi'rcient, and longitudinal and transversal dispersivity in a well
field. We also show the results from using a "parallel computing" scheme of obtaining

the flow field and transport at desired time steps.

We also show an example using macrodispersivity with values calculated from the
Gelhar and Axness (1983) article based on spectral representation of the medium. We
use derived longitudinal and transversal macrodispersion coefficients and input these
values local dispersivity values. We also introduce a heterogeneous medium to compare
with the averaged macrodispersion values. The derived input values for the
macrodispersivity are used to create an equivalent heterogeneous medium with local
longitudinal and transversal values. Then we compare and contrast the two results.

The advection-diffusion equation is rather difficult to solve because of the advective
(sharp-front) process that is being modeled. We compare and contrast using visual
examples the random walk method and modified method of characteristics. Random
walk provides great insight into the problem, but is discrete by nature. The discrete
nature of the method allows one to transfer back and forth from particles to plume
throughout a single simulation. The restrictions (i.e. spatial dimension or many release
particles) for the different methods are the driving force behind the new conceptual
modeling tool of hierarchical modeling used in IGW.

For stochastic flow theory, we first look at the governing flow equations, but with
random components. In other words, the hydraulic conductivity is represented as a
random field. The field or stochastic process replaces the average value parameter in the
original differential equation. We show examples of heterogeneity and the need for
stochastically representing the porous medium. We use the multiple models feature of
IGW and show the effects of different variability on a migrating plume.

Next we show a waste pond example of a stochastically represented heterogenous

medium with vertical fractures and sand lenses. The mathematical representation of flow

in the fractures is derived from the Darcy-Weisbach equation in fluid mechanics. The
fractures are modeled using the polyline feature of IGW.

We give an example of the Monte Carlo method (statistical sampling). The Monte
Carlo method was developed in the late 18005 [Zhang, 2002] and then was briefly
introduced to groundwater applications during the early 1960s by Warren and Price
(1961) and then by Shvidler (1964). But it wasn’t till the work of Freeze (1975) that the
Monte Carlo method became popular in groundwater theory. The Monte Carlo method is
based on probability distributions of multiple realizations. We start by discussing the
exploratory data analysis used with random sampling from a heterogenous field. The
random sample of points is referred to as scatterpoints in IGW. The exploratory data
analysis of the scatterpoints gives the usual moment statistics (i.e. median, mode, mean,
covariance, standard deviation, skewness, etc.) and other statistical charts such as the
histogram, probability distribution function, cumulative density firnction, and h-
scatterplot. We then opt to have conditional Monte Carlo simulations. Conditional
simulations are dependent on the field data and honor the data at the given locations,
unlike unconditional simulations. The scatterpoints are then used to form the
experimental semivariogram. From the variogram we obtain the range or correlation
scale, and sill (or variance). Both are used to recover the corresponding covariance
function. The covariance structure is the controlling function that relates any two points
in the domain and is used to construct equivalent realizations in the Monte Carlo
simulation. We give a thorough description of how to generate statistically equivalent
conditional realizations using a matrix decomposition method (Cholesky decomposition),

spectral approach, and briefly introduce the turning-bands method. And we also give the

formulas (recursive and experimental) for the post processing of the statistics for all the
realizations (ensemble mean, variance, etc.)

We then return to the wellfield example used to show the different contaminant
transport effects, but this time model it using stochastic parameters. We also show the
effects of a nonstationary medium (i.e. drawdown, source/sink boundaries, etc.) with
different distributions that are associated with the logarithmic hydraulic conductivity, the
head, and the concentration. We investigate the assumption that the concentration can be
represented as a Gaussian distribution.

Finally, we summarize the results for the given examples and discuss the ultimate goal
of bridging the gap between theory and real world applications. Through classroom
applications enabled to educate students and by presenting the theory associated with the
applications, we hope to bridge the gap in research literature as well as in traditional
groundwater curricula [Li and Liu, 2004]. We also emphasize the set-up of stochastic
representation for practitioners and the advantages of introducing probability-based
results. Moreover, we give insight into problems that can be considered two-dimensional
and allow us to decrease the computational time required for a three-dimensional model.

As the ultimate goal is to achieve solutions in a timely manner, we introduce future

work being done for later versions of IGW that greatly reduces computational time.

l A New Paradigm for Groundwater Modeling

Computational groundwater modeling is based on numerical solutions to partial
differential equations. The numerical solution is obtained for the flow and the transport
equations. But, there are drastic differences in "how" this is accomplished in the different
software packages available. The major difference between the solutions is whether the
results are found sequentially or in a parallel-fashion. A sequential process refers to first
obtaining the solution to the flow equation and then obtaining the solution to the transport
equation. The sequential process, as used in the old paradigm, causes many delays in that
the modeler cannot readily observe the results of the model. The interplay between the
geology, hydrology, chemistry, and management decision-making is lost in the sequential
process. A parallel process, as used in IGW, does not refer to the usual definition of
”parallel computing", but that at each time step both the flow and the transport equations
are solved. At each time step the flow and transport can be viewed, which visually
empowers the user and distinguishes IGW from other software packages. In combination
with the parallel computing, the “new paradigm” also uses an object-oriented structure.
The combination allows for the modeler to incorporate conceptual changes throughout
the modeling process and to view the results of the numerical solution in considerably

less time.

1.1 Beyond the “Black Box”

The new paradigm allows the user to control or steer the computation. The modeler
can choose the solver, the contaminant transport method, the cell dimensions, the time-
step, tolerances for the convergence criterion, starting head value, and virtually any input

that is applied to the underlying structure. This is powerful in that the modeler has

control of the data processing. IGW invites scientists and engineers to use the different
computational “blocks” to build both a proper conceptual model and numerical model.
Allowing the user to take part in the computational choices also allows the user to
understand the processing of the data, which is usually described as the “black box”.
IGW has defaults for each of the inputs above and uses these different values to solve the
model, if not provided. In addition, IGW invites the user to manage program execution

including the integrating, overlaying, and visualizing data and results.

1.2 IGW In Action

Much of the work presented in this paper is a testimonial to both the programming
structure of IGW and the visualization. We show with some examples considered novel
to civil engineering, the “new paradigm” for groundwater modeling in action. We show
throughout this work the underlying equations in order to incorporate the theory behind
the software package. We do this for completeness, but also to invite the modeler to go
beyond the “black box” and to see the benefits of understanding the “behind the scenes”
mathematical explanations to groundwater theory. This paper is to assist in
understanding the “new paradigm”, but also to assist in the construction of real-time

models.

2 Introduction to Ground Water Flow

The most firndamental determination in groundwater movement is the flow pattern.
Even when considering contaminant transport, the results depend on the flow model. We
begin by showing the mathematical equations that the software package Interactive
Ground Water (IGW) uses to represent the physical characteristics of the porous medium
and the flow pattern. We then consider applications of such equations in IGW.

In groundwater modeling, it is customary to obtain the solution to the governing flow
equation and then use the general flow pattern to obtain the solution to the transport
equation. Although we give a sequential presentation, the software package IGW
actually uses a “parallel-computing” scheme. In other words, the flow and transport are
computed, analyzed, and visualized at each time step. We begin by discussing the
theoretical background of the governing flow equations and the IGW implementation.
2.0 Governing Equation

For saturated groundwater flow, we combine the conservation of mass and Darcy’s

equation to obtain the following three-dimensional partial differential equation

a}; a a};
SSE:c—3x—;[Kij5x7]+n xeD (2.0-1)

where S, = specific storage [U],
I; = hydraulic head [L]
x,- = spatial coordinate [L]
K1,: hydraulic conductivity tensor [L/l']

r7 = the source/sink term [l/T].

Here and in subsequent equations we use indicial notation, with summation implied

over repeated indices. In (2.0-1) the indices 1' and j range from "one" to the number of
spatial dimensions M]. We assume that the head equation applies over the Nd

dimensional spatial domain D [Li et al., 2003]. In this volume, we refer to IGW 3 for 2D
and 2.5D modeling. In the next volume, we will explore the capabilities of the 3D

modeling version IGW 4 [Liu, 2003].

2.1 IGW Horizontal Model

Many groundwater flow situations are described in terms of flow patterns and
variations over a hallo-dimensional horizontal plane. This is representative for flow along

an aquifer layer. So we next average equation (2.0—1) in the vertical direction across the

aquifer depth 210p to zbm assuming that flow is mostly horizontal within the aquifer and

obtain
2, . z,
[p S, 27"": = fiv - (KVii)+ niiz (2.1-1)
zbot zbot

Throughout this work we often refer to the x, y, and z coordinates with repeated
indices and use x1, x2, and x3, respectively. We will reference the horizontal plane with
either the x or y coordinate, and the vertical plane with the z coordinate. Since most
engineering curricula stresses the x, y, and 2 directions in problem set-ups as does the user

interface in IGW as shown in equation (2.1-1), we occasionally write the equations using

the x, y, and z notation. On the other hand, much research is acquainted to the vector

representation and repeated indices of the x1, x2, and x3 directions. For the mathematical

10

derivations and the general equations we use this notation and also introduce the nabla (or
del) for gradient notation.
Equation 2.1-1 is a nonlinear partial differential equation, which can be written in the

IGW depth-averaged form as follows

Sagfihflmgjw (2.1-2)

where S(h) = storage coefficient [-]
Tyfh) = transmissivity [L2/T]
w = source/sink term [L/T]
h(x,y) = vertically averaged head [L].
We present a discussion for each of the different terms in (2.1-2). We first discuss the
storage term. Then we will discuss the transmissivity term and also provide an example
of anisotropy. Finally, we discuss the sources and sinks.

2.1.1 Storage Coefficient

Specific storage for a saturated aquifer is defined as the volume of water that a unit
volume of aquifer releases from storage under a unit decline in hydraulic head [Freeze
and Cherry, 1979]. So, for a confined aquifer, when the hydraulic head is reduced the
aquifer is compacted (i.e. the skeleton compresses the aquifer) and the water expands due
to the decrease in pressure. This dynamic is often referred to as the elasticity of the
medium. In an unconfined aquifer, the release from storage is the actual dewatering of
the soil pores. For the same yield as a confined aquifer, an unconfined aquifer requires

less head change over less area.

11

Depending on the head, IGW uses either the specific yield (unconfined aquifer) or the

specific storage (confined aquifer), as follows:

Sy, ['1 ifh<ztop

. 2.1-3
SS (zrop " zbot )a [’I if h 2 ztop ( )

502) = {

In the modeling domain, the aquifer type may vary. So, a built-in feature of IGW,
allows one to enter values for both the specific yield and the specific storage. IGW will

then internally decide what value to pick depending on the head in the aquifer. Default
values are used for specific storage (0.00001 m'l) and/or the specific yield (0.1), if no

value is specified.
We can also use the following equation [Freeze and Cherry, 1979] to estimate the

specific storage
SS = pg(a + nefl) (2.1-4)
where p = the density of the water,
g = the acceleration of gravity,
a = the compressibility of the aquifer skeleton,
ne = the effective porosity, and
,8 = the compressibility of the water (~ 4.6 x 10'10 m2/N).

This equation represents the physical properties of both the fluid and the porous medium.
Many books such as Freeze and Cherry (1979) and the IGW Online Help file also give
values for the specific storage and specific yield. In general, specific yield values are
orders of magnitude higher than specific storage. We can specify the storage coefficient

at known points throughout the domain, and allow IGW to interpolate between the points.

12

Any spatial parameter (i.e. storage coefficient, etc.) can be interpolated from known data
at given locations, which is usefirl when configuring the model with actual data points.

2.1.2 Transmissivity

Mathematically, we represent the transmissivity term as three different expressions as

follows
Ki,- (h — 2b,, ), [L2 /T] if h < 2,0,

Ta- ('1) = Kij(zrop — 21...). [L2 m if h 2 2,0, (2.1-5)
0 (aquifer locally dry) if h < 2,”,

We observe that the transmissivity is dependent on the hydraulic conductivity and the
aquifer thickness. We first discuss the hydraulic conductivity, and then we discuss the
aquifer thickness.

2.1.2.1 Conductivity

We define the conductivity matrix K as follows

K K
K=|: ‘1 12]. (2.1-6)
K21 K22

Ifthe medium is isotropic, then the principal conductivities K11 = K22 are the only

nonzero components of the matrix. IGW defines the tensor as follows

: Ki1+K22 +Kil “K'zz

 

 

K11 2 2 008(26), (2.1-7)
K' +K' K' —K'
K22 : 112 22 _ ll 2 22 005(26), (2.1-8)
and
K12 = K21: MSIIKZG) . (2.1-9)

13

Here K {1 and K '22 are defined locally and are the principal values along the preferential

directions. The parameter 6 represents the orientation of anisotropy. We label these
parameters in the upcoming IGW example of anisotropy in different zones. It is not
uncommon to discover two dominant directions of hydraulic conductivity that are at a
certain angle to each other in field studies [Kresic, 1997].
2.1.2.2 Aquifer Thickness

Besides characterizing the hydraulic conductivity, we also discuss how IGW computes
the aquifer thickness. In IGW we input the elevation of the surface, the top of the
aquifer, and the bottom of the aquifer. Ifthe head is less than the assigned top of the

aquifer, then the aquifer is considered to be an unconfined aquifer. Then the thickness
of the aquifer is defined to be h—zbo,. Ofien, for an unconfined aquifer the top of the

aquifer is assigned the same value as the surface elevation initially. For a confined

aquifer, the assigned top of the aquifer is less than or equal to the head in the aquifer. In
this case, the thickness of the aquifer is defined to be zmp-zbot.

We observe that IGW uses three different cases that apply to this equation, namely,
the unconfined aquifer, confined aquifer, and the dry region. The thickness of the aquifer
is defined using the type of aquifer (confined or unconfined). As this equation is solved
across the domain, IGW allows the type of aquifer to vary. In other words, IGW can
model different transitional regions that can actually change over time and space. In the
case where the aquifer type is inactive (i.e. transmissivity equals zero), it is similar to
MODFLOW and the designation of "dry cells". IGW has a default "re-wetting" scheme
that is dependent on the head around the dry area. Ifthe head rises above the aquifer plus

an additional threshold, the "dry cells" are assigned a head value that is a fraction of the

14

threshold (re-wetted). IGW has a default setting, in which the modeler is able to change
if needed.
Iterative Solution and Starting Head

For an unconfined aquifer, the governing equation is nonlinear. Thus, the head must
be solved iteratively as well as the transnrissivity. IGW uses an iterative form of
equation (2.1-2).

Form = 0,1,2,...

 

m+1 m+1
S(hm)ahaz = i[flj(hm)§L—]+w (2.1-10)

ax, 6x]-

where m = 0, is the starting head 11". We will discuss more about the starting head later in
this section.
Outer Iteration or Water Table Iteration

IGW uses an outer iteration or water table iteration, which is controlled by the
tolerances for the transmissivity values. It allows the user to choose the maximum
number of iterations and the relative error. Iterative solutions need an initial value and
since the iteration is dependent on the head value, a starting head value must be input into
the model. If we know the initial starting head, then it can be used as the initial head
value. The choice of starting head is very important, for instance, if the starting head
value is chosen to be less than the bottom of the aquifer, the program will crash. IGW
provides "intelligen " default values for the starting head so the aquifer is always wet
initially. Note: If the solution converges, it does not depend on the starting head value. If

the aquifer is initially described as "dry cells", then the solution will not converge. For a

steady-state solution, the default value for the starting head is the top of the aquifer, ho =

15

2,01,. For a transient solution (the initial head values are known) the default value for the

starting head is obtained from the previous time-step.
Inner Iteration or Matrix Solver

The inner iteration is a matrix equation that IGW internally solves for the head and is
based on the finite difference scheme. IGW has many different iterative methods
available that the user may choose from such as the Jacobi method, Gauss-Seidel method,
successive over relaxation method (SOR), nonstationary SOR, ADI, conjugate gradient,
algebraic multigrid (AMG), generalized minimum residual, etc. As seen before in the
transmissivity equation (2.1-5), there is three different classifications of the aquifer type
that are dependent on the resulting head values. This allows IGW to model extremely
complex situations with multiple aquifer types in a domain. IGW can model unconfined
and confined aquifers, transition regions and partially unconfined aquifers over space and
over time by this iterative scheme. The spatial parameters (i.e. conductivity, top, bottom,
surface elevation, etc.) can be interpolated from known data at given locations, which is
usefirl when configuring the model with actual data points.
2.1.2.3 Illustration: IGW Anisotropy

We now give an example of locally defined conductivity and anisotropy as modeled in
IGW. Anisotropy in the horizontal extent arises in groundwater modeling in a number of

places. For example, anisotropy is present when considering riverbeds and fi'actured

media. For the following two zone model, we entered a value of 50 m/day for K ’1 1, a

value of 2 for the anisotropy ratio K ’1 1/K ’22, and 0 =30 degrees (see Figure 2.1). The

overall conductivity was set equal to 50 m/day. We set constant-head boundaries of

16

0 meters and 5 meters on the left and right hand sides of the model, respectively, and no-

flow boundaries at the top and bottom.

llllllll

av I+v
K22325II/h)”
’4'-

-[ y .

Kil=50mllay

#

a'

lllllll

336.0 504.0

 

Steady Flow. Time Elapsed = 0 days (0.00 years)

Figure 2.1 We observe a change of flow direction with anisotropy.

2.1.3 Sources and Sinks

IGW has very strong capabilities in modeling sources/sinks, in particular, the tracking
of inflow and outflow to the source/sink. In this section, we first present the different
sources and sinks represented as non-head-dependent fluxes and then we address the
head-dependent fluxes. The explicit fluxes or non-head-dependent fluxes, such as
pumping and recharge, are assigned flux values. The head-dependent fluxes are rivers,

lakes, wetlands, drains, surface seepage, and evapotranspiration. We consider these to be

17

implicit fluxes in that the values are calculated based on the location of the hydraulic
head.
2.1.3.1 Explicit Flux: Non- Head-dependent

We now present the sources/sinks that are assigned explicit values and are not
dependent on the head.
Injection and Pumping Wells

IGW has the capability of modeling injection/pumping wells. The physical quantities
associated with a well are the magnitude of the pumping/injecting rate Qp, and the
location of the well. IGW uses a negative rate for pumping and a positive rate for
injecting. We use the Dirac delta function 5 , which has dimensions [l/L] to represent a

well as follows

W= Qp(t)5(x-xw)5(y-yw) (2.1-11)
Note that the location of the well is at (xw, yw). So, physically speaking, the well is on
when the Dirac delta frmction is O, which corresponds to the location of the well, and is
017 at any other location. Hence, IGW represents injection and pumping wells as point
sources or sinks. IGW allows the pumping/injecting rate Q, to be a general function of
time, periodic, or random (stochastic process over time).

Well screen in [G W

For the 2D model, the results do not depend on the well screen as long as the screen is
in the aquifer (i.e. located between 2,0,, and zbm). IGW allows for the screen interval b to
change based on the modeler’s input. The horizontal model uses a fully screened well, b

= 2,0,, - zbw. In the profile model of the 2.5 D solution, the default setting is the middle

18

one-third of the aquifer. We will discuss this in more detail in Section 2.3 IGW Profile

Model-2V2 D Model.

 

 

 

 

 

 

 

worn-i) :

...........

 

 

 

......

67211
Steady Flow. Time Elapsed - 0 days 9.00 years)

 

 

 

 

 

 

 

 

 

Figure 2.2 Pumping well located at (494.45, 361.43) with pumping rate 700 GPM.

Recharge

Recharge is the rate at which water infiltrates to the aquifer; it has units [L/T]. As
mentioned before, IGW has excellent capabilities in representing sources/sinks. Amongst
the other sources/sinks, the recharge parameter shows the versatility of IGW. Recharge
can be applied over any user-defined polygon. IGW requires the user to provide the
recharge values based on the user’s discretion. Recharge can be defined as a constant
value, random function over a polygon, interpolated scatterpoint data, and time-

dependently (transient description) depending on the data available. IGW allows the

19

modeler to define the recharge as deterministic, piecewise, time-dependent, or as a
stochastic process. In addition, the user can also use the recharge to calibrate the model
and/or interpolate to known recharge scatterpoints in a specific polygon. The
mathematical representation for the recharge is
w = 6‘ (2.1-12)

We give an example of a steady-state solution for a large domain. In this case, the
domain size is 792,000 feet x 594,000 feet, which is approximately 150 miles long by
102 miles wide. The constant head boundaries on the east and west sides are set to —5
meters. The recharge rate is set to 2 in/year and the overall conductivity in the model is

set to 80 m/day. We obtain the following results (see Figure 2.3)

20

 

Figure 2.3 Recharge zone assigned constant value of 2 in/yr (artificial recharge
value). (a) Shows conceptually the recharge basin. (b) Shows the horizontal model
set-up, K=60 m/day, East and West River with constant head values of -5 111.

2.1.3.2 Implicit Fluxes: Head-dependent Flux

The head-dependent flux arises for many of the different sources/sinks. We begin by
generally defining the head-dependent flux, and then we give a more specific definition
for each case. The head-dependent flux w is the product of the leakance L and the
difference between the head ho in the source and the head h immediately below or
adjacent to the source/sink.

w= L(h0 —h). (2.1.13)

21

We define the leakance L to be
=— (2.1-14)

where K ' is the vertical conductivity (or conductance) of the interface separating the
aquifer from the source, and d ' is the thickness of the interface. When entering the
leakance in IGW, it should reflect how well the aquifer communicates or is connected
with the surface water body. Sometimes the surface water body does not have any effect
on the aquifer, and then there are cases where the surface water body is the main source.
The leakance can be constant, spatially dependent, and/or time-dependent. IGW allows
the user to enter both K ’ and d' separately, the ratio K '/ d ' , or only K ' (d' is then
calculated by IGW).

We explore different cases in which the head-dependent flux is applied. IGW models
rivers, lakes, wetlands, drains, and surface seepage by the head-dependent flux.
Rivers and Lakes

IGW considers rivers and lakes to be surface water bodies. IGW uses the head in the
aquifer h versus the head in the surface water body ha and Darcy's Law across the
sediment interface for the head-dependent sources/sinks. In the following cases, we
explicitly write the leakance as K '/ d '. IGW has the following cases for sources/sinks:
Case 1: Stationg Surface Water Bodies

If the head in the aquifer is greater than the river bottom elevation Re (see Figure 2.4-

Figure 2.8) then we obtain

235% _h) ifthe (2.1-15)

22

If the head in the aquifer is less than the river bottom elevation (see Figure 2.9 and
Figure 2.10), then we obtain

w=_K_(h0—Re) ifh<Re (2.1-16)
d'(X)
Notice that this equation is simply gravitational flow using Darcy's Law across the

medium.
In equation (2.1-15) if h > ho, then we will obtain a negative value. This represents

that the aquifer (see Figure 2.4 and Figure 2.5) is losing water to the surface water body
and thus, the surface water body is a sink. When the surface water body penetrates the
aquifer (see Figure 2.4, Figure 2.5, and Figure 2.7), then IGW uses K ' as the vertical
conductivity of the sediments in the bottom of the source/sink, and d' is the thickness of
the sediments. If the surface water body does not penetrate the aquifer (see Figure 2.6
and Figure 2.8), then the user needs to provide the harmonic average of the conductivity
for the riverbed sediments and the medium (see Figure 2.6). We give the harmonic
average as follows:

, d'
K =————— (2.1-17)
d1 (12
_+___.
K1 K2

where d ' = (11 + d 2 and the subscripts represent the riverbed sediments and the separating

medium, respectively. And in the case where d' varies (i.e. top/bottom of the aquifer

varies), IGW calculates the variable thickness (see Figure 2.8 and Figure 2.10).

23

 

Profile Model 105 Prolile Model 10!)

Unconfined Aquifer

   

Figure 2.4 Surface water body gaining Figure 2.5 Surface water body
water in a confined aquifer. gaining water in an unconfined
aquifer.

Profile Model 105 Profile Model 105

Confined Aquifer

   

Figure 2.6 Surface water body losing Figure 2.7 Surface water body losing
water in a confined aquifer. water in an unconfined aquifer.

Prolile Model 105

 

Confined Aquifer .

Figure 2.8 Confined variable thickness losing aquifer.

   

Unconflned

Profile Model 105

 

Figure 2.9 Disconnected surface water body.

Confined with variable thickness

Profile Model 112

 

Figure 2.10 Disconnected surface water body.

Case 2: "Expanding" or "Contracting" Surface Water Bodies
Ifthe surface water body in the aquifer is time-dependent and the aquifer thickness

varies, then we obtain

w(x,t) = #02033) — h(x,t)) if h 2 R, (2.1-18)
d (x,t)
and
K' .
w(x,t)='—(h0(x,t)-Re) rfh <R, (2.1-19)
d (x)

for x e no) a {x : hort) 2 zsuyaceki.

25

 

We observe that the head in the surface water body can vary in time, which creates a
moving boundary.
Special Modeling Features of I GW

IGW has a special feature called scatterpoints that are assigned different elevations of
the river bottom. IGW also has a toggle button labeled "same as the surface" that assigns
the topographical surface to be used for the river bed elevations (assuming the surface
topography varies). IGW has the ability to model this complex dynamic by using the
toggle button labeled “transient” for the river/lake. This dynamic is very difficult to
represent, in that, we do not know the head until after seeing the solution. See Chapter 3
Section 3 Transient Lake Dynamics with a Moving Boundary for an example.
Drains

Drains are used to represent many different types of sinks i.e. wetlands, surface
seepage, and French drains. Drains are used to lower the water table and are used to
represent intermittent rivers, ephemeral streams, or surface water bodies with low
amounts of water (i.e. lakes with low water levels can be included). We have the
following three cases:
Else 1: Fixed Depth

If the head in the aquifer is less than or equal to the elevation head of the drain, then
we obtain

w = 0 if h < zdmm (2.1-19)
where dey'n is the elevation of the drain. Physically, this shows that no water is leaving

through the drain if the head is lower than the drain. We show a cross-section of a

rectangular curtain drain (see Figure 2.11).

26

If the head in the aquifer is greater than the elevation head of the drain (see Figure

2.12), then we obtain
K ' ,
W = 7, (h - zdrain) If h > 2.1.0... (2.1—20)

where K '/ d ' is the leakance defined by the separating layer of sediments. For a French
drain, the leakance represents how permeable the bedding and backfill would be in the
trench.

Equation (2.1-20) is the mathematical description for water leaving the aquifer
through the drain. We show an example of this below, where the drain elevation was set

to -2 m (see Figure 2.12).

 

 

Figure 2.11 Illustration of drain elevation placed above the water table. This figure
shows that no water is leaving the aquifer through the drain. There is no impact on
the flow system.

27

 

‘7.

——.‘»~’:

DRAIN ELEV; = -2 m

 

 

“‘1‘ , 9:. ‘.- '
rnnuwmz‘fixm was. .
_ wru- ~fi%““¢fl& 3": 9‘"

 

Figure 2.12 Profile model of French drain (drain elevation = -2 m). This figure
shows conceptually the water leaving the aquifer and the lowering of the water
table.

ii
I

lil
504.0

 

Figure 2.13 Shows the horizontal model set-up with the regional flow pattern
towards the drain.

Case 3: Variable Dgpth—Segps and Wetlands

Surface seepage and wetlands occur in locations that are not necessarily well defined.
It is ofien assumed that a wetland has a small depth so that the head in the wetland is

negligible.

28

For surface seepage and wetlands, the head in the aquifer can be time-dependent. We

represent this process as follows

L

W = dr(x) (h(X,t) _ zdrain) if h > zdfflifl

, (2.1-20)
0 If h S zdrain

for x e on) .=_ {x : h(t) 2 zd,a,.,,(x)}.
The area extent is defined as the region that could potentially contain surface seepage or

wetlands. IGW automatically defines the seeps area or potentially wet area in this region.
If we are modeling surface seepage, then we assign zdmm to the same as top elevation
and the surface becomes the drain and acts similar to the French drain (see Figure 2.13).
For wetlands, we assign 24mm to be a certain water level in the wetland (see Figure 2.14).

If the wetland is gaining or losing water from the aquifer, then the head is a function of
time and cannot be assumed to be constant. The leakance of the wetland is similar to the
leakance of surface water bodies, in that we consider the thickness of the wetland
sediments and conductivity through the sediments or separating medium. In the case of a
drain, IGW relies on the user to provide the proper value for the leakance.

We observe that the head can vary in time, which creates a moving boundary. Also
the sediment thickness may vary spatially. This dynamic is very difficult to represent, in
that, we do not know the head until after seeing the solution. In the case where we do not
know the boundary for the “drain”, we define a large region (sometimes the whole
domain) as a drain with the drain elevation set equal to the surface elevation. See
Chapter 3 Section 5: Transient Seepage Dynamics with a Moving Boundary for an

example.

29

 

 

Figure 2.14 Profile model of wetland region delineated by using the drain feature in
IGW. We set the drain elevation equal to the surface elevation. This figure shows
conceptually the water leaving the aquifer through the wetland.

 

Figure 2.15 Shows the horizontal model set-up with the regional flow pattern
towards the wetland.

Evapotranspiration

(This feature is part of IGW 4 and throughout this work we use IGW 3.)

30

Evapotranspiration occurs when the water table goes through the root system of plants.

IGW models this by the following

0 Ifh<zmin
W: E fl“;- r'z - <h<z (2.1-21)
t mrn max
Zmax-zmin
Et lfh>zmax

where E, is the evapotranspiration constant, h is the head in the aquifer, and 2mm is the

bottom of the root zone, and 2m is the top of the root zone.

 

rm

 

z(mirl)

 

Figure 2.16 Illustration of evapotranspiration zones.

We can also model evaporation by using a negative recharge value. But the above
evapotranspiration model is more physically-based because the system reflects the
dependence on the head.

2.1.4 Boundary Conditions
IGW allows no-flow (or flux), prescribed-head, and head-dependent boundary

conditions (also known as first, second, and third type boundary conditions).

31

2.1.4.1 Mathematical Descriptions for IGW Boundary Conditions
The IGW default boundary condition is no-flow. We first look at the default set-up of
IGW and then describe the other boundary conditions later.
1. No-flow boundary condition.
The no-flow default boundary condition is given as

22.
an

0 (2.1-22)
where n is the normal vector to the direction of flow.

We first look at the default set-up of IGW. When the user defines a polygon as the
modeling domain, IGW assigns a default boundary condition of no-flow (see
Figure 2.17). The no-flow boundary condition is a flux condition and can be

represented mathematically as

- T95 = 0, (2.1-23)
an

where T is the transmissivity. In Figure 2.17 we observe that the no-flow boundary

and the equipotentials which are orthogonal at the boundary.

32

 

 

 

 

 

Figure 2.17 No-flow default boundary set for modeling domain in IGW.

Next, we show how different boundaries can be formed using the sources/sinks and
include the polyline feature. Even though the default boundary of the domain is no-flow,
IGW provides user-defined sources/sinks that can be imposed as boundaries (see Figure
2.18). This allows the user to define a regional flow pattern, and allows the user to
experiment with different boundary conditions in order to best fit the actual flow data.

The next type of boundary condition that we discuss is the flux boundary condition.

2. Flux boundary condition.

The general form for the flux boundary condition is given as follows

:1? = my) (2.1—24)

33

3.

where y is a constant rate of the inflow/outflow of the boundary. The no—flow
boundary condition is a special case of the flux boundary condition. This flux
condition can be implemented through a recharge zone along the boundary in
IGW. IGW also uses this boundary condition for recharge zones in the following

form

6h
—T——=c x,t -b 2.1-25
an ( ) ( )

where e is the recharge (infiltration) rate [L/T], b is the width of the zone (see
Figure 2.18).

Mixed boundary condition.

The mixed boundary condition is the linear combination of the head-dependent

and the flux boundary condition.
ah
all 73 a 7 ( )

where a , fl , and 7 are constants.

IGW uses the following form of the general equation

6h K'
_ 725-T _ 71,0“) _ h) (2.1-27)

where K ' / d' is defined as the leakance (see Section 2.1.3), b is the width of the
zone, and ho is the specified head or stage. The mixed boundary condition can be
implemented in IGW as a head-dependent flux zone along the boundary.
Prescribed-head boundary condition.

The general form assigns the head to a constant value

h = 7(x,t) (2.1-28)

34

where y is a firnction of space and time. The prescribed-head boundary condition
can be implemented in IGW along a polyline.
IGW uses the following form of the general equation

h = ho (x,t) (2.1-29)
where ho is a specified head value.

The initial constant-head boundary can be used when adding stress far from the
boundary in the modeling domain. By using scattered data to delineate the head at the
boundary before adding internal stress such as pumping, the impact can be modeled so
the initial head is correct on the boundary. As a result, the flow pattern is different
internally.

We give an example where the different boundary conditions are placed around the

perimeter of the modeling domain (see Figure 2.18).

35

6h

   
 

NO FLOW T— = 0
6y
en
T— = 10d(2 — h)
mxsn BOUNDARY common a
i h f-

 

.._,.. ......

CONSTANT -HEAD BOUNDARY CONDITION
NOIJJflhIOD AHVGMIIOE Xfl'fl

 

 

‘ =—: <—: - .
NO-FLOW BOUNDARY CONDITION

 

Figure 2.18 Boundary conditions by sources and sinks and d represents the width
of the aquifer (i.e. for this case 40 meters (the default aquifer thickness».

We give the following table of how the boundary condition was formed (see
Table 2.1).

Condition IGW

stage -
leakance K'/d' =10

 

Table 2.1. Table of boundary conditions and IGW representations.

36

2.2 IGW Vertical Model

The IGW 2D version model can be used to model certain vertical dynamics. In this
model set-up we are only interested in the results in the xz plane. This allows us to obtain
results in the vertical direction without creating a full three-dimensional model.

The governing equation for the IGW Vertical Model is equation (2.1-2), except we
are considering the x 1x2 plane to represent the xz-plane, and b represents the thickness of

the vertical slice in which we are studying. We obtain the following governing equation

for a 2D vertical model

s. .,,6_h=_a_[,,...,_ar_].w (22.1)
t 6x

where S, =specific storage[l/L],
K [j = conductivity tensor [L/T],
w = source/sink term [UT], and
b = 2,0 p - zbm, thickness of vertical slice [L].
Notice that thickness of the vertical slice is considered to be mop-21W. In order for the

2D version model to use the same solution process as the horizontal model, we have to
make the vertical model set-up match the existing horizontal model set-up. The
horizontal model set-up has the following key components (see Figure 2.19):

A. The xy-plane is the modeling domain, we refer to this as the "model space".

E. The attached vertically averaged subsurface and associated parameters.

37

 
    
 

. r )VERTICALLY
,5 AVERAGED
CONFINE l

AQUIFER Z . ,, *
x ‘ ‘ ‘ - .. " m0?)

Figure 2.19 Confined aquifer in the horizontal model corresponds with the vertical
model set-up.

In order to show that the vertical model has an equivalent set-up, we must show
equivalent components to A and B listed above. For A, we simply change the y-
coordinate to the z-coordinate. Just as the horizontal model has boundary conditions on
all sides of the model, the vertical model requires similar boundary conditions. The key
difference is the top boundary condition, which can be very complex. This boundary
condition must be added to form an equivalent domain (see Figure 2.20).

For B, we no longer have a vertically average subsurface when considering the y-
direction. So, we must provide an equivalent relation for the y-direction using the
existing horizontal aquifer parameters. The equivalent relation is that we can create a
uniform panel that will essentially "freeze" the y-direction, since we are only interested in
the xz dynamics. The uniform panel is created by using the "confined aquifer" settings
for the horizontal model, which gives us an equivalent relation to B. We create a uniform
"confined aquifer" in the y-direction by using the default IGW settings (if needed, the

user may change these). For a confined aquifer, we use the IGW default value for the

38

specific storage 0.00001 m'l. For the aquifer thickness, the IGW default setting is the top
of aquifer, zmp, is set to —10 meters, and the bottom of the aquifer, zbm, is set to —50

meters.

But, in order to make sure that the numerical program always uses the "confined
aquifer" settings for the y-direction, 2,0,, must be assigned a value less than the minimum
water table elevation. The user must check that the water table elevation (known by the
user) is greater than the top of aquifer, zwp, which has a default of -1 0m. If the water
table elevation is less than 2,0,” the program thinks it is an "unconfined aquifer" and the

uniformity in the y—direction will not exist.
For example, in Figure 2.20 we show a water table that has a minimum elevation of 15

meters. So, this implies that z(top) must be set less than 15 meters.

39

UNIFORM
PANEL

25
20
15 r:—'f‘_=.,

10

0 MODEL SPACE

 

  

CONFINED

VAQUIFER

 

 

Figure 2.20 The vertical model set-up shows the unconfined aquifer from the
horizontal model “flipped” vertically to create a uniform panel. We then can model
both unconfined and confined aquifers in the vertical model explicitly. The vertical
dynamics will only be in the xz-plane.

2.2.1 Boundary Conditions

The 2D vertical model in IGW can have several different types of boundary

conditions. The sides of the model and the bottom of the model can have all three types

40

of boundary conditions: constant-head, flux (no-flow), or mixed boundary conditions.
Refer to the horizontal model boundary conditions for example boundary conditions.

Similarly, the boundary condition(s) along the top of the model can vary. We consider
the confined aquifer and unconfined aquifer scenarios. Along the top of the model, we
can have a fixed variable boundary condition (unconfined aquifer) or we can specify
different conductivity zones and/or a flux boundary condition (confined aquifer). We
show illustrations for the confined and unconfined aquifers below in Figure 2.21 and
Figure 2.22. We refer the reader to Chapter 3: Illustrative Applications of Groundwater
Theory Sections 3.1 and 3.2.

In the unconfined case, the water table boundary condition differs fiom the boundary

conditions discussed previously. The water table boundary condition is as follows

h = Zwatertable (x) (2-2'2)
The Toth solution recreated in Chapter 3 Section 3.2 gives an example of a firnction

1.186(1 for zwatertable(x) -

41

IMPERVIOUS ZONE
(INACTIVE)

    
 
 

  

CONSTANT -HEAD ZONE
h=h0

. r
\NO—FLOW: FLUX BC

6h
K 51-20 LEAKANCE: MIXED B.

  

 

Kx’g‘“0 —K-%=£(ho—h
SOURCE/SINK an d 0

SOURCE/SINK DRAIN ah

HORIZ. WELL Ky ”67 =

 

Steady Flow, Tlme Elapsed = 0 days (one years)

Figure 2.21 Confined aquifer set-up for the vertical model showing different
boundary conditions and sources/sinks.

CONSTANT-HEAD ZONE

 

 

 

h=h0
v
UNCONFINED AQUIFER ‘ . KNOWN BC
PRESCRIBED HEAD. h = zwatenam (x)
6h
Kx ._ = 0 /-
ax ah
SOURCE/SINK K x ._ = o
SOURCE/SINK DRAIN 5"
HORIZ. WELL K .E = 0
.V 6y
Steady Flow, Time Elapsed - 0 days (0.00 years)

 

Figure 2.22 Unconfined aquifer set-up for the vertical model showing different
boundary conditions and sources/sinks.

 

' For the leakance zone, one may use a prescribed head (as shown in the unconfined case) or may use a
mixed boundary condition (as shown in the confined case). Also one can just assign a different
conductivity value for that zone.

42

2.2.2 Sources/Sinks

For the vertical model, the sources/sinks term is represented in the dynamics of the
flow under the water table. The sources/sinks are similar to the horizontal model, but
their application may differ.
2.2.2.1 Implicit Flux: Head-dependent

The head-dependent sources/sinks (i.e. river/lakes, drains, and evapotranspiration) can
be used in the vertical model as drains, for example. We show examples of the drains in
Figure 2.21 and
Figure 2.22. We also show an example of a leakance zone that can be modeled using
either a prescribed head, mixed boundary condition, or by using a different conductivity.
2.2.2.2 Explicit Flux: Non- Head-dependent

The non- head-dependent sources/sinks (i.e. recharge, pumping/injection wells) can
also be used in the vertical model.

For example, we can model horizontal wells by placing a well in the vertical section.
A real-world application of horizontal wells is for landfill analysis (see Figure 2.21 and
Figure 2.22). Mathematically, we have the following expression for the horizontal well

w = Qp6((x — xw)(z — zw)) (2.2-3)
where (xw, 2,.) is the location of the well, Qp is the magnitude of the pumping/injecting
well, and 6' is the Dirac delta function. The placement of the well needs to be located
sufficiently below the water table, so that the effect of purnping/inj ecting will not impact
the water table (unconfined aquifer) or reduce the head substantially (confined aquifer).
The well operates similarly to the vertical well presented in the horizontal model

sources/sinks.

43

2.2.3 Water Table Generation and Boundary Conditions

We use the IGW polyline feature to create the appropriate water table for the
unconfined aquifer. The modeling domain has a default no-flow boundary. In order to
change the default boundary conditions, different sources/sinks or the polyline feature
may be used to create the desired condition.

To summarize, the IGW vertical model is extremely valuable in capturing the vertical
dynamics that are not available in the horizontal model. The vertical model is limited, in
that; IGW cannot be used to predict free-surface. We just assume that we know the
water table for this case. But, we consider the free surface problem in the IGW Profile-

2‘/2 D Model.

2.3 IGW Profile Model-2V2 D Model

In order to provide results that, in general, can give some insight into a three-
dimensional model, we must be able to model a free surface as well as capture cross-flow
through the vertical model discussed above. Without forming a complete three-
dirnensional model, IGW is able to give results for what is considered 2% D. A "2% D"
model obtains information for the cross-flow through the vertical model space from the
horizontal model, which uses a vertical average. Since the vertical dynamics are not
solved for explicitly, the model results are considered to be "2% D". IGW automatically
updates the vertical model and the results are seen in the "IGW Profile Model". In other
words, the IGW Profile Model internally couples the vertical model with the horizontal
model.

The user starts by developing a horizontal model and then slices the model to obtain

the desired cross-sectional vertical model. The water table in the vertical model is

automatically updated from the horizontal model at each time step. We are capable of
viewing the results in the IGW Profile Model at each time step.

The IGW Profile Model is different from the traditional MODFLOW profile model, in
that it does not assume that “all flow occurs parallel to and in the plane of the profile
[Anderson and Woessner, 1992].”

2.3.1 Governing Equation

We have the same two-dimensional horizontal governing equation (2.1-2). The profile
model (vertical model with cross-flow) is developed using the two-dimensional vertical
governing equation and boundary conditions. However, the cross-section does not have
to be aligned with the principal axes and the profile model allows for a variable water
table, unlike the vertical model.

The governing equation for the profile model is as follows:

5—62 =1 Ki]- -b—ah— + w (2.3-1)
at 561' 661'

where S = S,b specific storage coefficient [-] if the aquifer is confined,

K ,-- = conductivity tensor [L/T],

w = is the net flux source/sink term and cross-flow [L/T],
b = thickness of vertical slice (also referred to as the user-defined thickness DB in

this section) [L],
{I = direction along polyline [L], and
4‘2 = vertical direction [L].
Notice that the storage coefficient does not have the option of being the "specific

yield". The profile model is constructed using the same approach as the vertical model,

45

that is, creating a uniform panel (see Section 2.2). The key difference is that the water
table information comes directly fiom the horizontal model. In the horizontal model, the
aquifer type is determined and the specific storage or specific yield is chosen. The
storage coefficient is then part of the water table solution and automatically updated into
the profile model.

2.3.2 Sources/Sinks

The sources/sinks used in the profile model are directly related to the sources/sinks in
the horizontal model. IGW uses the same mathematical expressions as presented in
Section 2.1 Horizontal Model: Sources/Sinks section. Additionally, in the profile model,
IGW uses the source/sink term as part of the mass balance. Since it is coupled with the
horizontal model, the source/sink term represents the inflows and outflows into the slice.
It is assumed that the inflow/outflow is uniform across the depth, since the horizontal
model uses a depth-averaged governing equation (2.1-2). IGW has a built-in feature,
which bounds the profile on each side. The user has the option of defining the thickness
of this slice DB, which is also used as b in the vertical and profile models (see Figure

2.23).

46

 

Figure 2.23 Profile with IGW built-in bounds. This figure shows conceptually the
mass balance obtained from the cross-flow in the IGW profile model.

Traditional profile modeling only allows for variation in the profile plane. But as
mentioned earlier, the IGW Profile Model has an associated thickness and allows for flow
to not be aligned with the cross-section. Referring to Figure 2.23, IGW calculates the
cross-flow of the vertical slice. Since the flow is going into and perhaps, out of, the
vertical slice, it acts as a source or sink. The expression of the cross-flow source/sink is

the following

 

 

‘,.ho ho £ -
=K-}/h DB/2 +K%DB/2 [T] (2'32)

where K i1 is the conductivity at the squares in Figure 2.24, and h-l and h+1 are obtained
5

from the horizontal model. The head ho is the head in the profile model, which also

comes from the horizontal model.

47

Profile Model 109
9.72

 
   
     
    

./w.

m F-f-flIIZ 5391:; m £2215 ELIE

 

'3...

bw-r-l
.lar/Trfi'r-r

_4l_

5‘.

    

]l.blTTITT“f~.

..t_ _
- [‘H‘J_[‘._l r .
' T

 
  
 
 

  
  

A I

All

  

I
I
h

'I
/
0

.J‘ ' ‘l . \\ . g‘ . . \_
renal-r): l :l - 2, .\. .
q _ ._ \

 

 

 

 

Figure 2.24 Two-dimensional horizontal model with profile model of the cross-flow
in the IGW Profile Model-2.5D Model. We observe the coupling between the IGW
horizontal model and vertical model in this figure.

We show two wells (singularities) in the above example. The flux into the profile
model about the well is not represented accurately by the lateral flux because the physical
dynamics of the well are considered "radial". Although the flux quantities may need
adjusting, the model still provides insight into the well dynamics.

2.3.3 Boundary Conditions
In addition to equation (2.3-1), we need boundary conditions for the profile model,

which is similar to the vertical model (see Figure 2.25). The left and right boundary

conditions are defined by IGW, but depend on the profile chosen. If the profile extends

48

to the no-flow default boundary, then the boundary is no-flow. If the right or left
boundary does not extend to the no-flow boundary, then it is considered a constant-head
boundary. The constant-head boundary is an approximation because vertical variability
cannot be determined from the depth-averaged horizontal model. The top boundary is
the water table obtained from the horizontal model if the aquifer is unconfined. If the
aquifer is confined, the top boundary is the top of the aquifer. The bottom boundary is
no-flow by default.

The boundary conditions including the water table itself are obtained from the
horizontal model. The horizontal model automatically updates the profile model with
information pertaining to the water table. This, again, is done to model a free surface and
cross-flow in a vertical model.

No—Flow and Flux

The top boundary of the aquifer can be no-flow if the aquifer is confined (see the right

side of Figure 2.25). The mathematical expression for the no-flow boundary condition is

xfl=o
6n

We also see that the left side boundary condition is the default no-flow boundary
condition.
Constant-Head

In Figure 2.25, we observe that the water table and the right side of the profile are
constant-head (prescribed head) boundaries. The prescribed-head is used in the profile
model for the top boundary condition (water table) in the unconfined portion and has the

following mathematical representation

h : zwatertable (x) '

49

The constant head boundary to the right has the following mathematical representation
h = 7

where 7 is the averaged value from the horizontal model. On the lefi hand side of Figure
2.25, we see the constant-head zone obtained from the horizontal model. When placing a
constant-head boundary condition inside the no-flow default boundary, the boundary
condition for the model is considered to be the constant-head boundary condition and
would have the same mathematical expression as displayed above.
Mixed or Head-dependent

Similar to the vertical model, the leakance from a surface water body is considered to
be a mixed boundary condition (see Figure 2.25). The mixed boundary condition has the
following mathematical expression

ah K'
—K-—=—h—h
a d’(0 )

Head-dependent OI Mixed

Wenondilo from xy-modol

  
   

peouriueisuog

Constant head from
polyline (internal boundary)

Figure 2.25 Boundary conditions for the profile model corresponding to Figure
2.24.

50

We previously discussed the boundary conditions in the horizontal model and the
vertical model. We now show the mathematical expressions for the boundary conditions
for the profile model. In Figure 2.25, we observe that the riverbed is now aligned with
both principle axes. Since it now depends on two different directions because of the
coupling of the horizontal and vertical models, we use the head-dependent flux in both
directions. For example, the bottom of the riverbed would have the following boundary
condition

6h K'

-K _=_
262 d'

(h0 — h). (2.3.3)

The sides of the riverbed would have the following boundary conditions

6h K'
-K —=-—h -h. 2.3-4

Recall that K '/ d ' is the leakance defined for the riverbed sediments.

Note: The value for the leakance is not the same for the horizontal and vertical models.
In the horizontal model the leakance value will be less than in the profile model. IGW
has a correction factor that is used in the profile file to account for this difference. It
increases the leakance value to the corresponding value in the profile model.

2.3.4 Advantages and Limitations

In a sense, we have created another dimension by considering the inflows and
outflows of the slice. But there are some limitations in this approach. Previously, we
have mentioned that the vertical extent cannot include variability because it is obtained
from the horizontal mode], which uses a depth-averaged equation. Moreover, the
boundary conditions are approximate and may be different from no-flow or constant-head

conditions. Finally, if there are singularities (i.e. wells) in the profile, then the solution

51

will be inaccurate about the singularity. The solution is inaccurate because the model
uses the depth averaged Darcy's Law, which assumes uniform flow locally. The actual
physical dynamics close to the actual pumping of the well are not uniform and this
creates a discrepancy. Even though the solution is inaccurate, it still gives insight into the
problem.

The advantage of the profile model is that it gives a relative intuition into the vertical
dynamics while retaining the horizontal dynamics. The coupled horizontal and vertical
model domains form a pseudo three-dimensional model, but with far less computation
time than a three-dimensional model. In essence, we are exploiting the two-dimensional

modeling to accommodate the timely needs of solving groundwater problems.

2.4 Radially Symmetric Flow

The IGW 2D version can be used to model certain radially symmetric flow for either
confined or unconfined with a known water table (free surface is known). We are
capable of modeling this because the 3D governing equation reduces to a 2D equation
using radial symmetry. The model is in the rz-plane and we take advantage of the radial
symmetry. We obtain the pumping effect by defining the well screen region as a recharge
zone with a negative rate.

Before presenting the different results, we first show the mathematical derivation to
justify the use of the vertical model as the domain for the radially symmetric problem.

2.4.1 Mathematical Set-up.

We start with the governing flow equation written in the following form

v -(KVh) = S, % (2.4-1)

52

where h = f (6, z, r). Additionally, we use boundary conditions to implement the
pumping.
2.4.1.1 Boundary Conditions

The flux boundary condition at the well casing is

ah Q
° 6: = #— for Zscreen _bot S Zw S zscreen _top (2°4'2)
W screen
where both rw and 2,, are within the cylindrical domain of the well (see Figure 2.26 and

Figure 2.27). Let bsmen represent the well screen thickness (see Figure 2.27).

The boundary condition at the influence radius (outer radius) is considered to be no-

flow (see Figure 2.26 and Figure 2.27) and has the following mathematical representation

 

K @- = O (2.4-3)
6r
NO-FLOW
Influence
1‘
Well Casing
Flu BC I
PLAN VIEW

Figure 2.26 Plan view of well section with boundary conditions.

53

The top boundary condition for the well is similar to the vertical and profile model. If the
aquifer is confined, then we have a no-flow condition for the top of aquifer. The
mathematical expression for this is

K -2 = 0 (2.4.4)
6n

If the aquifer is unconfined, then a prescribed-head boundary condition is used for the
fixed water table. The prescribed-head boundary condition has the following

mathematical expression

h = Z watertable (r) (2-4'5)
2.4.1.2 Derivation of the 2D Radial Symmetric Equation
Expanding equation (2.4-1) with the divergence operator we obtain the equation in
cylindrical coordinates

m-.. .61:

1 6 l 6
:El’(KVh)rl+;3§l-’<Vhla+ 62 s at

(2.4-6)

We assume that the problem is radially symmetrical, which gives K =K(r,z) and h=h(r,z).

We next average about the 6 -direction or the circumference of the radial problem and

obtain
B-r 0-r
j [lib-(Kw), ]+ m + w]d6l = j Sire (2.4-7)
0 r 6r 6z 0 6t

We set 61 and 02 equal to O and 6 - r , respectively, and then average over the entire

circumference of 1ength6 - r. So then (2.4-7) becomes

6r 6r 6r 62 62 6r
(2.4-8)
6 6h 6h

+—K 6r—=S —6r
62 ”z 62 Sat

54

Next we let T; = K i- b, where b = 6r, and obtain

grflg.gr,g.§r,g.§z.rzzg=ssg_fa (2.4-9)
We can simplify the equation further, in that, we only assume that the conductivity varies
in layers. In other words, we assume that the aquifer is isotropic in the horizontal
direction and that the principal axes are aligned with the flow direction. So, the off-

diagonal terms in the transmissivity tensor are zero, and so we obtain

grflg+érzg=gga (2.4-10)
Equation (2.4-10) is in the same form as that used in the horizontal model and vertical
model for the xz-plane, which justifies our choice of using a vertical model. There is a
slight variation in the physical set-up of this problem in that the “thickness”, b, is defined
by an arc which increases with distance r from the well (see Figure 2.27). Thus, the
transmissivity increases with r, when r is not fixed. To account for this increase in
transmissivity, we use linear regression with scatterpoints to obtain the following values
in the model. We show multiple layers in the depiction below to emphasize that the user

can model multiple scenarios. For example, the user can place two well screens in

different aquifer systems or layers.

55

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

z 1)
Well
Casing . |
No-F low
No-Flow
No-Flow
No-F low
No-F low or Prescribed Head
. No-Flow
Flux Aquifer
4r
No-F low

 

Figure 2.27 Depiction of the radial well problem set-up using IGW 2D model. We
define the thickness of the aquifer b as the arc length for a given radius location.
We also give sample boundary conditions

56

3 Illustrative Applications of Groundwater Theory

We now introduce illustrative applications that directly apply to the theory discussed
in Chapter 2 for groundwater flow. The illustrative applications used throughout this
paper are to show both how the "new paradigm" applies to the IGW simulations and also
to provide guidance in defining a groundwater model. The illustrative applications for
groundwater flow show both the innovative structure of IGW and they also validate the
software solutions by solving routine problems in groundwater (i.e. seepage beneath a

dam and Toth's solution for regional and localized flow).

3.1 Groundwater Flow beneath a Dam

We utilize the strong numerical component of IGW to solve a very familiar problem in
fluid mechanics, namely, flow beneath a dam. This problem is notorious in engineering
problems in that it is important for the following reasons [Manna et al., 2003]:

(a) Determination of rate of settlement of a saturated compressible soil layer.

(b) Calculation of seepage through the body of dams and stability of slopes.

(c) Calculation of uplift pressures under hydraulic structure and their safety against

piping.

(d) Groundwater flow towards wells and drainage of soil.
We first consider the flow beneath the dam and then we add a sheetpile. Then, we
compare the seepage rate for each case. Definition sketches are provided for each case
(see Figure 3.1 and Figure 3.2). We provide a brief history of the type of dam chosen to

validate the dimensions used in this model.

57

Brief History of Triangular Dams

The design chosen predates the 1936 arch dams, and is representative of dams in the
early 1900’s such as the Elephant Butte dam (1912-1916), Friant Dam (1939 -— 1942), the
Lower Crystal Springs Dam (1888), and the Grand Coulee Darn (1941). The design of
the triangular gravity dams is attributed to the French engineer J. Augustin Tortene de
Sazilly (1812-1852). Other triangular dams were noted earlier in Mexico, but the

designers are unknown [Polaha].

 

fl=o ax

5" 5” l“ K-8.64 m/day

 

 

— 1000 foot
@ = 0
62

 

 

 

 

 

Figure 3.1 Conceptual drawing of the triangular dam set-up without sheetpile. This
figure shows the assumed boundary conditions and sources/sinks for the model.

58

Inactive zone

Shoetpilo
Inactive 10'" Constant Head-
20 foot

  
    

C ons tam H earl =
200 f eel

    
    

 

 

 

 

 

 

 

6h 6h
_ = 500 foot ,. — = 0
6x 0 K-8.64 mvday 6x
1000 feet —
fl = 0
62

Figure 3.2 Conceptual drawing of the triangular dam set-up with a sheetpile. This
figure shows the assumed boundary conditions and sources/sinks for the model.

In this case, the dam has a structural height of 300 feet, hydraulic height of 200 feet, crest
width of 20 feet, and base width of 250 feet. The IGW default depth (y-direction) is 40
meters for the entire model base. Most crest lengths corresponding to similar dimensions
are approximately 1500 feet long, which is approximately 10 times greater than the 40
meters default. By assuming that the reservoir upstream from the darn has a sufficiently
large volume, we can assume that the change in volume is negligible. And by assuming

that the discharging river is at a predefined fixed depth, we assume that the heads h 1 and

h 2 are constant (see Figure 3.1). This implies a steady-state system. Additionally, we
assume that we have a saturated homogenous porous medium with laminar flow and that

both the water and soil is incompressible. From these assumptions, we are able to use

59

Darcy’s Law. Another common assumption associated with the vertical model is that the
horizontal conductivity is approximately an order of magnitude greater than the vertical
conductivity. We account for this anisotropy by assigning the default ratio of horizontal
to vertical conductivity to 10. We obtain the following 2-D steady-state homogenous

anisotropic equation

6 6h 6 6h 6 6h 6 6h
——K— —K— —K— —K—=O.3.1-1
ax( rraxj+ax[ rzazl+az( zraxj+az[ 2282) ( )

where the conductivities are determined from the full conductivity tensor.

For each case, we have a set of boundary conditions. Due to the complexity of the set-
up we show the boundary conditions in the respective definition sketches (see Figure 3.1
and Figure 3.2). We assume that the porous medium consists primarily of clay or low
permeability material, which gives a conductivity of 10‘7 cm/s or 8.64 m/day [Freeze and
Cherry, 1979]. Chen when installing a gravity dam the subsurface is grouted to ensure
that the medium has uniform low conductivity with no fi‘actures. This would yield a
heterogeneous medium.

The IGW numerical set-up for these problems was based on a modeling domain of
1000 x 1000 fiz. The domain was then discretized using a 100 x 100 grid. For

visualization, we utilize the feature of particle tracking via selected points within the
model domain to show the flow lines of the simulation.

We obtained the following flow solutions for the two cases (see Figure 3.3 and Figure
3.4). In order to see the actual change in volume of the seepage under the dam for
interval of time, we computed the mass balance for the porous medium (see Figure 3.4
and Figure 3.6). We found that the flux was approximately 10.2 m3/day for the first case.

With the sheetpile, we obtained a flux of 8.6 m3/day. The results show a reduction in the

60

volume of seepage for a given time period. Recall that this is the seepage rate for a
default thickness of 40 m or 131 it. As mentioned before the average crest length is
approximately 1500 ft, we use a width of 1,310 ft for this size of this dam. So we need to
multiply the seepage rates by a factor of 10 to obtain the total seepage rate. The case
without a sheetpile yields a seepage rate of 102.0 ma/day, and the added sheetpile reduces

the seepage rate to 86.0 m3/day.

lllll

400.00 (fiat)

 

Figure 3.3 The numerical solution (flownet) to the seepage beneath the dam in IGW.
We use the IGW feature of particle tracking to view the flow lines beneath the dam.

61

Water Balance

 

 

1"A'
.‘p ‘

 

 

 

 

Figure 3.4 IGW calculates the water balance. The constant-head reservoir zone has
a negative flux (10.2 m3/day), which represents the seepage leaving the zone per the
default aquifer thickness (40 m). The constant-head stream zone has a positive flux
(10.2 m3/day), which represents the seepage entering the zone.

|l|ll

4001mm“)

4000 (an)

 

Figure 3.5 The numerical solution (flownet) to the seepage beneath the dam with an
added sheetpile in IGW. We use the IGW feature of particle tracking to view the
flow lines beneath the dam. The sheetpile causes a longer flow path reducing the
seepage rate.

62

Water Balance

 

 

 

:10 CI-leag

E _

‘é o

: __ deln de Out River Well Rch. Drain Stor. Error
:3

E-1O

 

Figure 3.6 IGW calculates the water balance with the added sheetpile. The
constant-head reservoir zone has a negative flux (8.6 m3/day), which represents the
seepage leaving the zone per the default aquifer thickness (40 m). The constant-
head stream zone has a positive flux (8.6 m3/day), which represents the seepage
entering the zone.

3.2 Toth Solution and Local, Intermediate, and Regional Flow Patterns.

The hierarchy of groundwater modeling starts with a regional flow model. The
regional model encompasses local models or models at a subscale. The flow patterns
obtained for these models are solutions to the elliptical partial differential equation or
Laplace’s equation. Toth (in 1962) derived the analytical solution to the boundary-value
problem representing steady-state flow in a vertical, two-dimensional, saturated,
homogenous, isotropic flow field with the water-table represented as a constant value
[Freeze and Cherry, 1979]. To obtain, the solution of the water table at a given incline,
he projected the solution from the constant case to the inclined water table. He also
considered the case where the water table is specified as a sine curve. And then similarly
to before, he extended the solution by projecting the equipotentials to a sine curve at a
certain angle of incline. The analytical solutions for the constant value water-table and

sine curve water-table configurations are given in the book Groundwater [Freeze and

63

Cherry, 1979]. Freeze and Witherspoon then expanded the model to numerical solutions
that can handle anisotropy and nonhomogenous domains [Fetter, 2001]. We begin this
study by showing the numerical solution of pseudo-Toth’s analytical solution.

3.2.1 Toth’s Solution

Suppose that we have an unconfined aquifer with recharge zones in high topographical
regions and discharge zones in low topographical locations. Toth shows that for a
shallow groundwater basin, only local flow systems will form. But for deeper
groundwater basins, intermediate and regional flow systems are present. We show these
results using both a shallow aquifer and a deep aquifer using the sinusoidal form of the
water table that was used in the analytical derivation. We also show an example that has
a greater incline and amplitude of the sinusoidal undulations.

Mathematical set-up. For steady-state flow in a vertical, two-dimensional, saturated,

homogeneous, isotropic flow, we obtain the following Laplace's equation

2 2
Lil = o (3.2-1)
6x2 622

3‘

With the following boundary conditions and water table given as the following

g—xll(0,z) =gg(L,z) = 0

62k (x,0) = 0 (3.2-2)
h(x)|z=water table = 20 + ax + a sin bx
Where c = tan a, a = a'/ cos a, andb = b'/cosa , a' being the amplitude of the sine curve
and b' being the frequency [Freeze and Cherry, 1979]. In each case, we input the water

table equation in (2.1-3) into MATLAB to obtain a series of points to input into IGW (see

Appendix B).

64

For the shallow aquifer, we used an elevation of 1000 feet, amplitude of 50 feet, and a
frequency of 0.01. We obtain the following results from IGW (see Figure 3.7). For a

shallow aquifer we have a depth to length ratio of 1:20.

 

 
 

‘- 111111

 

1‘4: 1 =:£
Surname“); ; c

Figure 3.7 Local flow pattern in a shallow aquifer obtained from Toth's solution.

 

 

 

 

 

 

 

 

For the deep aquifer, we used an elevation of 10,000 feet and the same amplitude and
frequency as the shallow aquifer. We now show the results for a deep aquifer, which has
depth to length ratio of 1:2. We observe the regional, intermediate, and local flow

patterns that arise in the deep aquifer.

65

 

I "’r . (Ctr w" y ‘s§é£> 2 my ‘5'
900000 . . $$59vw~ff§EEww
|. :' :"”g“‘“ 59’

4. ‘i

 

 

 

Figure 3.8 Regional flow pattern in a deep aquifer obtained from Toth's solution.
We use the IGW particle tracking feature to observe the regional flow lines. The
deep aquifer contains the local, intermediate, and regional flow patterns.

The next case shows the correlation between the depth of the local flow pattern and
the amplitude of the sinusoidal undulations. We use the same sinusoidal representation
of the water table with elevation 1000 feet, amplitude 200 feet, and frequency 0.01. In the
previous case, the local contours are more rounded. In this case, we see more “box-like”
flow lines that are slightly deeper than the previous case. In general, changes in the

topography with greater undulations cause different flow patterns.

66

 

 

! E , 5 ‘ - i
l ‘ ; ; - i
r . : '

5000.0(ree't)§ " 'ronoonfi - f 3 ‘rsononfi

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 3.9 Toth's solution with extreme undulations. We see a more box-like flow
pattern with greater undulations in a shallow aquifer. The flow line at the 10,000
feet location is the best representation of the box-like flow line in this example.
3.2.2 Flow Patterns with Heterogeneity-Layering

Freeze and Witherspoon expanded Toth’s solution for cases of anisotropy and
nonhomogenous domains. We begin by showing the flow in an isotropic system (K1,: K,
= l) with no vertical exaggeration. We consider the medium to consist of clean sand/
silty sand, which gives a conductivity of 10 m/day. We refer to this case as the “base

’9

case .

3 20000 3

 

Figure 3.10 IGW numerical solution to the isotropic system considered by Freeze
and Witherspoon. We see how a slight elevation difference causes the flow lines to
extend vertically downward.

67

In order to create an anisotropic medium, we set the horizontal to vertical conductivity

ratio to 10 or 10:1. But, since the modeling domain is the xz-plane, we enter the ratio as

Kx':Ky' = 10. This implies that the general trend in the flow is horizontally.

 

 
   

    

 

IllIII/[IlllllllllllllIlI/llllll/IIIll_Ill/1111111111111]!!!ll!(at.
EI-T'F’ ~7'7 -——‘T—";-— “we" ~..- . 2

__————____’_zd.u—I—"‘ @ ..........
WWW.I‘§WEEEEEEZZEEEE§

 

 

  

l ......

100-oorm)

 

Figure 3.11 IGW numerical solution to the anisotropic system considered by Freeze
and Witherspoon. The flow lines are not orthogonal to the head contours.

In some regions, anisotropy may not be present to a large extent, for example,
studying the flow in a confined aquifer. But, heterogeneity is almost unavoidable at the
large-scale. So the next few cases, we study different nonhomogeneous configurations.

The upper layer has conductivity equal to l m/day, and the lower layer has conductivity

equal to 10 m/day.

68

 

Figure 3.12 IGW numerical solution to the heterogeneous system considered by
Freeze and Witherspoon. The upper zone has a conductivity of 1 m/day and the
lower zone has a conductivity of 10 m/day.

The next case (see Figure 3.13), the upper layer has a conductivity of 100 m/day and
the lower layer has a conductivity of 1 m/day. We observe that the flow lines are parallel

to the low conductivity zone.

 

 

 

 

..............................

 

 

. Illlll
IIIIII

 

 

 

 

 

 

Figure 3.13 IGW numerical solution to the heterogenous system considered by
Freeze and Witherspoon. The upper zone has a conductivity of 100 m/day and the
lower zone has a conductivity of 1 m/day.

69

In the next case (see Figure 3.14), we have an upper conductivity layer of 10 m/day
and a lower conductivity layer of 100 m/day. The flow lines proceed to the high

conductivity region very rapidly.

 

 

 

9;“ =-=~';.g
h‘ d‘ g; - s 3 5 a ______ : l

=‘“ ‘ u ................. ".
‘i A. ___ ___________ 1‘ . ; ; , , ‘ ...... ;
5. ~ :10000(m)~:-u;3“’§-t:2mnu;.f.?1::;;;;;;

‘ 5 ‘- - Q o w v - - - - - -

 

 

 

 

 

Figure 3.14 IGW numerical solution to the heterogenous system considered by
Freeze and Witherspoon. The upper zone has a conductivity of 10 m/day and the
lower zone has a conductivity of 100 m/day.

Finally, we extend the heterogeneity to the horizontal direction by placing a low
conductivity zone in the middle of two high conductivity zones. The low conductivity
zone has a conductivity value of 10 m/day. The high conductivity zones each have a

conductivity of 100 m/day.

70

 

 

Figure 3.15 IGW numerical solution to the heterogenous system considered by
Freeze and Witherspoon. The left zone and right zone have a high conductivity of
100 m/day, and the middle zone has a low conductivity of 10 m/day.

3.3 Moving-Boundary Time-Dependent Lake Dynamics

For the next example, we consider a lake placed in the center of a rectangular area
with no-flow boundaries. We assume that each of the no-flow boundaries is based on a
groundwater divide. As seasonal variations cause the river levels to rise and fall, we
encompass this variation into the modeling domain. In other words, the governing flow
equation is now time-dependent or transient. We add an average constant recharge to the
entire domain. We classify this problem as a 2.5-D problem which uses the governing
equation (2.3-1) presented before for the profile. Recall that a 2.5-D model couples the
horizontal model with a vertical or profile model. Since the two domains are coupled the
water table is no longer restricted to be fixed in the profile model, but is allowed to vary
based on the horizontal model.

3.3.1 Source/Sink: Transient Dynamics of the Lake
In this problem, we define the entire domain as a river/lake. The bottom of the lake is

defined to be the same as the surface elevation. By using the IGW "transient" feature, the

71

stage in the lake is determined from the previous time-step. The lake boundary is
obtained by solving for the water table at the current time-step. This also delineates the
leakance area for the aquifer, where the sediment thickness d ’ is dependent on the
location of the lake boundary. We also added constant recharge to the entire domain for
this example.

This set-up allows the actual boundary of the lake to change. Other software packages
such as MODFLOW assign “river cells” that well-define the lake a priori to running the
simulations that usually delineate the lake or river. In this case, we exploit the parallel-
computing feature of IGW, in which, it updates the model at each time step and defines
the lake boundary at each time step. This allows for a lake to expand and reduce with
seasonal fluctuations. We next explore how the boundary is delineated at each time-step.

3.3.2 Moving-boundary-IGW delineation

The capability of IGW to update at each numerical step and each time step redefining
the lake boundary is novel of this software package. The boundary of the lake is assigned
a mixed boundary condition that is solved for at each time step. Recall from Section
(2.1.4) a mixed boundary condition is mathematically given as

ah K’ .
— T5;- _ 712010 (x,t) — h(x,t)) (3.3-1)

for x e no) a {x : h0(f) 2 zswface(x)l,
where n is the normal component to the flow direction, K '/ d ' is defined as leakance, b =
(h(x,t) - zbot) is the thickness of the unconfined aquifer at the boundary of the lake, ho is

the stage which is obtained fi'om the previous time-step f , and h is the head in the aquifer

72

solved at each time step. IGW determines the boundary of the internal lake at the current

time-step based on the location of the water table.

In this example, we define the entire domain as a river/lake with default leakance of

5 /day. The hydraulic conductivity throughout the region was assumed to be silty sand,

which is typical of riverbed sediments and was given a corresponding value of 0.001

cm/sec. We added 10 in/year of recharge to the entire aquifer. We then set the bottom of

the river/lake as the surface. The surface terrain was created using scatterpoints

(conditional interpolation, the curve goes through the designated scatterpoints). In this

case, we used the following scatterpoints (see Table 3.1) given in the form (surface

elevation, top elevation, bottom elevation) from left to right to create the surface in an

overall domain size of 1000 m x 750 m. We discretized the model using a 100 x 75 grid,

which produced approximately 10 m x10 m cells.

 

 

 

 

 

 

 

 

 

Scatterpoint Surface (m) Top (m) Bottom (m)
102 35 35 -50
103 45 45 -52
104 -15 -15 -45
105 2.5 2.5 -45
106 -8 -8 -42
107 25 25 -50
108 40 4O -55

 

 

 

 

Table 3.1 Input scatterpoint values to form topography.

The transient behavior of the river was simulated using the default transient feature in

IGW. The sinusoidal variation has a period of 360 days, amplitude of 2 m, and a random

fluctuation with correlation scale of 15 days and standard deviation of 1 m.

We provide three different flow solutions corresponding to the time of 10 days, 50

days, and 190 days (see Figure 3.16 to Figure 3.21). In the first set, we see that the

aquifer is discharging to the lake. We predict that at a later time the level of the lake

 

should rise. We also predict that the lake might be unified periodically. We display the

river depths throughout the simulation in the horizontal model.

l‘mfilr .‘llrnlrl lllh’

I
r
r
r
r
r
r
r
l
l
1
r
I
I

///W

Figure 3.16 Transient lake profile view at 10 days. We have two lakes in the figure
because of the low stage. This figure corresponds to the horizontal plan view in
Figure 3.17.

 

74

 

 

 

 

l? 'i f if Wmmfluuom'fiwt it

Figure 3.17 Transient lake plan view at 10 days. We observe two lakes using the
IGW feature of viewing the river/lake depth.

In the next set of figures (see Figure 3.18 and Figure 3.19), we observe that the stage
has increased. And as predicted we see the lakes merge together. We also note that the

groundwater is discharging to the lake causing a rise in the lake level.

75

l’rnlile Muriel lllR

 

Figure 3.18 Transient lake profile view at 50 days. The stage has increased to form
one lake. This figure corresponds to the horizontal plan view in Figure 3.19.

s

.o p 4 r I - s x ~

‘2001100 4001) - . ~ ~ ~ ~

mrmmmmuodmnum

 

Figure 3.19 Transient lake plan view at 50 days. We observe one lake using the
IGW feature of viewing the river/lake depth.

76

In the set of figures below (see Figure 3.20 and Figure 3.21), we see the seasonal
change of the lake by the head in the lake decreasing, which is typical of the dry season.
The lake starts to separate into two lakes, and there is a reduction in the aquifer discharge

to the lake.

Profile Model 108

an it

 

Figure 3.20 Transient lake profile view at 190 days. We start to see the lake level
decreasing almost forming two lakes. This figure corresponds to the horizontal plan
view in Figure 3.21.

77

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

\ x \ u u . a . a » r . .
\ 1| 0 n , I t ' ' ' ' '
\ r I I I I I

I I I l l I
(I I I f I
I I II I I I
I l r , ' '
. 'Il , r r
r 4
x . .
. . \ l a r
.\ \ x \ s
\ \ \ x \ t
\\ x \ \ \
\‘l \ \ \
I r \ i \ \
I II I r l \ \ t \ t . .

I I I I I u n . . \ . . .
l . \ t . . t . .
o a 600” n \ smn . . . .

W-annoewummmm ... ,— .~

Figure 3.21 Transient lake plan view at 190 days. The stage has decreased almost
forming to separate lakes.

3.4 Moving-Boundary Time-Dependent Seepage

Often during the wet season the water table will rise and meet the surface. The
occurrence is known as groundwater seepage. Groundwater seepage happens in various
locations and is similar to a wetland. Similarly to the lake simulation, we delineated
seepage boundaries that are larger than the actual seepage areas to allow for seasonal
fluctuations. We created a hummocky topography that represents natural occurring low

spots in a region. This is very typical for a densely forested region at the toe of a slope.

78

3.4.1 Source/Sink: Transient Dynamics of Seepage

In this problem, we define the entire domain as a drain with transient recharge. The
drain elevation is defined to be the same as the surface elevation. By using the IGW
"transien " recharge feature, we are able to model the temporal fluctuations in the water
table. The seepage boundaries are obtained by solving for the water table at the current
time-step. This also delineates the leakance area for the aquifer, where the sediment
thickness d' is dependent on the location of the seepage boundaries.

This set-up allows the actual boundary of the seepage areas to change. Other software
packages such as MODFLOW assign “drain cells” that well-define the seepage or
wetland area a priori to running the simulations that usually delineate seepage areas.
Similarly to the transient lake example, we exploit the parallel-computing feature of
IGW, in which, it updates the model at each time step and defines the seepage boundaries
at each time step. This allows for the seepage areas to expand and reduce with seasonal
fluctuations. We next explore how the boundary is delineated at each time-step.

3.4.2 Moving-boundary-IGW delineation

The capability of IGW to update at each numerical step and each time step redefining
the seepage boundaries is novel of this software package. The boundaries of the seepage
areas are assigned a mixed boundary condition that is solved for at each time step. Recall
from Section (2.1.4) a mixed boundary condition is mathematically given as

612 K '
4377120.... (x.y)-h(x,y,t)) (3.4-1)

for x e no) a {x : h(t) 2 zd,a,.,,(x)},
where n is the normal component to the flow direction, K '/ d’ is defined as leakance, b =

(h(x,y,t) - 21,0.) is the thickness of the unconfined aquifer at the boundary of the seepage

79

areas, zdrain is the surface elevation which varies across the domain, and h is the head in

the aquifer solved at each time step. IGW determines the boundary of the internal

seepage areas at the current time-step based on the location of the water table.

llllllllll

8

illilllllll

I l l I I l I l
I l I I I I I I
I I I I

g
5

liltilllllllli

‘

\

|\I£1

a II: ‘ K \\
/ If..." I \“~'\ \
Ix Ms}

I
200mm):,:4oon:- LI “soon

l\\\

Olill

 

Translont Flow, Tlme Elapsed = 20 days (0.05 years)

Figure 3.22 Seepage problem set-up, horizontal model of high water table causing
groundwater seepage. We use the IGW visualization feature of viewing the drain
locations.

The domain of this problem was set-up using the default of 1000 m x 750 m. The
scatterpoints used for this problem are given in Table 3.2 below. The hydraulic

conductivity is 0.01 cm/sec, which is considered clean sand. We added two profiles,

80

which cross diagonally, northwest to southeast, and northeast to southwest. We use the
profile models throughout this section to observe the fluctuations in the water-table. For
the high water-table we observe many areas of seepage in Figure 3.22. The water-table

along the two profiles is given in Figure 3.23 and Figure 3.24.

l'll'filt‘ \lmlt-l l l-

 

Figure 3.23 NE to SW Profile model of high water table causing groundwater
seepage.

l'rullln' .".lIrIlt'l lll’r

 

W/

81

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Scatterpoint Surface (m) Top (ml Bottom (m)
122 45 45 ~45
123 10 10 -50
124 15 15 -38
125 25 25 -32
126 45 45 ~35
127 6 6 -45
128 25 25 -35
129 -2 -2 -50
130 30 30 -55
131 50 50 -47
132 2 2 -55
133 40 40 -35
134 -4 -4 -45
135 35 35 -42
136 60 60 -45
137 60 60 -37
138 60 60 -36

 

 

 

 

Table 3.2 Scatterpoints used in the model, row-wise, to create the topography.

In the first set of figures, we observed a very high water table. The rising of the water
table is very typical in the US. Northwest region that it is often referred to as the high
winter water table. The high water table yields seepage in all the corresponding low
spots of the topography. We next see the fluctuation of the water-table representative of
the dry season or summer season. In Figure 3.25 we observe that the number of seepage

areas has been reduced. Also, the diameter of the seepage area has greatly been reduced.

82

 

I\\\
t

i
\
li‘lllOIOIIIJ

\
‘
-
~
-
..
-
'
-
+
-
-
.

l E

3

5 .

 

Figure 3.25 Reduced water table at 60 days later. We see that the number of
seepage areas have decreased with the reduction in the water table elevation.

83

l’r n-‘rlv Mm‘r’l 'l/
'67J2

- 7.?

-1288 . .4

*//

 

Figure 3.26 NE to SW Profile model of lower water table causing reduced
groundwater seepage.

lQ'lulllt‘ .‘llurlrl l lfr

 

Figure 3.27 NW to SE Profile model of lower water table showing reduced
groundwater seepage.

In the horizontal view, we observe that the equipotential lines are not interrupted by
the topography, since the water table is below the hummocky topography spots. The
effect of the water-table fluctuations is critical when representing the subsurface properly

as displayed in this example.

84

3.5 Pumping Well

For this simulation we first used the Theis equation for pumping in a confined aquifer
[Freeze and Cherry, 1979] to obtain a head solution for a given set of parameters. In this
case we set the well radius (area of impact) to be 26.25 feet (8 m). The well-casing was
assumed to have a 2 inch radius (0.05 m). The thickness of the confined aquifer was
assumed to be 5.9 feet (1.8 m). The specific storage varies radially. We assign a pumping
rate of -8.676 m3/day. In the vertical model we provide a negative recharge value of
-96.4 m/day for the area inside the middle layer of the well casing. We obtain the
following IGW program results (see Figure 3.28). We have chosen to show the well with
inactive layers so that the user is aware that this feature can also represent a multi-layered
well system. In the multi-layered system, certain zones can be used to represent different

layers in the aquifer or different layers in which to add contamination, for instance.

85

 

59.57

INACTIVE

 

 

 

INACTIVE

INACTIVE

 

Figure 3.28 Pumping well (radial symmetry) in a vertical model.

The boundary conditions for the IGW model simulation were no-flow boundaries at
the left and right, except where the well screen is located. The scatterpoints were used to
represent the radially symmetry of the pump. In this case, the thickness of the cylindrical
representation of the well radius was represented as scatterpoints. We then used linear
regression. We obtained a reasonable solution in that the head contours are denser
around the well casing and lessen as they extend further. In Appendix A, we show more

details of how the actual set-up of this problem was obtained.

86

4 Introduction to Contaminant Transport

4.1 Governing Equation for Contaminant Transport

We start with the following the governing equation for 3-D transport of a single
chemical constituent in groundwater, considering advection, dispersion, fluid

sinks/sources, equilibrium-controlled sorption, and first-order irreversible rate reactions

[Zheng and Bennet, 1995]; is

ac“: 6 ac? a . . b_
R —=— D--+D*l—— —— .C+ C - C+—C 4.1-1
d 6t 6x,-[( ‘1 6le axilv' l q‘ 3 “l n, l l l

where C = the solute concentration [M/L3],
B = the aquifer thickness [L],
ne = the effective porosity H,
v, = seepage or average pore velocity in the x,- direction [L/T],
Dy- = the dispersion coefficient tensor [Lz/T],
D* = the effective diffusion coefficient [LZ/T]

xi = the Cartesian coordinate [L],

C = the concentration of solute species adsorbed to solids [M/M],

pb = the bulk density of the solids [M/L3],

ry= the decay coefficient [l/T],

qs = the volume flow rate per unit volume of the source or sink [UT], and

Cs= the solute concentration in the source or sink fluid [M/L3].

The retardation factor is defined as

87

(4.1-2)

For the IGW two-dimensional models, the transport is assumed to be primarily in the
horizontal direction and that there is no vertical dispersion (i.e. fully mixed in the z

direction). Dividing through by the retardation factor Rd in equation (4.1.1) we obtain

f39+i-3-(V,C)= —1— a (D,j +D*)a—é+ q—S—C,— Wan—’15 (4.1-3)
at Rd 6x,- Rd ax—i 6x,- Rd Rd ne

We then vertically average equation (4.1-3) as follows

L—(ngtBC)+ +R——(ner ViC)— — 'R—— 16—6 ("88(Dij + D *)-:x—C] +
“lax ‘1 xi " (4.1—4)
—C, ——[n,BC+ pbBC]
+Rd

Then using the product rule and rearranging equation (4.1-4) we obtain the following

29+__a_€.= 1 i neB(Dy°+D*laC _ qs+VmeB+Rd_1_a£ C+ qus
at Rd ax,- neBRd ax,- 6xj neBRd BRd at neBRd

 

 

(4.1-5)
The parameter descriptions are the same as in equation (4.1-l) except that q, is the

volumetric flow rate per unit area of the source sink [UT] and C is the averaged solute
concentration [M/L3]. Note that the thickness of the aquifer is time-dependent, which is

necessary for the unconfined aquifer.

4.2 Discussion of the Transport Equation Terms

Similar to the governing flow equation, we provide a mathematical summary of each

of the terms and process in the governing transport equation (4.1-4). We begin by

88

discussing the processes of the transport equation, namely, advection, dispersion, decay,

and sources/sinks.

4.2.1 Advection

Advection is the movement of the particle or contaminant dependent only on the
velocity of the flow. IGW obtains the seepage velocity from the solution to the

governing flow equations at each time-step. Recall the seepage velocity is given as

K--
v,- = ——ll-fl- for i,j = l, 2 (4.2-1)
ne 6):,-

where h = hydraulic head [L],
x,- = the directional flow length [L],
72, = effective porosity [-], and
K I]: hydraulic conductivity tensor [L/T].
4.2.2 Diffusion
Diffusion in solutions is the process whereby ionic or molecular constituents move
under the influence of their kinetic activity in the direction of their concentration gradient
[Freeze and Cherry, 1979]. In porous media, we usually refer to the effective diffusion
coefficient, D“, which accounts for the longer paths of diffusion caused by the presence
of particles in the solid matrix and because of adsorption on the solids. The effective
diffusion coefficient becomes important for sufficiently small velocities (i.e. v << 1).
The major ions in groundwater (N 3+, K+, Mg2+, Ca2+, Cl', HCO3', SO42) have

diffusion coefficients, Df, in the range for diffusion range 1 x 10'9 to 2 x 10'9 mz/s at

89

25°C. For nonreactive chemical species the diffusion coefficient, Df, has values in the
range of l x 10'10 and 1 x 10'“. The effective diffusion coefficient is obtained by

13* = wa (4.2-2)

where D* = the effective diffusion coefficient [LZ/T]
a) = the empirical coefficient [-]
Df= the diffusion coefficient [Lz/T]
The process of diffusion is also sometimes referred to as selfldr'fi‘usion, molecular
diflusion, or ionic diffusion.
4.2.3 Mechanical Dispersion
We now look at the mixing and spreading process that is dependent on the velocity.
The mixing process is referred to as dispersion. Hydrodynamic dispersion is defined by
both dispersion (mechanical dispersion) and diffusion (for low velocities). Since, we
have already discussed diffusion in Section 4.2.2, we now discuss mechanical dispersion.
Mechanical dispersion is dispersion caused entirely by the motion of the fluid. The
process by which solutes are transported by the bulk motion of the flowing groundwater
is known as advection. There are three mechanisms that cause mechanical dispersion:
varying velocities in a pore channel, the different pore sizes, and tortuosity, branching,
and fingering of pore channels. In IGW the mechanical dispersion tensor [Bear, 1972] is
defined by the following components

2

D11=(aL —aT)l:1vT+aT|v| (4.2-3)
2
D22 = (0L —a7~ )r—jl-l-GTIVI (4.2-4)

9O

1312 = 021 = (0L "057)ng (4-2-5)

where 0:1,: the longitudinal dispersivity [L],
aT= the transverse dispersivity [L],

v = the seepage velocity vector [UT], and

Iv] = ‘lvlz + v; the normalized velocity vector [L/T].

Spreading of the solute in the direction of bulk flow is known as longitudinal
dispersion. Spreading in directions perpendicular to the flow is called transverse
dispersion [Freeze and Cherry, 1979]. Dispersivity is a characteristic of the porous
medium and is determined based on the pore sizes.

If the coordinate system coincides with the principal axes, then the dispersion tensor

reduces to
D11 = aL |v| (4.2-6)
1922 = aTlvl (4.2.7)
012 = 021 = 0 (42-8)

The dispersivity values (a L and aT) represent the variance in the flow direction due to
mechanical dispersion (mixing) within the soil matrix. The longitudinal dispersivity is
usually found to be at an order of magnitude higher than the transverse dispersivity.
When considering the vertical direction, the vertical dispersivity is considered to be
somewhat random. This is due to increased variance relative to distance.

4.2.4 Macrodispersion

91

Macrodispersion is the dispersion caused at the field-scale due to heterogeneity.
Previously, we have been discussing the concept of dispersion and how to obtain the
corresponding values based on laboratory or local characterization of the porous medium.
The concept of macrodispersion is that additional dispersion can occur in the field due to
heterogeneity. The mathematical expression for the macrodispersion tensor is given as

follows (assume the coordinate system coincides with the principal axes)

1)11 = A lel (4.2.9)
022 = AT|v| (4.2-10)
D12 = 1921: 0 (4.2-11)

where A L = the longitudinal macrodispersivity [L], and
A T = the transversal macrodispersivity [L].

In IGW, we can input the macrodispersivity A for each direction as seen in the above
equations (4.2-9 to 4.2-10) in place of the local dispersivity values in equations (4.2-6 to
4.2-8), in order to characterize macrodispersion in the model. The microscopic scale for
dispersivity is also known as the local dispersivity. And the macrodispersivity is
sometimes referred to as the global dispersivity.

There are different methods on determining the values for macrodispersivity. For
example, tracer studies conducted in the field are one method to obtain values for

macrodispersivity. We can also use analytical expressions for two-dimensional flow as

follows [Gelhar, 1986]

AL = a},1/ yz (4.2-12)

92

 

2
0' aL
AT: f 2 (14.351) (4.2—13)
87 “L

where a; is the variance of the logarithmic conductivity (i.e.f=1n K), It is the statistical

value related to the size of heterogeneity, y = q /J1K1 with q as the directional flux, JL as

the longitudinal gradient, and K L as the conductivity in the longitudinal direction. The
local dispersivities, a L and a7 , are also used in determining the macrodispersivity

values.

4.2.5 Radioactive Decay or Biodegradation

In the governing equation (4.1-4), we see the decay termfllf— [neBC + pbBC], which
d

represents the loss of mass of both the dissolved phase (C) and the sorbed phase (5 ). In
IGW we are only capable of modeling first order decay. IGW allows the modeler to
input the decay coefficient. The decay coefficient is usually given in terms of the half
life

1‘12— (4.2-14)

= t1/ 2
where 11/2 is the half-life of radioactive or biodegradable materials [Zheng, 1990].
4.2.6 Sorption
The process of sorption occurs on the surfaces of solids in which an electrical charge
is imbalanced and may be satisfied by adsorbing a charged ion. Some of the common
ions in groundwater were discussed in Section 4.2.2. The transfer by adsorption or other

chemical processes of contaminant mass from the pore water to the solid part of the

93

porous medium, while flow occurs, causes the advance rate of the contaminant front to be
retarded [Freeze and Cherry, 1979].
In IGW we are capable of modeling equilibrium sorption characterized by the linear

isotherm as follows
C = KdC (4.2-15)
where C = mass of solute on the solid phase per unit mass of solid phase [M],

Kd= the partitioning coefficient for the solute with soil [L3/M], and

C = the concentration of solute in solution [M 2/L3].
The natural range of values for the partitioning coefficient is from near zero to 103 mL/ g.
For Kd values that are orders of magnitude larger than 1, the solute is essentially

immobile [Freeze and Cherry, 1979]. The partitioning coefficient can be defined in a
variety of different ways. Most often the partitioning coefficient will be chosen to model

a particular species that has a corresponding empirical isotherm. The IGW Help file
provides different Kd values for selected elements and organic compounds.

4.2.7 Sources and Sinks

The source/sink term of the governing equation, q,C,, represents the solute mass

entering the model domain through sources or leaving the model domain through sinks
[Zheng, 1990]. For the IGW 2D model, contamination cannot enter the model through
sources/sinks (i.e. rivers/lakes, drains, etc.). The concentration source/sink must be

defined by initial conditions and/or boundary conditions. In essence for the IGW 2D
models, we can remove the qu, term from equation (4.1-4) and represent the

sources/sinks of the concentration through the initial and boundary conditions.

94

4.2.7.1 Initial conditions
The initial condition is given as
C (x, y,O) = f (x, y) (4.2-l6)
This represents an instantaneous source. In IGW the user will define a concentration
polygon, which represents an instantaneous source with a given concentration.
4.2.7.2 Boundary conditions used as sources/sinks
IGW does not allow concentration to enter through the external boundary or
unsaturated zone. So at the boundary the concentration is zero
C(x,y,t) = 0 x 6 6!) (4.2-17)
where 60 represents the boundary. However, IGW does allow for the concentration to
be defined by internal boundary conditions
C(x,y,t) = f(x,y,t) x,y 6 D1 ; D (4.2-l8)
where D 1 is a designated portion (or sub-domain) contained within the domain D. In

IGW the user will define a concentration polygon, which represents a continuous source
with a given concentration.
4.2.7.3 3D IGW Sources/Sinks
The 3D version of IGW allows for the contamination to enter the model via

sources/sinks. IGW allows for the contamination to enter the system by the following
processes:

1. Infiltration: the user provides the recharge and concentration values,

2. Surface water body: the user defines the bottom elevation, leakance, head and

concentration values,

3. Injection well: the user provides the injection rate and concentration value,

95

4. Constant-head: the user provides a polygon with a prescribed-head and
concentration value,

5. Specified internal boundary condition: the user provides a polygon with a
prescribed concentration at a certain time.

6. Specified initial condition: the user provides a polygon with a prescribed

concentration as an initial condition.

4.3 Numerical Methods for the Contaminant Transport Equation

The solution of the transport partial differential equation will often not possess an
analytical solution, unless under major simplifying assumptions. This leads to a
numerical approach of solving the transport equation. Much research in computational
fluid dynamics has been completed for the notorious advection-dispersion equation. The
equation is rather difficult to solve in that the equation can be advective-dominant or
diffusion-dominant. The equation is such that for advective-dominant flow, we
encounter a sharp-front process which is very difficult to solve accurately. We discuss a
couple of the different methods to obtain insight into the different options that are
available in IGW. IGW (2D and 2.5D modeling domains) is equipped with the following
methods: the modified method of characteristics (MMOC), and the random walk method
(RW). IGW (2D and 2.5D) allows the user to choose between MMOC and RW as the
primary transport solver. We refer the reader to the MT3D manual for other schemes,
such as the TVD scheme, used in the 3D version of IGW. In order to illustrate the
methodology, we use one-dimensional uniform flow. The methodology remains the

same in more complex cases.

96

4.3.1 Finite Difference (First Order Upwinding scheme)

We show the explicit first-order upwinding scheme for solving the transport equation
because it is conceptually easy to understand. This method is not used in IGW, but we
use it to address problems with stability.

We first consider the following equation, which represents only the advection portion

of the advection-dispersion equation

LC = ...-IE (4.3-1)
6t 6x
where V = 7;:— > 0,and C = C(x, t). Additionally, we have the initial condition
d
C(x,0) = F (x) . (4.3-2)

This is used to represent the initial size of the plume. Using the explicit first-order

upwinding scheme we obtain
."+1 _ .n c." _ {1
9.79. = -2. [4.27121] . (4.3-3)

Collecting like terms and solving for C?“ , we obtain

Cpl = __"iA’ (Cir! _ (31.11 )+ C," (4.3-4)

This method is conditionally stable for the Courant (CFL) number, C, = VS: .<.. l.

 

4.3.2 Stability

Rearranging equation (4.3-4), we obtain

C," +1 = C,C,"_1 + (1 - C,)C{’ (4.3-5)

97

For C, = 1 , the truncation error is zero. In this case, we would recover the exact
solution. For C, at 1 and 0 S C, s 1, the numerical scheme averages the concentration

values at the previous time step to obtain the concentration at the next time step (see

Figure 4.1). In this case, the concentration (represented as the triangle in Figure 4.1)

would be obtained by linearly interpolating between nodes C ."_1 and C,-". For C, > 1, the

1

concentration particle is located to the left of C ."_1 in the figure below. As the first order

1
upwinding scheme only uses the nodal concentrations at i-1 and i, the concentration of a
particle to the left of the i-I node will be obtained using extrapolation (see Figure 4.2).
This causes instability and so we obtain that the Courant number must be less than or
equal to one. Ideally we would like a robust method that would not restrict the time-step
or spatial-step. Later in the chapter we look at other methods that do not require this

restriction.

 

n+1 o 0 El 0

i-2 i-l i i-l-l

 

 

 

Figure 4.1 Finite difference grid that shows the location of a concentration particle
that would provide a Courant number in the range 0 s C, 51. This causes the first

order upwinding scheme to interpolate to predict the concentration at the next time
step n+1.

98

 

i-Z i-l i i+1

 

 

 

Figure 4.2 Finite difference grid that shows the location of a concentration particle
that would yield a Courant number in the range C, > 1. This causes the first order
upwinding scheme to extrapolate to predict the concentration at the next time step
n+1.

4.3.3 Modified Method of Characteristics (MMOC)

As mentioned before the user may choose to use the modified method of
characteristics (MMOC) as the primary transport solver in IGW. This method does not
have a stability criterion like the first order upwinding scheme and is considered to be
always stable.

Mathematical Set-up for MMOC

Let 6" be the location of the concentration particle at the previous time step. If .5 n is

located between the concentration nodes Cl."_1 and C," , then we use the same equation as
the first order upwinding scheme as follows

C,” = C,C,"_1 + (1 — C, )C,-" for o s C, s 1. (4.3-6)

If 4‘" is located between two arbitrary nodes, we obtain

99

1
C,"+ =C,Cl."_(m+1) +(1—C,)C,-”__m for C, >1 (4.3-7)

where m = LC,_] the floor of the Courant number. (The floor is the integer part of a
number, for example, [2.2] = 2 or [2.6] = 2)

The location of the concentration at the previous time-step is found by using backward
particle tracking. In the case of constant velocity backward particle tracking is
represented as follows

x,H = x" — u" (x, , y" )At (4.3-8)

y.-. = y. — V. (X.sy.)At (4-3-9)
where u, and v, are the x and y velocity components at a particular location, and At is the
time-step.

In general IGW uses the fourth-order Runga-Kutta method. This technique calculates
four different velocities along the particle path (one initial, two in the middle, and one in
the final position) and then uses a weighted average of the four as the velocity in the
tracking calculation [Paulson, 2002]. For velocities at non-nodal locations, IGW uses a
bilinear interpolation scheme to determine these velocities [Paulson, 2002]. Once the
location is found, IGW uses equations (4.3-3) and (4.3-7) to determine the concentration

at the next time step.

4.3.4 Particle Tracking

We have already discussed backward particle tracking used to locate the concentration
particle at a previous time step for the MMOC method. Note: In Chapter 2 and Chapter 3
related to groundwater flow, we used forward particle tracking to show the advection only

movement through the porous medium. Particle tracking can only model transport that is

100

caused by advection. We can release a collection of particles and see plume
characteristics based on advection only. IGW uses the following mathematical
expressions for forward particle tracking as follows

xn+1 = x" +un (xmyn )At (4.3-10)

yn+1 = y" +vn (xn,yn)At (4.3-11)
We extend the particle tracking equations with a random component to obtain the
mathematical expressions for the random walk method in solving the transport equation.

The random component allows us to model dispersion.

4.3.5 Random Walk

IGW uses the following mathematical expressions for tracking the particles using RW

xn+1 = x" + u, (x, ,yn )At + ,lZDLAtficosfl — ,l2DTAtgsin19 (4.3-12)
yn+1 = y" + v, (x, ,y, )At + ,l2DLAtrf sing + ,l2DTAtg cost? (4.3-13)

where ,lZD L Ar; and 420ng are Gaussian random variables with dependency on
the time-step At , a normal random number 5 = N (0,1) or 4' = N (0,1) , and the
longitudinal and transversal dispersion coefficients D1, and Dr, respectively. We obtain

0 by using the inverse tangent as follows

19 = tan‘1 h-
an

Also, notice that if the dispersion coefficients were set to zero, we would recover the
forward particle tracking equations (4.3-10) and (4.3-11).

We show an example of random walk in a heterogeneous medium that was
constructed deterministically in Chapter 5. A disadvantage of random walk is that it is

discrete rather than continuous by nature, which causes one to release many particles

lOl

19(104) in order to replicate a concentration plume. A small amount of particles can be

used if the location of the plume is important rather than the concentration magnitude in
the plume. The amount of particles released should correspond to the concentration of
the plume. In other words, if we have a high concentration, then we should have a high
density of particles.

A couple of advantages of random walk are that it gives perfect mass balance results
and is free of numerical dispersion. It allows the user to switch back and forth between
discrete particles and a concentration plume. It also preserves the tendency for the general
movement of particles to be fi'om high to low concentrations.
4.3.5.1 Comparison of RW and MMOC in IGW

We give an illustration comparing RW and MMOC using a concentration plume over
time in a homogeneous medium (see Figure 4.3). The model was discretized using 100 x
100 cells, which gives a cell dimension of 10.1 m x 10.2 m. We use 10,000 particles to
obtain the initial concentration of 500 ppm in the rectangular region below for the
random walk method and we defined the plume to be instantaneous. We also assigned
dispersivity values of l m for the longitudinal dispersivity and 0.1 m for the transversal
dispersivity. Since, the concentration is derived at each time step from the number of
particles in a particular zone; we can actually change back and forth between a plume and

particles at different time steps for the RW method.

102

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Random Walk
1120 days
00.00 -._.
600.00
, MMOC
400.00
1120 days
00.0(m)
400.09g.) 800.00

 

 

 

 

 

Steady Flow, Time Elapsed - 4240 days (11.62 years) ‘

 

Figure 4.3 Plume migration using Random Walk in a homogeneous medium versus

the plume migrating using MMOC.

In the first snapshot of 1120 days, we see how the RW particles form a distorted
plume with jagged edges. The jagged edges and even discontinuities are from the
Gaussian random motion of the particles. Also at 1120 days we see that numerical
dispersion is already present in the MMOC plume. At 4240 days we see the
concentration as a plume for the RW results, and for the MMOC results we see more

numerical dispersion.

103

4.3.5.2 IGW Improved Methods for 3D Modeling Domains

In IGW (3D modeling), there are different methods that are to improve the accuracy of
the solution. These methods include the total variation diminishing method (TVD) and
the MT3D schemes (refer to the MT3D manual for further information): method of
characteristics (MOC), the modified method of characteristics (MMOC), and the hybrid
method of characteristics (HMOC).

4.3.6 IGW Hierarchical Modeling

Beyond the different numerical methods, IGW introduces a new innovative modeling
application to obtain efficient and accurate solutions to the transport equation. This
innovate application is called hierarchical modeling. Previously, we have seen that the
accuracy for each numerical method depends on spatial discretization. If the cell
dimensions are too big, the solution will be highly inaccurate, regardless of the stability
of the method. We now turn to hierarchical modeling to obtain more accurate and more
efficient solutions to the transport equation.

Hierarchical modeling is a process of forming submodels that are finely discretized.
For contaminant transport, very fine grids are of key importance. But, the entire domain
does not need to be finely discretized, only the regions that may contain heterogeneity or
transport. The process starts by forming a parent model that will be discretized (m x 12
cells). Then children or regions of interest within the parent model will be discretized
more finely (i x k cells). The number of cells can actually be the same in the parent and
child model, but the child domain is a smaller area which creates a finer grid. This
process can continue until the actual boundaries (i.e. inside the well) restrict the domain.

We give a depiction of the process below (Figure 4.4).

104

 

 

 

 

 

 

 

 

 

 

7 4’ ‘1
Parent Model
m x n cells
H h h h h H

@ 4 1 2 3 4 [ 3 1314

it 41

a 1-

4————-Child Model

a j x k cells

0 F e a

 

 

 

 

 

Figure 4.4 Hierarchical model with boundary conditions.

The children and grandchildren in the model assume boundary conditions from the
parent. If the boundary is not a constant-head boundary (i.e. uniform along the side of
the child domain), the boundary is linearly interpolated from the parent model or model
preceding the child in the hierarchical diagram. In this case, we let H1, H2, H3, and H4
denote the nodes from the parent model and we show the interpolation between nodes for

the boundary condition at each of the child nodes h1, h2, kg, and M.
hi = H2 + fiéjzxiz-(xi — x1) (4.3-22)

where h,- = h(x,~).

105

5 Illustrative Applications of Contaminant Transport

The innovative applications of contaminant transport are numerical dispersion,
accuracy of hierarchical modeling, and effects of decay, partitioning, and dispersivity.
IGW provides additional numerical techniques for the necessity of minimizing truncation
error in a partial differential equation directly related to the grid resolution. It also is
novel in that the effects of different chemical properties of a pollutant can be studied in a

porous media.

5.1 Effects of decay, partitioning, and dispersion

We refer to the transport equation (4.1-4) and the terms related to the decay,
partitioning (part of retardation factor), and the dispersion coefficient. We begin by
showing an example that is only dependent on the advective movement in the flow. We
then add decay, partitioning, and longitudinal and transversal dispersivity. Finally, we
compare static profiles of the different chemical properties. For each of the cases, we
discretized the domain using 100 x 77 cells of size 21.6 m x 21.7 m (so there is some
numerical dispersion in each of the plume geometries).

5.1.1 Advection only

We consider the flow between two rivers with a well field. The North and South
Rivers are constant head boundaries with stages of zero and five feet, respectively. The
Town Well located nearest to the North River has a pumping rate of 500 GPM. Rural
Well 1 and Rural Well 2 each have a pumping rate of 100 GPM, and Rural Well 3 has a
pumping rate of 200 GPM. In this case, we have an overall hydraulic conductivity of 100

fi/day and effective porosity of 0.3. No chemical properties are associated with plume.

106

‘1
332
, 355
I”
3
2.
155
; m
: 1.“
‘ 17.
151
12?
'32
N
51
.5251:
jog

 

 

000 .0(l'eet) 4000.0 6000.0
Steady Flow, Time Elapsed = 1120 days (3.07 years)

Figure 5.1 Plume migration with no dispersivity, decay, or partitioning

5.1.2 Dispersion

In the next case (see Figure 5.2), we set the longitudinal dispersivity to be 100 m and
the transversal dispersivity to be 5 m. The values for local dispersivity are usually on the
order of centimeters. We used large dispersivity values to represent the macrodispersion
of the medium. We observe the plume extending further in the direction of flow than

perpendicular to the flow.

107

 
   

 

 
        
 

 

 

 

 

 

 

 

 

4000.00
VVelll W91” 05
92
n
zoooooaeet) a
1 a
u
_ ‘:
eRural i I
Well3 . .-.
Static Profile
2000.00(feet) 4000.00 6000.00
re. ; i " *

‘ 1"5 .. a .{ ~5
' J .' , "..
1: I . 1".ng‘ '1 s x‘ - .... H“... .»,.,‘t_, r, . . g

Figure 5.2 Longitudinal dispersivity = 100 m and Transversal dispersivity = 5 m.

5.1.3 Decay

In the next case (see Figure 5.3), we set a decay rate of 0.005 l/day. In this case, we
observe the plume decaying and we expect that it will eventually be negligible. Note that
the legend significantly changes for this example. The red zone actually represents 116.7
ppm afier 3 years, when we started with 500 ppm. In three years, the plume has already
decreased substantially.

The significance of this example is that we can model different remediation techniques
such as biodegradation or microbial remediation of a plume traveling with an associated

velocity.

108

 

 

Figure 5.3 Decay with lambda = .005 l/day at 3 years

5.1.4 Partitioning

For the last case (see Figure 5.4), we set the partitioning coefficient to 1 mI/g. The
partitioning coefficient reduces the advective movement of the plume greatly. We see
that only about half of the plume moves advectively, while the other half remains in the
original position. The partitioning or adsorption to the medium is having a substantial
effect on the plume migration. We observe that the concentration of the plume still has a

maximum of 500 ppm.

109

 

2W.w(feet)

Afflu‘al i '/ i

E

;.:g:easzz§sraasee*

        
   
   

....J

' Static Profile

 

 

 

 

f—
20w.m(feet) 4000.00 6N0“)

 

Steady Flow, Time Elapsed - 1120 days (3.07 years)

Figure 5.4 Partitioning coefficient Kd = 1 ml/g

110

Static Profile

 

soo , , .- . Advecfive Only

+ Partition

—*— Decay

—-— DI sperslo n

3
0

Concentration
8

 

01!) 020 0.40 0.60 0.80 1.00
Normalized X

Figure 5.5 Static profiles of different plume migration properties.

We observe the Gaussian distribution of the plume for the first case of advection only
(see Figure 5.5). When adding decay, we still obtain a Gaussian distribution, except the
peak concentration is greatly reduced over time. The partitioning (adsorption to the
medium) of the plume restricts the peak from migrating fi'om the initial location and has a
non-Gaussian distribution. Some of the plume moves due to advection and is represented
above (see Figure 5.5). The static profile of dispersion shows that the peak is reduced
and moves further than the advective flow. This is due to the longitudinal dispersivity set

to an exaggerated 100 m. We observe that the profile related to dispersion yields a

111

Gaussian-like plume that is shifted. But, the overall distribution is classified as non-
Gaussian.
To summarize, advection only and decay yield a Gaussian distribution. While the

partitioning and dispersion yield a non-Gaussian distribution.

5.2 Particle Tracking and Wellhead Protection Areas

In Chapter 2 and Chapter 3, we used particle tracking with continuous lines to visualize
the flow lines in the examples. In the examples, we used forward particle tracking.
Another feature of particle tracking is using backward particle tracking with continuous
lines to visualize wellhead protection areas (WHPA). In this case, we simply place
particles around the well and run the simulation backwards in time. We give an
illustration of this from an example that we will pursue in more detail in Chapter 4 (see
Figure 5.6). We have two rivers that differ in head such that the gradient is fi‘om the
South River to the North River. Then we added four pumping wells with varying
pumping rates. Notice that the simulation time is approximately 15 years. Different
jurisdictions (that have adopted a WHPA) require different lengths of time for the
WHPA. Also, we can identify whether or not the plume will migrate into a well or wells
during this time. It appears that a portion of the plume will enter Rural Well 1 in 15 years

(see Figure 5.6).

112

32:2:1' Rm~a11:...21

..... Welll'f...:: ::

2000.00(feet) IZZIZ:Z"1

:1“: Run] ..........

. .... Well3 ..IIIIIIIII.
reconvene-33333400000-

In order to obtain the effect of heterogeneity on solute transport we need a method that

Chapter 6 and Chapter 7.

5.3 Numerical Dispersion

 

 

 

 

113

.......
..........
..........
.......

........

Figure 5.6 Wellhead protection area using backward particle tracking.

encompasses the dispersion. Small-scale heterogeneities cause large-scale fluctuations
and in order to capture this phenomenon, we not only need to categorize the contaminants

appropriately but also the porous medium itself. We look at stochastic subsurface flow in

In Figure 5.7, we use hierarchical modeling to show that the grid resolution is directly
proportional to numerical dispersion. The parent model uses a grid of 100 x 77 cells with

each cell 21.6 m x 21.7 m. The child model uses a much finer grid of 100 x 120 cells

    

with each cell 4.5 m x 4.5 m. The parent model shows numerical dispersion around the
edge of the particles. But, in the child model, numerical dispersion is greatly reduced by

the finer grid (see Figure 5.8).

 

  
  
  

4000.00

2000.00(feet)

  

Melly ...... ..-.._:

2000.00fleet] . 4000.00 5000.00

 

Figure 5.7 Numerical dispersion directly related to grid resolution. We discretize
throughout the domain using a coarse grid. In hot spots of the domain (i.e.
contaminant transport areas), we discretize the submodel more finely.

114

 

 

 

 

Figure 5.8 Numerical dispersion with very fine grid. We use the large model to
capture the boundary conditions of the child model. We then use the submodel to
model regions of detailed interest.

The concentration plume in Figure 5.7 does not preserve the same shape from the fine
grid to the coarse grid. Using a more finely discretized grid, we see some dispersion in
Figure 5.8, but the shape of the plume is preserved. In summary, the coarse grid parent
model shows significant error attributed to numerical dispersion. On the other hand, in
the more finely discretized child model, the plume resembles the shape of the particle
tracking solution. In general, when modeling contaminant transport, hierarchical
modeling is more accurate and more efficient than finely discretizing the entire domain

and solving.

115

5.4 Recovery of Head Contours about a Well Using Hierarchical Modeling

The accuracy of the numerical model cannot always be observed by the output of the
parent model. In this case (see Figure 5.9), we pump at 50 GPM from Rural Well 1. The
contours from the drawdown curve are not seen in the parent model. So, we form a child
model around the well and then we obtain the contours from the drawdown. We form a
grandchild model (to the parent model), and we firrther improve the accuracy of the
model.

We visually are capable of observing the details in each of the models. From this, we
can see that the hierarchical models give a more accurate solution and depiction. It is
also computationally efficient, since only the areas of interest are finely discretized. As
mentioned before, we obtain the boundary conditions from the parent model of a

particular submodel.

116

 

 

 

 

Figure 5.9 In the parent model the drawdown curves are not showing at all. By
using hierarchical modeling and forming a child model about the well, IGW is able
to visually recover the drawdown curves.

6 Stochastic Groundwater Modeling Using IGW

Stochastic groundwater theory is developed, in part, to account for the heterogeneous
nature of the porous medium. Within the porous medium, the conductivity alone can
vary by many orders of magnitude in a very small distance. Incorporating heterogeneity
into the modeling domain is important because the effects of small-scale variability cause
large-scale impacts on both flow and transport. In the following illustration (see Figure
6.1), we see that the plume migrates differently depending on magnitude of variability.
The variability increases clockwise starting in the upper left hand corner. We observe

that small changes in the variability of the medium creates noticeable differences in the

fate and transport of contaminants.

 

    

Steady Flow. Time Elapsed -'1a1o days (4.96 a...)

Figure 6.1 Clockwise from upper left corner, we have mean conductivity = 20 m/day
with In K variance 1, 2, 3, and 4.

118

In IGW we investigate heterogeneity by the following:
1. Deterministic representation of heterogeneity, and
2. Stochastic approach to modeling heterogeneity, which includes the Monte Carlo

method

6.1 Deterministic Representation of Heterogeneity

Deterministic representations of heterogeneity are user defined features such as layers,
zones, polylines, or scatterpoints, which represent different materials. The direct input of
the features is considered "deteministic" in that the user explicitly represents a
preferential channel or material. This technique is useful for modeling features that
correspond to data such as conductivity values fiom pump tests, or developing
hypothetical situations such as fiactures beneath a waste pond.

Layers or zones are created by making polygons and assigning different conductivities
to each polygon. Polylines are used to model preferential channels and conductivities
along the line. In particular, polylines can be used to model a special type of preferential
channel, namely, fi'actures in IGW.

We give a brief overview of how IGW calculates the conductivity of an individual
fracture. Recall that the mere derivation of flow in a porous medium is based on
capillary flow. Since the capillary flow is the basis of flow through porous media, we
allow certain capillaries to represent fractures.

We assume that we have a stratified system, similar to that of different layers of
conductivity in the porous medium. But, the difference is that the conductivity in the

fracture is much greater than that in a stratified system with porous media for each of the

119

layers. In this case, the expression for the conductivity along the direction of the polyline

or fracture, K], has been derived fi'om flow through parallel plates and the cubic law is

b2
Kf =k07. (6.1-1)

where 7 is the specific weight of water, b is the fracture width, and k9 is dependent on

the direction of the channel and gravitational constant.

We then substitute this conductivity into the mean flow equation and obtain

b2
§=kg——J (6.1-2)
7

where J is the flow gradient. IGW provides the option of choosing the fracture width, b.
The hydraulic conductivity of the fracture can also be specified. In IGW the conductivity
across the polyline (perpendicular to the preferential flow channel) is the same as the
conductivity matrix. Later in this chapter, we give an example of a leaking waste pond.

Another feature that allows one to deterministically represent heterogeneity is
scatterpoints. Scatterpoints are used to input data or known values. Scatterpoints also
allow the user to create a continuous variation with the IGW interpolation schemes such
as inverse distance weighting, kriging, and regression
6.1.1 Layers, Zones, Polylines (preferential channels and fractures)

IGW has various features in order to input physical geometries. In this section we
present a couple of examples of deterministic heterogeneity modeled in IGW. We first
show an example of zone-based heterogeneity and then we show an example of fractures
beneath a waste pond.

Example of Zone-based heterogeneity

120

We present the models using the special feature of multiple models in IGW. Multiple
models are models that are exact copies of one another. These copies are then changed
individually to compare different results.

We define regions with conductivity 0.1 m/day and an overall conductivity of 20
m/day. The less conductive layers help in dispersing the contaminants.

A couple of deterministic ways in which the variability of the aquifer can be created is
in a fingering formation (top) or with lenses (bottom) (see Figure 6.2). We employ the
random walk transport solver for modeling the transport in this case. Recall that the
random walk method allows the user to view the migrating contaminants as a plume or

particles (as shown below).

 

 

 

 

 

 

Steady Flow,1"lme Elapsed = 5560 days (15.23 years)

Figure 6.2 Heterogeneity was formed deterministically by the lingering formation
and the lenses. We use the random walk transport solver to model the plume
migration. We see that the lingering and the lenses do tend to spread the
contamination.

121

The top figure above shows the plume separating due to the contrast in the high
conductivity region. As the group of particles moves through the fingers, some particles
move fi'om the high conductivity zone into the low conductivity zone as observed at the
edges of the low conductivity zone. In the bottom figure we observe the plume moving
around the different lenses. Observing more carefully, we see that some of the particles
have entered into the low conductivity lenses.

Example of Polylines used to Model Heterogeneity

Clay lined waste ponds and landfills used to be very common in the 1960’s, especially
in areas with an abundance of natural clay as in northern Ohio. It was assumed to be
impermeable. But fractures can form, causing the waste or leachate to readily flow
through the medium, sometimes into aquifers. In the worst-case scenario, a waste pond
could potentially leak into an aquifer used for drinking water.

Due to the difficulty involved in modeling fractures or determining the actual location
of the fractures, the subsurface is not usually modeled with a fractured medium.
Modeling fractures is critical when considering contaminant transport to the subsurface
especially since fractures are the main contributing pathway for flow in clay. In IGW,
the polyline feature represents the fractures explicitly.

In Figure 6.4, we give a hypothetical waste pond set-up and the assumed boundary
conditions. We introduce a stochastic or random medium for the general domain of the
model and will expand on this in later sections. We assign a value of 10'6 cm/sec (glacial
till) for the mean value in the random hydraulic conductivity field, and set the effective
porosity to 0.1 for the upper unconfined aquifer. The medium also has a correlation scale

ratio of 3:1 for vertical to horizontal. The bottom layer has a mean conductivity of

122

0.0001 cm/sec with effective porosity 0.15, which is characteristic of a limestone bedrock
aquifer. The sand lenses have a conductivity of 0.1 cm/sec (clean sand) and porosity 0.3.
And we use a fracture width of 0.01 m to define each of the fiacmres.

We would also like to mention the importance of the water table in the vertical model.
In this case, if the water table was not present the concentration would rapidly flow out of
the surface (i.e. concentration travels from high concentration regions to low
concentration regions). So, we impose a natural water table that is linear between the
waste pond and the lake and is constant from the waste pond to the right no-flow
boundary. We assign a no-flow boundary at the right by assuming that a groundwater
divide exists at that location. The waste pond is modeled as a continuous concentration
zone beneath the water table. Continuous was chosen instead of instantaneous; by
assuming the waste pond is large enough that the volume will not change or that the

waste pond is being continually filled by more waste.

       
 
  
  

WATER
CONSTANT TABLE
HEAD
h,‘ 550 ft hw I 570 R

 

__-_-'_«_,_,_...:;:2 ............ . 6h
SAND . -
LENSE ‘
NO-FLOW
CONSTANT
HEAD

580 it

 

UMEST ONE AOUIFE ',

\_fl=0

l
3
_ 2000 a ‘3‘ J

 

 

Figure 6.3 Waste pond set-up with boundary conditions.

123

We show the logarithmic hydraulic conductivity field with the fi‘actures below. Notice
that the fractures have high conductivity in comparison to the surrounding medium. The
sand lenses also show a similar high conductivity value. We expect waste to travel
through the fractures and then fill the sand lenses until other fi‘actures are reached, and

then eventually enter the limestone bedrock aquifer.

 

J
100.00 feet

 

400.00 feet

Figure 6.4 Hydraulic conductivity and illustration of correlation scales in the
media for the waste pond.

 

600.

 

”32:“ x”

\
100.00(feet) & V722?”
flk

400.com) 800.00 1200.00 1600.00
~ : Steady Flow, Tlme Elapsed = 54760 days (150.03 years)

        

 

 

Figure 6.5 Waste pond with leakage at 150 years.

We notice that the waste travels primarily through the fractures and the sand lenses.

' The pond leaks more on the west portion than on the east portion due to the water table.
The contours of the no-flow boundary of the water table can be seen in between the lake
and the waste pond as they are perpendicular with the water table. The water table
changes direction at the west comer of the waste pond and this creates an equipotential
that follows the contour of the pond relatively closely (unlike the east portion). So, even
though the medium has a very low mean conductivity, the flow pattern is affected by both
the fracture (in the direction of flow) and the flow being perpendicular to the
equipotentials.

Limitations

The deterministic approach is an averaged or smoothed version of the porous medium.
When applying different zones to create heterogeneity, the zones are averaged regions.
When using scatterpoints, the data is smoothed by inverse distance weighting, kriging, or
regression. In other words, the covering is averaged.

Besides averaging, there simply is not enough data to cover point to point variations.
Much of modeling is dependent on calibrating the numerical model with the experimental
data in a least squares sense. In other words, deterministic models are common from a
practitioner’s perspective because of the ability to calibrate the model with known values
from the field. But, this in itself is misleading, in that, one presumes that the model is
exactly replicating that of the real world. When, in fact, it does not represent the effect of

the variability.

125

In order to represent the variability, we use random fields or stochastic processes. In
the next section we explain the stochastic approach to modeling heterogeneity. We also

introduce the Monte Carlo method.

6.2 Stochastic Approach to Modeling Heterogeneity

The need for a stochastic approach to modeling heterogeneity arose from the
limitations of data We cannot map heterogeneity uniquely because of the limitations of
data. The data itself is only adequate for producing the mean, variance, and correlations.
From this we can obtain information about trends and the statistical measure of
heterogeneity. Uncertainty in the model is caused by the data limitation used in
quantifying spatial variability. The spatial variability is important because it can cause
large-scale impacts. In order to capture the large-scale fluctuations caused by variability,
we need to first be able to somewhat replicate the variable medium. We consider the
small-scale variability of hydraulic properties as being random [Gelhar, 1993]. Since we
already have well-developed equations for fluids with a continuum description, the
probabilistic description must be able to update the previous equations. In order to do so,
we encounter stochastic processes or random fields. Due to the many orders of
magnitude present in the porous medium, we often use the logarthmic conductivity in
place of the usual conductivity when stochastically modeling the porous medium. In
addition, the statistical properties of In K are desirable for characterizing the medium.
We give a brief derivation of the natural dependency on In K in the 1-D groundwater
flow equation as follows.

The 1-D steady-state groundwater flow equation (infinite domain) is as follows

126

2—Kfl=0. (6.2-l)
6x 6x

Expanding the equation using the product rule we obtain

2
§_K_§_I.I_+Ka h...

-—— -— 0 . 6.2-2
6x 6x 6x2 ( )

We assume that the hydraulic conductivity is strictly greater than zero and divide through

by K as follows

2
iélifl +M = o. (6.2-3)
K 6x 6x 6x2

We then observe the dependency on In K as follows

601092.614-
6x 6x 6x2

 

0. (6.2-4)

We usually refer to f = ln(K) (log conductivity) or sometimes f = 1n (T) (log
transmissivity) in the stochastic groundwater equations. IGW allows the user to display
the conductivity random field using the logarithmic hydraulic conductivity as we will
show later in this section.

We will discuss further the usage of the In K throughout this section. We first present
some of the basic components, terminology, and theory required to conceptualize a
stochastic groundwater model. At the very core of stochastic modeling is generation of a
random field or stochastic process.

6.2.1 Statistical Characterization

Using the sample data, we can obtain the mean, variance, and covariance. We now

give the mathematical representations for each of the sample statistics.

127

6.2.1.1 Sample Statistics: Mean, Variance, Covariance
Mean

The sample mean is calculated as

)7 = % in (6.2-5)
i=1
Variance
The variance is calculated by
2 1 N -
0' = —Z(xj — X (6.2-6)
N i=1
Covariance
1 N N _ _
Cf(h)=-fi22(xi—X xj—X) (6.2.7)

where x,- —xj =|h|.

6.2.2 Random Field Representation of Heterogeneity

A stochastic process is defined as a collection of random variables. A stochastic process
usually refers to time-dependent random variables. In particular, we refer to the
logarithmic hydraulic conductivity stochastic process (In K =1) as a spatially-dependent
collection of random variables written as f (x) where it stands for (x1,x2 ,...,xd )T where
d is the dimensionality. A random field can be generated using different statistical
distributions. In this chapter, we assume that the random field is jointly Gaussian

N ( ,u, 0'2 ) , which is categorized uniquely by its mean and variance or first and second

moments. We also assume that the random field is stationary. Stationary implies that the

mean (or expectation), variance, and covariance are not spatially dependent. If both

128

these criterion are met; that is, the process is jointly Gaussian and is stationary, then the
process is considered to be ergodic. If the random field is ergodic, then the spatial
statistics from each realization are equivalent to the ensemble statistics (or statistics from
the collection of realizations). There are two types of random fields that we further
investigate; that is, unconditional random fields and conditional random fields.
6.2.2.1 Unconditional Random Field

An unconditional random field uses the general statistics (i.e. mean, variance,
covariance), but does not actually honor the data values at the data point locations. We
first calculate the mean and variance and then show an arbitrary interpolation (or single
realization) between the data without honoring the data at each location. Then we use an
interpolating scheme such as kriging or inverse distance weighting to average between
the data points. We also account for the correlation scale about the data. We discuss this
further in Section 6.2.4 Generation of Different Realizations.
6.2.2.2 Conditional Random Field
Conditional simulation uses both the statistics of the data and the actual data points.
Hence, the statistics are “conditioned” on the data. The data is honored at each location
for a single realization. So when using the interpolating scheme, the realization will
actually match the data values at each data location. We also incorporate the correlation
scale into the data fitting scheme. Regions where data is lacking will be determined
randomly and will have the maximum variance. At the data point, we expect the variance

to be zero. We discuss this further in Section 6.2.4 Generation of Different Realizations.

129

We later present three methods of generating realizations both unconditionally and
conditionally, namely, matrix decomposition, spectral methods, and the turning bands
method (see Section 6.2.4).
6.2.2.3 Example of Different Realization Models created with the Same Statistics

We next show an example of different realizations created for the logarithmic
conductivity field that are statistically equivalent. We created these in IGW using the
"random" selection for conductivity and assigning a mean value of 20 m/day, variance of

4, and correlation scale of 30 meters. We observe that each realization shows a different

heterogeneous medium.

 

 

Figure 6.6 Each of the submodels have a mean conductivity value of 20 m/day, In K
variance of 4, and a correlation scale of 30 m. We observe each of the submodels
are different due to the "seed number" or random generator. Each of the
submodels are statistically equivalent.

130

6.2.3 Characterizing Aquifer Heterogeneity using Data

For a given set of data (i.e. random sampling over a polygon, field data, etc.), IGW
allows the modeler to obtain the PDF, CDF, histogram, and h-scatterplot for the data.
The h-scatterplot is useful, in that, the modeler is able to see how many correlation pairs
exist for the data. If one or two pairs exist, then the analysis will not be accurate. The
separation distance or lag can be adjusted to obtain the maximum number of lags.

The modeler is also able to view the PDF and CDF, which is a guideline for
determining whether or not the process is Gaussian. If the process is not Gaussian, the
mean and variance may not be good measure to characterize the data. Note that the log
normal distribution that we use for the log conductivity is characterized by the mean and
the variance. We refer to the log conductivity field as log normal or Gaussian throughout
this work.

When we characterize an aquifer using data, there are many possible outcomes due to
spatial variations of the medium between data points. IGW uses a log-normal
distribution for generating the conductivity random field. If the distribution is not log-
norrnal, the generated random field will need to be modified. For example, developing
conductivity zones where the scatterpoint data is log-normal. In order to evaluate
whether or not the distribution is Gaussian, we look at the sample statistics from the
exploratory data analysis (i.e. mean, median, mode, variance, and etc.) to help

characterize the data.

131

6.2.3.1 Exploratory Data Analysis

IGW allows the user to obtain the general statistics for a set of data. This is useful in
gathering insight into the nature of the medium. For a given data set IGW uses the
following statistics:

Mean

The sample mean is calculated as
N
2x, (6.2-8)

Variance

The variance is calculated by

0'2 =—1-%(x - 4? (6.2-9)
N 1

Median
The sample median is defined so that it is larger than half the values and smaller than the

other half. That is,

x, where] = (n + l) / 2, if n = odd
xm = , (6.2-10)
(x, +x1+1)/2 wherel=n/2, 1fn=even
Mode
The mode of a probability distribution function p(x) is defined as
120)!me = PMAX (6.2-11)

132

This implicit equation is solved by finding the maximum value of p(x,-) which gives

X mod = x j , X mod is the value of x where the probability distribution is at its maximum

value.
Average Deviation
The average deviation is calculated by
1 N —
Avg. Dev.= Félxi - X I (6.2-12)
Skewness

The skewness is calculated by

l N (x- - X}?
Skew = _ __1 6.2-13
N; a ( )
Kurtosis
The kurtosis is calculated by
N
Kurt— — {i0 ZL—x‘ —————:):X }—3 (6.2-14)
N _
i—l
PDF

IGW uses the following procedure to calculate the PDF:

1. Define Xm to be equal to MAX(x,, x2, . . . ,xg, . . ., xN) where the MAX fimction
extracts the maximum value in the associated set.

2. Define Xm to be equal to M1N(x1, x2, . . . ,xi, . . .,x~) where the MIN function extracts
the minimum value in the associated set.

3. Calculate Axas

133

Ax: Xmax -ern (6.2-15)

 

where M is the number of intervals (user specified).
4. Sum the total number of data in each interval [xj,xj+/] where j=1 , 2, . . .M. These

values are assigned to nj. Note that

M
in} =N (6.2-16)
F1
5. Calculate the PDF as
( ) "j (6 2 17)
x . = _ _ -
p 1 N
CDF
The CDF is derived from the PDF through
x
C(x) jp(x)dx (6.2-18)
Xmin
which is calculated numerically in IGW as
j
C(xj) z 202 p(x,-) (6.2-19)
i =1

Once the preliminary statistics have been calculated, we then try and categorize the
data. We categorize the data by fitting the data to a known variogram model or analytical
covariance ftmction.
6.2.3.2 Variogram analysis

We investigate the experimental and theoretical variogram in this section. The

experimental variogram is obtained fi'om the data points and then plotted using the

134

variogram function. The theoretical variogram best fits the data with a known variogram
function such as the exponential, hole-type, or Gaussian variogram.

Experimental Variogram

The experimental variogram is the graph that is most commonly used in applied
geostatistics to explore spatial interdependence. It contains information about the scale
of the fluctuations of the variable [Kitanidis, 1997]. In geostatistics, we refer to this as
the variogram or semivariogram (the semivariogram is used throughout this paper and
since we are only discussing it in this context we interchangeably use “variogram”). The
variogram by itself should not be used as a form of complete analysis. IGW uses the

following definition of the (semi-) variogram defined from the data as follows:
1 I I
70!) = 5 E [(f (x. y) - f (x ,y ”2] (62-20)
where f(x,y) and f( x', y') are the values at different locations of the random function f,

and h = "x - x'" . We refer to h as the separation distance. Some of the spatial variable

quantities that are of interest are: porosity, chemical concentration, and precipitation
[Kitanidis, 1997]. We use the IGW variogram to obtain information about the following
characteristics [Kitanidis, 1997]:
1. The presence of variability at the scale of the sampling span. This depends on the
behavior of the experimental variogram near the origin, i.e., at small separation
distances.
2. The presence of variability at a scale comparable to the sampling domain. This

depends on the behavior of the experimental variogram at large distances.

135

Theoretical Variogram

For now, we concentrate on the large-scale behavior of the system. The variogram at
distances comparable to the size of the domain determines whether the function is
stationary (wide-sense) or nonstationary. The graph of the variogram can easily
determine whether or not a certain parameter has a stationary or nonstationary behavior
(see Figure 6.7). In general, nonstationary behavior contains trends (or moving averages)

and is very prominent around sharp contrasts in the media (i.e. boundary conditions).

 

Stationary
$Hr———— ——

 

 

Range

Nugget

 

 

 

 

Figure 6.7 Variogram used to show the difference between a stationary and
nonstationary process. The figure also shows the sill, range, and nugget.

136

In the above figure (see Figure 6.7), we have the terms sill, range, nugget, stationary
and nonstationary. A stationary random function is bounded by its variance or sill as the
separation distance tends to infinity. In general, any points within this separation
distance (in the above figure, it is three separation distances) are positively correlated
with each other and this distance is given as the range. The nugget is the measurement
uncertainty scale not resolved by the data. In other words, if we were to examine two
points (i.e. conductivity values) at the same location, there would be a certain amount of

error. This we associate with the nugget effect.

With nonstationary behavior, we see that a “sill” is never obtained. So the
nonstationary sill is defined as the scale at which two measurements of the variable
become practically uncorrelated [Kitanidis, 1997]. We obtain the correlation scale or
range, a , fi'om the horizontal axis at the point where the sill or variance of the process is
obtained.

From equation (6.2-20), we see that the variogram depends on the nature of the
random firnction. The spatial function can be defined in a variety of ways. In this case,
we look at functions that can be defined by a Fourier series or summation of sines and
cosines. The connection between the spatial functions being defined by a Fourier series,
the nature of the covariance function (which is the Fourier transform of the spectral
density used in effective parameters), and the relation between covariance functions to
the variogram (recall the variogram is obtained directly from the data), provides a bridge
from the theory to real-world applications.

The variogram has direct relation to the covariance R(h) as follows

702) = —R(h) + R(O) (6.2-21)

137

The mean and covariance (first and second moments) of a random function are defined

by the following
m(x.y) =E[f(x,y)] (6.2-22)
R(h) = E[(f(x. y) — m(x.y)Xf(x'. y') — m<xzy')]
(6.2-23)
Note that for a separation distance of h = 0, we obtain the variance from the covariance

function. We refer to Kitanidis (1997) for the details between the stationary model and
the intrinsic isotropic model. In this paper, we work with the stationary model for the
variogram, which requires the sill to equal the variance, 7(00) = R(O) = 02. (The isotropic
intrinsic model does not require the sill to be the limit of the variogram as h —) 00. So,
the intrinsic model is useful in classifying the nonstationary model.)
6.2.3.3 IGW covariance models

We now give a few different covariance firnctions and their associated variogram
[Kitanidis, 1997].
Gaussian Model

For the Gaussian model we have

2 l “l
R(h) = 0' exp —-—2- (6.2-24)
L

2
702) = 0'2[1- exp[— fa] (6.2-25)
L

where (I2 > 0 and L > 0 are the two parameters of this model. The range, a , is defined

as the distance at which the correlation is 0.05, i.e. a z 7L / 4

138

Exponential Model

For this model, the covariance and variogram are given by

R(h) = 02 exp(- 1;] (6.2-26)

7(h) = 02(1- exp[- gm (6.2-27)

where the parameters are the variance 0'2 > 0 and the length parameter (or integral

scale) l> O. The range is a z 31.

Spherical Model
For the spherical model,
02 l--§—}-I—+l—,i for0 < h < a
R(h) = 2 a 2 a3 ’ - _ (6.2-28)
0, for h > a
2 3 h 1 h3
0' ————— , forOShSa
M) = 2 a 2 a3 , (6.2-29)

0' , forh>a

where the parameters are the variance 0'2 > 0 and the range a > 0.

Of the models mentioned thus far, the realizations of the random field for the Gaussian
model are the only realizations that are “smooth” enough to differentiate. The other
realizations are continuous, but not differentiable. The other noticeable difference
between the models is the exponential and spherical models each have a linear behavior

at the origin (as h —-> 0, 7(h) at h). But the Gaussian model has parabolic behavior at the

origin (ash -+ 0,7(h) oc hz).

139

6.2.4 Generation of Different Realizations

From the variogram, IGW obtains the necessary statistics for generating many
realizations that maintain the same statistical nature, namely, the variance, correlation
scale and covariance (obtained from the variogram). There are many methods of
generating the different realizations such as spectral methods, the Turning Bands method,
matrix decomposition (Cholesky decomposition), and Gaussian sequential simulation.
IGW allows the user to choose between the spectral algorithm, sequential Gaussian
simulation, and the turning bands algorithm. In addition, the simulation can be
considered unconditional or conditional. Unconditional simulation uses the statistics of
the data, but does not actually honor the data points. Conditional simulation uses both
the statistics of the data and the actual data points. Hence, the statistics are “conditioned”
on the data. We discuss four method of generating realizations, namely, matrix
decomposition, spectral method, turning bands method, and sequential Gaussian
simulation.
6.2.4.1 Basic Decomposition
The following mathematical derivation is taken from the SAS OnlineDoc Version 8
(1999). Additional comments and explanations have been provided for this mathematical
derivation. IGW does not use the matrix decomposition method. It is presented here in
its entirety because it is conceptually easy to understand.

Unconditional Simulation

We first generate a Gaussian random vector Z [n x1] with mean m and covariance C.

Let Z be represented as the following vector fixed at the locations s1, s2, . . ., S], and with

corresponding covariance matrix C

140

PZ(S1)- f C(O) C(Sr-Sz) C(SI—Sk)\
Z= z(fz) , and C: C(Sz‘Sr) C(O) C(S2-Sk) (6.2-30)
_Z(Sk)_ task—Sr) C(Sk ‘32) C(O) )

 

 

 

 

That is Z,- is a random variable, m,- is the mean of Z,, and Cg- is the covariance of Z,- and Z]-

for i, j = l, . . ., n. Then, we can obtain a realization of the vector fi'om the following
equation:
Z = m + Lu (6.2-31)

where u is an n x 1 vector of standard normal variates and L is an n x n matrix such that

LLT = C 6.2-32)
Then the lower triangular matrix L is obtained from the Cholesky decomposition.

The Cholesky decomposition is used only for positive definite symmetric matrices. Since
the correlation between Z,- and Z,- is the same as between 2,- and Z,, we have a symmetric

matrix. We know that the matrix is positive definite because the nature of the covariance

structure. That is, the matrix is diagonally dominant because the maximum correlation

value is between a point and itself. Since the matrix is diagonally dominant the

determinant is always positive. Ifthe matrix is not positive definite, then the method can

still be applied, but leads to a complex matrix L, which is impractical [Kreyszig, 1993].
For the given matrix equation

T

011 012 013 "'11 0 0 "'11 0 0

021 022 023 = "121 "122 0 "'21 "122 0 . (6-2-33)

031 “32 033 "’31 "’32 ”'33 "’31 "'32 "’33

we obtain the following formulas for the Cholesky decomposition [Kreyszig, 1993] as

follows

141

 

j-l
my =r/“jj-Zm3. j=z.-~,n
s=l
a .1 (6.2-34)
mj1=—1— 1:2,” ,n
mll

k—l
l .
mjk =— ajk— E mjsmks j=k+l,---,n
mkk 3:1

Once the lower triangular matrix L is found it is then multiplied by the normal random
vector u and summed with the mean to generate another Z random vector that preserves
the mean and covariance structure. This is repeated N times, where N is the number of
realizations specified.
Conditional simrfirtion

For a conditional simulation, the distribution of Z must be conditioned on the
observed data. The following explanation uses conditional distributions of multivariate

normal random variables.

Let X~N,,,(m, 2), where

X m 2 2
x:[ l],m=l: 1],and2=[ “ ‘2]. (6.2-35)
2321 2322

The subvectoer is kxl,X2 is nxl, 211 is kxk, 222 is nxn,and 212 =25] is
k x n , with k + n = m. The fill] vector X is partitioned into two subvectors X1 and X2,

and E is similarly partitioned into covariances and cross covariances.

With this notation, the distribution of X1 conditioned on X2 = x2 is M, (iii, 73), with

61 = m + 2122302 —m2). (6.2-36)

and

142

i = 211—2122552121. (6.2-37)
Refer to [Searle, 1971] for details. The correspondence with the conditional spatial

simulation problem is as follows. Let the coordinates of the observed data points be

denoted 31,32 ,- - - , 3”,, , with values 21, 32 ,- - - , 5". Let Z denote the random vector

 

 

"1(31 )
Z = 262) (6.2-38)
.201).
The random vector Z corresponds to X2, while Z corresponds to X1. Then
(Z [Z = 'i)~ N), (61,5) as in the previous distribution. The matrix
6 = C11 "C12C5iC21 (62-39)

is again positive definite, so a Cholesky factorization can be performed. Similarly to the
unconditional case, we obtain a lower triangular matrix L from the Cholesky
factorization. We then multiply L by the random normal vector u. and add iii to the
transformed vector. This is repeated N times, where N is the value specified for the
number of realizations.

6.2.4.2 Turning-Bands

The turning-bands method (TBM) for simulation of random fields was first presented by
Matheron (1973). Matheron named the method “turning-ban ” because of the
discretized segments (or bands) on each random line of the discrete character of the one-
dimensional generation [Bras and Rodriguez-Iturbe, 1985]. The lines turning also
correspond to the bands turning. We see applications of the method originating from the

Ecole des Mines de Paris (e.g., J oumel, 1974; Delhomme, 1979). The spectral method

143

although better for larger domains than the matrix decomposition, is less advantageous
than TBM for large number of points in the field. Instead of generating multiple fields
with equivalent statistical properties to a particular multidimensional field as done with

sampling fi'om the spectrum, TBM performs simulations along lines in space. Then for

each point in 91" , the corresponding weighted sum of values of the line process is
assigned.
We provide a brief introduction to the method as done by Bras and Rodriguez-Iturbe

(1985), which is taken from Mantoglou and Wilson (1982).
Consider a multidimensional stationary and isotropic process that is denoted by Z (-) .

Assume that the random field has zero mean and is normally distributed. (If not, we must
first provide a transformation.) We show a similar figure given by Mantoglou and

Wilson, 1982 for the two-dimensional example (see Figure 6.8). In this example, the

lines are generated fi‘om an arbitrary origin. The sampling of lines occurs at different 0,-

directions that are uniformly distributed between 0 and 2a . We then generate the mean
and covariance function C165), a one-dimensional process, along each line. The points
on the line are then projected to the region in the two-dimensional domain that is of
interest. We generate the one-dimensional process Z ,- (0) for the points along line 1'.
Figure 6.8 shows one point in space with location vector it”; its projected counterpoint on
line i is {M- and the value of the one-dimensional process at that point is Z ,- («f N,- ). We

can rewrite the projection using an inner product as follows:

Zi(:N.-)= Zi(xN we), (6.2-40)

144

where u,- is the unit vector on line i and the parentheses contain the inner product of the

vectors 1” and u,- [Bras and Rodriguez-Iturbe, 1985].

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 6.8 Schematic of field and turning-band lines (from Mantoglou and Wilson,
1982).

Once the one-dirnensional covariance C1 (5) has been generated for each of the L

lines, the point N of the two-dimensional region is assigned the simulated value

145

L
Zs(x~)= fie-(m ~u.~) (6.2-41)
i=1

This gives us the covariance of the one-dimensional process. We now need to show the
relation with the two dimensional covariance C S (v). From Mantoglou and Wilson

(1982), take two points of the field with location vectors x1 and x; as shown in Figure 6.9.

The simulated points are given by

ZS(xl)=%éZi(xl 41:) (62-42)
and
Z5(x2)=—1—ZL:Z-(x1-u).(6.2-43)
@121 J J

146

 

Figure 6.9 Projection of vector v onto a turning-bands line (from Mantoglou and
Wilson, 1981).

Recall that the one-dimensional process has zero mean. The field ZS also has zero mean

and the covariance between the points is

CS(XI’XZ)= E [71:22:01 "01%;21'61 '“j)]

(6.2-44)

2%: iE[Zi(x1'“i)Zj(x2'uj)]

i=1 j=l
Recall that the expectation of the product of two independent events is zero, unless the

events are identical. In equation (6.2-44), we can simplify the summation by noting that

147

the covariance is nonzero for i = j, and zero otherwise. Simplifying equation (6.2-44), we

obtain
1 L
CS(X1,X2)= 22512410 '“i)zi(xz '“il (6-2-45)
i=1

Note that the expression in the summation can be written as the covariance of single line

using the separation vector v as follows
E[Z,-(x1-u,-)Zi(x20ui)]=C1(v-uj), (6.2-46)
where C1 (v - u j )= C(g) is the covariance of the one-dimensional process, which we

assumed to be stationary and thus v = x2 — x1. Replacing the expression in the

summation with the covariance of the one-dimensional process gives
1 L
Cs (Krdiz)= ZZCI (V '“i)- (6-2-47)
i =1
Since the unit vector is uniformly distributed about the unit circle, we obtain
1 L
CS (v) = lim —2C1(v-u,-)= E[C1 (v - u)], (6.2-48)
L—)oo L i=1

where v = [V]. We write the continuous form of the expectation for equation (6.2-48) as

follows
E[C1(v - u)] = I: C1 (v - u) f (u)du , (6.2-49)

where f(u) denotes the probability density frmction of u and Q? the unit circle (for two-
dirnensional case) or unit sphere (three-dimensional case). The differential length or area
on Q? at the end of vector u is given by du. For the two-dimensional case, we have

f(u) = i (6.2-50)

27:.

148

Substituting the probability density function into equation (6.7-29) we obtain the

covariance function for the two-dimensional field as follows

CS(v)=2L (mi, C1(v-u)du. (6.2-51)

circle

This gives the relation between the covariance of the one-dimensional process and the
covariance for the two-dimensional field. In IGW we can obtain the variogram (and
covariance) for a given set of data points. IGW will best fit the data with a spherical,
exponential or Gaussian model. From this and using the best-fit variogram equation, we
can obtain a covariance model for the data. In practice we assume to know the

covariance function C(v) to be preserved during the simulation. So the field covariance
C S (v) is assigned the assumed covariance function. We then need to obtain the one-
dimensional covariance function C1(§)

The details can be seen in Bras and Rodriguez—Iturbe (1985). We give the relation
between the two-dimensional covariance function C(v) = C S (v) and the one-dimensional
covariance C1(v) along the turning bands in polar coordinates as follows

C S—(v) — :j— C1 (5) (6.2-52)
v2 _
where 5 = v sin 6 , d5 = vcos 316. This integral equation cannot be evaluated directly to
give an expression for C1(v) as a function of C (v) Mantoglou and Wilson (1982) give

examples of the exponential covariance and for the Bessel type of two-dimensional

covariance derived from the spectral density firnction.

149

Since IGW uses these different models, we show the exponential model, Bessel
model, and the hole-type model. These models are selected when forming a random
conductivity field.

Exponential Model
The two-dimensional model has a covariance firnction
C(v) = 025“" (6.2-53)

and the one dimensional covariance is given as

ale):011{fiance—Local}, (6.2-s4)

where 10 is a Bessel function of order zero, and L0 is a modified Struve fimction of order
zero (Abramowitz and Stegun, 1965).

Bessel Model

The two-dimensional model has a covariance function
C(v) = 02b VK1(b v) (6.2-55)

and the one-dimensional covariance is given a
C( _ 2 iiL-bé- _b: -_ -
1.9-0 1 2 51(1):) e E1( b5) , (6.256)

where Ei is the exponential integral function.
Hole-type Model

The hole-type model is used to fit a model of the exponential or Bessel processes.
Mantoglou and Wilson (1981) show that the two-dimensional covariance function in the
turning-bands method that relates to the one-dimensional hole covariance frmction is

given by

150

C(v)=0'2{Io(av)—L0(av)+av[11(av)—L1(av)-%]}, (6.2-57)

where 11 and L1 are Bessel and Struve functions of order 1 and a is the parameter of

C1 (4‘) in the following equation

C1(§) = 0'2 (1 — a§)e"a§ o s g < oo. (6.2-58)
Since this covariance firnction has similar behavior to both the exponential and Bessel
processes, we fit these processes by the hole-type covariance function.
Generation of realizations
The generation of realizations along the turning-bands lines uses a moving-average

process with the weighting function defined as

205(1—aw)e_aw wZ 0
0 w<0

f(w) = ( (6.2-59)

As in the spectral approach, the covariance function of the one-dimensional process or
along the bands is the discrete covariance function. We refer the reader to Bras and
Rodriguez-Iturbe (1985) for the full detailed explanation for conditional and
unconditional simulations.
6.2.4.3 Spectral Algorithm

The spectral method requires a little more background in understanding the method,
since the simulations and preserving of statistical properties is done in the frequency
domain, instead of in the spatial domain. We first give a general outline to the approach,

and then explain the necessary theory for each of the steps in accordance with IGW 3

151

User's Manual [Paulson, 2002]. The beginning of the spectral approach is obtaining the
covariance structure from the scatterpoint data created by IGW or entering field data.
Using the variogram analysis IGW best fits the data with a covariance model (i.e.
exponential, hole-type, etc.). It then transforms the covariance to the frequency domain
using the Fast Fourier Transform (F FT ). IGW is then able to easily generate different
stochastic processes that are statistically equivalent. We give the mathematical
descriptions of the processes below.

The stationary covariance structure between two spatial quantities from a zero-mean

2D stochastic process h(r1,r2) = h(r) is given as
R(rl , 22 ) = R(r) = 15[h(r1 + 2'] , r2 + 72 W0] (6.2-60)
where m is the conjugate function of h(r). The spectrum of the stochastic process,

S(k1,k2) = S(k), can be expressed as the FFT of the covariance function as follows
S(k) = ”R(‘r)exp[27ti(k - 1)]dr (6.2-61)

and
R(r)= 117:3- ”S(k)exp[— 272i(k-r)]dk (6.2-62)

We then write the above covariance and spectral relationship in their discrete form using

 

the discrete FFT
N1 -1 N2 -1 m n m n
Sum = Z Z Rm1.m2 exv[2m'( 13,1 + If, 2 H (6.2-63)
"11:01:12 =0 l 2
1 Nl—lNz—l .
le,m2 =W;N_2 Z anbnz exp[— 27110111an A1 +man2 A2)] (6'2'64)
I!

1=0n2=0

152

where N1, N2 = the 2D domain dimensions [L],
k = the frequency [LI],

2' = the spatial quantity [L], and

A = the grid spacing [L].
Since we will be employing this method to generate many logarithmic hydraulic
conductivity fields, we use f = In K values as the spatial quantities from the 2D stochastic
process. We now replace equation (6.2-64) with the notation for the log hydraulic
conductivity field.
If the fm is a stochastic process, the covariance between any two points m1 and m2 can be

written as follows
1 N1 -1N2 -1
fmbm2 =—-—— Z Eran,”2 exp[— 27zi(m1k,,lA1 +m2kn2A2)]. (6.2-65)
N1N2 _ _
n1 —0 n2 —0
The stochastic processes F in the frequency domain are random under the assumption
that fm is a stochastic process. Assuming that the spectral density F represents a wide-

sense stationary process, that is,

E[F,,l,,,2 ] = 0 , (6.2-66)
E[Fn,,, me] = 0 , where n.- ¢ m,- (6.2-67)

where I"? is the conjugate function of F. Taking the expectations of equation (6.2-65) and
multiplying by the conjugate 7 , we obtain the following relationship
Eanlmz Fn1,n2 ] = N1N2Sn1,n2 2 (6-2‘68)

which relates to the spectral density given earlier. So the objective is to first generate a

stochastic process, for example, a log hydraulic conductivity random field. From this we

153

can randomly sample the field and then IGW will find the best fit covariance structure for
the data. IGW then uses the FFT to obtain the random spectrum in the fi'equency domain.
This is where we generate many random spectrums that have the same statistical
properties.
6.2.4.4 Sequential Gaussian Simulation

The Sequential Gaussian Simulation can be found in the GS Library manual and we

refer the user to that manual.

6.2.5 Recursive Statistical Post-Processing

For each of the realizations created the mean head, log conductivity, and concentration
are obtained. The post-processing of all the realizations allows one to obtain the overall
statistics of the simulation such as the probability density function (PDF), the cumulative
distribution function (CDF), the mean, the variance, and the covariance. For the
collection or ensemble of all the realizations we use recursive formulas for the statistical
moments that are listed below.

The IGW field-based statistics are a little more challenging to describe in that they are
based upon the stochastic process or random field. We let the random field be noted
f(x, y). IGW uses recursive formulas for many realizations or stochastic processes (as
done with the Monte Carlo method), so that it only has to store the data for the current
realization. It would take up too much memory of the computer to store each realization
and the post-process the statistics. Recursive formulas are equivalent to the definitions

provided above, but are made efficient for multiple realizations.

154

6.2.6 Monte Carlo Flow and Transport Modeling

The Monte Carlo method is used to generate many different solutions that are
statistically equivalent and also have a component of randomness. This allows a model
to incorporate variability between data points.
6.2.6.1 Basic Idea

The basic idea for the Monte Carlo method is to generate many random In K
realizations. We then solve the flow equation for each realization and obtain flow
realizations. Once we obtain the flow realizations, we use this information to generate
the plume realizations by solving the transport equation.
6.2.6.2 New Information from Monte Carlo

After generating the many different likely scenarios or realizations, IGW computes the
ensemble means (minimum variance estimate) and confidence intervals (variances). We
are then able to analyze the statistics and obtain the uncertainty and risk-based
predictions (probability) for the flow and transport in a heterogeneous aquifer.
6.2.6.3 Monte Carlo Simulation Procedure

We go through the Monte Carlo Simulation step-by-step in Chapter 7: Illustrative
Applications of Stochastic Groundwater Theory. The first example in that section is the
Monte Carlo method.

The general procedure is as follows:

1. Create the model set-up with boundary conditions, sources/sinks, and monitoring
wells.

2. Create the conductivity random field by using the "Random" choice in IGW.

3. Run the model and then sample the field to obtain hypothetical data points. If field

data is available, input that data into a spreadsheet format.

155

4. Then choose the type of simulation to be completed using the data (i.e. conditional
simulation, unconditional, etc.)

5. Use the exploratory data analysis to see whether or not the data is Gaussian and
whether or not the data is stationary. In order, to determine stationarity a sufficient
amount of data points is necessary.

6. Select the Monte Carlo method as the type of solution and the number of realizations.
Also, IGW allows the user to choose the number of visual outputs for each realization.
7. Observe the post-statistics such as the concentration breakthrough curve at the
monitoring wells. The PDF, CDF, and stochastic process are also available for the log
conductivity, head, and concentration. Also, the means and variances can be obtained
graphically as well as computationally giving results of uncertainty, confidence intervals,
and risk prediction.

6.2.7 Recursive Formulas for Ensemble Statistics

Recursive Mean

The recursive mean is calculated by

_ "(k-l) _
[700 (x, y) = F (x, y)(: 1) + f(xs y) , k =1929H-Nk (62-69)

where k is the realization index number. The total number of realization, Nk, is set by the
user and can theoretically approach infinity.
Recursive Variance

The recursive variance is calculated by

 

(k—I) _ _ we 2
0.00 (x, y) = 0' (x, J’Xk 1) + [[(x9 y) F (x9 y)] k =1,2,"., Nk (6.2-'70)

156

Recursive Covariance

Considering a second random field, g(x, y), the covariance between point Pj(x, y) and

Pg(x, y) is expressed as

 

 

I I “(‘0
e e k—l

_ “(16) _ “(’0
+ [f(x.y) F (x.y)JLgl(cxo,yo) G 0000)] k =1.2.---.Nk

(6.2-71)

157

7 Illustrative Applications of Stochastic Groundwater Theory

7.1 Probability Model

In the 19605, Warren and Price (1961) and Shvidler (1964) presented applications of
the Monte Carlo method to flow in porous media. Although the Monte Carlo method was
seen at this earlier time, it wasn’t popular until Freeze (1975) used the Monte Carlo
method to study both steady-state and transient flow in a one-dimensional domain, where
the hydraulic conductivity is assumed to be statistically independent in adjacent blocks
[Zhang, 2002]. “Freeze’s (1975) work is arguably the beginning of the stochastic era for
flow in porous media [Zhang, 2002].”

We first explain conceptually the Monte Carlo approach and then we show the
mathematical set-up of the method. Conceptually, the Monte Carlo method is statistical
sampling. Statistical sampling is completed on different realizations (or paths) that are
statistically equivalent. The realizations are generated using a variety of methods, which
we explore next.

As mentioned earlier in Chapter 5, the statistical sampling approach to modeling
heterogeneity in groundwater flow allows us to use both random fields and field data.
Even though the data may be sparse, we can form a conditional simulation, which honors
the data at those points. In this innovative example, we first create a randomly generated
conductivity field with a mean of 100 fi/day, log conductivity variance of 1, and with a
correlation scale of 10 m. We then sample 150 data points (see Figure 7.1). From these
points, we use exploratory analysis of the scatterpoints to obtain additional statistical

information. We then choose a conditional simulation and show the associated

158

variogram with the points. The variogram (see Figure 7.2) shows that range is about 45
m and the sill is located at 0.61. The correlation scale is obtained from the range, and in
this case the correlation scale is 75 m. We also obtain that the log hydraulic conductivity
can be represented as a stationary process because the variogram is bounded with
increasing separation distance.

We slightly change the set-up from the previous contaminant transport model by
adding a lake in the center of the domain. This is done for illustrative purposes of
showing the effects of an internal boundary that creates nonstationarity within the model.

We find that nonstationarity changes the distribution of the plume.

159

 

\4

onitorin- Well 2

'l,

g.

,.

' 000.00 feet

 

‘ _.
‘0;

Figure 7.1 Flow field with representative scatterpoints used for a conditional Monte
Carlo simulation.

160

 

Model direction 1 (angle-90)
Experiment data direction 1

 

 

 

 

 

 

Variogram
9
WW
9 O Q ..
300 400 500 000 700
Distance (11)

Figure 7.2 Variogram of scatterpoints

IGW gives the best-fit solution for the range, sill, and nugget. The nugget is the value
given close to the origin. In this case, the nugget is assigned zero as the variogram is
predicted to go through the origin itself. This represents the small-scale variability. We
notice that for this particular sampling the data is bounded and so we can assume that the
process is stationary within the domain of sampling.

We obtain insight into an appropriate model for the data fi'om the variogram. In this
example, exponential model best fits the data with range 45 meters and variance 0.61.
Later this becomes critical in the choice of method for the generation of realizations. For
instance, the spectral method requires a covariance structure that will represent the data.
7.1.1.1 Different realizations and static profiles

To contrast the mean behavior, we show some of the results from individual
realizations that vary considerably from one to another. We also show the associated
static profile for each of the realizations shown below to give an idea of how the

concentration distribution changes considerably in the different realizations. The first

161

realization we show with the log conductivity. The remaining realizations we show just

the plume to better visualize the differences.

-
“'1
Monitorinv Well 2
Static Profl

Transient Flow, Time Elapsed - 0 days (0.00 years), In Realization 100

 

Figure 7.3 Realization 100 plume migration at the end of 4800 days.

162

 

 

[Static Profile Grapifl

 

 

 

500,

i

8.00

I; ’ \

Sam

H

C

8200

C

8100'
"0 100 200 300 m

Lengm (m)

 

 

 

Figure 7.4 Static profile from original location to Monitoring Well 1.

163

@m' Wel

 

 

 

2000.00 feet 4000.00 0000.00

Transient Flow, Time Elapsed = 4800 days (13.16 years), In Realization 104

Figure 7.5 Realization 104 plume migration at the end of 4800 days.

164

..r

m.\-.

 

 

 

 

,0, [Static Profile Graph]
C
I0 m
E
‘E 200
o
0
C
o 100
0
30 100 200 300 400
Length (m)

 

 

 

Figure 7.6 Static profile from original location to Monitoring Well 1.

 

 

 

 

Figure 7.7 Realization 104 plume migration at the end of 4800 days.

165

2 . a 'L': '1. .. . .5“ .22 ' . . . ‘ . . 4‘
W ~
. > - ' , I

 

 

 

 

a, [Static Profile Graph]
.§ m
E ..J
5...
o
C
0 so
0
Length (m)

 

 

 

Figure 7.8 Static profile from original location to Monitoring Well 1.

 

 

 

 

000.00 feet 000.00 00.00
[ Transienti-‘lowfl'lmeslapeed-m daysttatsyem), lnReallaallentW [

Figure 7.9 Realization 107 plume migration at the end of 4800 days.

166

 

Figure 7.10 Static profile from original location to Monitoring Well 1
The results above are all very different and yet, they are still statistically equivalent. We
need to be able to characterize the mean and variance for the ensemble of all realizations
to get an estimate of the most probable locations that the plume will travel. IGW
provides he user with a "means and variances" graph that is processed for all realizations.
7.1.1.2 Means and Variances

We next show the means and variances of the plume and conditional field after 100
realizations (see Figure 7.11). This is very import for practitioners, in that; one can
obtain a confidence interval to predicting the location of a plume migrating. One of the
major advantages of stochastic modeling is that multiple realizations can be used to show
an exhaustive process of many possibilities that could happen with equal probability.
The accumulation of all such realizations gives a mean behavior (ensemble mean) and

predicted variance. However, caution should be taken in that the means and variances

167

describe a Gaussian process. If the distribution of a parameter is not Gaussian, then the

means and variances may be inaccurate.

 

 

Figure 7.11 Means and Variances after 100 realizations.

Many statistics can be calculated once the Monte Carlo simulation is complete. The
process, probability density function, and cumulative density function are all obtainable
from the simulation for the different specified points of interest. In this case, we
specified two different monitoring wells. We will take a look at some of the interesting
distributions, but not all the distributions for each well. The objective in calculating the

different distributions is to show that we cannot assume that the head, in K, or

168

concentration distributions are Gaussian. In particular, internal boundary conditions have
an effect on the nature of the distribution.
7.1.1.3 Head distributions

We first look at the probability density frmction for the head at each of the monitoring
wells (see Figure 7.12 and Figure 7.13). The density functions for Monitoring Well 2
(located closest to the initial location of the plume) resembles that of a normal probability
distribution. But the density firnction for Monitoring Well 1 (located at the edge of the
pond) is slightly skewed (the peak is shifted to the right instead of being centered). The
primary reason for the change in the density function is the proximity to the pond

(internal boundary).

 

 

 

 

Figure 7.12 Head at Monitoring Well 2, away from internal boundary.

169

‘ Probability at Monitoring Well 1

 

Figure 7.13 Head at Monitoring Well 1, near the internal boundary.

7.1.1.4 Concentration distributions

We next look at the concentration density fimctions at the monitoring wells. In the
contaminant transport section, we discussed the concentration of the plume being
represented as a Gaussian distribution. We predict that the internal boundary condition
will change the distribution, and that the distributions away from the boundary conditions
should still represent a Gaussian distribution. We observe the following results (see
Figure 7.14 and Figure 7.15). Neither of the monitoring wells show a Gaussian
distribution, which means the means and variances will not give an accurate description
of the concentration. Monitoring Well 1 shows a first a peak and then greatly reduces,
similar to an exponential disu'ibution. Monitoring Well 2 actually shows two peaks (first

a smaller peak and then a slightly larger peak)

170

   

g, Probability of Miiiiiluriiiq Well 1

 

lr-‘IEII

Min WE
Men WEI—J
Medan WE
Mode WE
Ave Err WE

 

 

 

 

 

 

 

 

 

 

 

 

 

10000 2001» Sid
Me «p. 9:11; {$1. SWE
Kutosis WE]
i smkkanmcmmv' W 1" Pieces Process Cur-...»; Choice
[2 Heed m (CM '7 Ham [7 News”
- CDF 17 Mom [.7 Mom-Std
mam am ]

 

 

 

 

 

Figure 7.14 Concentration in a non- Gaussian distribution at lake boundary.

   

 

A Probability .il Moriiliiriiiq Wu" Z

 

Wm
Men [THE
Medm WE]
Mode WE]
AveEir WE

 

 

 

10000 310.00 300m
m

 

I101004E
0.1034 E
Ku‘tosia [1.1889 E

 

 

 

 

 

 

(“Process
GPDF
(“CDF

SeieclaPaiameteitoVisueize
r c .

‘ “‘l

 

 

Process Cur-re Choice
[7 Hefidion l7 MmStd
17 Mean I7 Mom-Std

 

 

 

 

 

Figure 7.15 Concentration in a non-Gaussian

l7l

distribution away from the lake.

We see that the results of the concentration distributions as the different wells are very
different. Yet, neither one of them is a Gaussian-like distribution. It poses the question
of what type of distribution does concentration have with variability, if any? We have
seen effects of decay, partitioning, and dispersivity in a homogeneous aquifer with the
distributions showing some similarity to the Gaussian distribution (whether shifted,
reduced peak, etc.) We would expect that chemical properties of the pollutant or even the
chemical interaction between both the medium and the pollutant to yield non-Gaussian
behavior.
7.1.1.5 Log conductivity distribution

We next observe the probability distribution functions for the logarithmic
conductivity. As noted above in the generation of realizations, we use a random normal
vector to generate the conductivities in between the data points. So one would expect
that the distributions of the log conductivity to resemble a Gaussian distributions. At
Monitoring Well 2 (furthest from the lake), we obtain a primarily Gaussian distribution
are 1n K =3.0, which corresponds to a conductivity value of 20.08 m/day. We used an

average value of 20 m/day to form the In K field.

172

l’rulidblllty at Monitorinri Well 2

. ___l
‘ -_..|
_J
.'__J
__i

 

Figure 7.16 Ln K distribution at Monitoring Well 2.

Observing the In K probability distribution firnction at Monitoring Well 1, which is
closer to the boundary of the lake, we find that the distribution still resembles a Gaussian
distribution. This is due to the generation of the field itself. Even though the region
around a lake would primarily have a higher conductivity, we used the same average
conductivity value of 20 m/day in the region near the lake. The reason is that the lake
does not firlly penetrate the aquifer, and we are capturing the heterogeneity of the aquifer.
If the lake were to have penetrated deep into the aquifer, then a deterministically defined
conductivity would have to be assigned to the region of the lake. We would then have to

resample the region of interest.

173

Probability at Monitoring Well 1

 

Figure 7.17 Ln K distribution at Monitoring Well 1 (near the lake).

In summary, we have seen that the distribution change for the head and concentration
when there is an internal boundary condition. This can easily be predicted in that the
flow of both the discharge from the river and the contaminant are interrupted by the lake
boundary. The conductivity, however, was not affected as long as the lake boundary
does not penetrate deep into the aquifer. If it does, then the conductivity zones will need
to be added to capture the change in soils.
7.1.1.6 Concentration breakthrough curves
We continue with our analysis by observing the concentration breakthrough curves at the
different monitoring wells. The concentration breakthrough curve shows an average
breakthrough time of about 2500 days (point of inflection) (see Figure 7.18). The mean
and standard deviation are given for the well. We also can visualize a particular

realization with respect to the mean and standard deviations. This is very powerful in

174

that we can determine a range of time when the concentration for the plume to enter the

well. Similar results show for Monitoring Well 1, but with a later breakthrough time.

llllll'l'ilill'xx-ll~1Iilll'llllllll-IUI‘“_‘

1”" _ ____.r _ ..

 

Figure 7.18 Concentration breakthrough curve at Monitoring Well 2.

175

8 Conclusion

The bridge between theory and applications begins with replicating the basic concepts of
groundwater flow, contaminant transport, and stochastic flow and transport, in a visual
numerical setting. The visual software that we used was the software package Interactive
Groundwater (IGW). We gave examples of local, intermediate, and regional flow by
recreating the Toth solution. Then we recreated the Freeze and Witherspoon results with
the added anisotropy and heterogeneity to the Toth solution. For contaminant transport,
we showed the affects of decay, partitioning, and dispersion in migrating plume. We also
showed the correlation between the grid resolution and numerical dispersion. Finally, we
showed the basic concepts of representing the porous medium as a random field
(stochastic flow). We showed the effects of small-scale variability on a migrating plume.
We also investigated using effective parameters (or average parameter values) versus
representing the heterogeneity of the medium using a random field.

We extended beyond the basic concepts with innovative examples such as flow
beneath a dam, transient lake with recharge, groundwater seepage using drains,
hierarchical modeling in a well field, and conditional Monte Carlo simulation using the
fast Fourier spectral approach.

Every piece of the numerical model starting from the hardware to the algorithms used
plays a critical role in the efficiency of a software package. We have shown in the above
examples yet another perspective of the computational bottleneck, namely, the
visualization of results. In this case, much of the inefficiency associated with the

visualization has been removed by the computational structure of IGW. This is not

176

present in other software packages such as GMS and MODF LOW. Although results can
be seen within minutes, there is still a need for efficient numerical algorithms so that the
interactive model becomes "continuous" in delivering results. The crux of this paper has
been to provide researchers, practitioners, and instructors with examples that will aid in
conceptualizing the dynamics of groundwater. We provide final comments regarding the

issues that are currently present as applied to researchers, practitioners, and instructors.

8.1 Research Issues

Some of the current issues that are present in groundwater research are the
computational bottleneck, representing nonstationary dynamics, and reducing numerical
error.

The application of spectral methods and perturbative methods in representing the
variable porous medium has helped considerably in resolving the high computational
time for modeling variability. But, yet much of the stochastic applications assume that
the fluctuations of the log conductivity field and hydraulic head values can be represented
by stationary processes. Yet, there are many instances when this cannot be assumed (i.e.
internal boundary conditions, trends within the medium). Li and McLaughlin (1991 ,
1995) established a nonstationary spectral method to better deal with these changes in the
medium. The theoretical results have proven to give better resolution in the problematic
regions, but the practicality of the application in the numerical generation of the solutions
was hindered by the computational time. The results can be verified with analytical
solutions that are complex and tedious, and only exist under simplifying assumptions.
Current research for IGW includes a new method that uses the stationary spectral method

to analyze the general behavior. Then in areas that are problematic (i.e. internal

177

boundary conditions), numerical modeling with fine grids are used to resolve the flow
pattern about the complex geometry of the domain. There is much progress being made
in exploring new theories and their application in groundwater modeling. The impacts
will only draw us closer to being able to model the real-world with a three-dimensional

computationally efficient software package.

8.2 Practitioners Issues

From the practitioner's perspective, groundwater modeling needs to be correlated with
the given data taken from well logs, pump tests, and other data collection techniques.

The underlying question that is still prominent amongst groundwater modelers is whether
or not deterministic or stochastic modeling should be the method of choice. Gelhar
(1993) addresses this same issue in his book Stochastic Subsurface Hydrology. The
answer is that it should depend on the objective of the project. In some instances, when
flow prediction is the only desired result, then deterministic modeling may be adequate.
When considering contaminant transport, then stochastic modeling may be appropriate.
So the issue at hand is determining when each type of modeling should be used.

Another issue that is still problematic is that the field can change drastically within a
few steps of each data collection point. In order to capture heterogeneity in a model, we
have to consider the Heisenberg uncertainty principle. That is, where the measurement
device interferes with the fundamental measurement. In every subsurface class, we study
the different ways to obtain the fundamental parameter values, for example, determining
conductivity of a certain soil type. From analyzing soil samples in the lab to tracer
studies in the field, the determination of the hydraulic conductivity lends itself to very

elaborate procedures. The bottom line is that the characterization of data incorporated

178

into groundwater model has many limitations and assumptions. It is impossible to both
physically and economically represent a project site with enough data.

This leads to using calibration techniques that often use well log data. When simply
modeling the general flow pattern in a basin, this method can give valid results. But as
mentioned before, if contaminant transport is included, then this may not be adequate.

Where do we go fi'om here? From a practical standpoint, software such as IGW needs
to be stressed to incorporate stochastic techniques with data points. Using "real-time
steering", which allows the scientist or engineer "to be an equal partner with the computer
in manipulating and maneuvering the 3D visual presentations of the modeling results [Li
and Liu, 2004]." It gives the user control over much of what has been considered a
"black box" in the past. For instance, the modeler can "interactively steer the
computation, control the execution sequence of the program, guide the evolution of the
subsurface flow and plume migration dynamics, control the visual representation of data
during processing, and dynamically modify the computational process during its
execution [Li and Liu, 2004]." In instances where the location of the plume is critical,
probability modeling should be used since it gives confidence intervals of the location of
the plume at a desired time. Moreover, with legal issues rising, especially in groundwater
lawsuits, visual software needs to be incorporated in the courtroom to educate the jurors,
judges, and lawyers. In particular, the placement of boundary conditions (and how the

model can vary) is critical in the outcome of the flow pattern.

8.3 Instructional Issues

The instructional issues that we face today stem from presenting both theory and

applications to students that have a wide variety in background. The actual subject

179

matter of groundwater flow uses partial differential equations and the theory is
considered to be at a graduate level because of the underlying equations. Darcy's Law is
the main equation that gets the majority of focus at the undergraduate level, and yet the
solution to that equation can lead to very complex theories. However, in training students
both at the graduate and advanced undergraduate levels, educational stress should be on
conceptualizing the physical dynamics of the system. This can be done in pointing out
the relationships between parameters, without actually solving the equations. Knowing
whether or not a software package is giving valid solutions is the most pressing issue
facing groundwater modeling in the classroom today.

Instructional software such as IGW lets students build a hypothetical real-world where
the physical dynamics can be singled out amongst the theory. The advances that need to
be made in groundwater modeling lie in asserting that students conceptualize many
models and be able to predict the results. IGW is a "sketch book" for very complex
dynamical groundwater situations. Students can develop their own hypothetical
situations and quickly see results. This will allow students to see the necessary
components of a model (i.e. boundary conditions, water table location for a vertical
model, sources/sinks, etc.). Finally, after forming an intuitive picture of groundwater
dynamics, advanced studies should be given in building a data dependent model. This
will provide the student with a better background to decide whether or not the model is

accurate.

180

9 References

Abramowitz, M. and I. A. Stegun (1965), Handbook of Mathematical Functions, Dover.

Afshari, S. (2003 ), Improved Finite Difference Methods for Modeling Two-Dimensional
Flow in General Anisotropic Media. Master’s Thesis, Department of Civil and
Environmental Engineering, Michigan State University.

Anderson, MP. and Woessner, W.W. (1992). Applied Groundwater Modeling;
Simulation of F low and Advective Transport. San Diego, California: Academic Press.

Bakr, AA, L.W. Gelhar, AL. Gutjahr, and JR. MacMillan (1978), Stochastic analysis
of spatial variability in subsurface flows, 1. Comparison of one- and three-dimensional
flows, Water Resour. Res. , 14(2), 263 —272.

Bear, J. (1979), Hydraulics of Groundwater. McGraw-Hill, Israel.

Bras, RL. and I. Rodriguez-Iturbe (1985), Random Functions in Hydrology, Addison-
Wesley.

Fetter, C.W. (2001). Applied Hydrogeology. Upper Saddle River, New Jersey: Prentice-
Hall.

Freeze, RA and Cherry, IA. (1979). Groundwater. Englewood Cliffs, New Jersey:
Prentice-Hall.

Freeze, RA and Witherspoon, RA (1966). Theoretical Analysis of Regional
Groundwater Flow: 1. Analytical and Numerical Solutions to the Mathematical Model.
Water Resources Research 2, no. 4:641-56.

Freeze, RA and Witherspoon, PA (1967), Theoretical Analysis of Regional
groundwater Flow: 2. Effect on Water-Table Configuration and Subsurface Permeability
Variation. Water Resources Research 3, no. 22623-34.

Gelhar, L.W. (1993), Stochastic Subsurface Hycb'ology. Englewood Cliffs, New Jersey:
Prentice-Hall.

Gelhar, L.W., and CL. Axness (1983), Three-dimensional stochastic analysis of
macrodispersion in aquifers, Water Resour. Res, 19(1), 161-180.

Gutjahr, A, B. Bullard, S. Hatch, and L. Hughson (1994), Joint conditional simulations
and the spectral approach for flow modeling, Stock Hydro]. Hydraul, 8(1), 79-108.

181

H. Liao, K.J. Paulson, S.G. Li, C.F. Ni (2003), IGW 3 Reference Manual, Department of
Civil and Environmental Engineering, Michigan State University.

Joumel, AG., and CI. Huijbregts (1978), Mining Geostatistics, Academic Press,
London.

Kitanidis, PK. (1997), Introduction to Geostatistics: Applications to Hydrogeology.
Cambridge, UK, Cambridge University Press.

Kreyszig, E. (1993), Advanced Engineering Mathematics (Seventh Edition), New York,
John Wiley & Sons, Inc.

Li, S.—G. and D. McLaughlin (1991), A nonstationary spectral method for solving
stochastic groundwater problems: Unconditional analysis, Water Resourc. Res, 27(7),
1589-1605.

Li, S.-G. and D. McLaughlin (1995), Using the nonstationary spectral method to analyze
flow through heterogeneous trending media, Water Resourc. Res, 31(3), 541-551.

Li, S.-G., D. McLaughlin, and H. Liao (2003), A computationally practical method for
stochastic groundwater modeling, Advances in Water Resources, 26, 1137-1148.

Li , S.-G. and Q. Liu (2004), Interactive Ground Water (IGW): An innovative digital
laboratory for groundwater education and research, Computer Applications in
Engineering Education, 1 1(4).

Li, S.-G. and Q. Liu (2004), A New Paradigm for Groundwater Modeling, In review by
Ground Water, October 2004.

Lunrley, IL, and HA Panofsky (1964), The Structure of Atmowheric Turbulence, John
Wiley & Sons, New York

Manna, M, AK. Bhattacharya, S. Choudhury, and SC. Maji (2003), Groundwater flow
beneath a sheetpile analyzed using six-noded triangular finite elements, The Institution of
Engineers (India): Technical Journals, Civil Engineering, 080309.

Mantoglou, A, and J .L.Wilson (1982), The turning bands method for simulation of
random fields using line generation by a spectral method, Water Resour. Res, 18, 1379-
1394.

Matheron, G. (1973), The intrinsic random functions and their applications, Adv. App].
Prob, 5, 438-468.

Mizell, S.A, AL. Gutjahr, and L.W. Gelhar (1982), Stochastic analysis of spatial

variability in two—dimensional steady groundwater flow assuming stationmy and
nonstationary heads, Water Resour. Res. , 18(4), 1053-1067.

182

 

‘ . .
r
.. D
.. . .
.J \

 

.
..l x
s
'
I!
§ . . .u ,
...
It‘
..
e
I‘
i:
l: a
a
'
l!
.
‘t.
a i
a
s .
o
‘1
u.
.
I
v
if. .
. A:
a r
..i.
. i
. ..

i . . it:
.
V
l A
\. . .n
. .
t. a .-
li .
b,
:4
. .
i. . u b
L
f
.
i .
y

 

Paulson, K]. and 8.6. Li (2002), Interactive Groundwater Users Manual, Department of
Civil and Environmental Engineering, Michigan State University.

Polaha, M.V. Cracking Dams .' An Interactive Web Site for K-12. [Online]
http://www.simscience.org/cracks/thesis.htm1

Potter, MC. and Wiggert, DC. (2002). Mechanics of Fluids. Pacific Grove, California:
Brooks/Cole.

Rubin, Y., and G. Dagan (1988), Stochastic analysis of boundaries effects on head spatial
variability in heterogeneous aquifers: 1. Constant head boundary, Water Resour. Res,
24(10), 1689-1697.

Rubin,.Y., and G. Dagan (1989), Stochastic analysis of boundaries effects on head spatial
variability in heterogeneous aquifers: 2. Impervious boundary, Water Resour. Res, 25(4),
707-712.

Shvidler, M.I. (1964), Filtration Flows in Heterogeneous Media, Consultants Bureau,
New York.

SAS OnlineDoc, Version 8 (1999), The computational and theoretical details of spatial
simulation, Cary, NC: SAS Institute Inc.

T oth, J. (1962). A theory of groundwater motion in small drainage basins in central
Alberta. J. Geophys. Res, 67, 4375-4387.

Toth, J. (1963). A theoretical analysis of groundwater flow in small drainage basins. J.
Geophys. Res, 68, 4795-4812.

Walpole, Myers and Myers. (1998) Probability and Statistics for Engineers and
Scientists. Sixth edition. Upper Saddle River, New Jersey: Prentice Hall, Inc.

Warren, J .E., and HS. Price (1961), Flow in heterogeneous porous media, Soc. Petrol.
Eng. J., 1,153-169,1961.

Zhang, D. (2002). Stochastic Methods for Flow in Porous Media. San Diego, California:
Academic Press. '

Zheng, C., (1990). M T 3D, A modular three-dimensional transport model for simulation
of advection, dispersion and chemical reactions of contaminants in groundwater systems,

Report to the US. Environmental Protection Agency, 170 p.

Zheng, C. and Bennet, G. D. (1995). Applied Contaminant Transport Modeling. New
York: Van Nostrand Reinhold.

183

10 APPENDICES

184

10.1 Appendix A

Appendix A is written as the “how-to” documentation of the IGW examples given in the
main text. They are to aid in both using the software as well as for instructional
purposes. The appendix is also written such that it can serve independently as a manual
for students to use.

We use the following steps as a general step-by-step guideline for the flow and transport
models problems. The stochastic examples are done separately:

1.

2.

9

Domain: Determine the size of domain needed and whether or not the units will
be in meters or feet.

Zones and Parameter Values: Construct the zones needed and assign parameter
values. To replicate the examples above we provide the coordinates for defining
each zone. For instructional purposes, we suggest that the instructor roughly
estimate the zones. IGW has been designed to give quick illustrations that allow
changes. Also, when defining the zones, IGW gives the option of calculating
zone mass balances.

Boundary Conditions: Construct boundary conditions, i.e. no-flow, constant
head, etc. This step can also include the generation of a particular flow pattern.
Flow from left to right can be constructed using two different stages in rivers at
the left and right boundaries, respectively.

Internal Sources/Sinks: Add internal sources/sinks (i.e. recharge zones,
lakes/rivers, drains). With the sources/sinks decide whether transient or steady-
state is appropriate

Run Flow Model/Add Particles: The determination of the flow pattern can be
done by discretizing and running the model. If also considering contaminant
transport in the model, continue with step 6. IGW does not require that the flow
and transport be'solved separately.

Contaminant Zones: Add contaminants either by constructing a zone that is
either an instanteous or continuous source. This is done similarly to the internal
sources and sinks. Or construct a line or zone of particles.

Profile Model: Add cross-sectional models by defining profile (polylines) along
paths of interest.

Transient or Steady-State: For particles and concentration, the model will run in
transient mode. For flow one must choose one of the processes (river, drain,
recharge, etc.) to be transient.

Grid Size: Discretize and run the model.

10. Static Profiles: Define static profiles along paths of interest.
11. Water Balance: View water balance results.

185

10.1.1 Seepage beneath a Dam.

Inactive zone

  
 
     

Shootpile

5186“" 1°” Constant Head-
20 feet

Constant Head=
200 feet

  
   

500 feet

K-8.64 m/day

 

1000 feet ..

 

 

 

 

 

Figure 10.1 Triangular dam with sheetpile

1. Domain: 1000 it x 1000 it.
2. Zones and Parameter Values: Constructed the following zones (with coordinates)
and assign parameter values:

Porous Medium (IGW currently only allows units of meters to be entered):

In order to keep units of feet, we convert feet to meters for entry of coordinates.
(0,0), (304.8, 0), (304.8, 152.4), (190.5, 152.4), (190.5, 137.2), (114.3, 137.2),
(114.3, 152.4), (0,152.4, end)

Parameter values: We assigned a conductivity value of 8.64 m/day.

Triangular Dam:

(114.3, 137.2), (114.3, 228.6), (120.4, 228.6), (190.5, 152.4), (190.5, 137.2, end)
This is an inactive zone.

Reservoir:

(0,152.4), (114.3, 152.4), (114.3, 213.4), (0,213.4, end)

186

Shallow Controlled Discharge:
(190.5, 152.4), (304.8, 152.4), (304.8, 158.5), (184.9, 158.5, end)

We also opted to calculate the zone mass balance for the porous medium. This allows us
to figure out the volume of seepage for a particular thickness.

3. Boundary Conditions: Boundary conditions were assumed to be no-flow for the left,
right, and bottom. The triangular dam is considered an “inactive” zone. The reservoir is
assigned a constant head boundary of 200 feet, and similarly, the shallow controlled
discharge assigned a constant head of 20 feet.

4. Internal Sources/Sinks: There were no internal sources or sinks for this example.

5. Run Flow Model/Add Particles: The flow pattern can be obtained at this step, but in
order to visualize the flow lines, we assign a polyline with particles evenly spaced. The
solution for this problem is a steady-state solution. Since we “activate” particles for
visualization, the solution is obtained using “transient” mode.

6. Contaminant Zones: We constructed a particle polyline by using the following
endpoints:
(0, 148.32) and (114.3, 148.32). We then assigned 12 particles along the line.

7. Profile Model: Since this is a cross-sectional or vertical model, we do not define a
profile model.

8. Transient or Steady-State: Since we have added particles for visualization, the
solution is obtained using “transient” mode.

9. Grid Size: We discretized the model with 200 x 200 cells, which gives a cell
dimension of 1.5 meters x 1.5 meters, approximately. When using particles whether for
visualization or for contamination, the model must be finely discretized.

10. Static Profiles: No static profiles were done for this model.

11. Water Balance: The zone mass balance can be seen at the end of the simulation.
IGW has a default thiclmess of 50 meters. As mentioned before the average crest length
is approximately 1500 it, we use a width of 1,640 ft for this size of dam. So we need to
multiply the seepage rates by a factor of 10 to obtain the total seepage rate. We
calculated the zone mass balance without the sheet piling and then with the sheet piling.
For the sheet piling, we added an “inactive” zone with the following coordinates and
repeated steps 1-11.

Sheetpile: (114.3, 137.2). (117.3, 137.2), (117.3, 106.7), (114.3, 106.7, end).

187

10.1.2 Toth Solution: Local Model and Regional Model.

1. Domain: 20,000 it x 2,000 ft (Local)
20,000 ft x 15,000 it (Regional)
2. Zones and Parameter Values: Constructed the following zones with coordinates:

Porous Medium (IGW currently only allows units of meters to be entered):
In order to keep units of feet, we convert feet to meters for entry of coordinates.
Values for this were generated using a MATLAB program (see Appendix B) and
are specific to the undulations (see text for analytical equation). In this case, we
have 20 = 1000 feet, u = 72/100 , a' = 100 and b' = .011 . The list of specific input
values is given as the 41 x 2 matrix T in Appendix B. For instructional purposes,
a roughly drawn figure will yield similar results.
Parameter values: We assigned a conductivity value of 25 m/day.

 

ll .11 I

 

Figure 10.2 Porous medium

Inactive zones:
The inactive zones were formed along the top of the porous medium to represent
the surface. The undulations actually represent the water table, not the surface. At
certain places the water table may meet the surface (i.e. surface seepage,
wetlands, etc.)

3. Boundary Conditions: Boundary conditions were assumed to be no-flow for the left,
right, and bottom. The “caps” or unsaturated zone is considered an “inactive” zone. In
this vertical model, the water table is the upper boundary condition. The water table is
formed using the polyline feature, and we used the values for matrix T in Appendix B
(evenly spaced points in the sinusoidal variation). IGW assumes that the water table
boundary is fixed in the vertical model.

4. Internal Sources/Sinks: There were no internal sources or sinks for this example.
5. Run Flow Model/Add Particles: The flow pattern can be obtained at this step, but in
order to visualize the flow lines, we assign a polyline with particles evenly spaced. The

solution for this problem is a steady-state solution. Since we “activate” particles for
visualization, the solution is obtained using “transient” mode.

188

6. Contaminant Zones: We constructed a particle polyline by lowering the vertical
coordinate of each undulation point by 30 meters (90 feet)). We then assigned 30
particles along the line.

7. Profile Model: Since this is a cross-sectional or vertical model, we do not define a
profile model.

8. Transient or Steady-State: Since we have added particles for visualization, the
solution is obtained using “transient” mode.

9. Grid Size: We discretized the model with 200 x 100 cells, which gives a cell
dimension of 30.63 meters x 6.16 meters, approximately. When using particles whether

for visualization or for contamination, the model must be finely discretized.

10. Static Profiles: No static profiles were done for this model.

189

10.1.3 Freeze and Witherspoon.

For this example, we use the general set-up given for the isotropic case. This contains the
general outline for the porous medium and the unsaturated zone. Then for the different
cases, we will reference the changes.

1. Domain: 300 meters x 100 meters

2. Zones and Parameter Values: Constructed the following zones (with coordinates):
(U = upper zone, L = lower zone)

Porous Medium:

Case 1: Isotropic (Kh=10 m/day; Kh:Kv=1 :1)

Case 2: Anisotropic (Kx=10 m/day; Kr:Ky=10:l)

Coordinates(Cases 1&2): (0, 0), (300, 0), (300, 52), (35, 50), (0, 48, end)

 

Surface

 

 

fi/// /// / / // /// //

Porous Medium

Figure 10.3 Freeze and Witherspoon Anisotropic Case 2

Case 3: Low upper conductivity zone. (K U=l m/day; KL=10 m/day)
Coordinates: Lower zone (0, 0), (300, 0), (300, 40), (0, 40, end)
Upper zone (0, 40), (300, 40), (300, 52), (35, 50), (0, 48, end)

 

i Surface
I

 

W// ////////////
UnerZone

    

Lower Zone

Figure 10.4 Freeze and Witherspoon Anisotropic Case 3

Case 4: Very high upper conductivity zone. (KU=100 m/day; KL=1 m/day)
Coordinates(Cases 4 & 5): Lower (0, 0), (300, 0), (300, 20), (0, 20, end)
Upper (0, 20), (300, 20), (300, 52), (35, 50), (0, 48, end)

190

Case 5: Very high lower conductivity zone. (K U=10 m/day; KL=1000 m/day)

 

Surface

 

_,_,—~7-7'77//////////////////////////
UpperZone

Figure 10.5 Freeze and Witherspoon Anisotropic Cases 4 and 5

 

Case 6: High, low, high conductivity zones horizontally.

(Kleft=100 m/day; Kmiddle=10 m/day; Kright=100 m/day)
Coordinates: Left zone (0, 0), (60, 0), (120, 50.6375), (35, 50), (0, 48, end)
Middle zone (60, 0), (100, 0), (200, 51.2375), (120, 50.6375, end)

Right zone (100, 0), (300, 0), (300, 52), (200, 51.2375, end)

/—Su rfoce

g //////
Left Zone Middle , Right Zone
Zone -; 3L

Figure 10.6 Freeze and Witherspoon Anisotropic Case 6

 

 

 

 

Inactive zones:
The inactive zones were formed along the top of the porous medium to represent
the surface.

3. Boundary Conditions: Boundary conditions were assumed to be no-flow for the lefi,
right, and bottom. The unsaturated zone is considered an “inactive” zone. In this vertical
model, the water table is the upper boundary condition. The water table was assigned
using the polyline feature and the following points:

Coordinates: (0, 48), (35, 50), (300, 52, end)
IGW assumes that the water table boundary is fixed in the vertical model.

4. Internal Sources/Sinks: There were no internal sources or sinks for this example.

191

5. Run Flow Model/Add Particles: The flow pattern can be obtained at this step, but in
order to visualize the flow lines, we assign a polyline with particles evenly spaced. The
solution for this problem is a steady-state solution. Since we “activate” particles for
visualization, the solution is obtained using “transient” mode.

6. Contaminant Zones: We constructed a particle polyline by lowering the surface
vertical component of each point by 2 meters. We then assigned 30 particles along the
line.

7. Profile Model: Since this is a cross-sectional or vertical model, we do not define a
profile model.

8. Transient or Steady-State: Since we have added particles for visualization, the
solution is obtained using “transient” mode.

9. Grid Size: We discretized the model with 100 x 85 cells, which gives a cell
dimension of 3.03 meters x 1.19 meters, approximately. When using particles whether
for visualization or for contamination, the model must be finely discretized.

192

10.1.4 Transient Lake.

1. Domain: 1000 meters x 750 meters

2. Zones and Parameter Values: Constructed the following zones (with coordinates):

Porous Medium (represent with one zone in the horizontal plane):

Used scatterpoints to create a topographical surface with a lake in cross-sectional
view. Note: The lake is formed by the topography, not the actually assigning of
cells. The following is a table of the scatterpoints used from left to right in the
domain. (We do not give the exact coordinates for each of the scatterpoints, but

this can be done for critical projects where data is known.)
Parameter Values: Conductivity = 0.001 cm/sec

 

\ . ‘3
‘ \
- I I i i ‘\. ’
\‘ \
m‘ -. \ I ‘( . 'l .
. .\ .
u - ~ ~ . . \'
400.00 . . . -
S 102 ~
/
. ., ,2
. . . . '7. . . ”I: .
rv—I'rr‘f'”: . ,/ .

 

. , 1 -
'f'200.00(m

 

 

Transient Flow. Time Elapsed = 10 days (0.03 years)

 

Figure 10.7 Lake in cross sectional view

193

 

 

 

 

 

 

 

 

Scatterpoint Surface (m) Top (m) Bottom (m)
102 35 35 -50
103 45 45 -52
104 -l 5 -15 -45
105 2.5 2.5 -45
106 -8 -8 -42
107 25 25 -50
108 40 4O -55

 

 

 

 

 

 

Table 10.1 Scatter points

3. Boundary Conditions: Boundary conditions were assumed to be no-flow for the left,
right, bottom and top. The unsaturated zone is considered an “inactive” zone.

4. Internal Sources/Sinks: The internal sources or sinks for this example include the
river/lake and recharge.
Parameter Values: Recharge = 10 in/yr

River = Transient

Sediment Conductivity = 0.1 m/day

River Bottom = same as surface

5. Run Flow Model/Add Particles: The flow pattern can be obtained at this step. In this
case, we do not use particles to depict the flow line. Since we defined the river/lake as
“transient”, the solution is obtained using “transient” mode.

6. Contaminant Zones: No particles used.

7. Profile Model: We defined a profile model. In the profile model, the water table was
obtained from the horizontal model. We created a profile model shown above (the line
through the scatterpoints).

8. Transient or Steady-State: The solution is obtained using “transient” mode.

9. Grid Size: We discretized the model with 100 x 75 cells, which gives a cell

dimension of 10.10 meters x 10.13 meters, approximately. When using particles whether
for visualization or for contamination, the model must be finely discretized.

194

10.1.5 Surface Seepage

1. Domain: 1000 meters x 750 meters
2. Zones and Parameter Values: Constructed the following zones (with coordinates):

Porous Medium (represent with one zone in the horizontal plane):

Used scatterpoints to create a topographical surface with a lake in cross-sectional
view. Note: The seepage regions are formed by the topography, not the actually
assigning of cells. The following is a table of the scatterpoints used from left to
right in the domain. (We do not give the exact coordinates for each of the
scatterpoints, but this can be done for critical projects where data is known)
Parameter Values: Conductivity = 0.001 cm/sec

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.. , ,.’f’.ZSPI33 221351" . ._ _
.-- -- .(Dm .. 400.00.— ...600.00......m0.00u—.-.1
mum Flew,11rns‘Elapssrd=20 am passions)

Figure 10.8 Scatter points for topography and seepage

195

Scatterpoint Surface (m) Bottom (m)
122 45
123 10
124 15
125 25
126 45
127 6
128 25
129 -2
130 30
131 50
132 2
133 40
134 -4
135 35
136 60
137 60
138 60

 

Table 10.2 Scatter points for topography

3. Boundary Conditions: Boundary conditions were assumed to be no-flow for the left,
right, bottom and top. The unsaturated zone is considered an “inactive” zone.

4. Internal Sources/Sinks: The internal sources or sinks for this example include the
drain and transient recharge.

Parameter Values: Recharge = Transient
Drain = Same as surface elevation

5. Run Flow Model/Add Particles: The flow pattern can be obtained at this step. In this
case, we do not use particles to depict the flow line. Since we defined the recharge as
“transient”, the solution is obtained using “transient” mode.

6. Contaminant Zones: No particles used.

7. Profile Model: We defined two profile models. In the profile models, the water table
was obtained from the horizontal model. The profile models were defined diagonally
across the domain as shown above.

8. Transient or Steady-State: The solution is obtained using “transient” mode.

9. Grid Size: We discretized the model with 50 x 37 cells, which gives a cell dimension
of 10.27 meters x 10.45 meters, approximately.

196

10.1.6 Pumping Well

1. Domain: 26.25 h x 52.49 ft (8 m x 16 m)
2. Zones and Parameter Values: Constructed the following zones (with coordinates)
and assign parameter values:
Middle Aquifer Layer: 1.8 m x 8 m (This zone needs to correspond with the well
screen since it is the only layer turned on).
Porous Medium (IGW currently only allows units of meters to be entered):
In order to keep units of feet, we convert feet to meters for entry of coordinates.

Parameter values: We assigned a conductivity value of 103.68 m/day in the
lower aquifer layer. The upper layers were assigned as “inactive”.

scatterpoint -> , d- scatterpoint

 

59.57

well 7 g ;
\ INACTIVE .

 

 

 

 

well
screen

INACT IVE

INACTIVE

INACT IVE .
scatterpoint * scatterpornt

Figure 10.9 Pumping well

197

3. Boundary Conditions: Scatterpoints were used to depict the thickness of the aquifer
expanding radially. The thickness was assigned by using the “top of aquifer” and
“bottom of aquifer” inputs. Scatterpoints at the edge of the model were given a thickness
value of 50.265 meters and in the center 0.314 meters. Boundary conditions were
assumed to be no-flow for the bottom and top of the lower aquifer. The cast and west
boundaries are also no-flow boundaries, except where the well screen is located.

4. Internal Sources/Sinks: We assigned the lower well zone a negative recharge value
of -96.4 m/day. The area of the zone is 0.09 m2 (0.05 x 1.8). For a pumping rate of -
8.676 m3/day.

5. Run Flow Model/Add Particles: The flow pattern can be obtained at this step. The
solution for this problem is a transient-state solution. Since we calculate the zone mass
balance over time, we run the model in transient mode.

6. Contaminant Zones: Particles were not used in this model.

7. Profile Model: Since this is a cross-sectional or vertical model, we do not define a
profile model.

8. Transient or Steady-State: The solution was obtained using the “transient” mode, so
that we can obtain the cumulative volume of the water pumped in the water balance.

9. Grid Size: We discretized the model with 161 x 161 cells, which gives a cell
dimension of 0.05 meters x 0.1 meters, approximately. When using particles whether for

visualization or for contamination, the model must be finely discretized.

10. Static Profiles: No static profiles were done for this model.

198

10.1.7 Well Head Protection Area

1. Domain: 7200 ft x 5400 ft.
2. Zones and Parameter Values: Constructed the following zones (with coordinates)
and assign parameter values:

Porous Medium (IGW currently only allows units of meters to be entered):
In order to keep units of feet, we convert feet to meters for entry of coordinates.

Plume zone: (914.4, 457.2), (944.88, 444.95), (975.36, 365.76), (944.88, 286.57),
(914.4, 274.32), (883.92, 286.57), (853.44, 365.76), (883.92, 444.95)

These points are located on an ellipse with semimajor axis =300fl and semiminor
axis = 200 ft, centered at (914.4, 365. 76).

South river zone:
(0, O), (O, 60.96), (2194.56, 60.96), (2194.56, 0, end)

North river zone:
(0, 1584.96), (0, 1645.92), (2194.56, 1645.92), (2194.56, 1584.96, end)

Parameter values:
Conductivity = 100 ft/day
Porosity = 0.3

3. Boundary Conditions: No flow for the lefi and right boundaries. Constant head
boundaries for the top and bottom boundary conditions.

199

 

 

 

 

fl}: mw.n&;__;;
4000.00 . ‘...‘.:.':..;;:;
.................... ,I- :.I-ulll‘.....2.
uraIW-ellllw Run“ [...Z: ....

1:21:21 w.u2 """

2000.00(fee 1:112: 111 i .........
i:::n u'alwell3 :‘i: ...... IIIIIIIII'
00000166 ----- 000.00 0.00'“

 

mew. 11m Elm-dam duel-1813M!)
Figure 10.10 Well protection area

4. Internal Sources/Sinks: We have the following pumping rates (point source/sink) for
the wells:

Town Well = -500 gpm

Rural Well 1 = -100 gpm

Rural Well 2 = -100 gpm

Rural Well 3 = -200 gpm

South River: Stage = 5 meters
Leakance = 5 /day
Bottom Elevation = -30 meters

North River: Stage = 0 meters
Leakance = 5 /day
Bottom Elevation = -30 meters

5. Run Flow Model/Add Particles: The flow pattern can be obtained at this step. The

solution for this problem is a transient solution, since we place particles around the wells.
We run the model backwards in time to obtain the wellhead protection area.

200

6. Contaminant Zones: Particles were placed around the well, using the default number.
7. Profile Model: We did not define a profile model.

8. Transient or Steady-State: The solution was obtained using the “transient” mode.

9. Grid Size: We discretized the model with 30 x 40 cells, which gives a cell dimension
of 51.44 meters x 52.25 meters, approximately. When using particles whether for
visualization or for contamination, the model must be finely discretized. We used a time

step of 10 days.

10. Static Profiles: No static profiles were done for this model.

201

10.1.8 Random Walk

1. Domain: 1000 it x 500 ft.

2. Zones and Parameter Values: Constructed the following zones (with coordinates)
and assign parameter values:

Porous Medium (IGW currently only allows units of meters to be entered):
In order to keep units of feet, we convert feet to meters for entry of coordinates.

Parameter values: Conductivity = 20 m/day

Particle Zone: (100, 200), (100, 300), (200, 300), (200, 100, end)

 

4760 days

 

 

.999

 

   

 

 

 

 

 

 

 

 

 

 

400 I] 0(m) . 800 .00

 

Figure 10.11 Random walk

4. Internal Sources/Sinks: The plume region is considered an “instantaneous” source.
We used polylines for the constant head boundaries on the left and right sides.

Left polyline = 2 meters
Right polyline = 0 meters

5. Run Flow Model/Add Particles: The flow pattern can be obtained at this step. The
solution for this problem is a transient-state solution, since we observe the
plume/particles moving.

6. Contaminant Zones: Particles were used by picking “Random Walk” as the
contaminant transport method.
Local dispersivity: Longitudinal = 1 meter
Transversal = 0.1 meter
Vertical = 0.0 meter

202

7. Profile Model: We did not define a profile model.

8. Transient or Steady-State: The solution was obtained using the “transient” mode.

9. Grid Size: We discretized the model with 50 x 25 cells, which gives a cell dimension
of 20.4 meters x 20.8 meters, approximately. When using particles whether for
visualization or for contamination, the model must be finely discretized.

10. Static Profiles: No static profiles were done for this model.

11. Water Balance: The zone balance was not calculated for this example.

203

10.1.9 Numerical Dispersion and Hierarchical Modeling

1. Domain: 7200 ft x 5400 ft.
2. Zones and Parameter Values: Constructed the following zones (with coordinates)
and assign parameter values:

Porous Medium (IGW currently only allows units of meters to be entered):
In order to keep units of feet, we convert feet to meters for entry of coordinates.

Plume zone: (914.4, 457.2), (944.88, 444.95), (975.36, 365.76), (944.88, 286.57),
(914.4, 274.32), (883.92, 286.57), (853.44, 365.76), (883.92, 444.95)

These points are located on an ellipse with semimajor axis =3 00 ft and semiminor
axis = 200 ft, centered at (914.4, 365.76).

 

South river zone:
(0, 0), (O, 60.96), (2194.56, 60.96), (2194.56, 0, end)

North river zone:
(0, 1584.96), (0, 1645.92), (2194.56, 1645.92), (2194.56, 1584.96, end)

Parameter values:
Conductivity = 100 ft/day
Porosity = 0.3

204

 

...... ' @W)

0Rurd Well 1 Rural

Wel 2

 
   
  
   
  

2000.00neeq ___,..._.~/

\8......j..

 

  

 

200mm") 4000.00 6000.00

 

Steady Flow, Time Elmer! -:2040 dive 15.5:Oeyeers)

Figure 10.12 Numerical dispersion

4. Internal Sources/Sinks: We have the following pumping rates (point source/sink) for
the wells:

Town Well = -500 gpm

Rural Well 1 = -100 gpm

Rural Well 2 = -100 gpm

Rural Well 3 = -200 gpm

South River: Stage = 5 feet
Sediment conductivity = 0.1 m/day
Sediment thickness = River bot elev. - aquifer top elev.
Bottom Elevation = -30 feet

North River: Stage = 0 feet
Sediment conductivity = 0.1 m/day
Sediment thickness = River bot elev. — aquifer top elev.
Bottom Elevation = -30 feet

205

5. Run Flow Model/Add Particles: The flow pattern can be obtained at this step. The
solution for this problem is a transient-state solution, since we observe the
plume/particles moving.

6. Contaminant Zones: We have one contaminant zone, which we represented as a
plume and with particles. We used MMOC for the plume, and forward particle tracking
for the particles. Both zones used the coordinates for the “plume zone” listed above.

7. Profile Model: We did not define a profile model.
8. Transient or Steady-State: The solution was obtained using the “transient” mode.

9. Grid Size: We discretized the parent model with 100 x 77 cells, which gives a cell
dimension of 21 .6 meters x 21.7 meters, approximately. We discretized the child model
with 100 x 120 cells, which gives a cell dimension of 4.5 meters x 4.5 meters,
approximately. When using particles whether for visualization or for contamination, the
model must be finely discretized. We then ran the model and viewed the sub-model
output of the child model, which is more accurate.

10. Static Profiles: No static profiles were done for this model.

11. Water Balance: The water balance was not calculated for this example.

206

10.1.10 Hierarchical Modeling

1. Domain: 7200 ft x 5400 ft.
2. Zones and Parameter Values: Constructed the following zones (with coordinates)
and assign parameter values:

Porous Medium (IGW currently only allows units of meters to be entered):
In order to keep units of feet, we convert feet to meters for entry of coordinates.

Plume zone: (914.4, 457.2), (944.88, 444.95), (975.36, 365.76), (944.88, 286.57),
(914.4, 274.32), (883.92, 286.57), (853.44, 365.76), (883.92, 444.95)

These points are located on an ellipse with semimajor axis =3 00 ft and semiminor
axis = 200 ft, centered at (914.4, 365.76).

South river zone:
(0, 0), (0, 60.96), (2194.56, 60.96), (2194.56, 0, end)

North river zone:
(0, 1584.96), (0, 1645.92), (2194.56, 1645.92), (2194.56, 1584.96, end)

Parameter values:
Conductivity = 100 ft/day
Porosity = 0.3

207

 

 

 

 

       
    

!"""""""'3
‘Rural Well 1 i
—--—--:.-————+.fiw6n 2/

   

 

2000.00 feet '

Figure 10.13 Hierarchical modeling

4. Internal Sources/Sinks: We have the following pumping rates (point source/sink) for
the wells:

Town Well = -500 gpm

Rural Well 1 = -50 gpm

Rural Well 2 = -100 gpm

Rural Well 3 = -200 gpm

South River: Stage = 5 feet
Sediment conductivity = 0.1 m/day
Sediment thickness = River bot elev. — aquifer top elev.
Bottom Elevation = -30 feet

North River: Stage = 0 feet
Sediment conductivity = 0.1 m/day
Sediment thickness = River bot elev. — aquifer top elev.
Bottom Elevation = -30 feet

208

 

5. Run Flow Model/Add Particles: The flow pattern can be obtained at this step. The
solution for this problem is a transient-state solution, since we observe the
plume/particles moving.

6. Contaminant Zones: We have one contaminant zone, which we represented as a
plume and with particles. We used MMOC for the plume, and forward particle tracking
for the particles. Both zones used the coordinates for the “plume zone” listed above.

Plume concentration = 550 ppm (instanteous source)

7. Profile Model: We did not define a profile model.

8. Transient or Steady-State: The solution was obtained using the “transient” mode.

9. Grid Size: We discretized the parent model with 100 x 77 cells, which gives a cell
dimension of 21 .6 meters x 21.7 meters, approximately. We discretized the child model
with 200 x 97 cells, which gives a cell dimension of 4.4 meters x 4.5 meters,
approximately. We discretized the grandchild model with 160 x 82 cells, which gives a
cell dimension of 2.8 meters x 2.9 meters. When using particles whether for visualization
or for contamination, the model must be finely discretized. We then ran the model and
viewed the sub-model output of the child model, which is more accurate.

10. Static Profiles: No static profiles were done for this model.

11. Water Balance: The water balance was not calculated for this example.

209

10.1.11 Advection only, decay, partitioning, and dispersion of a migrating plume

1. Domain: 7200 it x 5400 ft.
2. Zones and Parameter Values: Constructed the following zones (with coordinates)
and assign parameter values:

Porous Medium (IGW currently only allows units of meters to be entered):
In order to keep units of feet, we convert feet to meters for entry of coordinates.

Plume zone: (914.4, 457.2), (944.88, 444.95), (975.36, 365.76), (944.88, 286.57),
(914.4, 274.32), (883.92, 286.57), (853.44, 365.76), (883.92, 444.95)

These points are located on an ellipse with semimajor axis =3 00 fl and semiminor
axis = 200 fl, centered at (914. 4, 365.76).

South river zone:
(0, O), (0, 60.96), (2194.56, 60.96), (2194.56, 0, end)

North river zone:
(0, 1584.96), (0, 1645.92), (2194.56, 1645.92), (2194.56, 1584.96, end)

Parameter values:
Conductivity = 100 ft/day
Porosity = 0.3

210

 

    

2000.00(feet)

“annuarmrms

. :99

 

 

000.0(feet) 4000.0 6000.0

Steady Flow, Time Elapsed = 1120 days (3.07 years)

 

Figure 10.14 Advection only, decay, partitioning, and dispersion of a migrating
plume

4. Internal Sources/Sinks: We have the following pumping rates (point source/sink) for
the wells:

Town Well = -500 gpm

Rural Well 1 = ~100 gpm

Rural Well 2 = -100 gpm

Rural Well 3 = -200 gpm

South River: Stage = 5 feet
Sediment conductivity = 0.1 m/day
Sediment thickness = River bot elev. - aquifer top elev.
Bottom Elevation = -30 feet

North River: Stage = 0 feet
Sediment conductivity = 0.1 m/day
Sediment thickness = River bot elev. — aquifer top elev.
Bottom Elevation = -30 feet

211

5. Run Flow Model/Add Particles: The flow pattern can be obtained at this step. The
solution for this problem is a transient-state solution, since we observe the
plume/particles moving.

6. Contaminant Zones: We used MMOC for the numerical method of tracking the
plume.

Plume concentration = 500 ppm (instantaneous source)

Advection only: No additional parameter values required.
Decay: Decay = 0.005 /day
Partitioning: Partitioning coefficient = 1 mL/ g
Dispersion: Longitudinal = 100 m

Transversal = 5 m

7. Profile Model: We did not define a profile model.

 

8. Transient or Steady-State: The solution was obtained using the “transient” mode.

9. Grid Size: We discretized the parent model with 100 x 77 cells, which gives a cell
dimension of 21.6 meters x 21.7 meters, approximately. We only show the advection
case.

10. Static Profiles: Static profiles were obtained for each example. The data was then
exported to an excel file. A graph with the different effects was then created.

 

212

10.1.12 Deterministic heterogeneity using random walk and multiple models

Multiple Models: Create a template model (i.e. model that has the basic same
parameters and boundary conditions), then click on the “multiple models” button in IGW.
In this case we formed two models. Since random walk uses particles and plume
representation, we show two of the same multiple models (side by side) (see Figure
below).

1. Domain: 1000 m x 500 m.
2. Zones and Parameter Values: Constructed the following zones (with coordinates)
and assign parameter values:

Porous Medium (IGW currently only allows units of meters to be entered):
In order to keep units of feet, we convert feet to meters for entry of coordinates.

Parameter values: Conductivity (outside of lenses and fingering) = 20 m/day

Fingering: Arbitrary coordinates were chosen.
Conductivity (in fingering formation) = 0.1 m/day

Lenses: Arbitrary coordinates were chosen.
Conductivity (in lenses) = 0.1 m/day

Particle Zone: (100, 200), (100, 300), (200, 300), (200, 100, end)

 

(Q ‘1'.” ‘ '

”9‘2... W
’ .Qlfiw ,
dais-r" ‘ W

 

 

20000} ”\{\

l i | l, \,\~.l) \ '
i l 4%)? E l . m
Steady Flow, Time Elapsed - 5560 days (15 .23 years)

l)? :fi\ 10“” ”ff“
A

 

 

 

 

 

 

 

Figure 10.15 Deterministic heterogeneity using random walk and multiple models

213

4. Internal Sources/Sinks: The plume region is considered an “instantaneous” source.
We used polylines for the constant head boundaries on the left and right sides of each
model.

Left polyline = 2 meters
Right polyline = 0 meters

5. Run Flow Model/Add Particles: The flow pattern can be obtained at this step. The
solution for this problem is a transient-state solution, since we observe the

plume/particles moving.

6. Contaminant Zones: Particles were used by picking “Random Walk” as the
contaminant transport method.

7. Profile Model: We did not define a profile model.

8. Transient or Steady-State: The solution was obtained using the “transient” mode.
9. Grid Size: We discretized the model with 200 x 201 cells, which gives a cell
dimension of 5 meters x 5 meters, approximately. When using particles whether for

visualization or for contamination, the model must be finely discretized.

10. Static Profiles: No static profiles were done for this model.

214

 

10.1.13 Log hydraulic conductivity represented as a variable random field.

Multiple Models: Create a template model (i.e. model that has the basic same
parameters and boundary conditions), then click on the “multiple models” button in IGW.

In this case we formed four models.

1. Domain: 1000 m x 500 m.

2. Zones and Parameter Values: Constructed the following zones (with coordinates)

and assign parameter values:

Porous Medium (IGW currently only allows units of meters to be entered):
In order to keep units of feet, we convert feet to meters for entry of coordinates.

Parameter values:

Random Conductivity (starting in the upper left hand corner and proceeding
clockwise) = 20 m/day (given average value for each model)

Upper Left Corner

Random => Scale 1
Spectral Algorithm
,1, = 1y =10
(Variance) 0'2= l
Exponential Model
Angle = 90
Nugget = 0.01

Scale 2
Spectral Algorithm
A, = 1y = 100

(Variance) 02 = O
Exponential Model
Angle = 90
Nugget = 0.01

The overall variance is then obtained by summing the variances for each scale. In
this case (upper left hand corner), we have a variance of 1.

Upper Right Corner

Random => Scale 1
Spectral Algorithm
,1, = 1y =10

(Variance) 0'2 = 1
Exponential Model
Angle = 90
Nugget = 0.01

Scale 2
Spectral Algorithm
xix = 11y = 100

(Variance) 0'2: 1
Exponential Model
Angle = 90
Nugget = 0.01

The overall variance is then obtained by summing the variances for each scale. In
this case (upper left hand comer), we have a variance of 2.

Lower Right Corner
Random => Scale 1

Spectral Algorithm

Scale 2
Spectral Algorithm

215

lx=ly=10 lx=ly=100

(Variance) 02= 1 (Variance) 0'2= 2
Exponential Model Exponential Model
Angle = 90 Angle = 90

Nugget = 0.01 Nugget = 0.01

The overall variance is then obtained by summing the variances for each scale. In
this case (upper left hand comer), we have a variance of 3.
Lower Left Corner

Random => Scale 1 Scale 2
Spectral Algorithm Spectral Algorithm
Ax=ly=10 Ax=ly=100
(Variance) 0'2 = 1 (Variance) 0'2 = 3
Exponential Model Exponential Model
Angle = 90 Angle = 90
Nugget = 0.01 Nugget = 0.01

The overall variance is then obtained by summing the variances for each scale. In
this case (upper left hand comer), we have a variance of 4.

Particle Zone: (100, 200), (100, 300), (200, 300), (200, 100, end)

 

  

mono A m
’ " .~ , Steady Flow. Time Elapsed - 1810 days (4.96 years)-

 

Figure 10.16 Log hydraulic conductivity

216

4. Internal Sources/Sinks: The plume region is considered an “instantaneous” source.
We used polylines for the constant head boundaries on the left and right sides of each
model.

Left polyline = 0 meters
Right polyline = -2 meters

5. Run Flow Model/Add Particles: The flow pattern can be obtained at this step. The
solution for this problem is a transient-state solution, since we observe the
plume/particles moving.

6. Contaminant Zones: We picked 50 particles in the zone and used forward particle
tracking.

7. Profile Model: We did not define a profile model.

8. Transient or Steady-State: The solution was obtained using the “transient” mode.
9. Grid Size: We discretized the model with 200 x 100 cells, which gives a cell
dimension of 10.2 meters x 10.2 meters, approximately. When using particles whether

for visualization or for contamination, the model must be finely discretized.

10. Static Profiles: No static profiles were done for this model.

217

10.1.14 Macrodispersion versus local dispersion using effective parameters

Multiple Models: Create a template model (i.e. model that has the basic same
parameters and boundary conditions), then click on the “multiple models” button in IGW.
In this case we formed two models.

1. Domain: 1000 m x 500 m.
2. Zones and Parameter Values: Constructed the following zones (with coordinates)

and assign parameter values:

Porous Medium (IGW currently only allows units of meters to be entered):
In order to keep units of feet, we convert feet to meters for entry of coordinates.

Parameter values:

Upper Model:

Conductivity = 20 m/day

Porosity = 0.3

Local dispersivity (macrodispersion for this case):
Longitudinal = 20 meters
Transversal = 3.25 meters
Vertical = 0 meters

Lower Model:
Random Conductivity = 20 m/day (given average value for each model)
Random => Scale 1 Scale 2

Spectral Algorithm Spectral Algorithm

xix=ly=10 xix=2y=10

(Variance) 0'2= 1 (Variance) 0'2: 1

Exponential Model Anistropic Bell Model

Angle = 90 Angle = 90

Nugget = 0.01 Nugget = 0.01

The overall variance is then obtained by summing the variances for each scale. In
this case (lower model), we have a variance of 2.

Porosity = 0.3

Local dispersivity (macrodispersion for this case):
Longitudinal = 10 meters
Transversal = 1 meters
Vertical = 0 meters

218

 

Particle Zone (used in template model): (100, 200), (100, 300), (200, 300),
(200, 100, end)

suadyFlolemeElapeed-Mdayanzsm

 

Figure 10.17 Macrodispersion versus local dispersion using effective parameters

4. Internal Sources/Sinks: The plume region is considered an “instantaneous” source.
We used polylines for the constant head boundaries on the left and right sides of each
model.

Left polyline = 0 meters
Right polyline = -2 meters

5. Run Flow Model/Add Particles: The flow pattern can be obtained at this step. The
solution for this problem is a transient-state solution, since we observe the
plume/particles moving.

6. Contaminant Zones: We picked 50 particles in the zone and used forward particle
tracking.

7. Profile Model: We did not define a profile model.

8. Transient or Steady-State: The solution was obtained using the “transient” mode.

219

9. Grid Size: We discretized the model with 100 x 100 cells, which gives a cell
dimension of 10.1 meters x 10.2 meters, approximately. When using particles whether
for visualization or for contamination, the model must be finely discretized.

10. Static Profiles: No static profiles were done for this model.

220

 

10.1.15 Waste pond model with individual fractures

1. Domain: 2000 ft x 700 fi.
2. Zones and Parameter Values: Constructed the following zones (with coordinates)
and assign parameter values:

Porous Medium (IGW currently only allows units of meters to be entered):
In order to keep units of feet, we convert feet to meters for entry of coordinates.

Parameter values:

Porous medium zone (outside of sand lenses):

(39.624, 30.48), (82.296, 167.64), (112.776, 176.784), (213.36, 176.784), (228.6,
173.736), (243.84, 164.592), (304.8, 152.4), (365.76, 164.592), (381.0, 173.736),
(396.24, 176.784), (609.6, 176.784), (609.6, 30.48, end)

Random Conductivity = 10'7 cm/sec (given average value-glacial till)

Random => Scale 1 Scale 2
Spectral Algorithm Spectral Algorithm
Ax=ly=10 xix=zly=10
(Variance) 0'2= 1 (Variance) 0'2= 1
Exponential Model Anistropic Bell Model
Angle = 90 Angle = 90
Nugget = 0.01 Nugget = 0.01

The overall variance is then obtained by summing the variances for each scale. In
this case (lower model), we have a variance of 2.

Porosity = 0.4

Local dispersivity (macrodispersion for this case):
Longitudinal = 10 meters
Transversal = 1 meters
Vertical = 0 meters

Waste pond zone (contaminant zone):

(228.6, 173.736), (243.84, 164.592), (304.8, 152.4), (365.76, 164.592), (381.0,
173.736, end)

Lake: (0,0), (0,167.64), (82.296, 167,64), (39.624, 30.48), (3048,1524, end)

Sand lenses:
Upper left lens (ellipse centered at 700, 350):

221

 

(152.4, 106.68), (182.88, 113.28), (213.36, 114.3), (243.84, 113.28), (274.32,
106.68), (243.84, 100.08), (213.36, 99.06), (182.88, 100.08, end)

Lower Middle lens (ellipse centered at 800, 250):
(182.88, 76.2), (213.36, 82.80), (243.84, 83.82), (274.32, 82.80), (304.8, 76.2),
(274.32, 69.60), (243.84, 68.58), (213.36, 69.60, end)

Upper right lens (ellipse centered at 1200, 375):
(304.8, 114.3), (335.28, 120.90), (365.76, 121.92), (396.24, 120.90), (426.72,
114.3), (396.24, 107.70), (365.76, 106.68), (335.28, 107.70, end)

Conductivity (sand lenses) = 0.1 cm/sec
Porosity =

Limestone aquifer zone:
(0, 0), (30.48, 15.24), (39.624, 30.48), (609.6, 30.48), (609.6, 0, end)

Conductivity (limestone aquifer) = 104 cm/sec
Porosity = 0.15

 
   
   
 

 
 

WATER .
0323*“ TABLE V355”
hw-550+o.04(750-x)
11,- 550 ft Waste hw - 570 a
_, ‘ __l?_ond_ __
_- ‘ ‘ ’ """""""" -_ ' ah_
SAND / F RACTU R E %— -0
#i 580 e CONSTANT
HEAD
h, - 552 ft
k 29.
az 0
w 2000 e .

 

 

Figure 10.18 Waste pond model with individual fractures

4. Internal Sources/Sinks: The plume region is considered an “instantaneous” source.
We used polylines for the constant head boundaries on the left and right sides of each
model. We also use polylines to construct the fractures.

Left polyline (along lake boundary) = 550 feet

Right polyline (in limestone aquifer) = 552 feet

Fracture polylines: polyline feature in IGW choose "Model polyline as a fracture" in the
Hydraulic Conductivity section in the IGW interface and use the fracture width of 0.01 m

222

5. Run Flow Model/Add Particles: The flow pattern can be obtained at this step. The
solution for this problem is a transient-state solution, since we observe the
plume/particles moving.

6. Contaminant Zones: The waste pond is the contaminant zone. It has a concentration
of 100 ppm and is a continuous source.

7. Profile Model: We did not define a profile model.

8. Transient or Steady-State: The solution was obtained using the “transient” mode.
9. Grid Size: We discretized the model with 100 x 53 cells, which gives a cell
dimension of 4.1 meters x 4.2 meters, approximately. When using particles whether for

visualization or for contamination, the model must be finely discretized.

10. Static Profiles: No static profiles were done for this model.

223

 

10.1.16 Probability Model

1. Domain: 7200 it x 5400 ft.
2. Zones and Parameter Values: Constructed the following zones (with coordinates)
and assign parameter values:

Porous Medium (IGW currently only allows units of meters to be entered):
In order to keep units of feet, we convert feet to meters for entry of coordinates.

Plume zone: (914.4, 457.2), (944.88, 444.95), (975.36, 365.76), (944.88, 286.57),
(914.4, 274.32), (883.92, 286.57), (853.44, 365.76), (883.92, 444.95)

These points are located on an ellipse with semimajor axis =300fi and semiminor
axis = 200 ft centered at (914.4, 365.76).

South river zone:
(0, 0), (0, 60.96), (2194.56, 60.96), (2194.56, 0, end)

 

North river zone:
(0, 1584.96), (0, 1645.92), (2194.56, 1645.92), (2194.56, 1584.96, end)

Lake: (807.41, 923.32), (745.30, 879.53), (716.07, 832.03), (712.42, 788.29),
(737.99, 759.09), (814.21, 737.20), (891.44, 722.60), (946.24, 733.55), (1008.35,
744.50), (1081.41, 773.69), (1110.64, 817.49), (1121.60, 853.98), (1096.03,
897.77), (1052.19, 934.27), (960.85, 956.17), (880.48, 948.87, end)

Parameter values:
Conductivity 100 ft/day (average)
Porosity = 0.3

3. Scatterpoint data

For this step we must create the data using scatterpoints. So, we sample randomly 150
different points in a somewhat arbitrary sub-region of the overall model. This region is
the predicted area that the plume will travel.

 

A. We first create the random field for the entire domain.

 

Random Conductivity = 100 fi/day

Random => Scale 1 Scale 2
Spectral Algorithm Spectral Algorithm
Ax=ly=10 Ax=ly=100
(Variance) 0'2= l (Variance) 0'2: 0
Exponential Model Anistropic Bell Model
Angle = 90 Angle = 90

224

Nugget = 0.01 Nugget = 0.01

The overall variance is then obtained by summing the variances for each scale. In this
case, we have a variance of 1.

B. Then select the polygon for the sub-region.

C. Then, in the file "Utilitiy", we select “Random Sampling” with 150 data points.

D. IGW then returns a table (ScatterPoints.csv file) with the data points. We then save
the data in the same file directory as the main model.

B. We then re- open the main model and input the ScatterPoints.csv file for the sub-
region zone. By right clicking on a scatterpoint we can select "Switch List" to enable
conditioning on the data points, Then we can select the newly formed group of
conditional points and obtain another menu, which includes the different interpolatory
schemes and exploratory data analysis. (Notice in the figure below that the random
field in the sub-region is created from the scatterpoints in the sub-region area.)

 

 

Figure 10.19 Probability model

4. Exploratory Data Analysis: At this point we can view the exploratory data analysis
of the scatterpoints as well as view the experimental variogram. The exploratory data
analysis gives a variety of statistics for the scatterpoint data such as the maximum,
minimum, mean, median, mode, variance, standard deviation, skewness, coefficient of

225

variance, the lower quantile 25%, and the upper quantile 75%. It also show the
histogram, PDF, CDF, and h-scatterplot of the data points. The variogram gives insight
into whether or not the data (such as conductivity data in this case) is nonstationary or
stationary.

Model-direction 1 (angle-90)
0 Experimental data —direction1

 

Vario gram

as .
(L .. . O O .
as

 

0.4”

 

0.2:

 

 

a=..‘..7..‘..=s.‘...'... no
Dhtmcefll)
Figure 10.20 Variogram

5. Internal Sources/Sinks: The plume region is considered an “instantaneous” source.
We have the following pumping rates (point source/sink) for the wells:

Town Well = -500 gpm

Rural Well 1 = -100 gpm

Rural Well 2 = -100 gpm

Rural Well 3 = -200 gpm

South River: Stage = 5 feet
Sediment conductivity = 0.1 m/day
Sediment thickness = River bot elev. - aquifer top elev.
Bottom Elevation = -30 feet

North River: Stage = 0 feet
Sediment conductivity = 0.1 m/day
Sediment thickness = River bot elev. — aquifer top elev.
Bottom Elevation = -30 feet

6. Run Flow Model/Add Particles: The flow pattern can be obtained at this step, but for
there are a couple of extra steps for performing a Monte Carlo simulation. The solution
for this problem is a transient-state solution, since we observe the plume/particles
moving.

A. For a Monte Carlo simulation, we first chose a “conditional simulation” and the
“spectral approach” in the menu for the conditional data points.

B. Then under the “Solver” menu, we chose a Monte Carlo simulation with 100
realizations.

226

We also placed a static profile to obtain the concentration over distance in each of the
realizations.
The time step was set to 40 days for duration of 4800 days (13.15 years). We viewed
the results at 4800 days for each realization.
Specify that the “Means and Variances” should be calculated.
Add monitoring wells to obtain the different post statistics and breakthrough curves.
We selected the options to "Monitor head and concentration" and to "Monitor
probability distribution".

Monitoring Well 1 is located at (727.031, 733.547).

Monitoring Well 2 is located at (836.634, 547.423).

.5 .0

7’1!“

6. Contaminant Zones: We used MMOC for the numerical method of tracking the
plume.

Plume concentration = 500 ppm (instantaneous source)
7. Profile Model: We did not define a profile model.
8. Transient or Steady-State: The solution was obtained using the “transient” mode.

9. Grid Size: We discretized the model with 200 x 154 cells, which gives a cell
dimension of 10.7 meters x 10.8 meters, approximately. When using particles whether
for visualization or for contamination, the model must be finely discretized.

10. Results and Static Profiles: We obtained the following various results.

A. We obtained the “Means and Variances” graph, which shows the probability of
the plume passing through a particular point.

B. We obtained the experimental probability density firnction (PDF), the cumulative
distribution function (CDF), and the process for the hydraulic head, conductivity,
and concentration.

C. For the concentration, we obtained the concentration breakthrough curves at the
different monitoring wells.

D. After the 100 realizations were finished we ran a couple more simulations to show
interesting static profiles graphs from model.

227

 

10.2 Appendix B

Matlab code with associated sinusoidal curves.

clear all;

x=0:500:20000;

m=size(x)

u=pi/100;

a=100;

A=a/cos(u);

A_=2*A;

b=22/20000

B=b/cos(u);

C=tan(u);

zO=1000;

i=1 :41;
f(i)=zO+C*x+A*sin(B*x);
g(i)=zO+C*x+A_*sin(B*x);
plot(12*2.54*x./ 100,12*2.54*f./ 100)
hold on
plot(12*2.54*x./100,12*2.54*g./100,'r')
axis tight;

T=zeros(4l ,2)
U=zeros(41,2)
T(:,1)=12*2.54*x'./ 100;
T(:,2)=12*2.54*f./ 100;
U(:,1)=12*2.54*x'./ 100;
U(:,2)=12*2.54*g'./ 100;

T

U

228

 

 

450

Meters

Figure 10.21 Undulations for Toth examples

’1":

 

 

./-

Undulations for Toth Examples

1

T

1

 

1.06r003 "'
0 0.3048

0.1524
0.3048
0.4572
0.6096
0.7620
0.9144
1.0668
1.2192
1.3716
1.5240
1.6764
1.8288
1.9812
2.1336
2.2860
2.4384
2.5908
2.7432
2.8956
3.0480
3.2004
3.3528
3.5052
3.6576
3.8100
3.9624
4.1148
4.2672

0.3255
0.3416
0.3496
0.3486
0.3403
0.3287
0.3184
0.3141
0.3183
0.3312
0.3505
0.3719
0.3904
0.4020
0.4047
0.3992
0.3884
0.3769
0.3696
0.3701
0.3795
0.3966
0.4177
0.4380
0.4528
0.4594
0.4570
0.4479

1000

2000

3000
Meters

4000

229

5000

 

4.4196
4.5720
4.7244
4.8768
5.0292
5.1816
5.3340
5.4864
5.6388
5.7912
5.9436
6.0960

0.4361
0.4266
0.4235
0.4292
0.4434
0.4634
0.4846
0.5022
0.5 124
0.5 137
0.5070
0.4958

1.0e+003 *
0 0.3048

0.1524
0.3048
0.4572
0.6096
0.7620
0.9144
1.0668
1.2192
1.3716
1.5240
1.6764
1.8288
1.9812
2.1336
2.2860
2.4384
2.5908
2.7432
2.8956
3.0480
3.2004
3.3528
3.5052
3.6576
3.8100
3.9624
4.1 148
4.2672
4.4196
4.5720
4.7244
4.8768
5.0292
5.1816
5.3340
5.4864
5.6388
5.7912
5.9436
6.0960

0.3415
0.3687
0.3800
0.3732
0.35 19
0.3238
0.2986
0.2850
0.2887
0.3098
0.3436
0.3815
0.4137
0.4322
0.4328
0.4169
0.3905
0.3628
0.3435
0.3396
0.3537
0.383 1
0.4204
0.4562
0.4812
0.4894
0.4800
0.4569
0.4286
0.4047
0.3938
0.4004
0.4240
0.4591
0.4967
0.5272
0.5429
0.5406
0.5224
0.4952

230

   

1(11111131111111)11111111