Citation |

- Permanent Link:
- https://ufdc.ufl.edu/UF00100699/00001
## 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:
- 1999
- 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 )
## UFDC Membership |

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 |