Published online December 16, 2005

A Modified Mualem–van Genuchten Formulation for Improved Descriptionof the Hydraulic Conductivity Near Saturation
Marcel G. Schaap* and Martinus Th. van Genuchten
in which Ko is a matching point at saturation (cm d21), L
is an empirical pore-connectivity parameter, and Se iseffective saturation given by

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

ABSTRACTThe unsaturated soil hydraulic properties are often described using
Mualem–van Genuchten (MVG) type analytical functions. Recent
studies suggest several shortcomings of these functions near saturation, notably the lack of second-order continuity of the soil water
retention function at saturation and the inability of the hydraulicconductivity function to account for macroporosity. We present a
modified MVG formulation that improves the description of thehydraulic conductivity near saturation. The modified model introducesa small but constant air-entry pressure (hs) into the water retention
curve. Analysis of the UNSODA soil hydraulic database revealed an
optimal value of 24 cm for hs, more or less independent of soil texture.
The modified model uses a pressure dependent piece-wise linear correction to ensure that deviations between measured and fitted conductivities between pressure heads of 0 and 240 cm were eliminated. A
small correction was found necessary between 24 and 240 cm, and amuch larger correction was needed between 0 and 24 cm. An average
RMSE in logK of only 0.26 remained for a data set of 235 samples. The
resulting modified MVG model was found to have small systematic
errors across the entire pressure range. The modified model appears
well suited for large-scale vadose zone flow and transport simulations,
including inverse modeling studies.

MSe (h) 5

quire description of unsaturated soil hydraulicproperties over a wide range of pressure heads, h. The
hydraulic properties are often described using the poresize distribution model of Mualem (1976) for thehydraulic conductivity in combination with a water retention function introduced by van Genuchten (1980).
The soil water retention equation, u(h), is given by

5

ur 1

u s 2 urn
[1 1 jahj ]m
us

h#0

[1]

h.0

where u is the volumetric water content (cm3 cm23) atpressure head h (cm); ur and us are the residual andsaturated water contents, respectively (cm3 cm23); a(.0, in cm21) is related to the inverse of the air-entry
pressure; n (.1) is a measure of the pore-size distribution (van Genuchten, 1980); and m 5 1 2 1/n. The corresponding MVG hydraulic conductivity function, K(h), isK(Se ) 55

Ko SLe [1 2 (1 2 S1/m
)m ]2
e
Ko

h#0
h.0

[3]

Accurate measurements of the unsaturated hydraulicconductivity are generally difficult and time-consuming.
When the water retention parameters are known, itis therefore common practice to use Eq. [2] to estimate
K(Se) or K(h) while making assumptions about thevalues of Ko and L. Most commonly, a measured saturated hydraulic conductivity (Ks) is used for Ko, while L is
fixed at 0.5 (Mualem, 1976). When water retention andunsaturated hydraulic conductivity data are available, an
optimization may be performed for several or all of theparameters (Ko, L, ur, us, a, and n) contained in Eq. [1]and [2] (e.g., van Genuchten et al., 1991). This may bedone by optimizing the parameters in Eq. [1] and [2]separately, sequentially, or with a simultaneous procedure in which the objective functions for Eq. [1] and [2]are combined (Yates et al., 1992). Some or all of theparameters in Eq. [1] and [2] are sometimes estimatedsimilarly by means of inverse procedures from measured
laboratory and/or field measurements of the pressure
head, the water content, and/or fluid fluxes (e.g.,
Hopmans et al., 2002; Butters and Duchateau, 2002).
While Eq. [1] and [2] have found widespread use, theyalso have several limitations caused by the particularmathematical properties of Eq. [1] or by the use of default values for Ko and L. Luckner et al. (1989) and
Vogel et al. (2000), among others, showed that the retention curve lacks second-order continuity at saturation
when n , 2 (i.e., d2u/dh2 becomes discontinuous at
h 5 0). If combined with Eq. [2], the resulting conductivity function, K(h) will then exhibit a sharp dropat pressures just below saturation, especially whenn approaches its lower limit of 1.0 (Vogel et al., 2000).
For 1.0 , n , 1.3, the reduction in conductivity may
become so dramatic that it may cause numerical instabilities in simulations of near-saturated infiltration.
Vogel et al. (1985, 2000) solved this problem by forcing
the retention function to be constant at us between saturation (h 5 0) and some very small negative value of thepressure head, hs. While only minimally affecting the soil
water retention curve macroscopically, the modification
significantly reduced the nonlinearity in the hydraulicconductivity function near saturation. Vogel et al. (2000)

ANY VADOSE ZONE flow and transport studies re-u(h) 5u(h) 2 ur
u s 2 ur

[2]

USDA-ARS, George E. Brown, Jr. Salinity Lab., 450 W. Big Springs
Rd., Riverside, CA 92507. Received 11 Jan. 2005. *Corresponding
author (mschaap@ussl.ars.usda.gov).

Abbreviations: FMVG, Mualem–van Genuchten model with fitted
values for Ko and L; MMVG, modified Mualem–van Genuchten;
MVVG-A, modified Mualem–van Genuchten refitted with standard
retention parameters (hs 5 0); MVG, Mualem–van Genuchten;
TMVG, traditional Mualem–van Genuchten model with fixed values
for Ko and L.

Published in Vadose Zone Journal 5:27–34 (2006).
Original Researchdoi:10.2136/vzj2005.0005ª Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

27

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

28

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

showed that hs values as small as 21 or 22 cm produced
much more stable numerical simulations. They obtained
substantial improvements in the predicted hydraulic conductivity function of one clay soil when a value of 22 cm
was used. However, they did not investigate the optimalvalue for hs in terms of fitting errors of the retention or
the unsaturated hydraulic conductivity properties.
Another limitation of Eq. [2] involves the shape of thehydraulic conductivity function near saturation, especially of structured media (i.e., macroporous soils or unsaturated fractured rock), and how best to estimate thematching point, Ko. van Genuchten and Nielsen (1985)
and Luckner et al. (1989) pointed out that the matching
point, Ko, ideally should be located at a small negative
pressure head to avoid the effects of macropore flow thatcannot be captured with Eq. [2]. Schaap and Leij (2000)
and Schaap et al. (2001) confirmed that fixing Ko at Ks
leads to overpredicted hydraulic conductivity values atmost pressure heads. For a range of textures they foundthat the fitted Ko values were generally about one orderof magnitude smaller than the measured Ks values. In
addition, Schaap and Leij (2000) found that the fitted L
values were often negative, with an optimal value of 21.
Schaap et al. (2001) demonstrated that Eq. [2] with fitted Ko systematically underestimated hydraulic conductivity between 0 and 210 cm. However, as compared tothe traditional MVG model (with Ko 5 Ks and L 5 0.5),
their approach gave a much better description of the
conductivity at more negative pressures. This indicatesthat Eq. [2] needs to be modified further to describe theunsaturated hydraulic conductivity accurately in both
the wet and dry ranges. The main objective of this study
hence was to integrate the results found by Schaap andLeij (2000) and Schaap et al. (2001) with the model ofVogel et al. (2000). We investigate the optimal value of hs
for four broad soil textural groups in terms of the fitting
accuracy for both water retention and the unsaturated
hydraulic conductivity. While preserving the advantagesof the Vogel et al. (2000) model, the new model has
the potential to provide a much better description of theunsaturated hydraulic conductivity near saturation. The
performance of our modification (MMVG) will becompared to the traditional Mualem–van Genuchten
model with fixed (TMVG) and fitted (FMVG) values forKo and L.

THEORY
Modified Mualem–van Genuchten Model
Vogel et al. (2000) solved the conceptual andnumerical limitations of Eq. [1] and [2] by introducing
a small minimum capillary height, hs, and a fictitious(extrapolated) parameter um . us in Eq. [1]. Theirapproach maintains the physical meaning of us as ameasurable quantity, while the definition of effective
saturation Se is also not affected. The modified water
retention equation is given byuðhÞ 55

ur 1

u m 2 urn
(1 1 jahj )mus

h , hsh $ hs[4]

where
um 5 ur 1 (us 2 ur )(1 1 jahs jn)m[5]

Combination of Eq. [4] with Eq. [2] leads to themodified Mualem–van Genuchten (MMVG) model
12F(Se ) 2
LKo Se, h , hs
K(Se ) 5
12F(1)[6]Ko, h $ hs5

where
F(Se ) 5 [1 2 (Se )1=m ]m[7]

in which
Se (h) 5

u s 2 uru(h) 2 ur
Se (h) 5
u m 2 ur
u m 2 ur

[8]

This model reduces to Eq. [2] for hs 5 0. We refer to
Vogel et al. (2000) for a detailed discussion of the above
equations. They suggested that hs be about 21 or 22 cm,
but did not indicate whether other values would lead tobetter descriptions of observed retention and especiallyconductivity data. In this study we investigated whether
hs can be fixed to a single optimal value applicable todifferent soil textures. The advantage of having a fixed hs
is that it would avoid nonuniqueness problems by notadding another free hydraulic parameter to Ko, L, ur, us,a, and n, while keeping applications of Eq. [6] as simple
as possible. Vogel et al. (2000) showed that the above
modification affects the shape of the retention curve
only minimally relative to Eq. [1], but that the effects on
the unsaturated hydraulic conductivity can be very significant (their Fig. 2 and 3), especially for relatively finetextured soils when n approaches 1.0.

Macroporosity Modifications
Schaap and Leij (2000) showed that the fitted valuesof Ko in Eq. [2] were generally about one order ofmagnitude smaller than Ks, thus causing an apparent
discontinuity in the hydraulic conductivity when Ks
would be used for h 5 0 and Eq. [2] for h , 0 cm.
Conceptually, the discrepancy between Ks and Ko can be
explained by the presence of macropores or fracturesthat dominate the flow regime near saturation and micropores that control matrix flow. Matrix flow would thus
be active at all pressures, while macropore flow would
be dominant only near saturation and become negligible
at some relatively small negative pressure. This conceptualization suggests that Eq. [2] with fitted Ko and L
values should be used only for matrix flow. The challenge
to be addressed below is how best to account for theeffects of macroporosity by further modifying Eq. [6].
Previous attempts to incorporate the effects of macropore flow in the hydraulic conductivity generally focusedon composite, dual-porosity type functions in which
separate equations are used for the macropore andmicropore contributions to the hydraulic conductivity
function (e.g., Peters and Klavetter, 1988; Durner, 1994;
Mohanty et al., 1997; Köhne et al., 2001). Consistent withthis approach, we will use Eq. [6] to describe matrix flow,

29

www.vadosezonejournal.org

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

but adapt this equation for macropore flow nearsaturation. The modification we propose uses two
conductivity terms, Km(h) given by Eq. [6] for matrix
flow and a dimensionless pressure-dependent scale
factor A(h) to account for the increased conductivity
near saturation between a certain pressure head hm and
0. The conductivity function thus becomesK(h) 5 A(h)Km (h) h . hm
Km (h)
h # hm5[9]

Ideally, hm would be equal to hs because this would
simplify the resulting model. However, no compelling
physical reason exists for this since hs was introduced to
improve in an empirical way the numerical behavior of
Eq. [2] near saturation. We therefore will not assume apriori that hs and hm are equal.
In essence, A(h) introduces a pressure-dependentmatching point in Eq. [6]. In this study we assumed that
A(h) has the same general functional form for all
conductivity data sets, irrespective of soil texture. Thisavoids the necessity to parameterize A(h) for eachconductivity data set separately, which is difficult or
impossible given the scarcity of available near-saturated
conductivity data. The functional shape of A(h) will be
studied by first fitting Eq. [6] to conductivity data using
the first-order approximation that hm 5 hs. Because we
explicitly assume that Eq. [6] does not hold near saturation, it is imperative that this equation be fitted only
to pressure heads less than hs. Contrary to Eq. [2], pressure heads between hs and 0 cm are thus not used to fit
Eq. [6]. We will subsequently derive the shape of A(h)by studying the residuals between the fitted Km(h) andthe measured K(h) data. These residuals are expected to
exhibit considerable scatter since measured hydraulicconductivity data near saturation are subject to experimental difficulties, while different data sets will lead todifferent values of Ks, Ko, n, and L in Eq. [6]. To better
scrutinize the shape of A(h) we computed scaled
residuals according toR(h) 5log K(h) 2 log Km (h)log Ks 2 log Km (h)[10]This equation forces R(h) to be between 0 and 1 in theregion where Km(h) underestimates measured K(h)data, irrespective of the actual values for Km(h) and Ks.
A plot of R vs. h should indicate also whether or not hm
can be assumed equal to hs.
MATERIALS AND METHODS
The data for this study were taken from the UNSODA
database (Leij et al., 1996; Nemes et al., 2001) and are thesame as those used by Schaap and Leij (2000). We used 235
laboratory data sets that had at least six u–h pairs and at least
five K–h pairs. Samples with chaotic data or with limited
retention or conductivity ranges were omitted. Figure 1 provides the textural distribution of the samples and their classification in only four textural groups designated as Sands (100data sets), Loams (41), Silts (58), and Clays (36).
The parameters in Eq. [1], [2], [4], and [6] were fitted towater retention and unsaturated hydraulic conductivity data

Fig. 1. Textural distribution of samples used in this study and thedelineation of the four textural groups.

with the simplex or amoeba algorithm (Nelder and Mead,
1965; Press et al., 1988). The objective function for water
retention was taken as

O
NW

OW (p) 5(ui 2 ui 9)2

[11]

i51

where ui and ui9 are the measured and estimated watercontents, respectively, NW is the number of measured water
retention data points for each sample, and p is the parameter
vector {ur, us, a, n}. The parameter hs is not fitted but fixed forthe entire data set. For optimization of the unsaturated
hydraulic conductivity parameters we minimizedONK

OK (p) 5[log10 (Ki ) 2 log10 (Ki 9)]2

[12]

i51

where Ki and Ki9 are the measured and estimated hydraulicconductivity, respectively, Nk is the number of measured K(h)data points and p 5 {Ko, L}. Logarithmic values of Ki wereused in Eq. [12] to avoid a bias toward high conductivities in
the wet range. For Eq. [6] we used only conductivity data atpressure heads #hs cm, unless mentioned otherwise. Fitted Kovalues were constrained to be smaller than or equal to Ks. In
the case of Eq. [2], this produced 62 samples with Ko 5 Ks (outof 235), while for Eq. [6] the number was only 19. We believe
this constraint was hit more often for Eq. [2] because of its
shape problems near saturation. We did not carry out
simultaneous fits of the water retention and conductivity data.
The goodness of fit was quantified with the RMSE

sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
OW,K (p)RMSEW,K (p) 5
NW,K 2 np[13]

where NW,K is the number of water retention or hydraulicconductivity measurements and np is the number of parameters that were optimized. Results will be presented as averages
for each textural group as well as for the complete data set.
Because logarithmic values were used in the conductivity
optimizations, RMSEK results are dimensionless. One RMSEK
unit can be interpreted as a factor 10 error in K(h). To more

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

30

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

effectively compare errors for different hs values we computed
relative RMSE values by dividing RMSEW and RMSEK for
hs, 0 with the corresponding RMSE values for hs 5 0 (i.e.,the MVG model).
To study whether Eq. [2] and [6] underestimate oroverestimate observed hydraulic conductivities at particularpressure heads, we computed mean errors for 10 consecutive
pressure head ranges between 0, 23.2, 210, 232, 2100,
2320, 21000, 23200, 210 000, 232 000, and 2100 000 cm
according to

MEK 5

O

1 N
(log Ki 9 2 log Ki )
N i51

[14]

where N is the number of observations derived from multiple
samples in each pressure head range. As with RMSEK, MEK
values are dimensionless because logarithmic values of the
conductivities are used.

Parameter Optimization Procedure
The optimal value of hs for water retention and conductivity
was studied by comparing fitting errors of Eq. [4] and [6] for hs
values between 0 and 220 cm. The most optimal value of hs will
be used to compute R given by Eq. [10]. The results should yield
information about the functional shape of A(h) and whether or
not hm can be assumed equal to hs. We subsequently study
whether Eq. [9] (the MMVG model) will provide a better
description, in terms of RMSE and ME, than Eq. [2]. The latter

equation will be used in two versions, one with fixedparameters (i.e., Ko 5 Ks and L 5 0.5, after Mualem, 1976)
and one where Ko and L are optimized (using results from
Schaap and Leij, 2000). We also will assess whether the modified model can be simplified by directly using parametersfrom Eq. [1], rather than parameters from the more complicated Eq. [4].

RESULTS AND DISCUSSION
Optimal Value for hs
Figure 2A shows that relative RMSEW values for water
retention are almost constant for 25 , hs , 0 cm. This
indicates that the water retention curve within thispressure head range is macroscopically not affectedby hs. For hs less than 25 cm, the RMSEW values
increased for all soil textural classes, but especially for
the Sands. The differences between textural groupaverages of the fitted ur, us, and n parameters in Eq. [1]and [4] (Table 1) are small for an hs value of 24 cm (thevalue ultimately chosen). Fitted a values for Eq. [4] were
somewhat larger than those for Eq. [1]. This is probably
due to correlation between hs and a, since both act
as an air-entry pressure. The RMSEW values for Eq. [4](hs 5 24 cm) were virtually the same as those for Eq. [1],indicating that both equations fitted the retention data
equally well.

Fig. 2. Dependence of relative RMSE on values of hs for the four textural groups used in this study: (a) water retention, (b) hydraulic conductivity,and (c) the mean of water retention and conductivity.

31

www.vadosezonejournal.org

Table 1. Average retention parameters for each soil textural group as obtained with Eq. [2] and [4].
Eq. [4] (hs 5 24 cm)

Eq. [1]
Textureur

us

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

3

All
Sands†
Loams§
Silts¶
Clays#

cm cm
0.055
0.052
0.056
0.031
0.098

23

0.442
0.396
0.512
0.428
0.512

a

n

21

cm
0.02190.02630.04070.01200.0178

ur

RMSEW
3

1.64
2.23
1.19
1.38
1.30

0.01220.01220.01190.01410.0096

us
23

cm cm
0.051
0.051
0.046
0.039
0.088

a

n

RMSEW

1.63
2.23
1.18
1.38
1.29

cm3cm23
0.01240.01240.01210.01430.0095

21

0.439
0.393
0.506
0.426
0.509

cm
0.02450.02750.05370.01320.0204† Sand, Loamy Sand, Sandy Loam, and Sandy Clay Loam.
§ Loam and Clay Loam.
¶ Silty Loam and Silt.
# Clay, Sandy Clay, Silty Clay, and Silty Clay Loam.

The relative RMSEK values for the hydraulic conductivity (Fig. 2b) decreased sharply between hs 5 0 andhs 5 21 cm, especially for the Clays and Loams. The
reductions are caused by the exclusion of conductivities
with h . hs. This approach leads to more realistic Ko and
L values and produces better descriptions of the dry
parts of the hydraulic conductivity function. The improvement for the Sand group was not expected because
this group is normally associated with n values .2(Table 1). However, further inspection of the database
showed that 50 of the 100 samples of the Sand group hadn values,2. Thus, the Sand group also benefited fromthe improved description of the conductivity by using
Eq. [6] rather than Eq. [2]. For hs less than 21 cm, therelative RMSEK continued to decrease with lower hs
within the entire range between 0 and 220 cm (Fig. 2b).
The RMSEW and RMSEK graphs in Fig. 2a and 2bshow different optimal values for hs (i.e., minima in the
curves) for the water retention and conductivity data. To
arrive at a optimal value for hs we computed the averageof the relative RMSE values as (RMSEW 1 RMSEK)/2.
The optimum value of hs is now located between 23 and25 cm for the Sands, Loams, and Silts, while the Clayshave a relatively broad minimum with an optimum value
at approximately 28 cm (Fig. 2c). Since the data set as awhole had an optimum at 24 cm, we decided to use thisvalue for the remainder of our study. Use of Laplace’s
capillary law indicates that this pressure head of 24 cm
corresponds with a circular pore of about 0.04-cm radius.

Correction Near Saturation
Figure 3 shows scaled conductivity residuals, R, vs.
pressure head, h, for data from all 235 hydraulic conductivity data sets. Values for R are given as averages for1-cm pressure head intervals between 0 and 210 cm, 2cm intervals between 210 and 250 cm, and 5-cm
intervals for h less than or equal to 250 cm. The graphclearly shows that R decreases from near 1 at h 5 0 cm toapproximately 0 near h 5 240 cm. This indicates thatEq. [6] underestimates observed unsaturated hydraulicconductivities within this range, while below 240 cm nosignificant under- or overestimation was apparent. The
trend in these data suggests a sharp decrease in R nearsaturation and a slower decrease further away fromsaturation. A two-element piecewise linear function wasused to describe this pattern, with R 5 0 at h 5 240, achange in slope at h 5 24 cm, and R 5 1 at h 5 0 cm.
The change in slope was deliberately located at 24 cm tokeep the definition of Eq. [9] as simple as possible, andconsistent with the earlier derived value of 24 for hs.
With a least-squares minimization procedure we
determined that R at 24 cm should be 0.25 (i.e., thebreak in the fitted curve of Fig. 3), leading to thefollowing expression for R(h):0R(h) 5 0:2778 1 0:00694h
1 1 0:1875h5

h , 240 cm
240 # h , 24 cm [15]
4 # h # 0 cm

Fig. 3. Average scaled conductivities, R (Eq. [16]), for pressure heads between 0 and 2100 cm. The averages were calculated for pressure intervalsof 1 cm between 0 and 210 cm, 2 cm between 210 and 250 cm, and intervals of 25 cm between 250 and 2100 cm. The piece-wise linear curve
indicates the most optimal path through the data between 240 and 0 cm (Eq. [14]).

32

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

RMSEK values for this model (MMVG) and threeother models (TMVG, FMVG, and MVVG-A, modified
Mualem–van Genuchten refitted with standard retention parameters [hs 5 0] [MMVG-A], the latter to be
discussed further below) for all textural groups. We notethat the TMVG model has no degrees of freedom
whereas the FMVG and MMVG variants have two each(i.e., Ko and L). Fitted log(Ko) values for the MVVG
model were considerably lower than those found by
Schaap and Leij (2000) for the FMVG model (Eq. [2])assuming variable Ko and L. The smaller Ko values for
MMVG are due to the conceptual split between matrix
and macropore flow. Average fitted L values were
somewhat larger for the MMVG model than for the
FMVG formulation. However, the L values were stillnegative, which indicates that the new model does not
resolve problems with the physical interpretation of thepore interaction term SeL in Eq. [2] and [6] as discussed
previously by Kosugi (1999) and Schaap and Leij (2000).
The average RMSEK of the FMVG model was about 0.9
order of magnitude lower than the RMSEK of the
traditional TMVG model, which uses Ko 5 Ks and L 5
0.5 (Table 2). The RMSEK of the MMVG model is still
lower than for the FMVG model (0.261 vs. 0.410 for
FMVG), thus indicating that Eq. [16] provides a substantially better description of the unsaturated conductivity than Eq. [2].
We noted above that the RMSEW and the fitted
parameters values obtained with Eq. [1] and [4] differed
only marginally. This suggests that the MMVG model
can be used directly with retention parameters estimatedwith Eq. [1] (e.g., using the RETC code; van Genuchten
et al., 1991), rather than having to fit the more complicated Eq. [4] to observed retention data. Mean Ko and
L parameters for model MMVG-A are listed in Table 2.
This model uses Eq. [16] but assumes that all retention
parameters (i.e., ur, us, a, and n) were obtained with Eq.
[1]. The differences in averages of fitted parametervalues are very small, especially for logKo. Likewise,RMSEK values are only slightly larger for MMVG-A(0.263 vs. 0.261) for all 235 samples. This confirms thatEq. [16] can be used reliably with retention parameters
estimated with Eq. [1]. However, it is still necessary to
use Eq. [5] through [8] and Eq. [15] and [16] to computethe unsaturated hydraulic conductivity function.

Solving Eq. [10] for K(h) gives

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

K(h) 5[Ks
K m (h)

]

R(h)

[16]

Km (h)

where Ks is now the measured saturated hydraulicconductivity, R is calculated according to Eq. [15], andKm(h) is calculated using Eq. [6]. Note that thisapproach does not require an extra parameter, as did
the traditional TMVG model. However, in addition tounsaturated hydraulic conductivity data, a measured Ks
is now necessary to fit Ko and L with Eq. [16]. Although
K(h) is not constant between hs and 0 cm but increases
to account for macropore flow effects, Eq. [16] preserves the advantage over Eq. [2] in that it does nothave a very steep slope at h 5 0 cm. If Km(h) were
calculated according to Eq. [2], the numerical problems
found by Vogel et al. (2000) would most likely remain.
The need for a two-element piecewise linear correction in Fig. 3 is interesting because it may reflect thecontributions of three different sets of pores (0 to 24,
24 to 240, and less than 240 cm) to the hydraulic conductivity, rather than only macropores and micropores.
These three pore regions would comprise macropores,
mesopores, and micropores (e.g., Luxmoore, 1981), or
primary fractures, secondary fractures, and the soil orrock matrix (cf. Simunek et al., 2003). Having threeoverlapping pore regions would also give credence tothe use of multiporosity or multipermeability models
for preferential flow in structured media (e.g., Wilson
et al., 1992; Gwo et al., 1995; Hutson and Wagenet,1995), rather than only two regions as is common inmost dual-porosity and dual-permeability models (vanGenuchten and Sudicky, 1999; Simunek et al., 2003).
Still, we note that the deviations between 24 and240 cm in Fig. 3 could have been caused in part also
by mathematical or other limitations in Eq. [4] or [6]not specifically accounted for in our study.
Contrary to Eq. [6], which only holds for hydraulicconductivity data with h # hs, Eq. [16] was able to describe all conductivity data, including the near-saturated
data. We therefore decided to fit also the Ko and L
parameters in the Km term of Eq. [16] directly to all
available conductivity data for the entire range of available pressure heads. Table 2 lists fitted parameters andTable 2. Fitting results for different versions of the Mualem–van Genuchten (MVG) and modified Mualem–van Genuchten (MMVG)
models. The results for the traditional Mualem–van Genuchten model with fixed values for Ko and L (TMVG) and Mualem–van
Genuchten model with fitted values for Ko and L (FMVG) were derived from Schaap and Leij (2000). The MVVG model includes
correction near saturation, and the MVVG-A represents the same model but refitted with standard retention parameters (hs 5 0).
TMVG Eq. [1]log(Ks)All
Sands†
Loams§
Silts¶
Clays#

Log(Ko)cm d211.92
log(Ks)2.24log(Ks)2.03
log(Ks)1.70log(Ks)1.31
log(Ks)

L
0.5
0.5
0.5
0.5
0.5

MMVG Eq. [4] (hs 5 24 cm)

FMVG Eq. [1]RMSEK

log(Ko)

1.309
1.216
1.343
1.119
1.703

cm d21
1.03
1.29
1.42
0.82
0.26

† Sand, Loamy Sand, Sandy Loam, and Sandy Clay Loam.
§ Loam and Clay Loam.
¶ Silty Loam and Silt.
# Clay, Sandy Clay, Silty Clay, and Silty Clay Loam.

L
23.09
21.28
26.97
21.22
25.96

RMSEK

log(Ko)

0.410
0.395
0.398
0.403
0.481

cm d21
0.45
0.99
0.47
0.08
20.45

L
21.59
20.48
25.51
21.26
20.75

MMVG-A Eq. [1]RMSEK

log(Ko)

L

RMSEK

0.261
0.275
0.242
0.247
0.268

cm d21
0.45
0.98
0.49
0.08
20.50

21.36
20.40
24.50
21.10
20.85

0.263
0.276
0.254
0.247
0.265

33

www.vadosezonejournal.org

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Systematic Deviations vs. Pressure
Mean errors obtained with the TMVG, FMVG, and
MMVG models for nine pressure head classes are shownin Fig. 4. Negative ME values indicate underestimationof the conductivity, while the gray bars denote thenumber of samples in each pressure head class. Resultsfor MMVG-A are not shown since they were indistinguishable from those for MMVG. Notice that the
TMVG model overestimates conductivities between 0
and 2100 cm for the Sands, between 0 and 210 000 cm
for the Silts, and between 0 and 23000 cm for the Clays.
For the Loams, the TMVG model describes conductivity
well until a pressure of 2100 cm. Conductivities of theSands and Loams in particular are severely underestimated with the TMVG model at more negative
pressures. The FMVG model further underpredicts the
conductivity between saturation and 230 cm, thus confirming the findings shown in Fig. 3. Below 230 cm themean errors are nearly zero. The Sand and the Silts showa small underestimation at 210 000 cm with the FMVG
and MMVG models, which likely is a result of the small
number of samples in this pressure range (gray bars).
Because of the correction inherent in Eq. [15] and [16],the MMVG model as expected performed better thanthe FMVG model near saturation. A slight underestimation exists for the Sands and Silts at pressures at about
23 cm. For the Loams, the MMVG model overestimatedthe conductivity between 0 and 210 cm. For all samples, the MMVG model performed much better thaneither the TMVG model or the FMVG models, thus indicating that this model is able to describe unsaturated
conductivities correctly for the entire range of pressure
heads between 0 and 210 000 cm.

SUMMARY AND CONCLUSIONSA Mualem–van Genuchten model with improved
description of conductivity near saturation was presented. As previously proposed by Vogel et al. (1985,
2000), the model introduces a small but non-zero airentry pressure, hs, in the water retention curve. The
model was further modified to account for hydraulicconductivity matching points (Ko) that are much smaller
than the measured saturated conductivities (Ks). In thisstudy we found that the optimal value for hs was 24 cm.
We also found that the modified model still underpredicted K(h) between 0 and 240 cm. We were able
to correct for this underestimation by introducing apiecewise linear function that accounts for the effects ofmacropores and fractures near saturation. A small correction was needed between 24 and 240 cm (as causedby mesopores), while a stronger correction was necessary between 0 and 24 cm (as caused by macropores).
After refitting the Ko and L parameters, an average
RMSE of 0.261 remained for the database containing

Fig. 4. Mean errors (ME) for nine pressure head classes for the soil textural groups. The number of observations in each pressure class is indicatedwith the gray bars (right axis). Note that the scale on the right axis is not the same for all graphs.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

34

VADOSE ZONE J., VOL. 5, FEBRUARY 2006

235 samples. This error is lower than when the traditional Mualem-van Genuchten model was fitted (0.410)
with variably Ko and L (model FMVG), and much lower
than when default parameters were assumed for the
TMVG model (1.309). A plot of mean errors vs. pressure head showed that the MMVG has small systematicdeviations from measured hydraulic conductivities. The
modifications were found to be beneficial for all four
textural groups used in this study.
It is attractive to interpret the modified model interms of matrix and macropore flow, possibly augmented with flow through intermediate sized mesopores.
However, we emphasize here that the corrections werecompletely empirical and did not rely on statistical poresize distribution concepts such as those inherent in themodel of Mualem (1976). The modified model should be
used with some caution since hydraulic conductivity
measurements near saturation are generally quite unreliable. Other mathematical and conceptual limitationsin the Mualem–van Genuchten equations may well havefound their way into the derived conductivity corrections as well. Still, we believe that the new MMVG
formulation provides a substantial improvement in themacroscopic description of the hydraulic conductivity
function. The model seems to be especially suited forlarge-scale studies that require realistic simulations ofsaturated or near-saturated infiltration into soils. Sinceboth the retention and conductivity data could be
described well, the model may also be well suited forinverse modeling studies. However, because the modelin its present form does not explicitly define separate
matrix and macropore domains, it may not be well suited
for simulations of preferential flow.
ACKNOWLEDGMENTS
The authors gratefully acknowledge the support by NSF andNASA (EAR-9804902, EAR-0440024), and the ARO (39153–
EV). This paper is based in part also on work supported by
SAHRA (Sustainability of semi-Arid Hydrology and Riparian
Areas) under the STC Program of the National Science
Foundation, Agreement no. EAR-9876800.

REFERENCES
Butters, G.L., and P. Duchateau. 2002. Continuous flow method for
rapid measurement of soil hydraulic properties. I. Experimental
considerations. Available at www.vadosezonejournal.org. Vadose
Zone J. 1:239–251.
Durner, W. 1994. Hydraulic conductivity estimation for soils withheterogeneous pore structure. Water Resour. Res. 30:211–223.
Gwo, J.P., P.M. Jardine, G.V. Wilson, and G.T. Yeh. 1995. A multiplepore-region concept to modeling mass transfer in subsurface media.
J. Hydrol. (Amsterdam) 164:217–237.
Hopmans, J.W., J. Simunek, N. Romano, and W. Durner. 2002. Inverse
methods. p. 963–1008. In J.J. Dane and G.C. Topp (ed.) Methods ofsoil analysis. Part 4. SSSA Book Ser. 5. SSSA, Madison, WI.
Hutson, J.L., and R.J. Wagenet. 1995. A multiregion model describing
water flow and solute transport in heterogeneous soils. Soil Sci. Soc.
Am. J. 59:743–751.

Köhne, J.M., S. Köhne, and H.H. Gerke. 2001. Estimating the hydraulic functions of dual-permeability models from bulk soil data.
Water Resour. Res. 38(7)1121. doi:10.1029/2001WR000492.
Kosugi, K. 1999. General model for unsaturated hydraulic conductivity for soils with lognormal pore-size distribution. Soil Sci. Soc.
Am. J. 63:270–277.
Leij, F.J., W.J. Alves, M.Th. van Genuchten, and J.R. Williams. 1996.
The UNSODA Unsaturated Soil Hydraulic Database; User’s
manual, Version 1.0. EPA/600/R-96/095, National Risk Management Laboratory, Office of Research and Development, USEPA,
Cincinnati, OH..
Luckner, L., M.Th. van Genuchten, and D.R. Nielsen. 1989. A consistent set of parametric models for the two-phase flow of immiscible
fluids in the subsurface. Water Resour. Res. 25:2187–2193.
Luxmoore, R.J. 1981. Micro-, meso-, and macroporosity of soil. Soil
Sci. Soc. Am. J. 45:671–672.
Mohanty, B.P., R.S. Bowman, J.M.H. Hendrickx, and M.Th. van
Genuchten. 1997. New piecewise-continuous hydraulic functions
for modeling preferential flow in an intermittent flood-irrigatedfield. Water Resour. Res. 33:2049–2063.
Mualem, Y. 1976. A new model predicting the hydraulic conductivityof unsaturated porous media. Water Resour. Res. 12:513–522.
Nelder, J.A., and R. Mead. 1965. A simplex method for function
minimization. Comput. J. 7:308–313.
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. (Amsterdam) 251:151–162.
Peters, R.R., and E.A. Klavetter. 1988. A continuum model for water
movement in an unsaturated fractured rock mass. Water Resour.
Res. 24:416–430.
Press, W.H., B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling. 1988.
Numerical recipes in C. 1st ed. Cambridge Univ. Press, New York.
Schaap, M.G., and F.J. Leij. 2000. Improved prediction of unsaturated
hydraulic conductivity with the Mualem-van Genuchten model.
Soil Sci. Soc. Am. J. 64:843–851.
Schaap, M.G., F.J. Leij, and M.Th. Van Genuchten. 2001. Rosetta: A
computer program for estimating soil hydraulic parameters withhierarchical pedotransfer functions. J. Hydrol. (Amsterdam) 251:
163–176.
Simunek, J., N.J. Jarvis, M.Th. van Genuchten, and A. Gärdenäs. 2003.
Review and comparison of models for describing non-equilibrium
and preferential flow and transport in the vadose zone. J. Hydrol.
(Amsterdam) 272(1–4):14–35.
van Genuchten, M.Th. 1980. A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J.
44:892–898.
van Genuchten, M.Th., F.J. Leij, and S.R. Yates. 1991. The RETC codefor quantifying the hydraulic functions of unsaturated soils. Rep.
EPA/600/2–91/065. R.S. Kerr Environmental Research Laboratory,
USEPA, Ada, OK.
van Genuchten, M.Th., and D.R. Nielsen. 1985. On describing andpredicting the hydraulic properties of unsaturated soils. Ann.
Geophys. 3:615–628.
van Genuchten, M.Th., and E.A. Sudicky. 1999. Recent advances invadose zone flow and transport modeling. p. 155–193. In M.B.
Parlange and J.W. Hopmans (ed.) Vadose zone hydrology: Cutting
across disciplines. Oxford Univ. Press, New York.
Vogel, T., M. Cislerova, and M. Sir. 1985. Reliability of indirect estimation of soil hydraulic conductivity. Vodohosp Cas 33(2):204–224.
Vogel, T., M.Th. van Genuchten, and M. Cislerova. 2000. Effect of the
shape of the soil hydraulic functions near saturation on variablysaturated flow predictions. Adv. Water Resour. 24(2):133–144.
Wilson, G.V., P.M. Jardine, and J.P. Gwo. 1992. Modeling the hydraulicproperties of a multiregion soil. Soil Sci. Soc. Am. J. 56:1731–1737.
Yates, S.R., M.Th. van Genuchten, A.W. Warrick, and F.J. Leij. 1992.
Analysis of measured, predicted, and estimated hydraulic conductivity using the RETC computer program. Soil Sci. Soc. Am. J. 56:
347–354.