M
Z

EvaluaƟng InteracƟons between Groundwaterand Vadose Zone Using the HYDRUS-Based
Flow Package for MODFLOW

S

S

:V

Navin Kumar C. Twarakavi,* Jirka Šimůnek, and Sophia Seo
In the past, vadose zone processes in groundwater flow models were drama cally simplified (or even neglected) due toconstraints on computa onal resources. The one-dimensional unsaturated flow package HYDRUS, recently developed
for the groundwater model MODFLOW, was evaluated and compared with other contemporary modeling approaches
used to characterize vadose zone eﬀects in groundwater models. Being fully incorporated into the MODFLOW program,the HYDRUS package provides MODFLOW with recharge fluxes at the water table, while MODFLOW provides HYDRUS withthe posi on of the groundwater table that is used as the bo om boundary condi on in the package. The performance ofthe HYDRUS package was analyzed for three case studies of increasing complexity: (i) a one-dimensional infiltra on experiment; (ii) a two-dimensional water table recharge experiment; and (iii) a hypothe cal regional-scale groundwater flow
problem. The computa onal need and modeling eﬃciency of the HYDRUS package was compared with other relevant
MODFLOW packages (VSF, UZF1, and REC-ET). For smaller scale problems (up to two dimensions), the VSF process andthe HYDRUS and UZF1 packages performed comparably well in terms of modeling eﬃciency and simula on mes.
Because of the high computa onal demand, it was not feasible to use the VSF process on a typical personal computer
for the hypothe cal large-scale groundwater problem. The HYDRUS package provided a much more eﬃcient alternave to VSF for this large-scale groundwater problem and could be er account for vadose zone processes than the UZF1
and REC-ET packages. For large-scale groundwater problems, the HYDRUS package provides an op mal tradeoﬀ between
computa onal eﬀort and accuracy of model simula ons for coupled vadose zone–groundwater problems.
A: ET, evapotranspiration; REC-ET, Recharge–Evapotranspiration package; UZF1, Unsaturated Zone Flow package; VSF,
Variably Saturated Flow package.

W

through the variably saturated (vadose)zone is an important part of the hydrologic cycle
because it inﬂuences partitioning of water among various ﬂow
components. Depending on hydrologic, geologic, and soil characteristics, rain and snowmelt are partitioned at the land surface
into runoﬀ, inﬁltration, evapotranspiration (ET), groundwater
recharge, and vadose zone storage (Fig. 1). Water ﬂow in thevadose zone especially aﬀects the transfer rates between the land
surface and the groundwater table, which are two key hydrologic boundaries. Evaluation of almost any hydrologic process,therefore, requires that water ﬂow through the vadose zone isappropriately taken into account. Modeling of vadose zone ﬂow
processes, however, is a complex and computationally demanding
task that is often handicapped by the lack of data necessary tocharacterize the hydraulic properties of the subsurface environment. Consequently, vadose zone ﬂow processes have rarely been

properly represented in hydrologic models (Ward, 2002; Scanlon,
2002; Keese et al., 2005). For example, models that simulate surface and near-surface hydrology usually oversimplify the impact ofvadose zone ﬂow processes and rarely consider three-dimensional
regional groundwater ﬂow. Similarly, regional-scale groundwatermodels often simplify vadose zone ﬂow processes by calculating
groundwater recharge externally without proper considerationof changes in groundwater levels (e.g., Lorenz and Delin, 2007;
Shah et al., 2007; Uddameri and Kuchanur, 2007). To overcome
this frequent simpliﬁcation, there is an urgent need for methods
that can eﬀectively simulate water ﬂow through the vadose zone
in large-scale hydrologic models (Winter et al., 1998). This issue
is especially important for groundwater models.

N.K.C. Twarakavi and J. Šimůnek, Dep. of Environmental Sciences, Univ.
of California, Riverside, CA 92521; S. Seo, Colorado School of Mines, 1500
Illinois St., Golden, CO 80401. Received 25 Apr. 2007. *Corresponding
author (navink@ucr.edu).
Vadose Zone J. 7:757–768
doi:10.2136/vzj2007.0082© Soil Science Society of America
677 S. Segoe Rd. Madison, WI 53711 USA.
All rights reserved. No part of this periodical may be reproduced ortransmi ed in any form or by any means, electronic or mechanical,
including photocopying, recording, or any informa on storage andretrieval system, without permission in wri ng from the publisher.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

F . 1. A schema c showing the processes (including the key
vadose zone processes) aﬀec ng subsurface hydrology.
757

Most traditional attempts at characterizing vadose zone
ﬂow processes in groundwater models follow the water budget
or “residual” approach (Scanlon, 2002). Water budget methods
are based on the water budget equation. The water budget equation is a mathematical representation of the fact that all water
arriving at the water table leaves the system as groundwater ﬂow,
is discharged through sinks such as surface water bodies, is evapotranspirated, or is retained as storage (Scanlon, 2002). One may
indirectly estimate the water table recharge from the water budget
equation by measuring or estimating all other components in thewater budget. An example of the water budget approach is the useof the Recharge and Evapotranspiration (REC-ET) packages in
MODFLOW (Harbaugh et al., 2000), a modular three-dimensional ﬁnite-diﬀerence groundwater ﬂow model. MODFLOW
was developed by the U.S. Geological Survey and is one of themost widely used groundwater ﬂow models. It is obvious that
the REC-ET package tends to oversimplify and underestimate
the eﬀects of vadose zone ﬂow on groundwater. In spite of its
widespread use, the REC-ET approach suﬀers from the following
major limitations: (i) the method is not reliable for deeper watertables; and (ii) the applicability of this approach is questionable
for arid and semiarid regions where soil capillary pressures play adominant role in vadose zone ﬂow (Gee and Hillel, 1988; Lerner
et al., 1990; Hendrickx and Walker, 1997).
A more promising approach to properly represent vadosezone ﬂow processes in groundwater models involves couplinggroundwater and vadose zone models. A coupled model simulates
the eﬀects of near-surface hydrologic processes on groundwaterﬂow by linking a groundwater model with a selected vadose zone
model in space and time. The majority of currently available
vadose zone models are based on either the Richards equation(Richards, 1931) or the kinematic wave equation (Colbeck, 1972;
Smith, 1983; Smith and Hebbert, 1983). While the Richards
equation considers ﬂow due to both capillary and gravity forces,
the kinematic wave equation neglects capillarity and considers only gravity. In coupled models, the groundwater recharge
is calculated internally in the model based on existing surfacehydrologic conditions and water table levels. By simultaneously

considering surface meteorological conditions, water table levels,and the hydraulic properties of the vadose zone, coupled models
represent reality better than traditional approaches such as the
REC-ET package; however, evaluation of interactions between
the near-surface and groundwater ﬂow processes using coupledmodels has been a desirable but diﬃcult goal.
It is desirable to develop MODFLOW packages other than
REC-ET that would better account for processes in the vadosezone. By combining these packages with MODFLOW, not only
can the vast groundwater modeling capabilities of MODFLOW
be harnessed, but these new numerical packages can be quickly
distributed among the large number of MODFLOW users. The
virtues of a coupled vadose zone–groundwater model should be
evaluated based on the following criteria: (i) accuracy in representation of the physical processes that drive vadose zone ﬂow,(ii) usability for diﬀerent groundwater modeling scenarios, (iii)applicability to diﬀerent spatial and temporal scales, i.e., from lab
or ﬁeld to regional spatial scales and from hourly to decadal temporal scales, and (iv) applicability to diﬀerent meteorological andclimactic conditions, such as humid, arid, and semiarid regions.
Three MODFLOW packages accounting for processes in thevadose zone have been recently developed: the Variably Saturated
Flow (VSF) process (Thoms et al., 2006), the Unsaturated
Zone Flow (UZF1) package (Niswonger et al., 2006), and the
HYDRUS package (Seo et al., 2007). Table 1 lists the strengths
and weaknesses of selected currently available approaches thatincorporate vadose zone ﬂow into MODFLOW. It may be notedthat all currently available coupled modeling techniques havesome weaknesses and some strengths.
Among the available packages for MODFLOW, the VSF
process most robustly represents the vadose zone processes, as it
can consider all major variably saturated ﬂow processes as wellas their three dimensionality. The thorough consideration of thevadose zone ﬂow processes in the VSF process, however, makes it
computationally very demanding. On the other hand, the UZF1
and REC-ET packages radically simplify vadose zone processes.
As a result, these packages are computationally eﬃcient but may

T
1. A comparison of the Recharge–Evapotranspira on (REC-ET), Unsaturated Zone Flow (UZF1), and Variably Saturated Flow (VSF)packages for MODFLOW that incorporate the eﬀects of water flow in the vadose zone.
Characteris c
Descrip on of vadose zone
processes
Characteriza on of vadosezone flow processes
Vadose zone flow modelingcapabili es
Applicability across larger
spa al domains
Computa onal eﬀort
Range of groundwatermodeling scenarios
Applicability to arid andsemiarid areas

REC-ET
water balance at the land
surfacepoor

UZF1
one-dimensional kinema cwave equa on
be er

VSF
three-dimensional Richards
equa onbest

REC and ET only

ET, root water uptake, andinfiltra on

yes

yes

low
large†

moderatelarge†

ET, root water uptake,surface ponding, andinfiltra on
diﬃcult due to computa onal yes
eﬀort
very highmoderatelarge†
large†

not very applicable; needs
extensive calibra onvery applicable for deep
highly applicablewater tables, lessapplicable for shallow
water tables where
capillary rise may becomea significant process† Since it is based on MODFLOW.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008758

HYDRUS
one-dimensional Richards
equa onbe er than UZF1, not as goodas VSF
many capabili es as available
in HYDRUS-1D

highly applicablenot necessarily provide an accurate characterization of vadosezone ﬂow processes.
The objective of this study was to brieﬂy review the aforementioned approaches (i.e., the VSF process and the UZF1
and REC-ET packages) used to account for vadose zone ﬂow
in MODFLOW and then to compare them, using case studies of increasing complexity, with the HYDRUS package. The
HYDRUS package for MODFLOW was developed to provide
a balance between computational eﬃciency and accuracy. Being
one dimensional, the HYDRUS package signiﬁcantly simpliﬁes
the calculations but cannot consider water ﬂow in the vadosezone in multiple dimensions.
It is important to note that when calibrated against collectedﬁeld data, one may expect the REC-ET and UZF1 packages
to perform relatively well for many practical applications. This
study looked only at how accurately these approaches (i.e., VSF,HYDRUS, UZF1, and REC-ET) perform relative to each otherwhen ﬁeld-estimated soil hydraulic parameters or their literature
values are given. Calibration, which is not considered here, may
lead to estimates of soil hydraulic parameters that do not always
correspond to ﬁeld-estimated values.

more computationally demanding treatment of water ﬂow inboth the unsaturated and saturated zones. The large computational demand stems from the fact that the numerical solutionof the Richards equation requires much ﬁner discretization of
three spatial dimensions and smaller time steps than traditional
groundwater models. This seriously limits the applicability of
the VSF process for regional-scale groundwater ﬂow problems(domains >100 km2) (Thoms et al., 2006).
The Unsaturated Zone Flow Package

A number of researchers (e.g., Pikul et al., 1974; Refsgaard
and Storm, 1995; Niswonger et al., 2006) have proposed a simpler methodology that signiﬁcantly decreases the computational
demand without greatly compromising the eﬃciency of thecoupled modeling approach. The proposed approach involves
coupling a one-dimensional vadose zone ﬂow model with a threedimensional groundwater ﬂow model (such as MODFLOW).
Pikul et al. (1974) and Niswonger et al. (2006) noted that thisapproach probably provides the most efficient solution forgroundwater ﬂow models, especially for large-scale applications.
This approach, i.e., consideration of only one-dimensional verticalﬂow in the unsaturated zone and fully three-dimensional groundwater ﬂow, has been used, for example, in the MIKE SHE model(Refsgaard and Storm, 1995) and the UZF1 package (Niswonger
et al., 2006) for MODFLOW.
The UZF1 package couples a vadose zone ﬂow model basedon the numerical solution of the one-dimensional kinematicwave equation with MODFLOW. Unlike the Richards equation,
which considers both gravity and capillarity as driving forces for
ﬂow in the vadose zone, the kinematic wave equation considers
only gravity-driven ﬂow. The model, based on the kinematic wave
equation, relates water ﬂuxes directly to the degree of saturation.
The Brooks and Corey model (Brooks and Corey, 1964), one of
the commonly used models relating the moisture content, θ, tothe unsaturated hydraulic conductivity, K(θ), (and the ﬂux, q), isused in the UZF1 package:

MODFLOW Packages AccounƟng for Vadose
Zone Processes
The Recharge–Evapotranspira on Package

The Recharge (REC) and Evapotranspiration (ET) packages
can be used together to provide a simplistic characterization ofvadose zone processes in MODFLOW. The REC package is used
to simulate a speciﬁed downward recharge ﬂux across the top of
the model domain. The recharge ﬂux can be varied spatially andwith time. To estimate the volumetric ﬂow rates at the water table,
these ﬂuxes are simply multiplied by the horizontal area of thecells. The REC package is a simpliﬁed representation of vadosezone ﬂow and does not consider vadose zone ﬂow processes suchas storage and runoﬀ. The ET package in MODFLOW is used tosimulate the discharge of water to evaporation and transpiration.
In the ET package, a maximum evapotranspiration rate is supplied to the model as a function of space and time. To consider
the inﬂuence of the groundwater depth on evapotranspiration
rates, the ET package uses a user-deﬁned extinction depth. While
the ET package may be simple to use and understand, it may
tend to oversimplify the impact of vadose zone processes. Also,
the necessity to supply several rather empirical parameters (e.g.,
the extinction depth and maximum evapotranspiration rates)
increases the uncertainty of modeling results.

⎛ θ − θ r ⎞⎟ε⎟q = K (θ ) = K s ⎜⎜[1]⎜⎝ θ s − θ r ⎠⎟⎟where q is the water flux [L T−1], K(θ) is the unsaturated
hydraulic conductivity [L T−1] expressed as a function of thewater content θ (dimensionless), Ks is the saturated hydraulicconductivity (L T−1), ε is the Brooks–Corey exponent (dimensionless), and θs and θr are the saturated (porosity) and residual
water contents (dimensionless), respectively. The UZF1 package
considers evaporation and root water uptake (transpiration) byassuming that the water loss occurs instantaneously in the soil
proﬁle between the soil surface and a user-speciﬁed depth called
the ET extinction depth (see Fig. 1).
Application of the kinematic wave equation for vadosezone modeling has its own advantages and disadvantages. Many
researchers have debated whether or not variably saturated ﬂow
can be treated using the kinematic wave approach (e.g., Singh,2002). The applicability of the kinematic wave equation (such as
in the UZF1 package) to simulate vadose zone ﬂow depends inlarge part on the soil hydraulic properties, climatic conditions, andthe depth to the groundwater table. The UZF1 package requiresthat the unsaturated zone is homogeneous, which can signiﬁcantly limit applicability of the package. While in coarse-textured

The Variably Saturated Flow Process

Incorporating the numerical solution of the three-dimensional Richards equation into the groundwater ﬂow model is themost accurate way to represent the complex nature of physicalprocesses in the unsaturated part of the subsurface. An example
of such an approach is the VSF process (Thoms et al., 2006) for
MODFLOW. The VSF process solves the three-dimensional formof the Richards equation for the entire MODFLOW domain.
In the VSF process, the ﬁnite-diﬀerence MODFLOW domain
is expanded to include the variably saturated zone, and the“mixed form” of the Richards equation is used as the governing
equation. The VSF process thus oﬀers more rigorous but much
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008759

soils, deep vadose zone proﬁles, or humid climates gravity usually
dominates ﬂow in the unsaturated zone and thus the kinematicwave approach is applicable, in ﬁne-textured soils, proﬁles withshallow groundwater levels, or arid climates, neglecting capillary
forces may lead to signiﬁcant errors. Under such conditions, thekinematic wave equation may fail to describe the dominant ﬂow
processes. The UZF1 package can thus be applied mainly in situations where gravity-dominated water ﬂow occurs.

The ﬁnal HYDRUS package thus simulates only one-dimensional water movement in variably saturated porous media. The
Galerkin-type linear ﬁnite-element scheme is used in HYDRUS
to numerically solve the Richards equation.
Water uptake by plant roots has a great eﬀect on water inthe root zone. Root water uptake is represented in HYDRUS as
an extraction or sink term, S(h), that distributes the potential
transpiration across the root zone. Feddes et al. (1978) described
the sink term as

The HYDRUS Package

S ( h ) = α ( h )S p

The HYDRUS package (Seo et al., 2007) was developed for
the MODFLOW-2000 (Harbaugh et al., 2000) environment
to combine extensive modeling capabilities of both HYDRUS
and MODFLOW. The HYDRUS package incorporates into
the MODFLOW suite a vadose zone ﬂow model based on the
one-dimensional Richards equation. The package was developed
to consider the eﬀects of precipitation, inﬁltration, evaporation,
plant water uptake, soil moisture storage, and water accumulation at the ground surface and in the vadose zone. It is basedon the HYDRUS-1D program (Šimůnek et al., 2005, 2008),
which simulates one-dimensional water movement in the variably saturated zone.
In the coupled HYDRUS–MODFLOW system, vadose zone
and groundwater ﬂows are modeled using two separate governing
equations. Similarly to the UZF1 package, groundwater ﬂow in
MODFLOW is modeled by solving the following mass-conservation equation using a ﬁnite-diﬀerence approximation:[4]

suggesting that the root extraction, S, depends on the pressure
head, h, the potential root water uptake rate, Sp [T−1], and thestress response function α(h), which characterizes plant response
to water stresses.
The HYDRUS package for MODFLOW, similarly to the
UZF1 and REC-ET packages, does not take into account subsurface runoﬀ because of the one-dimensional nature of the package.
The impact of subsurface runoﬀ needs to be considered independently when these packages for MODFLOW are used. On theother hand, the HYDRUS package can consider surface runoﬀ.
The HYDRUS package for MODFLOW has an option wherein any
excess water on the soil surface can either accumulate there or be
immediately removed by surface runoﬀ (Šimůnek et al., 2005).
Spa al Discre za on∂h∂twhere Kx, Ky, and Kz are hydraulic conductivities [L T−1] in the
direction of the x, y, and z coordinates, respectively; h is the pressure head [L], W is the volumetric ﬂux per unit volume throughsources or sinks [T−1], Ss is the speciﬁc storage of the porous
material [L−1], and t is time [T].
In the HYDRUS package, vadose zone water ﬂow is describedmathematically using the modiﬁed Richards equation:

The eﬃciency of a coupled vadose zone–groundwater model
depends to a large extent on how these two submodels interact with each other in space and time. The MODFLOW model
uses the ﬁnite-diﬀerence approximation of the mass conservation
equation to simulate groundwater ﬂow. The groundwater modeling domain for MODFLOW is discretized into grids or blocks as
described in Harbaugh et al. (2000) and the number of vadosezone proﬁles may be as large as the number of rows and columnsof this ﬁnite-diﬀerence grid. Based on similarities in soil hydrology, topographical characteristics, and depth to groundwater,
the discretized MODFLOW domain can be divided into zones,
which comprise one or more cells of the MODFLOW model (Fig.
2). One HYDRUS soil proﬁle is then assigned to each of these

⎤
∂θ
∂ ⎡
∂h
=
[3]
⎢ K ( h ) − K ( h )⎥ − S⎥⎦∂ t ∂ z ⎢⎣∂zwhere θ is the volumetric water content (dimensionless), h is thesoil water pressure head [L], t is time [T], z is the vertical coordinate [L], S is the sink term usually accounting for root water
uptake [T−1], and K(h) is the unsaturated hydraulic conductivity[L T−1] as a function of h or θ. Note that the Richards equation
is highly nonlinear due to the dependence of the unsaturated
hydraulic conductivity, K(h), and the water content, θ(h), on thecapillary pressure head, h. The two most widely used approaches
representing these nonlinear relationships are the Brooks andCorey (Brooks and Corey, 1964) and van Genuchten–Mualem(van Genuchten, 1980) models. Both models are available in the
HYDRUS package.
The computer program HYDRUS-1D (Šimůnek et al., 2005)
was adapted and simpliﬁed for the HYDRUS package. The simpliﬁcation involved removal of subroutines simulating soluteand heat transport, hysteresis in the soil hydraulic functions, andboundary conditions that were irrelevant for the coupled model.

F . 2. A discre zed aquifer system in MODFLOW and two associated HYDRUS soil profiles. One HYDRUS soil profile is assigned toeach MODFLOW zone. Note that the discre za on of HYDRUS soil
profiles is much finer than that of the MODFLOW domain.

∂ ⎜⎛∂ h ⎞⎟ ∂ ⎜⎛∂ h ⎞⎟ ∂ ⎛⎜∂ h ⎞⎟⎟+Kx
Kz⎟ + ⎜⎜ K y
⎟ −W
⎟
⎟
⎜
⎜
∂x ⎝
∂x ⎠ ∂y ⎝∂ y ⎠⎟ ∂ z ⎝∂ z ⎠⎟[2]= Ss

www.vadosezonejournal.org · Vol. 7, No. 2, May 2008760

zones (Fig. 2). It is assumed that the HYDRUS
soil proﬁle adequately represents vadose zone ﬂow
for the entire zone.
Created vertical soil profiles are then discretized vertically into finite elements. The
ﬁnite-element mesh is constructed by splitting
the soil proﬁle into one-dimensional elementsthat are connected to each other at nodal points.
Once the ﬁnite elements are constructed, they
may not be changed during the simulation. Careshould therefore be taken to ensure that the depthof the soil proﬁle extends from the soil surface
to the deepest possible water table that may be
expected during the simulation. To ensure convergence of the numerical solution, ﬁnite-element
dimensions should be relatively small at locations
where sharp pressure head gradients are expected.
Such smaller elements are usually needed close to
the soil surface where meteorological factors cancause rapid changes in water content and pressure
head gradients, and at interfaces between diﬀerent soil horizons. Soil texture also needs to beconsidered during the discretization process. For
example, coarse-textured soils generally requireﬁner discretization than ﬁne-textured soils dueto the higher nonlinearity of their soil hydraulic properties. Once the spatial discretization of
the soil proﬁle is performed, the distribution of
diﬀerent soil materials in the proﬁle needs to be
described (Fig. 2) (Seo et al., 2007).
Time Discre za on

F . 3. Flowchart describing the coupled modeling approach used in the HYDRUS package for MODFLOW: (a) steps shown in gray correspond to the treatment of variably
saturated water flow (i, stress period; j, me step; k, soil profile number; nt, number ofme steps; ns; number of stress periods; np, number of HYDRUS profiles); (b) calculaons carried out by the HYDRUS package during one MODFLOW me step.

The computational eﬃciency of the coupled
HYDRUS–MODFLOW system is enhanced bysimulating vadose zone and groundwater ﬂows
at their own, often diﬀerent, time steps. Thisis needed because a proper treatment of the
Richards equation requires smaller time steps than those usuallyused in MODFLOW simulations. Figure 3 describes the coupling procedure used in the HYDRUS package. The two models(HYDRUS and MODFLOW) interact, i.e., exchange information about the groundwater recharge and the groundwater level,
only at the end of each MODFLOW time step, during which
HYDRUS may perform multiple time steps to simulate unsaturated zone ﬂow. MODFLOW receives the recharge ﬂux fromHYDRUS and calculates a new water table depth for the nexttime step. A new water table depth is calculated and assigned as
the pressure head bottom boundary condition in the HYDRUS
package for the next MODFLOW time step. The iteration procedure in the HYDRUS package is similar to that described in the
HYDRUS-1D manual (Šimůnek et al., 2005). See this reference
and Seo et al. (2007) for more details.

regional-scale groundwater ﬂow problem. While the HYDRUS,
UZF1, and VSF packages were used in the ﬁrst case study, only
the HYDRUS and UZF1 packages were applied in the secondcase study since results for the VSF process for this application
can be found in the literature. Finally, the REC-ET, UZF1, andHYDRUS packages were used in the third case study.
The one-dimensional Las Cruces inﬁltration experiment ofWierenga et al. (1991) was used ﬁrst to evaluate the eﬀectivenessof the HYDRUS, VSF, and UZF1 packages to simulate ﬂow
in the vadose zone without considering groundwater ﬂow. The
water table recharge experiment of Vauclin et al. (1979) was thenused to evaluate whether a combination of a one-dimensional
vadose zone module with a groundwater model can approximate
this obviously two-dimensional problem. Finally, a complex
regional-scale groundwater ﬂow problem was used to evaluate
the eﬀectiveness of diﬀerent vadose zone packages in accounting
for various processes in the vadose zone.

Case Studies
The performance of the HYDRUS package was analyzed
and compared with other vadose zone ﬂow packages in three casestudies: (i) the Las Cruces one-dimensional inﬁltration experiment (Wierenga et al., 1991), (ii) the two-dimensional water table
recharge experiment (Vauclin et al., 1979), and (iii) a hypotheticalwww.vadosezonejournal.org · Vol. 7, No. 2, May 2008

Case Study 1: One-Dimensional Infiltra on Experiment

The ﬁrst case study involved the one-dimensional inﬁltration
experiment at the Las Cruces trench site (Wierenga et al., 1991).
The experiment involved a comprehensive ﬁeld study, conducted761

sure heads at depth z, respectively. A soil proﬁle depth of 600 cm
was used and the simulation was run for 35 d.
Experimental results of the Las Cruces trench inﬁltration
experiment (Wierenga et al., 1991) are compared with results
simulated using the HYDRUS, VSF, and UZF1 packages in Fig.
4. Figure 4 shows the soil water content proﬁles for diﬀerent daysof the experiment and compares model predictions of the VSF,UZF1, and HYDRUS packages with the experimental data. As
expected, the HYDRUS package and the VSF process performed
similarly as they both solve the same Richards equation for
one-dimensional problems. The UZF1 package only slightly overpredicted water contents behind the wetting front. A comparisonof the time needed for the simulations by the VSF process andthe UZF1 and HYDRUS packages was done. It was noted thatthe computational demand of the VSF, UZF1, and HYDRUS
packages was similar for the one-dimensional case study.
The UZF1, HYDRUS, and VSF packages provided similar results for this one-dimensional inﬁltration experiment andneeded comparable computational times.

in southern New Mexico, the primary purpose of which was todevelop a data set for validating and testing numerical models.
For this purpose, the study site was heavily instrumented withneutron probes, tensiometers, and solute samplers for measuring water contents, pressure heads, and solute concentrations(Wierenga et al., 1991), respectively. More than 500 soil samples(undisturbed and disturbed) were taken at the experimental site
and analyzed in the laboratory for bulk density and to ﬁnd thesaturated hydraulic conductivity and the soil water retention
curve. The inﬁltration study involved application of water to a
4-m-wide area using closely spaced drips with an average surface
ﬂux of 1.82 cm d−1 for 86 d of the experiment. To reduce thedisruption of the experimental conditions by rain and evaporation, the irrigated area and its surroundings were covered by apond liner.
Wierenga et al. (1991) performed a one-dimensional simulation of the inﬁltration experiment using a numerical model basedon the ﬁnite-diﬀerence approximation of the Richards equation.
They considered a uniform soil proﬁle with an equivalent saturated hydraulic conductivity Ks of 270.1 cm d−1. The RETC code(van Genuchten et al., 1991) was used to analyze the retention
curve data for undisturbed and disturbed soils (>500 soil samples),resulting in the following retention curve parameter values (vanGenuchten, 1980): θs = 0.321, θr = 0.083, α = 0.055 cm−1, andn = 1.51. These values were then used by the HYDRUS package
and VSF process. A zero extinction depth was used in the UZF1
package, as no evaporation losses were considered. Morel-Seytoux
et al. (1996) developed equations describing the parameter equivalence between the Brooks–Corey exponent and van Genuchten
parameters. From these equations, the Brooks–Corey exponent
for the UZF1 package was estimated to be 6.92.
The one-dimensional simulation of the inﬁltration experiment was performed using MODFLOW with the HYDRUS, VSF,
and UZF1 packages. The ﬁnite diﬀerence mesh for MODFLOW
consisted of a one-cell grid. Initial pressure heads (hi = −100 cm)
in the soil proﬁle were the same as those used by Wierenga etal. (1991). While a constant water ﬂux was used as the upperboundary condition (q0 = 1.82 cm d−1), freedrainage was considered at the lower boundary. The implicit assumption in this boundary
condition is that the groundwater table is deep
enough so that it does not aﬀect ﬂow in the soil
proﬁle. The initial and boundary conditions interms of the water content, θ(z, t), are describedas follows:θ ( z ,0 ) = θ init ( z ) = θ [ hinit ( z ) ][5a]q ( z = 0, t ) = 1.82 cm d − 1[5b]θ (b, t ) = θ init (b )[5c]where b is the depth of the soil proﬁle, whichmust be large enough so that the wetting front
does not aﬀect the water content at the bottomof the soil proﬁle during the simulation, θinit(z)and hinit(z) are initial water contents and preswww.vadosezonejournal.org · Vol. 7, No. 2, May 2008

Case Study 2: Two-Dimensional Water Table Recharge Experiment

The HYDRUS and UZF1 packages were used to model thetwo-dimensional transient water table experiment of Vauclin etal. (1979). The same data set was used previously by Thoms etal. (2006) to evaluate the VSF process. See these references fortheir results.
The experimental setup consisted of a 6.0- by 2.0-m box
containing sandy soil. The initial water table elevation was 0.65 m
from the bottom. A constant ﬂux of q = 3.55 m d−1 was applied
across the center 1.0 m of the soil surface while the rest of thesurface was covered to prevent evaporation. Due to the symmetryof the experiment, only one half of the experiment was modeled and the model domain was thus 3.0 by 2.0 m. The initial
total head was set equal to 0.65 and the right boundary cells
were constrained to the initial water table position throughout
the 8-h simulation. The grid was discretized into uniform cellsof 0.1-m width and 0.05-m depth. For comparison purposes,

F . 4. A comparison of measured water contents and corresponding es mates calculated using the HYDRUS, Variably Saturated Flow (VSF), and Unsaturated Zone Flow(UZF1) packages for (a) Day 19 and (b) Day 35 for the Las Cruces trench experiment (datafrom Wierenga et al., 1991).
762

simulations were performed using the exact setup described in
the documentation of the VSF process. Only two soil proﬁlesrepresenting the soil directly below the recharge zone and therest of the transport domain were used in calculations with the
HYDRUS and UZF1 packages. The saturated hydraulic conductivity of 840 cm d−1 was used. The initial and boundary conditions
are described as follows:
Domain0 ≤ x ≤ 3 m, 0 ≤ z ≤ 2 m, 0 ≤ t ≤ 8 h
Δ x = 0.1 m, Δ z = 0.05 m, Δ t = 1 min
Initial conditionh ( x , z ,0 ) = 0.65 − z [in meters][6]

Boundary conditionq ( 0 ≤ x ≤ 0.5,2, t ) = 3.55 m d − 1
q ( x > 0.5,0, t ) = 0.00 m d − 1
q ( 0, z , t ) = 0.00 m d − 1
h ( 3, z , t ) = 0.65 − z
where q(x,z,t) is the ﬂux at spatial coordinates x and z at time t.
The Brooks–Corey exponent ε was set equal to 6.37 in the UZF1
model, based on the estimate from Carsel and Parrish (1988). A
zero ET extinction depth was assumed, as the experiment was
designed to minimize all evaporative losses. The HYDRUS package used retention curve parameters similar to those used for the
VSF process (Thoms et al., 2006; Vauclin et al., 1979), i.e., θs =0.30, θr = 0.01, α = 0.033 cm−1, and n = 4.1.
Figure 5 compares water tables simulated using the HYDRUS
and UZF1 packages with the experimental data. One may also
compare the performance of the VSF process by referring toThoms et al. (2006). Water tables calculated using the HYDRUS
package are similar to those simulated using the VSF processeven though the numerical solution of the Richards equation in
the HYDRUS package is limited to only the vertical direction. It
was observed that the one-dimensional nature of the vadose zone
modeling used in the HYDRUS package did not signiﬁcantly
aﬀect the correspondence of simulated results with experimental data. Note that while only vertical ﬂow was allowed in thevadose zone, horizontal ﬂow below the water table redistributed
recharged water and resulted in smooth water tables; however, acomparison of results calculated with the UZF1 package (Fig.
5b) with those obtained using the VSF and HYDRUS packages
shows that the UZF1 package marginally underestimated thedepth of the water table. This may be attributed to the kinematicwave approximation used in the UZF1. Also, the uncertainty insoil hydraulic parameters may be responsible for the diﬀerences.
The underprediction of the water table depth by the UZF1 package is relatively larger at later times (i.e., 8 h). The calibrated
UZF1 package would probably provide similar results to thoseby HYDRUS and VSF.
A comparison of computational times needed for the VSF,HYDRUS, and UZF1 showed that the UZF1 package required
the least computational eﬀort, followed by HYDRUS and VSF.
Simulation times required by the HYDRUS and UZF packages were,
however, signiﬁcantly smaller than required by the VSF process.
It can be concluded that for small-scale groundwater problems (up to two dimensions, such as in Case Studies 1 and 2)
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

F . 5. A comparison of water table posi ons calculated using the(a) HYDRUS and (b) Unsaturated Zone Flow (UZF1) packages withthe experimental data of Vauclin et al. (1979).

with downward ﬂow in the vadose zone, the UZF1 package has
accuracy similar to the HYDRUS package and the VSF process,at least for certain cases such as those considered here. The thirdcase study was designed to test the performance of the HYDRUS
package for a regional-scale groundwater ﬂow problem.
Case Study 3: Hypothe cal Regional-Scale Groundwater Problem

The third case study involved a hypothetical large-scale
groundwater ﬂow problem in a semiarid to arid region. The
geometry of the modeling domain (Fig. 6) was based on thetest example described in Prudic et al. (2004) and Niswonger etal. (2006). In this case study, we compared the eﬀectiveness of
the REC-ET, UZF1, and HYDRUS packages in characterizing
vadose zone processes at a regional scale. The VSF process was
not used here because of its extraordinary computational demand(Thoms et al., 2006) for such large-scale applications.
The model domain was designed to represent an alluvial basinwith loamy soils. Figure 6 shows the model domain and other key
characteristics for this hypothetical regional-scale groundwaterﬂow problem. The ﬂow domain was divided into uniform gridsof 1524- by 1524-m size. Two cells were assigned a general headboundary condition to simulate head-dependent ﬂux boundaries
to allow ﬂow in and out of the system. At the head-dependent763

F . 7. (a) Land surface eleva on, (b) depth to bedrock, and (c)water table depth at the beginning of the simula on for the hypothe cal regional-scale groundwater flow problem.
F . 6. Model domain, spa al distribu on of hydraulic conduc vi es
and specific yields, wells (red circles), and general head boundaries
for the hypothe cal regional-scale groundwater flow problem.

We assumed that the vadose zone consisted of loamy soils
throughout the model domain. The following Brooks–Corey(Brooks and Corey, 1964) model parameters were used in the
UZF1 package, as suggested for loam by Carsel and Parrish(1988): θs = 0.30, θr = 0.00, Ks = 3 × 10−5 cm s−1, inverse of the
air entry pressure α = 0.0897 cm−1, the pore size distribution
index n = 0.22, and the pore-connectivity parameter l = 1. An ET
extinction depth of 2.65 m (Shah et al., 2007) was used and the
Brooks–Corey exponent, ε, was calculated to be 12.09.
The HYDRUS package offers a variety of models forcharacterizing the water retention and unsaturated hydraulicﬂux boundary, water enters the model domain if the head in thecell is less than a certain user-deﬁned reference head and leaves
the model domain otherwise. The alluvium valley aquifer was
assumed to have greater hydraulic conductivity than the uplandareas. Figure 7 shows surface elevations, bedrock depth, and initialwater table depths in the study area. The geologic settings werevaried spatially to present a complex three-dimensional case.
The modeling time was divided into 12 equal stress periods,
each of which lasted 30.42 d. Except for the ﬁrst stress period,
they were modeled in MODFLOW in the transient mode with
15 time steps for each stress period. The ﬁrst stress period wasmodeled as steady state. The meteorological conditions wereassigned to represent a semiarid climate where potential evaporation rates are substantially higher than precipitation rates. Suchmeteorological settings require consideration of both downward
and upward water ﬂuxes in the soil proﬁle and provide, therefore,a good case study to compare the HYDRUS package to the UZF1
package. While Table 2 provides the base precipitation, potential
evaporation, and well pumping rates for the 12 stress periods,
Fig. 8 shows precipitation rate factors that were used to varyprecipitation rates throughout the ﬂow domain. While diﬀerent
precipitation rates were assigned for each stress period, the spatial
distribution of precipitation rates was considered to be the samefor all stress periods.

www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

T
2. Precipita on, poten al evapora on, and pumping rates
for a hypothe cal regional-scale groundwater flow problem.
Stress period

1
2
3
4
5
6
7
8
9
10
11
12

764

Precipita on rate

Pumping rate

Poten alevapora on rate

mm d−1
0.03
0.21
0.26
0.53
0.13
0.26
0.08
0.13
0.08
0.05
0.03
0.03

m3 s−1
−0.06
−0.06
−0.03
−0.03
−0.08
0.00
−0.06
0.00
−0.03
0.00
−0.08
0.00

mm d−1
1.32
1.32
1.32
1.32
1.32
1.32
1.32
1.32
1.32
1.32
1.32
1.32

F . 8. Zona on showing the spa al distribu on of precipita on
for the study area of the hypothe cal groundwater flow problem.
For any stress period, the actual precipita on rate in the zone is
obtained by mul plying the precipita on rates given in Table 2 by
the zone precipita on rate factors.

F . 9. MODFLOW zones used to define HYDRUS soil profiles in the
hypothe cal groundwater flow problem.
F . 10. Ini alwater contents
as a func on of
depth in HYDRUS
soil profiles represen ng selected
MODFLOW zones(colors correspondto zones in Fig. 9)
in the hypothe calgroundwater flow
problem.

conductivity of the soil. The van Genuchten–Mualem model(van Genuchten, 1980) was used in the HYDRUS package torepresent the unsaturated hydraulic conductivity and water content dependency on the capillary pressure. As suggested by Carsel
and Parrish (1988), the following van Genuchten parameters forloam were used: θs = 0.30, θr = 0.00, Ks = 3 × 10−5 cm s−1, α =0.036 cm−1, n = 1.56, and l = 0.5.
For both HYDRUS and UZF1 packages, it was essential toﬁrst create zones that represented relatively homogeneous units
with similar soil and hydrogeologic properties so that one soil proﬁle could be assigned to each zone. To create the zones, the fuzzy
c-means clustering algorithm of Dunn (1973) was used. The fuzzy
c-means algorithm is a method of clustering that allows one pieceof data to be split into a user-deﬁned number of clusters such thatthe data points in each cluster are as similar as possible. The fuzzy
c-means method was used independently of the MODFLOW–HYDRUS environment as it is not a part of it. The zones werecreated based on surface elevations, hydraulic conductivities,
initial water table heads, and locations of cells. The number ofclusters was chosen to be 20 since additional clusters did not signiﬁcantly improve uniformity within each cluster (MODFLOWzone). Figure 9 shows the zones used in the HYDRUS package
for the hypothetical vadose zone–aquifer interaction problem. The
same zones were used for the UZF1 package.
The HYDRUS soil proﬁles corresponding to each of thesezones were then created and discretized vertically into ﬁnite elements. Figure 10 shows initial soil water contents for selectedsoil proﬁles. One may note a large variation of initial soil watercontents between diﬀerent soil proﬁles even though the same soil
texture was used throughout the transport domain (i.e., loamy
soils). This was due to diﬀerences in other variables, such as thewater table depth and the elevation of the land surface. In otherpackages, such as REC-ET and UZF1, it is diﬃcult to take into
account the variability of soil proﬁles.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

One of the beneﬁts of the HYDRUS package is its ability
to estimate transient ﬂuxes in the groundwater table based on
the hydrologic conditions and water residence time in the soilcolumn without signiﬁcantly compromising the computational
requirements. Changes in the magnitude and direction of theﬂux can be caused by changes in water contents in the HYDRUS
soil proﬁles as a result of time-variable surface meteorologicalconditions and the position of the water table. Figure 11 shows
the groundwater ﬂux zones estimated for diﬀerent stress periods.
Note that, depending on various hydrologic and topological conditions, the HYDRUS package predicts both positive (downward)
recharge and negative (upward, capillary rise) discharge ﬂuxes.
Water table ﬂuxes at any cell are directly inﬂuenced by surface
inﬁltration, evaporation, and transpiration, as well as pumping
rates in and around a particular cell. During the initial stressperiods, the HYDRUS package predicted considerable upward
ﬂuxes, especially in cells where the precipitation rates were lower765

F . 11. Groundwater table fluxes (recharge vs. discharge) as predicted by the HYDRUS package at the end of Stress Periods (a) 3,(b) 6, and (c) 12 for the hypothe cal groundwater flow problem.

F . 12. Depths to the water table for the groundwater flow
problem at the end of Stress Period 12 as es mated by the (a)Recharge–Evapotranspira on (REC-ET), (b) Unsaturated Zone Flow(UZF1), and (c) HYDRUS packages.

(compare Fig. 11 with Fig. 8). Cells with deeper initial watertables were less aﬀected by evaporation than those with shallow
initial water tables. As the simulation time proceeded, the numberof cells with upward ﬂuxes decreased because of the inﬁltration
front movement toward the water table.
Figure 12 compares the ﬁnal water table depths estimatedby the three packages. Water table depths estimated by the
UZF1 and HYDRUS packages were comparable except that the
HYDRUS package consistently predicted marginally deeper watertables. Water table depths predicted using HYDRUS at the end of
the simulation were deeper by between 0 and 1.31 m than those
calculated with the UZF1 package. This is probably due to a moreaccurate characterization of capillary pressures and ﬂuxes in thevadose zone by HYDRUS. The HYDRUS package can consider
upward pressure gradients that cannot be simulated by either the
UZF1 or the REC-ET packages, both of which tend to predictgreater downward water table ﬂuxes.
Figure 13 shows water table depths at the end of diﬀerent
stress periods as a function of initial water table depths. Thisﬁgure indicates the impact of vadose zone ﬂow on model predictions of groundwater tables. Low capillary pressures in soil proﬁlesof arid zones often lead to slow upward water movement. Finalwater table depths predicted using the HYDRUS package deviate
from those predicted by the other packages, especially for cellswith deep initial water tables. One may attribute these diﬀerences
to the inﬂuence of capillary forces on vadose zone ﬂow, which

are fully considered only in the HYDRUS package. Groundwatertable depths calculated using the UZF1 package also deviate fromthose calculated using the REC-ET package, albeit in some cells
to a lesser extent than those predicted using the HYDRUS package. This is due to the delay in groundwater table ﬂuxes resulting
from the gravity-driven ﬂow considered by the UZF1 package
and neglected by the REC-ET package. One can infer that vadosezone processes resulting from diﬀerent soil hydraulic characteristics can be better represented using the HYDRUS package thanthe REC-ET and UZF1 packages, especially for soils, such as
medium- and ﬁne-textured soils, where capillary forces play animportant role in determining water ﬂuxes.
The computational time needed to perform calculations withthese vadose zone packages is of paramount importance because it
seriously aﬀects their applicability to regional-scale problems. The
UZF1, REC-ET, and HYDRUS packages were run on a 1 GB
RAM, 3.40 GHz Intel Pentium based personal computer. The
simulation of the hypothetical regional-scale groundwater model(Case Study 3) using the REC-ET package needed, as expected,
the least amount of time (5 s). The MODFLOW simulations with
the UZF1 and HYDRUS packages took approximately 20 and
26 s, respectively. While the computational demand of the UZF1
and HYDRUS packages is comparable, the HYDRUS package
provides a more comprehensive characterization of vadose zone
ﬂow processes.

www.vadosezonejournal.org · Vol. 7, No. 2, May 2008766

F . 13. Depth to the water table es mated using the Recharge–Evapotranspira on (REC-ET), Unsaturated Zone Flow (UZF1), and HYDRUS
packages at the end of Stress Periods (a) 3, (b) 6, (c) 9, and (d) 12 as a func on of the ini al water table.

Summary and Conclusions

application. Although previous research (e.g., Shah et al., 2007)
provides ET extinction depths for diﬀerent textural classes, determination of the ET extinction depth for practical cases can becumbersome. The presence of vegetation further complicatesthe determination of the ET extinction depth. In contrast, the
HYDRUS package and the VSF process use soil hydraulic andplant parameters, such as the van Genuchten and Feddes parameters, respectively, that are readily available for a wide variety of soil
textures and plants (Leij et al., 1996; Lilly, 1997). The HYDRUS
package also oﬀers a variety of models for describing soil hydraulic properties. The necessary parameter values required for thesemodels in HYDRUS may be estimated experimentally or using
pedotransfer functions (e.g., Schaap et al., 2001).
Another weakness of the current UZF1 package is that it simulates unsaturated ﬂow only for homogeneous vadose zones andit cannot consider multiple soil horizons with varying hydraulicproperties. User-speciﬁed layering can be easily accommodated
when using the HYDRUS package. The HYDRUS package thus
oﬀers a good alternative to the UZF1 package when these factors,e.g.., vegetation or multiple horizons, are signiﬁcant for a particular application. The HYDRUS package thus may be expected
to perform better for regional-scale groundwater problems

We evaluated the recently developed HYDRUS package
for MODFLOW and compared its performance with other
MODFLOW packages (REC-ET, UZF1, and VSF) that accountfor processes in the vadose zone. Based on the HYDRUS-1D
software, the HYDRUS package considers the eﬀects of inﬁltration, soil moisture storage, evaporation, plant water uptake,precipitation, runoﬀ, and water accumulation at the ground surface. The HYDRUS package was compared with other currently
available packages (VSF and UZF1) for case studies of varying
complexities. For smaller scale problems (up to two dimensions),
the VSF process and the UZF1 and HYDRUS packages perform
similarly. The VSF process provides the best accuracy due to its
thorough consideration of vadose zone processes and is most suitable for smaller scale problems. For three-dimensional models,the HYDRUS package demonstrated a signiﬁcant improvement
in modeling accuracy compared with the UZF1 package and asigniﬁcant decrease in computational demand compared with
the VSF process. This new package represents a promising toolaccounting for the eﬀects of vadose zone processes on groundwater levels and ﬂuxes.
Even though the UZF1 package provides comparable results
to the HYDRUS package, a number of diﬃculties arise during its
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008767

water recharge. Hydrogeol. J. 10:18–39.
Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. Rosetta: A computer
program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol. 251:163–176.
Seo, H.S., J. Šimůnek, and E.P. Poeter. 2007. Documentation of the HYDRUS
package for MODFLOW-2000, the U.S. Geological Survey modular
ground-water model. GWMI 2007-01. Int. Ground Water Modeling Ctr.,
Colorado School of Mines, Golden.
Shah, N., M. Nachabe, and M. Ross. 2007. Extinction depth and evapotranspiration from ground water under selected land covers. Ground Water
45:329–338.
Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2005. The HYDRUS-1D
software package for simulating the one-dimensional movement of water,
heat, and multiple solutes in variably-saturated media, Version 3.0. HYDRUS Softw. Ser. 1. Dep. of Environ. Sci., Univ. of California, Riverside.
Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2008. Development andapplications of the HYDRUS and STANMOD software packages and related codes. Vadose Zone J. 7:587–600 (this issue).
Singh, V.P. 2002. Is hydrology kinematic? Hydrol. Processes 16:667–716.
Smith, R.E. 1983. Approximate sediment water movement by kinematic characteristics. Soil Sci. Soc. Am. J. 47:3–8.
Smith, R.E., and R.H.B. Hebbert. 1983. Mathematical simulation of interdependent surface and subsurface hydrologic processes. Water Resour. Res.
19:987–1001.
Thoms, R.B., R.L. Johnson, and R.W. Healy. 2006. User’s guide to the Variably
Saturated Flow (VSF) process for MODFLOW. Techniques and Methods
6-A18. Available at http://pubs.usgs.gov/tm/2006/tm6a18/ (veriﬁed 15
Apr. 2008). USGS, Reston, VA.
Uddameri, V., and M. Kuchanur. 2007. Estimating aquifer recharge in Mission
River watershed, Texas: Model development and calibration using genetic
algorithms. Environ. Geol. 51:897–910.
van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.
van Genuchten, M.Th., S.M. Lesch, and S.R. Yates. 1991. The RETC code forquantifying the hydraulic functions of unsaturated soils. Version 1.0. U.S.
Salinity Lab., Riverside, CA.
Vauclin, M., J. Khanji, and G. Vachaud. 1979. Experimental and numerical
study of a transient, two-dimensional unsaturated–saturated water table
recharge problem. Water Resour. Res. 15:1089–1101.
Ward, S. 2002. Recharge and groundwater models: An overview. Hydrogeol. J.
10:110–120.
Wierenga, P.J., R.G. Hills, and D.B. Hudson. 1991. The Las Cruces trench site:
Characterization, experimental results, and one-dimensional ﬂow predictions. Water Resour. Res. 27:2695–2705.
Winter, T.C., J.W. Harvey, O.L. Franke, and W.M. Alley. 1998. Ground waterand surface water a single resource. Circ. 1139. USGS, Denver, CO.

with complex layering in the vadose zone and with alternating
recharge–discharge ﬂuxes.
The HYDRUS package currently simulates only water ﬂow
and is distributed as an open-source code. We intend to expand
the HYDRUS package to also simulate solute transport so that
the MODFLOW–HYDRUS code will produce concentrations
as a function of time that can be incorporated into the source
function for MT3D. This is expected to be especially useful forregional-scale studies involving nonpoint-source pollution.

References
Brooks, R.H., and A.T. Corey. 1964. Hydraulic properties of porous media. Hydrol Pap. 3. Colorado State Univ., Fort Collins.
Carsel, R.F., and R.S. Parrish. 1988. Developing joint probability distributionsof soil water retention characteristics. Water Resour. Res. 24:755–769.
Colbeck, S.C. 1972. A theory of water percolation in snow. J. Glaciol. 2:369–385.
Dunn, J.C. 1973. A fuzzy relative of the ISODATA process and its use in detecting compact well-separated clusters. J. Cybernetics 3:32–57.
Feddes, R.A., P.J. Kowalik, and H. Zaradny. 1978. Simulation of ﬁeld water use andcrop yield. Simulation Monogr. PUDOC, Wageningen, the Netherlands.
Gee, G.W., and D. Hillel. 1988. Groundwater recharge in arid regions: Review
and critique of estimation methods. Hydrol. Processes 2:255–266.
Harbaugh, A.W., E.R. Banta, M.C. Hill, and M.G. McDonald. 2000. MODFLOW-2000, the U.S. Geological Survey modular ground-water model
user guide to modularization concepts and the ground-water ﬂow process.
USGS, Denver, CO.
Hendrickx, J., and G. Walker. 1997. Recharge from precipitation. p. 19–38. In
I. Simmers (ed.) Recharge of phreatic aquifers in (semi-) arid regions. A.A.
Balkema, Rotterdam, the Netherlands.
Keese, K.E., B.R. Scanlon, and R.C. Reedy. 2005. Assessing controls on diffuse groundwater recharge using unsaturated ﬂow modeling. Water Resour.
Res. 41:W06010, doi:10.1029/2004WR003841.
Leij, F.J., W.J. Alves, M.Th. van Genuchten, and J.R. Williams. 1996. The UNSODA—Unsaturated Soil Hydraulic Database: User’s manual version 1.0.
Rep. EPA/600/R-96/095. USEPA Natl. Risk Manage. Res. Lab., Cincinnati, OH.
Lerner, D.N., A.S. Issar, and I. Simmers. 1990. Groundwater recharge: A guideto understanding and estimating natural recharge. Rep. 8. Int. Assoc. Hydrogeol., Kenilworth, UK.
Lilly, A. 1997. A description of the HYPRES database (Hydraulic Properties ofEuropean Soils). p. 161–184. In A. Bruand et al. (ed.) The use of pedotransfer functions in soil hydrology research. Proc. Worksh. of the Project“Using Existing Soil Data to Derive Hydraulic Parameters for Simulation
Modelling in Environmental Studies and in Land Use Planning,” 2nd, Orleans, France. 10–12 Oct. 1996. INRA, Orleans.
Lorenz, D.L., and G.N. Delin. 2007. A regression model to estimate regional
ground water recharge. Ground Water 45:196–208.
Morel-Seytoux, H.J., P.D. Meyer, M. Nachabe, J. Touma, M.Th. van Genuchten,
and R.J. Lenhard. 1996. Parameter equivalence for the Brooks–Corey andvan Genuchten soil characteristics: Preserving the eﬀective capillary drive.
Water Resour. Res. 32:1251–1258.
Niswonger, R.G., D.E. Prudic, and R.S. Regan. 2006. Documentation of theUnsaturated-Zone Flow (UZF1) package for modeling unsaturated ﬂow
between the land surface and the water table with MODFLOW-2005.
Techniques and Methods 6-A19. Available at http://pubs.usgs.gov/
tm/2006/tm6a19/ (veriﬁed 14 Apr. 2008). USGS, Denver, CO.
Pikul, M.F., R.L. Street, and I. Remson. 1974. A numerical model based oncoupled one-dimensional Richards and Boussinesq equations. Water Resour. Res. 10:295–302.
Prudic, D.E., L.F. Konikow, and E.R. Banta. 2004. A new Stream-Flow Routing(SFR1) package to simulate stream–aquifer interaction with MODFLOW2000. Open-File Rep. 2004-1042. USGS, Carson City, NV.
Refsgaard, J.C., and B. Storm. 1995. MIKE SHE. p. 809–846. In V.P. Singh(ed.) Models of watershed hydrology. Water Resour. Publ., Highlands
Ranch, CO.
Richards, L.A. 1931. Capillary conduction of liquids through porous mediums.
Physics 1:318–333.
Scanlon, B.R. 2002. Choosing appropriate techniques for quantifying groundwww.vadosezonejournal.org · Vol. 7, No. 2, May 2008768