Geoderma 306 (2017) 144–151

Contents lists available at ScienceDirect

Geodermajournal homepage: www.elsevier.com/locate/geoderma

Fractal-based models for the unsaturated soil hydraulic functionsa,⁎aM.A. Alfaro Soto , H.K. Chang , M.Th. van Genuchten
a
b

MARK

b

Departamento de Geologia Aplicada and Centro de Estudos Ambientais, UNESP - Universidade Estadual Paulista, Brazil
Federal University of Rio de Janeiro, UFRJ, Rio de Janeiro, Brazil

A R T I C L E I N F O

A B S T R A C T

Handling Editor: Morgan Cristine L.S.

Reliable estimates of the unsaturated soil hydraulic properties are needed in many research and engineering
projects. Fractal based approaches have become popular for describing the hydraulic functions. New models
based on fractal geometry are presented to describe the wetting phase relative permeability (unsaturated hydraulic conductivity) of porous media. Analytical expressions were derived using the fractal capillary pressure(water retention) model proposed by Alfaro Soto and Vilar (2006) and the statistical pore size distribution
models of Mualem (1976a) and Burdine (1953) for the unsaturated hydraulic conductivity. The newly derived
models provide predictions of the relative permeability to water, as well as of the fractal dimension for modelingthe pore structure of a soil. Results are compared with previous formulations by Brooks and Corey (1964), vanGenuchten (1980) and Kosugi (1996), as well as with experimental water retention and hydraulic conductivity
data of four soils having widely diﬀerent soil textures. The proposed models agreed closely with observed data.
The analytical models may be attractive alternative for estimating the unsaturated hydraulic conductivity whenno laboratory or ﬁeld measurements are available.

Keywords:
Relative permeability
Fractal geometry
Water retention curve
Hydraulic conductivity

1. Introduction
Estimates of the unsaturated soil hydraulic conductivity are neededin various branches of the earth sciences, including civil and petroleum
engineering, hydrogeology, soil science, vadose zone hydrology, andenvironmental and geotechnical engineering. Although the unsaturated
hydraulic conductivity can be measured using a range of laboratory andﬁeld tests (e.g., Dane and Topp, 2002), the required experiments are
generally too expensive, complex, and time-consuming for routine applications, especially when the conductivity is very low. Moreover,results can vary considerably if such factors as the wetting and drainage
history of the porous medium or speciﬁc characteristics of the ﬂuid
(water, air, non-aqueous phase liquid) or medium are not considered in
the tests. These various limitations have encouraged the use of indirect
methods involving both empirical equations as well as more processbased formulations to determine the unsaturated hydraulic conductivity, K, as a function of the pressure head, h, or the water content,
θ (Brooks and Corey, 1964; Durner, 1994; Vereecken et al., 2010).
Statistical pore-size distribution models have made it possible toestimate the K(θ) or K(h) functions based on the characteristics of thesoil water retention (or capillary pressure - saturation) curve, θ(h). A
range of models can be used for this, as exempliﬁed by theoreticalformulations by Purcell (1949), Childs and Collis-George (1950),
Burdine (1953), Millington and Quirk (1961) and Mualem (1976a),⁎Corresponding author.
E-mail addresses: alfaro@rc.unesp.br (M.A. Alfaro Soto), chang@rc.unesp.br (H.K. Chang).

http://dx.doi.org/10.1016/j.geoderma.2017.07.019
Received 1 December 2016; Received in revised form 1 July 2017; Accepted 17 July 20170016-7061/ © 2017 Elsevier B.V. All rights reserved.

among others. Using these theories, various analytical models havebeen obtained for the hydraulic conductivity functions, such as the
widely used models of Brooks and Corey (1964) and van Genuchten(1980).
Fractal geometry has been used also for describing the hydraulicproperties of unsaturated porous media. Tyler and Wheatcraft (1990)
determined the hydraulic conductivity functions using a fractal modelfor the water retention curve, along with the statistical pore-size distribution models of Burdine (1953) and Mualem (1976a). Rieu andSposito (1991) combined fractal descriptions of the water retention
curve and the particle size distribution to obtain the fractal dimension
necessary for calculating the hydraulic conductivity. Other contributions were fractal formulations of the porous medium where used in
combination with predictive pore-size distribution approaches for thepermeability are given by Fuentes et al. (1996), Xu (2004) and Li(2010), among others.
In this paper we propose an alternative fractal-based description of
the unsaturated soil hydraulic functions using the conductivity models
of Burdine (1953) and Mualem (1976a). Below we ﬁrst brieﬂy summarize various formulations for the hydraulic properties based on bothtraditional (mostly empirical) and fractal-geometric descriptions of thehydraulic properties. These descriptions are then compared with the
new fractal-based formulations.

Geoderma 306 (2017) 144–151

M.A. Alfaro Soto et al.

Kr (Se ) = Se2 [1 − (1 − Se1 m)m] (m = 1 − 2 n)

2. Existing models for the soil hydraulic properties

The Mualem-based equation for Kr(Se) is obtained similarly by
substituting Eq. (7) into Eq. (2) and assuming m = 1 − 1/n on permissible m and n values, to give

2.1. Traditional analytical formulations
A large number of empirical equations and statistical pore size
distribution models have been used for the hydraulic properties of
unsaturated media, including functions proposed by Brooks and Corey(1964), Campbell (1974), van Genuchten (1980), Fredlund et al.
(1994), Kosugi (1996), and Assouline and Tartakovsky (2001). A detailed overview of various models is provided by Leij et al. (1997). Here
we brieﬂy summarize only those formulations that relate immediately
to the work discussed in this paper.
Most of the predictive equations for the hydraulic conductivity arebased on the statistical pore-size distribution models of Burdine (1953)
or Mualem (1976a), which are given by

Kr (Se ) = Se0.5 [1 − (1 − Se1 m)m]2 (m = 1 − 1 n)

∫ h2dτ(τ ) ∫ h2dτ(τ )0Se (h) =1
erfc {ln(h h 0) σ 2 }2(10)in which ho and σ are parameters characterizing the log-normal poresize distribution, and erfc is the complementary error function. Kosugi(1996) combined Eq. (10) with both Eqs. (1) and (2) to obtain

Kr (h) =
(1)

0

(9)

Another set of equations was derived by Kosugi (1996) who proposed a water retention curve assuming a log-normal pore radius distribution, to obtain

1

Se

Kr (Se ) = Se2

(8)

1 2
Se erfc [ln(h h 0) σ 2 + σ 2 ]2(11)for the Burdine-based conductivity Kr(h) function, and

and
Se

1

∫

∫

dτ
⎡
Kr (Se ) = Se0.5 ⎢
h
(τ )
⎣0

0

2

dτ ⎤
h (τ ) ⎥⎦Kr (h) =(2)

Tyler and Wheatcraft (1990) were among the ﬁrst to determine thepermeability of a porous medium to water using fractal mathematics.
For this they combined a fractal model for the water retention curvewith the statistical models of Burdine (1953) and Mualem (1976a).
Many others have followed since, including Toledo et al. (1990), Rieu
and Sposito (1991), Giménez et al. (1997), and Xu (2004), and various
references therein. Most of the fractal-based functions are relatively
simple power functions of the pressure head, which, if combined with
Burdine's conductivity model (Eq. (1)), will yield equations for threedimensional media of the form

in which S (=θ/θs) is relative saturation, Sr (=θr/θs) is residual saturation, and θ, θs and θr are the volumetric, saturated and residual
water contents, respectively.
Most early models for the hydraulic properties used a simple power
expression for the water retention curve, exempliﬁed by the classicalmodel of Brooks and Corey (1964), which can be combined with Burdine's model (Eq. (1)) to yield the following set of hydraulic functions(4)

Kr (Se ) = Se3 + 2 λ

(5)

(13)

11 − 3D(S ) = S 3 − D

(14)

e

e

in which D is the fractal dimension. Eqs. (13) and (14) are identical to
the Brooks and Corey equations given by Eqs. (4) and (5) if the poresize distribution exponent λ is redeﬁned in terms of the fractal dimension as
(15)

λ=3−D

(6)

which has become a key parameters in many fractal-based soil andhydraulic property studies.
Eq. (13) can also by combined with Mualem's conductivity model(Eq. (2)), in which case the fractal based conductivity function becomes

Limitations often associated with Eqs. (4), (5) and (6) are the sharp
breaks in the Se(h) and Kr(h) curves at the air entry value, the absence ofan inﬂection point, and possibly poor ﬁts with retention data in the dry
range. Some of these limitations motivated van Genuchten (1980) to
introduce a smooth (continuous diﬀerentiable) water retention curve of
the form

Se (h) = [1 + (αh)n]−m

Se (h) = (ha h)3 − DKrin which ha (generally referred to as the air entry value or bubblingpressure) and λ (also known to as the pore size distribution index) areessentially empirical parameters.
Eq. (4) can be combined also with Mualem's model (Eq. (2)) to give

Kr (Se ) = Se5 2 + 2 λ(12)2.2. Fractal formulations for the hydraulic functions(3)

Se (h) = (ha h) λ2 ]} 2for the Mualem-based Kr(h) expression.

respectively, where Kr = K(h)/Ks is the relatively hydraulic conductivity, Ks is the saturated hydraulic conductivity, h is the pressure
head (for notational convenience assumed here to be positive for unsaturated conditions), and Se is eﬀective saturation of the wetting
phase:Se (h) = (θ (h) − θr ) (θs − θr ) ≡ (S (h) − Sr ) (1 − Sr )1 0.5
Se {erfc [ln(h h 0) σ 2 + σ
219 − 3D

Kr (Se ) = Se 6 − 2D

(16)

Eqs. (13), (14) and (16), with various assumptions, modiﬁcations
and extensions, have been derived or used in many studies (Tyler andWheatcraft, 1990; Toledo et al., 1990; Fuentes et al., 1996; Giménez
et al., 1997; Ghanbarian-Alavijeh et al., 2010, Xu and Dong, 2004;
Ghanbarian-Alavijeh and Hunt, 2012b, among others), with someproducing diﬀerent power functions, especially of Eqs. (14) and (16),while others interpreted D alternatively as the fractal dimension of thepore radius, pore surface or particle size distribution. Also, some studies
assumed applicability of the equations over the entire moisture domain(thus assuming that θr = 0), while others applied the equations over
only a limited part (e.g., between the air entry value and the pressure
head at the permanent wilting point (15,000 cm or some other value)).
Similarly as the traditional Brooks and Corey Se(h) curve, Eq. (13)has limitations in terms of the sharp break of the curve at the air entry(7)

in which α, m and n, similarly as ha and λ in Eq. (4), are also empiricalparameters.
Eq. (7) with only minimal restrictions on permissible values of mand n (i.e., n > 1, and m > 0) can be combined directly with the
Burdine and Mualem predictive equations to obtain very general but
relatively complicated mathematical expressions for the hydraulicconductivity functions (van Genuchten and Nielsen, 1985; Dourado
Neto et al., 2011). With the restriction m = 2 − 1/n on permissible
values of m and n, the Burdine-based equations for the water retention
and hydraulic conductivity functions are obtained by substituting Eq.
(7) into Eq. (1) and integrating to yield (van Genuchten, 1980)
145

Geoderma 306 (2017) 144–151

M.A. Alfaro Soto et al.

value, ha, and the function not exhibiting a inﬂection point often observed in ﬁeld data. For this reason several have derived alternative
fractal-based expressions for the hydraulic properties (e.g., Cihan et al.,2009; Li, 2010). One is by Li (2010) who derived relative permeabilities
using fractal modeling of the medium in combination with the predictive conductivity models of Purcell (1949) and Burdine (1953). The
water retention curve used in his study was of the form (Li, 2010)h = h max (1 − bSe )−1 λS [1 + (h ho)−λ]h > haS (h) = ⎧ r⎨
1
0≤h ≤ ha⎩where S(h) = θ/θs and Sr = θr/θs.
Alfaro Soto and Vilar (2006) derived Eq. (21), similarly as Tyler andWheatcraft (1990), by using a Menger sponge algorithm for the fractal
pore size distribution without any assumption about the geometry ofmass and the solid-pore interface. The equation incorporates two portions of the water retention curve, one associated with initial desaturation of macropores and represented by the region where the fractal
characteristics are evident (with the fractal dimension being within in
the range 0 ≤ D ≤ 3), while a second region assumes that the lower
limit of the water content when h → ∞ is equal to the residual watercontent, θr.
We consider Eq. (21) to be a semi-empirical equation since themodel is made to partially satisfy a fractal condition in that the fractal
portion is invariant upon dilation and has fractal scaling as suggested
by Tyler and Wheatcraft (1990). In other words, the system remains
identical under diﬀerent magniﬁcations (Adler, 1992) and has a hierarchical structure with larger units containing successively smaller andsmaller units (Giménez et al., 1997).
Fig. 1 shows S(h) plots based on Eq. (21) using hypothetical valuesof 0.23 for Sr and 0.5 m for ho, with D varying from 1.0 to 2.9. The plots
reﬂect the fact that for given values of ho (in this case 0.5 m) and Sr, theair entry value ha decreases as the fractal dimension increases. Please
note also that all curves go through the point (ho, 2Sr) = (0.5, 0.46),
which follows immediately from Eq. (21). In all cases, relative saturation decreases between pressure heads of ho and ha, but this would
occur over a wider pressure head range if the residual water content isallowed to decrease simultaneously when D increases. The curves resemble to some extent typical Brooks and Corey retention curves, except that they do not exhibit straight lines on complete log-log plots ofSe(h), but only in the dry pressure head range.
Fig. 2 shows similar plots for diﬀerent values of residual saturation(Sr), varying from 0.1 to 0.5, while ho and D are maintained at 0.50 m
and 1.0, respectively. The plots show not only that, as expected, relative
saturation at a given pressure head increases with an increase in Sr, butalso with an increase in the air entry value, ha, as implied by Eq. (22).
The simultaneous increase in Sr and ha is typical of the water retention
curves of relatively ﬁne-textured soils (e.g., Carsel and Parrish, 1988;
Schaap et al., 2001).
The ﬂexibility of Eq. (21) in describing water retention data ofvarious soils is demonstrated in Fig. 3. The soils studied were Lamberg
clay, Caribou silty loam, Rubicon sandy loam, and Gravelly sand G.E.9,
with data sets obtained from the catalogue of Mualem (1976b). Results
obtained with Eq. (21) are compared in Fig. 3 with those using themodels of Brooks and Corey (1964), van Genuchten (1980) and Kosugi

(17)

where

b = 1 − (ha h max )−λ(18)in which, as before, λ = 3 − D, and where ha and hmax are the lower
and upper cutoﬀ values for Eq. (17), taken by Li (2010) as the air entry
value and the pressure head at the residual water content, respectively.
Eq. (17) was previously used also by Bird (1998), as well as by Perfect(1999) who assumed the residual water content to be zero. The inverse
expression of Eq. (18) is given by (Bird, 1998; Perfect, 1999;
Ghanbarian-Alavijeh and Hunt, 2012a):−λ−λSe = (h max− h−λ ) (h max− ha−λ ) ha ≤ h ≤ h max(19)

Using Eq. (17), Li (2010) derived a model for the relatively hydraulic conductivity of the wetting phases using the predictive conductivity model of Purcell (1949), leading to(2+λKr = 1 − (1 − bSe ) λ) (1 − β )
2+λ
λ

(20)
−λ

where b = 1 − β, and β = (ha/hmax) . One may verify that
Kr(Se = 0) = 0, and Kr(Se = 1) = 1 in the above equations. Li (2010)
showed for two media (Berea sand and a rock sample) that Purcell's
model worked very well for the wetting phase.
These various results suggest that fractal geometry can be used todescribe the water retention curve, as well as of the hydraulic conductivity to water. These issues are further explored in the next section
by using fractal theory to describe the water retention and hydraulicconductivity functions.
3. Alternative fractal model for the hydraulic functions
In this section, the fractal model for the water retention curve developed by Alfaro Soto and Vilar (2006) is used to derive formulations
for the unsaturated hydraulic conductivity function using the approaches of Mualem (1976a) and Burdine (1953).
3.1. Fractal model for the water retention curve
Alfaro Soto and Vilar (2006) expanded upon previous work byconsidering the pore-size distribution instead of the particle size distribution, and introducing an equivalent water retention curve which
does not assume a gradual and continuous reduction in pore size. Theirwater retention model is of the form

1.E+06
1.E+04Pressure head (m)θ [1 + (h ho)−λ]h > haθ (h) = ⎧ r⎨
1
0≤h ≤ ha⎩(21)

where λ = 3 − D as before, θr is the residual water content, ho is ansemi-empirical parameter, and ha is the usual air entry value characterizing a break in the retention function near saturation typical of
the Brooks and Corey (1964) formulation. This parameter is a functionof ho and λ as follows:ha = ho (Sr (1 − Sr ))1 λ(Sr (1 − Sr ))(h ho)−λh > ha10 ≤ h ≤ ha

1.E+02
1.E+00
1.E-02
1.E-04
1.E-060.0(22)

In terms of eﬀective (Se) and relative (S) saturations, Eq. (21) can be
written also in the forms

Se (h) = ⎧
⎨
⎩

(23b)

0.2

0.4

0.6

0.8

1.0

1.2

Water Saturation (-)
D=1.0

D=2.5

D=2.8

D=2.9

Fig. 1. Capillary pressure head – relative saturation curves calculated with Eq. (22) assuming Sr = 0.23 and ho = 0.50 m, with D varying from 1 to 2.9.

(23a)
146

Geoderma 306 (2017) 144–151

M.A. Alfaro Soto et al.

Pressure head (m)1.E+04(1964), van Genuchten (1980) and Kosugi (1996). However, the plot
for Lamberg clay in Fig. 3a revealed lower values for R2, mostly becauseof less that optimal ﬁts of the very wet and dry sides of the curve. The
ﬁgure also shows that extrapolation of the ﬁtted van GenuchtenMualem and Kosugi retention curves to the dry end did lead to lower
residual water contents than the Brooks-Corey and our curves.

1.E+02

1.E+003.2. Fractal models for the hydraulic conductivity
1.E-02

1.E-04

0.0

0.2

0.4

0.6

0.8

1.0

We next used Eq. (21) or (23a) to derive predictive equations for the
relative permeability (Kr) using the models of Burdine (1953) and
Mualem (1976a) as given by Eqs. (1) and (2), respectively. Solving Eq.
(23a) for h(Se) gives

1.2

1

Water Saturation (-)
Swr=0.10

Swr=0.25

Swr=0.35

h = ho [Se Sr ] D − 3
Swr=0.50(24)

Substituting Eq. (24) into Eq. (1) and carrying out the integrations
yields the relative permeability model based upon Burdine's approach:Fig. 2. Water pressure head – saturation curves calculated with Eq. (22) assuming
ho = 0.5 m and D = 1.0, with Sr varying between 0.10 and 0.50.

11 − 3D

Kr (Se ) = Se 3 − D(1996), while Table 1 presents the ﬁtted parameters obtained using
nonlinear regression. The fractal dimension varied between 0.70 and
2.61. No apparent relationship existed between texture and the value ofD. However, all D values were less than the theoretical upper limit 3.0(Hunt, 2004).
The plots in Fig. 3 and the R2 values in Table 1 indicate good
agreement of Eq. (21) with the data. They show that results with Eq.
(21) are very similar to those based on the models of Brooks and Corey(25)

Similarly, the relative permeability of the wetting phase can becalculated using the model of Mualem (1976a) by substituting Eq. (24)into Eq. (2) and integrating to yield
19 − 5D

Kr (Se ) = Se 6 − 2D(26)

The remarkable result here is that Eqs. (25) and (26) exactly duplicate the Brook-Corey fractal based Kr(Se) expressions given by Eqs.
(14) and (16). This even though diﬀerent soil water retention models

Fig. 3. Experimental capillary pressure head – volumetric water content data and ﬁtted curves obtained with Eq. (22), as well as with the models of Brooks and Corey (1964), vanGenuchten (1980), Kosugi (1996), Toledo et al. (1990) and Fuentes et al. (1996), for four soils: a) Lamberg clay, b) Caribou silty loam, c) Rubicon sandy loam and d) Gravelly sand.

147

Geoderma 306 (2017) 144–151

M.A. Alfaro Soto et al.

Table 1Fitted parameters of the water retention functions of Alfaro Soto and Vilar (ASV), Brooks
and Corey (BC), van Genuchten-Mualem, van Genuchten-Burdine and Kosugi, given by
Eq. (20), Eq. (7) with m = 1 − 1/n, Eq. (7) with m = 1 − 2/n, and Eq. (10), respectively.
Retention

Parameter

Unit

Function
1

ASV

BC

VGM

VGB

K

θs
θr
ho
D
R2
θs
θr
ha
λ
R2
θs
θr
α
N
m = 1 − 1/n
R2
θs
θr
α
n
m = 1 − 2/n
R2
θs
θr
ho
σ
R2

(–)
(–)
(cm)
(–)
(–)
(–)
(–)
(cm)
(–)
(–)
(–)
(–)
(cm− 1)
(–)
(–)
(–)
(–)
(–)
(cm− 1)
(–)
(–)
(–)
(–)
(–)
(cm)
(–)
(–)

Soil2
a

b

c

d

0.544
0.026
997
2.50
0.962
0.531
0.000
2.6
0.460
0.973
0.533
0.000
0.177
1.79
0.442
0.998
0.528
0.000
0.262
2.60
0.231
0.994
0.536
0.000
12.0
1.207
0.999

0.442
0.195
108
2.61
0.990
0.438
0.242
66.7
0.569
0.996
0.440
0.306
0.009
3.05
0.672
0.999
0.440
0.299
0.010
3.53
0.434
1.000
0.440
0.312
135.7
0.587
0.999

0.382
0.155
80
0.70
0.999
0.381
0.156
66.7
2.310
0.999
0.383
0.177
0.011
6.88
0.855
0.990
0.382
0.174
0.012
6.97
0.713
0.992
0.382
0.178
90.7
0.264
0.990

0.324
0.071
140
1.87
0.994
0.312
0.071
47.6
1.126
0.989
0.322
0.081
0.015
2.96
0.662
0.999
0.320
0.078
0.017
3.58
0.442
0.998
0.322
0.084
84.7
0.645
0.998

Fig. 5. Wetting phase relative permeability (Kr) - relative saturation (S) curves calculatedwith Eq. (27), based on Mualem's pore size distribution model assuming the same values
for ho (0.5 m) and D (1.0) as in Fig. 2, with Sr varying between 0.1 and 0.5.

parameters as in Figs. 1 and 2 respectively (i.e., ho = 0.5 m, Sr = 0.23
and D varying between 1.0 and 2.9 in Fig. 4 and D = 1 and Sr varying
between 0.1 and 0.50 in Fig. 5). The plots in Fig. 4 show a considerable
decrease in Kr as the fractal dimension (D) increases, leading to muchsteeper curves. For the theoretical upper limit of 3.0 for D, the curve
would become a step function, being 1 at full saturation, and zero
elsewhere. This situation is very similar to the traditional van Genuchten-Mualem and van Genuchten-Burdine equations, which havelower limits of the van Genuchten n values of 1 and 2, respectively (vanGenuchten and Nielsen, 1985; Vogel et al., 2001).
Fig. 5 shows that, as expected, an increase in residual saturation (Sr)forces the wetting relative permeability curves to be restricted to increasingly smaller parts of relative saturation close to saturation.

1
ASV: Alfaro Soto and Vilar (2006); BC: Brooks and Corey (1964); VGM: van
Genuchten (1980)-Mualem; VGB: van Genuchten (1980)-Burdine; K: Kosugi (1996).
2(a) Lamberg clay; (b) Caribou silty loam; (c) Rubicon sandy loam; (d) Gravelly sand.

are involved (Eqs. 13 and 23a). Because of the diﬀerent retention
equations, the Kr(h) functions will be diﬀerent also for the two models.
The fact that the same Kr(Se) expressions are obtained is a consequenceof the assumed scaling properties of the Burdine and Mualem theories(Eqs. (1) and (2), respectively), which are forcing the predictive conductivity curves to end at Kr(Se) = 1 when Se = 1. In other words, themultipliers of the Se(h) expressions (i.e., of Eqs. 13 and 23a) can be
diﬀerent, but not the exponents of Se in the two retention equations.
Figs. 4 and 5 shows plots of the relative hydraulic conductivity (Kr)versus relative saturation (S) calculated with Eq. (27), using the same

3.3. Example calculations of the relative hydraulic conductivity
Eqs. (25) and (26) were used next to obtain the relative permeabilityof the wetting phase using the Burdine and Mualem approaches, respectively. The parameters Sr, ho and D in the equations (Table 1) were
estimated from the water retention curves shown in Fig. 3 using nonlinear least-squares. The accuracy of the ﬁts were quantiﬁed by meansof the root-mean-square deviation (RMSD) between the modeled andmeasured values. Values of RMSD close to zero provide better predictions, with larger values indicating greater deviations from measured
values. Table 2 summarizes the results of the RMSD values obtained foreach soil using the following equation:
Table 2
Root Mean Square Deviation (RMSD) values of predicted relative hydraulic conductivity
functions using the pore size distribution models of Mualem and Burdine.
Mualem(1976b)catalogue
numberSoil name

a

1003B

b

3301

c

3501

d

4135

Lambergclay
Caribousilty loam
Rubiconsandy loam
Gravelly
sand G.E.9

Soil

Fig. 4. Wetting phase relative permeability (Kr) - relative saturation (S) curves calculatedwith Eq. (27), based on Mualem's pore size distribution model assuming the same values
for Sr (0.23) and ho (0.5 m) as in Fig. 1, and with D varying between 1.0 and 2.9.

1Conductivity
model1

M
B
M
B
M
B
M
B

Model2
ASV

BC

VG

K

0.40
0.31
0.12
0.13
0.25
0.40
0.09
0.10

0.46
0.37
0.16
0.18
0.26
0.41
0.12
0.09

0.36
0.29
0.88
0.66
1.26
1.29
0.07
0.13

0.41
0.32
0.87
0.96
1.27
1.41
0.08
0.15

M: Mualem (1976a), B: Burdine (1953).
ASV: Alfaro Soto and Vilar (2006); BC: Brooks and Corey (1964); VG: van Genuchten(1980); K: Kosugi (1996).
2

148

Geoderma 306 (2017) 144–151

M.A. Alfaro Soto et al.

Fig. 6. Comparisons of unsaturated hydraulic conductivity curves, K(θ), obtained with Eq. (26), Brooks and Corey (1964), van Genuchten (1980) and Kosugi (1996), all assuming
applicability of Burdine's (1953) model, and experimental data for four soils: a) Lamberg clay, b) Caribou silty loam, c) Rubicon sandy loam and d) Gravelly sand.

RMSD =[log(Km) − log(Kp )]2 (n − 1)

In view of the results in Figs. 3, 6 and 7, we believe that the proposed fractal formulation for the water retention curve in combination
with the Mualem or Burdine permeability models provide an possible
alternative for describing the unsaturated soil hydraulic properties.
Additional comparisons with measured data of diﬀerent soil are clearlyneeded. The accuracy of the K(θ) predictions may well be a function ofsoil texture, with some studies suggesting that Burdine's model isequally or more accurate than Mualem's model for very coarse-textured
soils, while the reverse may be true for medium and ﬁne-textured
systems (van Genuchten and Nielsen, 1985; Haverkamp and Parlange,
1986; API, 2007).
Although the models presented here (Eqs. (25) and (26)) provide
reasonably accurate predictions of θ(h) and K(θ), we note that thepredictions are based on a retention model that still assumes a constantsaturation until the air-entry pressure. As mentioned earlier, this situation may not be realistic for certain ﬁeld soils having relatively
broad pore-size or particle-size distributions (van Genuchten andNielsen, 1985; others). This restriction will also aﬀect predictions of thehydraulic conductivity if the description of the retention curve nearsaturation is imperfect.
Also, our results indicate that Eqs. (25) and (26) tend to underestimate the measured conductivity data in the dry range, although thisalso occurred with the other models discussed here. Some of the discrepancy may have been due to diﬃculty in describing the water retention curve on the basis of a single fractal dimension (D) applied tothe entire curve. Although Eq. (21) is not linear, the model is based on aconstant fractal dimension, which appears to be a limitation of manysoils showing multifractal characteristics in their pore or particle size

(27)

where Km and Kp represent, respectively, the measured and predicted
relative permeability and n represents the number of samples. The log(K) values were used since the conductivity data follow a log-normal
distribution (Kolmogorov-Smirnov statistical test using a signiﬁcancelevel of 5% for all comparisons).
Fig. 6 shows predicted curves of the hydraulic conductivity, K(θ),
for the same four soils as in Fig. 3. The plots compare the proposed
Burdine-based permeability model (Eq. (25)) with the Burdine-based
equations by Brooks and Corey (Eq. (5)), van Genuchten (Eq. (8)) and
Kosugi (Eq. (11)), as well as with the experimental data as documented
by Mualem (1976b) where the saturated hydraulic conductivity were:2.25 × 10− 5 cm/s, 1.66 × 10− 4 cm/s, 2.83 × 10− 4 cm/s and2.78 × 10− 4 cm/s in soils: a, b, c and d, respectively.
The results in Fig. 6 indicate close agreement between the K predictions with the fractal-based model and the experimental data, as wellas those obtained with the other Burdine based models. The fractal
model agreed very well with the experimental data of especially soils b,c, and d. RMSD values for the fractal-based model were lowest in mostcases as shown by the entries in Table 2. The RMSD values suggest that,
overall for our soils, Eq. (24) performed best, followed by the models ofBrooks and Corey, van Genuchten and Kosugi.
Similar results as in Fig. 6 are shown in Fig. 7 for the Mualem-based
K(θ) predictions obtained using Eq. (26). The fractal model again performed better in most cases as reﬂected by the RMSD values in Table 2.
The van Genuchten and Kosugi models for Rubicon sandy loam (soil c)in particular produced relatively poor results using either the Burdine
or Mualem based expressions in Figs. 6 and 7, respectively.
149

Geoderma 306 (2017) 144–151

M.A. Alfaro Soto et al.

Fig. 7. Comparisons of unsaturated hydraulic conductivity curves, K(θ), obtained with Eq. (26), Brooks and Corey (1964), van Genuchten (1980), Kosugi (1996) and Toledo et al. (1990),
all assuming applicability of Mualem's (1976a) model, and experimental data for four soils: a) Lamberg clay, b) Caribou silty loam, c) Rubicon sandy loam and d) Gravelly sand.

using our fractal-based functions, Eq. (7) with independent m and nvalues may well be more accurate than most or all approaches, including our formulation. While the resulting hydraulic conductivity
models in combination with Eq. (7) are deﬁnitely more complicated,
the recent introduction of advanced mathematical packages such as
Mathematica, Matlab or Mathcad with their tool boxes of special functions (in particular, the incomplete beta and hypergeometric functions)should soften this limitation. An even more far-reaching approach
would be to implement form-free descriptions of the hydraulic properties as shown by Iden and Durner (2007).

distributions. Several studies (e.g., Grout et al., 1998; Posadas et al.,
2001; Wang et al., 2008) have pointed out that multifractal concepts
are generally very well suited for the particle size distribution of manysoils, especially when ﬁne-textured. Since pore size and particle size areclosely correlated, the same may be true also for the water retention
curve.
Multifractal models may be especially attractive for tropical soils,
such as latosols and oxisols, which are more heterogeneous in their poresize distribution than many temperate soils. Measurement protocols
and other factors may play a role also as explained well in a recent
paper by Bieganowski et al. (2013) for the particle size distribution. As
multifractal approaches have proliferated for the particle-size distribution, similar philosophies may be beneﬁcial also for the soil hydraulic properties. Our approach in this paper still uses one fractal dimension, but implemented in a nonlinear manner over the entire
retention function. In view of the shortcomings of having an air-entry
value in the retention formulation, as done also in the classical Brooks
and Corey model, perhaps a multi-fractal approach may well prove to
be more promising.
Alternatively one could also implement more ﬂexible but purely
empirical retention functions, such as Eq. (7) as originally used by van
Genuchten (1980) with independent m and n values. As shown by van
Genuchten and Nielsen (1985), the presence of independent m and nvalues assumes that the curvatures of the curve near saturation at the
air entry value will be independent of the curvature in the dry side of
the curve. van Genuchten and Nielsen (1985) showed that this form of
the retention function is always more accurate than the Books andCorey and restricted van Genuchten-Mualem and van Genuchten-Burdine functions. While we did not further investigate this problem here

4. Conclusions
In this paper we presented and tested two new models based on
fractal geometry designed to predict the relative permeabilities of thewetting phases in two-phase systems for diﬀerent types of soils. The
models assume a more realistic fractal based description of the water
retention curve, in conjunction with the statistical pore-size distribution
models of Burdine (1953) and Mualem (1976a). The proposed models
are diﬀerent from those based on purely empirical formulations of thewater retention curve that do not necessarily reﬂect our understandingof the interaction between water and soil. The proposed formulations
showed close agreement with existing models for the wetting permeabilities, as well as with measured experimental data for four diﬀerent
soil types. An advantage of the models is their applicability over the fullrange of pore sizes in a soil. They appear useful tool for studying anddescribing the hydraulic properties of soils and rocks.

150

Geoderma 306 (2017) 144–151

M.A. Alfaro Soto et al.

Ghanbarian-Alavijeh, B., Hunt, A.H., 2012b. Unsaturated hydraulic conductivity inporous media: percolation theory. Geoderma 187–188, 77–84.
Ghanbarian-Alavijeh, B., Liaghat, A., Huang, G.-H., van Genuchten, M.Th., 2010.
Estimation of the van Genuchten soil water retention properties from soil texturaldata. Pedosphere 20 (4), 456–465.
Giménez, D., Perfect, E., Rawls, W.J., Pachepsky, Y., 1997. Fractal models for predicting
soil hydraulic properties: a review. Eng. Geol. 48 (1–3), 161–183.
Grout, H., Tarquis, A.M., Wisner, M.R., 1998. Multifractal analysis of particle size distributions in soil. Environ. Sci. Technol. 32, 1176–1182.
Haverkamp, R., Parlange, J.-Y., 1986. Predicting the water retention curve from particle
size distribution: I. Sandy soils without organic matter. Soil Sci. 142, 325–339.
Hunt, A.G., 2004. Comparing van Genuchten and percolation theoretical formulations of
the hydraulic properties of unsaturated media. Vadose Zone J. 3, 1483–1488.
Iden, S.C., Durner, W., 2007. Free-form estimation of the unsaturated soil hydraulicproperties by inverse modeling using global optimization. Water Resour. Res. 43,
W07451. http://dx.doi.org/10.1029/2006WR005845.
Kosugi, K., 1996. Lognormal distribution model for unsaturated soil hydraulic properties.
Water Resour. Res. 32 (9), 2697–2703.
Leij, F.J., Russell, W.B., Lesch, S.M., 1997. Closed-form expressions for water retention
and conductivity data. J. Hydrol. 35 (5), 848–858.
Li, K., 2010. More general capillary pressure and relative permeability models fromfractal geometry. J. Contam. Hydrol. 111 (1–4), 13–24.
Millington, R.J., Quirk, J.P., 1961. Permeability of porous solids. Trans. Faraday Soc. 57,
1200–1206.
Mualem, Y., 1976a. A new model for predicting the hydraulic conductivity of unsaturated
porous media. Water Resour. Res. 12, 593–622.
Mualem, Y., 1976b. A catalogue of the hydraulic properties of unsaturated soils. In:
Research Project Report No 442. Haifa, Technion, Israel Institute of Technology.
Perfect, E., 1999. Estimating soil mass fractal dimensions from water retention curves.
Geoderma 88, 221–231.
Posadas, A.D., Gimenez, P.D., Bittelli, M., Vaz, C.M.P., Flury, M., 2001. Multifractal
characterization of soil particle-size distributions. Soil Sci. Soc. Am. J. 65,
1361–1367.
Purcell, W.R., 1949. Capillary pressures, their measurement using mercury and the calculation of permeability therefrom. J. Pet. Technol. 1, 39–46.
Rieu, M., Sposito, G., 1991. Fractal fragmentation, soil porosity and soil water properties.
I. Theory. Soil Sci. Soc. Am. J. 55, 1231–1244.
Schaap, M.G., Leij, F.J., van Genuchten, M.Th., 2001. Rosetta: a computer program forestimating soil hydraulic parameters with hierarchical pedotransfer functions. J.
Hydrol. 251, 163–176.
Toledo, P.G., Novy, R.A., Davis, H.T., Scriven, L.E., 1990. Hydraulic conductivity ofporous media at low water content. Soil Sci. Soc. Am. J. 54, 673–679.
Tyler, S.W., Wheatcraft, S.W., 1990. Fractal processes in soil water retention. Water
Resour. Res. 26 (5), 1047–1054.
Vereecken, H., Weynants, M., Javaux, M., Pachepsky, Y., Schaap, M.G., van Genuchten,
M.Th., 2010. Using pedotransfer functions to estimate the van Genuchten–Mualem
soil hydraulic properties: a review. Vadose Zone J. 9, 795–820.
Vogel, T., van Genuchten, M.Th., Cislerova, M., 2001. Eﬀect of the shape of the soil
hydraulic functions near saturation on variably-saturated ﬂow predictions. Adv.
Water Resour. 24 (2), 133–144.
Wang, D., Fu, B., Zhao, W., Hu, H., Wang, Y., 2008. Multifractal characteristics of soil
particle size distribution under diﬀerent land-use types on the Loess Plateau, China.
Catena 72, 29–36.
Xu, Y., 2004. Calculation of unsaturated hydraulic conductivity using a fractal model forthe pore-size distribution. Comput. Geotech. 31, 549–557.
Xu, Y., Dong, P., 2004. Fractal approach to hydraulic properties in unsaturated porous
media. Chaos, Solitons Fractals 19, 327–337.

Acknowledgements
This work was supported by the Fundação para o Desenvolvimento
da Unesp – FUNDUNESP, São Paulo, Brazil.
References
Adler, P.M., 1992. Porous Media: Geometry and Transports. Butterworth and Heinemann,
Stoneham, MA.
Alfaro Soto, M.A., Vilar, O.M., 2006. Evaluation of a pore fractal model for the predictionof soil water retention curve. In: Miller, G.A., Zapata, C.A., Houston, S.L., Fredlund,
D.G. (Eds.), Unsaturated Soils 2006. ASCE Geotechnical Special Publ. 147 2. pp.
2441–2452. http://dx.doi.org/10.1061/40802(189)208.
API, 2007. Light Non-aqueous Phase Liquid Distribution and Recovery Model (LDRM).
Publication 4760. V.1,2. American Petroleum Institute, Washington, DC.
Assouline, S., Tartakovsky, D.M., 2001. Unsaturated hydraulic conductivity function
based on a fragmentation process. Water Resour. Res. 37 (5), 1309–1312.
Bieganowski, A., Chojeckib, T., Sochan, A., Lamorski, K., 2013. Methodological aspects offractal dimension estimation on the basis of particle size distribution. Vadose Zone J.
12 (1). http://dx.doi.org/10.2136/vzj2012.0064.
Bird, N.R.A., 1998. Comment on “An improved fractal equation for the soil water retention curve” by E. Perfect et al. Water Resour. Res. 34, 929–930.
Brooks, R.H., Corey, A.T., 1964. Hydraulic properties of porous media. In: Hydrology
Paper 3. Colorado State Univ., Fort Collins, CO.
Burdine, N.T., 1953. Relative permeability calculation size distribution data. Trans. Am.
Inst. Min. Metall. Eng. 198, 71–78.
Campbell, G.S., 1974. A simple method for determining unsaturated conductivity fromsoil moisture retention data. Soil Sci. 117, 311–314.
Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water
retention characteristics. Water Resour. Res. 24 (5), 755–769.
Childs, E.C., Collis-George, N., 1950. The permeability of porous materials. Proc. R. Soc.,
Ser. A 201, 392–405.
Cihan, A., Perfect, E., Tyner, J.S., 2009. Predicting relative permeability from water retention: a direct approach based on fractal geometry. Water Resour. Res. 45,
W04404. http://dx.doi.org/10.1029/2008WR007038.
Dane, J.H., Topp, G.C., 2002. Methods of soil analysis, part 4, physical methods. In: SSSA
Book Series 5. Soil Science Soc. Am., Madison, WI (2002).
Dourado Neto, D., de Jong van Lier, Q., van Genuchten, M.Th., Reichardt, K., Metselaar,K., Nielsen, D.R., 2011. Alternative analytical expressions for the general van
Genuchten-Mualem and van Genuchten-Burdine hydraulic conductivity models.
Vadose Zone J. 10, 618–623. http://dx.doi.org/10.2136/vzj2009.0191.
Durner, W., 1994. Hydraulic conductivity estimation for soils with heterogeneous pore
structure. Water Resour. Res. 32 (9), 211–223.
Fredlund, D.G., Xing, A., Huang, S., 1994. Predicting the permeability function for unsaturated soils using the soil-water characteristic curve. Can. Geotech. J. 31,
533–546.
Fuentes, C., Vauclin, M., Parlange, J.-Y., Haverkamp, R., 1996. A note on the soil-water
conductivity of a fractal soil. Transp. Porous Media 23, 31–36.
van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898.
van Genuchten, M.Th., Nielsen, D.R., 1985. On describing and predicting the hydraulicproperties of unsaturated soils. Ann. Geophys. 3 (5), 615–628.
Ghanbarian-Alavijeh, B., Hunt, A.H., 2012a. Comments on “More general capillarypressure and relative permeability models from fractal geometry” by Kewen Li. J.
Contam. Hydrol. 140-141, 21–23.

151