Advances in Water Resources 24 (2001) 133±144
www.elsevier.com/locate/advwatresEect of the shape of the soil hydraulic functions near saturation onvariably-saturated ¯ow predictions
T. Vogel a, M. Th. van Genuchten b,*, M. Cislerova ca

Department of Hydraulics and Hydrology, Czech Technical University, Thakurova 7, Prague 16629, Czech RepublicbGeorge E. Brown, Jr., Salinity Laboratory, USDA, ARS, 450 West Big Springs Road, Riverside, CA, USAc
Department of Irrigation, Drainage and Landscape Engineering, Czech Technical University, Thakurova 7, Prague 16629, Czech Republic
Received 8 May 1998; received in revised form 10 June 2000; accepted 24 June 2000

Abstract
Relatively small changes in the shape of the soil water retention curve near saturation can signi®cantly aect the results ofnumerical simulations of variably saturated ¯ow, including the performance of the numerical scheme itself in terms of stability andrate of convergence. In this paper, we use a modi®ed form of the van Genuchten±Mualem (VGM) soil hydraulic functions toaccount for a very small, but non-zero minimum capillary height, hs , in the soil water retention curve. The modi®ed VGM model iscontrasted with the original formulation by comparing simulation results for in®ltration in homogeneous soils assuming bothconstant pressure and constant ¯ux boundary conditions. The two models gave signi®cantly dierent results for in®ltration in ®netextured soils, even for hs -values as small as )1 cm. Incorporating a small minimum capillary height in the hydraulic properties leads
to less non-linearity in the hydraulic conductivity function near saturation and, because of this, to more stable numerical solutions of
the ¯ow equation. This study indicates an urgent need for experimental studies that assess the precise shape of the hydraulicconductivity curve near saturation, especially for relatively ®ne-textured soils. For one example we found considerable improvementin the predicted conductivity function when a value of )2 cm for hs was used in the modi®ed VGM model. Ó 2000 Elsevier Science
Ltd. All rights reserved.
Keywords: Soil hydraulic properties; Unsaturated ¯ow; Numerical solutions

1. Introduction
Numerical models simulating the movement of waterand chemicals in variably saturated soils are increasingly
used in a variety of research and engineering projects.
Unsaturated ¯ow is usually described with the Richards'
equation which, for one-dimensional vertical ¯ow, isgiven by


oh h
o
oh

K h
K h ;
1ot
oz
oz
where h is the volumetric water content, h the soil waterpressure head, K the hydraulic conductivity, z the distance (positive downward), and t is the time. The use of
Eq. (1) requires knowledge of two soil hydraulic functions: the soil water retention curve, h h, and the un*

Corresponding author. Tel.: +1-909-369-4847; fax: +1-909-3424964.
E-mail addresses: vogel@fsv.cvut.cz (T. Vogel), rvang@ussl.ars.usda.gov (M.T. van Genuchten), cislerova@fsv.cvut.cz (M. Cislerova).

saturated hydraulic conductivity function, K(h).
Accurate measurement of the unsaturated conductivity
is generally cumbersome, costly, and very time-consuming. Consequently, many attempts have been madeto develop indirect methods which predict the conductivity function from the more easier measured water
retention curve [2,5,10,16,]. Most of the predictive conductivity models are based on the assumption of having
an ideal capillary medium characterized by a certain
pore-size distribution function. Although not necessarily
valid for all soils, this assumption is widely accepted asan eective working hypothesis [16].
A variety of empirical equations have been used todescribe the soil water retention curve. One of the mostpopular equations is the power-law function initially
proposed by Brooks and Corey [1]. This equation leads
to an air-entry value, ha , in the soil water retention
curve above which the soil is assumed to be saturated.
The Brooks and Corey equation was modi®ed by van
Genuchten [26] to enable a more accurate description ofobserved soil hydraulic data near saturation, especially

0309-1708/00/$ - see front matter Ó 2000 Elsevier Science Ltd. All rights reserved.
PII: S 0 3 0 9 - 1 7 0 8 ( 0 0 ) 0 0 0 3 7 - 3

134

T. Vogel et al. / Advances in Water Resources 24 (2001) 133±144

for undisturbed and many ®ne-textured soils. As shownby van Genuchten and Nielsen [27], Vogel et al. [29] andVogel and Cislerova [30], among others, the choice of
the analytical model for h h can signi®cantly aect thepredicted K(h) function obtained with one of the statistical pore-size distribution models. One reason for thisis that the predicted conductivity function is extremely
sensitive to small changes in the shape of the retention
curve near saturation. This sensitivity is a major cause of
the sometimes signi®cant dierences between predicted
K(h) functions obtained with the Brooks and Corey andvan Genuchten retention functions. The dierences are
specially important for ®ne-textured soils which can
exhibit extreme non-linearity in K(h) close to saturation
when van Genuchten's equations are used [27,28]. The
dierences are generally much less severe for coarsetextured soils.
As will be shown further in this paper, the presence ofhighly non-linear K(h) relationships near saturation canalso substantially impact the performance of numericalsolutions of Eq. (1) in terms of the accuracy, stability,
and rate of convergence of the invoked numerical
scheme. While numerical solutions for solving the variably saturated ¯ow equation have been signi®cantly
improved in recent years [4,6,9,14], most improvements
focused on numerical problems associated with the in®ltration of water in very dry coarse-textured soils[6,7,25]. By comparison, less attention has been paid tonumerical problems encountered when simulating in®ltration in ®ne-textured soils.
In this paper, we focus on the non-linear descriptionof the soil hydraulic properties near saturation. Specifically, we address two closely related issues: (1) how
uncertainty in the shape of the hydraulic functions,notably the hydraulic conductivity, near saturationtranslates into uncertainty in predicting soil water contents during in®ltration, and (2) how extremely smallchanges in the shape of the water retention curve nearsaturation can impact the stability and convergence ofnumerical solutions of Eq. (1). Results will be presented
in terms of van Genuchten's [26] expressions for the soil
hydraulic properties, as well as modi®cations thereof.

2. Description of the unsaturated soil hydraulic properties
The predictive model of Mualem [16] for the relative
hydraulic conductivity function, Kr (h), may be written in
the form:
 Z Se2 , Z 1
2
dx
dx
l
;
2
Kr S e   S eh x
0 h x
0where Kr  K/Ks , K is the unsaturated hydraulic conductivity, Ks the saturated hydraulic conductivity, l thepore connectivity parameter usually assumed to be 0.5

following Mualem [16], and Se is the eective saturation
given bySe h h h hr
;
hs hr

3

where h is the volumetric water content, and hr and hs are
the residual and saturated water contents, respectively.
Substituting van Genuchten's [26] expression for the soil
water retention curve, i.e.,hr  1hs ahhrn m ; h < 0;4
h h 
h P 0:
hs ;into Eq. (2) leads to the following equation for the unsaturated hydraulic conductivity

Ks Kr h; h < 0;5
K h 
h P 0;
Ks ;
where


Kr Se   Sel 1

Se1=m

1

m  26

in which n and a are empirical shape parameters (withn > 1), and m  1 1=n. The restriction m  1 1=n isnecessary to permit direct integration of Eq. (2) once theinverse of Eq. (4) is substituted into Eq. (2). Much more
complicated expressions for Kr (h) result when theparameters m and n are assumed to be mutually independent [28]. We will refer to Eqs. (4)±(6) as the original
VGM model.
To increase the ¯exibility of the VGM model in describing retention and hydraulic conductivity data nearsaturation, the following modi®cations in Eqs. (4)±(6)were suggested by Vogel et al. [29] and Vogel andCislerova [30]:
hr  1hm ahhnr m ; h < hs ;h h 
7
hs ;h P hs ;Ks Kr h; h < hs ;8K h Ks ;h P hs ;where

Kr Se   Sel

1
1

F Se   1

Se1=m

Se 

hs
hm

hr
Se
hr

F Se 
F 1
m

;

2
;

9
10
11

in which as before m  1 1=n and n > 1. The parameter hs is denoted here as the minimum capillary height.
Notice from Fig. 1 that the parameter hs in the modi®edretention function is replaced by a ®ctitious (extrapolated) parameter hm P hs . This approach maintains
the physical meaning of hs as a measurable quantity,
while the de®nition of the eective degree of saturation,
Se , is also not aected. The minimum capillary height,

T. Vogel et al. / Advances in Water Resources 24 (2001) 133±144

135

should have some non-zero minimum capillary height(hs < 0) associated with the largest pore of the soil. For
example, the equation of capillary rise (e.g., [8, p. 23])shows that water in a capillary of radius 1 mm shouldhave a capillary elevation of approximately 15 mm.
While it may appear here that we are trying to assign
physical signi®cance to hs , this coecient in reality isprobably best treated pragmatically as an additional®tting parameter that enables a more precise descriptionof the retention and hydraulic conductivity curves nearsaturation.
3. Shape of the hydraulic functions near saturation

Fig. 1. Schematics of the soil water retention and hydraulic conductivity functions as described with the modi®ed VGM model.

hs , becomes zero when hm  hs , in which case Eqs. (7)±(11) reduce to the original VGM model.
Plots of the modi®ed VGM model in Fig. 1 suggest
that, at least phenomenologically, the parameter hs
serves the same purpose as the air-entry value, ha , in thesoil hydraulic equations of Brooks and Corey [1], that is
to introduce a break in the functions at some negative
pressure head. This is only partially true. Whereas
Brooks and CoreyÕs air entry value is usually somewhere
between )10 and )1000 cm depending upon soil type(e.g., [24]), hs in our study will be much smaller in absolute value (e.g., )1 to )2 cm) so that the retentionfunction will retains its more-or-less S-shaped form
typical of the VGM model (this will be further discussed
later). The parameter hs is a direct consequence of thesoil system being a capillary medium which always

We now address the question as to how well theoriginal and modi®ed VGM equations can describe thesoil hydraulic functions close to saturation. The presence in Eq. (7) of a non-zero minimum capillary height
has little or no eect on the retention function, but
strongly aects the predicted shape of the unsaturated
hydraulic conductivity function. This is demonstrated
by a comparison of Figs. 2 and 3 which, for the original
and modi®ed VGM models, respectively, show variations in the soil hydraulic characteristics for a series ofn-values, while keeping all other parameters constant(a  0:05 cm 1 ; hs  0:40 and hr  0:10). Notice that the
Kr (h) function in Fig. 2 exhibits an abrupt drop at
saturation when n becomes less than about 1.5. The
abrupt near-instantaneous and highly non-linear decrease in Kr (h) is caused by the 1=h term in Eq. (2)which goes to in®nity at Se  1, and hence leads to a
singularity in the original VGM conductivity equation.
The magnitude of the decrease at h  0 depends not
only on the parameter n, but to a lesser extent also on a,which may be viewed as a scaling factor for h. The
extreme non-linearity occurs only when 1 < n < 2.
Close inspection of Figs. 2(b) and (c) shows that the soil
water capacity function, C h  dh=dh, and the relative
conductivity function, Kr (h), both change their shapenear saturation when n changes from n < 2 to n > 2.
The slopes of these two functions change from 1 forC(h) and 1 for Kr (h) when n < 2 to some non-zero ®nite values when n  2, and to zero when n > 2. Thiseect of the value of n on the shape of the retention
curve near saturation was discussed previously by
Luckner et al. [12] and Nielsen and Luckner [20]. They
showed that the restriction n > 1 guarantees only ®rstorder continuity in h h at saturation, while n > 2 ensures also second-order continuity in h h, and hence®rst-order continuity in C(h).
A comparison of Figs. 2 and 3 shows that the introduction of a very small but non-zero value of the minimum capillary height, hs , in the VGM model leads tosigni®cant dierences in the predicted K(h) functions forn values close to 1 (i.e., for many ®ne-textured and/or

136

T. Vogel et al. / Advances in Water Resources 24 (2001) 133±144

Fig. 2. Plots of the soil hydraulic properties near saturation as a function of the parameter n (original VGM model).

undisturbed soils characterized by relatively broadpore-size distributions), and to far less dramatic or eveninsigni®cant dierences when n is large (i.e., for morecoarse-textured, disturbed and/or compacted soils having relatively narrow pore-size distributions). By comparison, corresponding dierences between the retention
curves for the two models (Figs. 2(a) and 3(a)) are
negligible, and de®nitely within the range of routine
measurement errors. However, notice from Fig. 3(b)that the modi®ed VGM soil water capacity curve is now
discontinuous at h  hs . Also, Fig. 3(c) shows that themodi®ed hydraulic conductivity curve now no longer
exhibits the abrupt decrease when the pressure head
becomes equal to hs < 0. This result is due to the factthat Mualem-type unsaturated conductivity functions
are obtained by evaluating the integrals in Eq. (2) from 0
to Se (or 1), which is equivalent to a pressure head range
from 1 to h (or 0). When hs < 0, the singularity in
Eq. (2) at h  0 is avoided and integration will always
lead to a only modestly increasing K(h) curve near saturation. In other words, the slope K 0 h of the conductivity function remains ®nite when h reaches any hs < 0,as opposed to becoming in®nite when h ! hs  0 forn < 2 as shown in Fig. 2(c).
The eects of having a small non-zero hs on the hydraulic properties is illustrated here with one example.

Fig. 4 shows measured water retention and hydraulicconductivity data for Beit Netofa clay as obtained by
Rawitz [23] and cited in Mualem [17]. The data are described in terms of both the original and modi®ed VGM
equations. The data set is well suited for our purposes
here because several unsaturated conductivity data
points were measured close to saturation. Hydraulicparameters for Eqs. (4) and (7) are given in Table 1. The
parameters a and n were optimized using least-squares[28]. The parameter hm was adjusted so that the value of
hs in (7) would become equal to )2 cm. Note that, as
expected, the original and modi®ed VGM retention
curves resulting are visually indistinguishable, whereasdierences in the predicted hydraulic conductivity
functions are signi®cant and manifested over the entire
range of pressure heads (i.e., not only near saturation).
From Figs. 4(a) and (b), it is evident that the modi®ed
VGM method for this clay soil leads to a much better
agreement between the predicted and measured conductivity curves.
4. Numerical simulations
The eects of changes in the shape of the soil hydraulic functions near saturation on simulated pressure

T. Vogel et al. / Advances in Water Resources 24 (2001) 133±144

137

Fig. 3. Plots of the soil hydraulic properties near saturation as a function of the parameter n (modi®ed VGM model with hs  )2 cm).

head pro®les were further analyzed by numericallysimulating in®ltration into an initially unsaturated
homogeneous soil column. Results were obtained with
the HYDRUS (version 5.0) one-dimensional Galerkin®nite element computer program documented by Vogelet al. [31]. The solution invokes mass lumping of thetime derivative in the ¯ow equation, and incorporates
the mass-conservative numerical scheme of Celia et al.
[4]. We considered in®ltration in a 100-cm long soilcolumn that initially was assumed to be in equilibriumwith an imposed pressure head of )1000 cm at thebottom of the column. Calculations were performed
using nodal spacing of 0.5 and 1.0 cm, depending upon
the non-linearity of the simulated problem, and an
adaptive time stepping algorithm as described by Vogelet al. [31]. Three combinations of boundary conditions
were tested: (1) a zero pressure head at the bottom of
the column combined with a zero ¯ux at the top of thecolumn, leading to upward in®ltration against gravity;(2) a zero pressure head at the top of the column
combined with a zero ¯ux at the bottom, leading to
downward in®ltration, and (3) constant ¯ux in®ltration
at the top of the soil pro®le, combined with a zero ¯ux
at the bottom. Calculations were carried out for manysoil types using n values between 1 and 1.5. Additional
simulations were carried out for upward in®ltration

using the soil hydraulic data listed in Table 2, as well
as by varying n while keeping the other hydraulicparameters constant. The data in Table 2 were obtained
by Carsel and Parrish [3] by analyzing a large databaseof observed soil water retention data. We believe thatthe entries in Table 2 provide reasonable estimates of
the VGM hydraulic parameters for dierent soil textural groups. As an illustrative example we present here
detailed results only for the clay soil (Soil no. 12) whose
n-value is 1.09. The hydraulic properties of this soil are
displayed in Fig. 5 in terms of the original VGM formulation with hs  0, and the modi®ed VGM modelwith hs  2 cm. Note again the insigni®cant dierences between the two h h curves, and the large differences between the Kr (h) curves over the entire rangeof h-values.
The simulation results for upward and downwardin®ltration are presented in Figs. 6 and 7, respectively,
while those for constant ¯ux in®ltration assuming asurface water ¯ux, q0 of Ks /2 are given in Fig. 8. The®gures compare results obtained with the original VGM
formulation with a zero minimum capillary height(hs  0), with calculations obtained with the modi®ed
VGM model assuming hs  2 cm. Figs. 6 and 7 show
signi®cant dierences in the position of wetting front
between the original and modi®ed hydraulic models.

138

T. Vogel et al. / Advances in Water Resources 24 (2001) 133±144

Fig. 4. Comparison of measured soil hydraulic data for Beit Netofa clay with calculated curves using: (i) the original VGM equations assuming azero minimum capillary height, hs  0; (ii) the modi®ed VGM with hs  )2 cm.

Table 1
Estimated hydraulic parameters for Beit Netofa clay
hr

hs

hm

a

n

Method

0.010
0.010

0.450
0.450

±
0.450080.001730.00170

1.170
1.172

Original VGM
Modi®ed VGM

The calculated in®ltration rates are also quite dierent.
These dierences are consistent with the dierent shapesof the original and modi®ed VGM hydraulic conductivity functions in Fig. 5.
Fig. 7(c) shows severe oscillations in the calculatedin®ltration rate for the hs  0 case when downward in®ltration with a constant pressure head boundary condition was simulated. Some oscillations also developedin the calculated pressure heads for the constant ¯uxboundary condition simulation (Fig. 8). Numerical instabilities associated with these oscillations are a directconsequence of the non-linear shape of the hydraulicconductivity function near saturation. For example, as
the in®ltration process shown in Fig. 7 for hs  0 proceeds and the top part of the soil column approaches
saturation, the pressure head goes to zero and the surface ¯ux slowly approaches the value of Ks . At this time
the numerical scheme becomes aected by the extremely
non-linear shape of the conductivity function near saturation (Fig. 5), and starts to produce oscillations.
When hs  2 cm, on the other hand, the conductivity
function decreases much more gradually and no oscillations occur. When downward in®ltration with a constant ¯ux boundary condition is simulated (Fig. 8), theT. Vogel et al. / Advances in Water Resources 24 (2001) 133±144

139

Table 2Average values of VGM soil hydraulic parameters for 12 major soil textural groups according to Carsel and Parrish, 1988 [3]Texturehr (m3 /m3 )hs (m3 /m3 )a (1/cm)n

Ks (cm/day)(1) Sand(2) Loamy sand(3) Sandy loam(4) Loam(5) Silt(6) Silt loam(7) Sandy clay loam(8) Clay loam(9) Silty clay loam(10) Sandy clay(11) Silty clay(12) Clay

0.045
0.057
0.065
0.078
0.034
0.067
0.100
0.095
0.089
0.100
0.070
0.068

0.430
0.410
0.410
0.430
0.460
0.450
0.390
0.410
0.430
0.380
0.360
0.380

0.145
0.124
0.075
0.036
0.016
0.020
0.059
0.019
0.010
0.027
0.005
0.008

2.68
2.28
1.89
1.56
1.37
1.41
1.48
1.31
1.23
1.23
1.09
1.09

712.80350.20106.10
24.96
6.00
10.80
31.44
6.24
1.68
2.88
0.48
4.80

Fig. 5. Soil hydraulic properties of a clay soil (Soil no. 12 in Table 2) assuming hs  0 and hs  )2 cm.

pressure head gradient again goes to zero as the top partof the column ®lls with water and the water ¯ux becomes equal to the imposed surface ¯ux, q0  Ks =2. The
steep part of the conductivity function close to Ks /2 now
again causes small pressure head oscillations (Fig. 8(c)).
By contrast, no oscillations occur for upward in®ltration(Fig. 6); this because the pressure head gradient in the

more saturated bottom part of the column goes to +1
instead of zero. The upward water ¯ux rate through thebottom of the column decreases quickly to values outside of the most non-linear portion of the Kr (h) curve in
Fig. 5, leading to numerically more stable simulations.
The pressure head distributions in Fig. 6 eventuallywould reach equilibrium.

140

T. Vogel et al. / Advances in Water Resources 24 (2001) 133±144

Fig. 6. Calculated upward in®ltration into a clay soil column (Soil no. 12) subject to a zero pressure head at the lower boundary.

To assess the relative importance of the eects of hs
on water ¯ow in dierent soils, we next carried out aseries of calculations for upward in®ltration assuming
for all cases minimum capillary heights, hs , of both 0 and)2 cm. The total volume of in®ltrated water after one
day of in®ltration (denoted by I and I for the original
and modi®ed VGM models) were calculated for all soil
types listed in Table 2. Fig. 9 gives a plot of the relativechange, I  I=I, in the in®ltrated amount of water asa function of the parameter n. Each pair of simulations
resulted in an asterisk in Fig. 9 (one for each soil in
Table 2). Relative changes in in®ltration were also calculated by ®xing the parameters a  0:01 cm; hs  0:40;
hr  0:10; Ks  10 cm day 1 , but letting the parametern vary between 1.0 and 3.0. The entire sequence wasrepeated once more with a  0:1 cm 1 . The results in
Fig. 9 indicate that a hyperbolic type function accurately
bounds the most sensitive areas of the numerical algorithm. The ®gure shows that dierences in the amount ofsimulated in®ltration resulting from the original andmodi®ed VGM models can reach more than one hundred percent for soils having n values close to 1.

5. Discussion and conclusions
In this study, we showed that extremely small changes
in the shape of the soil water retention curve near saturation can lead to signi®cantly dierent curves for theunsaturated hydraulic conductivity function when incorporated into the predictive Kr (h) model of Mualem[16]. The modi®ed conductivity curves can lead to remarkably dierent simulation results for ®ne-textured
soils. These ®ndings were shown to hold for both upward and downward ponded in®ltration, and downward
constant-¯ux in®ltration. The dierences are a consequence of the extreme sensitivity of MualemÕs model tosmall changes in h h close to saturation. The sensitivity
is relatively low for n > 2, but increases dramaticallywhen n approaches 1.0, i.e., for most ®ne-textured and/or certain undisturbed soils. The sensitivity also depends
to some extent on the value of scaling parameter a(Fig. 9).
We found that the non-linear nature of the original
VGM conductivity function with hs  0 can cause
serious numerical convergence problems for n values

T. Vogel et al. / Advances in Water Resources 24 (2001) 133±144

141

Fig. 7. Calculated downward in®ltration into a clay soil (Soil no. 12) subject to ponded in®ltration at the surface.

close to 1. Even if convergence is achieved, numericaloscillations may develop in the in®ltration rate for aconstant pressure head boundary condition, and in thecalculated pressure head when a constant ¯ux boundary
condition is imposed. The convergence and oscillation
problems were eectively removed by using the modi®ed
VGM model with a small non-zero minimum capillaryheight, hs .
This study indicates an urgent need for very precise
measurements of the unsaturated hydraulic conductivity function close to saturation for ®ne-textured
soils (i.e., for n values less than about 1.3). Theoreticalconductivity models, such as the Brooks±Corey andVGM models, are increasingly used in numericalmodels predicting subsurface water ¯ow and chemicaltransport. While the equations have been validated
against observed data for a large number of soils,
most or all of these have been medium-or coarse-textured soils [11,16,18,19,28]. The soils in many cases
were also repacked in the laboratory. In contrast,
there is a lack of validation for ®ne-textured soils,
mostly because of a conspicuous absence of reliable

experimental data in the literature. This is unfortunate
since, as shown in this paper, the greatest sensitivitiesof the predictive theories to changes in the retention
curve occur for such ®ne-textured soils. Current predictions of the hydraulic conductivity function using
the original VGM model may well be unreliable for nvalues close to 1.0. One example soil, Beit Netofa clay,
was used in this paper to illustrate the type of dierences that can be expected when the original VGM is
modi®ed to include a minimum capillary height, hs .
For this soil we used a value of )2 cm for hs . Preliminary analyses of the UNSODA unsaturated soil
hydraulic property database [11,19] involving some 235
combined retention and unsaturated hydraulic conductivity data sets showed that: (1) hs values between 0(including zero) and )5 cm produce essentially thesame h h function; (2) the use of any hs < 0 (i.e.,
excluding zero) has a very signi®cant eect on thecalculated K(h) function as compared to the original
MVG model; (3) the eect of changing hs on the K(h)curve of a particular soil is about the same over therange 10 < hs < 1 cm.

142

T. Vogel et al. / Advances in Water Resources 24 (2001) 133±144

Fig. 8. Calculated downward in®ltration into a clay soil (Soil no. 12) subject to ¯ux-controlled in®ltration at the soil surface (constant ¯ux,
q  Ks =2).

We emphasize here that our study focused on unstructured, mostly non-swelling soils. Additional complications arise for structured (macroporous or
fractured) media associated with the presence of decayedroot channels, earthworm channels, natural interaggregate pores, rock fractures, and/or drying cracks.
Alternative formulations may be needed for such
structured media, possibly using sums of several VGM
type or other models [13,15,21,22,32]. Like the introduction of slightly negative hs value, bimodal or multimodel functions lead to signi®cant changes in K(h) nearsaturation, but generally at slightly lower pressure heads(e.g., )3 to)10 cm, or lower, depending upon the type ofmacropores of fractures present). Conceptually, thereasons for such changes are of course completelydierent for structured as compared to ®ne-textured
porous media.
Based on this study, as well as ongoing work, we
recommend to always use the modi®ed VGM modelwith a small value (e.g., )2 cm) for hs when no precisely measured unsaturated hydraulic conductivity

data are available and the parameter n (as estimatedfrom retention data only) is less than about 1.3. Whenan optimization algorithm is used to determine the
VGM-model parameters (either by ®tting retention
data or by solving the inverse problem by ®tting thehydraulic parameters to measured time±space variations in pressure heads, water contents and/or ¯uxes),
we recommend to use an additional optimized parameter hm (Eq. (7)) which during the optimization should
remain slightly larger than hs (e.g. hm  hs  0:0001).
Alternatively, the following approximation may be
used for hm :h m  hr  hs

n m

hr 1  2a  ;12which follows from Eq. (7) with hs  2 cm. Our recommendation for hm is based in part on conceptualissues (including the Beit Netofa Clay example), and inpart on more pragmatic considerations to remove theextreme non-linearity in the predicted K(h) function, andhence to minimize numerical solution problems.

T. Vogel et al. / Advances in Water Resources 24 (2001) 133±144

Fig. 9. Calculated relative increase in the amount of in®ltrated water
when using a non-zero capillary height, hs , of )2 cm as compared tohs  0 (the asterisks represent results for the soils listed in Table 2; solidlines were obtained using two values of a in the VGM model andvarying the value of n between 1 and 3).

Acknowledgements
This work was supported in part by the Grant
Agency of the Czech Republic through the Academy of
Sciences of grant No. A3060001.
References[1] Brooks RH, Corey AT. Hydraulic properties of porous media.
Hydrology Paper No. 3, Fort Collins, Colorado: Colorado State
University; 1964; p. 27.
[2] Burdine NT. Relative permeability calculations from pore-size
distribution data. Petrol Trans Am Inst Min Eng 1953;198:71±7.
[3] Carsel RF, Parrish RS. Developing joint probability distributionsof soil water retention characteristics. Water Resour Res
1988;24(5):755±69.
[4] Celia MA, Bouloutas ET, Zarba R. A general mass-conservative
numerical solution for the unsaturated ¯ow equation. Water
Resour Res 1990;26(7):1483±96.
[5] Childs EC, Collis-George N. The permeability of porous materials. Proc Roy Soc 1950;201:392±405.
[6] Huang K, Mohanty BP, van Genuchten MT. A new convergencecriterion for the modi®ed Picard iteration method to solve thevariably saturated ¯ow equation. J Hydrol 1996;167:69±91.
[7] Hills RG, Hudson DB, Porro I, Wierenga PJ. Modeling onedimensional in®ltration into very dry soils, 1. Model development
and evaluation. Water Resour Res 1989;25(6):1259±69.
[8] Kirkham D, Powers WL. Advanced soil physics. New York:
Wiley±Interscience; 1972. p. 534.
[9] Kirkland MR, Hills RG, Wierenga PJ. Algorithms for solving
Richards' equation for variably saturated soils. Water Resour Res
1992;28(8):2049±58.

143[10] Leij FJ, Russell WB, Lesch SM. Closed-form expressions forwater retention and conductivity data. Ground Water
1996;35(5):848±58.
[11] Leij FJ, Alves WJ, van Genuchten MT, Williams JR. The
UNSODA Unsaturated Soil Hydraulic Database. User's Manual,
Version 1.0. EPA/600/R-96/095, Cincinnati, OH: National Risk
Management Laboratory, Oce of Research and Development,
US Environmental Protection Agency; 1996. p. 103.
[12] Luckner L, van Genuchten MT, Nielsen DR. A consistent set ofparametric models for the ¯ow of water and air as immiscible¯uids in the subsurface. Water Resour Res 1989;25(10):2187±93.
[13] Mallants D, Tseng P-H, Toride N, Timmerman A, Feyen J.
Evaluation of multimodal hydraulic functions in characterizing aheterogenous ®eld soil. J Hydrol 1997;195:172±99.
[14] Milly PCD. A mass-conservative procedure for time stepping inmodels of unsaturated ¯ow. Adv Water Resour 1985;8(1):32±6.
[15] Mohanty BP, Bowman RS, Hendrickx JMH, van Genuchten MT.
New piecewise-continuous hydraulic functions for modelingpreferential ¯ow in an intermittent ¯ood-irrigated ®eld. Water
Resour Res 1997;33:2049±63.
[16] Mualem Y. A new model for predicting the hydraulic conductivityof unsaturated porous media. Water Resour Res 1976;12:513±22.
[17] Mualem Y. A catalogue of the properties of unsaturated soils.
Research Report 442. Technion, Haifa, Israel: Israel Institute ofTechnology; 1976.
[18] Mualem Y. Modeling the hydraulic conductivity of unsaturated
porous media. In: van Genuchten MT, Leij FJ, Lund LJ, editors.
Proc. Int. Workshop on Indirect Methods for Estimating theHydraulic Properties of Unsaturated Soils. Riverside: Universityof California; 1992. p. 15±36.
[19] Nemes A, Schaap MG, Leij FJ. The UNSODA unsaturated soil
hydraulic property database, version 2.0. Available at http://
www.ussl.ars.usda.gov/MODELS/unsoda.htm, 2000.
[20] Nielsen DR, Luckner L. Theoretical aspects to estimate reasonable initial parameters and range limits in identi®cation procedures for soil hydraulic properties. In: van Genuchten MT, Leij
FJ, Lund LJ, editors. Proc. Int. Workshop on Indirect Methods
for Estimating the Hydraulic Properties of Unsaturated Soils,
Riverside: University of California; 1992. p. 147±60.
[21] Othmer H, Diekkruger B, Kutilek M. Bimodal porosity andunsaturated hydraulic conductivity. Soil Sci 1991;152:139±50.
[22] Peters RR, Klavetter EA. A continuum model for water movement in an unsaturated fractured rock mass. Water Resour Res
1988;24:416±30.
[23] Rawitz E. The in¯uence of a number of environmental factors on
the availability of soil moisture (in Hebrew), Ph.D. thesis,
Rehovot, Israel: Hebrew University; 1965.
[24] Rawls WJ, Ahuja LR, Brakensiek DL. Estimating soil hydraulicproperties from soils data. In: van Genuchten MT, Leij FJ, Lund
LJ, editors. Proc. Int. Workshop on Indirect Methods for
Estimating the Hydraulic Properties of Unsaturated Soils, Riverside: University of California; 1992. p. 329±340.
[25] Ross PJ. Ecient numerical methods for in®ltration using
Richards' equation. Water Resour Res 1990;26(2):279±90.
[26] van Genuchten MT. A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci Soc Am J
1980;44:892±8.
[27] van Genuchten MT, Nielsen DR. On describing and predicting
the hydraulic properties of unsaturated soils. Ann Geophys
1985;3:615±28.
[28] van Genuchten MT, Leij FJ, Yates SR. The RETC code forquantifying the hydraulic functions of unsaturated soils. Report
No. EPA/600/2-91/065, R. S. Kerr Environ. Res. Laboratory, US
Environmental Protection Agency, Ada, Oklahoma; 1991. p. 85.
[29] Vogel T, Cislerova M, Sir M. Reliability of indirect determinationof soil hydraulic conductivity. Vodohosp Cas 1985;33(2):204±24(in Czech).

144

T. Vogel et al. / Advances in Water Resources 24 (2001) 133±144[30] Vogel T, Cislerova M. On the reliability of unsaturated hydraulicconductivity calculated from the moisture retention curve. Transport Porous Media 1988;3:1±15.
[31] Vogel T, Huang K, Zhang R, vanGenuchten MT. The HYDRUS
code for simulating one-dimensional water ¯ow, solute transport,

and heat movement in variably-saturated media, version 5.0,
Research Report 140, Riverside, California: US Salinity Laboratory; 1996. p. 131.
[32] Zurm
uhl T, Durner W. Modeling transient water and solutetransport in a biporous soil. Water Resour Res 1996;32:819±29.