WATER RESOURCES RESEARCH, VOL. 38, NO. 12, 1323, doi:10.1029/2001WR001233, 2002

Intercode comparisons for simulating water balance of surficial
sediments in semiarid regions
Bridget R. Scanlon,1 Marty Christman,2 Robert C. Reedy,1 Indrek Porro,3 Jirka Simunek,4and Gerald N. Flerchinger5
Received 6 February 2002; revised 23 July 2002; accepted 23 July 2002; published 31 December 2002.

[1] Near-surface water balance modeling is often used to evaluate land-atmosphereinteractions, deep drainage, and groundwater recharge. The purpose of this study was tocompare water balance simulation results from seven different codes, HELP, HYDRUS1D, SHAW, SoilCover, SWIM, UNSAT-H, and VS2DTI, using 1–3 year water balancemonitoring data from nonvegetated engineered covers (3 m deep) in warm (Texas) andcold (Idaho) desert regions. Simulation results from most codes were similar andreasonably approximated measured water balance components. Simulation of infiltrationexcess runoff was a problem for all codes. Annual drainage was estimated to within ±64%by most codes. Outliers result from the modeling approach (storage routing versus
Richards’ equation), upper boundary condition during precipitation, lower boundary
condition (seepage face versus unit gradient), and water retention function (van Genuchten
versus Brooks and Corey). A unique aspect of the code comparison study was theability to explain the outliers by incorporating the simulation approaches (boundary
conditions or hydraulic parameters) used in the outlying codes in a single code andcomparing the results of the modified and unmodified code. This approach overcomesthe criticism that valid code comparisons are infeasible because of large numbers ofdifferences among codes. The code comparison study identified important factors for
INDEX TERMS: 1866 Hydrology: Soil moisture;
simulating the near-surface water balance.
1818 Hydrology: Evapotranspiration; 1875 Hydrology: Unsaturated zone; 1836 Hydrology: Hydrologic
budget (1655); KEYWORDS: hydrologic budget, soil moisture, unsaturated zone, water balance modeling
Citation: Scanlon, B. R., M. Christman, R. C. Reedy, I. Porro, J. Simunek, and G. N. Flerchinger, Intercode comparisons for
simulating water balance of surficial sediments in semiarid regions, Water Resour. Res., 38(12), 1323,
doi:10.1029/2001WR001233, 2002.

1. Introduction[2] The purpose of this study was to compare the performance of seven different codes (HELP, HYDRUS, SoilCover, SHAW, SWIM, UNSAT-H, and VS2DTI) for
simulating the near-surface water balance. These codes werechosen because they can simulate flow in response tometeorological forcing, they represent different modelingapproaches (storage routing and Richards’ equation), they
are fairly well documented, they have been widely used andtested, and most are in the public domain. The simulations
were applied to nonvegetated, engineered covers in semiarid
regions in warm (Chihuahuan Desert of Texas) and cold
(Snake River Plain, Idaho) deserts. Although simulation of1
Bureau of Economic Geology, University of Texas at Austin, Austin,
Texas, USA.
2
GeoSyntec Consultants, Austin, Texas, USA.
3
Idaho National Engineering and Environmental Laboratory, Idaho Falls,
Idaho, USA.
4
George E. Brown Jr. Salinity Laboratory, Riverside, California, USA.
5
U.S. Department of Agriculture, REE-ARS-PWA-WMR, Boise, Idaho,
USA.

Copyright 2002 by the American Geophysical Union.
0043-1397/02/2001WR001233$09.00

nonvegetated systems may be considered limiting, theadded complexities of simulating vegetation and evapotranspiration, and the lack of field measurements of vegetation parameters and processes (interception, leaf areaindex, rooting depth, transpiration) would greatly complicate the code comparison analysis. Application of thesecodes to sites using detailed field-monitored water balancedata is valuable for evaluating how well the codes simulate
field measurements. A unique aspect of this code comparison study is that it was extended beyond the traditional
comparison of code results to attribution of differences tospecific processes or parameters by incorporating the different processes or parameters into a single code in order toreproduce the results in other codes. This approach helps toovercome the criticism that valid code comparisons are
infeasible as a result of too many differences among codes
because it reduces the comparison to single issues. Uniquecontributions of this study include number of codes being
compared, detailed field monitoring of the water balance,
length of the monitoring record for simulation (1– 3 yr), anddifferent climate settings evaluated (warm and cold desert).
Although water balance-monitoring data from engineered
covers were used for input and comparison with simulation
results, results of this study are not limited to evaluation of
the engineered covers but should be applicable to general

59 - 1

59 - 2SCANLON ET AL.: INTERCODE COMPARISON

water balance modeling of nonvegetated systems. The
approaches used for comparing codes could also be usefulin any future code-comparison studies.
1.1. Background[3] Why are simulations of the water balance of surficial
sediments important? Near-surface water balance modelingis used to evaluate land-atmosphere interactions, estimate
groundwater recharge and contaminant transport, and evaluate the performance of engineered covers for waste containment. Partitioning of precipitation into the various waterbalance components is used for evaluating the regional andglobal water cycle [Milly, 1994; Milly and Dunne, 1994].
Recharge estimates based on water balance modeling are
critical for quantification of water resources and for assessing aquifer vulnerability to contamination. Many previous
modeling studies used to estimate groundwater recharge arebased on 1-D water balance codes that simulate drainage
below the root zone [Rockhold et al., 1995; Kearns andHendrickx, 1998]. Such 1-D models can be applied to areasof different climate, vegetation, and land use to determine
the impact of these factors on recharge and to develop
recharge maps [Hatton, 1998]. The impact of land-use
change can also be evaluated using these codes [Pierce etal., 1993]. The 1-D models should be appropriate in areas
where lateral flow processes are not important, i.e. noinfiltration excess runoff and low topography [Dawes etal., 1997; Hatton, 1998].
[4] Numerical water balance models are often used to
predict the effectiveness of engineered covers in minimizing
infiltration into underlying waste [Fayer et al., 1992; Khireet al., 1997; Andraski and Jacobson, 2000]. Sensitivity
analyses have been conducted to determine important factors controlling subsurface flow, such as time discretizationof precipitation input, vegetation, soil depth and texture,
hydraulic parameters, liquid and vapor flow, and hysteresis[Fayer and Gee, 1992; Meyer, 1993; Rockhold et al., 1995;
Stothoff, 1997; Kearns and Hendrickx, 1998; Andraski andJacobson, 2000]. Recent advances in computer technology,
more computationally efficient codes, and availability of
input data on climate and hydraulic properties online makelong-term simulations of the near-surface water balancemuch more feasible. Weather generators, such as USCLIMATE and GEM [Richardson, 2000] can be used todevelop long-term climate records for the simulations.
Pedotransfer functions can be used to estimate hydraulicparameters from information on soil texture [Schaap andLeij, 1998; Schaap et al., 1998].
[5] Codes used for simulating the near-surface waterbalance vary in complexity from bucket codes [Schroeder
et al., 1994; Flint et al., 2000] and simple, analytical codes[Kim et al., 1996; Simmons and Meyer, 2000] to more
complex, numerical codes [Fayer et al., 1992; Scanlon andMilly, 1994; Khire et al., 1997; Stothoff, 1997]. Bucket codes
approximate the storage capacity of the root zone by using abucket or reservoir. Deep drainage occurs when precipitation
exceeds runoff and evapotranspiration. Simple 1-D bucket
models have been used to estimate groundwater rechargeover large areas [Hatton, 1998; Flint et al., 2000] and toevaluate the importance of episodic recharge [Lewis andWalker, 2002]. Those interested in global climate change
have been developing simplified analytical solutions ofRichards’ equation in models of near-surface processes toincorporate into global climate models and to assess largescale water balance and ecohydrologic issues [Milly, 2001;
Rodriguez-Iturbe et al., 2001]. More complex numericalsolutions of Richards’ equation are used for water balance
modeling where detailed information on hydraulic parameters is available [Fayer et al., 1992; Flerchinger et al., 1996;
Khire et al., 1997]. The computational requirements ofnumerical models are much greater than those of analyticalmodels.
[6] The reliability of model results can be evaluated (1) bycomparing results obtained from different codes and (2) bycomparing model results with field monitored water balancedata. Programs such as nuclear waste disposal have included
a wide variety of code comparison studies (e.g., INTRACOIN, 1981 – 1986; HYDROCOIN; 1984 – 1989; andINTRAVAL, 1987 – 1993] for evaluating geosphere transportmodels for radioactive substances in the subsurface [Larssen
et al., 1995]. The Society of Petroleum Engineers conducts
an annual study titled Comparative Solution Project thatcompares the performance of several simulators by using
benchmark problems [Christie and Blunt, 2001] (http://
www.spe.org/csp/). These studies provide valuable insights
into important modeling approaches that affect simulations ofsubsurface flow.
[7] Many code-comparison studies are based on hypothetical input data sets (e.g., Society of Petroleum Engineers,
INTRAVAL, Phase 1). Many previous modeling studies havebeen hypothetical and did not compare model results withfield data [Stothoff, 1997; Kearns and Hendrickx, 1998;
Khire et al., 2000]. However, some code-comparison studies
include field data for model input and for comparison withsimulation results (INTRAVAL, Phase 2). One of the primarysources of field measurements of near-surface water balanceis provided by monitoring of engineered covers for waste
containment [Fayer et al., 1992; Khire et al., 1997; Wilson etal., 1999]. The Alternative Cover Assessment Program(ACAP) has been set up by the Desert Research Institute(DRI, Nevada), Pacific Northwest National Laboratory (Hanford, Washington), and the Environmental Protection
Agency (EPA) to provide detailed water balance monitoringof engineered covers in different types of climate, vegetation,and soils [Wilson et al., 1999]. Limited code-comparison
studies have been conducted using water balance-monitoringdata from engineered covers [Nichols, 1991; Berger et al.,
1996; Khire et al., 1997; Ogan et al., 1999; Wilson et al.,1999]. Including field data for model input and comparison
with model results should strengthen modeling studies.
1.2. Code Descriptions[8] The codes evaluated in this study include Hydrologic
Evaluation of Landfill Performance (HELP, version 3[Schroeder et al., 1994]; http://www.wes.army.mil/el/
elmodels), HYDRUS-1D (version 3.0 [Simunek et al., 1998];
http://www.ussl.ars.usda.gov/MODELS/HYDR1D1.HTM),
SoilCover (version 4.1 [Wilson, 1990; Wilson et al., 1994;
GEO2000, 1997]; http://www.geo2000. com), Simulation of
Heat and Water (SHAW, version 2.4 [Flerchinger and Saxton,
1989; Flerchinger et al., 1996]; http://www.nwrc.ars.usda.
gov/models/shaw/index.html), Soil Water Infiltration andMovement model (SWIM, version 2 [Verburg et al., 1996];
http://www.clw.csiro.au/products/swim),UNSAT-H (version
3.0 [Fayer, 2000]; http://hydrology.pnl.gov/resources.asp),

59 - 3SCANLON ET AL.: INTERCODE COMPARISON
Table 1. Attributes of Various Codes Compared in this Studya
Model

PE Input

PE Calc.

UNSAT-HHYDRUS-1D
SHAW
SoilCover
SWIM
VS2DTIHELP

XXPenman*PM
Penman*Snow

RunoffLower b.c.

Water Retention

K(q,h) Function

I.C.

GUI

VG, BC
VG, BCBC*
Fredlund
VG, BC
VG, BCBC

MualemMualemBurdineMualemMualemMualemBurdinehq or h
h
h
h
h
q

X
X

SCS

UG
UG, S
UG
h or q(c)
UG, S
UG, S
UG

X

X
X
X
Penman

X

X
X
X

a

PE, potential evaporation; calc., calculated internally in the code; b.c., boundary condition; GUI, graphical user interface; SCS, Soil Conservation
Service curve number; UG, unit gradient; S, seepage face; BC, Brooks and Corey; VG, van Genuchten; BC*, Brooks and Corey with zero residual watercontent (equivalent to Campbell equation); q, water content; h, pressure head; PM, Penman Monteith; Penman*, modified Penman equation [Doorenbos
and Pruitt, 1977]. HYDRUS-1D and VS2DTI were modified to distribute PE daily using a sine function similar to that used in UNSAT-H and described byHillel and van Bavel [1976].

and Variably Saturated 2 Dimensional Transport Interface(VS2DTI, version 1.1 [Healy, 1990]; http://water.usgs.gov/
software/vs2di.html). Graphical user interfaces are available
for most codes (Table 1). SHAW has an interface for data
input but does not have a postprocessor. Detailed descriptionsof the codes can be found in the user’s manuals listed above.
It is impossible to describe all the attributes of the various
codes in this paper; however, many of the attributes related
to water balance modeling of HELP, HYDRUS, SoilCover,and SHAW were described by Wilson et al. [1999].
[9] The code comparison focuses on various componentsof the water balance equation
P þ Irr  E  R0  D ¼ S

ð1Þwhere P is precipitation, Irr is irrigation, E is actual
evaporation, R0 is surface runoff, D is drainage, and S ischange in water storage. All codes evaluated in this study
simulate subsurface water flow using Richards’ equation,with the exception of HELP, which uses a storage-routing or
water balance approach. SHAW uses the Green and Ampt
approach to simulate infiltration but uses Richards’ equation
to simulate redistribution. All of the codes simulate liquid
flow; however, SoilCover and SWIM also simulate isothermal vapor flow, and HYDRUS-1D, SHAW, and UNSATH simulate isothermal and thermal vapor flow. All codes use
some type of adaptive time stepping approach that allows thetime step size to increase when the code converges readily
and decrease when there are convergence problems considering a number of prespecified constraints. Initial time step
size is generally a fraction of a second, and time step size canincrease to a prespecified maximum (e.g., hr, day).
[10] HELP and SHAW are the only codes that simulatesnow accumulation and melt. If air temperature is below 0,
precipitation accumulates as snow. Snowmelt is simulated
when the temperature rises above freezing. Runoff is simulated explicitly in HELP by using an empirical relationshipbetween precipitation and runoff (curve number) developedby the U.S. Soil Conservation Service [U.S. Department ofAgriculture, 1985]. In all the other codes, runoff occurs
when precipitation intensity exceeds infiltration capacity ofthe soils. Water that has infiltrated can move up by evaporation or down as a result of gravity or matric-potential
gradients. The upper boundary condition can be specified
using meteorological data from which the codes can calculate potential evaporation (PE) internally using original or amodified version of the Penman [Penman, 1948] equation(Table 1). Alternatively, precipitation and PE can be input

directly into many codes (Table 1). PE is not used in SHAW;
instead the energy budget equation is solved to simulate
actual evaporation.
[11] The codes use a variety of different approaches tosimulate evaporation. Potential evaporation is controlled byatmospheric conditions whereas actual evaporation is limitedby the rate at which soil can transmit water upward to the land
surface. HELP uses a two-stage modified Penman energybalance method [Ritchie, 1972]. Stage1 evaporation is equal
to potential evaporation (PE) and continues until the flux
reaches a value equal to the hydraulic conductivity at a matricpotential head of 1 m. Stage 2 evaporation follows and is
calculated as a function of the number of days since stage 1
ended and approaches 0 when the wilting point (153 m) isreached. Evaporation is assumed to occur within the evaporative zone depth specified by the user. HYDRUS-1D and
UNSAT-H allow evaporation to occur at the potential rate
when the head at the surface node is between 0 and aprespecified lower value. When the head reaches the lower
bounding value, the boundary condition changes from a
constant flux (PE) to a constant head, and evaporation iscontrolled by the rate at which water can be transmitted to thesoil surface. SWIM and SoilCover calculate evaporation (E)as a sink term using the following equation [Campbell,
1985]:
E ¼ PERHs  RHa
1  RHa

ð2Þ

where RH is relative humidity and the subscripts represent
soil surface (s) and atmosphere (a). The evaporation rate isequal to PE when RHs is 1 and is equal to 0 when RHs is equal
to RHa. Simulating evaporation as a sink allows infiltrationand evaporation to occur simultaneously, and evaporation
can be accounted for at the time discretization used internallyin the code rather than that of the precipitation input(generally daily or hourly). The value of RHs is calculated
as a function of matric potential at the surface node. The
SWIM user manual [Verburg et al., 1996] indicates that thesimulations are not very sensitive to variations in RHa andsuggests a value of 0.5. SoilCover also includes a more
complicated approach for calculating evaporation that was
not used in this study [Wilson, 1990]. Evaporation in VS2DTI
is computed as an upward flux driven by the pressure
potential gradient between soil and atmosphere:
E ¼ Kuha  hL

ð3Þ

59 - 4SCANLON ET AL.: INTERCODE COMPARISON

where Ku is the unsaturated hydraulic conductivity of thesurface soil, ha is the pressure potential of the atmosphere, h is
the pressure potential of the top node, and L is the distance
between the top node and the soil surface. VS2DTI
documentation [Lappala et al., 1987] suggests a value of1,000 m for ha and also states that evaporation is not very
sensitive to variations of ha. The user must also enter a record
for PE to serve as the upper limit for evaporation (i.e., if thecalculated upward flux is greater than PE, evaporation is set
equal to PE).
[12] The lower boundary condition available in most codes
includes a unit gradient option, which allows water to drain
when it reaches the boundary (Table 1). Many codes also
include a seepage face option, which only allows drainage
when the boundary becomes saturated (HYDRUS-1D,
SWIM, and VS2DTI). This boundary condition is required
to simulate lysimeters that do not have wicks to allow
drainage under unsaturated conditions. Wickless lysimeters
require saturation for drainage to occur. SoilCover only
simulates a constant pressure or constant water content
lower-boundary condition. Initial conditions for most codes
consist of matric potentials, except for the HELP code, which
requires water content data.
[13] Many codes (HYDRUS-1D, SWIM, UNSAT-H, andVS2DTI) include multiple analytical functions for water
retention and unsaturated hydraulic conductivity. Most codes
include the van Genuchten water retention function [vanGenuchten, 1980] and the Mualem hydraulic conductivity
function [Mualem, 1976]. SHAW uses Campbell’s water
retention function [Campbell, 1974, 1985], which is similar
to the Brooks and Corey function [Brooks and Corey, 1966]
except that the residual water content is zero, and the Burdine
hydraulic conductivity function [Burdine, 1953]. HELP uses
porosity, wilting point, field capacity, and saturated hydraulicconductivity data to calculate the Brooks and Corey water
retention and Burdine hydraulic conductivity functions. SoilCover requires tabulated water retention data and saturated
hydraulic conductivity, which are fitted to water retention andunsaturated hydraulic conductivity functions developed byFredlund et al. [1994] and Fredlund and Xing [1994].
1.3. Site Descriptions[14] The Texas site is located near Sierra Blanca, which isabout 120 km southeast of El Paso, Texas. The site is within
the Chihuahuan Desert of Texas. Long-term (30-yr) meanannual precipitation at Sierra Blanca is 320 mm. The site
consists of heavily instrumented engineered covers that were
installed in the summer of 1997. The surface dimensions of
the engineered cover are 34  17 m. Subsurface water
movement was restricted primarily to the upper 2 m of the
profile, which consists of 0.3 m of topsoil (sandy clay loamwith 15% gravel), underlain by 1.7 m of compacted native
material (sandy clay loam) (Figure 1). This material is in turnunderlain by a capillary barrier that consists of sand andgravel. The various components of the water balance equation are monitored, except evaporation, which is calculatedby difference. The monitoring system consists of an onsite
meteorologic station that monitors precipitation, wind speed,
relative humidity, temperature, and solar and net radiation; asurface-runoff system that includes interception drains,
underground storage tanks and pressure transducers; watercontent monitoring using time domain reflectometry (TDR)

in the topsoil, 10 vertical neutron probe access tubes to 2 mdepth, vertical profiles of heat dissipation sensors (depths:
0.15, 0.3, 0.6, 0.9, 1.2, 1.5, and 2 m) to monitor matricpotentials; and a 3-m-deep, wickless lysimeter to monitor
drainage. TDR signals below the topsoil were generally
attenuated because of high salinity and high clay content.
Water storage of the engineered cover was calculated byintegrating the water content data averaged from all 10
neutron probe access tubes with depth. The site was irrigatedin August and September 1998 to establish vegetation byusing a 0.3-m grid drip-irrigation system. Seedlings of five
perennial bunchgrass species were transplanted on a nominal
0.8-m grid. The vegetation is expected to have a negligible
effect on the water balance in September and October 1998
because time was insufficient for roots to be established. The
application rate was approximately 17 mm/hr. Irrigationwater was applied between 2:00 and 5:00 a.m. and resulted
in an average of 5.7 mm per application. Daily irrigations
from 5 – 31 August resulted in 150.9 mm of water applied,
and alternate day irrigations from 3 – 27 September resulted
in 74.8 mm of water applied.
[15] The Idaho site is located at the Idaho National Engineering and Environmental Laboratory in southeastern Idaho[Porro, 2001] on the Snake River Plain. Long-term (40-yr)
mean annual precipitation is 221 mm. The site consists of aconcrete structure containing 10 cells, each of which is 3 m 
3 m  3 m (four walls and a floor). Replicates of a monolithic
soil cover and a capillary barrier cover were constructed in
the cells. Data from only one of the monolithic soil cover cells
are used in this study. The texture of the soil is silt loam, andthe upper 0.15 m of the profile has 25% gravel by volume
mixed with the silt loam soil to reduce wind erosion(Figure 1). Daily climate data were recorded at a NOAA
station 11 km northeast of the test facility and include daily
precipitation, solar radiation, temperature, wind speed, andrelative humidity. Various components of the water balance
were monitored, including soil water storage, using TDR
probes at approximately 0.2-m depth intervals, and drainage
from the bottom of the cell. Matric potential was monitored
by tensiometers at 0.2-m intervals from 0.4 to 1.4 m depth
and at 1.8, 2.6, and 3 m depths. Sidewalls on the test cells
prevent run-on and runoff. Actual ET was calculated from thewater balance equation (1) for both sites.

2. Model Input[16] Parameters for the various models were chosen sothat the simulations were approximately equivalent. All
codes were set up to simulate liquid flow only, except
SoilCover, which has no option for deselecting isothermal
vapor flow.
2.1. Texas Site[17] Simulations were conducted for the 1998 water year
(October 1997 through September 1998) when the site wasunvegetated. PE was calculated externally using net radiation
in a Penman-Monteith equation [Monteith and Unsworth,1990] and was used in the HYDRUS-1D, SoilCover, SWIM,UNSAT-H, and VS2DTI input (Table 1). Basic meteorological data (precipitation, wind speed, relative humidity, temperature, and solar radiation) were used in the HELP and
SHAW input. An SCS curve number of 94 was specified by
HELP (sandy clay loam soil, 2% slope, 17-m slope length)

59 - 5SCANLON ET AL.: INTERCODE COMPARISON

Figure 1. Texture profiles and initial matric potentials (y) for the (a) Texas and (b) Idaho sites.
for runoff. The lower boundary condition was specified as a
seepage face in codes that include this option (HYDRUS-1D,
SWIM, VS2DTI) and as a unit gradient (HELP, SHAW andUNSAT-H) or a constant matric-potential maintained at theinitial value (260 cm) (SoilCover) in codes that do not.
[18] The soil profile (3.05 m) was divided into six layers
representing the different materials (Figure 1). A total of 103
nodes were used to represent the profile, with nodal spacing
ranging from 0.2 cm at the soil surface and 2 cm at materialinterfaces, to a maximum of 15 cm within materials. The
original SHAW code allowed only 50 nodes and was modified to increase the number of nodes. SoilCover allows amaximum of 100 nodes, and the grid for this code differs
slightly from that used in other codes. All soil layers in the
profile were specified as vertical percolation layers in HELP.
HELP internally divided the Texas profile into 14 segments,
with the upper 7 located in the 45-cm-thick evaporative zone
(default evaporative zone thickness for the code). Initialconditions for the simulations were based on linear interpolation of matric potentials measured by the heat dissipation
sensors on 1 October 1997. The corresponding water contents were estimated for HELP using the water retentionfunction.
[19] Water retention data for soils in the upper three soil
layers were based on laboratory measurements using a hanging water column (0 to 0.02 MPa [2 m]) and a pressure-plate apparatus (0.01 to 0.5 MPa [1 to 50 m]). The
measurements of layer 1 (the topsoil) were made on the finesoil fraction and the saturated water content of this layer was
modified to account for the 15% gravel mixed with the sandy
clay loam assuming that the gravel holds negligible water(Table 2). Water retention data and saturated hydraulic conductivity (Ks) for the underlying soil layers were obtained
from various sources (layer 4 [UNSODA 2642]; layer 5[Rockhold et al., 1993]; layer 6 [UNSODA 4650]). The
UNSODA database was described by Leij et al. [1996].
Residual water contents were set to zero for the van Genuchten water retention functions (Table 2). To determine van
Genuchten shape parameters (a and n), van Genuchten water
retention functions were fitted to the measured q and h datausing Solver in Microsoft Excel (Table 2). Detailed water
retention data generated from the van Genuchten function
were used to fit Brooks and Corey water retention parameters
using Solver (Table 2, BC hb and BC l) for the SHAW code
(Figure 2).
[20] Initial simulations were conducted using representative Ks values from existing databases for layers 1 and 2(sandy clay loam [Schaap and Leij, 1998]). Field measurements of Ks in the topsoil gravelly sandy clay loam (layer 1)
using a Guelph permeameter were similar to the databasevalues, and the field measurements were used in the final
simulations (Table 2). In contrast, Ks of the subsoil com-

Table 2. Model Input Parameter Values: Layer Thickness, Particle Size Distribution, Saturated Hydraulic Conductivity (Ks) and
Saturated Water Content (qs) for All Codes; Water Content at Field Capacity (qfc) and at Wilting Point (qwp) for the HELP code; Residual
Water Content (qr) and van Genuchten Shape Parameters (a, n) for HYDRUS-1D, SWIM, UNSAT-H, and VS2DTI; and Brooks andCorey Air Entry Pressure Head (hb) and l Parameters for SHAW a
Site
Texas
Texas
Texas
Texas
Texas
Texas
Idaho
Idaho
Idaho
a

Layer Thickness, m Percent Gravel, Sand Silt, Clay
1
2
3
4
5
6
1
2
3

0.30
1.70
0.30
0.30
0.30
0.15
0.15
2.85
0.10

15, 48, 19, 18
0, 55,18, 27
0, 89, 3, 8
63, 28, 4, 5
gravel0, 92, 7, 1
25, 16, 40, 19
0, 22, 53, 25
gravel

Ks

qs

qfc

qwp

qr

41
20
639
10
159840
587
94
43
30240

0.45
0.35
0.40
0.14
0.51
0.38
0.36
0.47
0.42

0.230
0.280
0.160
0.110
0.0014
0.040
0.115
0.237
0.032

0.090
0.150
0.030
0.060
0.0001
0.002
0.043
0.019
0.013

0.000
0.000
0.000
0.000
0.000
0.000
0.035
0.015
0.005

VG a, cm1 VG n BC hb, cm BC l
0.027
0.010
0.020
0.007
10.95
0.050
0.036
0.005
4.930

1.276
1.167
1.464
1.188
1.722
1.774
1.601
2.090
2.190

Layer 3 in the Idaho profile is not real but was used in SHAW and UNSAT-H simulations to approximate a seepage face.

24.94652.07226.62073.8360.085
12.28011.200109.70
0.900

0.240
0.138
0.368
0.191
0.706
0.605
0.321
0.685
0.706

59 - 6SCANLON ET AL.: INTERCODE COMPARISON

Figure 2. Example water retention functions for (a) layer
1 at the Texas site and (b) layer 2 at the Idaho site. Curves
represent HELP (H), Brooks and Corey with zero residual
water content (BC) (used in SHAW), van Genuchten (VG),
and SoilCover (SC) water retention curves.
pacted sandy clay loam (layer 2) measured in the laboratory
using a constant head permeameter [ASTM D2434, 1994]
and a flexible wall permeameter [ASTM D5084, 1990] andin the field using a Guelph permeameter (Soilmoisture
Equipment Corp., Santa Barbara, CA) were much lower
than the database value of 9.6 cm/d. Simulations using thesemeasured Ks values did not match the measured waterbalance as well as the database values for layer 2. These
differences may be related to scaling issues as suggested by
the order of magnitude increase in hydraulic conductivity
between laboratory core measurements and field measurements using a Guelph permeameter. Therefore, in the codecomparison analysis we used the original database Ks values
for layer 2.
2.2. Idaho Site[21] Simulations were conducted for the period 21 July
1997 through 31 October 1999. PE calculated by UNSAT-H
was used as input to the other codes that had a PE input option(Table 1). Meteorologic data were input into HELP and
SHAW. The precipitation data were manipulated to accountfor snow accumulation and melting in all codes except
SHAW, which calculates these processes internally. Therefore, in all simulations, precipitation assumed to be snow
accumulated on days with an average temperature below
32F (0C), and was released as snowmelt on warmer days
using the degree-day method [Mockus, 1972; Magnuson,
1993]. PE was set equal to zero when the soil was consistently frozen below 20 cm depth on the basis of temperaturemonitoring data because soil moisture in frozen soil reduces
hydraulic conductivity and inhibits evaporation, and snow on
the surface acts as a sealant between the soil and theatmosphere. A similar approach was used by Fayer et al.
[1992] to approximate a snow cover. The lower boundary

condition was specified as a seepage face in codes that
include this option (HYDRUS-1D, SWIM, VS2DTI). A
unit-gradient lower boundary condition and a 10-cm-thickgravel layer at the base (Table 2, layer 3) were used in SHAWand UNSAT-H to approximate a seepage face although theengineered cover did not include a gravel layer at the base,unlike the Texas site. HELP simulates downward flow usinga unit gradient, which also applies to the lower boundary. A
constant matric potential (5 cm) lower boundary was
specified for SoilCover simulations. Initial conditions forthe simulations were based on linear interpolation of matricpotentials measured using tensiometers. The corresponding
water contents were calculated for HELP using the water
retention function.
[22] The soil profile (3.00 m) was divided into two layers(Figure 1b). An additional 10 cm of gravel was included at
the base for codes that did not include a seepage face. A totalof 120 nodes were used to represent the profile in all codes
except SoilCover, which has a maximum of 100 nodes. Nodal
spacing was 0.2 cm at the soil surface, 2 cm at materialinterfaces, and a maximum of 24 cm within materials. HELP
internally divided the Idaho profile into 11 segments, with theupper 7 located in the 30-cm-thick evaporative zone (default
evaporative zone thickness for code).
[23] Previous modeling studies using laboratory-measured water retention and Ks data in UNSAT-H simulations
did not adequately reproduce the measured water balancecomponents at the site [Porro and Martian, 1999; Porro,
2001]. The water application rate in the initial wetting phase(16.4 – 19.2 cm/d) generally exceeded the laboratory-measured Ks values (layer 1, 5.9 cm/d; layer 2, 8.9 cm/d). Porroand Martian [1999] used trial-and-error calibration withUNSAT-H to estimate the hydraulic parameters. In thisstudy we used van Genuchten water retention and Ks values
generated by the UNSAT-H calibration as input to othercodes that used van Genuchten parameters or to generate
other water retention functions for other codes (Table 2 andFigure 1). Because UNSAT-H was used for calibration, theperformance of this code could not be evaluated for the
Idaho site.

3. Results and Discussion[24] Various simulation results are shown in Tables 3
and 4 and in Figures 3 through 5. For ease of comparison,
the different simulations have been numbered in Tables 3
and 4, and these numbers are included in parentheses whennecessary for clarification.
3.1. Texas Site[25] Measured and simulated water balances for the Texas
site were compared for 1 October 1997 through 30 September
1998.
3.1.1. Runoff[26] A total of 6 cm of runoff was measured at the site.
Runoff occurred primarily in October 1997 following an
intense precipitation event (2.1 cm of precipitation in 15min; 1.7 cm runoff) and in August and September 1998,
when the soil cover was irrigated to establish vegetation (22.3
cm irrigation; 3.5 cm runoff). Runoff was underpredicted by
all codes when daily precipitation was input. HELP and
SoilCover predicted 1.5 and 0.9 cm of runoff, respectively,
whereas all other codes predicted zero runoff when using

59 - 7SCANLON ET AL.: INTERCODE COMPARISON
Table 3. Measured and Simulated Water Balance Components for the Texas Site (October 1997 Through September 1998)a
Code

Time

Measured
(1)
(1)
(2)
(2)
(3)
(3)
(4)
(4)
(5)
(5)
(6)
(6)
(7)
(8)
(9)
(10)
(11)
(12)
(13)

UNSAT-H (UG)
UNSAT-H (UG)
HYDRUS-1D (S)HYDRUS-1D (S)SHAW (UG)
SHAW (UG)
SoilCover
SoilCover
SWIM (S)SWIM (S)VS2DTI (S)VS2DTI (S)HELP (UG)
HYDRUS-1D (UG)
HYDRUS-1D (S, hyst.)
UNSAT-H (VS2DTI)UNSAT-H (HYDRUS-1D)
UNSAT-H (vapor)
UNSAT-H (hysteresis)

d
h
d
h
d
h
d
h
d
h
d
h
d
d
h
d
d
h
h

R0

E

0.0

32.6

0.0
0.0
0.0
0.0
0.1
0.0
0.0
4.4
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0

29.7
29.9
35.0
33.4
31.5
33.1
34.6
25.7
34.1
33.6
17.9
29.9
21.5
35.0
33.4
18.0
34.4
30.1
30.7

RMSE (E)S

RMSE (S)

4.1
2.0
2.0
1.6
1.2
0.8
0.7
1.4
3.5
1.4
1.4
7.1
2.0
5.2

6.7
6.6
1.6
3.3
4.8
3.5
2.1
6.6
2.6
3.0
18.8
6.8
14.2
1.2
3.0
18.4
2.0
6.4
5.8

D
0.0

1.7
1.7
1.7
1.1
0.7
1.0
1.4
1.0
1.5
1.5
7.1
1.9
4.5

0.3
0.3
0.0
0.0
0.3
0.2
0.0
0.0
0.0
0.0
0.0
0.0
0.9
0.4
0.4
0.3
0.3
0.3
0.3

a
R0, surface runoff; E, evaporation; RMSE, root mean square error; S, change in water storage, D, drainage; d, daily, h, hourly precipitation input; S,
seepage face and UG, unit gradient lower boundary conditions, and UNSAT-H (vapor) simulates isothermal vapor flow. The water balance components (E, S,
and D) represent results from applying net precipitation (42.7 cm precipitation minus 6 cm of runoff = 36.7 cm net precipitation). Potential evaporation =164 cm. Significant outliers are shown in bold. Code 1 (code 2): simulation using code 1 adjusted to represent code 2.

daily precipitation input. The effect of precipitation intensity
on simulated runoff was tested by using daily and hourly
precipitation in all codes that include this option (all codes
except HELP) and 15-min precipitation in HYDRUS-1D,
SWIM, and VS2DTI, which are the only codes that allow
input on any time interval. SoilCover simulated 6.2 cm ofrunoff using hourly precipitation input, which is similar to themeasured runoff (6.0 cm); however, there was little difference (<0.1 cm) in simulated runoff between daily and hourly
precipitation input for all other codes. Using 15-min precipitation input, SWIM simulated 0 runoff and HYDRUS-1Dand VS2DTI simulated 1 cm runoff (October 1997) but did
not simulate the runoff during irrigation, which was evenmore intense than represented by the 15-min precipitation
input. The temporal resolution of the precipitation measurements was 15 min from October 1997 through June 1998.
Measurements were made at much higher resolution (10 to15 s) from June through October 1998. Additional simulations with actual intensities for this time period resulted in
negligible (0 to 0.1 cm) increases in runoff. The inability tosimulate runoff may result from spatial focusing of irrigation
at the drip centers (0.3 m grid) and uncertainty in thehydraulic conductivity of the surficial sediments. The simulation results generally demonstrate the difficulties of simulating infiltration-excess runoff. Similar results were found ina study by Meyer [1993], who showed that use of daily versus
hourly precipitation underpredicted runoff by a factor of 5 at ahumid site in South Carolina characterized by high-intensity,
short-duration rainfall. One option to overcome this problem
may be to calibrate the models to runoff by decreasing thehydraulic conductivity of the surficial sediments. The
hydraulic conductivity of the topsoil had to be decreasedby an order of magnitude to simulate 6 cm of runoff at thesite using UNSAT-H.
[27] Underprediction of runoff in most codes results inmore infiltration into the system, which in turn affects theremaining terms in the water balance equation. To better

compare the simulated water balance of all codes, runoff(6 cm) was subtracted from precipitation (42.7 cm), and net
precipitation (36.7 cm) was input to the simulations. On
days when runoff was measured, the precipitation applied at
the soil surface was reduced by the amount of the measured
runoff. Using net precipitation, all codes simulated zero
runoff, except SHAW (0.2 cm; daily precipitation input) and
SoilCover (4.4 cm, hourly precipitation input) (Table 3).
Simulation of runoff using SoilCover with daily precipitation input does not suggest that this code is more reliable
than the others because SoilCover also simulated runoff
when net precipitation was input.
3.1.2. Drainage, Evaporation, and Water Storage
Change[28] No drainage was measured at the site and most codes
simulated zero drainage (Table 3). Small amounts of drainage(0.9 cm) were simulated using HELP, SHAW, and UNSATH at the start of the simulation. These drainage amounts may
reflect the wet initial conditions and nonequilibrium of theinitialization of the simulated system (profile matric potentials and boundary conditions), the simulation technique(storage routing versus Richards’ equation), and/or the useof unit gradient rather than seepage face lower boundary
conditions. HELP uses a storage routing approach and only
considers gravitational gradients and ignores matric potential
gradients, which are often upward in semiarid regions.
Therefore water flows vertically downward under unit-gradient conditions, provided that the water content in theevaporative zone is above the wilting point. The representation of the lower boundary as a unit gradient rather than aseepage face may also contribute to the excess drainage
simulated by all three codes. The effect of the lower boundary
condition was evaluated by rerunning the HYDRUS-1D
simulation with a unit gradient rather than a seepage facelower boundary condition, which resulted in 0.4 cm ofdrainage (simulation 8; Table 3), similar to the values
simulated by SHAW and UNSAT-H.

59 - 8SCANLON ET AL.: INTERCODE COMPARISON

Table 4. Measured and Simulated Annual Water Balance Components for the Idaho Sitea
Pre-WY 98
P 75.9 PE 45.8

(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
(9)
(10)
(11)
(12)
(13)
(14)
(15)
(16)
(17)
(18)

WY 98
P 24.2 PE 118.0

WY99
P 19.8 PE 129.3

Code

Time

E

S

D

E

S

D

E

S

D

Measured
UNSAT-H (UG, G), VG, Mualem
HYDRUS-1D (S), VG, Mualem
SHAW (UG, G), BC, qr = 0, Burdine
SoilCover (const. h), (Fredlund)
SWIM (S), VG, Mualem
VS2DTI (S), VG, Mualem
HELP (UG), BC
HYDRUS-1D (UG, G)HYDRUS-1D (UG)
HYDRUS-1D (hysteresis)
UNSAT-H (UG)
UNSAT-H (BC, Burdine, qr = 0)
UNSAT-H (BC, Burdine, qr 6¼ 0)
UNSAT-H (BC, Mualem, qr 6¼ 0)
UNSAT-H (VS2DTI)UNSAT-H (HYDRUS-1D)
UNSAT-H (1) (vapor)
UNSAT-H (1) (hysteresis)

d
d
d
d
d
d
d
d
d
d
d
d
d
d
d
d
d
d

5.2
5.2
10.2
9.3
9.9
9.5
3.5
10.0
10.3
9.3
10.2
4.3
6.8
6.6
9.2
3.4
9.8
5.2
5.0

48.0
48.7
48.0
48.5
49.7
51.1
51.3
7.3
47.9
11.7
46.6
10.4
48.5
46.9
44.5
48.9
48.5
48.7
40.8

22.7
22.1
17.7
18.1
16.2
15.3
21.2
58.7
17.7
78.3
19.1
82.0
20.6
22.4
22.3
23.6
17.7
22.1
30.1

15.4
13.1
18.7
19.9
16.8
16.1
10.0
15.4
19.0
14.3
18.5
9.8
19.8
19.4
26.7
8.3
16.9
13.4
12.0

0.5
0.8
0.8
0.7
1.2
0.2
0.7
1.9
1.0
13.4
0.7
11.1
1.4
1.3
4.7
0.9
0.5
0.7
0.2

8.3
10.3
6.3
5.0
8.7
7.9
13.6
7.0
6.3
23.4
6.5
25.6
5.9
6.1
2.3
15.0
7.8
10.1
12.5

11.7
10.1
13.5
22.8
11.2
10.4
9.0
13.6
13.7
8.9
13.3
6.7
16.4
16.2
23.7
7.9
11.7
10.4
9.7

0.8
2.6
1.3
3.6
1.3
1.3
2.8
2.6
1.2
1.3
3.2
5.3
5.1
5.1
7.1
2.8
2.3
2.7
3.0

8.9
12.3
7.7
0.6
10.0
10.8
13.7
8.8
7.3
12.3
9.7
18.5
8.5
8.8
3.2
14.7
10.5
12.1
13.2

a
Pre-WY98 (12 July to 30 September 1998); WY, water year; P, precipitation; PE, potential evaporation; d, daily precipitation input; E, evaporation; Schange in soil-water storage; and D, drainage; UG, unit gradient lower boundary condition; G, gravel to approximate a seepage face; S, seepage face lower
boundary condition; qr, residual water content, BC, Brooks and Corey water retention function; Burdine and Mualem are hydraulic conductivity functions.
Simulation 1 (UNSAT-H) was the calibration run used to generate hydraulic parameters for all the other codes. Simulation 8 was conducted to determine
whether a seepage face could be approximated by a unit gradient lower boundary condition with a gravel layer at the base of the profile. Simulations 9 and11 evaluate the use of a UG lower boundary condition. Simulations 10 and 18 examine the impact of hysteresis. Simulation 12 was conducted to replicate
the SHAW simulation (3). Simulation 13 evaluates the impact of the residual water content (by comparing with 12); simulation 14 examines the effect of
the hydraulic conductivity function. Simulation 15 was set up to replicate the upper boundary conditions from VS2DTI in UNSAT-H and simulation 16
replicates that of HYDRUS-1D. Simulation 17 examines isothermal vapor flow. All units in cm. Noteworthy simulation results are shown in bold. Code1(code 2): simulation using code 1 adjusted to represent code 2.

[29] The main process in the water balance at this site ispartitioning of water between evaporation and water storagechange. These two components are inversely related; i.e., if
evaporation is overpredicted, then water storage change isunderpredicted, and vice versa. The accuracy of the simulations was evaluated by calculating the root mean squareerror (RMSE) between simulated and measured storage andevaporation change. The time periods over which thechanges were evaluated correspond with the neutron probe
water content measurement intervals (monthly) (Figure 3):
RMSE ¼n
1X
ðXm  Xs Þ2i
n i¼1

!0:5
ð4Þ

where X represents evaporation or water storage change, n is
the number of time periods evaluated (15), and the

Figure 3. Time series of daily precipitation and irrigation input and measured and simulated water
storage change for the Texas site. Results from the various codes listed in order in the legend from left toright correspond to the order on the graph from top to bottom on the right. Simulation numbers inparentheses correspond to Table 3. The main outliers are VS2DTI with daily precipitation input (V(d))
and HELP (H). Results from all the other codes are similar and cannot be readily distinguished. Irrigation
was applied in August and September 1998. Error bars represent 2s uncertainties in the water content
measurements.

SCANLON ET AL.: INTERCODE COMPARISON

59 - 9Figure 4. Time series of cumulative evaporation for the Texas site. Results from the various codes listed
in order in the legend from left to right correspond to the order in the graph from top to bottom on theright. Simulation numbers in parentheses correspond to Table 3. The main outliers are VS2DTI with daily
precipitation input (V(d)), HELP (H), and SoilCover with hourly precipitation input (SC(h)). Resultsfrom all other codes are similar and cannot be distinguished readily. Error bars represent 2s uncertainties
in the precipitation, water storage, and drainage measurements.
subscripts represent measured (m) and simulated (s) values.
Accurate simulation of water storage change is importantbecause large storage increases may ultimately result indrainage.
[30] The above comparisons ignore any uncertainties in
the measurements. Uncertainties in the storage measurements were estimated considering counting uncertainties
with the neutron probe (±0.0001 cm), neutron probe calibration slope error (±0.01 cm), and spatial variability in watercontent based on the 10 neutron probe access tubes (±1.3 cm,
2s). Uncertainties in calculated evaporation incorporate
uncertainties in precipitation and irrigation input and water
storage. Errors result from partial filling of the tipping bucket(0.9 cm) and premature tipping as a result of high intensities (0.5 cm), and from irrigation application related to theflowmeter ((±0.1 cm). Some of these errors cancel and resultin a net error of +0.5 cm. Therefore the uncertainty in theevaporation estimates is 1.3, +1.8 cm.
[31] All codes except HELP and VS2DTI adequately
simulated evaporation and storage change, as indicated by
low RMS errors (Table 3 and Figures 3 and 4). Some codes
underestimated while other codes overestimated the measured water balance values. Differences between simulated
and measured annual evaporation ranged from 3% to 9% of
the measured value using daily precipitation input, whereas
differences between simulated and measured water storage

Figure 5. Time series of daily precipitation and applied irrigation and measured and simulated drainage
at the Idaho Site. Drainage curves were restarted on 1 October each year. Results from the various codes
listed in order in the legend from left to right generally correspond to the order in the graphs from top tobottom. Simulation numbers in parentheses correspond to Table 4. The plot was irrigated from 21– 25
July 1997 until breakthrough occurred. Daily irrigation amounts for this period are not shown for clarity
and were 5.4, 18.2, 17.3, 20.1, and 9.6 cm. HELP simulated 58 cm of drainage in response to thisirrigation (not shown to scale). The main outliers are HELP (H) for the first few months, SHAW (SH) for
1999, and VS2DTI (V) for 1999 and 2000. Results from all other codes are similar and cannot readily be
distinguished from measured (M) values.

59 - 10SCANLON ET AL.: INTERCODE COMPARISON

change ranged from 17% to 63%, excluding HELP and
VS2DTI. This comparison of measured and simulated
values underscores the difficulties of accurately estimatinga small number, such as water storage change or drainage,
versus a larger number, such as evaporation, because
percentage errors are amplified when a small number isestimated.
3.1.3. Explanation of Differences Between
Model Results[32] Two main outliers in the simulation results were
HELP and VS2DTI (Figure 3). RMS errors for these codes
ranged from 4.5 to 7.1 cm when daily precipitation data
were input (Table 3). Both codes underestimate evaporationand overestimate water storage change. The overestimationof water storage change by HELP is consistent with overestimation of drainage by HELP discussed previously. The
overestimation of water storage change by VS2DTI is dueto the approach it uses to assign upper boundary conditions
during precipitation. Specifically, VS2DTI sets PE to zero
during days with precipitation. When daily precipitation is
input, PE is set equal to zero for the entire day. Thisapproach underestimates evaporation and overestimateswater storage change (Table 3). The effect of the upperboundary condition is greatly reduced when actual precipitation is input at an hourly interval because the duration of
precipitation is generally 1 to 2 hr and PE can occur
throughout the rest of the day.
[33] Differences approaches for simulating the upperboundary condition with daily precipitation input are more
appropriate than that used by VS2DTI. HYDRUS-1D subtracts PE from precipitation and applies net precipitation or
net evaporation. This approach is also used in SWATRE[Belmans et al., 1983] and results in slight overestimation of
evaporation and underestimation of water storage change
when daily precipitation is input. Examples of otherapproaches include that used by UNSAT-H, in which
precipitation is applied at a specified rate (default value 1
cm/hr) until all precipitation is applied and evaporation isallowed to occur throughout the rest of the day. Simulation
results using daily and hourly precipitation input are almost
identical for UNSAT-H (Table 3), suggesting that thisapproach works well for approximating natural conditions
with daily precipitation input when there is no runoff.
SoilCover and SWIM simulate infiltration and evaporation
at the same time by simulating evaporation as a sink term
(equation 2). Results from SWIM and SoilCover are similar
to those of HYDRUS-1D. The SoilCover simulation resultswith hourly precipitation input differ from the others
because runoff was overestimated.
[34] Effects of the upper boundary condition were further
evaluated by running a single code (UNSAT-H) repeatedly,using a different upper-boundary condition for each run.
UNSAT-H simulations that were set up to approximate theVS2DTI approach (PE = 0 on rain days; simulation 10) andthe HYDRUS-1D approach (net precipitation or net PE
applied; simulation 11) generally reproduced the results from
VS2DTI and HYDRUS-1D, confirming the explanation forthe different simulated water balance components (Table 3).
Many previous studies have used daily input for precipitation
data [Fayer and Gee, 1992; Rockhold et al., 1995; Kearns
and Hendrickx, 1998], and some studies have used monthly
precipitation input [Wilson et al., 1999]. The approach usedby the various codes for simulating the upper boundary
condition when daily precipitation is input is critical because(1) only daily meteorological data are available for manysites, (2) most weather generators are limited to generating
weather data at a minimum time resolution of one day[Richardson, 1981], and (3) the feasibility of conducting
long-term (30- to 100 yr) simulations is greatly facilitatedby using daily, as opposed to shorter term, meteorologicaldata.
[35] The seasonal distributions in soil-water storage andevaporation were compared with those calculated from themeasured water balance (Figures 3 and 4). Measured water
storage increased in December 1997 because of snowmelt
and gradually decreased in the spring and early summer.
Increases in water storage in late summer resulted fromsummer rain, and large increases in August and September
resulted from irrigation and rain. Most codes generally overestimated water storage in the spring and summer, which is
attributed to underestimation in evaporation. Similar results
were found by Fayer et al. [1992] and Khire et al. [1997].
Fayer et al. [1992] improved the simulation of seasonalevaporation by decreasing the pore interaction term in the
Mualem hydraulic conductivity function from 0.5 to 0.0,
which increased evaporation in the spring and summer andnegligibly affected evaporation in the winter. The increase inevaporation was attributed to higher unsaturated hydraulicconductivity (Ku) when the pore interaction term was zero.
The impact of increased Ku on the water balance parameters
appears counterintuitive. One would expect that increasing
Ku would increase the drainage; however, Stothoff [1997]
conducted a series of simulations to demonstrate that increasing Ku allowed water to be transmitted up to the surface more
readily, which maintained a moist surface over longer time
periods that resulted in increased evaporation. These results
may be restricted to semiarid regions. Reducing the poreinteraction term in the UNSAT-H simulations of the Texas
site increased the unsaturated hydraulic conductivity by up toa factor of 3. However, this increase in unsaturated hydraulicconductivity only resulted in a negligible increase in evaporation; therefore additional studies are required to evaluate
simulation of evaporation.
3.2. Idaho Site[36] Measured and simulated water balances for the Idahosite were compared for the initial period (21 July 1997
through 30 September 1997; pre-WY98) and water years
1998 and 1999 (WY98, WY99). UNSAT-H was used tocalibrate the hydraulic properties [Porro and Martian,1999]; therefore the performance of this code was notevaluated for this site. Only daily meteorological data wereavailable for this site. Because the lysimeter cells wereenclosed on all sides, runoff could not occur. SHAW
simulated 4.9 cm of runoff in WY98 when the soil wasfrozen. A ponding depth of 10 cm was used to prohibit thiscode from simulating runoff.
[37] Most drainage from the profile (22.7 cm) occurred inJuly and August 1997 in response to irrigation in July
(Figure 5). Drainage also occurred in March/April each year
(5 – 6 cm) in response to spring snowmelt. Additional drainage (2 –3 cm) occurred in May/June each year in response tolong periods (10 d) of high precipitation. Uncertainties inthe drainage measurements are considered to be less than 1%

SCANLON ET AL.: INTERCODE COMPARISON

based on uncertainties in the tipping bucket and pressuretransducer data. Drainage toward the margins of the lysimeter
was collected separately to ensure that there was no preferential flow along the walls of the lysimeter. Simulated annual
drainage was within ±64% of measured drainage for mostcodes (simulations 1, 2, 4, 5, and 6; Table 4). The two main
outliers were HELP, which overestimated drainage by afactor of 2.6 for the initial irrigation period, and SHAW,
which generally underestimated drainage, particularly in
WY99 (factor of 14.8). The time series of drainage indicatesthat simulated drainage generally occurred much earlier than
measured drainage each year (Figure 5). The difference in thetiming between the simulated and measured drainage isattributed to inability of the codes to simulate reduction in
infiltration and in hydraulic conductivity due to soil freezing.
[38] Measured water storage increased by 48 cm after the
initial irrigation in July 1997 (pre-WY98 in Table 4). Measured water storage changes in the following water years wereminimal (i.e., WY98 and WY99 in Table 4). Uncertainties in
the water storage measurements could not be determinedbecause the TDR probes were not calibrated. The Topp
equation [Topp et al., 1980] was used to estimate water fromthe TDR data. Most codes generally simulated the annual
water storage change, with the exception of HELP, which
greatly underestimated the water storage change because it
overestimated the drainage after the initial irrigation. With the
exception of HELP, all codes simulated measured water
storage change during the first few months (July throughSeptember 1998) to within ±7%. Much larger percentage
deviations occurred in WY98 and WY99 because the measured storage values were very low.
[39] Most codes generally overestimated evaporation, withthe exception of VS2DTI, which underestimated annual
evaporation by a factor of as much as 1.5 (Table 4). SHAW
overestimated evaporation by a factor of 1.9 during water
year 1999.
3.3. Explanation of Differences Between Model Results[40] Additional simulations were conducted to further
evaluate the discrepancies in simulated results among thevarious codes. The main outlier is the HELP simulation
results for the first time period (pre-WY98) (Figure 5). The
primary difference between HELP and the other codes is theuse of storage routing versus Richards’ equation to simulate
flow. HELP uses a unit gradient to simulate downward flowin the profile when the matric potential exceeds the potential
at the wilting point. The unit gradient also applies to the lower
boundary, which forces water to drain as it reaches the lower
boundary. However, the most appropriate lower boundary
condition to represent the lysimeter is a seepage face because
it does not have wicks that would allow drainage underunsaturated conditions. The seepage face requires water tobuild up at the lower boundary until saturated conditions
develop and then drainage occurs. HYDRUS-1D, SWIM,
and VS2DTI include a seepage face option and a unit gradient
lower boundary condition with a gravel layer was used in
UNSAT-H and SHAW to approximate a seepage face.
[41] To evaluate the impact of the different lower boundary conditions, the various options were implemented in
HDYRUS-1D, i.e., a seepage face (simulation 2), unit
gradient and gravel layer (simulation 8), and unit gradient
only (simulation 9) (Table 4). Simulation results based on a

59 - 11

unit gradient and a gravel layer (8) are similar to those basedon a seepage face (2), indicating that the combination of aunit gradient lower boundary condition and a gravel layer
can be used to approximate a seepage face. In contrast,
simulation results based on a unit gradient alone (without agravel layer, 9) differ markedly from those based on a
seepage face (2). The increased drainage simulated by
HYDRUS-1D with a unit gradient lower boundary condition (9) is similar to that based on HELP (6); however, themagnitude of the drainage values simulated by HYDRUS1D exceed those simulated by HELP by as much as 19.6 cm
in pre-WY98 (Table 4). Simulation results based on
UNSAT-H with a unit gradient lower boundary condition(no gravel, 11) were similar to those based on HYDRUS-1D(9). The combination of a gravel layer and a unit gradient
lower boundary approximates a seepage face because thecapillary barrier effect of the gravel requires water to build
up almost to saturation on the gravel layer before breakthrough can occur which is similar to the saturated conditions in the seepage face. A gravel layer could not be usedin HELP to approximate a seepage face because it only
simulates flow in response to gravitational forces and not in
response to matric potential (capillary and adsorptive)
forces and cannot simulate a capillary barrier. The resultsof these simulations indicate that use of a unit gradient tosimulate a wickless lysimeter can greatly overestimate drainage from the lysimeter. The large differences in drainage
obtained using different lower boundary conditions raise thequestion of whether wickless lysimeters are suitable formonitoring near-surface water balance or performance ofengineered covers. Unless a capillary barrier underlies theengineered cover, wickless lysimeters may overestimate
water storage and, therefore underestimate drainage.
[42] The second notable outlier from the Idaho site
simulations is the SHAW simulation results. SHAW generally overestimated evaporation and underestimated water
storage change and drainage, particularly during WY99.
The upper boundary condition for SHAW assumes that all
precipitation occurs in the first hour and evaporation issimulated for the remaining 23 hr. This boundary conditionis similar to that for UNSAT-H, which assumes precipitation
occurs at a default rate of 1 cm/hr until all precipitation is
input. An additional UNSAT-H simulation was run in which
all precipitation was applied in the first hour and the results
were identical to those in which precipitation was applied at
a rate of 1 cm/hr. One of the differences between SHAW andthe other codes is that SHAW uses the Campbell water
retention function (= Brooks and Corey water retentionfunction with zero residual water content) and the Burdine
hydraulic conductivity function. Although UNSAT-H simulations incorporating similar retention and hydraulic conductivity functions as SHAW (simulation 12) were unable
to replicate SHAW simulation results exactly, they demonstrate the sensitivity of water balance simulations to water
retention functions because of their impact on unsaturated
hydraulic conductivities. UNSAT-H simulations set up to
replicate SHAW (simulation 12) differ from the UNSAT-H
base case (simulation 1) in these three factors: Brooks andCorey versus van Genuchten water retention function,
residual water content = 0 versus > 0, and Burdine versus
Mualem hydraulic conductivity function. Additional simulations were conducted to evaluate the impacts of each of

59 - 12SCANLON ET AL.: INTERCODE COMPARISON

these factors separately (Table 4). The Brooks and Coreyfunction (simulation 14) results in much greater annual
evaporation (by as much as a factor of 2) and lower
drainage (by as much as a factor of 3.6) than the van
Genuchten function (simulation 1). These differences are
attributed to the higher unsaturated hydraulic conductivity
that results from using the Brooks and Corey versus the van
Genuchten water retention functions (Figure 6). As discussed earlier, increasing unsaturated hydraulic conductivity
increases upward flow rather than drainage in these settings,
similar to the results of Stothoff [1997]. Varying the residual
water content from zero to the values for the different
materials (0.005 to 0.035 m3/m3; Table 2) had a negligible
impact on the simulation results (12 versus 13) because it
did not vary the unsaturated hydraulic conductivity. The
different hydraulic conductivity functions also impacted
simulation results (13 versus 14). The Burdine function(simulation 13) resulted in lower annual evaporation (asmuch as a factor of 1.5) and higher drainage (as much as afactor of 2.8) relative to the Mualem function (simulation
14) because of the lower unsaturated hydraulic conductivity
resulting from the Burdine function (Figure 6). This trend isin the opposite direction to the differences in evaporationand drainage between SHAW (Burdine function) and
UNSAT-H (Mualem function); therefore the main factor
controlling the difference between SHAW and UNSAT-H isthe water retention function.
[43] As at the Texas site, simulation results from the Idahosite demonstrate the importance of the approach used tosimulate the upper boundary. VS2DTI overestimates annual
drainage by as much as a factor of 1.6 (Table 4 and Figure 5),
which is similar to the overestimation of water storage changesuggested by the Texas site data. The impact of the upperboundary condition was evaluated by conducting UNSAT-H
simulations, with the upper boundary condition specified asin VS2DTI (PE = 0 on days with precipitation; simulation 15).
These UNSAT-H simulations generally reproduced theresults from the VS2DTI simulations, particularly in
WY98. A similar approach was used to evaluate the impactof the upper boundary condition used in HYDRUS-1D. The
UNSAT-H simulations that applied net precipitation or PE on
precipitation days (similar to HYDRUS-1D) generally reproduced the simulation results from HDYRUS-1D, particularly
during pre-WY98 (Table 4, simulation 16). These results
indicate that variations in the upper boundary condition can
generally explain the differences in simulation results among
the codes HYDRUS-1D, UNSAT-H, and VS2DTI andemphasize the impact of the approach used to simulate theupper boundary condition on the simulation results.
3.4. Implications for Water Balance Modeling[44] The results of this study have important implications
for modifying existing, or developing new, water balance
codes and for applying such codes to simulate water balances. Results from this study suggest that the Richards’
equation-based codes are more appropriate for simulating
the near-surface water balance than those based on storage
routing, such as HELP. The storage-routing approachassumes that gravity is the only driving force in water
movement, and the approach approximates upward flowusing a depth zone for evaporation. Although the storagerouting approach accurately simulated drainage in the latter 2

Figure 6. (a) Unsaturated hydraulic conductivities and (b)ratios of hydraulic conductivities calculated from van
Genuchten (VG) and Brooks and Corey (BC) water
retention functions and from Mualem (M) and Burdine(B) hydraulic conductivity functions. The various ratios areused to show differences in unsaturated hydraulic conductivity on the basis of different water retention functions(BC versus VG), different hydraulic conductivity functions(M versus B), or both.
years at the Idaho site, it greatly overestimated drainage
during the first few months when the site was irrigated andoverestimated water storage change at the Texas site. Previous studies also demonstrate that this approach generally
overestimates drainage [Berger et al., 1996; Khire et al.,
1997; Wilson et al., 1999]. The greater accuracy of the
Richards’ equation based models, the increased computational efficiency of these codes in recent years, recent
advances in computer technology, and the availability of
online hydraulic parameter data make the use of Richards’
equation based codes more feasible.
[45] This study demonstrated the difficulties of simulating
infiltration-excess runoff because of high intensity precipitation and irrigation inputs and uncertainties in hydraulicconductivities. The results suggest that it may be infeasible torepresent actual precipitation intensity and to measure
hydraulic conductivity accurately enough to reproduce measured runoff. An alternative approach may be to calibrate themodel to simulate runoff by varying the hydraulic conductivity of the shallow soils. Use of daily precipitation input hada large impact on partitioning of water between evaporation,water storage, and drainage because of the different
approaches used by the various codes for simulating theupper boundary condition. Many previous studies have useddaily precipitation input to predict near-surface water balance[Fayer and Gee, 1992; Rockhold et al., 1995; Kearns andHendrickx, 1998] and some studies have used monthly
precipitation input [Wilson et al., 1999]. While some codes(SoilCover, SWIM) simulate evaporation as a sink, allowing
evaporation and infiltration to occur at the same time, many
codes simulate these processes separately. Setting PE to zero

SCANLON ET AL.: INTERCODE COMPARISON

on rain days underestimates evaporation (VS2DTI), whereas
subtracting PE from precipitation and applying net precipitation or net PE overestimates evaporation (HYDRUS-1D)(Tables 3 and 4). Applying precipitation at a specified rate
and allowing evaporation to occur during the remainder ofthe day seemed to approximate hourly precipitation data best
as shown by results from the UNSAT-H simulations of the
Texas data (Table 3). Therefore researchers who use daily
precipitation as input should understand how the code distributes precipitation and PE and should assess the impact ofdaily data by comparing simulation results with those basedon hourly or shorter precipitation input, if possible. Previous
studies by Stothoff [1997] noted that use of daily precipitation
input underestimated infiltration by as much as a factor of 3,
relative to hourly precipitation input using the BREATH code
(evaporation simulated as a sink, no measured runoff)
because spreading precipitation input throughout the day
allowed more evaporation to occur. Schemes to disaggregate
daily precipitation to smaller time steps have been developed
for select locations [Egbuniwe, 1975; Econopouly et al.,
1990; Socolofsky et al., 2001], but parameters need to be
developed for other regions. Such disaggregation could help
simulation of water partitioning more accurately among thevarious water balance components.
[46] Although conducting sensitivity analysis to hydraulicparameters was not the primary objective of this study, someof the simulation results demonstrated sensitivity to hydraulicparameters. Model input parameters for the Idaho site werebased on calibration using the UNSAT-H code with van
Genuchten water retention functions and Mualem hydraulicconductivity functions; therefore codes that use similarfunctions were expected to calculate water balances that
better match measured values than the other codes. Resultsfrom simulations of the Idaho site in this study showed thatuse of the Campbell water retention function (Brooks andCorey with zero residual water content, simulation 12) andthe Burdine hydraulic conductivity function (SHAW, simulation 3) resulted in increased evaporation and reduced
drainage relative to the van Genuchten water retention andMualem hydraulic conductivity functions (UNSAT-H, simulation 1). The main difference between the codes resulted
from different water retention functions. The Brooks andCorey water retention function resulted in much higher
unsaturated hydraulic conductivity relative to the van Genuchten function (Figure 6), which translates into higher
evaporation and lower drainage (Table 4).
[47] Simulation results from this study showed muchgreater sensitivity to water retention and hydraulic conductivity functions than previously demonstrated. Previous
studies show that the impact of using different hydraulicconductivity functions (Mualem versus Burdine) in a sandy
soil resulted in = 10% variation in the simulated waterbalance [Fayer and Gee, 1992]. The Fayer and Gee study
also showed that reducing residual-water content from 0.05
m3/m3 to 0.00 decreased drainage by 27%. Andraski andJacobson [2000] showed that use of different water retentionfunctions (van Genuchten [1980] and Rossi and Nimmo[1994]) resulted in differences in simulated water potentials
in the shallowest nodes but had little impact on the simulated
water balance components.
[48] Previous sensitivity analyses by Fayer and Gee[1992] indicated that inclusion of isothermal vapor flow

59 - 13was significant (reduced simulated drainage by as much as65%). Vapor flow was not simulated using any of the codesin this study except SoilCover. Simulations conducted withand without isothermal vapor flow using UNSAT-H weresimilar (simulations 1 and 12, Table 3; 1, 17, Table 4),
indicating that isothermal vapor flow was not important for
simulating the water balance of these sites. Fayer et al.
[1992] indicated that hysteresis was critical in simulating
drainage in a capillary barrier engineered cover. Most codes
evaluated in this study do not simulate hysteresis, with the
exception of HYDRUS-1D and UNSAT-H; therefore theimpact of hysteresis could not be evaluated using all codes.
Additional simulations were conducted with HYDRUS-1Dand UNSAT-H to evaluate hysteresis. Water retention functions for the wetting curves were estimated by assuming thatthe van Genuchten a parameter for wetting is 2 that fordrying and allowing 10% for entrapped air during rewetting.
Results from these simulations indicate that hysteresis has anegligible impact on the simulation results for the Texas site
(simulations 9 and 13; Table 3). Hysteresis had little impact
on the HYDRUS-1D simulations of the Idaho site (simulation 10; Table 4); however, hysteresis resulted in higher
drainage, particularly in pre WY 98 when UNSAT-H wasused (simulation 18; Table 4). This result is consistent withthe results from Fayer et al. [1992].
[49] Use of the appropriate lower boundary condition isimportant in accurately simulating water balance of lysimeters. Flury et al. [1999] also point out the importance of thelower boundary condition on solute transport in wickless
lysimeters because the soil has to be saturated before drainage
can occur. Wickless lysimeters are most accurately represented by a seepage face; however, some codes do not include
a seepage-face option. Many previous modeling studies ofengineered covers monitored by wickless lysimeters used aunit gradient lower boundary condition [Fayer and Gee,
1992; Khire et al., 1997]. Results from this study indicatethat the unit-gradient lower boundary condition could greatly
overestimate drainage from a wickless lysimeter that does notcontain an overlying gravel layer (Table 4). Simulations of
the Idaho site demonstrate that models using Richards’
equation can approximate a seepage face by simulating athin layer of gravel with a unit-gradient lower boundary
condition. Because many studies try to accurately simulate
drainage from lysimeters on the order of millimeters, specifying the most appropriate boundary condition is critical. The
impact of the seepage face versus unit gradient lower boundary condition only applies to wickless lysimeters, and not
the natural system. Unit gradient lower boundary conditions
should work well for the natural system.
[50] The ultimate goal of most water balance modelingstudies is to predict drainage. The drainage term is generally
much smaller than many of the other terms in the water
budget in semiarid regions, particularly precipitation andevaporation. In contrast, drainage is generally much higher
in humid regions (e.g., 6 to 51 cm; 1985 – 1994, Coshochton, Ohio) and more easily simulated [Wilson et al., 1999].
In addition, the drainage term generally accumulates errors
in all other terms in the water budget. Therefore it isimportant to consider the uncertainty in simulating various
water balance components, such as runoff, evaporation, andwater storage change, when estimating uncertainty in drainage estimates.

59 - 14SCANLON ET AL.: INTERCODE COMPARISON[51] In addition to providing insight on important featuresof codes for simulating near-surface water balance, theintercode comparison also provided information on thecomputational efficiency of various codes. Most codes tested
in this study required less than 1 min to run a 1-year
simulation; however, longer times were required by SHAW
(10 min), SoilCover (1 hr), and VS2DTI (8 min) on a
personal computer (Pentium 4, 1.4 GHz CPU, 256 MB
RAM). A variety of factors impact the computational efficiency of codes. Internal tabulation of water retention functions can result in reduced computational times by as much asa factor of 3 because table lookup (HYDRUS-1D, SWIM) ismuch faster than calculating results from an analytical
function.
[52] Code comparisons should be conducted regularly sothat more can be learned about code performance andimportant factors in simulating the near-surface water balance can be assessed. Because source codes are generally notprovided and only executable versions of codes are generally
available, it is often difficult to determine how a code
simulates a particular process. The proliferation of codes
and the inability of user manuals to address all potential
applications increase the need for code comparison studies.
To facilitate future model testing and code-comparison studies, a database should be developed of water balance monitoring data for different climate and soil conditions.

4. Conclusions[53] This study demonstrates the variability in simulated
water balance components using a variety of different
codes (HELP, HYDRUS-1D, SHAW, SoilCover, SWIM,UNSAT-H, and VS2DTI) on the basis of field monitored
data from engineered covers at warm (Texas) and cold
desert (Idaho) sites and provides some indication of the
expected reliability of simulated water balances. Simulation results from most codes were similar and generally
reproduced measured water balance components at the
Texas and Idaho sites. Both sites consisted of unvegetatedloam soil.
[54] Simulation of infiltration-excess runoff was a problem for all codes, underscoring the difficulties of representing actual precipitation intensities and of measuring
hydraulic conductivity of surficial sediments (as shownby the data from Texas). Drainage is the most criticalparameter for evaluation of contaminant transport, engineered covers for waste containment, and groundwater
recharge. Drainage could be estimated to within ±64%by most codes. Outliers for the various simulations could be
attributed to the following factors: (1) the modelingapproach, i.e., water storage routing versus Richards’ equation, (2) the upper boundary condition during precipitationand time discretization of precipitation input, (3) water
retention functions (i.e., van Genuchten versus Brooks
and Corey), and (4) the lower boundary condition (i.e.,unit gradient versus seepage face).
[55] The water storage routing approach does not seemto adequately represent the flow system in semiarid
regions. By assuming that gravity is the only driving force
and ignoring matric-potential gradients that are often
upward in semiarid regions, downward flow is generally
overestimated and ultimately results in overestimation ofdrainage.

[56] The approach used to simulate the upper boundary
condition during precipitation is crucial when precipitation
is input on a daily or larger time step. Setting PE to zero on
rain days (VS2DTI) greatly underestimated evaporation andoverestimated drainage. Subtracting PE from precipitation
and applying net precipitation or net PE on a daily basis(HYDRUS-1D) had a much lesser impact on simulation
results. The best approach is to disaggregate daily precipitation and apply it at a specified rate, allowing PE to occur
throughout the rest of the day, as shown by the UNSAT-H
simulations.
[57] The impact of water retention functions was demonstrated at the Idaho site, where increased unsaturated
hydraulic conductivity based on the Brooks and Coreyfunctions relative to the van Genuchten functions resulted
in overestimation of evaporation and underestimation ofdrainage. In contrast, the input value of residual watercontent (0 for Campbell function versus > 0 for Brooks
and Corey) had little impact on simulation results.
[58] The most appropriate lower boundary condition forsimulating wickless lysimeters is a seepage face. Simulations
using HYDRUS-1D demonstrated that this boundary condition could be approximated by simulating a thin bottom
layer of gravel with a unit gradient boundary condition incodes that use Richards’ equation but do not include a
seepage face option. However, use of a unit-gradient lower
boundary condition alone greatly overestimated drainage.
This study demonstrates the usefulness of conducting intercode comparisons to evaluate the reliability of water balance
simulations and to determine important factors controlling
water balance simulation results.

[59] Acknowledgments. This project was funded by the U.S. Environmental Protection Agency under contract IAC (CR827598-01-0). Publication of this manuscript was authorized by the Director, Bureau of
Economic Geology, The University of Texas at Austin. The authors gratefully acknowledge the assistance of Kelley Keese in running many of thesimulations. Rick Healy (USGS) provided help with VS2DTI, and Mike
Fayer provided valuable insights on UNSAT-H simulations and waterbalance modeling. The authors greatly appreciate extremely helpful andthorough reviews provided by three anonymous reviewers.

References
Andraski, B. J., and E. A. Jacobson, Testing a full-range soil-water retention function in modeling water potential and temperature, Water Resour.
Res., 36, 3081 – 3090, 2000.
ASTM D2434, Standard test method for permeability of granular soils(constant head), pp. 191 – 195, Am. Soc. Testing Materials, Philadelphia,
Pa., 1994.
ASTM D5084, Standard test method for measurement of hydraulic conductivity of saturated porous materials using a flexible wall permeameter,
pp. 63 – 70, Am. Soc. Testing Materials, Philadelphia, Pa., 1990.
Belmans, C., J. G. Wesseling, and R. A. Feddes, Simulation model of thewater balance of a cropped soil: SWATRE, J. Hydrol., 63, 271 – 286, 1983.
Berger, K., S. Melchior, and G. Miehlich, Suitability of Hydrologic Evaluation of Landfill Performance (HELP) model of the US Environmental
Protection Agency for the simulation of the water balance of landfill
cover systems, Environ. Geol., 28, 181 – 189, 1996.
Brooks, R. H., and A. T. Corey, Properties of porous media affecting fluid
flow, J. Irrig. Drainage Div. Am. Soc. Civ. Eng., 72(IR2), 61 – 88, 1966.
Burdine, N. T., Relative permeability calculation size distribution data,
Trans. Am. Inst. Min. Metall. Pet. Eng., 198, 71 – 78, 1953.
Campbell, G. S., A simple method for determining unsaturated conductivity
from moisture retention data, Soil Sci., 117, 311 – 314, 1974.
Campbell, G. S., Soil Physics with BASIC: Transport Models for Soil-Plant
Systems, Elsevier Sci., New York, 1985.
Christie, M. A., and M. J. Blunt, Tenth SPE comparative solution project: A
comparison of upscaling techniques, SPE Reservoir Eng. Eval. Eng., 4,308 – 317, 2001.

SCANLON ET AL.: INTERCODE COMPARISON
Dawes, W. R., L. Zhang, T. J. Hatton, P. H. Reece, G. T. H. Beale, and I.
Packer, Application of a distributed parameter ecohydrological model(TOPOG_IRM) to a small cropping rotation catchment, J. Hydrol.,191, 67 – 89, 1997.
Doorenbos, J. and W. O. Pruitt, Guidelines for predicting crop waterrequirements, FAO Irrig. Pap. 24, 2nd ed., Food and Agric. Organ.,Rome, 1977.
Econopouly, T. W., D. R. Davis, and D. A. Woolhiser, Parameter transferability for a daily rainfall disaggregation model, J. Hydrol., 118, 209 –
228, 1990.
Egbuniwe, N., Generation of hourly rainfall, Water Resour. Bull., 11, 706 –
713, 1975.
Fayer, M. J., UNSAT-H version 3.0: Unsaturated soil water and heat flow
model, theory, user manual, and examples, Rep.13249, Pac. Northwest
Natl. Lab., Richland, Wash., 2000.
Fayer, M. J. and G. W. Gee, Predicted drainage at a semiarid site: Sensitivity to hydraulic property description and vapor flow, in Indirect Methods for Estimating the Hydraulic Properties of Unsaturated Soils, editedby M. T. van Genuchten, F. J. Leij, and L. J. Lund, pp. 690 – 700, Univ.
of Calif., Riverside, 1992.
Fayer, M. J., M. L. Rockhold, and M. D. Campbell, Hydrologic modelingof protective barriers: Comparison of field data and simulation results,
Soil Sci. Soc. Am. J., 56, 690 – 700, 1992.
Flerchinger, G. N., and K. E. Saxton, Simultaneous Heat and water modelof a freezing snow-residue-soil system, 1, Theory and development,
Trans. ASAE, 32, 565 – 571, 1989.
Flerchinger, G. N., C. L. Hanson, and J. R. Wight, Modeling evapotranspiration and surface energy budgets across a watershed, Water Resour.
Res., 32, 2539 – 2548, 1996.
Flint, A. E., L. E. Flint, J. A. Hevesi, F. D’Agnese, and C. Faunt, Estimationof regional recharge and travel time through the unsaturated zone in arid
climates, in Dynamics of Fluids in Fractured Rock, Geophys. Monogr.
Ser., vol. 122, edited by B. Faybishenko, P. A. Witherspoon, and S. M.
Benson, pp. 115 – 128, AGU, Washington, D. C., 2000.
Flury, M., M. V. Yates, and W. A. Jury, Numerical analysis of the effect of
the lower boundary condition on solute transport in lysimeters, Soil Sci.
Soc. Am. J., 63, 1493 – 1499, 1999.
Fredlund, D. G., and A. Xing, Equations for the soil-water characteristic
curve, Can. Geotech. J., 31, 521 – 532, 1994.
Fredlund, D. G., A. Xing, and S. Huang, Predicting the permeability function for unsaturated soils using the soil-water characteristic curve, Can.
Geotech. J., 31, 533 – 546, 1994.
GEO2000, SoilCover user’s manual version 4.01, Unsaturated Soils Group,
Univ. of Saskatchewan, Saskatoon, Saskatchewan, Canada, 1997.
Hatton, T., The basics of recharge and discharge, part 4, Catchment scale
recharge modeling, 19 pp., Commonwealth Sci. Indust. Res. Organ.,
Collingwood, Victoria, Australia, 1998.
Healy, R. W., Simulation of solute transport in variably saturated porous
media with supplemental information on modifications to the U.S. Geological Survey’s computer program VS2D, U.S. Geol. Surv. Water Resour. Invest. Rep., 90-4025, 125 pp., 1990.
Hillel, D., and C.H.M. van Bavel, Simulation of profile water storage as
related to soil hydraulic properties, Soil Sci. Soc. Am. J., 40, 807 – 815,
1976.
Kearns, A. K. and J. M. H. Hendrickx, Temporal variability of diffuse
groundwater recharge in New Mexico, Tech. Completion Rep. 309, 43
pp., N. M. Water Resour. Res. Inst., Las Cruces, 1998.
Khire, M. V., C. H. Benson, and P. J. Bosscher, Water balance modeling ofearthen final covers, J. Geotech. Geoenviron. Eng., 123, 744 – 754, 1997.
Khire, M. V., C. H. Benson, and P. J. Bosscher, Capillary barriers: Design
variables and water balance, J. Geotech. Geoenviron. Eng., 126, 695 –
708, 2000.
Kim, C. P., J. N. M. Stricker, and P. J. J. F. Torfs, An analytical framework
for the water budget of the unsaturated zone, Water Resour. Res., 32,
3475 – 3484, 1996.
Lappala, E. G., R. W. Healy, and E. P. Weeks, Documentation of computer
program VS2D to solve the equations of fluid flow in variably saturated
porous media, U.S. Geol. Surv. Water Resour. Invest. Rep., 83-4099, 184pp., 1987.
Larssen, A., K. Pers, K. Skegius, and B. Dverstorp, International INTRAVAL project to study validation of geospheric transport models for performance assessment of nuclear waste disposal, phase 2 summary report,
118 pp., Swed. Nucl. Power Inspectorate, Stockholm, 1995.
Leij, F. J., W. J. Alves, M. T. van Genuchten, and J. R. Williams, The
UNSODA unsaturated soil hydraulic database, user’s manual, version
1.0, EPA/600/R-96/095, Natl. Risk Manage. Lab., Off. of Res. andDev., Cincinnati, Ohio, 1996.

59 - 15

Lewis, M. F., and G. R. Walker, Assessing the potential for episodic recharge in south-western Australia using rainfall data, Hydrogeol. J., 10,229 – 237, 2002.
Magnuson, S. O., A simulation study of moisture movement in proposed
barriers for the subsurface disposal area, EGG-WM-10974, 64 pp., Idaho
Natl. Eng. Lab., Idaho Falls, 1993.
Meyer, P. D., Application of an infiltration evaluation methodology to ahypothetical low-level waste disposal facility, NUREG/CR-6114, PNL
8842, 28 pp., Nucl. Regulatory Agency, Rockville, Md., 1993.
Milly, P. C. D., Climate, soil water storage, and the average annual waterbalance, Water Resour. Res., 30, 2143 – 2156, 1994.
Milly, P. C. D., A minimalist probabilistic description of root zone soil
water, Water Resour. Res., 37, 457 – 463, 2001.
Milly, P. C. D., and K. A. Dunne, Sensitivity of the global water cycle to thewater-holding capacity of land, J. Clim., 7, 506 – 526, 1994.
Mockus, V., Hydraulic effects of land use and treatment, in Estimation of
Direct Runoff From Snowmelt, SCS National Engineering Handbook,sect. 4, Hydrology, U.S. Dep. of Agriculture, Soil Conserv. Serv., Washington, D. C., 1972.
Monteith, J. L. and M. H. Unsworth, Principles of Environmental Physics,
Edward Arnold, London, 1990.
Mualem, Y., A new model for predicting the hydraulic conductivity ofunsaturated porous media, Water Resour. Res., 12, 513 – 521, 1976.
Nichols, W. E., Comparative simulations of a two-layer landfill barrier
using the HELP version 2.0 and UNSAT-H version 2.0 computer codes,
Doc. PNL-7583, Pac. Northwest Lab., Richland, Wash., 1991.
Ogan, B. D., G. V. Wilson, W. H. Albright, G. W. Gee, M. J. Fayer, and
S. Rock, Sensitivity analysis and validation of numerical models used in
the design of alternative landfill covers, paper presented at Annual Meeting, Soil Sci. Soc. of Am., Salt Lake City, Utah, 31 Oct. to 4 Nov. 1999.
Penman, H. L., Natural evaporation from open water, bare soil and grass,
Proc. R. Soc. London, Ser. A, 193, 120 – 145, 1948.
Pierce, L. L., J. Walker, T. I. Dowling, T. McVicar, T. J. Hatton, S. W.
Running, and J. C. Coughlan, Hydroecological changes in the MurrayDarling Basin, part 3, A simulation of regional hydrological changes, J.
Applied Ecol., 30, 283 – 294, 1993.
Porro, I., Hydrologic behavior of two engineered barriers following extreme
wetting, J. Environ. Qual., 30, 1 – 13, 2001.
Porro, I., and P. Martian, Hydrologic modeling of engineered barriers following a wetting test, Agronomy Abstracts, Am. Soc. of Agron., Madison, Wisc., 1999.
Richardson, C. W., Stochastic simulation of daily precipitation, temperature, solar radiation, Water Resour. Res., 17, 182 – 190, 1981.
Richardson, C. W., Data requirements for estimation of weather generation
parameters, Trans. Am. Soc. Agric. Eng., 43, 877 – 882, 2000.
Ritchie, J. T., A model for predicting evaporation from a row crop withincomplete cover, Water Resour. Res., 8, 1204 – 1213, 1972.
Rockhold, M. L., M. J. Fayer, and P. R. Heller, Physical and hydraulicproperties of sediments and engineered materials associated with grouted
double-shell tank waste disposal at Hanford, Report PNL-8813, Pac.
Northwest Lab., Richland, Wash., 1993.
Rockhold, M. L., M. J. Fayer, C. T. Kincaid, and G. W. Gee, Estimation ofnatural ground water recharge for the performance assessment of a lowlevel waste disposal facility at the Hanford site, PNL-10508, Battelle Pac.
Northwest Natl. Lab., 1995.
Rodriguez-Iturbe, I., A. Porporato, F. Laio, and R. Rodolfi, Plants in watercontrolled ecosystems: Active role in hydrological processes and response to water stress, 1, Scope and general outline, Adv. Water Resour.,24, 695 – 706, 2001.
Rossi, C., and J. R. Nimmo, Modeling of soil water retention from saturation to oven dryness, Water Resour. Res., 30, 701 – 708, 1994.
Scanlon, B. R., and P. C. D. Milly, Water and heat fluxes in desert soils, 2,
Numerical simulations, Water Resour. Res., 30, 721 – 733, 1994.
Schaap, M. G., and F. J. Leij, Database-related accuracy and uncertainty ofpedotransfer functions, Soil Sci., 163, 765 – 779, 1998.
Schaap, M. G., F. J. Leij, and M. T. van Genuchten, Neural network
analysis for hierarchical prediction of soil hydraulic properties, Soil
Sci. Soc. Am. J., 62, 847 – 855, 1998.
Schroeder, P. R., C. M. Lloyd, P. A. Zappi, and N. M. Aziz, The Hydrologic
Evaluation of Landfill Performance (HELP) Model, user’s guide forversion 3, EPA/600/168a, U.S. Environ. Prot. Agency Risk ReductionEng. Lab., Cincinnati, Ohio, 1994.
Simmons, C. S., and P. D. Meyer, A simplified model for the transient water
budget of a shallow unsaturated zone, Water Resour. Res., 36, 2835 –
2844, 2000.
Simunek, J., M. Sejna, and M. T. van Genuchten, The Hydrus-1D software
package for simulating the one-dimensional movement of water, heat,59 - 16SCANLON ET AL.: INTERCODE COMPARISON

and multiple solutes in variably-saturated media, user’s manual, version
2.0, U.S. Salinity Lab., Agric. Res. Serv., Riverside, Calif., 1998.
Socolofsky, S. E., E. Adams, and D. Entekhabi, Disaggregation of daily
rainfall for continuous watershed modeling, J. Hydrol. Eng., 6, 300 – 309,
2001.
Stothoff, S. A., Sensitivity of long-term bare soil infiltration simulations tohydraulic properties in an arid environment, Water Resour. Res., 33,
547 – 558, 1997.
Topp, G. C., J. L. Davis, and A. P. Annan, Electromagnetic determination ofsoil water content: Measurement in coaxial transmission lines, Water
Resour. Res., 16, 574 – 582, 1980.
U.S. Department of Agriculture, National Engineering Handbook, sect. 4,
Hydrology, U.S. Gov. Print. Off., Washington, D. C., 1985.
van Genuchten, M. T., A closed-form equation for predicting the hydraulicconductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44, 892 – 898, 1980.
Verburg, K., P. J. Ross, and K. L. Bristow, SWIMv2.1 user manual, Div. Rep.
130, Commonwealth Sci. Indust. Res. Organ., Canberra, Australia, 1996.
Wilson, G. V., W. H. Albright, G. W. Gee, M. J. Fayer, and B. D. Ogan,
Alternative cover assessment project, Phase 1 report, final technical report prepared by Desert Res. Inst., Nev., and Pac. Northwest Natl. Lab.,
Wash., Contract 44-0000-2038, 202 pp., U.S. Environ. Prot. Agency,
Washington, D. C., 1999.

Wilson, G. W., Soil evaporative fluxes for geotechnical engineering problems, Ph.D. thesis, Univ. of Saskatchewan, Saskatoon, Saskatchewan,
Canada, 1990.
Wilson, G. W., D. G. Fredlund, and S. L. Barbour, Coupled soil-atmosphere modeling for soil evaporation, Can. Geotech. J., 31, 151 – 161,
1994.




M. Christman, GeoSyntec Consultants, 106 E. 6th Street, Suite 800,
Austin, TX 78701, USA. (mchristman@geosyntec.com)
G. N. Flerchinger, U.S. Department of Agriculture, REE-ARS-PWAWMR, Boise, ID 83712, USA. (gflerchi@nwrc.ars.usda.gov)
I. Porro, Idaho National Engineering and Environmental Laboratory,
P.O. Box 1625, Idaho Falls, ID 83415, USA. (ixp@inel.gov)
R. C. Reedy and B. R. Scanlon, Bureau of Economic Geology,
University of Texas at Austin, 10100 Burnet Road, Austin, TX 78758,
USA. (robert.reedy@beg.utexas.edu; bridget.scanlon@beg.utexas.edu)
J. Simunek, George E. Brown Jr. Salinity Laboratory, 450 W. Big
Springs Road, Riverside, CA 92507, USA. ( jsimunek@ussl.ars.usda.gov)