EXPERIMENTAL AND NUMERICAL ANALYSIS OF GROUNDWATER ?
SURFACEWATER INTERACTION IN HORIZONTAL AND
SLOPING UNCONFINED AQUIFERS
Except where reference is made to the work of others, the work described in this thesis is
my own or was done in collaboration with my advisory committee. This thesis does not
include proprietary or classified information.
__________________________
Gopal Chandra Saha
Certificate of Approval:
________________________ ________________________
Joel G. Melville T. Prabhakar Clement, Chair
Profesor Profesor
Civil Engineering Civil Engineering
________________________ ________________________
Xing Fang George T. Flowers
Associate Professor Dean
Civil Engineering Graduate School
EXPERIMENTAL AND NUMERICAL ANALYSIS OF GROUNDWATER ?
SURFACEWATER INTERACTION IN HORIZONTAL AND
SLOPING UNCONFINED AQUIFERS
Gopal Chandra Saha
A Thesis
submitted to
the Graduate Faculty of
Auburn University
in Partial fulfillment of the
requirements for the
degree of
Master of Science
Auburn, Alabama
August 10, 2009
iii
EXPERIMENTAL AND NUMERICAL ANALYSIS OF GROUNDWATER ?
SURFACEWATER INTERACTION IN HORIZONTAL AND
SLOPING UNCONFINED AQUIFERS
Gopal Chandra Saha
Permission is granted to Auburn University to make copies of this thesis at its discretion,
upon the request of individuals or institutions and at their expense. The author reserves
all publication rights.
Gopal Chandra Saha
Signature of Author
August 10, 2009
Date of Graduation
iv
VITA
Gopal Chandra Saha, son of Santosh Kumar Saha and Niva Rani Saha, was born
on July 1
st
, 1979. In March 2003, he graduated from Bangladesh University of
Engineering and Technology (BUET) with a Bachelor?s of Science degree in Civil
Engineering. After graduation, he worked in a construction company for one year. Then
he did his Master?s of Science degree in Water Resources Engineering and Management
from Stuttgart University, Germany in December 2006. In January 2008, he started
coursework to pursue his Master?s of Science degree in Civil Engineering in Auburn
University under the guidance of Prof. T. Prabhakar Clement. He married Pinki Saha,
daughter of Utkott Saha and Jyotti Saha, on December 6
th
, 2008.
v
THESIS ABSTRACT
EXPERIMENTAL AND NUMERICAL ANALYSIS OF GROUNDWATER ?
SURFACEWATER INTERACTION IN HORIZONTAL AND
SLOPING UNCONFINED AQUIFERS
Gopal Chandra Saha
Master of Science, August 10, 2009
(B.Sc., BUET, Bangladesh, 2003)
84 Typed Pages
Directed by T. Prabhakar Clement
We investigated the interaction between the water moving from an unconfined
aquifer and a surfacewater body. This interaction can occur during or after a rainfall
event, which can cause the groundwater table to change. We developed a laboratory-scale
aquifer model, which is bounded by a stream and a water divide. The model simulates
both horizontal and sloping aquifer conditions. We conducted both steady-state and
transient experiments. First, we set the water table at a steady state using a constant
recharge rate and then the aquifer was allowed to drain. The experimental data were
vi
recreated using a numerical model, which simulated the observed water table variations.
The transient water table data indicated highly non-linear changes. It was found that
under similar recharge conditions, water table in the horizontal unconfined aquifer is
higher than that in the sloping unconfined aquifer (2.03 degree slope). Also, the
horizontal unconfined aquifer required more time to drain. The numerical model was able
to predict the observed groundwater table data. The proposed numerical model is a useful
tool for managing groundwater-surfacewater management problems.
vii
ACKNOWLEDGEMENTS
First of all, I would like to thank my advisor Prof. T. Prabhakar Clement for his
enormous guidance, suggestion and inspiration in completing the thesis successfully. His
constructive criticism has helped me to organize the work. I would also like to thank
Professors Joel Melville and Xing Fang for serving as my committee members, and for
their support. Special thanks to Jagadish, Shyam, Sunwoo and Eric for their help in doing
lab work.
I would also like to show my deepest gratitude to my parents and my family
members for encouraging me during my whole degree program. This research was, in
part, funded by the Auburn University WRC-Wolf Bay Basin Project # 101002-107311-
700-2055.
viii
Style manual or journal used Auburn University Graduate School Guide to Preparation of
Master?s Thesis, Water Resources Research Journal.
Computer software used Microsoft Office 2003 (Microsoft Word, Microsoft Excel),
Visual Basic.
ix
TABLE OF CONTENTS
LIST OF FIGURES ........................................................................................................... xi
LIST OF TABLES........................................................................................................... xiv
CHAPTER 1. INTRODUCTION .......................................................................................1
1.1 Background ......................................................................................................1
1.2 Objectives .........................................................................................................2
1.3 Organization of thesis .......................................................................................3
CHAPTER 2. LITERATURE REVIEW .............................................................................4
2.1 Unconfined hillslope flow .................................................................................4
2.2 Derivation of the governing flow equation......................................................10
2.3 Review of previous works related to groundwater-surfacewater Interactions.12
2.3.1 Review of analytical solutions to the governing flow equation ..........13
2.3.2 Review of numerical solutions.............................................................19
2.3.3 Review of experimental studies...........................................................21
CHAPTER 3. EXPERIMENTAL METHODS .................................................................23
3.1 Column experiments ........................................................................................23
3.2 Darcy experiments ...........................................................................................27
x
CHAPTER 4. DEVELOPMENT AND TESTING OF A NUMERICAL MODEL .........30
4.1 Numerical solution strategy .............................................................................30
4.2 Initial and boundary conditions .......................................................................33
4.3 Code verification using published datasets......................................................34
4.3.1 Comparison against Kim et al (2001) horizontal aquifer solution......34
4.3.2 Comparison against Verhoest et al (2002) sloping aquifer solution...36
4.3.3 Comparison against Sanford et al (1993) solution..............................38
CHAPTER 5. LABORATORY DATA AND MODELING RESULTS...........................42
5.1 Steady state experiments...................................................................................42
5.2 Modeling steady-state water table profiles .......................................................48
5.3 Transient experiments.......................................................................................51
5.3.1 Modeling the transient head boundary at the left side of the tank.........51
5.3.2 Transient solution for horizontal unconfined aquifer ............................54
5.3.3 Transient solution for sloping unconfined aquifer.................................57
5.3.4 Accuracy of the numerical model..........................................................58
5.4 Correlation between the sloping unconfined aquifer with the horizontal
unconfined aquifer ............................................................................................59
CHAPTER 6. SUMMARY, CONCLUSIONS AND RECOMMENDATIONS ..............61
REFERENCES ..................................................................................................................64
APPENDIX........................................................................................................................68
xi
LIST OF FIGURES
Figure 1: The schematic presentation of sloping bed unconfined aquifer...........................5
Figure 2: An arbitrary point in the porous medium of sloping bed unconfined aquifer.....6
Figure 3: Elevation of the water table in sloping bed unconfined aquifer...........................8
Figure 4: Horizontal bed unconfined aquifer under constant recharge............................10
Figure 5: Sloping bed unconfined aquifer under constant recharge ..................................10
Figure 6: 2-D control volume in an unconfined aquifer ...................................................11
Figure 7: Column experiment setup in the laboratory .......................................................24
Figure 8: Darcy experiment setup in the laboratory ..........................................................28
Figure 9: Comparison Kim et al. (2001) analytical model (linearized) results against our
proposed model (non-linear) results for a horizontal unconfined aquifer
example problem...............................................................................................35
Figure 10: Comparison between analytical model results of Kim et al. (2001) and the
linearized numerical model for a horizontal unconfined aquifer example
problem ............................................................................................................36
Figure 11: Comparison Verhoest et al. (2002) analytical results against linear and non-
linear numerical models results for a sloping unconfined aquifer example
problem ............................................................................................................37
xii
Figure 12: Comparison Sanford et al. (1993) analytical model results against our
proposed model (non-linear) results for a horizontal unconfined aquifer
example problem...............................................................................................39
Figure 13: Comparison Sanford et al. (1993) analytical model results against our
proposed model (non-linear) results for a sloping unconfined aquifer
(theta = 0.089 degree)example problem ...........................................................40
Figure 14: Comparison Sanford et al. (1993) analytical model results against our
proposed model (non-linear) results for a sloping unconfined aquifer
(theta = 0.137 degree)example problem ...........................................................41
Figure 15: Laboratory tank with recharge generator used in this study for horizontal
unconfined aquifer ..........................................................................................44
Figure 16: Laboratory tank with recharge generator used in this study for sloping
unconfined aquifer (i = 2.03 degree)...............................................................44
Figure 17: A schematic laboratory set up for tank experiment..........................................45
Figure 18: Experimental Steady-state water table profiles in horizontal unconfined
aquifer under different recharge rates...............................................................46
Figure 19: Experimental Steady-state water table profiles in sloping unconfined aquifer
(i=2.03degree) under different recharge rates..................................................47
Figure 20: Modeling steady-state water table profiles using analytical equation in
horizontal unconfined aquifer...........................................................................49
Figure 21: Modeling steady-state water table profiles using numerical model in
sloping unconfined aquifer (i = 2.03 degree)...................................................50
xiii
Figure 22: Modeling transient head boundary (left hand side of the tank) in horizontal
unconfined aquifer ............................................................................................52
Figure 23: Modeling transient head boundary (left hand side of the tank) in sloping
unconfined aquifer (i = 2.03 degree).................................................................53
Figure 24: Transient behavior of water table in horizontal unconfined aquifer (k = 90
cm/min (1296 m/day),
y
S = 0.35 and w = 5.4 cm/min at the beginning (steady
state))................................................................................................................55
Figure 25: Transient behavior of water table in horizontal unconfined aquifer (k = 90
cm/min (1296 m/day),
y
S = 0.35 and w = 4.42 cm/min at the beginning
(steady state))...................................................................................................56
Figure 26: Transient behavior of water table in sloping unconfined aquifer (i = 2.03
degree, k = 90 cm/min (1296 m/day),
y
S = 0.35 and w = 4.42 cm/min at the
beginning (steady state))................................................................................57
Figure 27: Comparison of experimental steady state water table profile between
horizontal and sloping unconfined aquifers under same condition .................60
Figure 28: Pressure at the bottom of a sloping bed unconfined aquifer ...........................68
xiv
LIST OF TABLES
Table 1: Results of Column experiment ............................................................................24
Table 2: Darcy velocities and Reynolds numbers in column experiments........................26
Table 3: Results from Darcy experiments .........................................................................28
Table 4: Darcy velocities and Reynolds numbers in Darcy experiments ..........................29
Table 5: Parameters used in Sanford et al. (1993) analytical model .................................38
Table 6: Recharge rates and constant head boundary conditions for both horizontal and
sloping unconfined aquifers (i = 2.03 degree) .....................................................43
Table 7: Mass balance errors for the water table profiles generated by the numerical
model for the horizontal unconfined aquifer ........................................................59
Table 8: Mass balance errors for the water table profiles generated by the numerical
model for the sloping unconfined aquifer (i = 2.03 degree) .................................59
1
CHAPTER 1
INTRODUCTION
1.1 Background
Interaction between groundwater and surfacewater systems is a common
phenomenon observed in nature. This interaction can occur during or after a rainfall
event. When rain water infiltrates into the ground it raises the water table which will then
induce groundwater drainage to a surface water body. For sustainable groundwater
management, it is important to understand the dynamics of groundwater table movement
during and after recharge events.
Groundwater discharge contribution to stream flow is known as groundwater
runoff or base flow. Traditionally it is assumed that during precipitation, stream flow is
generated primarily from surface runoff; while during dry period, stream flow is
generated primarily from groundwater flow (Hall, 1968, Norum et al., 1968). However,
recent studies have shown that groundwater flow (so called old water) can dominate
stream flow even during a precipitation event (Shanley et al., 2002). Typically, base flow
2
is a good indicator of aquifer storage characteristics within a groundwater basin
(Meyboom, 1961).
In order to estimate base flow, a rating curve of base flow can be made by
plotting mean groundwater table within a basin against stream flow during dry periods
(Schneider, 1961). During or after a rainfall event in a small basin the groundwater table
could rise rapidly and contribute to the base flow. The variations in base flow, with time
during and after rainfall event within a basin is indicated by a recession curve (Chow et
al., 1988). Hence, recession curve is a good measure of drainage rate from a basin
(Werner et al., 1951).
Groundwater table response to a recharge event may depend on the topography of
the aquifer bottom. It is necessary to understand this response because sloping hillslopes
are the basic landscape elements in many catchments (Torch et al., 2003). In addition,
hillslope subsurface flow is one of the key elements that control the hydrology of upland
watershed under both wet and dry conditions (Brutsaert, 1994).
1.2 Objectives
We investigated the hydraulics of unconfined flow near a surface water body
during and after a rainfall event under both horizontal and sloping aquifer conditions. We
developed a laboratory-scale aquifer model to conduct multiple experiments. A numerical
model, which is based on a fully implicit finite difference scheme with Picard iteration,
was used to model the experimental results.
3
1.3 Organization of thesis
This thesis is organized into six chapters including this introduction chapter.
Chapter 2 provides the literature review and derivation of the governing equations for
both horizontal and sloping-bed unconfined aquifers. It also summarizes previous works
related to groundwater-surfacewater interactions. Chapter 3 describes all the experimental
methods employed in this study. Chapter 4 details the development of a numerical model
and the validation of the model results based on published analytical and numerical
results. Chapter 5 discusses the modeling of laboratory data, and the transient model
solution for both horizontal and sloping unconfined systems. Chapter 6 summarizes the
conclusion of this study and provides some recommendation for future research in this
area. An alternate analysis of depth-averaged Darcy flux in sloping bed unconfined
aquifer is listed in the appendix.
4
CHAPTER 2
LITERATURE REVIEW
2.1 Unconfined hillslope flow
Changes in the groundwater table in an unconfined aquifer after a recharge event
has been studied by many researchers. Boussinesq (1877) first developed the general
depth averaged flux equation for unconfined flow in a sloping aquifer based on depth
averaged Darcy?s law. Figure 1 shows the schematic presentation of a sloping bed
unconfined aquifer.
In Figure 1, x is along the bed,
x is the Cartesian coordinate, and i is the slope. The
relationship between x and
x is given by the following equation:
cos
x
i
x
= (1)
5
Figure 1: The schematic presentation of sloping bed unconfined aquifer.
In Figure 2, an arbitrary point (x, Z) is considered in the porous medium of sloping bed
unconfined aquifer. The elevation of point (x, Z) = z
cos sinzZ ix i=+ (2)
From Darcy?s law,
H
vK
x
?
=?
?
(3)
Where v is Darcy flux or Darcy velocity or specific discharge
L
T
? ?
? ?
? ?
, K is the hydraulic
conductivity
L
T
??
??
??
, and
H
x
?
?
is hydraulic gradient.
h
i
x=0
X=L
1
h
2
h
D
z
x
x
i
6
Figure 2: An arbitrary point in the porous medium of sloping bed unconfined
aquifer.
Total head is defined as:
2
2
pv
H z
g?
= ++ (4)
where
p
?
is pressure head,
2
2
v
g
is velocity head, and z is elevation head. Since the
groundwater velocity (v) is typically a very small, we can neglect the velocity term.
Hence,
p
Hz
?
=+ (5)
Darcy velocity in Z direction
z
H
vK
Z
?
=?
?
x
Z
z
Z
(x,Z)
i
7
z
P
vK z
Z ?
???
=? +
??
?
??
(6)
If the flow is parallel to the bed, 0
z
v = . After substituting 0
z
v = in (6) we get,
0
P
Kz
Z ?
???
?+=
??
?
??
(7)
After differentiating we get
1
0
Pz
ZZ?
??
+=
??
(8)
From Figure 2 (or differentiating (2)) we can obtain cos
z
i
Z
?
=
?
; and using this we get,
1
cos 0
P
i
Z?
?
+=
?
(9)
Darcy velocity in x direction,
x
P
vK z
x ?
???
=? +
??
?
??
1
x
P z
vK K
xx?
??
=? ?
??
(10)
From Figure 2,sin
z
i
x
?
=
?
; and using this we can get,
1
sin
x
P
vK Ki
x?
?
=? ?
?
(11)
8
x
Z
Z
i
h
h = h(x)
Figure 3: Elevation of the water table in sloping bed unconfined aquifer.
The elevation of the water table function (h = h(x)) in sloping bed unconfined aquifer
is shown in Figure 3. Using this figure we can develop the following boundary condition
for the water table:
at Z h= , 0P =
From equation (8)
cos
P
i
Z
?
?
=?
?
(12)
After integrating
cos ( )P iZ f x?=? + (13)
9
Using the water table boundary condition we get,
() cosf xih?=
After substituting f(x) in (13) we get,
cos cosP iZ ih? ?=? + (14)
Substituting P in equation (11) we get,
(cos cos) sin
x
K
viZihK
x
??
?
?
=? ? + ?
?
cos sin
x
Kh
viKi
x
?
?
?
=? ?
?
cos sin
x
h
vKi i
x
?
??
=? +
??
?
??
(15)
This also shows that ()
xx
vvZ?
Depth-averaged Darcy flux can be defined as:
x
qvh
?
= (16)
where q
?
is depth averaged flux along the x direction
2
L
T
? ?
? ?
? ?
. After substituting (15) in (16)
we get the general depth averaged flux equation for unconfined flow in a sloping aquifer
as:
()
cosqKhhiz
x
?
?
=? +
?
cos
hz
qKh i
x x
?
??
??
=? +
??
??
??
cos sin
h
qKh ii
x
?
?
??
=? +
??
?
??
(17)
10
2.2 Derivation of the governing flow equation
The Boussinesq equation for a sloping unconfined aquifer with a constant
recharge rate can be formulated by combining equation (17) with the continuity equation.
Schematic of the horizontal and sloping bed unconfined aquifers considered in this study
are shown in Figures 4 and 5.
h
w
i
X=0
X=L
1
h
2
h
h
w
2
h
1
h
X=0
X=L
Figure 4: Horizontal bed unconfined aquifer under constant recharge
Figure 5: Sloping bed unconfined aquifer under constant recharge
11
Now, let us consider a 2-D control volume for groundwater occurring in these systems, as
illustrated in Figure 6.
Figure 6: 2-D control volume in an unconfined aquifer
Applying mass balance over this control volume for a system with a constant
recharge rate (w), we get the following conservation equation:
xxx yyy y
qq ytqq xtwxytSxyh?? ?? ?
?? ??
+? +?
???????+??+??=??
????
(18)
where
y
S is the specific yield and w is the recharge rate
L
T
? ?
? ?
? ?
. After dividing by x yt?? ??
and considering 0, 0, 0x y and t?? ?? ?? we get the following 2-D depth-averaged
flow equation,
x?
x
q?
?
xx
q?
?
+?
y
q?
?
yy
q?
?
+?
y?
12
y
hqq
Sw
txy
??
???
=? ? +
???
(19)
For one dimensional flow, the above equation reduces to,
y
hq
Sw
tx
?
??
=? +
??
(20)
Inserting (17) in (20) and considering homogeneous and isotropic values of
hydraulic conductivity and specific yield values, we get the following one-dimensional
Boussinesq equation for a sloping unconfined aquifer:
cos sin
y
hhh
SKih w
txxx
???? ??
??
=++
????
??
??
??
(21)
For horizontal (i = 0) unconfined aquifer, the above equation reduces to,
y
hh
SKh w
txx
?????
??
=+
????
???
??
??
(22)
2.3 Review of previous work related to groundwater-surfacewater interactions
An approximate analysis of groundwater flow resting on an impermeable bed was
initiated by Dupuit (1863). The common assumptions used in that study were: 1) the
stream or drain boundaries completely cut through the entire depth of the aquifer and 2)
groundwater flow occurs in a homogeneous, isotropic formation with hydraulic
properties that remain constant. Furthermore, these analyses assume the following
13
assumptions known as the Dupuit-Forchheimer assumptions (as summarized in Reddi,
2003):
a) For small changes in the slope of line of seepage, the hydraulic head is independent of
depth.
b) The hydraulic gradient causing flow is equal to the slope of the water table.
2.3.1 Review of analytical solutions to the governing flow equation
Since equation (21) is a non-linear equation, it does not have a general analytical
solution. Therefore, the equation is often simplified and solved by introducing additional
approximations. The most common approximation used is linearization. Boussinesq
(1877) suggested a simple linearization method that can be used when the water table
variation is small compared to the aquifer depth.
There are several analytical solutions for horizontal and sloping unconfined flow.
Some of them assume constant recharge rate and others assume time-varying recharge
rates. Since we used a constant recharge rate in our study, here we only review analytical
solutions that considered a constant recharge rate. Marino (1974) derived the following
analytical equation for transient water profiles in a horizontal aquifer.
22 2 2
0
00
2()(1)
(,) 1 (1)4 (1)4
22
nn
wDt l n x l n x
hxt h i erfc i erfc
S
KDt KDt
SS
??
==
??
??
??
??
++?
?= ? ? ? +? ?
??
??
??
??
(23)
14
He used the following initial and boundary conditions to develop the above
analytical equation:
Initial condition:
22
0
0hh?= at t = 0 and 0x?
Boundary conditions:
22
0
0hh?= at x= 0 and 0t >
22
0
()0hh
x
?
? =
?
at
2
l
x= and 0t >
where
0
h is the initial depth of saturation of the aquifer, h is the height of the water table
above the base of the aquifer after the incidence of recharge, l is the length of the aquifer
which is bounded by two fixed head boundaries,
2
()ierfcy is the second repeated integral
of the error function of argument y, h(x, t) is the hydraulic head at a particular time and
space, K is the hydraulic conductivity and D is the aquifer depth, w is the recharge rate,
and S is the specific yield.
Melville et al. (1987) derived steady-state dimensionless analytical equations
for different scenarios of hillslope flow problem. They used the following boundary
conditions to get those analytical equations:
Boundary conditions: 0
dh
dx
= (no flow) at x= 0 and 0t >
1
hh= at x= L and 0t >
When F = 0.25,
1
2
2
0
ln ln 1
2
2
X
HXX
H
X
HHH
H
??
??
??? ?
=??+
????
??? ?
??? ?
??
???
(24)
15
where
2
R
F
KT
= , R is vertical infiltration or recharge rate, K is hydraulic conductivity,
tanT ?= , ? is bed slope,
x
X
L
= , x is the coordinate parallel to the bed, L is the length
of aquifer,
2
hC
H
LS
= , cosC ?= , sinS ?= and h is the vertical depth of the water table.
When 0.25F <
()
()
1
2
2
0
2114
1
ln ln ln 1
21 4
2114
X
F
HXXH
F
XHHH
F
F
H
????
???
??
??
? ?
??
?? ??
??
=??+
? ???
?? ?? ??
???
?? ??
? ???
? ?
?+?
??
????
(25)
When 0.25F >
1
2
2
1
0
2
1
ln tan ln 1
2
41
41
X
HXXH
F
XHHH
F
F
H
?
?
??
????
?
????
??
? ?
??
?? ?? ??
??
=? ??+
? ???
?? ?? ?? ??
???
?? ?? ??
? ???
??? ?
?
??
??
??
??
(26)
Sanford et al. (1993) derived analytical equation for transient water profiles in
both horizontal and sloping aquifer. They considered a system with the following initial
and boundary conditions.
Initial condition: h = h (0, 0) at t = 0 and 0x?
Boundary conditions: h = h (0, t) at x = 0 and 0t >
0
dh
dx
= (no flow) at x = L and 0t >
16
For horizontal aquifer, the solution is:
()
2
2(0,)
2
hLxxht
?
=?+
(27)
where L is the length of the aquifer, h (0, t) is water table at the boundary (where x=0),
and
3
3
(0, )
i
VV
htL
Lsw
?
???
=?
??
??
where
i
V is the initial volume of water contained in the rectangular region, V is the
cumulative outflow volume at any time, s is the specific yield, and w is the width of the
system.
For sloping aquifer, the transient water profile can be predicted by the following
analytical equation:
()
2
2tan(0,)
2
hLxxx ht
?
?=??+
(28)
where ? is slope of the impermeable layer, and
2
3
31
(0, ) tan
2
i
VV
htL L
Lsw
? ?
???
=?+
??
??
Fetter (1994) derived the steady-state analytical equation for a horizontal bed
unconfined aquifer using the following boundary conditions:
Boundary conditions:
1
hh= at x= 0 and 0t >
2
hh= at x= L and 0t >
17
The steady-state analytical solution is:
(){}
(){ }
22
12
22
1
hhx
wL xx
hh
LK
?
?
=? +
(29)
where h is the hydraulic head in the unconfined aquifer, w is the recharge rate, L is the
length of the aquifer, K is the hydraulic conductivity,
1
h and
2
h are left and right hand
side boundary heads, respectively. The above equation was used in this study to analyze
steady-state data.
Kim et al. (2001) derived an analytical equation for transient water profiles in
horizontal unconfined aquifer using the following initial and boundary conditions:
Initial condition:
0
()hhx= at t = 0 and 0x?
Boundary conditions:
1
hh= at x= 0 and 0t >
2
hh= at x L= and 0t >
The derived analytical solution becomes as:
2
2
()
21
1
1
(,) sin
22
n
t
L
n
n
hhx Lnx
hxt x h Ce
LL
?
?
? ??
??
?
?
=
?
??
=? + + + +
??
??
?
(30)
where h(x, t) is the hydraulic head at a particular time and space, L is the length of the
aquifer,
1
h and
2
h is left and right hand side boundary heads, respectively.
w
S
? = where w is the recharge rate, and S is the specific yield.
KD
S
? = where K is the hydraulic conductivity and D is the aquifer depth.
18
2
21
01
0
2
sin
22
L
n
hhx Lnx
Chh x d
LLL
?? ?
??
???
??
=?+?+ ?
????
??
??
?
Verhoest et al. (2002) also developed an analytical equation for computing the
transient state profiles of both horizontal and sloping aquifers. They used the following
initial and boundary conditions to derive those analytical solutions.
Initial condition: h = D at t = 0 and 0x?
Boundary conditions:
1
D
F
y
= at x = 0 and 0t >
2
D
F
y
= at x = L and 0t >
where D is the depth of the soil layer, F is the Laplace transform of the water table depth,
y is the Laplace variable,
11
DDH= ? , and
22
DDH= ? ,
1
H is the water table at x =0,
2
H is the water table at x =L, L is the length of the aquifer.
For the sloping aquifer case, the analytical equation is:
( )
( )
()( )
()
()()
12
1
22
12
2
22 2 2 2
1
2
22
2
2 2
22 2 2
1
2 sinh( )
2 sinh( ) 2
21
sin exp
21
sin exp
1
aL x
n
aL
ax
n
n
ax aL
n
n
eLNafKHH ax
Nx
hH
afK aL afK
nDH DHe
nx n
eaKt
aL n L L
nLNe e
nx n
aKt
LL
fK a L n
?
??
?
?
??
?
??
?
?
=
?
?
=
+???
??
=? + +
??? ?
? ?
??
?+
? ???
+
??
? ?
??
??
?
??
??
+?+
??
??+?
?
?
?
??
??
(31)
19
where N is the recharge rate, f is the specific yield of the formation,
2
U
a
k
=? ,
an
sinki
U
f
= ,
coskpD i
K
f
= , k is the hydraulic conductivity, p is the linearization
constant, and i is the slope of the aquifer bed.
For horizontal unconfined aquifer i=0, U = 0, and a = 0 and equation (31) can be
simplified as:
()()
()( )
()
12
21
1
1
2
22 22
233 2
1
21
2
211
sin exp sin exp
n
n
n
n
DH DH
xH H NxL x
hH
LfK
LN
nx n nx n
KtKt
LL fKn LL
?
?? ??
?
?
=
?
=
??
??? ?
??
??
=+ + +
??
??
?? ??
??
?? + ?
?? ??
?? ??
?
?
(32)
2.3.2 Review of numerical solutions
Marino (1975) presented a digital computer model that simulated the response of
an unconfined aquifer to changes in stream levels. The model is based on the one-
dimensional Boussinesq equation for horizontal unconfined aquifer and was solved using
a predictor-corrector scheme that used non uniform grid spacing. The numerical scheme
was found to be unconditionally stable.
20
Beven (1981) solved non-dimensional form of one dimensional Boussinesq
equation for a sloping unconfined aquifer with constant recharge. The solution method
used was the implicit finite difference approximation.
Sims (1986) presented another approach for solving hillslope flow problem. He
developed a numerical model, which was based on explicit finite difference scheme to
solve the nonlinear problem. The major disadvantage of that model was that
discretization step must be small enough to generate a stable solution.
Serrano (1998) presented a numerical model for modeling transient stream/aquifer
interactions in an alluvial valley aquifer. The model is based on the one-dimensional
Boussinesq equation for horizontal unconfined aquifer which was solved using a
decomposition method.
Parlange et al. (2001) developed a numerical model based on the one-dimensional
Boussinesq equation for horizontal unconfined aquifer. The equation was solved using a
finite element program written in PDE2D, a general purpose partial differential equation
solver. PDE2D uses a second-order accurate backward Euler scheme and an adaptive
time stepping scheme to generate transient solutions.
Verhoest et al. (2002) presented a numerical model and compared the results
against a transient analytical solution. The numerical scheme used the Crank-Nicholson
approximation and a finite element mesh to solve the one-dimensional linearized form of
Boussinesq equation. The finite element scheme employed a piecewise-linear Lagrange
basis function and a piecewise-uniform weighting function.
21
Rocha et al. (2007) developed a numerical model that considered two-
dimensional linearized Boussinesq equation for a horizontal aquifer receiving a constant
recharge. The equation was solved using the MODFLOW employing the block-centered
finite difference approach.
2.3.3 Review of experimental studies
Vaculin et al. (1979) conducted laboratory experiment with recharge in a three-
dimensional soil slab (6 m ? 5 cm ? 2 m). The slab was packed with fine rive sand
between two perspex walls supported by a frame resting on impervious horizontal
boundary. Both the sides of the slab were connected to constant head reservoirs. At the
soil surface, a constant flux q = 0.148 m/hr was applied over a width of 1 m in the center.
The remaining soil surface was covered to prevent evaporation losses. During the
infiltration, water content was measured using gamma rays, and water pressure heads at
different points were measured by 20 tensiometers, each connected to its own transducer
mounted on one of the perspex walls. Another transducer measured outflow volume from
the constant head reservoirs.
Sanford et al. (1998) conducted a series of experiments in a three-dimensional
glass-sided flume (245 cm ? 30 cm ? 30 cm). The flume was filled with a sand-gravel
mixture to a depth of 23 cm. Several small monitoring wells were placed in the inside
wall of the flume at various locations to measure the water table heights. One side of the
flume was forced with a constant head boundary (x = 0), and other side was a no flow
22
boundary. Drainage holes were installed on the downslope end of the flume (x = 0) at
various heights above the impermeable layer. At the beginning of the experiment, the
flume was filled with water to a saturated depth, designated as h (0, 0), when the flume
was in the horizontal position. Then the upslope end was raised to the desired height and
water was added at the upslope end to maintain a constant water depth as close to h (0, 0)
as long as possible until steady-state condition was reached. At steady state, water table
heights were recorded in various wells. Then the drainage holes at h = 0 (x = 0) were
opened and water addition was stopped to produce a sudden drawdown. The water table
height in each well and the cumulative outflow were recorded periodically until the
drainage stopped. The same procedure was applied for several slopes.
Kim et al. (2001) performed laboratory experiments using a three-dimensional
physical unconfined aquifer model (67 cm ? 15 cm ? 50 cm) that used acrylic plates to
visualize the water table levels from outside. Seven manometers were placed at a regular
spacing in the front wall to measure water table heights. The model consisted of an
unconfined aquifer region bounded by two reservoirs. The unconfined aquifer region was
packed with uniform fine sand materials. Two different flow experiments were completed
using two different recharge conditions. In the rising head experiment, the water table
was allowed to rise by applying a constant recharge. In the falling head experiment, the
water table was allowed to drain until a steady-state condition was reached. Both the
rising and falling head experiments were performed for two different types of boundaries,
one with equal heads, and the other with unequal heads. Transient water table heights
were recorded at each manometer at a regular time intervals for all cases.
23
CHAPTER 3
EXPERIMENTAL METHODS
3.1 Column experiments
A schematic column experimental setup is shown in Figure 7. The column was
3-cm diameter and 60-cm long. The column was filled with uniform glass beads (mean
diameter 1.1 mm) up to 36.9 cm from the bottom of the column. A constant head
reservoir was used to provide constant flow rate for the column experiment. Flow rate
and water levels at the column inlet and at exit boundaries were measured under steady-
state conditions. Three column experiments with three different flow rates were
conducted. Table 1 provides a summary of the experimental results.
24
Figure 7: Column experiment set up in the laboratory
Experiment
no.
1
h (cm)
2
h (cm) Q
(
3
/mincm )
Calculated K
(cm/min)
Average K
(cm/min)
1 49.1 44.3 88 95.68
2 47.9 44.3 69 99.36 96.07
3 55.5 44.3 200 93.19
Table 1: Results from column experiments
h1
h2
Q
Holding
stand
Porous
medium
Screen
25
The hydraulic conductivity values shown in Table 1 are computed based on the
following analysis:
According to continuity equation we know,
QAq= (33)
where Q= flow rate
3
L
T
??
??
??
, A = area of the medium
2
L? ?
? ?
, and q = Darcy flux
L
T
??
??
??
.
From Darcy?s law, we can write
h
qK
x
?
=?
?
(34)
Equation (34) can be written as
12
hh
qK
x
?
=
?
(35)
where K is the hydraulic conductivity, x? is the thickness of porous medium layer,
1
h
and
2
h are head levels in the column at the top and bottom boundaries, respectively.
By substituting (35) in (33) we get,
2
12
4
hhd
QK
x
? ?
=
?
2
12
4 x
KQ
dh h?
?
=
?
(36)
where d is the diameter of the column
[ ]
L .
Hydraulic conductivity of glass beads can be calculated using the Hazen method
(Hazen 1911). This method uses the following empirical equation,
2
10
()KCd= (37)
26
where K is hydraulic conductivity (cm/sec),
10
d is the effective grain size (cm), and C is a
coefficient based on particle size and packing of medium. For well sorted clean coarse
sand, C is 120 to 150. Considering this coefficient for glass beads, we get hydraulic
conductivity 1.45 cm/sec (87.12 cm/min) to 1.815 cm/sec (109.9 cm/min).
Darcy velocity and Reynolds number (R) for different flow rates in column
experiments were calculated. Darcy velocity was calculated using equation (35).
Reynolds number was calculated using the following equation (Bear 1978),
qd
R
?
= (38)
where q = Darcy flux
L
T
??
??
??
, d = diameter of the porous medium
[ ]
L , ? = kinematic
viscosity of fluid
2
L
T
??
??
??
. For room temperature, kinematic viscosity of water is
2
6
10
sec
m
?
.
The calculated values of Darcy velocities and Reynolds numbers in column experiments
are tabulated in Table 2.
Flow Rate, Q
(
3
/mincm )
Darcy flux
(cm/min)
Reynolds
Number
88 12.45 2.28
69 9.76 1.79
200 28.29 5.19
Table 2: Darcy velocities and Reynolds numbers in column experiments
27
3.2 Darcy experiments
Several Darcy experiments were done in a three-dimensional flow tank (115 cm
? 2.5 cm ? 60 cm), which was filled by uniform glass beads (mean diameter 1.1 mm). A
constant head reservoir was used to provide constant flow rate at the left hand side of the
tank for the Darcy experiment. Water levels at the left and right hand side of the tank, and
outflow from the right hand side of tank were measured under steady-state conditions.
Three Darcy experiments with three different flow rates were conducted. A schematic
Darcy experimental setup is shown in Figure 8. Table 3 provides a summary of the
experimental results. Table 4 provides the calculated values of Darcy velocities and
Reynolds numbers in Darcy experiments.
28
Figure 8: Darcy experiment set up in the laboratory
Experiment
no.
1
h (cm)
2
h (cm) Q
(
3
/mincm )
Calculated K
(cm/min)
Average K
(cm/min)
1 23.7 23.2 107 85.6
2 24.4 23.5 200 88.9 90
3 25.6 23.8 422 93.8
Table 3: Results from Darcy experiments
115 cm
60 cm
2.5 cm
Glass tube
Glass beads
Wire mesh
29
Flow Rate, Q
(
3
/mincm )
Darcy flux
(cm/min)
Reynolds
Number
107 1.98 0.363
200 3.63 0.66
422 7.43 1.36
Table 4: Darcy velocities and Reynolds numbers in Darcy experiments
From Tables 2 and 4, it is seen that Reynolds numbers are less than 10.
Experimental results have found that Darcy?s law is valid when the Reynolds number is
less than 1 to 10 (Lindquist 1933. Rose 1945, Bear 1978).
30
CHAPTER 4
DEVELOPMENT AND TESTING OF A NUMERICAL MODEL
4.1 Numerical solution strategy
A fully implicit finite difference scheme with a Picard iteration scheme for the
non-linear terms was used to solve the following sloping unconfined flow equation.
cos sin
y
hhh
SKih w
txxx
???? ??
??
=++
????
??
??
??
(39)
The nonlinear spatial terms in the above equation (right hand side) can be written as:
() ()
()
11
2
11
cos sin
cos
22
sin
2
mm mm
ii ii
ii ii
ii
hh
Kih i w
xx x
hh hh
Ki hh hh
x
Kihh
w
x
+?
+?
??? ??
??
++=
??
??
?? ?
??
??
??
??
?? ??
??
?? ?
??
?? ??
?? ??
?? ??
?
??
???
++
??
?
??
(40)
31
In the above equation, m is the Picard iteration level, x? is grid size,
i
h and
1i
h
+
are
the water level at the current (m + 1 )-th Picard iteration level at the i-th node and (i +1)-
th node, respectively. The first term in the equation is approximated using the forward
difference approximation, and water levels (h) are averaged using the arithmetic mean of
Picard updated head values of the adjacent nodes. The second term is approximated using
the central difference approximation.
By using
()
1
1
2
mm
ii
avg
hh
h
+
+
=
( )
1
2
2
mm
ii
avg
hh
h
?
+
=
We can write equation (40) as:
()(){}
()
11 2 1
11
2
cos sin
cos
sin
2
avg i i avg i i
ii
hh
Kih i w
xx x
Kihhhhhh
Kihh
w
xx
+?
+?
??? ??
??
++=
??
??
?? ?
??
??
??
?? ?
???
+ +??
??
??
??
??
(41)
The temporal term (left hand side term) in the equation (39) can be expressed
using backward-Euler approximation as:
()
ii
yy
hph
h
SS
tt
?
?
=
??
(42)
32
where
i
ph is the water level at the previous time level, and t? is the time step. Using
finite-difference expressions for the spatial and temporal terms of equation (39), we get
the following equation:
()
()(){ }
()
11 2 1
2
11
cos
sin
2
avg i i avg i i
ii
y
ii
Kihhhhhh
hph
S
tx
Kihh
w
x
+?
+?
??
?? ?
?
=? ?
??
??
??
???
++
??
?
??
(43)
After rearranging all the terms, we get the following equation:
( ) ( ) ( )
2112 11avg i avg avg i avg i i
hhhhhhhph? ????
?+
?+ + ++ +?? = + (44)
where
(tan)
2
x i
?
?
=
,
2
()
(cos)
y
Sx
Kt i
?
?
=
?
2
()
(cos)
wx
Ki
?
?
=
Equation (44) can be written as:
11iii
ah bh ch d
?+
++ = (45)
where,
2
,
avg
ah ?=? +
12avg avg
bh h ?= ++
1
,
avg
ch ?=? ?
i
dph? ?= +
33
In equation (45), all the unknowns are in the left hand side and all the known are
in the right hand side. Using this equation, the approximate water level for the new time
level at (m +1)-th Picard level is obtained after solving the linear matrix problem. The
matrix equation (45) can be written in the following form:
[]{} {}A hd=
(46)
where [ ]A is the square matrix which contains all the coefficients (except d), { }h is the
unknown water levels, and { }d is based on the known water levels at the previous time
level. At every Picard step, a, b, c vales are updated. Once the convergence is established
the computation is moved to the next time level by updating the { }d vector.
4.2 Initial and boundary conditions
For initial conditions, water tables before applying to recharge or the steady-state
water table condition after applying a constant recharge rate were used. For boundary
conditions, constant head, known transient head, and no flow conditions were used.
Equation (45) was appropriately modified to implement these boundary conditions. For
no flow condition, cos sin 0
h
Kh i i
x
?
??
?+=
??
?
??
was used at the right hand side boundary in
both horizontal and sloping bed unconfined aquifers, since cos sin
h
qKh ii
x
?
?
??
=? +
??
?
??
.
34
4.3 Code Verification using published datasets
4.3.1 Comparison against Kim et al. (2001) horizontal aquifer solution
Before using our proposed numerical algorithm in our study, we verified the code
against other numerical and analytical results that are published in journal articles. First,
we compared our numerical results to the analytical results of Kim et al. (2001). The
parameters used by Kim et al. (2001) model are
1
h = 14.5 cm,
2
h = 14.6 cm, theta or i= 0
degree (horizontal bed), aquifer depth (D) = 16 cm, hydraulic conductivity (K) = 6.41
cm/min, specific yield (
y
S ) = 0.33, recharge (w) = 1.96 cm/min, ?x = 1 cm, ?t = 0.0001
min, length of aquifer (L) = 47 cm, and total simulation time to reach the steady state was
3.3 minutes. It is important to note that the Kim et al. model used a liner approximation to
simplify the model. They assumed the aquifer depth was large enough compared to the
water table variations and substituted the nonlinear water table term (h) by aquifer depth
(D) in one-dimensional Boussinesq equation for horizontal unconfined aquifer.
We found that the numerical results matched at the beginning, but with the
passage of time some deviations were observed. At the steady state, the peak water table
predicted by Kim et al. (2001) is about 0.29 cm higher than that our proposed non-linear
model (Figure 9). But when we modified our nonlinear code by using aquifer depth (D)
as the nonlinear water table term (h) in one-dimensional Boussinesq equation for
horizontal aquifer to implement the Kim et al. (2001) linear approximation, the numerical
results matched well with the analytical results with very little deviation (see Figure 10).
35
10
12
14
16
18
20
22
0 5 10 15 20 25 30 35 40 45 50
Distance (cm)
W
a
t
e
r
Ta
bl
e
(
c
m
)
t=0.3 min(model) t=0.3 min (kim)
t=0.8 min (model) t = 0.8 min (kim)
t=1.8 min (model) t = 1.8 min (kim)
t=3.3 min (model) t = 3.3 min (kim)
Figure 9: Comparison Kim et al. (2001) analytical model (linearized) results against
our proposed model (non-linear) results for a horizontal unconfined
aquifer example problem
36
8
10
12
14
16
18
20
22
0 5 10 15 20 25 30 35 40 45 50
Distance (cm)
W
a
t
e
r
Tabl
e
(
c
m
)
t=0.3 min(mod) t=0.3 min(ana)
t=0.8 min (mod) t= 0.8 (ana)
t=1.8 min(mod) t=1.8 min(ana)
t=3.3 min(mod) t=3.3 min(ana)
Figure 10: Comparison between analytical model results of Kim et al. (2001) and the
linearized numerical model for a horizontal unconfined aquifer example
problem
4.3.2 Comparison against Verhoest et al. (2002) sloping aquifer solution
We compared our numerical model results against the results of Verhoest et al.
(2002). The parameters used in Verhoest et al. (2002) model are
1
h = 0.5 m,
2
h = 1.5 m,
theta or i= 2 degree, aquifer depth (D) = 2 m, hydraulic conductivity (K) = 0.001 m/s,
37
specific yield (
y
S ) = 0.34, recharge (w) = 3 mm/hr, and the total simulation time to reach
the steady state was 200 days. Figure 11 compares Verhoest et al. (2002) analytical
results, which is based on linearization, against the linearized and non-linearized
numerical results. From the figure, it can be observed that Verhoest et al. (2002) model
yielded higher peak (about 0.6 m more) than the non-linear numerical model under
similar conditions.
0
0.5
1
1.5
2
2.5
3
0 10203040506070809010
Distance (m)
Water Tabl
e
(m)
t=200 days(Non-linear Num.
model)
t=200 days(Linear num model)
t=200 days (Verhoest et al(2002)
analytical)
Figure 11: Comparison Verhoest et al. (2002) analytical results against linear and
non-linear numerical models results for a sloping unconfined aquifer
example problem
38
But, when we modified our numerical code to implement the linearization
procedure, the numerical model results matched well with Verhoest et al. results (2002)
(see Figure 11).
4.3.3 Comparison against Sanford et al. (1993) solution
We also compared our numerical model results against the analytical results of
Sanford et al. (1993) for three different cases. In Sanford et al. (1993) the right hand side
was set to no-flow boundary condition and the left hand side was set to transient head
boundary condition. The parameters used in Sanford et al. (1993) analytical model are
shown in Table 5, and the comparison of Sanford et al. (1993) analytical results and our
proposed model results are shown in Figures 12, 13, and 14. It is to be noted here that the
analytical model of Sanford et al. (1993) was developed on equations 22 and 23, and the
left side boundary of the experimental tank was suddenly dropped removing the cork
from the bottom hole (h = 0) only.
Case Slope
(degree)
h(0,0)
cm
K (cm/s)
y
S
1 0 22.3 0.76 0.238
2 0.089 22.2 0.78 0.238
3 0.137 21.8 0.64 0.238
Table 5: Parameters used in Sanford et al. (1993) analytical model
39
0
5
10
15
20
25
0 50 100 150 200
Distance (cm)
H
e
i
ght
(
c
m
)
IC(O min)
Sanford (1.6 min)
proposed (1.6 min)
sanford(5.5 min)
proposed (5.5 min)
sanford(12 min)
proposed(12 min)
sanford (33 min)
proposed(33 min)
Figure 12: Comparison Sanford et al. (1993) analytical model results against our
proposed model (non-linear) results for a horizontal unconfined aquifer
example problem
40
0
5
10
15
20
25
0 50 100 150 200
Distance (cm)
H
e
i
ght
(
c
m
)
IC(O min)
Sanford (0.67 min)
proposed (0.67 min)
sanford(1.5 min)
proposed (1.5 min)
sanford(3 min)
proposed(3 min)
sanford (5 min)
proposed(5 min)
Figure 13: Comparison Sanford et al. (1993) analytical model results against our
proposed model (non-linear) results for a sloping unconfined aquifer
(theta = 0.089) example problem
41
0
5
10
15
20
25
0 50 100 150 200
Distance (cm)
Heig
ht (cm)
IC(O min)
Sanford (0.33 min)
proposed (0.33
min)
sanford(2 min)
proposed (2 min)
sanford(3 min)
proposed(3 min)
sanford (5 min)
proposed(5 min)
Figure 14: Comparison Sanford et al. (1993) analytical model results against our
proposed model (non-linear) results for a sloping unconfined aquifer
(theta = 0.137) example problem
From the Figure 12, we found Sanford et al. (1993) analytical results in
horizontal aquifer matched well against our proposed numerical model results. For the
sloping aquifer, the results matched quite well, except near the no-flow boundary (see
Figures 13 and 14).
42
CHAPTER 5
LABORATORY DATA AND MODELING RESULTS
5.1 Steady state experiments
A three-dimensional flow tank (115 cm ? 2.5 cm ? 60 cm) which was connected
to a constant head reservoir and a constant recharge generator was used in this study. The
flow tank was made of plexi glass plates hence the water level can be recorded from
outside the tank. Uniform glass beads (mean diameter 1.1 mm) were used as the porous
medium to construct an unconfined aquifer. Multiple 6 mm diameter glass tubes filled
with red-colored water were used to delineate the location of the water table in the
aquifer. Digital photos of the actual tank fitted with recharge generator in both horizontal
and sloping bed unconfined aquifers are shown in Figures 15 and 16 respectively. A
schematic of the laboratory setup is shown in Figure 17.
Experimental works conducted with water flow from the constant head reservoir
to laboratory tank through a recharge generator. In all the steady-state experiments,
constant head boundary condition was used for the left side of the tank, and no flow
43
boundary condition was used for the right side of the tank. During infiltration, water table
in unconfined aquifer was raised and when it became steady-state, the outflow from the
right hand side of the tank was collected using a beaker for 2 minutes. The flow was
measured by dividing the collected outflow volume by 2 minutes. The sloping aquifer
was made by lifting the left hand side of the tank placing a wooden block of 3.8 cm
height below the bottom. Observations from multiple steady-state experiments conducted
with different recharge rates under horizontal and sloping aquifers are shown in Figures
18 and 19 respectively. Recharge rates were designed in such a way that there was no
ponding. The recharge rates and constant head boundary conditions for both horizontal
and sloping aquifers are tabulated in Table 6.
Table 6: Recharge rates and constant head boundary conditions for both horizontal
and sloping unconfined aquifers
Slope (degree) Constant
Boundary
condition (cm)
Recharge
rate
(cm/min)
0 22.0 2.9
0 22.75 4.2
0 25.0 5.4
2.03 23.8 4.42
2.03 24.3 4.83
2.03 26.0 5.55
44
Figure 15: Laboratory tank with recharge generator used in this study
for horizontal unconfined aquifer
Figure 16: Laboratory tank with recharge generator used in this
study for sloping unconfined aquifer (i =2.03 degree)
45
Figure 17: A schematic laboratory set up for tank experiment
115 cm
60 cm
2.5 cm
107cm
No flow
boundary
Glass tubes
Glass beads
Wire mesh
Constant head
boundary
Recharge
46
0
5
10
15
20
25
30
35
40
2040608010120
Distance (cm)
W
a
t
e
r
Tabl
e Hei
ght
,
h (
c
m
)
exp(w=2.9 cm/min)
expt(w=4.2 cm/min)
expt(w=5.4 cm/min)
No Flow
Boundary
Figure 18: Experimental steady-state water table profiles in horizontal unconfined
aquifer under different recharge rates.
47
0
5
10
15
20
25
30
35
40
0 20 40 60 80 100 120
Distance (cm)
Water Tabl
e Height, h (cm)
expt (w=5.55 cm/min)
expt(w=4.83 cm/min)
expt (w=4.42 cm/min)
No flow
Boundary
Figure 19: Experimental steady-state water table profiles in sloping unconfined
aquifer (i = 2.03 degree) under different recharge rates.
48
5.2 Modeling steady-state water table profiles
Steady-state water table profiles observed were modeled using analytical and
numerical models. For the horizontal aquifer data, we used the analytical equation (29),
which was derived by Fetter (1994). The hydraulic conductivity estimated from Darcy
experiments was used in analytical equation (29), and the resulting analytical model
results matched well with the experimental steady-state results for different recharge rates
(see Figure 20).
For the sloping aquifer, steady-state water table profiles were simulated using the
numerical model. The same numerical modeling results matched reasonably well with the
experimental results (see Figure 21).
49
0
5
10
15
20
25
30
35
40
0 2040608010120
Distance (cm)
W
a
t
e
r
Ta
bl
e
H
e
i
ght
,
h (
c
m
)
expt(w=4.2 cm/min)
ana(w=4.2 cm/min)
expt(w=5.4 cm/min)
ana (w=5.4 cm/min)
exp(w=2.9 cm/min)
ana(w=2.9 cm/min)
No flow
Boundary
Figure 20: Modeling steady-state water table profiles using analytical equation in
horizontal unconfined aquifer
50
0
5
10
15
20
25
30
35
40
0 2040608010120
Distance (cm)
W
a
t
e
r
Ta
bl
e
H
e
i
ght
,
h (
c
m
)
expt (w=5.55 cm/min)
model(w=5.55 cm/min)
expt(w=4.83 cm/min)
model(w=4.83 cm/min)
expt (w=4.42 cm/min)
model (w=4.42 cm/min)
No Flow
Boundary
Figure 21: Modeling steady-state water table profiles using numerical model in
sloping unconfined aquifer (i = 2.03 degree)
51
5.3 Transient experiments
Multiple transient experiments were conducted in both horizontal and sloping
aquifers. For transient experiments, at first water tables was raised to a steady-state level
by applying a constant recharge rate, and then the system was allowed to drain after
removing the recharge. During transient experiments, water table at the left boundary
changed with time. Therefore, we had to model the constant head boundary condition as a
known transient head boundary. All the transient experiments were simulated using
known transient head boundary condition for the left side of the tank, and no flow
boundary condition for the right side of the tank. The specific yield (
y
S ) was estimated
by fitting the numerical model results to the transient experimental results. Water table
readings were recorded at regular time intervals using a digital camera. All the
experimental data sets and the modeling results are summarized in the following section.
5.3.1 Modeling the transient head boundary at the left side of the tank
Transient head boundary for both horizontal and sloping aquifers was described
using a best fit curve of the observed transient left hand boundary data (Figures 22 and
23). These data were collected by sampling the digital images at every half-minute from
an initial steady-state condition to fully drained condition. Multiple steady-state
experiments were done using different recharge rates in horizontal aquifer. Here only one
case with a constant recharge of 5.4 cm/min is shown as an example. At first the observed
52
transient left hand boundary data were plotted against time, and then a best-fit curve was
drawn for that profile of the experimental result. The best fit curve is given by equation
(39), which was used to model the left hand side transient head boundary condition (i.e. x
= 0 cm) in horizontal aquifer. Figure 22 shows the best fit model results against observed
head data.
432
( ) 0.0061 0.1356 1.1147 4.0671 25.348ht t t t t=?+?+ (47)
h = 0.0061t
4
- 0.1365t
3
+ 1.1147t
2
- 4.0671t + 25.348
R
2
= 0.99
18
19
20
21
22
23
24
25
26
0246810
Time (min)
Left boundar
y
w
a
ter l
e
vel (cm
)
h vs t (LHS)
Poly. (h vs t (LHS))
Figure 22: Modeling transient head boundary (left hand side of the tank) in
horizontal unconfined aquifer.
53
Similarly, the equation for transient head boundary condition at the left hand
side in sloping aquifer was developed. Figure 23 also shows the best fit model results
with observed head data. The equation of that best fit curve for recharge rate of 4.42
cm/min is:
432
( ) 0.0713 0.6172 1.5002 0.3683 23.875ht t t t t=? + ? + + (48)
h = -0.0713t
4
+ 0.6172t
3
- 1.5002t
2
- 0.3634t + 23.875
R
2
= 0.9921
18
19
20
21
22
23
24
25
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
Time (min)
L
e
f
t
b
oun
dar
y
w
a
t
e
r
l
evel
(
c
m
)
h vs t (LHS)
Poly. (h vs t (LHS))
Figure 23: Modeling transient head boundary (left hand side of the tank) in
sloping unconfined aquifer (i = 2.03 degree)
54
5.3.2 Horizontal unconfined aquifer- transient experimental data and modeling results
Figure 24 shows the transient experimental data and numerical results when the
water table was allowed to drain from the steady-state condition (initial condition) with a
constant recharge rate of 5.4 cm/min in the horizontal unconfined aquifer system. The
specific yield (
y
S ) was estimated as 0.35 by fitting these numerical model results to the
transient experimental results. Then this estimated value was used in the numerical model
to fit the numerical results with another experimental results (see Figure 25 where
recharge rate= 4.42 cm/min at the steady-state).
The numerical model results show good match with the experimental results at
different time intervals. The results show that the water table drops fast at the beginning
because of high hydraulic gradient, but drops slowly at later stage due to low hydraulic
gradient. Kim et al. (2001) made similar observation. This indicates a non-linear drainage
pattern for water-table movement in the unconfined aquifer.
55
10
15
20
25
30
35
40
0 2040608010120
Distance (cm)
W
a
t
e
r Tabl
e Hei
g
ht
,
h (
c
m
)
mod(st) exp(st) mod(0.25 min) exp(0.25 min)
mod(1 min) exp(1 min) mod(2 min) exp(2 min)
mod(4 min) exp(4 min) mod( 8 min) exp(8 min)
No flow
Boundary
Figure 24: Transient behavior of water table in horizontal unconfined aquifer (k =
90 cm/min (1296 m/day),
y
S = 0.35 and w = 5.4 cm/min at the beginning
(steady state))
56
10
15
20
25
30
35
40
0 2040608010120
Distance (cm)
W
a
t
e
r
Ta
bl
e
H
e
i
g
ht
,
h
(
c
m
)
mod(st) exp(st)
mod(0.25 min) exp(0.25 min)
mod(1 min) exp(1 min)
mod(2 min) exp(2 min)
mod(5 min) exp(5 min)
No flow
Boundary
Figure 25: Transient behavior of water table in horizontal unconfined aquifer (k =
90 cm/min (1296 m/day),
y
S = 0.35 and w = 4.42 cm/min at the beginning
(steady state))
57
5.3.3 Sloping unconfined aquifer- transient experimental data and modeling results
Figure 26 also shows the transient experimental data and numerical results when
the water table was allowed to drain from the steady-state condition (initial condition)
with a constant recharge rate of 4.42 cm/min in the sloping unconfined aquifer (i = 2.03
degree).
0
5
10
15
20
25
30
35
2040608010120
Distance (cm)
Water Table Height, h (cm)
exp(st) mod(st)
exp(0.5 min) mod(0.5 min)
exp(1 min) mod(1 min)
exp(3.5 min) mod(3.5 min)
No flow
Boundary
Figure 26: Transient behavior of water table in sloping unconfined aquifer ((i = 2.03
degree, K = 90 cm/min (1296 m/day),
y
S = 0.35 and w = 4.42 cm/min at
the beginning (steady state))
58
The numerical model results also show good agreement with the experimental
results. The experimental result shows that the water profile dropped very rapidly
because of sloping effect. It is observed that the water table dropped fast at the beginning
and slowly at the later stage, similar to the horizontal aquifer case. This indicates a non-
linear pattern in the water-table movement in unconfined aquifer.
5.3.4 Mass balance accuracy of the Numerical results
Our numerical model solved a set of non-linear equations using the Picard
iterative scheme. The model employed 1x? = cm, and 0.00025t? = min to simulate the
transient profiles. The mass balance errors were less than 1%, which are tabulated in
Tables 7 and 8 for horizontal and sloping unconfined aquifers respectively.
Elapsed Time (min) Mass balance error (%)
0.25 0.46
0.5 0.5
1 0.57
2 0.65
4 0.66
8 0.68
Table 7: Mass balance errors for the water table profiles generated by the
numerical model for the horizontal unconfined aquifer
59
Elapsed Time (min) Mass Balance Error (%)
0.5 0.57
1 0.6
3.5 0.6
Table 8: Mass balance errors for the water table profiles generated by the
numerical model for the sloping unconfined aquifer (i = 2.03 degree)
5.4 Correlation between the sloping unconfined aquifer with the horizontal unconfined
aquifer
Figure 27 compares the steady-state head profiles observed in the horizontal and
sloping unconfined aquifers under similar conditions. The figure shows that the peak
head values observed in horizontal and sloping (2.03 degree slope) unconfined aquifer
occurred at the no Flow boundary. The peak in the horizontal unconfined aquifer is
approximately 1.9 cm higher than that in the sloping unconfined aquifer under identical
recharge rate of 4.42 cm/min. But the total draining time in sloping unconfined aquifer is
approximately 1.5 minute less than that in horizontal unconfined aquifer because of the
influence of the gravity effect in sloping aquifer.
60
10
15
20
25
30
35
40
0 2040608010120
Distance (cm)
W
a
t
e
r
Tabl
e H
e
i
ght
,
h (
c
m
)
Steady state (2.03 degree slope)
Steady state (Horizontal)
No flow
Boundary
Figure 27: Comparison of experimental steady state water table profile between
horizontal and sloping unconfined aquifers (i = 2.03 degree) under same
condition
61
CHAPTER 6
SUMMARY, CONCLUSIONS AND RECOMMENDATIONS
We developed a laboratory-scale aquifer model to conduct steady-state and
transient experiments in both horizontal and sloping unconfined aquifers. The
experimental results were recreated using a numerical model. A three-dimensional flow
tank (115 cm ? 2.5 cm ? 60 cm) connected to a constant head reservoir and a constant
recharge generator was used to conduct steady- state and transient experiments in both
horizontal and sloping unconfined aquifers. The hydraulics conductivity of the porous
medium was measured using 1-D column and 2-D Darcy experiments. The average
hydraulic conductivity was estimated to be 90 cm/min. At first, multiple steady-state
experiments were conducted for different recharge rates. For steady-state experiments,
constant head boundary at the left hand side and no flow boundary at the right hand side
of the tank were used. Experimentally measured steady-state water table profiles were
modeled using analytical and numerical models. Horizontal aquifer data was described
using an analytical model. For sloping aquifer, a numerical model was developed to
simulate the experimental results. During transient experiments, it was observed that the
water table at the left hand boundary changed with time. Therefore, the left hand
62
boundary was modeled as a known transient head boundary. The transient head boundary
variations observed for both horizontal and sloping aquifers were modeled by using a
best fit curve of the observed head data. This data was collected by sampling the digital
images at every half-minute from an initial steady-state condition to fully drained
condition. All the transient experiments were simulated using transient head boundary
condition for the left boundary, and no flow condition for the right boundary.
The transient experiments and numerical simulations were done for known transient
head boundary condition. The transient results show that water table drops faster at the
beginning than that in the later stage because of high hydraulic gradient at the beginning.
These results indicate the non-linear nature of water-table drainage patterns in unconfined
aquifer system. It also shows that the peak water table in the horizontal unconfined
aquifer is higher than that in the sloping unconfined aquifer (2.03 degree slope) under the
same recharge condition, hence more time required for draining in horizontal unconfined
aquifer. The novel experimental setup developed in our laboratory also allowed us to
understand the interactions among an unconfined aquifer, a stream, and a water divide
during or after rainfall event, when the water level in the near by stream varies. The
numerical model predicted better groundwater table data than that in linearized analytical
model. The proposed numerical model is a useful tool for predicting the variations in
groundwater table during or after rainfall event.
For future application, the model needs to be further validated using a field scale
stream-aquifer interaction dataset. However, there are some challenges in the field for
characterizing the hydraulic properties (e.g. hydraulic conductivity). In the field,
63
hydraulic conductivity variations will be non-homogenous and anisotropic. This would
require a careful analysis of various averaging procedures used for interpretation the field
data.
64
REFERENCES
Bear J. 1978. Hydraulics of groundwater. McGraw-Hill , Inc. NY.
Beven K. 1981. Kinematic subsurface strormflow. Water Resources Research. 17(5):
1419-1424.
Boussinesq MJ. 1877. Essai sue la theorie des eaux courants. Memories de L? Academic
des Sciences de L? Institute de France 23(1): 252-260.
Brutsaert W. 1994. The unit response of groundwater outflow from a hillslope. Water
Resources Research 30(10): 2759-2763
Chow VT, Maidment DR, and Mays LW. 1988. Applied Hydrology, McGraw-Hill,
New York.
Dupuit J. 1863. Etudes theoriques et pratiques sur le movement des eaux dans les canaux
decouverts et a travers les terrains permeables, ed. 2, Dunod, Paris.
Fetter CW. 1994. Applied Hydrogeology. Prentice Hall: Englewood Cliffs, NJ; 691.
Hall FR. 1968. Base-flow recession- A review, Water Resources Research, v. 4, pp. 973-
983.
Hazen A. 1911. Discussion: Dams on sand foundations. Transactions, American Society
of Civil Engineers 73:199.
65
Kim DJ, and Ann MJ. 2001. Analytical solutions of water table variation in a horizontal
unconfined aquifer: Constant recharges and bounded by parallel streams.
Hydrological Processes 15: 2691-2699.
Lindquist E. 1933. On the flow of water through porous soil. Premier Congres des grands
barrages (Stockholm), 81-101.
Marino MA. 1974. Rise and decline of the water table induced by vertical recharge.
Journal of Hydrology. 23: 289-298.
Marino MA. 1975. Digital simulation model of aquifer response to stream stage
fluctuation. Journal of Hydrology. 25: 51- 58.
Melville JG, and PN Sims. 1987. An analysis of hillside seepage. Ground Water. V. 25,
Number 6.
Meyboom P. 1961. Estimating ground-water recharge from stream hydrographs. Journal
of Geophysical Research. V. 66, pp. 1203-1214.
Norum DI, and Luthin JN. 1968. The effects of entrapped air and barometric fluctuations
in the drainage of porous mediums. Water Resources Research. V. 4, pp. 417-
424.
Parlange JY, Stagnitti F, Heilig A, Szilgyi J, Parlange MB, Steenhuis TS, Hogarth WL,
Barry DA, and Li L. 2001. Sudden drawdown and drainage of a horizontal
aquifer. Water Resources Research. 37(8):2097-2101.
Rocha D, Feyen J, and Dassargues A. 2007. Comparative analysis between analytical
approximations and numerical solutions describing recession flow in
unconfined hillslope aquifers. Hydrogeology Journal 15:1077-1091.
66
Rose HE. 1945. An investigation into the laws of flow of fluids through beds of granular
materials. Proceedings of the Institute of Mechanical Engineers 153:141-148.
Reddi LN. 2003. Seepage in soils: principles and applications. John Wiley and Sons.
New York.
Sanford WE, Parlange JY, and Steenhuis TS. 1993. Hillslope drainage with sudden
drawdown: closed form solution and laboratory experiements. Water Resources
Research. 29(7):2313-2321.
Schneider R. 1961. Correlations of Groundwater Levels and Air Temperatures in the
Winter and Spring in Minnesota, U.S. Geological Survey Water-Supply
Paper1539- D, pp.14.
Serrano SE, and Workman SR. 1998. Modeling transient stream/aquifer interaction with
the non-linear Boussinesq equation and its analytical solution. Journal of
Hydrology. 206: 245- 255.
Shanley JB, Kendall C, Smith TE, Wolock DM, and McDonnell JJ. 2002. Controls on old
and new water contributions to stream flow at some nested catchments in
Vermont, USA. Hydrological Processes 16: 589-609.
Sims PN. 1986. Hillside seepage modeling. Master?s Thesis. Auburn University.
Troch PA, Paniconi C, and Loon EV. 2003. Hillslope-stroage Boussinesq model for
subsurface flow and variable source areas along complex hillslopes: 1.
Formulation and characteristic response. Water Resources Research. (11),
1316, doi: 10.1029/2002WR001728.
67
Vaculin M, Khanji D, and Vachaud G. 1979. Experimental and numerical study of a
transient, two-dimensional unsaturated-saturated water table recharge problem.
Water Resources Research. 15(5), 1089-1101.
Verhoest NEC, Pauwels VRN, Troch PA, and Troch FPD. 2002. Analytical solution for
transient water table heights and outflows from inclined ditch-drained terrains.
Journal of Irrigation and Drainage Engineering 128(6):358-36
Werner PW, and Sundquist KJ. 1951. On the ground-water recession curve for large
watersheds, International Association of Scientific Hydrology Publication, 33
pp. 202-212.
68
APPENDIX
An Alternate Analysis of Depth-averaged Darcy Flux in Sloping Bed Unconfined
Aquifer:
Figure 28: Pressure at the bottom of a sloping bed unconfined aquifer
In Figure 28, the weight of fluid phase is
f
wnghxy?=??
(I)where ? is the density of water
3
M
L
? ?
? ?
? ?
, g is gravitational acceleration
2
L
T
??
??
??
, h is
i
cos
f
wi
f
w
x?
h
y?
69
height of water table perpendicular to the sloping bed, n is porosity, x? is the length [L]
along the aquifer bottom, and y? is the width [L] of the aquifer.
Component of this weight perpendicular to the sloping bed, cos
f
wi, is equal to
cosnghxy i? ??
Pressure acting on perpendicular to the sloping bed,
cos cos
ff
wiwi
F
P
A Axy
== =
? ?
(II)
where F is the force acting perpendicular to the sloping bed
2
ML
T
? ?
? ?
? ?
, and
Anxy=??
is
the area normal to the sloping bed
2
L? ?
? ?
. After substituting equation (II) in equation (III),
we get
cos
cos
nghxy i
P gh i
nxy
?
?
??
==
??
(IV)
Total head is defined as:
2
2
pv
H z
g?
=+ + (V)
where
p
?
is pressure head,
2
2
v
g
is velocity head, and z is elevation head. Since the
groundwater velocity (v) is typically a very small, we can neglect the velocity term.
Hence,
p
Hz
?
=+ (VI)
Substituting (IV) in (VI), we get
70
cosgh i
Hz
g
?
?
=+
cosH hiz=+ (VII)
From Darcy?s law,
H
vK
x
?
=?
?
(VIII)
Where v is Darcy flux or specific discharge
L
T
? ?
? ?
? ?
, K is the hydraulic conductivity
L
T
? ?
? ?
? ?
,
and
H
x
?
?
is hydraulic gradient.
Depth-averaged Darcy flux can be defined as:
qvh
?
= (IX)
where q
?
is depth averaged flux along the x direction
2
L
T
? ?
? ?
? ?
. After substituting equations
(VII) and (VIII) in equation (IX) we get the general depth averaged flux equation for
unconfined flow in a sloping aquifer as:
(cos )qKhhiz
x
?
?
=? +
?
cos
hz
qKh i
x x
?
??
??
=? +
??
??
??
cos sin
h
qKh ii
x
?
???
=? +
??
?
??
(X)