Journal of Hydrology 588 (2020) 125041

Contents lists available at ScienceDirect

Journal of Hydrology
journal homepage: www.elsevier.com/locate/jhydrol

Research papers

Simple functions for describing soil water retention and the unsaturated
hydraulic conductivity from saturation to complete dryness

T

Rudiyantoa, , Budiman Minasnyb, , Ramisah M. Shaha, Budi I. Setiawanc, ,⁎Martinus Th. van Genuchtend,e,
⁎

⁎

⁎

a

Program of Crop Science, Faculty of Fisheries and Food Science, Universiti Malaysia Terengganu, 21030 Kuala Nerus, Terengganu, Malaysia
School of Life & Environmental Sciences, Sydney Institute of Agriculture, The University of Sydney, Sydney, NSW 2006, AustraliacDepartment of Civil and Environmental Engineering, IPB University, Bogor 16680, IndonesiadDepartment of Earth Sciences, Utrecht University, 3584CB Utrecht, NetherlandseCenter for Environmental Studies, CEA, São Paulo State University, UNESP, Rio Claro, SP 13506-900, Brazilb

ARTICLE INFOABSTRACTThis manuscript was handled by Corrado
Corradini, Editor-in-Chief, with the assistanceof Juan V. Giraldez, Associate Editor

Many applications require accurate soil hydraulic functions that can describe the dynamics of water in variably
saturated systems from full saturation to complete dryness. Current functions considering zero water content at
the dry end involve the summation of capillary and adsorption functions. This study proposes relatively simple
yet complete soil hydraulic functions from full saturation to oven dryness that can be incorporated into existingmodels. The proposed water retention model is constructed by scaling the saturation function of Fredlund andXing (1994) to describe capillary and adsorptive water retention as well as zero water content at oven dryness.
To describe the hydraulic conductivity due to capillary and film flow, the model of Wang et al. (2018) is used.
We further incorporated hysteresis in the water retention function of Wang et al. (2018). Testing on 241 hydraulic property data sets covering a wide range of soil textures from sands to clays, shows good agreement withlaboratory observations, especially for coarse-textured soils. The proposed model accounts for the following
processes: hysteretic capillary and adsorptive water retention, hysteretic hydraulic conductivity as a function ofpressure head, non-hysteretic hydraulic conductivity as a function of water content, the effects of entrapped air,and closed scanning curves. This paper successfully derives the parameters of the main wetting curve from themain drying curve. Using the proposed functions for soil water flow modelling should lead to improved dynamic
flow predictions, especially for drying soil conditions.

Keywords:
Soil hydraulic functions
Water retention
Unsaturated hydraulic conductivity
Hysteresis
Soil moisture
Dry porous media

1. Introduction
The unsaturated soil hydraulic properties involve the soil water
retention curve describing the relationship between the pressure head,h, and the volumetric water content, θ, and the hydraulic conductivity
relating the pressure head or water content, θ, with the hydraulicconductivity, K. The hydraulic properties control a range of vadosezone hydrologic processes such infiltration, redistribution, evaporation,root water uptake, and contaminant transport. Conceptual descriptions
and numerical simulations of variably-saturated water flow require arobust definition of the soil hydraulic functions (Durner and Flühler,
2005; Jury and Horton, 2004; Nielsen et al., 1986).
Various functional forms have been developed over the years to
describe the soil hydraulic properties, as summarized by Leij et al.,

1997; Assouline and Or, 2013; Kosugi et al., 2002, among others. Most
functions describe the soil water content in terms of effective saturation, which traditionally is defined as a scaled value between the saturated water content and some residual water content at the dry end of
the curve (e.g., Brooks and Corey, 1964; Kosugi, 1996; van Genuchten,1980). A major assumption of these retention functions is that the watercontent does not go below some residual value (generally non-zero, but
not always) when the medium dries out and the pressure head (h) approaches -∞ (e.g., Peters, 2013). However, water contents are known tobecome zero at oven dryness. Another feature is evident when the water
retention function is combined with a hydraulic conductivity model(e.g., Mualem, 1976), in which case the predicted K(θ) function is valid
only in the capillary region. As a result, these models are generally
unable to describe both water retention and the hydraulic conductivity

Corresponding author.
E-mail addresses: rudiyanto@umt.edu.my (Rudiyanto), budiman.minasny@sydney.edu.au (B. Minasny), budindra@apps.ipb.ac.id (B.I. Setiawan),rvangenuchten@hotmail.com (M.T. van Genuchten).
⁎https://doi.org/10.1016/j.jhydrol.2020.125041
Received 20 January 2020; Received in revised form 28 April 2020; Accepted 1 May 2020
Available online 08 May 20200022-1694/ © 2020 Elsevier B.V. All rights reserved.

Journal of Hydrology 588 (2020) 125041

Rudiyanto, et al.

Shape parameter [L-1]
Shape parameter α of the drying curve [L-1]
Shape parameter of α of the main wetting curve [L-1]
Relative saturation function [-]Relative saturation of the drying scanning curve for calculating K(h)Γkw*(h) Relative saturation of the wetting scanning curve for calculating K(h)θVolumetric water content [L3 L-3]
θ(h)Soil water retention curve [L3 L-3]
θd*(h)
Drying scanning curve
Wetting scanning curveθw*(h)
θr
Residual water content [L3 L-3]
Saturated water content [L3 L-3]θsdθsSaturated water content of the main drying curve [L3 L-3]
θ sw
Saturated water content of main wetting curve [L3 L-3]
FX
Fredlund and Xing (1994) model
FXWJD Modified FX model by Wang et al. (2018)
MDC
Main drying curve
MWC
Main wetting curve
PDI
Model by Peters (2013), Peters et al. (2013) and Iden andDurner (2014)
R2
Coefficient of determination
Root mean square error for fitted water retention curve
RMSEθ
RMSElog(K) Root mean square error for fitted log hydraulic conductivity curve
RMSS
Proposed soil hydraulic model in this study
SCE-UA Shuffled complex evolution global optimization method
SDC
Scanning drying curve
SWC
Scanning wetting curvevG
van Genuchten (1980) soil hydraulic model

Nomenclature

Α
αd
αw
Γ(h)
Γkd*(h)Correction factor for effective saturation [-]Cf(h)hPressure head [L]h0Pressure head at zero water content [L]haAir entry value [L](hΔwd, θΔwd)Reversal point of drying scanning curve in term watercontent(hΔdw, θΔdw)Reversal point of wetting scanning curve in term watercontent(hΔwd, SeΔwd) Reversal point of drying scanning curve in terms of effective saturation(hΔdw, SeΔdw) Reversal point of wetting scanning curve in terms ofeffective saturation
K(h)Unsaturated hydraulic conductivity function [L T−1]
Kc(h)Capillary conductivity function [L T−1]
Kf(h)Film conductivity function [L T−1]
Relative hydraulic conductivity [-]Kr(h)Ks
Saturated hydraulic conductivity [L T−1]
Capillary conductivity at saturation [L T−1]
Ksc
Ksf
Film conductivity at saturation [L T−1]
L
Tortuosity/connectivity parameter [-]m
Shape parameter [-]n
Shape parameter [-]Sc
Capillary saturation function [-]Adsorption saturation function [-]Sa
Se
Effective saturation [-]Sed*(h)Drying scanning curve in terms of effective saturation
Sew*(h) Wetting scanning curve in terms of effective saturation
Sekd*(h) Drying scanning curve for calculating K(h)Sekw*(h) Wetting scanning curve for calculating K(h)Correction factor [-]Xmover the complete moisture range, especially at very low water contents(Peters, 2013; Weber et al., 2019). Modelling soil from saturation todryness is essential as plant and ecosystem responds differently in dry
and wet climate conditions (Tiktak and Bouten, 1992). In addition,modelling water conditions in the dry region becomes more criticalwhen soils are subject to extreme drying events, certainly in arid andsemi-arid regions, with concomitant worsening conditions as drought
events are becoming more frequent recently (Samaniego et al., 2018).
Many different approaches have been used to describe soil water
retention at and near complete dryness when the water content goes tozero. One major approach modifies the water retention function as asummation of capillary and adsorption related functions where adsorption function decreases, often linearly on a log-scale of the pressure
head, toward zero water content at the dry end (Campbell & Shiozawa,
1992; Fayer & Simmons, 1995; Fredlund & Xing, 1994; Khlosi et al.,2006; Morel-Seytoux & Nimmo, 1999; Ross et al., 1991; Rossi &Nimmo, 1994). A superposition between capillary flow and adsorptionfunctions can then be used also to describe the hydraulic conductivity
functions (Peters, 2013; Rudiyanto et al., 2015). Closely related modelsof this type were developed also by Lebeau & Konrad (2010), Liao et al.
(2018), Wang et al. (2016), Weber et al. (2019), and Zhang (2011).
Another popular approach is a multiplication approach, where thesaturation index is multiplied by a correction factor to force this variable to go to zero at the dry end (Fredlund and Xing, 1994; Wang et al.,2016). To describe the hydraulic conductivity, Wang et al. (2018) useda similar approach by introducing a scaling parameter Se into the tortuosity term of the standard van Genuchten-Mualem (vG-M) formulation (van Genuchten, 1980). Compared to the first approach, thissecond approach is more straightforward since it only modifies existingmodels and does not introduce more parameters into existing soil hydraulic functions such as those by van Genuchten (1980). Wang et al.

(2018) found that their approach generally yielded equal or better
agreement with observations compared to the summation approach.
Another important characteristic of the water retention curve is soil
water hysteresis in terms of the relationship between pressure head andwater content not being unique for wetting and drying. Several empirical formulations have been introduced to describe hysteretic water
retention in the van Genuchten (1980) functions (e.g., Huang et al.,2005; Kool & Parker, 1987; Parker & Lenhard, 1987). More recently,
Rudiyanto et al. (2015) implemented the Parker and Lenhard (1987)
hysteretic model into the PDI model formulated by Iden and Durner(2014), Peters (2013), and Peters et al. (2014) model. Hysteresis in themodels of Rudiyanto et al. (2015); and Parker and Lenhard (1987) only
occurs for capillary water retention, but not for adsorptive water in thedry retention region. This is in contrast to several studies indicating thathysteresis also occurs in the adsorptive range (e.g., Arthur et al., 2015;
Davis et al., 2009; Globus and Neusypina, 2006; Ning and Morteza,
2015; Prunty & Bell, 2007).
The objective of this study is to propose a simple yet complete soil
hydraulic functions that describe the soil water characteristic from saturation to oven dryness. The proposed model involves the capillary
and adsorption regions as well as has a zero water content at ovendryness. The hydraulic conductivity function includes capillary and film
flow. The model also incorporates hysteresis in both capillary and adsorptive water retention, as well as entrapped air, while assuming
closed scanning curves. Finally, the performance of the proposed model
is compared with the existing formulation by Wang et al. (2018) and
the PDI model by Iden and Durner (2014).

2

Journal of Hydrology 588 (2020) 125041

Rudiyanto, et al.

2. Model development
2.1. Soil water retention
The soil water retention model is formulated in terms of the relative
saturation function, Γ, being a variable between 0 and 1 and defined interms of pressure head (h). The model for Γ according to van Genuchten(1980), abbreviated as vG here, is given by(h) = (1 + | h|n ) m(1a)while according to Fredlund and Xing (1994), abbreviated as FX, thefunction is(h) = (ln[exp(1) + | h|n ]) m

(1b)

where h is the pressure head, and α, n, and m are shape parameters. The
difference between the vG and FX models for Γ(h) is illustrated graphically in Fig. 1a. With increasing values of |h|, the FX model for Γ(h)still decreases gradually, while the vG model reaches a constant (zero)
value much earlier (Wang et al., 2018). The gradually decreasing final
part of the FX curve may be interpreted as being a characteristic of theadsorption region of water retention, reason why the FX formulation isused in this study instead of the vG model. However, the FX curve neverreached a value of zero at realistic h values in the dry side (generally
thought to be about -107 cm). To obtain a zero water content at ovendry conditions, effective saturation, Se may be defined as a rescaled
variable of Γ as proposed by (Iden and Durner, 2014):

Se (h) =(h )1(h 0 )(h 0 )

Fig. 2. Distribution of soils used in this study according to the USDA soil texture
triangle. The number in each texture class refers to the number of data set.

2.2. Unsaturated hydraulic conductivity
Following Wang et al. (2018), the relative hydraulic conductivity, Kr
is calculated using(2)

where h0 is the pressure head where the water content becomes zero.
Following Schneider and Goss (2012), the value of h0 in our study is set
at −6.3 × 106 cm. The use of Eq. (2) hence results in Se values between
zero and unity within the interval (0 to h0), as shown in Fig. 1a. Thus,
the soil water retention function is given by(h ) =s Se (h )Kr (h) = Se L 1

(1

)

1 1
m

1 2
n

(4)

where Γ, n, and m are from Eq. (1b), Se is from Eq. (2), and L is aparameter that accounts for the tortuosity and connectivity of the flow
system. As shown in Fig. 1b, with increasing log(h), the term [1-(1Γ)1/m]1-1/n]2 decreases gradually, while the SeL term decreases rapidly.
The product of these two terms results in two stepwise curves, whichmodel the hydraulic conductivity as a result of capillary and film flow,
respectively. The first part of the curve is similar to the Kr function of(3)

where θs is the saturated water content. This simple function, based on
Eqs. (1b) and (2), will be used to describe water retention in both thecapillary and adsorption regions from saturation to the dry end.

Fig. 1. Illustration of effective saturation and relative hydraulic conductivity curves generated with the proposed RMSS model (using the parameters
α = 0.022 cm−1, n = 3, m = 0.7, and L = 3.5) and the vG model (using the parameters α = 0.016 cm−1, n = 2.9, m = 0.655, and L = 0.5). Note that Kr for vGMualem uses Eq. (4) with Se = Γ from Eq. (1a).
3Journal of Hydrology 588 (2020) 125041

Rudiyanto, et al.

Fig. 3. Observed and fitted soil hydraulic properties: soil water retention (left) and hydraulic conductivity (right) curves for the RMSS, FXWJD and PDI models
covering 11 different soil textures.

the vG-Mualem model while the second part is larger than the vGMualem model (Mualem, 1976; van Genuchten, 1980). The unsaturated
hydraulic function, K(h) is now given by

K (h) = Ks Kr (h)where Ks is the saturated hydraulic conductivity.
4(5)

Journal of Hydrology 588 (2020) 125041

Rudiyanto, et al.

Fig. 3. (continued)

The difference between the proposed model (Eqns. (1b)-(5)) and themodel by Wang et al. (2018) is that the current model uses the same
effective saturation function of Se (Eq. (2)) for both water retention andthe hydraulic conductivity. The model of Wang et al. (2018) uses two

separate Se functions (as explained further in section 3.2). In addition,
the proposed model, here called the RMSS model, accounts for hysteresis, which will be described in section 2.3. Also, it is important tonote that we used here the assumption by Wang et al. (2018) that m and5

Journal of Hydrology 588 (2020) 125041

Rudiyanto, et al.

Fig. 3. (continued)

6

Journal of Hydrology 588 (2020) 125041

Rudiyanto, et al.

Fig. 3. (continued)

since θsd > θsw, the Se and Γ values for the unsaturated hydraulicconductivity of the MWC will be scaled using the ratio between θsw andθsd to obtain a non-hysteretic K(θ) function (see Eqs. (14) to (16)).

n in Eqs. (1a) and (4) are independent and not related to each otherthrough the original vG-Mualem restriction that m = 1–1/n. Theirapproach avoids the use of far more complicated functional forms for K(h), such as incomplete beta or hypergeometric functions (Dourado
Neto et al., 2011; van Genuchten and Nielsen, 1985).

2.3.2. Scanning curves
Drying (DSC) and wetting (WSC) scanning curves always start froma reversal point. DSC curves start at the reversal point, (hΔwd, θΔwd), fromwetting to drying and proceed to a drying reversal point at (hΔdw, θΔdw).
Their values are transformed into (hΔ, SeΔ) using the following equation
to obtain both (hΔwd, SeΔwd) and (hΔdw, SeΔdw):

2.3. Hysteresis
2.3.1. Main curves
The main hysteretic loop consists of a main drying curve (MDC)
starting from saturation and ending in the dry end, and a main wetting
curve (MWC) starting from the dry end and ending at the saturation
point. The parameters for the MDC are given by [θsd, αd, n, m], and theparameters for MWC are [θsw, αw, n, m] (Kool and Parker, 1987; Parker
and Lenhard, 1987; Rudiyanto et al., 2015). Note that the n and mparameters are identical for both curves. Due to entrapped air, the saturated water content of the MDC is generally larger than the saturated
water content of the MWC, i.e., θsd > θsw. The drying and wetting
curves for the unsaturated hydraulic conductivity use the same Ks and L
values, but their Se and Γ values are derived from the MDC and MWC
retention curves given by Eqs. (9) and (10), respectively. Note that

Se =

(6)

d
s

The drying scanning curve in terms of Se from the starting point at(hΔwd, SeΔwd) to the endpoint at (hΔdw, SeΔdw) is subsequently linearly interpolated as follows (Parker and Lenhard, 1987; Rudiyanto et al.,2015):Se d (h) =[Se wd[Se d (h wd )Se dw ][Se d (h)Se d (hdw )]Se d (hdw )] + Se Sedw

(7)

Se (hΔwd) and Sed(hΔdw) are calculated using Eq. (2) with h = hΔwd andd7Journal of Hydrology 588 (2020) 125041

Rudiyanto, et al.

Fig. 4. Boxplots of RMSE for three models: RMSS, FXWJD, and PDI models as a function of soil texture groups: (a) RMSE for water retention in terms of water content(θ) and (b) RMSE of log(K).

h = hΔdw, respectively. Note that the parameter set from the main drying
curve is used here and that the reversal point at (hΔwd, θΔwd) and the end
point at (hΔdw, θΔdw) should be known to predict the DSC. The scanning
curve is then re-transformed back into water content term usingd(h ) =

d d
s Se

(h)

d
k

[1

d (h )]0d (h )
0(10)

The wetting scanning curves are calculated in the same way. A WSC
starts at some reversal point, (hΔdw, θΔdw) from drying to wetting, and thenproceeds to the wetting endpoint, (hΔwd, θΔwd). Values are transformed
into (hΔ, SeΔ) form using the following equation to obtain both (hΔdw,
SeΔdw) and (hΔwd, SeΔwd).

(8)

For the unsaturated hydraulic conductivity, which needs values of
Se and Γ, Se is calculated using

Se kd (h) = Se d (h)(h) = Se d (h)

Se =

(9)

while Γ is calculated using Eq. (2):

w
s

(11)

The wetting scanning curve from the starting point at (hΔdw, SeΔdw) to8

Journal of Hydrology 588 (2020) 125041

Rudiyanto, et al.

Fig. 5. Distribution of L values of the (a) RMSS, (b) FXWJD and (c) PDI models for various soil textural classes. Black vertical lines show median values.

9Journal of Hydrology 588 (2020) 125041

Rudiyanto, et al.

Fig. 6. Distribution of L values of the RMSS and FXWJD models for data sets from all 11 soil textures. (a) L value ranging 0.1 to 10 and (b) 0.1 to 8. Black vertical lines
show median values.

the endpoint at (hΔwd, SeΔwd) is now interpolated using

Se w (h) =[Se dw[Se w (hdw )Se wd ]
[Se w (h)Se w (h wd )]Se w (h wd )] + Se wd

w (h ) = S we(h ) =w w
s Se

We used a large number of datasets to test the accuracy of the
proposed water retention and hydraulic conductivity functions. A totalof 221 datasets were taken from the UNSODA database (Nemes et al.,2001). We limited ourselves to laboratory datasets subject to drying andcontaining at least five pairs for both water retention and the hydraulicconductivity. Another 11 datasets were taken by Peters (2013),
Minasny and Field (2005), Mualem (1976), Pachepsky et al. (1984) andRudiyanto et al. (2015). Fig. 2 shows that all soil textures were represented, with some preponderance of the silt loam, loam, and sandy
loam groups.
The hysteretic model was tested using water retention data of

w

s

(14)

and Γ is calculated using
w
k

(h) =

w
w (h ) sds(16)

3.1. Soil hydraulic property data sets

For unsaturated hydraulic conductivity at wetting scanning curve,
the Se is calculated using the following formula:

Se kw (h) = Se w (h) sd

w (h )]0w (h )
0

3. Datasets and model comparisons

(13)

(h )

[1

Note that both Se and Γ are scaled using the ratio between θsw andθsd to account for the possible effects of entrapped air.

(12)

Sew(hΔdw) and Sew(hΔwd) are calculated using Eq. (2) with h = hΔdw andh = hΔwd, respectively. Note again that the parameters from the main
wetting curve are used here and that the reversal point at (hΔdw, θΔdw) andend point (hΔwd, θΔwd) should be known to predict WSC. The equation isthen re-transformed into water contents using
w

(h )

(15)

with
10

Journal of Hydrology 588 (2020) 125041

Rudiyanto, et al.

wetting scanning (WSC) curves were generated in term of Se (Eq. (17))
using the same formulations as for our RMSS (i.e., Eqs. (6, 7 and 8) andEqs. (11, 12 and 13), respectively).
Unsaturated hydraulic conductivities of the drying scanning curves,which need values of Se and Γ, Se was calculated using Eq. (9), while Γ
was obtained using Eq. (17), which can be rearranged as followsdk(h ) =Se d (h)Cf (h)(18)

Hydraulic conductivities of the wetting scanning curve were similarly
obtained using Se as given by Eq. (14), while Γ was calculated from Eq.
(15):
w
k

(h) =

Sew (h)Cf (h)(19)

Note that, similarly as our RMSS model, both Se and Γ are scaled using
the ratio between θsw and θsd to account for the possible effects of entrapped air. To differentiate with the original FXWJD model, we referred to the modified model with hysteresis at the extended FXWJD
model.
Fig. 7. Sensitivity analysis of the value of L of the RMSS and FXWJD models
using a total of 4160 K(h) points from 232 soil samples covering 11 soil textures.

3.2.3. PDI model
The PDI model describes soil water retention as a weighted sum of acapillary saturation function, Sc, and an adsorption saturation function,
Sa:

Tottori dune sand (Rudiyanto et al., 2015). The data included the main
drying (MDC) and wetting (MWC) curves, as well as scanning loops forboth drying and wetting. The unsaturated hydraulic conductivity forthe main drying retention curve was obtained using the evaporation
method as outlined by Schindler and Müller (2006). In addition, we
evaluated possible relationships between the drying and wetting parameters αd and αw used data from 24 soils covering various soil textures.
These data sets were taken from Rudiyanto et al. (2015), Vachaud andThony (1971), Mualem (1976a), Yang et al. (2004), Huang et al.
(2005), Nemes et al. (2001), Zhou (2013), and Viaene et al. (1994).

(h ) = ( s( )ln (1 + )
ln

h0
hr

(h 0 )(h 0 )(21)(h) = (1 + | h|n ) m(22)

Sc (h) =(h)
1

in whichwhere α and n are empirical shape parameters as before, and m = 1– 1/
n. For the adsorption term, Sa(h) is formulated as

3.2.1. FXWJD model
The proposed model was compared to two existing models that
apply to the full soil moisture range: the modified FX model, or FXWJD,
by Wang et al. (2018) and the PDI model by Peters (2013), Iden andDurner (2014) and Peters (2014). The FXWJD model defines effective
saturation for water retention, Se as:

Se (h) = 1(20)r Sa (h )where θs and θr are the saturated and residual water contents, respectively. The capillary term Sc(h) in Eq. (20) is defined as a scaling vGtype saturation function, Γ, given by

3.2. Model comparisons

h
1+ h
r

r ) Sc (h ) +Sa (h) =

Xm 1

h
ha
h
ln 1 + 0
hr

ln 1 +

forh

ha
(23)

1forh > hawith1(h) = Cf (h) (h)Xm = 1(17)

where Cf(h) is the correction factor, hr = -1500 cm, being the pressure
head corresponding to the residual water content, h0 is the same as in
Eq. (2), and Γ is the relative saturation function given by Eq. (1b). The
FXWJD water retention function is calculated using θ(h) = θs Se(h)where θs is the saturated water content and Se(h) is effective saturationaccording to Eq. (17). For the unsaturated hydraulic conductivity, theFXWJD model defines effective saturation as given by Eq. (2), which isthen used to the relative unsaturated hydraulic conductivity using Eq.
(4).

ln(2)

(

h

ln 1 + h 0

a

)

(24)

in which h0 is the pressure head at zero water content, and ha is an air
entry value approximated by the value of 1/α in vG’s retention model.
We note that since a smoothing factor in the PDI model does not affectsignificantly the shape of the soil water retention curve (Iden andDurner, 2014; Wang et al., 2018), this factor is not used here in theadsorption term.
The PDI hydraulic conductivity function is defined as the sum ofcapillary, Kc(h), and film, Kf(h), conductivities:(25)

K (h) = K c (h) + Kf (h)3.2.2. Hysteresis in the FXWJD model
The original FWXJD model did not account for hysteresis. We implemented our hysteresis concepts (section 2.3) into the main drying(MDC) and wetting (MWC) curves by using different values only for theparameters αd < αw, while the remaining parameters were assumed toremain the same (i.e., θsd ≥ θsw, n and m for water retention, and Ks and
L for the hydraulic conductivity). The drying scanning (DSC) and

The capillary conductivity function Kc(h) is calculated using

K c (h) = Ksc ScL 1

1

1
m

1

0m

1

m 2

(26)

where Ksc is capillary conductivities at saturation, Γ is given by Eq.
(22), Sc is calculated using Eq. (21), and L is the usual pore tortuosity/
11

Journal of Hydrology 588 (2020) 125041

Rudiyanto, et al.

Fig. 8. Hysteresis in the RMSS model. Observed and fitted main loops and predicted scanning curves of water retention at (a) logarithmic and (b) linear scales. The
main and scanning loops of the hydraulic conductivity versus (c) pressure head and (d) volumetric water content corresponding to the main and scanning loops ofwater retention in Fig. 7a.

connectivity parameter (generally assumed to be 0.5). The film conductivity function Kf(h) in Eq. (25) is given by

Kf (h) = Ksf

h0ha

The same procedure was followed for the PDI model by optimizing
three parameters [θs, α, n]. Subsequently, the hydraulic conductivity
data were fitted in log form by optimizing Ks and L for the RMSS andFXWJD models and Ksc, Ksf, and L for the PDI model, while otherparameters were fixed at the optimized water retention parameters
values. We hence now minimized RMSElog(k) given by:

1.5(1 Sa)(27)

where Ksf is the film conductivity at saturation, Sa is calculated using
Eq. (23) and ha = 1/α in van vG’s retention model. Table A1 of the
Appendix summarizes the above three models, as well as the traditionalvG formulation.

M

RMSElog(k ) =First, the RMSS and FXWJD equations were fitted to the water retention data by optimizing four parameters [θs, α, n, m] in terms ofminimizing the RMSEθ as follows
N1( data, i
N i=1

predicted, i )log(Kpredicted, i ))2(29)

Hysteresis was evaluated in the same manner using the MDC and
MWC datasets in their entirety by fitting simultaneously now 6 parameters [θsd, αd, n, m, θsw, αw] by minimizing RMSEθ.
The shuffled complex evolution (SCE-UA) global optimization
method (Duan et al., 1994) in the Hydromad R package (Duan et al.,1994) was used to optimize the parameters following the above procedures to minimize the defined objective functions. Constraints wereset for the parameters n > 1.1, m > 0.1, and L < 10. For samples
having no retention measurements near saturation, the saturated water

3.3. Parameter optimization

RMSE =

1(log(K data, i)
M i=1

2

(28)
12

Journal of Hydrology 588 (2020) 125041

Rudiyanto, et al.

Fig. 9. Hysteresis in the extended FXWJD model. Observed and fitted the main loop and predicted scanning curves of water retention at (a) logarithmic and (b) linear
scales. The main and scanning loops of the hydraulic conductivity versus (c) pressure head and (d) volumetric water content corresponding to the main and scanningloops of water retention in Fig. 8a.

content was fixed at the measured value.

broader particle size distribution and fine-textured soils, the RMSS
model had a larger RMSEθ compared with the two other models. This
discrepancy is mainly due to the defined shape of the RMSS function, as
shown in Fig. 3d-i and k. At pressure heads beyond the air entry region,
the shape of the RMSS model is concave, while the shape of the curvefor the finer-textured soils is linear with increasing values of log(|h|).
The other two models were better in describing these types of soils thanthe proposed RMSS model since they have two extra terms that makethem more flexible. The FXWJD model involves the product between asaturation function and a correction function, while the PDI model is asummation of capillary and adsorption functions (See Table A1 in Appendix).
The boxplot of RMSE for log(K) shows that our RMSS formulation
performed similarly as the FXWJD model for the hydraulic conductivityof all textural classes. Close agreement between the observed and fitted
data is apparent in both the capillary and film regions (see Fig. 3). In
contrast, the PDI model showed close agreement only for the coarsetextured soils with their distinct capillary and film contributions. For

4. Results and discussion
4.1. Model comparisonsFig. 3 compares the optimized water retention and hydraulic conductivity functions with the observed data of 11 representative samples.
Results for all 232 soils are available in the supplementary materials(Figs. S1–S11 and Tables S1–S3). Distributions of the model error are
summarized in Fig. 4 in the form of RMSE boxplots for θ and log(K).
The RMSEθ boxplot for the water content as a function of soil textural class shows that the proposed RMSS model performed quite similar as the two existing models (FXWJD and PDI models) in describing
water retention of coarse-textured soils (the sand, loamy sand, sandy
loam, and silt loam textures). All three models describe the entire rangeof water retention from saturation to dryness and capture both thecapillary and adsorption regions very well. However, for soils with a13

Journal of Hydrology 588 (2020) 125041

Rudiyanto, et al.

Table 1
Optimized main drying and wetting parameters of the RMSS model obtained from hysteretic water retention data of 24 soils.
Soils

Sources of dataθsd
(cm3 cm−3)αd(cm−1)

n
(-)

m
(-)

θsw
(cm3 cm−3)αw(cm−1)

Fitted Main Drying Curve
RMSE
R2

Fitted Main Wetting Curve
RMSE
R2

Tottori dune sand
Sand
Avondale clay loam
Gilham dune sandGuelph loam
White silica sandIda silt loam 15 cm
Ida silt loam > 15 cm
Fine sandRubicon sandy loam
Fine sandsand 4890
clay loam 1301
loam 4710
sand 2310
sandy clay loam 1300
silty clay loam 2020
silty clay loam 2021
silty clay loam 2022
silty clay loam 4700
sand 1410
silt loam 4870
Houston sandCoarse sand(Rudiyanto et al., 2015)(Vachaud and Thony, 1971)(Watson et al., 1975)(Mualem, 1976b)(Mualem, 1976b)(Huang et al., 2005)(Mualem, 1976b)(Mualem, 1976b)(Maqsoud et al., 2012)(Mualem, 1976b)(Yang et al., 2004)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Zhou, 2013)(Viaene et al., 1994)

0.375
0.363
0.365
0.303
0.527
0.354
0.548
0.531
0.271
0.383
0.408
0.303
0.367
0.406
0.348
0.386
0.664
0.640
0.670
0.396
0.359
0.487
1.000
0.259

0.040
0.021
0.010
0.037
0.021
0.013
0.010
0.015
0.051
0.019
0.048
0.009
0.005
0.031
0.028
0.017
0.009
0.004
0.009
0.019
0.047
0.010
0.065
0.030

8.855
10.0001.100
10.000
1.720
6.597
1.100
1.179
7.760
3.946
5.723
6.480
1.100
1.701
10.000
1.100
1.100
1.103
1.100
1.241
10.000
1.104
6.517
4.701

0.801
0.651
0.141
0.336
0.100
0.397
0.370
0.454
0.915
0.100
1.000
0.675
0.655
0.890
1.000
0.439
0.100
0.100
0.100
1.000
0.769
0.100
1.000
1.000

0.338
0.354
0.365
0.294
0.433
0.345
0.499
0.492
0.216
0.363
0.397
0.299
0.360
0.385
0.299
0.388
0.653
0.640
0.670
0.399
0.355
0.500
1.000
0.250

0.064
0.044
0.015
0.064
0.033
0.034
0.013
0.018
0.079
0.041
0.102
0.018
0.017
0.032
0.054
0.032
0.010
0.010
0.010
0.028
0.072
0.053
0.142
0.050
Mean
SD

0.014
0.011
0.009
0.012
0.010
0.014
0.007
0.007
0.007
0.023
0.010
0.005
0.008
0.006
0.018
0.010
0.015
0.008
0.018
0.013
0.012
0.018
0.044
0.012
0.013
0.008

0.013
0.015
0.015
0.009
0.033
0.013
0.013
0.008
0.012
0.022
0.021
0.007
0.022
0.007
0.018
0.008
0.018
0.036
0.031
0.005
0.012
0.007
0.054
0.013
0.017
0.012

0.995
0.990
0.943
0.988
0.991
0.992
0.996
0.997
0.988
0.953
0.998
0.997
0.941
0.995
0.995
0.924
0.894
0.894
0.871
0.848
0.994
0.979
0.989
0.990
0.964
0.045

0.991
0.982
0.901
0.996
0.913
0.993
0.987
0.996
0.984
0.966
0.990
0.996
0.560
0.993
0.996
0.876
0.870
0.792
0.654
0.988
0.995
0.999
0.980
0.984
0.933
0.114

Table 2
Hydraulic RMSS model parameters (αw and θsw = θsd) of the main wetting curve as predicted from the optimized main drying curve of 24 soils.
Soils

Tottori dune sand
Sand
Avondale clay loam
Gilham dune sandGuelph loam
White silica sandIda silt loam 15 cm
Ida silt loam > 15 cm
Fine sandRubicon sandy loam
Fine sandsand 4890
clay loam 1301
loam 4710
sand 2310
sandy clay loam 1300
silty clay loam 2020
silty clay loam 2021
silty clay loam 2022
silty clay loam 4700
sand 1410
silt loam 4870
Houston sandCoarse sandSources of data

(Rudiyanto et al., 2015)(Vachaud and Thony,
1971)(Watson et al., 1975)(Mualem, 1976b)(Mualem, 1976b)(Huang et al., 2005)(Mualem, 1976b)(Mualem, 1976b)(Maqsoud et al., 2012)(Mualem, 1976b)(Yang et al., 2004)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Zhou, 2013)(Viaene et al., 1994)

θsd

αd

n

m

θsw

αw

Fitted Main Drying Curve

(cm3
cm−3)

(cm−1)

(-)

(-)

(cm3
cm−3)

(cm−1)

RMSE

R2

RMSE

R2

0.374
0.363

0.037
0.022

9.614
10.000

1.000
0.553

0.374
0.363

0.067
0.039

0.005
0.009

0.998
0.991

0.030
0.062

0.988
0.937

0.362
0.305
0.525
0.357
0.542
0.522
0.266
0.382
0.410
0.303
0.368
0.360
0.350
0.383
0.662
0.638
0.653
0.390
0.359
0.489

0.012
0.034
0.020
0.011
0.010
0.017
0.051
0.015
0.050
0.009
0.005
0.028
0.028
0.013
0.007
0.004
0.007
0.028
0.045
0.011
0.062
0.029

1.666
10.000
1.841
5.802
1.447
1.767
10.000
9.658
5.195
6.001
1.100
3.056
10.000
1.101
1.100
1.100
1.100
1.181
10.000
1.100
9.301
7.512

0.100
0.543
0.100
0.684
0.231
0.195
0.704
0.100
0.974
0.852
0.634
0.462
1.000
0.624
0.259
0.100
0.100
0.572
0.959
0.100
1.000
0.903

0.362
0.305
0.525
0.357
0.542
0.522
0.266
0.382
0.410
0.303
0.368
0.360
0.350
0.383
0.662
0.638
0.653
0.390
0.359
0.489
0.000
0.256

0.021
0.062
0.036
0.020
0.018
0.030
0.093
0.026
0.090
0.016
0.009
0.050
0.050
0.024
0.013
0.007
0.013
0.051
0.081
0.019
0.112
0.052
Mean
SD

0.006
0.004
0.009
0.003
0.003
0.003
0.006
0.007
0.009
0.004
0.008
0.004
0.018
0.010
0.015
0.008
0.015
0.012
0.006
0.018
0.021
0.001
0.009
0.005

0.959
0.998
0.992
0.999
0.999
0.999
0.990
0.994
0.997
0.998
0.942
0.997
0.995
0.920
0.887
0.893
0.857
0.824
0.997
0.980
0.997
1.000
0.967
0.052

0.029
0.029
0.075
0.055
0.037
0.027
0.020
0.045
0.031
0.016
0.032
0.043
0.046
0.008
0.036
0.023
0.037
0.018
0.051
0.038
0.087
0.025
0.037
0.019

0.929
0.995
0.930
0.875
0.981
0.986
0.981
0.903
0.995
0.982
0.509
0.974
0.971
0.882
0.890
0.770
0.670
0.972
0.956
0.986
0.958
0.954
0.916
0.115

0.256

soils with a broader particle size distribution (generally n < 1.3), the
PDI model showed a worse fit due to poor performance of the Mualem
model for these soils (Iden et al., 2015; Vogel et al., 2000; Wang et al.,2018), as well as unclear separation between capillary and film regions.

Predicted Main Wetting Curve

4.2. Parameter L
The parameter L can have an important effect on predictions of theunsaturated hydraulic conductivity function using Mualem’s model(Mualem, 1976) since it affects the bend or slope of the curve, especially at low water contents. L is frequently interpreted as a pore tortuosity and/or connectivity factor (Mualem and Dagan, 1978). Based on
14

Journal of Hydrology 588 (2020) 125041

Rudiyanto, et al.

Table 3
Optimized main drying and wetting parameters of the FXWJD model obtained from hysteretic water retention data of 21 soils.
Soils

Sources of dataθsd
(cm3 cm−3)αd(cm−1)

n
(-)

m
(-)

θsw
(cm3 cm−3)αw(cm−1)

Fitted Main Drying Curve
RMSE
R2

Fitted Main Wetting Curve
RMSE
R2

Tottori dune sand
Sand
Avondale clay loam
Gilham dune sandGuelph loam
White silica sandIda silt loam 15 cm
Ida silt loam > 15 cm
Fine sandRubicon sandy loam
Fine sandsand 4890
clay loam 1301
loam 4710
sand 2310
sandy clay loam 1300
silty clay loam 2020
silty clay loam 2021
silty clay loam 2022
silty clay loam 4700
sand 1410
silt loam 4870
Houston sandCoarse sand(Rudiyanto et al., 2015)(Vachaud and Thony, 1971)(Watson et al., 1975)(Mualem, 1976b)(Mualem, 1976b)(Huang et al., 2005)(Mualem, 1976b)(Mualem, 1976b)(Maqsoud et al., 2012)(Mualem, 1976b)(Yang et al., 2004)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Zhou, 2013)(Viaene et al., 1994)

0.376
0.367
0.360
0.303
0.507
0.355
0.548
0.531
0.273
0.371
0.407
0.304
0.380
0.400
0.355
0.386
NA
NA
NA
0.400
0.357
0.480

0.039
0.021
0.017
0.036
0.015
0.012
0.008
0.013
0.050
0.015
0.048
0.009
0.019
0.029
0.028
0.025
NA
NA
NA
0.021
0.047
0.011
0.065
0.030

8.838
10.0004.790
10.000
6.963
6.288
1.100
1.158
7.436
10.000
6.014
6.380
1.100
1.796
10.000
1.343
NA
NA
NA
1.266
8.866
1.941
6.808
5.112

0.851
0.710
0.127
0.479
0.214
0.572
0.702
0.772
0.985
0.255
1.000
0.769
0.407
0.964
1.000
0.469
NA
NA
NA
1.000
0.842
0.343
1.000
1.000

0.338
0.351
0.354
0.294
0.433
0.344
0.499
0.492
0.217
0.370
0.397
0.300
0.367
0.380
0.291
0.383
NA
NA
NA
0.399
0.350
0.460
1.000
0.231

0.064
0.043
0.032
0.063
0.029
0.032
0.011
0.015
0.078
0.048
0.103
0.017
0.050
0.030
0.053
0.043
NA
NA
NA
0.030
0.071
0.033
0.143
0.045
Mean
SD

0.014
0.011
0.006
0.011
0.012
0.013
0.007
0.007
0.007
0.017
0.009
0.005
0.006
0.006
0.018
0.009
NA
NA
NA
0.013
0.016
0.010
0.044
0.010
0.012
0.008

0.013
0.015
0.004
0.008
0.012
0.012
0.013
0.008
0.012
0.019
0.022
0.007
0.019
0.008
0.017
0.008
NA
NA
NA
0.005
0.010
0.012
0.056
0.012
0.014
0.011

0.260

45 soil datasets, Mualem (1976a) found an optimal value of 0.5 for L. In
another study, Schaap and Leij (1999) found the optimal value of L in
the vG-Mualem model to be −1, based on an analysis of 235 laboratory
samples from the UNSODA database. The negative L value implies anincrease in hydraulic conductivity at low water contents (see Fig. 1 in
Schaap and Leij, 1999) as a compensation for the absence of film flow.
Thus, the negative L value may indicate that the vG-Mualem model has
an inadequate model structure to produce an appropriate shape of thehydraulic conductivity curve.
Fig. 5 shows the distribution of L values of the RMSS, FXWJD and
PDI models, for various soil texture classes. The sand, sandy loam,
loamy sand, silty clay loam, and clay loam textures show a unimodal
distribution of L for both models, while the other textures have a bimodal distribution. Unimodal (or narrow) distributions indicate that L
could be identified well, while bimodal (or high variance) distributions
suggest L to be highly uncertain. The unimodal distribution of L values
indicates a median L value between 1 and 3.8 for the RMSS model andbetween 1.1 and 3.5 for the FXWJD model. For the PDI model, themedian L value for sand was close to 3.5; while L values for the remaining soils were mostly between 7 and 10. The median L value as afunction of texture class can be used as an good estimate in the absenceof hydraulic conductivity data.
The distribution of L values when all soil samples are combined isshown in Fig. 6a. The plots show that both models have the same multi/
bimodal distribution with a dominant distribution ranging between 0.1
and 8, with median L values of about 4.4 for the RMSS model, 3.1 for
the FXWDJ model and 10 for PDI model. When samples with L value
above 8 are excluded, given their higher uncertainty, the distributions
are given in Fig. 6b. Results now show that the median value of L isabout 3.6 for RMSS and 3.0 for the FXWJD model. The median L value
found in this study for FXWJD was slightly lower than the value of 3.5
recommended by Wang et al. (2018). Interestingly, the L = 3 value we
obtained is identical to the theoretical value of a power law of the film
conductivity as a function of the film water content (Langmuir, 1938;
Tokunaga, 2009; Wang et al., 2013). Thus, the term of SeL of Kr (Eq. (4))
could be interpreted as the contribution of film flow while the term [1(1-Γ1/m)1-1/n]2 may be interpreted as the contribution of capillary flow
to the hydraulic conductivity.

0.995
0.991
0.963
0.990
0.988
0.994
0.997
0.997
0.987
0.995
0.997
0.997
0.964
0.995
0.995
0.926
NA
NA
NA
0.838
0.991
0.995
0.991
0.992
0.980
0.037

0.992
0.983
0.995
0.996
0.986
0.992
0.987
0.996
0.984
0.956
0.989
0.996
0.661
0.991
0.994
0.877
NA
NA
NA
0.987
0.995
0.998
0.978
0.980
0.967
0.075

We further performed a sensitivity analysis of the L parameter on
the predicted K(h) curves for both the RMSS and FXWJD models. The
water retention parameters and Ks were fixed based on the parameter
optimization results (section 4.1), while L was treated as an independent variable ranging from 0 to 10. RMSElog(K) values were calculated using a total of 4160 predicted and observed K(h) data points
from 232 soil samples covering 11 soil textures. The continuous curves
of RMSElog(K) as a function of L for the RMSS and FXWJD models arepresented in Fig. 7. Results show that RMSElog(K) had a minimum of
0.73 at L = 4 for the RMSS model, and FXWJD a minimum of 0.71 at
L = 3, which is in line with the median L values. These values areslightly higher than RMSElog(K) values of the best fitted data, whichwere 0.52 and 0.41 for the RMSS and FXWJD formulations, respectively.
4.3. Prediction of hysteresis
Fig. 8a,b shows observed and fitted main drying and wetting water
retention curves (MDC and MWC) for the RMSS model, both on a log
scale (from saturation to the dry end) and on a linear scale (from saturation to a pressure head of −100 cm). Similar results for the extended FXWJD model are shown in Fig. 9a,b. The observed and predicted scanning loops (drying and wetting scanning curves, DSC and
WSC) are also shown in Figs. 8 and 9. The DSC initially moves from areversal point on the MWC, i.e., (h, θ) = (-17.5, 0.227), to an end point(-36.5, 0.115) of the MDC, followed by a WSC back from the end point
to the initial reversal point on the MWC. Both models show similar
results, with the fitted main curves and predicted scanning loops both
agreeing well with the observed data, except in the dry range where
they underestimate the main curves. Fig. 8c,d show the hydraulicconductivity curves as a function of pressure head and water content,
respectively, corresponding to the main and scanning loops of water
retention in Fig. 8a. Similar hydraulic conductivity results for the extended FXWJD model are shown in Fig. 9c,d, corresponding to thewater retention data in Fig. 9a. The K(h) curves of both models show
two sub-curves reflecting the contributions of capillary flow at relatively high pressure heads and film conductivities at the lower pressure
heads in the dry range.
15

Journal of Hydrology 588 (2020) 125041

Rudiyanto, et al.

Rudiyanto et al. (2015) showed better agreement between observed
and fitted data in the dry region of Tottori dune sand.
Another concern relates to air entrapment during rewetting. The
RMSS and extended FXWJD models produce unclosed water retention
curves at saturation as shown by the data. The naturally saturated watercontent of the main drying curve is known to be generally larger than
the saturated water content of the main wetting curve due to air entrapment during rewetting in the capillary region. By comparison, the
RMSS and extended FXWJD scanning curves both follow a closed loop
between the drying and wetting scanning curves, in agreement withobserved data. This feature was also shown in our previous study(Rudiyanto et al., 2015).
The main and scanning K(h) branches showed hysteretic behavior inboth the capillary and film conductivity regions, which reflects hysteresis in capillary and adsorption water retention. These result are
different from our previous study (Rudiyanto et al., 2015) in which we
assumed hysteretic K(h) curves for only the capillary conductivity and a
non-hysteretic K(h) for the film conductivity. In terms of K(θ), thisshows that the RMSS model does not lead to hysteresis in the main andall scanning branches over the entire range of water contents. This is inaccordance with the often assumed non-hysteretic nature of K(θ) (e.g.,
Kool and Parker, 1987; Mualem, 1986). Fig. 9d shows that the extended
FXWJD model also produces almost no hysteresis in K(θ).
Finally, we note that the RMSS and extended FXWJD hydraulicparameters for hysteresis in the MDC and MWC are only different when
θsd > θsw and αd < αw, while the remaining parameters (n, m, Ks, andL) are assumed to be identical. These constraints lead to far fewer hydraulic parameters. When there is no entrapped air, θsd will of course
become identical to θsw. Commonly, water retention and hydraulicconductivity data are measured using drying processes, in which caseθsd, αd, n, m, Ks and L can be obtained from the drying curve, while αw(and possibly θsw is due to air entrapment issue) are still unknown. To
predict αw, a linear relationship between αd and αw could be establishedfrom the fitted main hysteretic loops of 24 soils for RMSS, and 21 soils
for the extended FXWJD involving 7 different soil textures. The data aresummarized in Tables 1 and 3 for the RMSS model and extended
FXWJD models, and further pertain to Figs. S12 and S13 in supplementary materials for the RMSS and extended FXWJD, respectively.
Results in Fig. 10 indicate a linear relationship (αw ≈ 1.8 αd) more orless independent of soil texture for both RMSS (R2 = 0.85) and theextended FXWJD (R2 = 0.84). The linear relationship was found to bestatistically significant (p < 0.001).
To test the applicability of the obtained linear relationship between
αw and αd (i.e., αw = 1.8 αd), we used the relationship to predict themain wetting scanning curves from the fitted main drying curves of 24
soils. Since we assumed that θsw is unknown, θsw = θsd. Predicted andobserved main wetting curves are shown in Figs. S12 and S13 of the
supplementary materials for the RMSS and the extended FXWJD, respectively, while the optimized parameters, RMSEθ, and the goodnessof fit (R2) for the RMSS and extended FXWJD models are listed inTables 2 and 4 respectively. The data show that the mean RMSEθ of thepredicted main wetting curve for both RMSS and the extended FXWJD
models are quite similar (0.037), slightly larger than the RMSEθ valuesof the fitted main wetting curve (0.017 for RMSS and 0.014 for theextended FXWJD). The R2 values of the predicted main wetting curvefor both RMSS and the extended FXWJD model were relatively large(R2 > 0.9), thus suggesting that the approximation αw = 1.8 αd ifs
useful for generating hysteresis loops for both the RMSS and extended
FXWJD models when relevant data are not available. The value of 1.8 isconsistent with values of 1.8 or 2.0 often assumed for hysteretic water
retention behavior (e.g., Luckner et al., 1989).

Fig. 10. Linear relationships between αd and αw as obtained from 24 and 21
soils covering 7 different soil textures for the (a) RMSS and (b) extended FXWJD
models, respectively.

We now briefly summarize several aspects of the RMSS and extended FXWJD models related to hysteresis. One concerns the degree ofhysteresis in water retention in both the capillary and adsorptive regions. The RMSS and extended FXWJD models both lead to hysteretic
capillary water retention when the water content drops and rises rapidly (h > -50 cm) along the main drying and wetting curves. The
main drying and wetting curves merge to a zero water content at thedry end. This suggests that hysteresis also occurs in the adsorptive region. Although beyond the observed data in Fig. 8a and 9a, previous
studies (Arthur et al., 2015; Davis et al., 2009; Globus and Neusypina,
2006; Ning and Morteza, 2015; Prunty and Bell, 2007) reported thathysteresis could also occur in the adsorption region. Our proposed
RMSS formulation and the extended FXWJD model are different fromour previous model (Rudiyanto et al., 2015), which could only produce
hysteretic water retention in the capillary region and assumed nonhysteretic behavior in the adsorptive region. Still, the formulation by

5. Conclusions
In summary, the proposed (RMSS and extended FXWJD) models are
relatively simple and straightforward in that they can be incorporated
16

Journal of Hydrology 588 (2020) 125041

Rudiyanto, et al.

Table 4
Hydraulic FXWJD model parameters (αw and θsw = θsd) of the main wetting curve as predicted from the optimized main drying curve of 21 soils.
Soils

Sources of data

θsd
(cm3
cm−3)

αd
(cm−1)

n
(-)

m
(-)

θsw
(cm3
cm−3)

αw
(cm−1)

Fitted Main Drying Curve
RMSE
R2

Predicted Main Wetting Curve
RMSE
R2

Tottori dune sand
Sand

(Rudiyanto et al., 2015)(Vachaud and Thony,
1971)(Watson et al., 1975)(Mualem, 1976b)(Mualem, 1976b)(Huang et al., 2005)(Mualem, 1976b)(Mualem, 1976b)(Maqsoud et al., 2012)(Mualem, 1976b)(Yang et al., 2004)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Nemes et al., 2001)(Zhou, 2013)(Viaene et al., 1994)

0.374
0.363

0.037
0.021

9.867
10.000

1.000
0.632

0.374
0.363

0.067
0.038

0.005
0.008

0.998
0.991

0.029
0.062

0.987
0.936

0.351
0.305
0.518
0.358
0.542
0.522
0.266
0.381
0.409
0.302
0.373
0.364
0.350
0.397
NA
NA
NA
0.382
0.359
0.467
1.000
0.256

0.013
0.034
0.017
0.011
0.008
0.013
0.051
0.014
0.050
0.009
0.025
0.026
0.028
0.027
NA
NA
NA
0.018
0.045
0.010
0.062
0.029

7.171
10.000
2.697
5.774
1.427
1.651
10.00010.000
5.350
7.512
1.862
2.807
10.000
1.100
NA
NA
NA
1.261
10.000
3.692
9.577
7.411

0.121
0.621
0.368
0.768
0.595
0.593
0.756
0.327
1.000
0.732
0.248
0.694
1.000
0.561
NA
NA
NA
1.000
0.983
0.224
1.000
0.955

0.351
0.305
0.518
0.358
0.542
0.522
0.266
0.381
0.409
0.302
0.373
0.364
0.350
0.397
NA
NA
NA
0.382
0.359
0.467
1.000
0.256

0.024
0.061
0.031
0.020
0.014
0.024
0.092
0.025
0.090
0.016
0.044
0.047
0.050
0.048
NA
NA
NA
0.032
0.081
0.018
0.112
0.051
Mean
SD

0.002
0.003
0.005
0.003
0.003
0.003
0.006
0.004
0.008
0.005
0.006
0.004
0.019
0.009
NA
NA
NA
0.011
0.006
0.003
0.022
0.001
0.007
0.005

0.996
0.998
0.998
0.999
0.999
1.000
0.990
0.997
0.998
0.997
0.967
0.997
0.995
0.928
NA
NA
NA
0.855
0.997
0.999
0.997
1.000
0.986
0.034

0.007
0.027
0.070
0.055
0.037
0.027
0.021
0.042
0.031
0.015
0.029
0.042
0.045
0.009
NA
NA
NA
0.017
0.050
0.041
0.088
0.025
0.037
0.020

0.986
0.995
0.962
0.875
0.980
0.985
0.981
0.903
0.995
0.985
0.562
0.973
0.972
0.859
NA
NA
NA
0.985
0.955
0.991
0.957
0.953
0.942
0.095

Avondale clay loam
Gilham dune sandGuelph loam
White silica sandIda silt loam 15 cm
Ida silt loam > 15 cm
Fine sandRubicon sandy loam
Fine sandsand 4890
clay loam 1301
loam 4710
sand 2310
sandy clay loam 1300
silty clay loam 2020
silty clay loam 2021
silty clay loam 2022
silty clay loam 4700
sand 1410
silt loam 4870
Houston sandCoarse sandTable A1
Synopsis of the RMSS, FXWJD and PDI functions used in this study, as well as of vG-Mualem equations.
Soil hydraulic model
RMSS (this study)

Soil water retention(h) =Se (h) =Hydraulic conductivitys Se (h )
(h)
(h0 )
1
(h0 )

K (h ) = Ks Se l 1(1

K (h) = Ks Senl 1

(1

Notes
1 2
1 1 n
m

Parameters: θs, α, n, m, Ks, L
Constant: h0 = − 0.63 107 cm

1 2
1 1 n
m

Parameters : θs, α, n, m, Ks, L
Constants:
hr = − 1500 cm
h0 = − 0.63 107 cm
L = 3.5 (in this study, L is treated as a parameter))(h) = (ln[exp(1) + | h|n ]) m
FXWJD(Wang et al., 2018)(h) =s Se (h )ln 1 +Se (h) = 1hhrhln 1 + 0
hr

Sen (h) =

(h)
1

)

(h0 )
(h0 )

(h) (h ) = (ln[exp(1) + | h|n ]) m
PDI

(Iden and Durner, 2014; Peters, 2014, 2013)(h) = ( sr ) Sc (h) +r Sa (h)K (h) = K c (h) + Kf (h )

(h0 )
(h0 )
(h) = (1 + | h|n ) m

Sc (h ) =(h)1

K c (h ) = Ksc SecL 1

Sa (h)

Kf (h) = Ksf

h
hra
h0
ln 1 +
haln 1 +

=

Xm 1

for h

1

( )
h0
ha

1
1

1 m 2
m
1
0m

Parameters: θs, θr, α, n, Ksc, L, Ksf (m = 1–1/n)Constants:
h0 = − 6300000 cm
ha = 1/α1.5(1 Sa)hafor h > h a1

Xm = 1
van Genuchten –Mualem (van Genuchten, 1980)ln(2)hln 1 + 0ha(h) = ( sr ) Se (h ) +Se (h) = (1 + | h|n ) mr

K (h ) = Ks Se l 1

immediately into existing flow models without requiring additionalparameters. The RMSS model describes the complete soil water retention and hydraulic conductivity curves, while accounting for capillary
and adsorption water retention, a zero water content in the very dry
region, capillary and film conductivity, and hysteresis. Our tests on 232
soils covering soil textural classes from sand to clay show that themodel works well especially for coarse-textured soils having relatively
narrow particle size distributions (such as sands, loamy sands, sandy
loams). The tortuosity/connectivity parameter L in the hydraulic

1 m 2
m

(1 Se )

Parameters: θs, θr, α, n, Ks, L(m = 1–1/n)

conductivity function was found to be a function of soil texture. The
median value of L was shown to be equal to 3.6 for our tested soils. We
further accounted for hysteresis in both the capillary and adsorptive
regions of the water retention curve, as well as for hysteresis in both thecapillary and film regions of the hydraulic conductivity as a function of
the pressure head, but not as a function of the water content. Hysteresis
was modeled using closed scanning loops. To generate hysteresis, theproposed RMSS and extended FXWJD models only required two additional parameters: θsw and αw. When no data are available of the main17

Journal of Hydrology 588 (2020) 125041

Rudiyanto, et al.

wetting curve, one can still model hysteresis assuming no entrapped air(in which case θsd = θsw), with αw to be approximated using αw = 1.8
αd. Changes in the hydraulic conductivity due to vapor flow (Peters,
2013; Rudiyanto et al., 2015; Saito et al., 2006) can be incorporated asneeded in the RMSS model. We did not consider vapor flow in this study
since the range of vapor conductivity is beyond the observed databasethat we used. We conclude that the inclusion of the proposed functions
in soil water flow modeling efforts can improve dynamic flow simulations from saturation to dry conditions of especially coarse-textured
soils.
All of the soil water retention and hydraulic conductivity data usedin this study are available for download at https://doi.org/10.7910/
DVN/URFRM9. The hysteretic water retention data are available fordownload at https://doi.org/10.7910/DVN/J2V6HH.

by constraining capillary bundle models to a maximum pore size. Adv. Water Resour.
85, 86–92. https://doi.org/10.1016/j.advwatres.2015.09.005.
Jury, W., Horton, R., 2004. Soil Physics, 6th Edition. John Wiley & Sons, New York.
Khlosi, M., Cornelis, W.M., Gabriels, D., Sin, G., 2006. Simple modification to describe thesoil water retention curve between saturation and oven dryness. Water Resour. Res.
https://doi.org/10.1029/2005WR004699.
Kool, J.B., Parker, J.C., 1987. Development and evaluation of closed-form expressions forhysteretic soil hydraulic properties. Water Resour. Res. 23, 105–114. https://doi.org/
10.1029/WR023i001p00105.
Kosugi, K., 1996. Lognormal distribution model for unsaturated soil hydraulic properties.
Water Resour. Res. 32, 2697–2703. https://doi.org/10.1029/96WR01776.
Kosugi, K., Hopmans, J., Dane, J.H., 2002. 3.3. 4 Parametric Models. Power 1–20.
Langmuir, I., 1938. Repulsive Forces between Charged Surfaces In Water, And The Cause
Of The Jones-Ray Effect. Science 80, 430–432. https://doi.org/10.1126/science.88.
2288.430.
Lebeau, M., Konrad, J.-M., 2010. A new capillary and thin film flow model for predicting
the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 46.
https://doi.org/10.1029/2010WR009092.
Leij, F.J., Russell, W.B., Lesch, S.M., 1997. Closed-form expressions for water retention
and conductivity data. Groundwater 35, 848–858. https://doi.org/10.1111/j.17456584.1997.tb00153.x.
Liao, K., Lai, X., Zhou, Z., Zhu, Q., Han, Q., 2018. A simple and improved model for
describing soil hydraulic properties from saturation to oven dryness. Vadose Zo. J.
17, 1–8. https://doi.org/10.2136/vzj2018.04.0082.
Luckner, L., Van Genuchten, M.T., Nielsen, D.R., 1989. A consistent set of parametricmodels for the two-phase flow of immiscible fluids in the subsurface. Water Resour.
Res. 25, 2187–2193. https://doi.org/10.1029/WR025i010p02187.
Maqsoud, A., Bussière, B., Aubertin, M., Mbonimpa, M., 2012. Predicting hysteresis of thewater retention curve from basic properties of granular soils. Geotech. Geol. Eng. 30,
1147–1159. https://doi.org/10.1007/s10706-012-9529-y.
Minasny, B., Field, D.J., 2005. Estimating soil hydraulic properties and their uncertainty:
The use of stochastic simulation in the inverse modelling of the evaporation method.
Geoderma 126, 277–290. https://doi.org/10.1016/j.geoderma.2004.09.015.
Morel-Seytoux, H.J., Nimmo, J.R., 1999. Soil water retention and maximum capillary
drive from saturation to oven dryness. Water Resour. Res. 35, 2031–2041. https://doi.org/10.1029/1999WR900121.
Mualem, Y., 1986. Hydraulic conductivity of unsaturated soils: Prediction and formulas,
in: Klute, A. (Ed.), Methods in Soil Analysis. Part 1, 2nd Edition. Agron. Monogr. 9.
ASA and SSSA, Madison, WI, WI, pp. 799–824.
Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated
porous media. Water Resour. Res. 12, 513–522. https://doi.org/10.1029/
WR012i003p00513.
Mualem, Y., 1976b. A catalogue of the hydraulic properties of unsaturated soils.
Mualem, Y., Dagan, G., 1978. Hydraulic conductivity of soils: unified approach to thestatistical models. Soil Sci. Soc. Am. J. 42, 392–395. https://doi.org/10.2136/
sssaj1978.03615995004200030003x.
Nemes, A., Schaap, M.G., Leij, F.J., Wösten, J.H.M., 2001. Description of the unsaturated
soil hydraulic database UNSODA version 2.0. J. Hydrol. 251, 151–162. https://doi.
org/10.1016/S0022-1694(01)00465-6.
Nielsen, D.R., Th. Van Genuchten, M., Biggar, J.W., 1986. Water flow and solute transport
processes in the unsaturated zone. Water Resour. Res. 22, 89S-108S. https://doi.org/
10.1029/WR022i09Sp0089S.
Ning, L., Morteza, K., 2015. Mechanisms for soil-water retention and hysteresis at highsuction range. J. Geotech. Geoenvironmental Eng. 141, 4015032. https://doi.org/10.
1061/(ASCE)GT.1943-5606.0001325.
Pachepsky, Y.A., Shcherbakov, R.A., Varallyay, G., Rajkai, K., Pachepsky, Y., Scherbakov,
R., Varallyay, G., Rajkai, K., 1984. On obtaining soil hydraulic conductivity curves
from water retention curves (in Russian). Pochvovedenie 10, 60–72.
Parker, J.C., Lenhard, R.J., 1987. A model for hysteretic constitutive relations governing
multiphase flow: 1 Saturation-pressure relations. Water Resour. Res. 23, 2187–2196.
https://doi.org/10.1029/WR023i012p02187.
Peters, A., 2014. Reply to comment by S. Iden and W. Durner on “Simple consistentmodels for water retention and hydraulic conductivity in the complete moisture
range”. Water Resour. Res. 50, 7535–7539. https://doi.org/10.1002/
2014WR016107.
Peters, A., 2013. Simple consistent models for water retention and hydraulic conductivity
in the complete moisture range. Water Resour. Res. 49, 6765–6780. https://doi.org/
https://doi.org/10.1002/wrcr.20548.
Prunty, L., Bell, J., 2007. Soil water hysteresis at low potential. Pedosphere 17, 436–444.
https://doi.org/10.1016/S1002-0160(07)60053-8.
Ross, P.J., Williams, J., Bristow, K.L., 1991. Equation for extending water-retention
curves to dryness. Soil Sci. Soc. Am. J. 55, 923–927. https://doi.org/10.2136/
sssaj1991.03615995005500040004x.
Rossi, C., Nimmo, J.R., 1994. Modeling of soil water retention from saturation to ovendryness. Water Resour. Res. 30, 701–708. https://doi.org/10.1029/93WR03238.
Rudiyanto, Sakai, M., van Genuchten, M.T., Alazba, A.A., Setiawan, B.I., Minasny, B.,
2015. A complete soil hydraulic model accounting for capillary and adsorptive water
retention, capillary and film conductivity, and hysteresis. Water Resour. Res. 51,
8757–8772. https://doi.org/10.1002/2015WR017703.
Saito, H., Šimůnek, J., Mohanty, B.P., 2006. Numerical analysis of coupled water, vapor,and heat transport in the vadose zone. Vadose Zo. J. 5, 784–800. https://doi.org/10.
2136/vzj2006.0007.
Samaniego, L., Thober, S., Kumar, R., Wanders, N., Rakovec, O., Pan, M., Zink, M.,
Sheffield, J., Wood, E.F., Marx, A., 2018. Anthropogenic warming exacerbates
European soil moisture droughts. Nat. Clim. Chang. 8, 421–426. https://doi.org/10.
1038/s41558-018-0138-5.

CrediT authorship contribution statement
Rudiyanto: Conceptualization, Methodology, Software, Formal
analysis, Investigation, Data curation, Writing - original draft,
Visualization. Budiman Minasny: Conceptualization, Writing - review& editing. Ramisah M. Shah: Supervision. Budi I. Setiawan:Supervision. Martinus Th. van Genuchten: Writing - review & editing.
Declaration of Competing Interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence the work reported in this paper.
Appendix A. Supplementary data
Supplementary data to this article can be found online at https://doi.org/10.1016/j.jhydrol.2020.125041.
References
Arthur, E., Tuller, M., Moldrup, P., Jensen, D.K., De Jonge, L.W., 2015. Prediction of clay
content from water vapour sorption isotherms considering hysteresis and soil organic
matter content. Eur. J. Soil Sci. 66, 206–217. https://doi.org/10.1111/ejss.12191.
Assouline, S., Or, D., 2013. Conceptual and parametric representation of soil hydraulicproperties: a review. Vadose Zo. J. 12https://doi.org/10.2136/vzj2013.07.0121.
vzj2013.07.0121.
Viaene, P., Vereecken, H., Diels, J., Feyen, J., 1994. A statistical analysis of six hysteresis
models for the moisture retention characteristic. Soil Sci. 157.
Brooks, R., Corey, A., 1964. Hydraulic properties of porous media. Hydrol. Pap. Color.
State Univ.
Campbell, G.S., Shiozawa, S., 1992. Prediction of hydraulic properties of soils using
particle-size distribution and bulk density data. Indirect Methods for Estimating theHydraulic Properties. of Unsaturated Soils.
Davis, D.D., Horton, R., Heitman, J.L., Ren, T., 2009. Wettability and hysteresis effects onwater sorption in relatively dry soil. Soil Sci. Soc. Am. J. 73, 1947–1951. https://doi.
org/10.2136/sssaj2009.00028N.
Dourado Neto, D., de Jong van Lier, Q., van Genuchten, M.T., 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 Zo. J. 10, 618–623. https://doi.org/10.2136/vzj2009.0191.
Duan, Q., Sorooshian, S., Gupta, V.K., 1994. Optimal use of the SCE-UA global optimization method for calibrating watershed models. J. Hydrol. 158, 265–284. https://doi.org/10.1016/0022-1694(94)90057-4.
Durner, W., Flühler, H., 2005. Soil Hydraulic Properties. Encycl. Hydrol. Sci., Major
Reference Works. https://doi.org/doi:10.1002/0470848944.hsa077c.
Fayer, M.J., Simmons, C.S., 1995. Modified soil water retention functions for all matricsuctions. Water Resour. Res. https://doi.org/10.1029/95WR00173.
Fredlund, D.G., Xing, A., 1994. Equations for the soil-water characteristic curve. Can.
Geotech. J. 31, 512–532. https://doi.org/10.1139/t94-061.
Globus, A.M., Neusypina, T.A., 2006. Determination of the water hysteresis and specificsurface of soils by electronic microhygrometry and psychrometry. Eurasian Soil Sci.
39, 270–277. https://doi.org/10.1134/S1064229306030057.
Huang, H.-C., Tan, Y.-C., Liu, C.-W., Chen, C.-H., 2005. A novel hysteresis model in unsaturated soil. Hydrol. Process. 19, 1653–1665. https://doi.org/10.1002/hyp.5594.
Iden, S.C., Durner, W., 2014. Comment on “Simple consistent models for water retention
and hydraulic conductivity in the complete moisture range” by A. Peters. Water
Resour. Res. 50, 7530–7534. https://doi.org/10.1002/2014WR015937.
Iden, S.C., Peters, A., Durner, W., 2015. Improving prediction of hydraulic conductivity

18

Journal of Hydrology 588 (2020) 125041

Rudiyanto, et al.
Schaap, M.G., Leij, F.J., 1999. Improved prediction of unsaturated hydraulic conductivity
with the Mualem-van Genuchten model. Soil Sci. Soc. Am. J. 64, 843–851. https://doi.org/10.2136/sssaj2000.643843x.
Schindler, U., Müller, L., 2006. Simplifying the evaporation method for quantifying soil
hydraulic properties. J. Plant Nutr. Soil Sci. 169, 623–629. https://doi.org/10.1002/
jpln.200521895.
Schneider, M., Goss, K.-U.U., 2012. Prediction of the water sorption isotherm in air dry
soils. Geoderma 170, 64–69. https://doi.org/10.1016/j.geoderma.2011.10.008.
Tiktak, A., Bouten, W., 1992. Modelling soil water dynamics in a forested ecosystem. III:
Model description and evaluation of discretization. Hydrol. Process. 6, 455–465.
https://doi.org/10.1002/hyp.3360060407.
Tokunaga, T.K., 2009. Hydraulic properties of adsorbed water films in unsaturated porous
media. Water Resour. Res. 45. https://doi.org/10.1029/2009WR007734.
Vachaud, G., Thony, J.-L., 1971. Hysteresis during infiltration and redistribution in a soilcolumn at different initial water contents. Water Resour. Res. 7, 111–127. https://doi.org/10.1029/WR007i001p00111.
van Genuchten, M., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 8, 892–898. https://doi.org/10.
2136/sssaj1980.03615995004400050002x.
van Genuchten, M.T., Nielsen, D.R., 1985. On describing and predicting the hydraulicproperties of unsaturated soils. Ann. Geophys. 3, 615–628.
Vogel, T., van Genuchten, M.T., Cislerova, M., 2000. Effect of the shape of the soil hydraulic functions near saturation on variably-saturated flow predictions. Adv. Water
Resour. 24, 133–144. https://doi.org/10.1016/S0309-1708(00)00037-3.

Wang, Y., Jin, M., Deng, Z., 2018. Alternative model for predicting soil hydraulic conductivity over the complete moisture range. Water Resour. Res. 54, 6860–6876.
https://doi.org/10.1029/2018WR023037.
Wang, Y., Ma, J., Guan, H., 2016. A mathematically continuous model for describing thehydraulic properties of unsaturated porous media over the entire range of matricsuctions. J. Hydrol. 541, 873–888.
Wang, Y., Ma, J., Zhang, Y., Zhao, M., Edmunds, W.M., 2013. A new theoretical modelaccounting for film flow in unsaturated porous media. Water Resour. Res. 49,
5021–5028. https://doi.org/10.1002/wrcr.20390.
Watson, K.K., Reginato, R.J., Jackson, R.D., 1975. Soil water hysteresis in a field soil. Soil
Sci. Soc. Am. J. 39, 242–246. https://doi.org/10.2136/sssaj1975.
03615995003900020010x.
Weber, T.K.D., Durner, W., Streck, T., Diamantopoulos, E., 2019. A modular framework
for modeling unsaturated soil hydraulic properties over the full moisture range.
Water Resour. Res. 55, 4994–5011. https://doi.org/10.1029/2018WR024584.
Yang, H., Rahardjo, H., Leong, E.-C., Fredlund, D.G., 2004. Factors affecting drying andwetting soil-water characteristic curves of sandy soils. Can. Geotech. J. 41, 908–920.
https://doi.org/10.1139/t04-042.
Zhang, Z.F., 2011. Soil water retention and relative permeability for conditions fromoven-dry to full saturation. Vadose Zo. J. 10, 1299–1308. https://doi.org/10.2136/
vzj2011.0019.
Zhou, A.-N., 2013. A contact angle-dependent hysteresis model for soil–water retention
behaviour. Comput. Geotech. 49, 36–42. https://doi.org/10.1016/j.compgeo.2012.
10.004.

19