Citation
Variably saturated three-dimensional rainfall driven groundwater pumping model

Material Information

Title:
Variably saturated three-dimensional rainfall driven groundwater pumping model
Creator:
Dogan, Ahmet, 1968- ( Dissertant )
Motz, Louis H. ( Thesis advisor )
Place of Publication:
Gainesville, Fla.
Publisher:
University of Florida
Publication Date:
Copyright Date:
1999
Language:
English
Physical Description:
xii, 254 leaves : ill. ; 29 cm.

Subjects

Subjects / Keywords:
Aquifers ( jstor )
Boundary conditions ( jstor )
Evaporation ( jstor )
Groundwater ( jstor )
Hydrological modeling ( jstor )
Modeling ( jstor )
Pressure head ( jstor )
Rain ( jstor )
Simulations ( jstor )
Soils ( jstor )
Civil Engineering thesis, Ph. D ( lcsh )
Dissertations, Academic -- Civil Engineering -- UF ( lcsh )
Groundwater flow -- Mathematical models ( lcsh )
Hydrogeology -- Computer simulation ( lcsh )
Lowry Lake ( local )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Abstract:
Many water-resource management and environmental quality problems require a better understanding of the complete hydrological process, which can be described only by using a variably saturated groundwater flow model. A new saturated/unsaturated three-dimensional rainfall-driven groundwater-pumping model has been developed to understand the response of a variety of hydrogeologic systems to both natural and anthropogenic impacts. This model was designed to simulate all of the important elements of the hydrological cycle other than the runoff and seepage processes as accurately as possible using physically based assumptions and equations. The uniqueness of the model is its hydrological and hydrogeological completeness such that the model runs using rainfall and climatologic data and calculates the three-dimensional pressure distribution over the entire hydrogeologic domain. The model also calculates the potential evapotranspiration for given climatological data. In the model, the greatest effort has been devoted to an improved conceptualization of the unsaturated zone, which is a crucial part of the hydrological system in a groundwater basin. This is because the unsaturated zone plays an important role in many hydrological processes such as recharge to groundwater, surface-groundwater interaction, solute transport, and evapotranspiration. Recent advances in modeling variably saturated flow are incorporated into the model. ( , )
Thesis:
Thesis (Ph. D.)--University of Florida, 1999.
Bibliography:
Includes bibliographical references (leaves 190-205).
Additional Physical Form:
Also available on the World Wide Web; Adobe Acrobat Reader required to view and print PDF file.
General Note:
Printout.
General Note:
Vita.
Statement of Responsibility:
by Ahmet Dogan.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
43846112 ( OCLC )
002531101 ( AlephBibNum )
AMP7014 ( NOTIS )

Downloads

This item has the following downloads:


Full Text











VARIABLY SATURATED THREE-DIMENSIONAL RAINFALL DRIVEN
GROUNDWATER PUMPING MODEL

















By

AHMET DOGAN


A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA


1999



























To my beloved and devoted father, Mehmet, the most honorable person in my little
world, who was struggling with lung cancer although he had never smoked at all.
Unfortunately, he passed away and walked to his beloved God on April 13, 1999, while
his only son was away from home to complete his Ph.D. study. I believe that my father is
at rest in peace now because he lived such a beautiful life to help others, and to comfort
others just for the sake of all mighty God.















ACKNOWLEDGMENTS


I wish to express my deep and sincere gratitude to Dr. L. H. Motz, my supervisor

and mentor during my long Ph.D. study, for his helpful support and wise guidance during

this study. Special appreciation is extended to Dr. K. Hatfield for his help, support, and

helpful technical discussions. I would also like to thank Dr. W. D. Graham for her

guidance and feedback. I appreciate Dr. K. L. Campbell for his valuable help and

guidance, especially about evapotranspiration. Particular appreciation is also extended to

Dr. R. J. Thieke. He is one of the most enthusiastic teachers in our department, and I

learned a lot in his class. I would like to thank Dr. R. W. Healy, the author of the model

code VS2D, for his comments and for providing me the most current version of VS2D. I

also give thanks to God for giving me the opportunity to meet and study with these great

people at the University of Florida and finish my Ph.D. study.

I would like to express my deep appreciation to my beloved late father Mehmet,

who passed away on 13th April 1999. He sacrificed a lot to raise me and to support my

education. He was a great man in my life, and I will try to follow in his footsteps to be a

good man. I thank my devoted mother Ayse from the bottom of my heart for her prayers,

unbelievable continuous support, and encouragement. I will never forget my beloved

wife Havva's help and support. She was always there to help me and support me any

time, anywhere. She sacrificed a lot to provide me comfort and a good study

environment during this Ph.D. study. My special appreciation is also extended to my









little son Mehmet and my baby girl Ayse Hilal. I forgot my troubles and found peace of

mind when I was playing with them. Finally, special thanks are given to my sisters

Selime and Bedia for their love and prayers.















TABLE OF CONTENTS

page

A C K N O W L E D G M E N T S ................................................................................................. iii

LIST OF TABLES ................................ .. .......... ............................ viii

LIST OF FIGURES ............................. .. .......... .............................. ix

A B S T R A C T ....................................................................................................................... x i

CHAPTERS

1 IN TR OD U CTION .......................................................................... .... ....... ...............

2 L ITE R A TU R E SU R V E Y .......................................... ..... ......................... ...............6...
Historical Development of Groundwater Hydrology and Hydraulics...................6...
R research in Saturated Flow ....................................... ....................... ...............8...
U nsaturated F low Studies ....................................... ....................... ............... 16
V ariably Saturated Flow Studies....................................................... ................ 23
Available Hydrologic Computer Models ................ ...................................29

3 DERIVATION OF THE VARIABLY SATURATED GROUNDWATER
F L O W E Q U A T IO N ................................................................................... ............... 4 3
General Three-Dimensional Saturated-Unsaturated Groundwater Flow
E q u atio n ........................................................................................................... 4 3
C onceptualization ................ .............. ........................................... 43
C continuity Equation ... .... ................ ............................................... 45
Storag e T erm ...................................................................... ... ... ............ . 4 7
D arcy-B uckingham Equation.................................................... .................. 49
Governing Equation (Modified Richards' Equation)..................................52
In the saturated zone: ............................. .. ............... .. .......... ..... 53
In the unsaturated zone: ..................................... .................................. 53
H ydraulic C onductivity ......................................... ......................... ................ 55
Sink/Source T erm .... ... .. ..................................................................... 60
D eterm nation of Evapotranspiration.......................................... ................ 63
Estimation of input parameters for PET calculations ..............................66
Determination of transpiration (or root water uptake)..............................70
E v ap oration ....................................................... .................... .... . ........ 7 8
Pum ping and R charge W ells ..................................................... ................ 80


v









D rains, Sinkholes, and Springs ................................................... 81
B oundary C condition s ........................................................................ ................ 82
Specified Flux B oundary Condition ........................................... ................ 82
Specified H ead Boundary Condition .......................................... ................ 83
V ariable B oundary Condition ....................... .......................................... 83
R iver B oundary ............. .. ............... ............................................... 85
G general H ead B oundary .............................................................. ................ 86

4 MATHEMATICAL MODEL DEVELOPMENT AND NUMERICAL
SOLUTION OF THE MODIFIED RICHARDS EQUATION ....................................88
C onceptualization of the M odel .............................................................................89
Spatial D iscretization ....................................... .. ........ .......... .. .............. ... 91
T em poral D iscretization ............................................ .............. ................ 92
Finite-difference Formulation of the Governing Equation ...............................93
Mixed Form of Richards Equation and Modified Picard Iteration
S c h e m e ............................................................................................................. 9 8
B boundary Conditions ................. ........................................................ 103
Prescribed head boundaries........................................... 104
Prescribed flux boundaries........................................... 105
R iver boundary ...................................................... ............... .......... ....... 111
Overland flow and ponding............... ..........................1...12
Rainfall and evaporation boundaries...... .... .................................. 112
D ew watering of a Confined Aquifer ......................................... ............... 114
Iteration L evels......................................................................................... 114
Conductance Term s (CN i+1/2,j,k) .............................. ... ............. ............... 115
Matrix Equation Solver (Preconditioned Conjugate Gradient Method
{P C G M }) ............................................................................................................. 1 1 8

5 VERIFICATION OF THE MODEL...... ........ ...... ......................124
E x am p le 1 ............................................................................................................ 12 4
E x am p le 2 ............................................................................................................ 12 8
E x am p le 3 ............................................................................................................ 13 3
E x am p le 4 ............................................................................................................ 13 8
E x am p le 5 ............................................................................................................ 14 1
E x am p le 6 ............................................................................................................ 14 3
E x am p le 7 ............................................................................................................ 14 6

6 APPLICATION OF THE MODEL................................... ......................150
Application of the Model to a Two-Dimensional Infiltration and
Evapotranspiration Problem ................... .. ......... ....... ..... .................. 150
Application of the Model to an Unconfined Sand Aquifer Pumping Test........ 155

7 APPLICATION OF THE MODEL TO A FIELD CONDITION IN NORTH
CEN TR AL FLORID A ................................................................... ............... 162
D description of the Study Area...... ............ .......... ..................... 162
L o c atio n ......................................................................................................... 1 6 2









C lim a te ...........................................................................................................1 6 4
G e o lo g y ..........................................................................................................1 6 4
Groundw ater Hydrology ....... ........... ............ ...................... 166
A application of the M odel ................. ........................................................ 168
Selection of the M odel Area ...... .......... ........ ...................... 168
B boundary Conditions ................. ........................................................ 168
M eteorological D ata..................................... ....................... ............... 170
E vapotranspiration ................. ............................................................ 171
L ak es ........................................................ . .... .... ......... ............... 17 2
Three-Dimensional Discretization ........... ....................... 174
Tw o-D im ensional D iscretization.......................................... ................... 176
Description of Input Parameters for the Two-dimensional Simulation
of the M odel .......................................................................................... . 178
M odel R results ........................................................................................ . 179

8 SUMMARY AND CONCLUSIONS .......... ...... ......................185
Applicability Lim stations of the M odel ....... ......... ...................................... 189
F u tu re S tu d y ........................................................................................................ 18 9

LIST OF REFEREN CES ..................... ................................................................ 192

APPENDICES

A FORTRAN CODE OF VARIABLY SATURATED THREE-DIMENSIONAL
RAINFALL DRIVEN GROUNDWATER PUMPING MODEL............................208

B INPUT FILES FOR THE MODEL SIMULATION IN THE UECB.......................230

B IO G R A PH IC A L SK E T C H ...........................................................................................256















LIST OF TABLES


Table page
2.1 Summary of selected saturated-unsaturated flow models......................................31
5.1 Param eters used for exam ple 1 ......................................................... 127
5.2 Param eters used for exam ple 2 ......................................................... 132
5.3 Param eters used for exam ple 3 ......................................................... 137
5.4 Param eters used for exam ple 4 ......................................................... 140
6.1 Param eters used for the V S2D problem ............................................... ............... 153
6.2 Parameters for the unconfined aquifer pumping problem.....................................160
7.1 Geologic layers in the Upper Etonia Creek Basin (based on Motz et al., 1993).......165
7.2 Hydrogeologic units of the Upper Etonia Creek Basin (based on Motz et al.,
1 9 9 3 ) ................. ... ..................... ........ ... ................................................ . ......... 1 6 7
7.3 Regional and Local Rainfall Data During the Simulation Period ...........................171
7.4 Lake Barco Pan Evaporation Coefficients....... .......... ....................................... 172
7.5 L ake stages in the m odel dom ain ......................................................... ............... 174
7.6 Parameters used for the model application in the UECB area...............................178
8.1 Sum m ary of new m odel ................................................ ...... .......... ............. 188
B. 1 Isoil matrix for material properties of the model domain in hydrologic
simulation of UECB, where, 1: Upper Floridan Aquifer (limestone); 2, 3, 4:
Confining Unit (Hawthorn Group); 5: Surficial Aquifer (sand); and 0: no
m material ............................................................................................................. . 2 3 0
B.2 Ibound matrix for the boundary properties of the model domain in hydrologic
simulation of UECB, where, 1: Active cell; 0: inactive cell; -2: fixed head
cell for Crystal Lake, -3:fixed head cell for Magnolia Lake, 9: general head
boundary cell, 7: rainfall and evapotranspiration boundary cell............................232
B.3 Meteorological data for the period September 1, 1994-August 31, 1995 ..............234
B.4 Initial pressure heads (m) and geometric elevations in the model domain for
hydrologic simulation of the UECB....... ....... ...................... 242















LIST OF FIGURES


Figure page
3.1 Conceptualization of hydrologic system ............................................... ................ 43
3.2 R representative unit volum e of an aquifer .............................................. ................ 44
3.3 Flow chart describing the principle sink/source terms in the model........................60
3.4 Flow chart for actual transpiration calculations in the model...............................62
3.5 Flow chart for the evapotranspiration calculations (Fares, 1996)...............................65
3.6 Schematic of the plant water stress response function, ar(h) (Feddes et al.,
1978). Water uptake below hi(air entrainment pressure, saturation starts) and
above h4 (wilting point) is set to zero. Between h2 and h3 water uptake is
maximum. The value of h3 varies with the potential transpiration rate Tp .............75
3.7 Water stress function as a function of pressure head and potential
transpiration (Jensen, 1983). ....................... .............. ............ ... ................ 77
4.1 Schematic representation of the physical components and the interaction
am o n g th em .................. ....................................................................................... 9 0
4.2 V ertical discretization of the m odel ....................................................... ................ 92
4.3 Flow into and out of cell i, j, k ................................................ ......... ................ 94
4.4 Diagram for calculation of vertical conductance in case of semi-confining
units........................ ......... ............... ......... ..................... 118
4.5 PCG m ethods (Schm it and Lai, 1994) .............................................. .................... 121
5.1 Comparison of the numerical model with results of Paniconi et al. (1991)............128
5.2 Comparison of the numerical model with the analytical solution of Srivastava
and Yeh (1991)................ . .. .. ............................................. 133
5.3 Comparison of the numerical model with the analytical solution of Srivastava
and Y eh (1991) for layered soils .................................................... .... ............... 136
5.4 Comparison of the numerical model with experimental results of Vauclin et
al. (1979) ................ .. . .. ....... ........... .. ............... 139
5.5 Three-dimensional model domain description for example 5 ................................141
5.6 Water-table elevations resulting from 3-D simulation of example 5 at y = 0
for various s tim e values .................................................................. ..... .............. 142
5.7 Three-dimensional water-table recharge. Water-table elevations are shown at
the end of 8-hr rainfall of 0.148 m /hr................................................ ................. 143
5.8 Three-dimensional water-table recharge and pumping. Water-table elevations
are shown at the end of a 4-hr rainfall of 0.148 m/hr and 6.25 m3/hr pumping.......144
5.9 Three-dimensional pumping from water-table. Water-table elevations are
shown at the end of a 4-hr pumping period at the rate of 6.25 m3/hr.....................144
5.10 Three-dimensional recharge to the water-table. Water-table elevations are
shown at the end of a 4-hr injection period at a rate of 6.25 m3/hr........................145









5.11 Three-dimensional recharge to the water table. Water-table elevations are
shown at a cross-section in the x-z plane at j = 1 for different time values.
Injection well is located at (i, j) = (15, 15) at k = 2, 3, 4, 5, 6 with the rate of
6 .2 5 m 3/h r......................... .. ........................ ................................................ .. 14 5
5.12 Problem definition sketch for example 7. ...... .............. ................................. 147
5.13 Water-table position at the steady-state condition for example 7 .........................148
5.14 Three-dimensional view of the water-table for example 7. ...............................149
6.1 Description of the problem of Lappala et al. (1987)...................... ...................151
6.2 Comparison of the results of VS2D and the current model ................................155
6.3 Cross section for the unconfined aquifer pumping problem...................................157
6.4 Comparison of the pumping test results of Nwankwor et al.(1992) and the
current m odel results ............... ............................ .... .. .. .. .. ..... .......... . ........... 16 1
7.1 September 1994 water table map in the UECB and the location of the model
dom ain (Source: Sousa, 1997). ...... ............. .............. ..................... 163
7.2 Topographic surface of the m odel area ..................................................................... 164
7.3 The model boundaries and September 1994 water table map. ................................169
7.4 L ake levels in the m odel dom ain. ......................................................... ............... 173
7.5 Horizontal discretization of the three-dimensional model domain........................ 175
7.6 Vertical discretization of the two-dimensional model domain...............................177
7.7 The model results versus the observed data in the well C520 during the period
of Sepem ber,1994-Septem ber, 1995 .................................... ............... ............... 181
7.8 The Rainfall and Evapotranspiration components in the model area. .....................182
7.9 Total head contours at tim e =150 days ................. ............................... ................ 183
7.10 Moisture content profiles at different times of the simulation............................. 184















Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy

VARIABLY SATURATED THREE-DIMENSIONAL RAINFALL DRIVEN
GROUNDWATER PUMPING MODEL

By

Ahmet Dogan

December 1999

Chairman: Louis H. Motz
Major Department: Civil Engineering

Many water-resource management and environmental quality problems require a

better understanding of the complete hydrological process, which can be described only

by using a variably saturated groundwater flow model. A new saturated/unsaturated

three-dimensional rainfall-driven groundwater-pumping model has been developed to

understand the response of a variety of hydrogeologic systems to both natural and

anthropogenic impacts. This model was designed to simulate all of the important

elements of the hydrological cycle other than the runoff and seepage processes as

accurately as possible using physically based assumptions and equations.

The uniqueness of the model is its hydrological and hydrogeological completeness

such that the model runs using rainfall and climatologic data and calculates the three-

dimensional pressure distribution over the entire hydrogeologic domain. The model also

calculates the potential evapotranspiration for given climatological data. In the model,









the greatest effort has been devoted to an improved conceptualization of the unsaturated

zone, which is a crucial part of the hydrological system in a groundwater basin. This is

because the unsaturated zone plays an important role in many hydrological processes such

as recharge to groundwater, surface-groundwater interaction, solute transport, and

evapotranspiration.

Recent advances in modeling variably saturated flow are incorporated into the

model. The model simulates the hydrogeologic system by solving the nonlinear three-

dimensional mixed form of the Richards equation using the modified Picard iteration

scheme and preconditioned conjugate gradient method. A new convergence criterion is

used for faster convergence in the iterations. The model treats the complete subsurface

regime as a unified whole, and flow in the unsaturated zone is integrated with saturated

flow in the underlying unconfined and confined aquifers. The model has the capability to

simulate pumping from the aquifer and artificial recharge. A transpiration and an

evaporation model are integrated into the groundwater flow equation as sink terms.

Input data for the model are the hydrogeologic and geometric properties of the

flow domain, meteorological data, vegetative cover, and soil moisture characteristics.

The output is in the form of groundwater heads, moisture contents, and actual

evapotranspiration. The model has been verified against other model results from the

literature.















CHAPTER 1
INTRODUCTION

The management and control of our water resources requires the conception,

planning, and execution of designs to make use of water without causing harm to the

environment. Approximately forty percent of the water used for all purposes in the

United States is derived from groundwater sources (Heath, 1983). Groundwater is a vital

and very vulnerable source of water supply. The main source of recharge for the

groundwater is precipitation, which may move through the soil directly to the

groundwater or it may enter surface-water bodies such as rivers, streams, lakes, and

wetlands and percolate from these water bodies to the groundwater. Interception,

depression storage, evapotranspiration, and soil moisture capacity must be satisfied

before any large amount of water can percolate to the groundwater. Precipitation can

supply large quantities of water for groundwater easily in such places where sandy soils,

flat topography, and high water-tables occur, i.e., in Florida. The surface-water and the

groundwater are strongly interrelated, and the use of one source may affect the water

available from the other source. Both surface water and groundwater should be

considered together in plans for water-resources development.

The groundwater-surface-water interaction process involves infiltration,

evapotranspiration, runoff, and seepage between streams and aquifers. A surface-water

model or a groundwater model alone cannot accurately simulate this process. Instead, a

complete hydrological system model is required, which can simulate the rainfall-runoff









relation, evapotranspiration, unsaturated flow, saturated flow, seepage, and pumping from

the aquifer.

There are two main reasons to develop and rely upon hydrologic models, i.e., to

understand why a flow system is behaving in a particular observed manner and to predict

how a flow system will behave in the future. In addition, models can be used to analyze

hypothetical flow situations in order to gain generic understanding of that type of flow

system. The first step in studying a groundwater system is to develop a conceptual

model, which describes the real hydrogeologic system. After conceptualization of the

real system, a mathematical model is developed to solve some form of the basic equations

of groundwater flow. Mathematical models can be classified as analytical or numerical

models, depending on the solution technique. Analytical models can be solved rapidly,

accurately, and inexpensively. Numerical models sometimes must be used when there is

a very complex hydrogeologic system where hydrogeologic and hydraulic properties vary

within the model area. Numerical solutions to the groundwater flow equations require

that the equations be reformulated using one of the approximation techniques, e.g., finite-

difference, finite element, or the method of characteristics. The requirements of water

resources planning have made numerical model simulations of the hydrologic response of

groundwater basins a very important technique. Successful resolution of many water-

resources management and environmental quality problems is possible through a better

understanding of the hydrological processes that take place near the ground surface, i.e.,

in the unsaturated, or vadose, zone.

A new saturated/unsaturated three-dimensional rainfall-driven groundwater-

pumping model, described in this dissertation, has been developed to understand better









groundwater level fluctuations and help to make reasonable groundwater policies. The

model was designed to simulate all of the important elements of the hydrological cycle as

accurately as possible in a manner that all assumptions and equations are physically

based. The uniqueness of the model is its three-dimensional hydrological and

hydrogeological completeness and better conceptualization of the unsaturated zone. The

unsaturated zone is a crucial part of the hydrological system in a basin. It plays an

important role in many modeling applications, e.g., for recharge estimation, surface-

groundwater interaction, solute transport, and evapotranspiration calculations. Therefore,

in the model, the main emphasis is given to simulation of the unsaturated zone, the

infiltration process, evapotranspiration, and the root water uptake process.

The model utilizes the finite-difference technique to solve the three-dimensional

form of the variably saturated groundwater flow equation. The finite-difference grids can

be generated as variable or constant size. The upper boundary in the model is at ground

surface, and the upper boundary conditions are determined using soil and meteorological

data. The upper boundary condition can be a positive flux boundary (i.e., before ponding

occurs) or a fixed head (i.e., after ponding occurs) during a rainfall event. It can be a

negative flux boundary or a fixed head boundary during the evaporation process. The

model treats the complete subsurface regime as a unified whole, and the flow in the

unsaturated zone is integrated with saturated flow in the underlying unconfined and

confined aquifers. This is achieved by solving the complete three-dimensional nonlinear

Richards equation (1931) throughout the whole flow domain. The model allows

modeling of heterogeneous and anisotropic geologic formations. It has the capability to

simulate anthropogenic effects such as pumping from an aquifer and artificial recharge.









A plant root water uptake (transpiration) model and an evaporation model are included as

sink terms in the groundwater flow equations. The model also includes a module to

calculate the potential evapotranspiration values for a given location and climatologic

data based on the Priestly and Taylor (1972) equation.

The required input data for the model are hydrogeologic and geometric properties

of the flow domain, meteorological data, vegetative cover, and soil type data for the

calculation of evapotranspiration, rainfall data, and soil-water characteristics. The output

provides groundwater heads in terms of pressure head, moisture-content profile in the

unsaturated zone, actual evapotranspiration, and exchange of water between surface-

water and groundwater systems.

A groundwater setting in north Florida was selected as an example of the model's

application. Florida has a unique hydrogeologic character with its flat topography, heavy

subtropical rainfalls, wetlands, high water-tables, and sandy soils, which cause significant

interactions between groundwater and surface-water systems. Florida's continuing

population growth has resulted in an increasing demand on the water supply. This

increasing demand mainly will be met using the state's groundwater resources. However,

excess usage of groundwater for public water supply, irrigation, and industry has led to

negative impacts, including saltwater intrusion, the lowering of lake and wetland water

levels, and the reduction of spring flow and stream flow. This problem is especially true

for north Florida.

Using the deterministic, fully distributed, physically based integrated hydrological

model, the effects of human interventions and effects of naturally varying recharge in the






5


form of rainfall patterns on the hydrological cycle can be determined. Using this model, a

more informed basis for policy and decision-making can be created.















CHAPTER 2
LITERATURE SURVEY


Historical Development of Groundwater Hydrology and Hydraulics


Although groundwater has been used since early times, an understanding of the

origin of groundwater as related to the hydrologic cycle was established only in the later

part of the seventeenth century. Several incorrect hypotheses explaining the occurrence

of groundwater were given by such early Greek philosophers and historians as Homer

(about 1000 BC), Anaxagoras and Herodotus (fifth century BC), Plato (427-347 BC), and

Aristotle (384-322 BC). Plato thought that one huge underground cavern in the earth was

the source of all rivers and that water flowed back from the ocean to this cavern.

Surprisingly, however, Plato's opinion includes an accurate description of the hydrologic

cycle (Baker and Horton, 1936). The Roman philosophers followed the Greek teachings

and contributed little to the subject. These hypotheses were unquestioned until the end of

the seventeenth century. The Roman architect Marcus Vitrivius (15 BC-58 AD) was

probably the first in the recorded history to have a correct grasp of the hydrologic cycle.

He realized that the mountains receive a large amount of water from melting snow that

seeps through the rock strata and emerges as springs at lower elevations. Al-Biruni (973-

1048) accurately explained the mechanics of groundwater movement as well as the

occurrence of natural springs and artesian wells "on the principle of water finding its own

level in communicating channels" (Kashef, 1986). Bernard Palissy (1509-1589) is









recognized as the first in modem history to explain the hydrologic cycle, the origin of

springs, and the relationship between wells and rivers (Cap, 1961). The first field

measurements were made by Pierre Perrault (1608-1680). He studied evaporation and

capillary rise and measured the rainfall and runoff of the upper drainage basin of the

Seine River in France (De Wiest, 1965). The findings of Perrault were verified several

years later by Edm Maiotte (1620-1684), whose report appeared in 1686 after his death.

Outstanding documents on the subject of artesian wells were written in 1715 by Antonio

Vallisnieri, President of the University of Padua, Italy (De Wiest, 1965).

In the nineteenth century, quantitative measurements were initiated by Darcy

(1856) and supplemented by the analytical work of Dupuit (1863), Thiem (1906), and

Forchheimer (1898). This work stimulated groundwater research in the twentieth century

and shifted groundwater hydrology from a descriptive subject to a more rigorous

analytical science (Kashef, 1986).

There have been three revolutions in physical hydrogeology: the historic set of

experiments carried out by Darcy (1856) in Dijon, France; the transient well-hydraulics

analysis by C.V. Theis in 1935; and the introduction of large digital computers in the

early 1960s. Darcy developed an empirical law on which nearly all subsequent work has

been based, and Theis developed a methodology for both the in-situ measurement of

hydrologic properties of geologic formations and the prediction of the response of

groundwater systems to pumping. Digital computers provide the means for assessment of

groundwater resources on a regional scale within the context of the full hydrologic cycle

(Freeze and Back, 1983).









Research in Saturated Flow


Two- and/or three-dimensional water flow through saturated porous media has

been known in its steady-state form since the work of Forchheimer (1898) in the late

nineteenth century. His understanding was based on an analogy with the heat-flow

equation. Theis invoked the same analogy in 1935 in presenting a solution to the

transient form of the flow equation, although he did not present the fundamental

differential equation itself. Since the movement of fluids in geological materials can be

understood based on treating fluid flow as a process mathematically analogous to heat

conduction in solids, the working mathematical model for the transient groundwater flow

is the partial differential equation of heat conduction, originally proposed by Fourier.

Fourier's theory was published in 1822 with additional works of Laplace, Lagrange,

Monge, and Lacroix. Darcy was aware of the studies of Fourier, Ohm, and Poiseuille and

made use of them in his work (Narasimhan, 1998a).

During the latter half of the nineteenth century, Boussinesq, Dupuit, Forchheimer,

Adolph Thiem, and Gunther Thiem made important contributions to the development of

the science. Dupuit (1863) developed a linear constitutive law, similar to Darcy's, based

on hydraulic theory rather than experimental evidence. He also produced the first formal

mathematical analysis of a groundwater hydraulics problem, that of radial flow toward a

pumped well in an unconfined aquifer. The assumptions invoked in his analysis, namely,

that the hydraulic gradient is equal to the slope of the water-table and that it is invariant

with depth, have come to be known as the Dupuit assumptions, and methods based on

these assumptions are still in wide use today.









Chamberlin (1885) is generally recognized as initiating the science of

hydrogeology in the United States. He outlined the seven prerequisites for artesian flow

and described the hydrogeologic properties of water bearing beds in his 1885 report. If

Gauss was the "prince of mathematicians," then surely Forchheimer was the prince of

groundwater hydraulics" (Freeze and Back, 1983). Forchheimer (1898) was the first to

note the analogy between groundwater flow and heat flow, and he was the first to use the

Laplace equation in the description of steady-state groundwater flow. He clarified the

Dupuit assumptions and recognized that steady-state flow in unconfined aquifers under

the Dupuit assumptions would obey the Laplace equation with respect to the square of the

hydraulic head rather than the hydraulic head itself.

Dupuit's formula for the discharge from a well in an unconfined aquifer required

advanced knowledge of the radius of the zone of influence at steady-state. Adolph Thiem

carried out extensive field investigations to clarify the controls on the radius of influence

in 1870. His son Gunther Thiem (1906) was the first to recognize that Dupuit's equations

could be applied at any two points on the profile of the cone of depression around a well.

This realization led to the first inverse application of a solution to the steady-state flow

equation and, hence, to the first use of pumping tests as a practical tool for in-situ

measurement of the hydraulic properties of geologic formations.

During the last part of the nineteenth century, nearly the same important

developments were duplicated in the United States because of the poor interchange of

information between Europe and the United States. C. S. Slichter of the U. S. Geological

Survey, working twelve years after Forchheimer and apparently unaware of his work,

utilized the same heat-flow literature to arrive at the Laplace equation and flow-net









construction. Another important contribution of Slichter was the investigation of the

physical significance of hydraulic conductivity, which was treated only as an empirical

constant by Darcy. He identified the geometric and viscous drag components of hydraulic

conductivity.

In the evolution of the ideas pertaining to the flow of fluids in geological media,

the period 1920-1940 must rank as truly remarkable (Narasimhan, 1998a). Oscar

Meinzer of the U.S. Geological Survey was one of the most famous hydrogeologists

during the early decades of the twentieth century in the United States. His two classic

water-supply papers (Meinzer, 1923 and 1928) are still reprinted and widely used today

(Freeze and Back, 1983). His major contribution to the understanding of the physics of

groundwater flow came in his 1928 paper on the compressibility of the artesian aquifers

wherein he invoked the effective stress principle. Meinzer recognized that the water in a

confined aquifer supports part of the overlying load and that aquifers compact when fluid

pressure is decreased. Although Terzaghi (1925) developed the basic concept of

effective stress in a laboratory soil column, Meinzer's realization that the same concept

applied to aquifers was a major breakthrough.

Weber (1928) made a successful attempt to analyze the unsteady flow of water to

a fully penetrating gravity well in an unconfined aquifer for the first time. In the 1930s,

the results of mathematical and experimental studies in the petroleum reservoir

engineering field were utilized by researchers in the groundwater field. Muskat (1934)

presented a detailed analysis of transient flow of compressible fluids in oil and gas

reservoirs. In the field of groundwater hydrology, Theis (1935) set up and obtained a

solution to the parabolic equation of groundwater flow similar to that of Muskat (1934).









He verified the credibility of his model by applying it to Wenzel's Grand Island,

Nebraska field data from an unconfined aquifer. Wenzel (1942) brought Theis' theory

into practical use by publishing a table of the exponential function values that appeared in

Theis' solution. Theis' work has been a milestone in groundwater hydrology and his

model is still used frequently today. Theis was careful in his paper to spell out the

assumptions on which his method was based, i.e., it applies to an idealized aquifer

configuration. The history of the subsequent development of the methodology of aquifer

hydraulics is largely a history of the systematic removal of his assumptions one by one.

Jacob (1946) extended Theis' method to heterogeneous media when he published

a paper on radial flow to a leaky aquifer, which opened up a new area of research relating

to multiple aquifer systems in groundwater hydrology and petroleum engineering.

The auger-hole methods and piozemeter methods were pioneered by Kirkham and

coworkers (Kirkham, 1946; Luthin and Kirkham, 1949; van Bavel and Kirkham, 1948).

These methods improved the estimation of the hydraulic conductivity of the saturated soil

below the water-table and are still being used.

Boulton (1954) pioneered the analysis of unconfined aquifers. He investigated the

transient flow of water to a well in an unconfined aquifer. Instead of solving the highly

complex flow process in the unsaturated zone embodied in Richards' equation, Boulton

simplified the effect of the unsaturated zone by introducing an empirical constant that

accounted for the delayed yield from the storage. As an approximation, he assumed that

the drainage from the unsaturated zone was an exponential function of time. The

resulting governing equation was solved for potentials within the saturated domain, while

yet approximately accounting for the contribution from the unsaturated zone by means of









a time dependent source term. His model still continues to be used by hydrogeologists

with minor modifications and extensions (Narasimhan, 1998a).

The effects of anisotropy and heterogeneity of the aquifers on flow were

investigated by Maasland (1957). Maasland also outlined the relationships between

stratified heterogeneous systems and homogeneous anisotropic systems in his paper.

During the early 1960s, doubts were expressed about the validity of Jacob's

development of the groundwater flow equation. The doubts were centered around the

fact that the effective stress laws he invoked assumed that only vertical stress existed. A

full analysis should have dealt with the interaction between a three-dimensional stress

field and a three-dimensional fluid flow field. Hydrogeologists discovered that Biot

(1941, 1955), a physicist working in a petroleum research institute, had already coupled a

three-dimensional stress field with the fluid -flow field. His work was interpreted in terms

of hydrogeological notation by Verruijt (1969) and placed in the context of earlier

groundwater development. In the mean time, De Wiest (1966) improved the Jacob

equation with respect to the variation of hydraulic conductivity with effective stress but

not with respect to the storage side of the equation. Cooper (1966) clarified the

relationship between the development of the flow equation in fixed coordinates and

deforming coordinates. Cooper concluded that Jacob's equation was correct for almost

all practical applications. Cooper and a group of hydrogeologists led by him made many

contributions to groundwater hydrology. These include interpretation of data from a slug

test (Cooper et al., 1967), analysis of transient pressure data from an anisotropic aquifer

(Papadopulos, 1965), transient flow of water to a well of large diameter (Papadopulos and

Cooper, 1967), and the response of a well to seismic waves (Cooper et al., 1965).









The study of leaky aquifers was pioneered by Jacob and his student Hantush.

Hantush (1960) considered the effects of aquitard storage in his leaky aquifer solution.

Hantush (1964) provided a comprehensive summary of developments related to leaky

aquifers as well as other aquifer configurations in his paper "Hydraulics of Wells". Toth

(1963) produced a set of analytical solutions to the steady-state boundary value problem

representing regional flow in a vertical profile. Neuman and Witherspoon (1969)

presented a complete solution that considers both water released from storage in the

aquitard and drawdowns in the hydraulic head in the unpumped aquifer.

A significant research milestone of the 1960s was the development of numerical

models. The era of the digital computer had started and computer development was

advancing with incredible rapidity. The digital computers provided the possibility of

solving transient flow problems in complex geological systems, which are impossible to

solve in closed form solutions. The early numerical solutions were based on the finite-

difference method and the method of relaxation, both of which were known before the

advent of computers. Stallman introduced finite-difference concepts into the

hydrogeological literature in 1956. Much later, Nelson (1968) used the finite-difference

method to study the inverse problem studies. The finite-element method (Clough, 1960),

which was initially designed for solving structural engineering problems, was soon

adapted to solve steady-state and transient problems of groundwater flow (Javandel and

Witherspoon, 1968).

Remson et al. (1965) helped popularize the computer modeling approach by

developing a steady-state computer model to predict the effects of a proposed surface-

water reservoir on the heads in an unconfined regional aquifer. Freeze and Witherspoon









(1966) presented numerical solutions that allowed consideration of more complex water-

table configurations and geologic environments.

By the early 1970s, computer simulation of aquifers was widely used in water

resources evaluations. This advance resulted largely from the development,

documentation, and availability of two aquifer simulation programs, the first by Pinder

and Bredehoeft (1968) of the U.S. Geological Survey, and the second by Prickett and

Lonnquist (1971) of the Illinois State Water Survey. The U.S. Geological Survey model

has been continually updated over the years. The first attempt to create a complete

hydrologic response model was made by Freeze and Harlan in 1969. Freeze (1971), who

was one of the pioneer numerical modelers, developed a three-dimensional variably

saturated transient groundwater flow model. His model was in finite-difference form and

used the line successive over relaxation method to solve the governing equation. The

finite-difference models of Freeze (1971) and Cooley (1971) however are not robust

because they incur numerical instabilities and convergence difficulties (Clement et al.,

1994). In the late 1970s, research emphasis was shifted from resources development

issues to environmental issues pertaining to chemical contamination. Since the

contaminant transport path typically goes through the unsaturated zone, soil scientists and

agricultural engineers began to investigate unsaturated soil characteristics and flow

processes in the vadose zone. Most of the researchers focused on unsaturated-saturated

flow problems. Neuman (1973), Brandt et al. (1971), and Haverkamp et al. (1977) are

among those researchers.

In the 1980s, topics such as leaky aquifers and unconfined aquifers gradually

receded from researchers' focus of attention. Interest steadily grew in characterizing flow









processes in the vadose zone, which is the path between wastes deposited at the land

surface and the water-table at depth. In the latter part of the 1980s, the motivation was to

develop better computer models and to search for better numerical techniques to solve

governing nonlinear partial differential equations. Advancements in computer technology

eased the researchers' job and motivated them to attempt to solve more complex,

challenging, and time consuming groundwater problems. Parallel to these advancements,

studies on numerical solution techniques increased rapidly. New numerical methods

were developed and applied in models. The boundary integral method (Liggett and Liu,

1983) and the analytic element method (Strack, 1989) were relatively new techniques that

were applied in models. Many sophisticated groundwater models were developed in the

late 1980s. The most well known of these groundwater models, MODFLOW, was created

by McDonald and Harbaugh (1988) of U.S. Geological Survey. MODFLOW is still

widely used by hydrogeologists.

In the 1990s, more challenging problems begun to be dealt with. Attempts were

made to couple variably saturated flow models, root water uptake models, and

groundwater models to simulate the complete hydrological process. With the help of

high-speed computers, hydrogeologists started modeling the surface-water groundwater

interaction process, and surface-water flow models were coupled with the groundwater

models. MODFLOW and BRANCH models were coupled by Swain and Wexler (1992)

of U.S. Geological Survey in 1992 to simulate non-steady river flow interaction with

groundwater in a successful coupled model referred to as MODBRANCH. Yeh et al.

(1996) developed a three-dimensional finite-element saturated unsaturated flow and

transport model. During the 1990s, contaminant transport and consequently unsaturated









saturated flow studies became very important because of increasing environmental

awareness and multi million dollar support of government agencies such as the U.S.

Department of Energy (DOE), U.S. Nuclear Regulatory Commission (NRC), U.S.

Environmental Protection Agency (EPA), and others. Contaminant transport is beyond

the scope of this dissertation. The studies for unsaturated flow and variably saturated

flow problems are described in the subsequent sections.

Another very important development was the introduction of Geographic

Information Systems (GIS) to water resources research in the 1990s. With the help of

GIS methods and the graphical interface programs, numerical models became very user-

friendly in terms of input data and post processing of the output data. This new technique

provided an interactive environment in which model grids, spatially referenced to a base

map, can be generated on the computer screen and the model results can be seen on the

screen immediately. It provides the capability for modelers to create, apply, and revise

groundwater models quickly and in a way never possible before. The first example of

this type of model is GWZOOM by Yan and Smith (1995), who created a system based

on GIS that works interactively with MODFLOW. Application of GIS to groundwater

problems is a very rapidly growing research area today, and it will be one of the primary

interests of researchers in the 21st century.

Unsaturated Flow Studies


The unsaturated zone in the hydrologic cycle transmits water falling or ponded on

the land surface to underlying groundwater or temporarily stores water near the land

surface for plant use. The first researchers who dealt with the unsaturated zone were soil









physicists. Later, agricultural engineers investigated the behavior of the unsaturated zone

around the plant root zone above the water-table. Starting with Terzaghi (1925), civil

engineers and geotechnical engineers became interested in the unsaturated zone to deal

with seepage and ground settlement problems.

The first research about unsaturated flow dates back to the early twentieth century.

This was conducted by Edgar Buckingham (1867-1940), who was a physicist at the

Physical Laboratory of the Bureau of Soils, U.S. Bureau of Agriculture. His theoretical

and experimental studies on the dynamic movement of soil gases and soil moisture led to

a major contribution to the foundation of soil physics. His first paper was published in

1904, but his major contribution to unsaturated flow research was his second paper,

which was published in 1907. This paper reported the results of studies on the movement

of soil moisture. Based on the works of Fourier and Ohm, Buckingham rigorously

defined the concept of capillary potential and proposed that the steady flux of moisture

through an unsaturated soil is directly proportional to the gradient of the potential, with a

coefficient of proportionality being a function of capillary potential. The mathematical

form of this statement was much the same as that of Darcy's law, except that the

parameter of proportionality was recognized by Buckingham to be a function of capillary

potential. It is remarkable that Buckingham, who was probably unaware of Darcy's work

(Sposito, 1987) gave a theoretical basis for Darcy's empirical law and extended the law to

the unsaturated zone. Buckingham provided a paradigm and unified the flow process in

the unsaturated and saturated zones. Some soil physicists persuasively argue that the

phrase "Darcy-Buckingham's law" should be used in place of Darcy's law (Narasimhan,

1998b). Buckingham appears to be the first researcher to address the transient movement









of water in the subsurface, and he is also widely known for developing the dimensional

analysis "pi theorem" (Buckingham, 1914).

At about the same time, Green and Ampt (1911) proposed a simple approximation

to quantify the vertical infiltration of water into an unsaturated soil. The Green and Ampt

idealization assumes that a sharp, piston-like zone of saturation advances with time as

water infiltrates into a soil. This approximation is still widely used.

Gardner and Widtsoe (1921) attempted to quantify the unsteady moisture

movement in unsaturated soils in terms of a transient diffusion equation analogous to

Fourier's transient heat conduction equation. They did not achieve satisfactory agreement

between experiment and theory, because they did not account for the dependency of

hydraulic conductivity on capillary potential suggested a decade earlier by Buckingham.

They tried to fit experimental data to a linear partial differential equation, when in fact a

nonlinear parabolic equation should have been used. In 1924, Terzaghi experimentally

studied the deformation of water-saturated clays and established the relationship among

external stresses, pore-water pressure, and deformation. Although his paper is classified

under the soil mechanics discipline, he proceeded to write down and solve the equation

for transient movement of water in a one-dimensional clay column by analogy with the

heat conduction equation (Narasimhan, 1988a). Tensiometers had become well

developed by the efforts of Willard Gardner and his coworkers in the late 1920s. Gardner

et al. (1922), in an abstract, published the first reference to the tensiometer (Narasimhan,

1998a), an instrument that has played a vital role in the evaluation of modern soil physics.

Because of the tensiometer, routine measurements of moisture-content and its relation to

capillary pressure had become possible (Richards, 1928). Combining Buckingham's









(1907) work on the equation of water motion in unsaturated soils with the newly available

soil moisture retention curves, Richards (1931) formally wrote down, for the first time,

the nonlinear partial differential equation describing transient flow of water in unsaturated

soils. He defined the moisture capacity as the slope of the moisture-content versus

capillary pressure curve. The Richards equation remained unsolved for nearly two

decades because of its nonlinearity. Klute (1952), Philip (1955), and others began to

obtain solutions for the Richards equation under highly simplified conditions using

numerical methods in the early 1950s.

Richards et al. (1956) demonstrated a method by which the hydraulic conductivity

function could be estimated in the field by measuring the depth profile of gauge pressure

head as well as moisture-content as a function of time during the redistribution of soil

moisture immediately following an infiltration event. In this experiment, the soil

moisture distribution was measured rather laboriously by the gravimetric method.

Gardner and Kirkham (1952) used neutron scattering to estimate quantitatively the soil

moisture, and this process developed into a workable field neutron probe, which is very

useful in measuring the soil moisture profile. During the 1960s, the field method of

Richards and his coworkers was improved by other researchers by utilizing the neutron

probe.

Gardner (1957) found that there is an exponential relationship between the

hydraulic conductivity and gauge pressure head over a limited range of gauge pressure.

This relationship made it possible for him to solve the Richards equation. The works of

Gardner (1957) and Philip (1955) continue to influence present day research relating to

hydraulic characterization of unsaturated soils in the field. Veihmeyer and Hendrickson









(1955) presented a comprehensive literature review about the relationship between

transpiration and the soil moisture. In 1957, Philip published his famous "Theory of

Infiltration" series in which he developed an infiltration equation based on the Richards

equation and Klute's (1952) equations for finite and infinite soil profiles. In 1957,

Gardner developed several steady-state solutions of the unsaturated moisture equation

with the application to evaporation from a water-table. Gardner (1960) did another study

about the plant root and soil moisture interaction and its dynamic behavior. The relation

of the root distribution to water uptake as a function of soil suction and water availability

was described by Gardner (1964). Brooks and Corey (1964) developed analytical

expressions to define the relationship between unsaturated hydraulic conductivity and soil

moisture-content based on the statistical predictive conductivity model of Burdine (1953).

Brooks and Corey (1964) obtained fairly accurate results with their equations.

In the 1970s, research continued to find better and more effective solutions to the

Richards equation to describe the infiltration process in the soil. Ripple et al. (1972)

investigated the relationship between bare soil evaporation and a high water-table.

Meteorological and soil equations of water transfer were combined in order to estimate

approximately the steady-state evaporation from bare soil under conditions of a high

water-table. Warrick (1975) described a one-dimensional linearized analytical solution to

the moisture flow equation for arbitrary input using simplified boundary conditions.

Mualem (1976) derived a new model for predicting the hydraulic conductivity based on

the soil-water retention curve and the conductivity at saturation. Mualem's derivation

leads to a simple integral formula for the unsaturated hydraulic conductivity, which

makes it possible to derive closed-form analytical expressions, provided suitable









equations for the soil-water retention curves are available, van Genuchten (1980)

developed a closed-form equation for predicting the unsaturated hydraulic conductivity

using Mualem's (1976) model. van Genuchten's closed form expression is still used

widely by hydrogeologists.

Haverkamp et al. (1977) developed one-dimensional moisture-content based

numerical models to solve the Richards equation for infiltration. Six different numerical

models were developed and compared to each other in terms of numerical errors and

computer time requirements. Those models were verified using experimental values and a

comparison was made between one of the six models and a calculated analytical solution

of Philip (1957) and Parlange (1971).

During the 1980s, groundwater contamination came into sharp focus because of

leaky gasoline tanks and other industrial wastes. Several government agencies gave large

financial support to contaminant transport studies, which require unsaturated flow

analysis. This motivation increased the number of investigations concerned with vadose

zone hydrology. The use of numerical models for simulating fluid flow and mass

transport in the unsaturated zone became increasingly popular. A lot of effort was made

to develop these kinds of models. A comprehensive list of numerical codes for single-

phase (water) flow in the vadose zone is given by Stephensen (1995).

Knoch et al. (1984) developed a one-dimensional physically based computer

model for predicting direct recharge to groundwater. Yeh et al. (1985) presented the

results of a stochastic analysis of unsaturated flow in heterogeneous soils. Results of their

stochastic theory for flow in heterogeneous soils were compared with experiments and

field observations. Effects of anisotropy on recharge, irrigation and surface runoff









generation, and waste isolation were discussed in the paper. Broadbridge and White

(1988) presented an analytical solution for a nonlinear diffusion-convection model

describing constant rate rainfall infiltration in uniform soil and the application of the

solutions (White and Broadbridge, 1988). Hills et al. (1989) developed a one-

dimensional model for infiltration into very dry soil. Numerical instability and

convergence problems caused by a sharp wetting front in very dry soil conditions were

dealt with. A two-step Crank-Nicholson procedure was developed, in which the first step

estimates the material properties and the second step uses the temporal averages of these

properties to calculate the unknown pressure head or moisture-content.

In the 1990s, more complex geometric situations such as layered and

heterogeneous soils, and variably saturated, multidimensional problems were

investigated. Warrick and Yeh (1990) presented a one-dimensional solution for a layered

soil profile. Ross (1990) developed efficient numerical methods for infiltration using the

Richards equation, proposing the "Advancing Front Method" as a better method for

infiltration redistribution and drainage problems. Celia et al. (1990) developed a general

mass-conservative numerical solution for the unsaturated flow equation. They reported

that the pressure based form of the Richards equation generally yields poor results and is

characterized by large mass balance errors. Conversely, mass is perfectly conserved in

numerical solutions based on the mixed (head and moisture-content) form of the Richards

equation.

Yeh and Harvey (1990) investigated various approaches for determining effective

conductivity values for heterogeneous sands and compared them to laboratory

measurements. Warrick et al. (1990, and 1991) and Broadbridge and Rogers (1990)









presented analytical solutions for steady and transient infiltration processes by solving the

Richards equation. An analytical solution for one-dimensional, transient infiltration in

homogeneous and layered soils was developed by Srivastava and Yeh (1991) using

exponential functions describing hydraulic conductivity, pressure head, and moisture-

content. Paniconi et al. (1991) evaluated iterative and non-iterative methods for the

solution of the Richards equations. They presented four different non-iterative solution

techniques and concluded that a second order accurate two-level "implicit factored" non-

iterative technique is a good alternative to iterative methods. Barry et al. (1993)

developed a class of exact analytical solutions for the Richards equations. Tracy (1995)

developed a three-dimensional analytical solution for unsaturated flow using a simplified

boundary condition. Chang and Corapcioglu (1997) presented a study on the effects that

roots have on water flow in unsaturated soils and included a root distribution model.

Variably Saturated Flow Studies


The necessity of modeling variably saturated flow was brought about by drainage

problems, which include both saturated and unsaturated flow. The first, relatively simple

variably saturated flow research was done in the field of agricultural engineering in the

late 1950s. These studies were limited to the one-dimensional drainage problems (e. g.,

Day and Luthin, 1956).

Starting in the late 1960s, rapidly developing computer technology motivated

researchers to analyze the entire hydrologic cycle. Sophisticated numerical methods and

high-speed computers made modeling of the entire geologic formation possible, from

ground surface to the impermeable bottom of a confined aquifer. The solution of the









highly nonlinear governing equation with very complex boundary conditions and

heterogeneous geologic formations combined with the transient nature of the problem,

required using very powerful high-speed computers.

Finite-difference approximations were widely used in several early studies. Most

of these variably saturated models were one-dimensional (e. g., Freeze and Harlan, 1969).

Transient numerical models that integrate the saturated and unsaturated zones

were pioneered by Rubin (1968), who developed the first multidimensional variably

saturated model. He developed a model using Darcy's flow equation (but it was actually

the Richards equation) for two-dimensional, transient water movement in a rectangular

unsaturated or partly unsaturated soil domain. He used alternating direction and implicit

difference methods. His paper was followed by several other two-dimensional

applications to various problems, i.e., Hornberger et al. (1969), Taylor and Luthin (1969),

Verma and Brutsaert (1970), and Cooley (1971). These studies used the Laplace equation

in the saturated zone and are thus limited to near-surface flow in homogeneous,

incompressible, and unconfined aquifers. All these studies were for small regions with

simplified boundary configurations. In late 1960s, only Jeppson (1969) considered a

variably saturated flow regime on a basin wide scale but with the restriction of steady-

state condition.

In 1971, a remarkable three-dimensional transient, saturated-unsaturated flow

model was developed by Freeze (1971). The model was designed as a regional model

applicable to any groundwater basin. The complete subsurface regime was treated as a

unified whole by solving the variably saturated flow equation in the unsaturated zone and

the saturated flow equation in the underlying unconfined and confined aquifers. This









was the first complete three-dimensional hydrologic response model. Jacob's (1940)

equation, as clarified by Cooper (1966), as a saturated equation and the Richards (1931)

equation as an unsaturated equation were combined and solved in terms of the pressure

head in deforming coordinates to take into account the compressibility of the formation.

In the same year, Cooley (1971) developed a finite-difference method for unsteady flow

in variably saturated porous media. He applied his model to a single pumping well

successfully. According to Wise et al. (1994), the Freeze (1971) and Cooley (1971)

models are not robust because they incur numerical instabilities and convergence

difficulties.

In most applications, the pressure-based form of the variably saturated flow

equation is used. Celia et al. (1990) and Kirkland (1991) claimed that the numerical

solution of the pressure-based Richards equation has poor mass-balance properties in the

unsaturated zone.

Brutsaert (1971) developed a two-dimensional model by solving the mixed form,

i.e., the pressure head and moisture-content based, Richards equation, using the finite-

difference method. Neuman (1973) developed a numerical model (UNSAT2) similar to

Rubin (1968), but it was a finite-element model. Narasimhan et al. (1978) created

TRUST, which is pressure-based for variably saturated flow and moisture-content based

for unsaturated flow problems. An integrated finite-difference method with a mixed

explicit -implicit time stepping procedure was used in TRUST.

In the 1980s, research studies on groundwater and solute transport increased

greatly. In this period, many remarkable research papers were published and the most

well known groundwater models MODFLOW and FEMWATER were created. van









Genuchten (1980) presented his closed form equation describing the relationship between

the hydraulic conductivity and pressure head. Yeh and Ward (1980) developed the two-

dimensional finite element variably saturated code called FEMWATER, which was

updated as a three-dimensional finite element model by Yeh and Cheng (1994) as

3DFEMWATER. Cooley (1983) presented his two-dimensional variably saturated finite

element model. In that paper, a new method for locating the position of a seepage face

was presented. Huyakorn et al. (1984 and 1986) developed two- and three-dimensional

finite element variably saturated flow models, respectively. In those models, the lower

and upper (LU) decomposition method, Newton-Raphson method, Picard iteration

schemes, and strongly implicit procedures were used. Voss (1984) of the U.S. Geological

Survey developed a two-dimensional finite-element simulation model called SUTRA for

saturated-unsaturated, fluid density-dependent groundwater flow with energy transport.

van Genuchten and Nielsen (1985) modified the original van Genuchten (1980) formula

to allow a non-zero value of specific moisture capacity in the saturated zone, which

makes the formula very useful for variably saturated flow models. Allen and Murphy

(1986) presented a variably saturated model that solved the mixed form of the Richards

equation using the finite-element method and the Gauss elimination technique during

iterations.

Kuiper (1986) compared seventeen different methods for solution of the

simultaneous nonlinear finite-difference approximating equations for groundwater flow in

a water-table aquifer in three-dimensions. The best methods were found to be those using

Picard iteration implemented with the preconditioned conjugate gradient method.









Lappala et al. (1987) developed a computer model, VS2D, for solving problems of

variably saturated, single-phase flow in porous media in two dimensions. Nonlinear

boundary conditions treated by the model included infiltration, evaporation, seepage

faces, and water extraction by plant roots. Subsequently, Healy (1990), who was one of

the co-authors of VS2D, added a solute transport capability to the VS2D and named the

new model VS2D. In 1986, another comprehensive saturated unsaturated flow model

named SHE was developed in Europe by a multinational group consisting of the Danish

Hydraulic Institute, SOGREAH of France, and the British Institute of Hydrology (Abbott

et al., 1986 a, b). SHE was revised and developed as MIKE SHE in 1995 by Refsgaards

and Storm (1995) of the Danish Hydraulic Institute. MIKE SHE is a complete hydrologic

system model that simulates overland and channel flow, snowmelt, evapotranspiration,

and saturated-unsaturated flow.

In the 1990s, research interests were focused on developing more numerically

stable, fast converging, and more accurate numerical methods to solve complex,

nonlinear partial differential equations. Another area of interest was finding exact

analytical solutions to the nonlinear Richards equation. During this period, a large

number of papers were published dealing with mathematical solution techniques,

coupling of surface-water and groundwater techniques, and GIS application techniques in

groundwater hydrology.

Kirkland et al. (1992) presented a successful and efficient example of a finite-

difference solution to two-dimensional, variably saturated flow problems. However, the

objective of Kirkland et al. was the development of competitive numerical procedures to

solve infiltration problems in very dry soils. Thus, they did not take into account the









effects of specific storage in their fundamental flow equation, and, consequently, their

model cannot be used to model accurately a wide variety of variably saturated flow

problems (Clement et al., 1994).

Clement et al. (1994) presented an algorithm for modeling variably saturated flow

in the two-dimensional finite-difference form. The mixed form of Richards equation was

solved in finite-differences using a modified Picard iteration scheme to determine the

temporal derivative of water content. Wise et al. (1994) presented a sensitivity analysis

of the same variably saturated flow model to soil properties. They showed that the

location of the phreatic surface and height of the seepage face are functions of the

capillary forces exerted in the vadose zone.

Yan and Smith (1994) attempted to integrate the South Florida Water

Management Model (SFWMM) (MacVicar et al., 1983) and MODFLOW to simulate the

hydrogeologic system of south Florida. They presented the algorithm of the proposed

model that is conceptualized and constructed with a reasonable level of detail regarding

the simulation of surface-water movement, groundwater movement, and interactions

between the surface-water and groundwater systems. The movement of water outside of

the aquifer is simulated using SFWMM, and the water movement within the aquifer

system is simulated using MODFLOW. The models are linked by processes that include

recharge, infiltration, changes in soil moisture in the unsaturated zone, evapotranspiration

from the unsaturated and saturated zones, and flow between surface-water bodies and the

aquifer as recharge or discharge. No further action has been taken to create the real

model, and it has remained as a proposal (personal communication with Yan, August,

1996). Their evapotranspiration and infiltration formulation was implemented in MIKE









SHE as an alternative to the complex ET modules of MIKE SHE to use in South Florida

Water Management District projects (Yan et al., 1998).

Clement et al. (1996) compared modeling approaches for steady-state unconfined

flow. The Dupuit-Forchheimer, the fully saturated flow, and the variably saturated flow

models were compared for problems involving steady-state unconfined flow through

porous media. The variably saturated flow model was found to be the most

comprehensive of the three.

Parallel to those studies mentioned above, a vast amount of research exists

concerning evapotranspiration calculations. Although evapotranspiration is an important

element of the saturated-unsaturated flow model, the theory of evapotranspiration and

development of the evapotranspiration models are beyond the scope of this dissertation.

In this study, one of the most widely accepted and physically based evapotranspiration

models, i.e., the Priestly-Taylor (1972) model, was selected for use in the numerical

model.

Available Hydrologic Computer Models


In this section, widely used, numerical saturated-unsaturated groundwater flow

models that provided insight and guidance in the development of the model in this study

are discussed briefly (Table 2.1). MODFLOW is included in this discussion because of

its very well organized modular structure and multi-layer aquifer simulation capability in

the saturated zone, although it is a saturated flow model only.

MODFLOW. MODFLOW is a three-dimensional finite-difference ground-water

flow model (McDonald and Harbaugh, 1988). It has a modular structure that allows it to









be modified easily to adapt the code for a particular application. Many new capabilities

have been added to the original model of 1988.

MODFLOW simulates steady and unsteady flow in an irregularly shaped flow

system in which aquifer layers can be confined, unconfined, or a combination of confined

and unconfined. Flows from external stresses such as flow to wells, areal recharge,

evapotranspiration, flow to drains, and flow through river beds can be simulated.

Hydraulic conductivities or transmissivities for any layer may differ spatially and be

anisotropic. However, they are restricted to having the principal direction aligned with

the grid axes. The storage coefficient may be heterogeneous, and the model requires

input of the ratio of vertical hydraulic conductivity to distance between vertically adjacent

block centers (vertical conductance). Specified head and specified flux boundaries can be

simulated.

In MODFLOW, the ground-water flow equation is solved using the finite-

difference approximation. The flow region is considered to be subdivided into blocks in

which the medium properties are assumed to be uniform. In the vertical direction, zones

of varying thickness are transformed into a set of parallel "layers". The associated matrix

problem can be solved by choosing one of several solver routines that are available, i.e.,

the strongly implicit procedure, slice successive over relaxation method, and

preconditioned conjugate gradient method. Mass balances are computed for each time

step and as a cumulative volume from each source and type of discharge.

In order to use MODFLOW, initial conditions, hydraulic properties, and stresses

must be specified for every model cell in the finite-difference grid. The primary output is










hydraulic head. Other output includes the complete listing of all input data, drawdowns,

and water-budget data.



Table 2.1 Summary of selected saturated-unsaturated flow models

Title Developer Distinct features Limitations
McDonald and 3-D finite-difference, distributed, No unsaturated flow
MODFLOW Harbaugh, 1988 saturated groundwater flow model, modeling, requires user
Boussinesq equation is solved, supplied ET and recharge
values.
Yeh et al., 1996 3-D finite-element, distributed, No root water uptake
saturated-unsaturated flow and feature, requires user
FEMWATER transport model. Pressure based supplied net rainfall,
Richards equation is solved, infiltration capacity, and
ET.
Skaggs, 1980 Lumped model, approximate Not distributed, no 3-D
methods are used. Variably saturated, simulation, and no
DRAINMOD designed for drainage and irrigation confined aquifer or
problems pumping.
Voss, 1984 2-D finite-element, distributed, No 3-D simulation, lack of
variably saturated flow and transport ET and surface-water
SUTRA model. General variably saturated groundwater interaction.
density dependent flow equation is
solved.
Lappala et al., 2-D finite-difference, distributed flow No 3-D simulation, poor
1987 and Healy, and transport model. Pressure based surface-water groundwater
VS2DT 1990 form of general variably saturated interaction.
flow equation is solved.
Bloom et al., Modified VS2DT with coupled ET No 3-D simulation.
1995 calculations and surface-water
WETLANDS groundwater interaction modules. 2-
D Richards equation is solved.
Abbott et al., 3-D finite-difference, distributed, Unsaturated zone is 1-D
1986 a, b, and saturated-unsaturated model. ET vertical. Difficulties exist
MIKE SHE Refsgaard and calculations, surface-water in coupling 3-D saturated
Storm, 1995 interactions are considered. Richards and 1-D unsaturated zone
and Boussinesq equations are solved, modules.
FRAC3DVS Therrien, and 3-D finite-element and finite- No ET calculations or
Sudicky, 1996. difference variably saturated flow surface water groundwater
and transport model for fractured interaction.
porous media. Mixed form of the 3-
D Richards equation is solved.
HYDRUS Vogel et al., 1-D finite-element variably saturated No 2-D or 3-D simulations.
1996. flow and solute and heat transport
model. Mixed form of the 1-D
Richards equation is solved using a
new convergence criterion.









The main limitation of the model is that it lacks the capability to simulate flow in

the unsaturated zone. Although MODFLOW can simulate evapotranspiration, recharge,

areal recharge, and river-groundwater interaction, some of the model parameters have to

be provided by the user to compensate for the lack of unsaturated zone simulation. These

parameters includes the maximum evapotranspiration rate, net recharge to the water-

table, etc. Therefore, the model may give erroneous results because of the user defined

parameters.

FEMWATER. FEMWATER is a 3D flow and contaminant transport finite-

element density-driven coupled or uncoupled model used to simulate both saturated and

unsaturated conditions (Yeh et al., 1996). FEMWATER was formed by combining two

older models, 3DFEMWATER (flow) and 3DLEWASTE (transport), which were written

by Yeh and Cheng (1994), into a single coupled flow and transport model.

The governing equation for flow is the three-dimensional Richards equation

modified to include a density dependent flow term and to consider consolidation of the

aquifer. There are four types of iteration methods for solving the linearized matrix

equations of the governing equation, i.e., successive point and block iteration, polynomial

preconditioned conjugate, and incomplete Choleski preconditioned conjugate gradient

methods.

The model requires identification of material properties representing the

hydrogeologic and transport characteristics of soil contained within the model. The

moisture-content, relative conductivity, and water capacity versus pressure head curves

should be supplied to the model. Initial conditions and four kinds of boundary conditions,

namely Dirichlet, Neumann, Cauchy, and variable boundary conditions, can be assigned









by selecting nodes or element faces. Features such as wells, constant head, and no-flow

boundaries can be defined. Transient data (such as recharge or well pumping), which is

typically available in hydrograph form, can be input and edited graphically. These data

can then be interactively assigned to a single element or a series of elements. The model

outputs are pressure head, moisture-content, Darcy velocity, and concentration values at

each node.

FEMWATER has no capability to calculate the evapotranspiration losses, which

have to be supplied to the model separately. The model does not calculate the rainfall-

runoff process. Therefore, the net precipitation has to be supplied to the model as a

boundary condition. Vegetative cover, precipitation, and evaporation interactions are not

considered in the model. The evapotranspiration losses calculated outside of the model

can be applied only at the ground surface as a boundary condition, and transpiration

losses below the ground surface around the root zone are not considered.

DRAINMOD. DRAINMOD is a lumped hydrologic simulation model developed

by Skaggs (1980). The model simulates the hydrology of poorly drained, high water-table

soils on an hour-by-hour, day-by-day basis for long periods of the climatological record

(e.g., 40 years). The model predicts the effects of drainage and associated water

management practices on water-table depths, the soil-water regime, and crop yields.

Initially, DRAINMOD was used as a research tool to investigate the performance

of a broad range of drainage and subirrigation systems and their effects on water use, crop

response, land treatment of wastewater, and pollutant movement from agricultural fields.

The specific objectives of DRAINMOD are to simulate the performance of water-table

management systems and to simulate lateral and deep seepage from the field.









DRAINMOD was developed for field-sized units with parallel subsurface drains. Most

of the hydrologic components considered in the water balance are formulated in the

model. DRAINMOD uses approximate methods to characterize the rates of infiltration,

drainage, evapotranspiration, and the distribution of soil water in the profile instead of

using numerical solutions to nonlinear differential equations. However, the estimates

provided by this approximate method are comparable to exact methods. A general flow

chart for DRAINMOD is given by Skaggs (1980). The following are the model

components:

1. Precipitation (hourly data is suited);

2. Infiltration: the Green-Ampt equation is used to compute infiltration;

3. Surface drainage: the average depth of depression storage, which must be

satisfied before runoff;

4. Subsurface drainage: the rate of subsurface-water movement into drain tubes or

ditches;

5. Subirrigation;

6. Evapotranspiration;

7. Soil water distribution; and

8. Rooting depth.

The input data can be summarized as follows: soil property inputs; hydraulic

conductivity; soil-water characteristics; drainage volume water-table depth relationship;

upward flux; Green-Ampt equation parameters; trafficability parameters; crop input data;

drainage system parameters; surface drainage; and effective drain radius. Outputs are:

yearly rankings of parameters such as number of working days and relative yields; daily









and monthly summaries of many of the output parameters; relative yield; and daily water-

table depths and drainage volumes for each year simulated.

Sensitivity analyses were conducted for different soils and water management

systems of North Carolina. The results are presented in the DRAINMOD reference report

(Skaggs, 1980). These results indicate that errors in the hydraulic conductivity (K) have

the greatest effect on predicting the water-table depth and water content of the soil

profile. DRAINMOD model was tested for use in humid and semi-arid climatic regions.

Skaggs (1980) reported the following limitations of this model:

1. DRAINMOD should not be applied to situations that are widely different from

conditions for which it was developed, without further testing; and

2. The field should have parallel subsurface drains.

In addition to those limitations, the model is not distributed so that heterogeneity

in the field can not be simulated. DRAINMOD uses approximate methods, which

requires extensive efforts to find appropriate coefficients for different geological

formations. The model was not designed to model confined aquifers, and lateral

movement of moisture in the unsaturated zone was not considered in the model.

SUTRA. SUTRA is a two-dimensional finite-element simulation model for

saturated-unsaturated, fluid-density-dependent ground-water flow with energy transport

or chemically-reactive single-species solute transport (Voss, 1984). SUTRA can be used

for areal and cross-sectional modeling of saturated ground-water flow systems and for

cross-sectional modeling of the unsaturated flow zone. SUTRA can also be used to

simulate solute transport to model natural or man-induced chemical species transport

including processes of solute sorption, production, and decay, and it may be applied to









analyze ground-water contaminant transport problems and aquifer restoration designs. In

addition, solute transport simulation with SUTRA may be used for cross-sectional

modeling of saltwater intrusion in aquifers in near-well or regional scales.

The model employs a two-dimensional hybrid finite-element and integrated-finite-

difference method to approximate the governing equations. These equations describe the

two interdependent processes that are simulated: (1) fluid density-dependent saturated or

unsaturated ground-water flow, and (2) transport of a solute in the ground water and solid

matrix of the aquifer.

Important limitations of the program are: it is two-dimensional, which is not

convenient to regional modeling of aquifers; and it does not have any capability to

calculate evapotranspiration losses and overland flow. The model's main emphasis is

solute transport and variable density flow, and it is not designed to simulate the complete

hydrological system such as extensive pumping from multiple aquifers that are highly

interactive with a river system and overland flow.

VS2DT. VS2DT solves problems of water and solute movement in variably

saturated porous media. The origin of the VS2DT is VS2D developed by Lappala et al.

(1987). Healy (1990) added a transport module and renamed the model as VS2DT. The

finite-difference method is used to approximate the flow equation, which is developed by

combining the law of conservation of fluid mass with the Darcy-Buckingham equation

and the advection-dispersion equation. The model can analyze problems in one and two

dimensions with planar or cylindrical geometries. There are several options for using

boundary conditions that are specific to flow under unsaturated conditions. There are

infiltration with ponding, evaporation, plant transpiration, and seepage faces. Solute









transport options include first-order decay, adsorption, and ion exchange. Extensive use

of subroutines and function subprograms provides a modular code that can be easily

modified for particular applications.

For the flow equation, spatial derivatives are approximated by the central-

difference method. Time derivatives are approximated by a fully implicit backward-

difference scheme. Nonlinear conductance terms, boundary conditions, and sink terms

are linearized implicitly using previous iteration step values. The relative hydraulic

conductivity is evaluated at cell boundaries by using full upstream weighting, the

arithmetic mean, or the geometric mean of values from adjacent cells. Saturated hydraulic

conductivities are evaluated at cell boundaries by using distance-weighted harmonic

means. Nonlinear conductance and storage terms can be represented by algebraic

equations or by tabular data.

The model requires initial conditions be input in terms of pressure heads or

moisture-contents for flow simulations and concentrations for transport simulations.

Hydraulic and transport properties of the porous media are also required. Flow

simulations require values for saturated hydraulic conductivity and for relative hydraulic

conductivity and moisture-content as functions of pressure head. Transport simulations

require values for dispersivity and molecular diffusion. Simulation results can be output

in terms of pressure head, total head, volumetric moisture-content, velocities, and solute

concentrations.

The main shortcomings of the model are that it can not simulate three-

dimensional problems. Also, it does not consider stream flow-groundwater interaction, it









does not consider interception losses from rainfall, and it cannot simulate multi-aquifer

systems.

WETLANDS. WETLANDS is a mathematical model for one- or two-

dimensional water flow and solute movement in variably- saturated multi-layered porous

media featuring optional surface-water bodies (ponds) and multiple root zones (Bloom et

al., 1995). The Richards equation and the advective-dispersive equation are solved

numerically using a finite-difference approximation, and the interaction of water levels in

the ponds with the surrounding soils are continually and dynamically adjusted.

A Priestly-Taylor model is used to simulate evapotranspiration by one to three

plant species. The controlling parameters can be specified to track seasonal variation in

sunlight, temperature, and rainfall. Solute transport can be affected by plant uptake,

passive sinks, or a variety of sorption phenomena as well as by water transport in the soil

and ponds.

WETLANDS is a finite-difference model using the strongly implicit method to

simulate water and solute transport. WETLANDS is a descendant of VS2DT, but it has

been modified to support simulations of shallow systems featuring ponds (or lakes, rivers,

etc.) with coupled evapotranspiration of multiple plant species. Output can consist of

matrices of head and solute concentration, temporal traces of stream flow and solute

concentration, mass balance monitoring, and data-sets depicting the water-table at any

given time.

The limitations of WETLANDS are that it is a two-dimensional model, and it

cannot simulate multi-aquifer systems or sink/sources such as pumping, drains, and

springs.









MIKE SHE. MIKE SHE is a distributed, physically based, three-dimensional,

finite-difference saturated unsaturated hydrological system model (Refsgaard and Storm,

1995). The MIKE SHE model was derived from the SHE model (Abbott et al., 1986a

and 1986 b). The model is applicable to a wide range of water resources problems related

to surface-water and ground-water management, contamination, and soil erosion. It is

designed as a modular structured model so that it can be easily modified or expanded. It

has the following components: evapotranspiration (ET), unsaturated zone flow (UZ),

saturated zone flow (SZ), overland and channel flow (OC), and irrigation module (IR)

components.

Different time scales can be used for different flow processes throughout the

simulation. For example, smaller time steps can be used in the unsaturated zone than are

used in the saturated zone. However, the UZ, ET, and OC modules use identical time

scales. This feature saves computer memory and enables faster simulations.

The UZ module plays an important role in the MIKE SHE, because all the other

components depend on boundary data from the UZ module. The flow is one-dimensional

vertically in UZ module. The governing equation for flow is the one-dimensional form of

the Richards equation. The UZ includes root extraction for the transpiration process,

which is explicitly incorporated in the equation by sink terms. The integral of the sinks

over the entire root zone depth gives the total amount of actual evapotranspiration. If the

root zone is homogeneous in certain regions, then the UZ calculations are only performed

in one representative column within those regions, and then lumped together for each

homogeneous region. The relationship between the moisture-content and pressure head

and hydraulic conductivity is a necessary input to the UZ module.









The ET module uses meteorological and vegetative input data to predict the total

evapotranspiration and net rainfall amounts. In the calculation of net rainfall amounts, the

processes of interception by the canopy, drainage from the canopy, evaporation from the

canopy surface, evaporation from the soil surface, and plant root water-uptake are

considered. An evapotranspiration model developed by Kristensen and Jensen (1975) is

used in the ET module.

The OC model is designed to route the excess ponded water on the ground-surface

towards the river system. The exact route and quantity is determined by the topography

and flow resistance as well as the losses due to evaporation and infiltration along the flow

path. Both the overland flow and the channel flow are modeled by approximations of the

St. Venant equations.

The SZ component of MIKE SHE calculates the saturated subsurface flow in the

catchment by solving the quasi-three-dimensional Boussinesq equation. MIKE SHE

allows for a fully three-dimensional flow in a heterogeneous aquifer with shifting

condition between unconfined and confined conditions. The spatial and temporal

variations of the hydraulic head are described by the nonlinear Boussinesq (1868)

equation and solved numerically by an iterative finite-difference technique. Successive

over relaxation and the preconditioned conjugate gradient solution techniques are

available in the model. In structure and flexibility, the SZ module is similar to

MODFLOW.

There is a difficulty in the linkage between the unsaturated zone and the saturated

zone in MIKE SHE because the UZ and SZ components run parallel, and thus they are

not solved in an integrated form. This difficulty has been solved by using an iterative









procedure based on mass balance calculations for the entire column including horizontal

flows in the saturated model. Because of the one-dimensional structure of the UZ

module, horizontal moisture movement in the unsaturated zone cannot be simulated in

the model. The model uses two different governing equations in the unsaturated and

saturated zones although one equation, the Richards equation, could be used for both

zones.

GMS (Groundwater Modeling System). GMS is not a groundwater model but,

instead, it is a groundwater modeling environment developed by the Department of

Defense (Yeh et al., 1996). GMS integrates and simplifies the process of groundwater

flow and transport modeling by seamlessly integrating all the tools needed for a

successful study. GMS supports the following models: MODFLOW, MODPATH,

MT3D, FEMWATER, SEEP2D, and RT3D.



FRAC3DVS. FRAC3DVS was developed by Therrien and Sudicky (1996). It is a

three-dimensional finite-element model that simulates saturated-unsaturated groundwater

flow and solute heat transport in porous or discretely fractured porous media. Galerkin

finite-element or finite-difference schemes can be selected. A conjugate-gradient like

solver is used to solve the systems of equations, and a full Newton-Raphson iteration

scheme is used to linearize the non-linear mixed form of the Richards equation.

HYDRUS. HYDRUS is a one-dimensional variably saturated groundwater flow

and solute and heat transport model developed by Vogel et al. (1996). It solves the mixed

form of the Richards equation using a new convergence criterion (Huang et al., 1996) to

speed up the iterative solution process. The model allows hysteresis to occur in both the






42


soil-water retention and hydraulic conductivity functions. It can simulate root-water

uptake, and it also includes heat transfer and heat movement simulations.















CHAPTER 3
DERIVATION OF THE VARIABLY SATURATED GROUNDWATER FLOW
EQUATION


General Three-Dimensional Saturated-Unsaturated Groundwater Flow Equation


Conceptualization

A general three-dimensional saturated-unsaturated groundwater flow equation can

be derived by considering the hydrological events and parameters that are depicted in

Figure 3.1. Using the concept of conservation of mass in a hydrological system such as

the one shown in Figure 3.1, the governing equation for groundwater flow can be derived.


Figure 3.1. Conceptualization of hydrologic system.









The conservation of mass concept considers that the sum of the inputs to the

system minus the sum of the outputs from the system is equal to change of mass in the

system per unit time. The mathematical description of the conservation of mass can be

described in terms of a unit volume taken from an interior location in the groundwater

system (Figure 3.2).






(qz)out


(q,)out



Ax
-- (qx)out
(qx)in Az (q)out




(qy)in




(qz)in


Figure 3.2 Representative unit volume of an aquifer.









Continuity Equation

The general groundwater flow equation is developed based on the mass continuity

(mass conservation) equation. The mass continuity equation can be written for a unit

volume of an aquifer (Figure 3.2) as




I- W= dS (3.1)
dt



where I is the mass inflow rate in the x, y, and z directions [MT-1], 0 is the mass outflow

rate in the x, y, and z directions [MT-1], W is a sink/source term representing the mass of

water injected into or removed from the aquifer per unit time [MT-1], and dS/dt is the

change in mass storage (S) per unit time [MT-1].

The mass inflow rates at x, y, and z in the x, y, and z-directions respectively are



I = (pqx)dz dy (3.2a)

Iy = (p qy)dx dz (3.2b)

Iz = (pqz)dx dy (3.2c)



where p is the fluid (water) density[ML-3], and q is the specific discharge, i.e., the Darcy

flux [LT-1].

Similarly, the mass outflow rates at x+dx, y+dy, and z+dz in the x, y, and z

directions (approximately by Taylor series expansion) respectively are









Ox = (p qx)dz dy +-(p q)dz dy dx
Ox

O = (p qy)dx dz + (p qy)dx dz dy
Oyy

Oz = ( qz)dx dy + -(p qz)dx dy dz
Oz



where dx dy dz is the volume of the unit representative element.



The right hand side of equation (3.1) can be written as



dS 8
dt -[(p 0) dx dy dz]
dt at


(3.3a)


(3.3b)


(3.3c)













(3.4)


where 0 is the volumetric moisture-content of the medium[L3L-3].

If equations (3.2a),(3.2b),(3.3c), and (3.4) are substituted into equation (3.1) and

rearranged, then,


L(P qx) + (pqy) + (pqz) +
ax ay z


(p 0)
Bt


(3.5)


w
where w = mass of water injected or removed from a unit volume of the aquifer
dx dy dz

per unit time [ML-3T-1].









Equation (3.5) is the continuity equation. There is only one equation and there are

five unknowns, i.e., the three components of the Darcy flux (qx, qy, qz), and 0 and p.

Thus, it is necessary to formulate four more equations to solve equation (3.5) for the five

unknowns.

Storage Term

The storage term on the right-hand side of the equation (3.5) can be expanded by

defining Sw as the saturated fraction of the porous medium and substituting this parameter

into the storage term in equation (3.5). This yields




(Sw) S+ p S +n qS (3.6)
at at at at




where r" is the porosity and Sw = -


Using the chain rule for differentiation, equation (3.6) can be rewritten in terms of

the pressure head h = p/y [L] as



O(p Sw) BSw Oh 8O Oh Bp Oh
t = pn +pSw+W +, SW (3.7)
8t Oh 8t Oh 8t Oh 8t



The first term of equation (3.7) accounts for the change in fluid storage due to a

change in the volumetric water content. It actually describes the effects of draining and

filling the pores. The first term can be redefined as









=- -=C(h) (3.8)
Oh Oh



where C (h) [L-1] is the slope of the moisture retention curve. This term is called the

specific moisture capacity, and it expresses the volume of water released per unit volume

of unsaturated zone for a unit decrease in pressure head h.

The second term in equation (3.7) accounts for the change in fluid storage due to

the compressibility of the solid matrix:




apg (3.9)
Oh



where ca is the solid matrix compressibility [LT2M-1], and g is the acceleration of gravity

[LT-2].

The third term accounts for the change in fluid storage due to fluid

compressibility:



=P p2g (3.10)
Oh



where P3 is fluid compressibility [LT2M'-]. Equations (3.8) through (3.10) can be

combined and substituted into equation (3.7) to give









=(po) p[C(h) + SSs]h (3.11)
dt dt



where Ss is the specific storage [L-1] defined as



SS = p g(Uc + q 3) (3.12)



The specific storage represents the volume of water released per unit volume of

aquifer per unit decline in pressure head. Equation (3.11) is the second equation of the

five equations necessary to solve for the five unknowns in Equation (3.5). This equation

is based on the assumptions that the aquifer and water are slightly compressible in the

saturated confined zone but incompressible in the unsaturated zone and in the unconfined

saturated zone. Therefore, Ss approaches zero in the unsaturated and the unconfined

aquifer zones because of its dependency on the compressibility of the solid matrix and

fluid. On the other hand, although C(h) can have significant values in the unsaturated

zone, it has very small values approaching zero in the saturated zone. This is because

C(h) is the slope of the moisture retention curve, which is zero in the saturated zone, i.e.,

the moisture-content is constant in the saturated zone.

Darcy-Buckingham Equation

Darcy developed his well-known formula for saturated flow conditions, and

Buckingham developed nearly the same relationship for unsaturated flow conditions.

Combining those two formulas results in the Darcy-Buckingham equation, which is used

as the flux equation for both saturated and unsaturated zones (Narasimhan, 1998b).









In Darcy's equation, the flux is linearly proportional to the hydraulic gradient, and

the proportionality constant is defined as the hydraulic conductivity. In Buckingham's

equation, in the unsaturated zone, the proportionality constant is not linear and the

hydraulic conductivity is a function of both the pressure head and the medium properties

of the unsaturated zone.

Defining the hydraulic head (or total head), H, as



H = h+z (3.13)



where h is the pressure head [L] and z is the elevation head [L], the Darcy-Buckingham

equation can be written as



OH
q = -K (h) H1 (3.14)




where q is the specific discharge [LT-1] in the x, y, and z directions, Kij (h) is the

hydraulic conductivity [LT-1] in the x, y, and z directions, and 1 represents the unit

distances in the x, y, and z directions.

Hydraulic conductivity is not only a function of the porous medium but also of the

fluid properties. Hubert (1956) pointed out that hydraulic conductivity is directly

proportional to the square of the mean grain size diameter (d2) and the specific weight of

the fluid (pg) and inversely proportional to the fluid viscosity ([t) (Bear, 1972). Together

with Darcy's original observation and dimensional analysis, the hydraulic conductivity









can be expressed as K = Cd2pg/[t. The term Cd2 is a property of the soil itself, and it is

called the intrinsic permeability, k. The coefficient C in the intrinsic permeability (k)

represents the grain-size distribution, the sphericity and roundness of the grains, and the

nature of their packing. The hydraulic conductivity K is written as K = kpg/[t.

K(h) is function of the pressure head in the unsaturated zone, but it is constant and

equal to saturated hydraulic conductivity in the saturated zone, i.e., Kij (h)= Ks. Some

typical values of Ks can be found on page 29 in Freeze and Cherry (1979).

In general, in a three-dimensional flow field, the hydraulic conductivity tensor

could have nine components. However, the hydraulic conductivity tensor is symmetric

such that Kxy(h) = Kyx (h), Kx (h) = Kz (h), and Kyz(h) = Kzy(h), and thus it has only six

components:




Kxx K Kxz
K KX KX Kyy (3.15)
KXZ Kyz Kzz



If the principle axis of anisotropy is aligned with the principle axis of flow, then

only three non-zero hydraulic conductivity terms remain, i.e., Kxx (h), Kyy (h), and Kzz (h).

Thus, the hydraulic conductivity tensor becomes





K = 0 K 0 (3.16)
0 0 KZZ









Governing Equation (Modified Richards' Equation)

Substituting equation (3.16) into the Darcy-Buckingham equation (3.14) in the x,

y, and z directions results in:


qx =-Kx(h) O
Ox

q, = -Ky(h) aO
ay

aH
qz =-Kz(h) O
Oz


ah
-Kx(h)
Ox

-K (h)
By

-Kz(h)(- + 1)
Oz


(3.17a)


(3.17b)


(3.17c)


These equations (3.17 a-c) are the third, fourth and fifth equations of the five

equations required in order to solve the five unknowns of equation (3.5).

Assuming that water density does not vary spatially, such that-- = 0, and
ax,

substituting the Darcy-Buckingham equation (equation (3.17)) and the general saturated-

unsaturated storage term (equation (3.11)) into the continuity equation (equation (3.5))

gives the following form of the three-dimensional Richards equation:




O(Kx(h)( O)) +(Ky(h)( )) +(Kz(h)(h + 1)) h
F+ +--- Qem =[C(h)+Ss] (3.18)
Ox + y + hz Qe Bt
Ox O Ozxt Ch) S,,s ]a t(3.18)









where Qext is a volumetric source or sink term, which is obtained by dividing w by the

density of water, Qe~ = [L3L 3T 1].
p

Equation (3.18) is the general three-dimensional saturated-unsaturated flow

equation that is called the "modified Richards equation" due to the inclusion of the

saturated zone, which is achieved by modifying the general storage term on the right hand

side of equation (3.5).

The expressions for C(h) and K(h) are both highly nonlinear, which makes the

solution of the governing equation very difficult and complex.

In the saturated zone:

C(h) = 0;

K(h) = Ks = constant;

Sw = 1; and

h hair entry .

Therefore, in the saturated zone, equation (3.18) becomes




O(K() )) -(Ky( )) (Kz( + 1)) Sh
S+ + Oz +-5 -S a h (3.19)
ox ay az t




In the unsaturated zone:

C(h) a 0;

C(h) >> S, Ss;









Swr = --<1;


h
K(h) = function of the pressure head.

Therefore, in the unsaturated zone, equation (3.18) becomes




O(Kx(h)( )) a(Ky(h)( )) I(Kz(h)( Oh+l))
x + y + z + = C(h) a- (3.20)
ox Sy 8z dt




The right side of the equation (3.20) describes the effects of draining and filling

the pores in the unsaturated region. Thus, expressing this concept in terms of the

temporal change in moisture-content would be more appropriate than expressing it in

terms of pressure head. In other words, the term C(h)(ah/Ot) = (dO/dh)(ah/Ot) can be

written more appropriately in its original simpler form, i.e., 6O/8t.

Using the term 80/&t in equation (3.20) converts the pressure-based modified

Richards equation into a mixed form of the modified Richards equation. Celia et al.

(1990) showed that the modified Picard iterative procedure for the mixed form of the

Richards equation is fully mass conserving in the unsaturated zone. By contrast, the

conventional pressure-based, backward Euler finite-difference formulations exhibit poor

mass-balance behavior according to Clement et al. (1994) and Celia et al. (1990). The

reason for this is that the discrete analogs of 80/8t and C(h) ah/8t are not equivalent even

though the time derivative of the moisture-content, 80/0t, is equal to C(h)ah/at, which is a









mathematically valid approximation (Clement et al., 1994). This inequality is amplified

owing to the highly nonlinear nature of the specific capacity term, C(h). Using the

modified Picard iteration method eliminates this problem by approximating directly the

temporal term 80/8t with its algebraic analog (Clement et al., 1994). The algebraic

approximation of the temporal term (09/at) and the modified Picard iteration method are

described in detail in Chapter 4 of this study.




Hydraulic Conductivity


The hydraulic conductivity, K, is constant with respect to time and equal to the

saturated hydraulic conductivity, Ks, in the saturated zone. In this study, a relative

hydraulic conductivity term, Kr, is used. Kr is the ratio of the unsaturated hydraulic

conductivity to the saturated hydraulic conductivity, K(h)/ Ks and thus K = Kr Ks.

The hydraulic conductivity in the unsaturated zone is defined as a function of the

pressure head, which can be derived from moisture-retention (or moisture characteristic)

curves, h versus 0. Several researchers have developed relationships between K(h) and

moisture retention curves. Measuring pressure head and moisture-content, h versus 0, is

easier than measuring pressure head versus K(h), so therefore h versus 0 relationships are

very useful in determining unsaturated hydraulic conductivity values.

The specific moisture-content C(h) is defined as the slope of the moisture

retention curve. It can be found by taking the derivative of the moisture-content with

respect to pressure head, h, or









C(h) =-- (3.21)
dh



Three different moisture-content-pressure head-hydraulic conductivity algebraic

relationships are used in this model study, i.e., the Brooks and Corey (1964) equations,

the van Genuchten and Nielsen (1985) relations, and the general power formula.

Brooks and Corey method. The Brooks and Corey (1964) equations are




Se = -Or = when h < ha, and (3.22)
Os -Or h


Se = --r 1 when h > ha (3.23)




where Se is the effective saturation, Or is the residual water content, and Os is the saturated

moisture-content, which is generally equal to the porosity of the formation (qr), and k is a

pore size distribution index that is a function of soil texture. The term ha is the bubbling

(or air entry) pressure head, equal to the pressure head required to desaturate the largest

pores in the medium, and it generally is less than zero.

The hydraulic conductivity is defined as




Kr K(h) h when h < ha; and (3.23)
Ksat ha

Kr = 1 when h > ha (3.24)









The specific moisture capacity C(h) can be calculated from


C(h) -(n Or) h k h


C(h) = 0 when h > ha


when h < ha, and


(3.25)


(3.26)


van Genuchten and Nielsen method. van Genuchten and Nielsen (1985)

developed a closed-form equation for hydraulic conductivity as a function of the pressure

head using the moisture retention curve:


for h < 0, and


(3.27)


(3.28)


K, = K(h) ( +3 5m/21+P)m -P 2

K,


where h3 = ha is air entry (or bubbling) pressure head[L], and n is a fitting
ha

parameter in the moisture retention curve, or m =1-1/n.

This closed form equation can be obtained by applying the fitting curve technique

to measured, or experimental, moisture-content-pressure head data. The moisture-

pressure head data generally fit the following equations:


Se = Or (1+ Op)- n
Os Or


if h < 0, and


(3.29)










Se = -r 1
0s Or


if h>0


(3.30)


where Or is the residual water content, and Os is the saturated moisture-content, which

generally equals the porosity of the formation (qr).

For the moisture-content relations, Paniconi et al. (1991) modified van Genuchten

and Nielsen's relation (equation (3.28) and (3.29)) in the form


0(h) = Or + (Os Or)(1 + )-m for h < h0, and

O(h) =Or+(s Or)(+0)-m+Ss(h -ho) for h>ho


(3.31)

(3.32)


where Ss is the value of specific storage for the pressure head h that is greater than the air

(h0Yn
entry pressure, Po = ho and ho is a parameter determined on the basis of continuity
requirement s imposed on S

requirements imposed on Ss, which implies that


s =(n- 1)(s r)h n-1
hs= ( + )m+1 h=ho


For a given value of Ss, equation (3.33) can then be solved for ho.

The specific moisture capacity C(h) can be calculated from


(3.33)


n-I

hs n(1+ p)m+1


when h < h0, and


(3.34)









C(h) = 0 when h > ho. (3.35)



For Ss = 0 and ho= 0, equations (3.29-35) revert to their original form in van

Genuchten and Nielsen (1985).

General power formula. A general power formula also can be used if there is a

moisture-content-pressure head (h-O) data set. The hydraulic conductivity K can be

described as a function of the effective saturation, Se:



K(h) n
Kr- K Se (3.36)
Ks




where Se - or and n is a parameter that has to be estimated by calibration. As a
0s -Or

guideline, the exponent n is usually relatively small for sandy soils (between 2 to 5) and

larger for clayey soils (between 10 to 20). The value of n influences the percolation rate in

the soil and thereby influences the actual evaporation rate.

In this method, any kind of soil moisture retention curve (0-h) can be used, but the

data should be supplied in a tabular form to the model. The model calculates the

intermediate values using interpolation methods.

The specific moisture capacity can be calculated using tabular values of the soil

moisture retention curve with the following formula:




C(hm+/2) = C- m+1 m (3.37)
( h h m+1 -hm










Sink/Source Term


The volumetric sink/source term (Qext) [L3T'1L-3] in equation (3.18) represents

the volume of water removed or injected per unit time from a unit volume of soil due to

sinks such as root water uptake in the unsaturated zone and pumping from wells and flow

from drains in the saturated zone. It is a source term in case of artificial recharge or

injection. Qext can be expressed as Qext = Wr + Ww + Wd where Wr represents root water

uptake, Ww represents well recharge or discharge, and Wd represents a drain. The

components of the sink/source terms of this study are briefly described in figure 3.3.


Source/sink terms in the model


fI
Pumping or Evapotranspiration Drains and
Injection Springs


Inputs Calculate PET from Input the drain head
flowrates and Priestly-Taylor eq'n or and conductance
application nodes from pan evaporation between drain and
application nodes values(see figure 3.5). surrounding cells.


Divide specified ei Separate Ep and Tp Calculate flow to
flow rates from PET as PET = Ep drains and divide
volume of the cell. + Tp based on LAI them by volume of
(equation 3.68). the cells.
I
I I
Calculate actual Ea, Calculate actual
and apply it to the transpiration (root
land surface as a water uptake) and
negative flux distribute it to the
boundary condition, root zone (see
figure 3.4).

Fiu(ec dModified Richards Equation n--



Figure 3.3 Flow chart describing the principle sink/source terms in the model.









Although evaporation and rainfall can be considered as sink and source terms,

respectively, they are treated in this study as upper boundary conditions. This is because

these processes occur on the land surface. Rainfall is applied to the land surface as a flux

boundary condition, and evaporation is separated from evapotranspiration by a procedure

that is described in the following sections. Then, it is also treated as a flux boundary

condition on the land surface. It could be in the form of soil evaporation or direct

evaporation from surface-water bodies such as lakes and rivers. The transpiration, or root

water uptake, is considered a sink term and applied to the cell nodes in the unsaturated

zone where it is occupied by roots.

The main part of the sink/source term in the unsaturated zone is the transpiration

(or root water uptake) process. To calculate the actual transpiration, the following steps

are taken (see figure 3.4):

1. Potential evapotranspiration (PET) is determined using one of the two

methods:

a. Pan evaporation method; or

b. Physically based equations (i.e., Priestly and Taylor Equation).

2. Potential evaporation (Ep) and potential transpiration (Tp) are determined from

PET as follows:

a. PET = Ep + Tp;

b. Ep = PET*exp(-0.4 LAI) where LAI = leaf area index; and

c. Tp = PET Ep

3. The actual transpiration (TA) is determined using two different options:

a. The method of Feddes et al. (1978); or










b.The method of Lappala et al. (1987), i.e., VS2D.

4. The actual evaporation (EA) is determined using the method of Lappala et al.

(1987). The actual evaporation is not treated as a sink term but it is

considered in the top boundary condition as a negative flux, and the modified

Richards equation takes EA into account at the top boundary.

The procedure outlined above is briefly summarized in figure 3.4, and it is also

described in detail in the following sections.



Actual Transpiration, Ta or
root water uptake calculations:


Tp is calculated from
Tp = PET(1-exp(-0.4 LAI)).



The root distribution
and root growth functionWp(z,t)
are determined from equations
(3.75) and (3.76).



The water stress response
function a3- (h) is determined
from equation (3.77) (see figure
3.6).


The actual transpiration is calculated
by multiplying Wp(z,t) and ar(h).
Equation 3.78 is the final form of
root water uptake function.


< ModifiedRichards Equation:


Figure 3.4 Flow chart for actual transpiration calculations in the model.









Determination of Evapotranspiration

Evapotranspiration in this study is considered to be the combination of

transpiration and evaporation. Potential evapotranspiration (PET) calculations can be

carried out using two different options. The first option is the easiest and the most

empirical one, i.e., the pan evaporation technique. If there are not enough available

climatologic data to calculate PET using one of the physically based equations, then the

pan evaporation method can be used. The pan evaporation method requires daily

measured pan evaporation values and a pan coefficient. The PET is calculated as



PET = Cpan*Epan (3.38)



where Cpan is pan coefficient, which is generally equal to approximately 0.7, and Epan is

the measured pan evaporation in [L/T].



The second method of PET calculations is physically based, i.e., based on the

energy conservation method. Although there are many equations in the literature to

calculate PET, the Priestly-Taylor equation was selected for use in this study. The

Priestly and Taylor (1972) equation is a derivative of the Penman equation (Fares, 1996).

It is advantageous among the others because it requires the least amount of input data.

Despite the empirical nature of its proportionality factor cU, the Priestly-Taylor equation is

based upon physical theory, and it reduces input data requirements (Buttler and Riha,

1989). It is also a simplified form of the Penman equation and is most reliable in humid

climates where an aerodynamic component has been deleted and the energy term









multiplied by a constant ca (Jensen et al., 1989). The Priestly-Taylor equation can be

written as




PET =ca (Rn -G) (3.39)
k(A+y)



where PET is potential evapotranspiration in mm t-1 (t is time unit), Rn is net solar

radiation [W m-2], A is the gradient of saturation vapor pressure-temperature curve

evaluated at the air temperature Ta, X is latent heat of vaporization [J m-3], G is soil heat

flux [W m-2], y is the psychrometric constant (0.067 kPa C-1) and ca is an empirical

parameter that depends on the nature of the surface, the air temperature, and time of day

and which varies from 1.05 to 1.38 (Viswanadham et al., 1991). Values of a are

generally between 0.6 and 1.1, according to Spittlehouse and Black (1981). Priestly and

Taylor (1972) obtained a mean value of ca equal to 1.26 for an extensive wet surface in

the absence of advection. Jury and Tanner (1975) showed that ca increases with heat

advection from surrounding areas and suggested a procedure for adapting the Priestly-

Taylor equation to such conditions (Fares, 1996).

The Priestly-Taylor equation requires values for six input parameters: the Priestly-

Taylor coefficient (c), slope of the saturation vapor pressure curve for water (A), net

radiation (Rn), soil heat flux (G), psychrometric constant (y), and latent heat of

vaporization (k). Fares (1996) developed an analytical model to calculate the above

parameters using available meteorological data, i.e., maximum and minimum daily

temperatures for a given location, day of the year, altitude and latitude of the location, and









albedo coefficient of the surface. If the albedo coefficient is equal to 0.05, the calculated

value of PET will be the potential evaporation directly from a free water surface.

However, if the albedo coefficient is equal to 0.15 (i.e., for pine trees), the result will be

the PET for pine trees. The flow chart of the calculations of PET using the Priestly-

Taylor equation is presented in figure 3.5.


Figure 3.5 Flow chart for the evapotranspiration calculations (Fares, 1996).









Estimation of input parameters for PET calculations

All the time units in PET calculations are in day. If it is desired, during the model

simulation daily values of PET can be converted into hourly values by assuming that

there is a sinusoidal distribution of the PET process based on a daily cycle in which PET

reaches its maximum value around noon time and reaches its minimum value around

midnight.

Calculation of the daily net radiation (Rn). Daily total values of Rn (MJ m-2 t-1)

can be determined from daily total incoming solar radiation, Rs (MJ m-2 t-1), and outgoing

thermal or long-wave radiation if they are not measured in the field. The following

relationship was proposed by Penman (1948) and modified by Wright(1982) to estimate

Rn:




R,, = (1 J)R, axk Tmmk l(a 0.139e aea + b (3.40)




where P3 is albedo (or reflectivity) coefficient of the surface, Rs is daily total incoming

solar radiation, c is the Stephan-Boltzman constant (4.903x10-9 MJ m-2 d-1 K-4), Tmak and

Tmink are the maximum and minimum daily air temperatures (K), respectively, ed is

saturation vapor pressure (kPa) at the dewpoint temperature (which is taken as the

minimum daily temperature), Rso (MJ m-2 d-') is daily total clear sky short wave radiation,

and a,, a, b are empirical coefficients.









Accepted values for the albedo coefficient (P3) are 0.05 for water surfaces; a range

of 0.15 to 0.60 for bare soil surfaces; 0.25 for most agricultural crops; and 0.1 for forests

(Fares, 1996). Wright (1982) estimated the a, empirical coefficient using



a1 = 0.26 + 0.1 e(00154 (J180))2 (3.41)



where J is the day of the year (1-365). If there are few clouds (Rs/Rso > 0.7), then a =

1.126 and b = -0.07; if it is cloudy (Rs/Rso< 0.7), then a = 1.017 and b = -0.06. Clear sky

short wave radiation (Rso) can be estimated from the Jensen et al. (1989) relationship:



Rso = 0.75 RA (3.42)



where RA is extraterrestrial radiation (MJ m-2 d-'). Although solar radiation (Rs) can be

measured using sophisticated instruments, the estimation of Rs using RA is also possible

using the Doorenbos and Pruitt (1977) equation:




Rs = 0.25+0.5N RA (3.43)




where n is the number of actual bright sunshine hours and N is the maximum possible

sunshine hours for that location. RA can be calculated using the following set of

equations developed by Duffied and Beckman (1980):









RA = 37.58 dr {ws sin(4) sin(6) + cos(l) cos(6) sin(ws)} (3.44)



where 4) and 1 are the longitude and latitude of the location in radians respectively (-E,

+W, -S, +N), dr is the relative distance of the earth from the sun, 6 is declination in

radians, and Ws is sunset hour angle in radians, which can be calculated as




dr = 1+ 0.033 cos 27T5} (3.45)


6 = 0.4093 sin27 284 365 (3.46)
365

ws = Arc cos(- tan(4) tan(6)) (3.47)



The average daily soil heat flux was approximated by Wright and Jensen (1972) as



G= (Ta -Tp) cs (3.48)



where Ta is average daily air temperature (C) at the height z, and Tp is the average daily

temperature (C) at that height for the previous three days. Parameter Cs is the general

heat conductance for the soil surface (Allen et al., 1989). G can be neglected if it makes a

very small contribution to the PET.

The slope of the saturated vapor function(A) can be calculated by taking the

derivative of the saturated vapor pressure (ed) equation with respect to temperature T

(Tetens, 1930) :











(16.78 T +117 (3.49)
ed = exp T+237.3 ( )


A = 098 ed (3.50)
(T + 237.3)2



The psychometric constant(kPa oC-') can be calculated as follows:



C P
y = (3.51)




where Cp is the specific heat of moist air at constant pressure (1.Olxl0-3 MJ kg-' C-1), P is

atmospheric pressure (kPa), ; is the ratio of molecular weights of air to water (0.622), and

k is the latent heat of vaporization (MJ kg-1), which can be calculated using the Harrison

(1963) relationship:



S= 2.5 ( 2.361x10-3 ) Ta (3.52)



where Ta is the air temperature in oC. The atmospheric pressure at a given altitude,

assuming a constant temperature lapse rate, can be calculated based on Burman et al.

(1987) as










Po To Y(z zo)
Y To



where Po and To are known atmospheric pressure (kPa) and absolute temperature (OK) at

elevation zo (m), and P is the desired pressure estimate at elevation z. The parameter Q is

the assumed constant adiabatic lapse rate. R is the specific gas constant for dry air (

286.9 Jkg-1 oK-1. Allen et al. (1989) suggested a Y value of 0.0065 K m-1 for saturated air

and of 0.01 K m-1 for dry air. The gravitational acceleration g is equal to 9.806 m s-2

Reference values for Po, To, and zo are set to those for the standard atmospheric pressure

definition at sea level, which are 101.3 kPa, 288 OK, and 0 m respectively.

Using all the above relationships, i.e., equations from (3.40) through (3.53), the

PET in equation(3. 39) can be calculated.

Determination of transpiration (or root water uptake)

Water is extracted from the unsaturated zone through plant roots. There are two

approaches to calculating root water uptake. The first one considers properties of a single

root (microscopic approach), and the second one considers the integrated properties of the

entire root system macroscopicc approach). In this study, the macroscopic approach is

followed.

The first step in root water uptake calculations is to determine the potential

transpiration rate. Then the next step is to find out how much of this potential

transpiration rate can actually occur under the restriction of soil and available moisture-

content conditions.









Potential transpiration (Tp) can be calculated as a fraction of the potential

evapotranspiration (PET) as a function of leaf area index (LAI) of the soil surface

(McCarthy and Skaggs, 1992; Fares, 1996; McKenna and Nutter, 1984). Developed by

Ritchie (1972) and modified by McKenna and Nutter (1984) the potential evaporation

from soil surface, Ep, and potential transpiration, Tp, can be calculated as follows:



Tp =PET Ep (3.67)

Ep = (PET) exp(-0.4 LAI) (3.68)



where LAI is the leaf area index. This term is defined as the ratio of total area of leaves

to the area of ground surface, and it can vary through the year depending on the type of

vegetation.

To calculate the actual amount of water taken up by a root system, a root water

uptake sink term, Wr, which represents the volume of water taken up by the roots per unit

volume of the soil in unit time [L3 L-3 T-'], was defined by Feddes et al.(1978) as



Wr(h)= ar(h)Wp(z) (3.69)



In equation (3.69) Wp (z) is the potential water uptake sink term [L3 L-3 T-'],

which is a function of depth and root density. It can be defined as the maximum possible

water uptake in favorable conditions, i.e., sufficient moisture-content around the root

zone. Wp stands for the distribution of the potential transpiration throughout the entire

root zone. The water stress response function ar (h) in equation (3.69) determines the









degree of restriction of the potential transpiration based on the available soil moisture-

content and potential transpiration rate.

Wp can be calculated using a root density distribution function and potential

transpiration rate. In the literature, many root distribution and water uptake functions can

be found (Molz, 1981). In this study, three different distribution functions were

considered.

The first root distribution model (Feddes et al., 1978) assumes a uniform

distribution of water uptake throughout the root zone, or




Wp = TP (3.70)
Zr



where Zr is the bottom of the root zone depth, and Tp is the potential transpiration [LT-1].

The second root distribution model (Prasad, 1988) is a linearly decreasing water

uptake model starting with a maximum value at the top and zero at the bottom of the root

zone, or




Wp(z) 1 (3.71)
Z r Zr



where z is the current depth and Zr is the bottom of the root zone depth.

The third root distribution model was developed for this study by modifying the

logarithmic root distribution model of Jensen (1983). His original relationship was











log W, (z) = logRo Rd z (3.72)




where Ro is a value of Wp at soil surface, and Rd is a parameter dependent on the crop and

soil. Equation (3.72) can be written as




W(z)= R0 = Cdz (3.73)
10Rdz



where Cd is a crop and soil coefficient, which has a value in the range of 0.1 < Cd < 1.

This third method is very flexible, and it can be applied to various vegetation

cover scenarios by changing the Cd values. If there is a uniform root distribution, then a

Cd value close to 1.0 is chosen. If there is a linearly decreasing crop distribution, then a

Cd value between 0.5 and 0.8 is chosen. Finally, if there is denser root distribution at the

top and much less root distribution at the bottom, i.e., there is grass and some bushes and

several pine trees in a unit area of the soil, then a Cd value less than 0.5 is chosen.

Hansen et al. (1976) gave a few characteristic values for Rd to calculate Cd values. In this

study, Cd values are assumed to be constant in time although the root depth can vary in

time. Integrating Wp(z) along the root zone gives the potential transpiration:



ZroCddz- R (CdZr -1)= T (3.74)

o In Cd









Solving equation (3.74) for Ro and substituting that into equation (3.73) gives the

following expression for the potential root water uptake function:




Wp (z) = d TP Cdz (3.75)




If Zr is a function of time, i.e., the annual vegetation with varying root depth, then

it can be calculated using Borg and Grimes (1986) relationship:



Zr (t) = ZT(0.5 + 0.5 sin[3.03(t / tT) -1.46]) (3.76)



where tT, and ZT are time to plant maturity and maximum rooting depth to be achieved at

t = tT respectively.

The water stress response function ar(h) is a prescribed dimensionless function of

the soil water pressure head (0< ar < 1). A schematic plot of the stress response function

used by Feddes et al. (1978) is shown in Figure3.6.

In this function, the wilting point is defined as the minimum moisture-content (or

corresponding pressure head) at which plant roots cannot extract any more water from the

surrounding soil.











ar(h)

1.0
= 5 mm/d
T | \p= 1 mm/d
-I



I I I (negative)
0.0 I I
h, h2 hg h3 h4 h
Absolute matrix head h

Figure 3.6 Schematic of the plant water stress response function, ar(h) (Feddes et al.,
1978). Water uptake below hi(air entrainment pressure, saturation starts) and above h4
(wilting point) is set to zero. Between h2 and h3 water uptake is maximum. The value of
h3 varies with the potential transpiration rate Tp.



In this study, a water response function ar(h) was developed by modifying the

relationship suggested by Kristensen and Jensen (1975):




C3

ar, (h) = Ta = 1- h h (3.77)
T h -hWP




where Ta, and Tp are the actual and potential transpiration respectively, and hfc and hwp are

pressure heads at field capacity and at wilting point of the soil around the root zone,

respectively. Field capacity is defined as the moisture-content (or corresponding pressure

head) at which gravitational drainage ceases. h is the pressure head around the root zone,

and C3 is a parameter greater then T p. The parameter C3 defines the shape of the water

stress function. For example, if C3 is equal to Tp, then the water stress function will

change linearly between 0 at hwp and 1 at hfc. The relationship in equation (3.77) is shown

in figure 3.7.








If equation (3.75) and equation (3.77) are substituted in equation (3.69) and

written in terms of pressure head, then the actual root water uptake model is obtained as




W,(h,z,t)= I- hf h TZ Cdz (3.78)




The integration of the root water uptake function, equation (3.78), along the root

zone gives the actual transpiration rate. The actual transpiration is restricted according to

the available moisture in the root zone, soil type, root distribution type, and the potential

evapotranspiration rate.

The model has an option to calculate the actual transpiration using a method based

on VS2D (Lappala et al., 1987). This method is relatively easy to use and requires the

input of the root pressure, which is the pressure applied by the roots on the surrounding

soil to extract water, and the root activity function. The root water uptake (wr) [T-1] for

each cell having a root zone is calculated using equation 3.79:



wr = (Ksat)hrzKr(h) r(z, t) (hroot h) (3.79)



where r (z, t) is the root activity function of depth and time [L-2], hroot is the pressure head

in the root zone [L] for the entire root system, and (Ksat)hrz is the average horizontal

conductivity (equal to (Kx + Ky)/2.0).











ar(h)


1.0 ---.--
I

0.5 increasing
STh


I I >
0.0 0.5 1.0 hfc- h'P



Figure 3.7 Water stress function as a function of pressure head and potential transpiration
(Jensen, 1983).


The total extraction by roots in a given column of cells can be calculated from



Qr (Wr dx dy dz)j (3.80)



where J is the total number cells in the column where roots are present.

If water is freely available to the plants, it is possible that a flux from the soil that

is larger than the potential transpiration rate may be computed using equations 3.79 and

3.80. Consequently, if the calculated value of Qr is larger then Tp, the value Of Wr is

T T
adjusted by the ratio T such that (wr) ij k= (Wr) ij k.
Qr Qr

The root activity function r (z, t) is defined as the length of roots in a given

volume of soil divided by that volume. This function is assumed to vary linearly between









the root activity at the bottom of the root zone and the root activity at top of the root zone.

Therefore, the root activities at the bottom and top of the root zone and the root depth as a

function of time need to be provided to the model.

Evaporation

Evaporation can take place from open water bodies or from soil surface, and it can

vary according to whether there is vegetation cover or not. Potential evaporation, Ep, is a

fraction of potential evapotranspiration (PET) and can be calculated using equation 3.68.

If there is an open water body such as lakes, rivers, etc., then Ep will be equal to PET by

equating LAI to zero in equation 3.68. Actual evaporation is determined by considering

the amount of rainfall intercepted by vegetation and the available moisture in soil.

Some researchers have reported that interception and transpiration are related and

transpiration will not start before the intercepted water is dried out by direct evaporation

(Hansen et al., 1976). To the contrary, Jensen (1979) and Van der Ploeg and Benecke

(1981) claim that both transpiration and evaporation from intercepted water can occur

simultaneously.

Interception (I) is the amount of water held by vegetation leaves during rainfall.

Jensen (1979) proposed that the maximum interception storage of the crop, Im, is linearly

related to the leaf area index, LAI:



Im= Cint LAI when P> Im (3.81a)

Im= P when P< Im (3.81b)









where Im is the maximum interception capacity [L], P is precipitation [L], and Cint is an

interception coefficient that can be taken as 0.2 for forest and for regions that have high

trees (Rutter et al., 1975) and 0.05 for short vegetation and agricultural crops (Jensen,

1979) when precipitation is measured in mm.

The demand for potential evaporation is met by intercepted water if it is sufficient.

If the intercepted water does not satisfy Ep, available water in the soil is used to satisfy the

Ep demand. The Ep demand is satisfied as long as the soil medium can conduct water to

the soil surface at a rate equal to Ep rate. As the soil near the surface becomes drier, then

soil evaporation is reduced to below the potential value. A physically based relationship

of Lapalla et al. (1987) is used for the prediction of actual evaporation from soil. The

relationship is described as



Ea = (Ksat)zKr SRES(hat htop) for Ea < Ep (3.82 a)

and

Ea= E for Ea>Ep (3.82b)



where Ea is the actual soil evaporation [LT-1], Ep is the potential soil evaporation [LT-1],

harm is the pressure potential of the atmosphere [L], htop is the pressure head at first node

on the land surface [L], and SRES is the surface resistance [L-'].

The atmospheric pressure potential hatrm can be calculated using the Kelvin

equation (Lappala et al., 1987):










hatm =RT ln(h) (3.83)
MWg



where R is the ideal gas constant [ML2T-2 K-Mol-1], T is absolute temperature [oK], Mw

is the mass of water per molecular weight [M Mol-1], hr is relative humidity [Lo], and g is

the gravitational acceleration [LT-2 ].

SRES can also be calculated, assuming that atmospheric pressure is applied to the

land surface, and the pressure at the top cell is applied to the node of the first cell.

Therefore, SRES will be the part of the conductance term between the top cell node and

the land surface such that



SRES = [2.0/DZtop] K, /(Ksatz)top (3.84)



where DZtop is the thickness of the first top cell, Kc is the saturated hydraulic conductivity

of the crust material, and Ksatz is the saturated hydraulic conductivity of the first cell on

the land surface.

The calculated actual evaporation is introduced as a negative flux boundary on the

land surface as a boundary condition.



Pumping and Recharge Wells

Wells are used for withdrawing water from the saturated zone of the aquifer or

adding water to the aquifer. The well discharge rate (Q) [L3T-1] should be specified for

each cell of the saturated zone. Negative values of Q are used to indicate a pumping well,









while positive values of Q indicate a recharge well. The specified Q value for each cell is

converted to a sink/source term Ww (Ww is part of main sink/source term Qext as

discussed above) by dividing the total Q by the number of cells. The node containing the

well itself is considered to be outside of the model, and six surrounding nodal blocks are

treated with the appropriate side as a flux boundary (Freeze, 1971). Such an approach

does not provide an exact duplication of flow conditions near the well, but it prevents the

well from becoming unsaturated immediately and unrealistically. For example, if a well

is open to five unconfined aquifer cells with uniform grid size (dx dy dz), then Ww for

each individual cell surrounding the cell containing the well itself will be calculated as

Ww = (Q/(5 dx dy.dz))/6 [T-']. If the top cell becomes unsaturated during water-table

drawdown, Ww is adjusted for the remaining four cells.

Drains, Sinkholes, and Springs

Drains are treated as specified head sink terms in this model. Drain heads are

specified for each cell in the saturated zone. Drainage flow occurs only in the saturated

zone if the water-table is above the position of the drains. The drain head should be

specified for each cell containing a drain. When the hydraulic head in the cell drops

below the specified drain head, the drainage flow ceases. It will start again if the water-

table rises above the specified drain head. Drainage flow is calculated based on the head

difference between the calculated hydraulic heads of the cell and the specified drain head.

The head difference is multiplied with the drain conductance that represents the resistance

to the flow because of material around the drain, the number of the holes in the drainpipe,

and converging streamlines in the immediate vicinity of the drain. After drain discharge

(Qd = Conductance Head Difference) is calculated, it is converted into a drain sink term









Wd (which is part of total sink term W) by dividing the calculated discharge by the total

volume of the drain cell such that Wd = Qd/(dx dy dz). The drain conductance is a

lumped parameter describing all of the head losses between the drain and the region of

the cell.

Direct recharge to an aquifer through a sinkhole can be treated as a negative drain

(source term). Springs are treated the same as drains (sink terms) by assigning a proper

conductance and a specified discharge head to each spring.




Boundary Conditions


Boundary conditions on all six sides of the flow domain must be known prior to

solving the governing groundwater flow equation. Typically, three types of boundary

conditions can be described along the boundaries: specified flux (Neumann); specified

pressure head (Dirichlet); and variable (between Neumann and Dirichlet) boundary

conditions. In each boundary condition, the prescribed values either can be constant or

they can vary with time.

Specified Flux Boundary Condition

Specified flux boundary conditions can be used to describe the rainfall

(infiltration), evaporation, and seepage processes. These processes are treated as sink or

source terms at the boundary element faces. No-flow boundaries and impervious

boundaries with zero flux also can be classified in this category. A specified boundary

condition can be defined formally as









q1/2(h)= qp(Xb,Yb, Zb,t) (3.85)



where q1/2 (h) is the flux at the boundary, and qp is the specified flux (evaporation or

infiltration rates) at the boundary nodes.

A detailed mathematical description of the boundary conditions in terms of the

finite-difference method and iteration procedures is presented in chapter 4 of this study.

Specified Head Boundary Condition

Specified head boundary conditions can be defined if there is a constant head

water body such as a lake, river, etc. The specified head boundary condition can be

written formally as



h= hd (Xb,yb,Zb,t) (3.86)



where h is the pressure head at the boundary node, and hd is the specified head at the

boundary coordinate of (xb, yb, Zb).

Variable Boundary Condition

Variable boundary conditions are used to describe the evaporation from the soil

surface and infiltration due to rainfall. These two hydrologic events have two stages such

that in one stage they are described as a flux boundary and in the other stage they are

treated as a constant specified head boundary condition. Variable boundary conditions

are called "variable" because they vary between flux boundary and specified head









boundary conditions during the simulation depending on the potential evaporation, the

conductivity of medium, and the availability of water such as rainfall.

Infiltration is a two-stage procedure. In the first stage, all rainfall enters the

system at the applied rate as long as the conductivity and sorptive capacities of the

medium are not exceeded. If these capacities are exceeded, water ponds on the surface

and infiltration decreases asymptotically to a value equal to saturated hydraulic

conductivity of the medium (Lappala et al., 1987). Rubin (1966) and Smith (1972)

presented extensive discussions of the ponding process and reported that surface runoff

cannot occur until ponding has begun. Using this concept, a maximum ponding depth

value (hpmax) is assigned for each top surface cell in the flow domain to handle the

rainfall-runoff-infiltration procedure. At land surface, two boundary conditions are

possible:



ql/2(h)= pp((Xbb,Yb,bt) if t
and

h = hpmax(xb,yb, Zb,t) if t> tp (3.87b)



where q1/2 (h) is the vertical flux at the boundary, pp is the rainfall rate flux, and tp is time

of ponding, which is determined during the iterations.

Evaporation is also a two-stage procedure at the boundary nodes, and it is

analogous to precipitation. It is dependent on both the potential evaporative demand of

the atmosphere and the ability of the porous medium to conduct water to the surface.

During the first stage of evaporation, there is an outward flux boundary at the surface.









This continues as long as water is conducted to the top layer soil where the soil moisture-

content is greater than the residual moisture-content, which can also be described as (hmin)

in terms of pressure head. In the second stage, evaporation ceases because of a lack of

moisture, and boundary conditions are set to hmin minimum pressure so that no more

water can be extracted from the soil. These two stages are described mathematically as



qj/2(h)= Ea(XbYb,Zb,h,t) if hl/2 > hmin (3.88a)

and

h= hmin(Xb, Yb, Zb,t) if hi/2 < hmin (3.88b)



where Ea is the actual evaporation, hl/2 is the pressure head at the face of the boundary

cell, and hmin is the minimum pressure at the residual water content level.

River Boundary

In the model, one-dimensional steady-state open-channel flow can interact with

the underlying porous medium. It is necessary to specify the river heads for each time

step. Actually the river boundary acts as a specified sink/source term in the top boundary.

The flow exchange between the river and the porous media can be calculated using

Darcy's equation based on the head gradient between the river and the underlying porous

media and the hydraulic conductance of the river bed material:




qr = Cr Hr (3.87)
DZr









where qr is the aquifer river exchange per unit length of the river [L3T-1L-1], Cr is the

conductance term for the river bed (calculated by averaging the river bed hydraulic

conductivity and the hydraulic conductivity of the first cell beneath the river bed), Hr is

the river head, H is the hydraulic head in the first cell beneath the river bed, and DZr is

the distance between the river bottom and the first cell beneath the river bed.

At each time step, using the prescribed river heads, Hr, and heads of the porous

medium H from the previous iteration level, the flow exchange qr between the river and

underlying groundwater system is calculated using equation (3.87). In the new iteration

level, the calculated qr is used as a specified flux boundary condition for top cells having

river segments superimposed on them. This procedure is followed for each river cell for

every time step.



General Head Boundary

A general head boundary can be used to provide a flow into or out of an active

cell at a boundary from an external source far from the boundary. That flow is calculated

as a function of the head difference between the active cell and the external source and

the conductance between the external source and the boundary cell. Functionally, this

type of boundary is similar to a drain or a river boundary (McDonald and Harbaugh,

1988). To specify a general head boundary condition, the head at the external source

(Hext) and the distance (Xghb) between the source and the boundary cell (which is used to

calculate the conductance between the boundary and the external source) need to be

specified. The flow between the external source and the boundary cell is calculated using

equation 3.88.






87





qghb K (Hext H) (3.88)
Xghb



where qghb is flux into or out of the boundary cell [LT-1], K is the hydraulic conductivity

between the boundary cell and the external source [LT-1], Xgbh is the distance between the

boundary cell and the external source [L], Hext is the specified head of the external source

[L], and H is the head at the boundary cell.















CHAPTER 4
MATHEMATICAL MODEL DEVELOPMENT AND NUMERICAL SOLUTION OF
THE MODIFIED RICHARDS EQUATION

As described in this chapter, a new variably saturated rainfall-driven groundwater-

pumping model has been developed. The model simulates a hydrogeologic system by

solving the nonlinear, three-dimensional form of the modified Richards equation

continuously throughout the whole flow domain from ground surface to the impervious

bottom of the lowest layer. The finite-difference method with a variable finite-difference

grid is used to solve the governing equation. The upper boundary in the model is at land

surface, and the upper boundary conditions are determined using soil and meteorological

data. The model treats the complete subsurface regime as a unified whole, and the flow

in the unsaturated zone is integrated with saturated flow in the underlying unconfined and

confined aquifers. The model allows modeling of heterogeneous and anisotropic

geologic formations. A plant root water uptake (transpiration) model and an evaporation

model are included in the governing flow equation as a sink term and boundary condition,

respectively, in the model.

In the previous chapter, the partial differential form of the three-dimensional

modified Richards equation was developed. Boundary conditions and sink/source terms

were mathematically described. This new numerical model was developed based on

those partial differential equations. The model was mathematically conceptualized and

developed using finite-difference approximation methods. The resulting finite-difference