Click
Here

WATER RESOURCES RESEARCH, VOL. 45, W04412, doi:10.1029/2008WR006903, 2009

for

Full
ArticleUsing pedotransfer functions in vadose zone models
for estimating groundwater recharge
in semiarid regions
Tiejun Wang,1,2 Vitaly A. Zlotnik,1 Jirka Šimunek,3 and Marcel G. Schaap4
Received 6 February 2008; revised 19 January 2009; accepted 4 February 2009; published 14 April 2009.

[1] Process-based vadose zone models are becoming common tools for evaluating spatial

distributions of groundwater recharge (GR), but their applications are restricted bycomplicated parameterizations, especially because of the need for highly nonlinear andspatially variable soil hydraulic characteristics (SHCs). In an attempt to address the
scarcity of field SHC data, pedotransfer functions (PTF) were introduced in earlier
attempts to estimate SHCs. However, the accuracy of this method is rarely questioned inspite of significant uncertainties of PTF-estimated SHCs. In this study, we investigated theapplicability of coupling vadose zone models and PTFs for evaluating GR in sand andloamy sand soils in a semiarid region and also their sensitivity to lower boundary
conditions. First, a data set containing measured SHCs was used in the simulations. A
second data set contained correlated SHCs drawn from the covariance matrix of the firstdata set. The third SHC data set used was derived from a widely used PTF. Althoughstandard deviations for individual parameters were known for this PTF, no covariancematrix was available. Hence, we assumed that the parameters of this PTF were
uncorrelated, thereby potentially overestimating the volume of the parameter space.
Results were summarized using histograms of GR for various sets of input parameters.
Under the unit gradient flow lower boundary condition, the distributions of GR for sandand loamy sand significantly overlap. Values of GR based on mean SHCs (or GR*)generally lie off the mode of the GR distribution. This indicates that the routinely usedmethod of taking GR* as a regional representation may not be viable. More importantly,
the computed GR largely depends in a nonlinear fashion on the shape factor n in thevan Genuchten model. Under the same meteorological conditions, a coarser soil with alarger n generally produces a higher GR. Therefore, the uncertainty in computed GR
is largely determined by the uncertainty in estimated n by PTFs (e.g., mean and standard
deviation). Under the constant head lower boundary condition, upward soil moisture
flux may exist from the lower boundary. Especially for regions with shallow water tables
where upward flux exists, choosing an appropriate lower boundary condition is more
important than selecting SHC values for calculating GR. The results show that the
distribution of GR is less scattered and GR is more intense if the constant head lower
boundary is located at deeper depths.
Citation: Wang, T., V. A. Zlotnik, J. Šimunek, and M. G. Schaap (2009), Using pedotransfer functions in vadose zone models for
estimating groundwater recharge in semiarid regions, Water Resour. Res., 45, W04412, doi:10.1029/2008WR006903.

1. Introduction[2] Analysis of groundwater sustainability in semiarid
areas relies on the knowledge of spatial distributions ofgroundwater recharge (GR) [Scanlon et al., 2006]. Althoughvarious approaches have been developed on the basis of1

Department of Geosciences, University of Nebraska at Lincoln,
Lincoln, Nebraska, USA.
2
Now at School of Natural Resources, University of Nebraska at
Lincoln, Lincoln, Nebraska, USA.
3
Department of Environmental Sciences, University of California,
Riverside, California, USA.
4
College of Agriculture and Life Sciences, University of Arizona,
Tucson, Arizona, USA.
Copyright 2009 by the American Geophysical Union.
0043-1397/09/2008WR006903$09.00

different physical, chemical, mathematical, and isotopic
techniques [Lerner et al., 1990; Hendrickx and Walker,
1997; Scanlon et al., 2002], estimation of GR still remains
one of the most difficult tasks in hydrology [National
Research Council, 2004], as GR may substantially varyacross landscapes because of spatial variations in soil texture,
vegetation, climate forcing, and hydrologic boundary conditions. Compared to methods that yield an areal-averaged
GR (e.g., catchment water balance method), hydrodynamic
process-based vadose zone models are becoming common
tools for evaluating GR and its spatial distribution [Keese etal., 2005; Small, 2005; Nolan et al., 2007]. However, the useof vadose zone models is usually restricted by complicatedparameterizations, especially because of lack of in situ soil
hydraulic characteristics (SHCs). Instead of measuring in situ

W04412

1 of 12

W04412WANG ET AL.: USING PEDOTRANSFER FUNCTIONS IN VADOSE ZONE MODELS

W04412

Figure 1. Location of the Barta Brothers Ranch Site (BBRS) in the Nebraska Sand Hills (contoured
area indicates the Sand Hills).
SHCs, estimation techniques have been developed such as
pedotransfer functions (PTF) [Wösten et al., 2001]. On thebasis of regression analysis of existing soil databases, PTFs
convert easily measurable soil properties (e.g., soil texture,
porosity, and bulk density) to SHCs. With the help of PTFs,
vadose zone models have been used for evaluating spatial
distributions of GR at large scales [Nolan et al., 2003, 2007;
Keese et al., 2005; Small, 2005; Anuraga et al., 2006].
[3] Using class-averaged SHCs for different soil textures[Schaap et al., 2001] and a vadose zone model, Small[2005] analyzed impacts of soil texture and climatic forcing
on GR in the southwestern United States. The study showedthat GR can occur from infrequent storms in coarser soils
when the ratio of precipitation (P) to potential evapotranspiration (ETp) exceeds 0.4, but cannot occur in finer soils
when P/ETp < 0.8. Keese et al. [2005] also used mean SHCs
estimated by the Rosetta program [Schaap et al., 2001] toevaluate GR in Texas and assessed impacts of climate,vegetation, and soil texture on GR. The study showed that
GR is higher in bare and homogeneous sand and can be
reduced by factors of 2 – 11 for soils whose textures varywith depth. On the basis of a similar method [Keese et al.,2005; Small, 2005], Anuraga et al. [2006] found that GR in
south India is more affected by soil texture than by land use.
[4] Previous studies that coupled vadose zone models and
PTFs have shown the significant impact of SHCs on evaluating GR under varying vegetation and climatic conditions.
However, the use of mean SHCs for evaluating GR, which isthen routinely taken as a regional representation, might bequestionable because of significant uncertainties associatedwith PTF-estimated SHCs [Schaap and Leij, 1998; Schaap
et al., 2001]. Nolan et al. [2007] estimated GR in the eastern
United States using both the unit gradient Darcian methodwith mean SHCs taken from PTFs and a chloride massbalance method. The study showed that the GR estimated by
the Darcian method exhibits higher spatial variability because of the uncertainty in estimated SHCs [Nolan et al.,2007]. Faust et al. [2006] also pointed out that GR estimatedat basin scales depends on the choice of PTFs. It is thus
imperative to assess the impact of the uncertainty of PTF-

estimated SHCs on GR, especially when using only mean
SHCs for GR evaluations [Keese et al., 2005; Small, 2005].
[5] In most studies, the unit gradient flow condition is
generally assumed at the lower boundary of the vadose zone
[Nolan et al., 2003, 2007; Keese et al., 2005; Small, 2005].
This boundary implies that water that passes downward
across the boundary does not return to the simulated
domain. However, such an assumption may not hold inregions of shallow water tables, as studies [Gillham, 1984;
Wu et al., 1996; Batelaan et al., 2003] have shown
significant impacts of water table on GR through processes
such as capillary rise or groundwater evapotranspiration.
[6] This study investigates the impact of the uncertainty of
PTF-estimated SHCs on computing GR by varying soil texture
and lower boundary condition. To be consistent with previous
studies, groundwater recharge (GR) is defined here as theamount of water leaving the simulated domain (i.e., lower
boundary). A commonly used vadose zone model, Hydrus-1D[Simunek et al., 2005] was chosen because its accuracy has
been verified by analytical techniques [Zlotnik et al., 2007].
Soil types of sand and loamy sand were selected. For each soil
type, three data sets with the van Genuchten parameters wereused. The first data set contained measured SHCs and isreferred to as the measured data set. The second data set
contained correlated SHCs that were drawn from the covariance matrix of the first data set. The third data set was derived
from a widely used PTF [e.g., Schaap et al., 2001]. Althoughstandard deviations for individual parameters were known forthis PTF, no covariance matrix was available, and therefore theparameters of the generated SHCs were assumed to be
uncorrelated in this study. The hydrometeorological andphenomenological data needed for the simulations were measured at the University of Nebraska’s Barta Brothers Ranch
Site (BBRS) in the eastern Nebraska Sand Hills.

2. Materials and Methods
2.1. Study Area
[7] This study was performed as a part of the Grassland
Destabilization Experiment (GDEX) at the Barta Brothers

2 of 12

W04412WANG ET AL.: USING PEDOTRANSFER FUNCTIONS IN VADOSE ZONE MODELS

Figure 2. Daily precipitation (P), potential evaporation(Ep), and potential transpiration (Tp) at the GDEX site forthe year 2004.
Ranch Site (BBRS) (Figure 1), which investigates theecological and geomorphic stability of the Nebraska Sand
Hills from an interdisciplinary perspective [Wang et al.,2008]. The 58,000 km2 Nebraska Sand Hills is the largest
native grassland-stabilized sand dune field in the Western
Hemisphere [Loope and Swinehart, 2000]. The high infiltration capacity of sandy soils makes this area an important
groundwater recharge source for the Ogallala aquifer. The
mean annual temperature at the site is 8.1°C and the meanannual P is 57.6 cm; the mean annual ETp is 136.0 cm[Szilagyi et al., 2005] and the mean annual humidity index(e.g., P/ETp) is 0.42. About 90% of the landscape at the
BBRS is composed of upland dunes and dry interdunal areas
covered by native grasslands, while the remaining 10%
consists of wet meadows. Ten circular plots, each 120 m indiameter, were constructed at the site, five of which were
equipped with meteorological instruments for evaluatingwater and energy balances. Surface soils (top 10 cm) are
sandy on ridges (average 94.4% sand) and only slightly lesssandy in swales (average 91.2% sand). Beneath the 10 cm
depth, contents of sand range from 95% to 97% regardless of
topographic positions [Wang et al., 2008].
2.2. Flow Model Description[8] The 1-D Richards equation, describing water flow insaturated and unsaturated porous media, can be formulated as

 

@q
@
@h
¼
KðhÞ
 KðhÞ  S
@t @x
@x

ð1Þ

where q is volumetric water content; t is time; x is a spatial
coordinate; h is water pressure head; S is root water uptake;and K is unsaturated hydraulic conductivity. A commonlyused model, Hydrus-1D [Simunek et al., 2005] was utilized
for solving equation (1) and computing GR.
[9] The Mualem – van Genuchten model [van Genuchten,
1980; Mualem, 1976], which analytically describes therelations among q, h, and K was used:qðhÞ ¼ qr þ ðqs  qr Þ=½ð1 þ jahjn Þm ; h < 0qðhÞ ¼ qs ; hKðhÞ ¼ KS

ð2Þ

0
Sel

½1  ð1  Se1=m Þm 2

ð3Þ

W04412where qr is residual soil moisture content; qs is saturated soil
moisture content; a, n, and l are shape factors: a is inversely
related to the air entry pressure, n is a measure of pore size
distribution, and lumped parameter l accounts for pore
tortuosity and connectivity; m = 1  1/n; KS is saturated
hydraulic conductivity; and Se = (q  qr)/(qs  qr). The
details of the three data sets used in this study are given in
Section 3.
[10] The atmospheric upper boundary condition for vadosezone is standard [Neuman et al., 1974]. This boundary
condition depends on climatic conditions (i.e., P and
potential evaporation Ep) and can switch from a prescribed
flux to a prescribed head when limiting pressure heads are
exceeded. Surface runoff occurs when P exceeds soil
infiltration capacity or soil becomes saturated. Two typesof lower boundary conditions were included: unit hydraulicgradient flow condition and constant head condition. For theunit gradient flow boundary, the soil profile depth was 5 m
with 201 nodes between the surface and bottom. For theconstant head boundary, the soil profile depth was 10 mwith 500 nodes between the surface and bottom. Constantheads were set to be 0, 5, and 7 m at the lower boundary,which correspond to zero constant heads at depths of 10, 5,
and 3 m below the surface, respectively. To minimize theeffect of initial conditions, the simulations were repeated atleast 10 times until the soil moisture profiles were equilibrated with climatic forcings [Keese et al., 2005; Small,2005].
[11] The root water uptake was simulated according to themethod of Feddes et al. [1978], which distributes potential
transpiration (Tp) over the root zone on the basis of root
density distribution and reduces Tp to actual transpiration on
the basis of soil matric potential. The root parameters on
pasture were used to determine the Feddes function[Wesseling, 1991]. Potential evapotranspiration ETp was
calculated using the Penman-Monteith equation [Allen etal., 1998] and then partitioned into Tp and Ep on the basis ofBeer’s law [Ritchie, 1972]:
Ep ¼ ETpek LAITp ¼ ETp  Ep

ð4Þ

ð5Þ

where k is an extinction coefficient set to be 0.5 and LAI isleaf area index.
[12] Hydrometeorological data for the year 2004 at the
GDEX site were utilized, including P, air temperature,humidity, solar radiation, soil heat flux, and wind speed.
P was 49.3 cm/a or 85.6% of the annual mean and ETp was
160.6 cm/a or 118% of the annual mean. The 2004 humidity
index was 0.31, which was lower than the annual mean of
0.42 in this area. The daily P, Tp, and Ep are plotted in
Figure 2. Note that low Tp rate is due to the low LAI in the
study area. In 3 m deep soil cores analyzed for root biomass,
60 –70% of the total root mass occurred in the top 20 cm
and 85– 90% occurred in the top 50 cm. On the basis of the
field observations [Wang et al., 2008], an exponential
distribution was fitted to describe the root density distribution, which is consistent with recent literature reviews on theoccurrence of root biomass in native grasslands [Jackson

3 of 12WANG ET AL.: USING PEDOTRANSFER FUNCTIONS IN VADOSE ZONE MODELS

W04412

W04412

Table 1. Pearson Correlation and Covariance Matrices for the Measured Data Seta
Sand

qr
qs
log10 alog10 nlog10 KSl

Loamy Sandqrqs

log10 alog10 nlog10 KS

l

qr

qs

log10 alog10 nlog10 KS

l

0.000
0.000
0.002
0.002
0.0050.238
0.005
0.0003
0.006
0.005

0.037
0.338
0.007
0.144
0.056

0.215
0.0249
0.105
0.081
0.077

0.064
0.173
0.652
0.506
0.450

0.097
0.068
0.121
0.230
0.409
-

0.001
0.007
0.005
0.013
0.004

0.299
0.010
0.0020.028
0.068

0.5106
0.434
0.0420.034
0.1340.761
0.1664
0.69600.018
0.027

0.533
0.690
0.159
0.175
0.103

0.030
0.3300.123
0.0520.056-

aNumbers in bold are the Pearson correlation.

et al., 1996; Schenk and Jackson, 2002]. LAI was measured at the site for partitioning ETp.

3. Data Sets on Soil Hydraulic Characteristics[13] In this study, three SHC data sets for sand and loamy
sand were used. The first data set included measured SHCs
and is referred to as the measured data set. In comparison,
the other two were generated, with one containing correlated
SHCs randomly generated from the measured SHC data set
and the other one containing uncorrelated SHCs derived
from the database of Schaap et al. [2001].
[14] The first data set with measured SHCs was derived
from the UNSODA database (see details in work by Nemes
et al. [2001]) and was previously used by Schaap and Leij
[2000] and Schaap and van Genuchten [2005]. From thisdatabase, we derived 51 samples for sand and 19 samples
for loamy sand with all samples having a complete set ofmeasured hydraulic parameters (e.g., qr, qs, a, n, KS, and l)that are needed to compute GR.
[15] The second data set with correlated SHCs wasderived from the measured data set by computing the classaverage values for qr , qs, log10 a, log10 n, log10 KS, and l forsand and loamy sand. A covariance matrix for theseparameters was computed for each of the classes (Table 1),
which allows us to account for parameter uncertainty as
well as correlations among the parameters. Under theassumption that the means and covariance matrix for qr ,qs, log10 a, log10 n, log10 KS, and l are an adequate statisticaldescription of the parameter distribution for underlying
classes of sand and loamy sand, new samples can be drawn
with a Monte Carlo procedure to obtain a more complete
distribution of parameter combinations. To this end wefollowed the procedure by Carsel and Parrish [1988]. We
generated 200 random samples with correlated parameters
by multiplying the upper triangular (Choleski) decomposition [Press et al., 1988] of the covariance matrix with
vectors containing six N(0,1) randomly distributed numbers. The result was then added to the mean parameters ofeach class. Means and standard deviations of the generated
parameters as well as a visual inspection of all bivariate
distributions showed that the generated distributions overlapped the distributions of the measured data set. However,this exercise led to some l with ‘‘very’’ negative values,
which can lead to unreasonable K() relationships that may
show increasing conductivities with decreasing water contents [Durner et al., 1999; Schaap and Leij, 2000]. Rather
than constraining l to yield monotonic K() relations, we set
l to be 0.5 [Mualem, 1976] for simplicity because thenumerical experiment showed less sensitivity to l in our
analyses.
[16] With the above Monte Carlo procedure, the notion of
PTF is extended from the traditional class mean parameters
to a procedure that also includes parameter uncertainty andcorrelations with the objective to yield (correlated) distributions of parameters as a prediction [Carsel and Parrish,
1988]. To compare with a more traditional approach, classmeans and standard deviations (s) published by Schaap etal. [2001] (Table 2) were also used in this study. Anotherreason for picking the database of Schaap et al. [2001] isbecause it has been widely used in various modelingpurposes related to vadose zones, such as soil moisture
movement and groundwater recharge [Keese et al., 2005;
Small, 2005; Miller et al., 2007], contaminant transport[Dann et al., 2006], and root water uptake [Demirkanli etal., 2008; Segal et al., 2008], etc. Given the enormousvolume of data used, the database of Schaap et al. [2001]
gives reasonable ranges of variability in SHC parameters.
No complete covariance matrix exists for this PTF and forthis reason we assume that the parameters are uncorrelated.
GR was then calculated on the basis of combinations of(mean  s), mean, and (mean + s) for five selectedparameters (i.e., qs, log10 KS, log10 a, log10 n, and l), whichleads to total 35 = 243 simulations. In order to reduce thenumber of simulations, qr = 0.01 was set as it is unimportantin the GR analysis (see Figures 5a and 6a). Although someof the combinations might not occur in natural conditions,the shape of GR distributions largely depends on meanand s of the factor n instead of how the SHC data set isformulated as shown in the following section, which does
not affect our conclusions.
[17] Table 2 shows mean values and standard deviationsof the SHCs for each of the three data sets. Differences are
small among those three data sets, especially for themeasured and the correlated data sets, which confirms thatthe distributions of measured and generated parameters are
similar. The difference between the uncorrelated data set
and the other two data sets is relatively higher because thedifferent data sources they are based on [see Schaap andLeij, 1998; Schaap et al., 2001]. The largest differences
are found for log10 Ks of sand, and for log10 n and log10 Ks ofloamy sand.
[18] It should be noted that the convergence of simulations with higher n values is slow (e.g., n > 3). This isbecause the water content decreases dramatically within arelatively narrow pressure head interval from 0 cm to about
500 cm, and then stays close to the residual water content
for lower pressure heads. Consequently, the hydraulic

4 of 12WANG ET AL.: USING PEDOTRANSFER FUNCTIONS IN VADOSE ZONE MODELS

W04412

W04412

Table 2. Average Values of Soil Hydraulic Parameters for the Three Data Setsa
Sand
N

qr

qs

log10 a log10 n log10 KS

0.048
0.379 1.512 0.464(0.036) (0.046) (0.308) (0.223)
Correlated data set
200 0.056
0.383 1.517 0.498(0.032) (0.042) (0.305) (0.220)
Uncorrelated data set 308 0.053
0.375 1.450 0.500(0.029) (0.055) (0.250) (0.180)
Measured data set

51

Loamy Sand

2.544
(0.732)2.599(0.719)2.810(0.590)lN0.149 19(1.532)
0.174 200(1.42)
0.903 201(0.486)qr

qs

log10 a log10 n log10 KS

0.061
0.418 1.493 0.308(0.041) (0.069) (0.364) (0.174)
0.065
0.417 1.518 0.309(0.036) (0.070) (0.340) (0.156)
0.049
0.390 1.460 0.240(0.042) (0.070) (0.470) (0.160)

2.275(0.619)2.405(0.598)2.020(0.640)l
0.680(3.146)0.216(2.838)
0.874(0.594)aStandard deviations are given in parenthesis; N is the number of samples in each data set.

capacity has a narrow peak within this pressure headinterval and is more or less zero below it. The hydraulicconductivity also decreases dramatically within this interval
and is almost zero below it. The combination of these threesoil hydraulic functions often causes problems with convergence of the numerical solution for a dry condition where alarge change in pressure heads is accompanied by a minimalchange in water contents. Certain numerical measures (e.g.,
reduce the maximum pressure head at the surface) weretaken for overcoming this problem.

4. Results and Discussions[19] These simulations yielded important conclusions on
using process-based vadose models for evaluating GR. In
the following sections, histograms are systematically used

to summarize numerous simulations of GR and to illustrateeffects of various parameters and processes on GR.
4.1. Groundwater Recharge Under Unit Gradient Flow
Lower Boundary Condition
4.1.1. Impact of Soil Texture on Groundwater
Recharge[20] Under the unit gradient flow condition, the distributions of computed GR and ETa from the measured SHC data
set are plotted in Figure 3. Because of the high infiltration
capacity of sand and loamy sand and the use of a spin-up
method (e.g., repetition of the climatic data), surface runoffand change in soil moisture storage are usually negligible in
the water balance (e.g., <1%), and thus are not analyzed
here. Also note that three simulations for both sand andloamy sand are removed from Figure 3 because of low

Figure 3. Distributions of groundwater recharge (GR) for (a) sand and (b) loamy sand and actual
evapotranspiration (ETa) for (c) sand and (d) loamy sand based on the measured data set under the unit
gradient flow boundary condition. The asterisk represents the results based on mean SHCs. Bin size is 5 cm/a.
5 of 12

W04412WANG ET AL.: USING PEDOTRANSFER FUNCTIONS IN VADOSE ZONE MODELS

W04412

Figure 4. Distributions of groundwater recharge (GR) under the unit gradient flow lower boundary
condition: (a) sand and (b) loamy sand from the correlated data set and (c) sand and (d) loamy sand fromthe uncorrelated data set. The asterisk represents the results based on mean SHCs. Bin size is 5 cm/a.
negative l values, which led to errors of mass balanceexceeding 10%, compared to the rest with errors less than
<0.5%.
[21] GR for sand (GRS) ranges from 0.06 to 32.8 cm/a(e.g., <0.01% and 66.5% of the total P, respectively) and
ETa for sand (ETaS) ranges from 16.4 to 49.3 cm/a (e.g.,
33.3% and >99.9%, respectively). GR for loamy sand
(GRLS) ranges from 0.02 to 25.7 cm/a (e.g., <0.01% and52.1%, respectively) and ETa for loamy sand (ETaLS) ranges
from 28.4 to 48.9 cm/a (e.g., 57.6% and 99.4%, respectively).
Although GR in coarser soils (e.g., GRS) is generally
believed to be higher [Keese et al., 2005; Small, 2005],
Figure 3 clearly shows that the distributions of GR and ETa
overlap for sand and loamy sand. This is because soil
texture is solely defined by particle size distribution; otherfactors (e.g., bulk density, particle arrangement, and organic
matter) that determine SHCs also influence GR. Therefore,
the practice of using estimated mean SHCs would not be
sufficient to embody the whole range of GR, and so caution
is needed when computing GR on large scales withestimated mean SHCs [Keese et al., 2005; Small, 2005].
Additionally, the distributions of corresponding GR and
ETa look more or less symmetric because of negligible
surface runoff and change in soil moisture storage.
[22] GR estimated from mean SHCs (denoted as GR*) is
usually used as a regional estimate of GR [e.g., Nolan et al.,2003, 2007; Keese et al., 2005; Small, 2005]. On the basisof mean SHCs of the measured data set, GR* and ET*a, also
from mean SHCs, are given in Figure 3. GR*S and GR*LS
are 22.4 and 10.2 cm/a or 45.4% and 20.7% of the total P,
respectively. ET*
aS and ET*aLS are 26.8 and 39.0 cm/a or
54.4% and 79.1% of the total P, respectively. On the basis of
the results, GR*LS is lower than GR*S, which is consistent
with the results from Keese et al. [2005] and Small [2005].
However, Figure 3 also reveals that both GR*
S and GR*LS lie offfrom the peaks on the histograms (Figure 3), which correspond to the highest frequency of the computed GR. Particularly, GR*S occurs in neither of the two peaks (Figure 3a).
To further illustrate this, the averages of GR were calculated
from the distributions in Figure 3, with averages of 13.8 cm/afor sand and 9.5 cm/a for loamy sand. It is obvious that thedifference between the average GR and GR*S for sand islarger than that for loamy sand, and GR*S is thus lessrepresentative in a relative sense. This is due to the larger
standard deviation of the parameter n for sand (Table 2),
the reason of which will be explained in the following
section. In general, mean SHCs possess higher occurrence
probabilities; however, from Figure 3, using mean SHCs
does not imply that GR* also corresponds to the mode of
the GR distribution, which is due to the nonlinear relationships between GR and SHCs.
[23] Results from the PTF-generated correlated and uncorrelated SHC data sets are plotted in Figure 4. Because
the distributions of corresponding GR and ETa look also
symmetric for the correlated and uncorrelated SHC data

6 of 12

W04412WANG ET AL.: USING PEDOTRANSFER FUNCTIONS IN VADOSE ZONE MODELS

sets, the results of ETa are thus not presented. For thecorrelated data set, GRS ranges from 0.3 to 30.7 cm/a(e.g., <0.01% and 62.3% of the total P, respectively), and
GRLS ranges from 0.2 to 27.5 cm/a (e.g., <0.01% and55.8%, respectively). For the uncorrelated data set, GRS
ranges from 2.5 to 33.4 cm/a (e.g., 5.1% and 67.7%,
respectively), and GRLS ranges from 0.07 to 19.1 cm/a(e.g., 0.14% and 38.7%, respectively). In general, resultsfrom the measured and generated data sets exhibit similar
ranges and patterns of GR. GR*S and GR*LS are 21.5 and10.7 cm/a (e.g., 43.6% and 21.7% of the total P, respectively) for the correlated data set and 20.3 and 2.2 cm/a(e.g., 41.2% and 4.5%, respectively) for the uncorrelated
data set. In comparison to GR*S, GR*LS from the PTFgenerated data sets are very different. Furthermore, both
GR*S correspond to one of the two peaks in the histograms,
whereas GR*LS of the uncorrelated data set is in the clusteredpeak and lower than that of the correlated data set that lies
off from the peak (Figure 4). The difference between GR*LS
values is due to different soil databases used to generatemean SHCs [Schaap and Leij, 1998; Schaap et al., 2001],
which is also supported by the conclusion of Faust et al.
[2006] that estimated GR depends on the choice of PTFs.
The simulation results strongly suggest that the method ofusing PTF-estimated mean SHCs to calculate GR* as aregional representative may be inappropriate, particularly
for sand, as estimated GR* depends on the choice of PTFs
and original soil databases used to generate mean SHCs
as well.
[ 24 ] On the basis of a water balance method, theNebraska Natural Resources Commission [1986] estimatedthat GR accounted for 20 to 30% of annual precipitations in
the Sand Hills. Chen and Hu [2004] used inverse vadosezone modeling and found that GR ranged from 18% to 30%
of annual precipitations at a site similar to the GDEX site in
the Sand Hills. Also using a water balance method, Szilagyi
et al. [2005] estimated that mean GR was about 10% of
annual precipitations in the region surrounding the GDEX
site. In comparison, given that sand is prevalent at the
GDEX site [Wang et al., 2008], all calculated GR*S (>40%of the total P) are higher than the above results, even though
the humidity index of 0.31 in 2004 was already lower than
the annual mean of 0.42 in this area. Interestingly, the GR ratethat accounts for 10% of the mean annual precipitation[Szilagyi et al., 2005] is closer to the left peak on thehistograms (Figures 3a, 4a, and 4c), whereas GR*S is closer
to the right one. This also demonstrates that the practice ofusing mean SHCs to estimate GR is not appropriate, at leastin the study area.
4.1.2. Sensitivity Analysis of Groundwater Recharge to
Soil Hydraulic Characteristics[25] Figures 3 and 4 show that the shape of GR distributions varies significantly for different soil textures. In
general, the distribution of GRS exhibits a bimodal pattern,whereas the one of GRLS exhibits a unimodal pattern. The
standard deviations of GR from each data set are given in
Table 3, which also reveals that the distributions of GRS are
more variable than those of GRLS. To investigate the soil
parameters that mostly control the shape of GR distribution,
relationships between GR and SHCs (e.g., qr, qs, a, n, KS,and l) derived from the measured data set are plotted in
Figures 5 and 6. Note that the results from the PTF-

W04412

generated data sets exhibit a similar pattern in Figures 5
and 6 and therefore are not analyzed here.
[26] Of the van Genuchten parameters, qr, qs, and KS havephysical meanings, whereas the other three parameters areshape factors: a is inversely related to air entry pressure, n
is a measure of pore size distribution, and lumped parameterl accounts for pore tortuosity and connectivity. Strikingly,
compared to the shape factors, the distribution of GR is notvery sensitive to qr, qs, and KS for both sand and loamy sand(Figures 5a, 5b, 5e, 6a, 6b, and 6e), even though those threeparameters vary significantly in the measured data set. The
results indicate that qr, qs, and KS are relatively unimportantfor determining GR in semiarid regions, possibly because of
low soil moisture content in these regions. This conclusion,
however, may not hold in humid regions.
[27] The shape of GR distribution is instead mostly
controlled by the factor n (Figures 5d and 6d). Figures 5d
and 6d show that GR increases with increasing n values. In
general, a coarser soil with larger pore sizes possesses ahigher n value, which leads to a higher GR. Physically andmathematically, the factor n controls GR mainly throughaffecting soil surface evaporation. As n becomes larger, thedecrease in hydraulic conductivity with decreasing soil
moisture becomes stronger (equation (3)), which prevents
upward capillary flow of soil moisture toward the evaporating surface because of low hydraulic conductivity in thesurface layer. Soil moisture only a few centimeters below
the surface becomes unavailable for evaporation when thesurface layer is dry and n is large, which reduces evaporation and therefore increases GR. Additionally, the factorl has a similar effect as the factor n has on soil moisture flux
(equation (3)), but to a much lesser degree. On the basis of
inverse modeling and Monte Carlo simulations, Pollacco etal. [2008] also found that the factor n is crucial for
estimating GR under the unit gradient flow condition as
well as qs. However, we found that qs is unimportant for
determining GR in our study.
[28] It is important to note that the relationship between
GR and n is highly nonlinear (Figures 5d and 6d). When n
is approximately less than 4, GR increases rapidly withincreasing n, then GR starts to level off when n becomes
greater than 4, mostly because of limited precipitations(Figures 5d and 6d). The nonlinearity between GR and nreveals the reason why GR* is different from the GR withthe maximum occurrence in the GR distribution (Figure 3),
as soil samples with larger n values have more weights fordetermining GR*. This also causes the tendency of GR* tobe higher than the most frequent GR (Figure 3). In conclusion, the spatial variability in n is extremely important for
determining an areal-averaged GR and its spatial distribution in semiarid regions. The high uncertainty in n estimatedby PTFs, particularly for sand [Schaap et al., 1998], makes
the use of PTF-estimated mean SHCs in vadose zone
models less reliable.
4.2. Groundwater Recharge Under Constant Head
Condition[29] Most studies that used vadose zone models toestimate GR assumed the unit gradient flow at the lower
boundary of the vadose zone [Nolan et al., 2003, 2007;
Keese et al., 2005; Small, 2005]. However, this boundary
condition may not be valid for regions where water tables
are shallow because groundwater may return to the atmo-

7 of 12

W04412WANG ET AL.: USING PEDOTRANSFER FUNCTIONS IN VADOSE ZONE MODELS

W04412

Figure 5. Sensitivity analysis of groundwater recharge (GRS) and actual evapotranspiration (ETaS) of
sand to the van Genuchten parameters: (a) qr, (b) qs, (c) a, (d) n, (e) Ks, and (f) l. Pluses represent resultsfrom mean SHCs.
sphere through either capillary rise or groundwater evapotranspiration [Gillham, 1984; Wu et al., 1996; Batelaan etal., 2003]. Chen and Hu [2004] used a vadose zone model
to study the impact of groundwater on evapotranspiration ata similar site in the NSH. The authors found that the annual
evapotranspiration from 1998 to 2000 is 7– 21% higher
when the groundwater table is included in the model.
Therefore, to study the impact of lower boundary conditions
on GR, the measured and uncorrelated data sets were used
to compute GR. This task was not performed for thecorrelated data set because previous simulation results
showed similar patterns of the GR distributions for themeasured and correlated data sets. Constant heads were
set to be 0, 5, and 7 m at the lower boundary, whichcorrespond to zero constant heads at depths of 10, 5, and 3

m below the surface, respectively. The simulation results are
summarized in Figure 7.
[30] Figure 7 shows that a significant fraction of GR is
negative, which implies an upward soil moisture flux fromthe lower boundary. This obviously differs from the results
based on the unit gradient flow lower boundary condition
that enforces positive GR values. As explained previously,
the factor l has a similar effect on GR as the factor n has.
Additional simulations show that if the negative l values inboth measured and uncorrelated data sets are replaced by
0.5 as Mualem [1976] suggested, the negative GR values
would become ‘‘less’’ negative or positive. Schaap and Leij
[2000] noted the enormously wide range of l values in theliterature, and if l becomes more negative, hydraulic conductivity may start increasing with decreasing soil moisture

8 of 12

W04412WANG ET AL.: USING PEDOTRANSFER FUNCTIONS IN VADOSE ZONE MODELS

W04412

Figure 6. Sensitivity analysis of groundwater recharge (GRLS) and actual evapotranspiration (ETaLS) of
loamy sand to the van Genuchten parameters: (a) qr, (b) qs, (c) a, (d) n, (e) Ks, and (f) l. Pluses represent
results from mean SHCs.
content, which is physically unrealistic and can cause
unacceptable errors in mass balance as explained previously.
If a water table is present, selection of l value may become
another important issue. However, because of the wide
range of l values [Schaap and Leij, 2000], there is still noconsensus in the literature on how to select l. Hence, this
aspect requires additional study.
[31] Figure 7 also shows that GR with zero constant headat 10 m is higher and its distribution is less scattered thanthose with zero constant heads at 3 and 5 m depths.
Mathematically, a deeper constant head has less influenceon the upper atmospheric boundary and simulated domain,
allowing more soil moisture to pass downward across theroot zone. Table 4 shows GR* under different lower
boundary conditions. GR*S and GR*LS generally are more

intense under deeper constant head conditions. By varying
the lower boundary conditions, the variability in GR*S issmaller than the one in GR*LS. Therefore, the degree of the
impact of the lower boundary on the GR distribution
depends on soil texture as well. For regions with shallow
Table 3. Standard Deviations of Calculated Groundwater
Recharge for Different SHC Data Sets Under the Unit Gradient
Flow Boundary Condition
SHC Data Set

Sand (cm/a)

Loamy Sand (cm/a)

Measured
Uncorrelated
Correlated

9.15
9.38
8.68

8.31
4.90
5.28

9 of 12WANG ET AL.: USING PEDOTRANSFER FUNCTIONS IN VADOSE ZONE MODELS

W04412

W04412

Figure 7. Distributions of groundwater recharge (GR) under the constant head conditions: (a) sand and(b) loamy sand from the measured data set and (c) sand and (d) loamy sand from the uncorrelated data
set. Bin size is 5 cm/a.
water depths, as mimicked by the constant head condition inthis study, choosing an appropriate lower boundary condition is more important than selecting values of SHCs toestimate GR.

5. Conclusions[32] In this study, we assessed the applicability of coupling process-based vadose zone models and pedotransfer
functions for soil hydraulic characteristics to estimate
groundwater recharge in a semiarid region, especially theviability of using PTF-estimated mean SHCs (e.g., soil
hydraulic characteristics) to estimate GR [Nolan et al.,2003, 2007; Keese et al., 2005; Small, 2005]. Field hydrometeorological and physiological data sets collected in the
Nebraska Sand Hills were used to drive the vadose zone
model. Three SHC data sets containing the van Genuchten
parameters for both sand and loamy sand were used,
including a measured, a PTF-generated correlated, and a
PTF-generated uncorrelated data set.
[33] Under the unit gradient flow lower boundary condition, the simulation results show that there is a significant
overlap in the GR distributions for sand and loamy sand.
More importantly, GR* that is estimated from mean SHCs
generally lies off from the most frequent range of GR.
Therefore, the practice of using mean SHCs to estimate GR
is not sufficient, particularly for evaluating GR on large
scales [Nolan et al., 2003, 2007; Keese et al., 2005; Small,2005]. It is also found that GR mostly depends on the shapefactor n. Coarser soils with larger n values produce higher
GR; however, this relationship is highly nonlinear. Becauseof this nonlinear relationship, GR* tends to be higher thanmost simulated GR and therefore does not correspond to
GR distribution peaks. The spatial variability in n isimportant for determining an areal-averaged GR and its
spatial distribution in semiarid regions. The high uncertainty
in estimated n by PTFs, particularly for sand [Schaap et al.,1998], makes the use of PTF-estimated mean SHCs invadose zone models less reliable. When water tables arepresent as represented by the constant head conditions inthis study, choosing an appropriate lower boundary condition is more important than selecting values of SHCs toestimate GR. Upward soil moisture fluxes may exist from
Table 4. Groundwater Recharge Rate for Different Lower
Boundary Conditionsa
Sand (cm/a)
Lower Boundary
Condition
Unit gradient flow
Constant head at 3 m
Constant head at 5 m
Constant head at 10 m

Loamy Sand (cm/a)

Measured Uncorrelated Measured Uncorrelated
22.4
22.6
22.6
22.4

20.3
15.6
20.0
20.6

10.2
7.1
10.1
10.2

2.2
9.6
0.3
2.1

a
The results are based on mean SHCs of the measured and uncorrelated
data sets.

10 of 12

W04412WANG ET AL.: USING PEDOTRANSFER FUNCTIONS IN VADOSE ZONE MODELS

water tables, which is largely controlled by the factor n as
well as by the factor l. However, further studies are neededon the criteria for selecting l values. In general, GR is higher
and its distribution is less scattered under deeper constanthead boundary conditions, the degree of which also depends
on soil texture.
[34] Acknowledgments. This study was supported by the NSF Biocomplexity program grant DEB 387 0322067 (V. A. Zlotnik and T. Wang),
the NSF Hydrologic Sciences Program grants EAR-0609982 (V. A. Zlotnik)
and, partially, EAR-0737945 (M. G. Schaap). The authors acknowledge theinitial technical review by N. D. Smith (UNL) and valuable suggestions from
Editor B. Berkowitz, Associate Editor F. Laio, and three anonymous
reviewers that led to significant improvements of this work. The hydrometeorological data sets and soil data sets can be obtained from thecorresponding author upon request.

References
Allen, R. G., L. S. Pereira, D. Raes, and M. Smith (1998), Crop evapotranspiration: Guidelines for computing crop water requirements, Irrig.
Drain. Pap. 56, U. N. Food and Agric. Organ., Rome.
Anuraga, T. S., L. Ruiz, M. S. Mohan Kumar, M. Sekhar, and A. Leijnse(2006), Estimating groundwater recharge using land use and soil data: A
case study in South India, Agric. Water Manage., 84, 65 – 76,
doi:10.1016/j.agwat.2006.01.017.
Batelaan, O., F. De Smedt, and L. Triest (2003), Regional groundwater
discharge: Phreatophyte mapping, groundwater modelling and impact
analysis of land-use change, J. Hydrol., 275, 86 – 108, doi:10.1016/
S0022-1694(03)00018-0.
Carsel, R. F., and R. S. Parrish (1988), Developing joint probability distributions of soil water retention characteristics, Water Resour. Res.,
24(5), 755 – 769, doi:10.1029/WR024i005p00755.
Chen, X., and Q. Hu (2004), Groundwater influences on soil moisture andsurface evaporation, J. Hydrol., 297, 285 – 300, doi:10.1016/j.jhydrol.
2004.04.019.
Dann, R. L., M. E. Close, R. Lee, and L. Pang (2006), Impact of data
quality and model complexity on prediction of pesticide leaching,
J. Environ. Qual., 35, 628 – 640, doi:10.2134/jeq2005.0257.
Demirkanli, D. I., F. J. Molz, D. I. Kaplan, and R. A. Fjeld (2008), A fullytransient model for long-term plutonium transport in the Savannah riversite vadose zone: Root water uptake, Vadose Zone J., 7, 1099 – 1109,
doi:10.2136/vzj2007.0134.
Durner, W., B. Schultze, and T. Zurmül (1999), State-of-the-art in inverse
modeling of inflow/outflow experiments, in Characterization and Measurements of the Hydraulic Properties of Unsaturated Porous Media,
edited by M. T. van Genuchten, F. J. Leij, and L. Wu, pp. 661 – 681,
Univ. of Calif., Riverside, Calif.
Faust, A. E., T. P. A. Ferré, M. G. Schaap, and A. C. Hinnell (2006), Can
basin-scale recharge be estimated reasonably with water-balance models?, Vadose Zone J., 5, 850 – 855, doi:10.2136/vzj2005.0109.
Feddes, R. A., P. J. Kowalik, and S. P. Neuman (1978), Simulation of Field
Water Use and Crop Yield, John Wiley, New York.
Gillham, R. W. (1984), The capillary fringe and its effect on water-table
response, J. Hydrol., 67, 307 – 324, doi:10.1016/0022-1694(84)90248-8.
Hendrickx, J., and G. Walker (1997), Recharge from precipitation, in
Recharge of Phreatic Aquifers in (Semi-) Arid Areas, edited by I. Simmers,
pp. 19 – 98, A. A. Balkema, Rotterdam, Netherlands.
Jackson, R. B., J. Canadell, J. R. Ehleringer, H. A. Mooney, O. E. Sala, andE. D. Schulze (1996), A global analysis of root distributions for terrestrial
biomes, Oecologia, 108, 389 – 411, doi:10.1007/BF00333714.
Keese, K. E., B. R. Scanlon, and R. C. Reedy (2005), Assessing controls on
diffuse groundwater recharge using unsaturated flow modeling, Water
Resour. Res., 41, W06010, doi:10.1029/2004WR003841.
Lerner, D. N., A. S. Issar, and I. Simmers (1990), Groundwater recharge, aguide to understanding and estimating natural recharge, Rep. 8, 345 pp.,Int. Assoc. of Hydrogeol., Kenilworth, U. K.
Loope, D., and J. Swinehart (2000), Thinking like a dune field: Geologic
history in the Nebraska Sand Hills, Great Plains Res., 10, 5 – 35.
Miller, G. R., D. D. Baldocchi, B. E. Law, and T. Meyers (2007), An
analysis of soil moisture dynamics using multi-year data from a networkof micrometeorological observation sites, Adv. Water Resour., 30, 1065 –
1081, doi:10.1016/j.advwatres.2006.10.002.
Mualem, Y. (1976), A new model for predicting the hydraulic conductivityof unsaturated porous media, Water Resour. Res., 12(3), 513 – 522,
doi:10.1029/WR012i003p00513.

W04412National Research Council (2004), Groundwater Fluxes Across Interfaces,
Natl. Acad. Press, Washington, D. C.
Nebraska Natural Resources Commission (1986), Policy issue study on
integrated management of surface water and groundwater, 45 pp.,
Lincoln, Nebr.
Nemes, A., M. G. Schaap, F. J. Leij, and J. H. M. Wösten (2001), Description of the unsaturated soil hydraulic database UNSODA version 2.0,
J. Hydrol., 251, 151 – 162, doi:10.1016/S0022-1694(01)00465-6.
Neuman, S. P., R. A. Feddes, and E. Bresler (1974), Finite element simulation of flow in saturated-unsaturated soils considering water uptake byplants, third annual report, Proj. A10-SWC-77, Hydraul. Eng. Lab.,
Technion, Haifa, Israel.
Nolan, B. T., A. L. Baehr, and L. J. Kauffman (2003), Spatial variability ofgroundwater recharge and its effect on shallow groundwater quality insouthern New Jersey, Vadose Zone J, 2, 677 – 691.
Nolan, B. T., R. W. Healy, P. E. Taber, K. Perkins, K. J. Hitt, and D. M.
Wolock (2007), Factors influencing ground-water recharge in the eastern
United States, J. Hydrol., 332, 187 – 205, doi:10.1016/j.jhydrol.2006.
06.029.
Pollacco, J. A. P., J. M. S. Ugalde, R. Angulo-Jaramillo, I. Braud, and
B. Saugier (2008), A linking test to reduce the number of hydraulicparameters necessary to simulate groundwater recharge in unsaturated
soils, Adv. Water Resour., 31, 355 – 369, doi:10.1016/j.advwatres.2007.
09.002.
Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling (1988),
Numerical Recipes in C: The Art of Scientific Computing, Cambridge
Univ. Press, Cambridge, U. K.
Ritchie, J. T. (1972), Model for predicting evaporation from a row crop withincomplete cover, Water Resour. Res., 8(5), 1204 – 1213, doi:10.1029/
WR008i005p01204.
Scanlon, B. R., R. W. Healy, and P. G. Cook (2002), Choosing appropriatetechniques for quantifying groundwater recharge, Hydrogeol. J., 10, 18 –39, doi:10.1007/s10040-001-0176-2.
Scanlon, B. R., K. E. Keese, A. L. Flint, L. E. Flint, C. B. Gaye, W. M.
Edmunds, and I. Simmers (2006), Global synthesis of groundwater recharge in semiarid and arid regions, Hydrol. Processes, 20, 3335 – 3370,
doi:10.1002/hyp.6335.
Schaap, M. G., and F. J. Leij (1998), Database-related accuracy and uncertainty of pedotransfer functions, Soil Sci., 163, 765 – 779, doi:10.1097/
00010694-199810000-00001.
Schaap, M. G., and F. J. Leij (2000), Improved predictions of unsaturated
hydraulic conductivity with the Mualem – van Genuchten model, Soil Sci.
Soc. Am. J., 64, 843 – 851.
Schaap, M. G., and M. T. van Genuchten (2005), A modified Mualem – vanGenuchten formulation for improved description of the hydraulic conductivity near Saturation, Vadose Zone J., 5, 27 – 34, doi:10.2136/
vzj2005.0005.
Schaap, M. G., F. J. Leij, and M. T. van Genuchten (1998), Neural network
analysis for hierarchical prediction of soil hydraulic properties, Soil Sci.
Soc. Am. J., 62, 847 – 855.
Schaap, M. G., F. J. Leij, and M. T. van Genuchten (2001), ROSETTA: A
computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions, J. Hydrol., 251, 163 – 176, doi:10.1016/
S0022-1694(01)00466-8.
Schenk, H. J., and R. B. Jackson (2002), The global biogeography of roots,
Ecol. Monogr., 72(3), 311 – 328.
Segal, E., T. Kushnir, Y. Mualem, and U. Shani (2008), Water uptake andhydraulics of the root hair rhizosphere, Vadose Zone J., 7, 1027 – 1034,
doi:10.2136/vzj2007.0122.
Simunek, J., M. T. van Genuchten, and M. Sejna (2005), The HYDRUS-1D
software package for simulating the one-dimensional movement of water,
heat, and multiple solutes in variably saturated media, version 3.0,
HYDRUS Software Ser. 1, 270 pp., Dep. of Environ. Sci., Univ. of Calif.,
Riverside, Calif.
Small, E. E. (2005), Climatic controls on diffuse groundwater recharge insemiarid environments of the southwestern United States, Water Resour.
Res., 41, W04012, doi:10.1029/2004WR003193.
Szilagyi, J., E. F. Harvey, and J. F. Ayers (2005), Regional estimation of
total recharge to ground water in Nebraska, Ground Water, 43, 63 – 69,
doi:10.1111/j.1745-6584.2005.tb02286.x.
van Genuchten, M. T. (1980), A closed-form equation for predicting thehydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44,892 – 898.
Wang, T., V. A. Zlotnik, D. Wedin, and K. D. Wally (2008), Spatial trends
in saturated hydraulic conductivity of vegetated dunes in the Nebraska
Sand Hills: Effects of depth and topography, J. Hydrol., 349, 88 – 97,
doi:10.1016/j.jhydrol.2007.10.027.

11 of 12

W04412WANG ET AL.: USING PEDOTRANSFER FUNCTIONS IN VADOSE ZONE MODELS

Wesseling, J. G. (1991), Meerjarige simulatie van grondwaterstroming voor
verschillende bodemprofielen, grondwatertrappen en gewassen met het
model SWATRE, report, DLO-Staring Cent. Wageningen, Netherlands,
Rap., 152, 63 pp.
Wösten, J. H. M., Y. A. Pachepsky, and W. J. Rawls (2001), Pedotransfer
functions: Bridging the gap between available basic soil data and missing
soil hydraulic characteristics, J. Hydrol., 251, 123 – 150, doi:10.1016/
S0022-1694(01)00464-4.
Wu, J., R. Zhang, and J. Yang (1996), Analysis of rainfall-recharge relationships, J. Hydrol., 177, 143 – 160, doi:10.1016/0022-1694(95)02935-4.
Zlotnik, V. A., T. Wang, J. L. Nieber, and J. A. Simunek (2007), Verification of numerical solutions of the Richards equation using a traveling wave

W04412solution, Adv. Water Resour., 30, 1973 – 1980, doi:10.1016/j.advwatres.
2007.03.008.




M. G. Schaap, College of Agriculture and Life Sciences, University of
Arizona, 524 Shantz Building, Tucson, AZ 85721, USA. (mschaap@cals.
arizona.edu)
J. Šimunek, Department of Environmental Sciences, University ofCalifornia, Geology Building 2320, Riverside, CA 92507, USA. (jiri.
simunek@ucr.edu)
T. Wang, School of Natural Resources, University of Nebraska at
Lincoln, 521 Hardin Hall, Lincoln, NE 68583, USA. (twang3@unl.edu)
V. A. Zlotnik, Department of Geosciences, University of Nebraska at
Lincoln, 214 Bessey Hall, Lincoln, NE 68588, USA. (vzlotnik1@unl.edu)

12 of 12