Received: 9 March 2021

Accepted: 9 September 2021

DOI: 10.1002/vzj2.20166
Vadose Zone Journal

ORIGINAL RESEARCH ARTICLE

Parasite inversion for determining the coefficients andtime-validity of Philip’s two-term infiltration equation
Parakh Jaiswal1
Jirka Šimůnek4

Jan Vanderborght3
Yifu Gao1
Mehdi Rahmati2,3Harry Vereecken3
Jasper A. Vrugt1,5

1 Dep. of Civil and Environmental

Engineering, Univ. of California, Irvine, CA
92697, USA
2 Dep. of Soil Science and Engineering,
Faculty of Agriculture, Univ. of Maragheh,
Maragheh, Iran
3 Forschungszentrum Jülich, Institute ofBio- and Geosciences: Agrosphere (IBG-3),
Jülich, Germany
4 Dep. of Environmental Sciences, Univ. of
California, Riverside, CA 92521, USA
5 Dep. of Earth System Science, Univ. of

California, Irvine, CA 92697, USA
Correspondence
Harry Vereecken, Forschungszentrum Jülich,
Institute of Bio- and Geosciences: Agrosphere (IBG-3), Jülich, Germany.
Email: h.vereecken@fz-juelich.de
Assigned to Associate Editor Nima Shokri.

Abstract
Many different equations have been proposed to describe quantitatively onedimensional soil water infiltration. The unknown coefficients of these equations characterize soil hydraulic properties and may be estimated from a n record, {𝑡̃𝑖 , 𝐼̃𝑖 }𝑛𝑖=1 ,of cumulative infiltration measurements using curve fitting techniques. The two-term√
infiltration equation, 𝐼(𝑡) = 𝑆 𝑡 + 𝑐𝐾s 𝑡, of Philip has been widely used to describe
measured infiltration data. This function enjoys a solid mathematical–physical underpinning and admits a closed-form solution for the soil sorptivity, S [L T−1/2 ], and multiple, 𝑐 [−], of the saturated hydraulic conductivity, 𝐾s [L T−1 ]. However, Philip’s
two-term equation has a limited time validity, 𝑡valid [T], and thus cumulative infil̃ 𝑡̃), beyond 𝑡 = 𝑡valid will corrupt the estimates of S and 𝐾s . This
tration data, 𝐼(
paper introduces a novel method for estimating S, c, 𝐾s , and 𝑡valid of Philip’s twoterm infiltration equation. This method, coined parasite inversion, use as vehicle
Parlange’s three-parameter infiltration equation. As prerequisite to our method, we
present as secondary contribution an exact, robust and efficient numerical solutionof Parlange’s infiltration equation. This solution admits Bayesian parameter estimation with the DiffeRential Evolution Adaptive Metropolis (DREAM) algorithm andyields as byproduct the marginal distribution of Parlange’s β parameter. We evaluate
our method for 12 USDA soil types using synthetic infiltration data simulated with
HYDRUS-1D. An excellent match is observed between the inferred values of S and𝐾s and their “true” values known beforehand. Furthermore, our estimates of c and𝑡valid correlate well with soil texture, corroborate linearity of the 𝑐(β) relationshipfor 0 ≤ 𝑡 ≤ 𝑡valid , and fall within reported ranges. A cumulative vertical infiltrationof about 2.5 cm may serve as guideline for the time-validity of Philip’s two-term
infiltration equation.

Abbreviations: BIC, Bayesian information criterion; DREAM, DiffeRential Evolution Adaptive Metropolis; MAP, maximum a posteriori density; MVG,
Mualem–van Genuchten; SWIG, Soil Water Infiltration Global
This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided
the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.
© 2021 The Authors. Vadose Zone Journal published by Wiley Periodicals LLC on behalf of Soil Science Society of America
Vadose Zone J. 2021;e20166.
https://doi.org/10.1002/vzj2.20166

wileyonlinelibrary.com/journal/vzj2

1 of 27

2 of 27

1

JAISWAL ET AL.

Vadose Zone Journal

INTRODUCTION AND SCOPE
Core Ideas

The topic of water infiltration into unsaturated soils has
attracted the attention of vadose zone hydrologists early
on. This process lends itself for detailed experimentation,mathematical–physical analysis, and spatiotemporal modeling and is central to hydrologic studies by water resourceand hydraulic engineers, hydrogeophysicists, hydrogeologists, plant scientists, agronomists, and ecologists (Philip,
1969b). Infiltration not only determines how much rainfall(irrigation) will accumulate on the soil’s surface or runoff
elsewhere following topographical gradients, but it also exerts
control on the soil moisture status of the root zone, which, inturn, regulates processes such as the absorption of soil waterby plant roots, evaporation, sensible heat fluxes, and groundwater recharge (Liu et al., 2011). Infiltration will also govern the transport of pollutants through the vadose zone, slope
stability, and land subsidence in areas that rely heavily on
groundwater pumping as primary water source (Chen et al.,2014). Farmers in the Central Valley of California, for example, rely heavily on infiltration to maintain an adequate soil
moisture status of the root zone in pursuit of an optimal crop
yield (Hellin & Schrader, 2003; Schoups et al., 2005).
After the topsoil has absorbed sufficient rainfall (or irrigation water), the water moves downward in the soil profile
under the influence of gravity and/or capillary action, soaking
and/or filling up the pore space, replenishing the root zone,
and possibly seeping into rocks through cracks. Infiltrability
is a term native to soil physics and hydrology that defines
the maximum rate at which rain or irrigation water can beabsorbed by a soil under given conditions (Horton, 1933).
The infiltration rate under natural conditions is not necessarily
readily deducible from soil moisture measurements, certainly
at larger spatial scales, and infiltration modeling with the aidof measurable quantities is of fundamental importance.
In the past decades, numerous mathematical–physicalapproaches have been developed to describe water infiltration into the vadose zone (Assouline, 2013). These models
or functions can be classified as empirical, semiempirical and(quasi-) mechanistic. Empirical infiltration models are mathematical functions that have been crafted by trial and error
with the sole purpose of matching laboratory or field measured cumulative infiltration, 𝐼(𝑡) [L], data, where t [T]denotes time. The resulting black-box function offers no
insights into the underlying data-generating process nor provides evidence why the function would accurately portray then vector of measured cumulative infiltration data. Examples
include the infiltration equations of Kostiakov (1932), Huggins and Monke (1966), modified Kostiakov (Smith, 1972),
and Collis-George (1977). As the parameters of empiricalinfiltration models lack a theoretical foundation, it is typically difficult, if not impossible, to relate their fitted values
to physical soil properties. Semiempirical (gray-box) infiltra-∙ We describe a method that can determine the time
validity of infiltration equations.
∙ This method infers as well the hydraulic propertiesof unsaturated soils.
∙ We present a robust and efficient numericalsolution of Parlange’s three-parameter infiltration
equation.
∙ Posterior distribution determines how well parameters are defined by calibration to data.
∙ Two and a half centimeters of cumulative infiltration is a good proxy for time validity Philip’s twoterm equation.

tion models conserve mass with a flux–concentration (storage) relationship that relates the infiltration rate, i [L T−1 ],to the cumulative infiltration, I, which is usually hypothesized. Examples include the models of Horton (1941), Holtan(1961), Singh and Yu (1990), and Grigorjev and Iritz (1991).
Then, the last group of quasi-deterministic (white-box) infiltration models adopt a reductionist approach and describe
the downward and/or horizontal movement of water usingmass balance principles coupled with a macroscopic law
for soil moisture flow. This use of first principles offers arobust protection against invalid assumptions and conventions and includes the mechanistic models of Green and Ampt
(1911), Philip (1957, 1969a), Mein and Larson (1971, 1973),
Smith (1972), Braester (1973), Stroosnijder (1976), Brutsaert (1977), Smith and Parlange (1978), Haverkamp et al.
(1994), Salvucci and Entekhabi (1994), Basha (1999, 2000,
2002), Zhu and Mohanty (2002), Chen and Gallipoli (2004),
and Valiantzas (2010). These models vary in complexity anddegree of sophistication and differ in their assumptions about
infiltration source type (point, circular, strip, and rectangular)
and location (above- or belowground), the shape of the infiltration front, soil heterogeneity, dominant processes and flow
dynamics, the hydraulic functions, and initial and boundary
conditions.
The mathematical–physical description of the infiltration
process has benefited tremendously from the theories andheavy lifting of Dr. John Philip (1927–1999). While he has
been recognized internationally for his many contributions to
the topics of water, gas, and energy movement, his researchon infiltration is arguably the most acclaimed worldwide. Dr.
Philip’s advanced mathematical training and deep intuition
about physical processes enabled him to make fundamental
advances to the theory and practice of water infiltration. In aclassic treatise on one-dimensional infiltration, Philip (1957)
describes the cumulative vertical infiltration, I(t), as a finite

JAISWAL ET AL.

Vadose Zone Journal(convergent) series of d > 3 expansion terms:
1

3Philip (1957) postulated that coefficients a3 , a4 , . . . , ad , ofEquation 1 must satisfy the following constraint:𝑑𝐼(𝑡) = 𝑎1 𝑡 2 + 𝑎2 𝑡 + 𝑎3 𝑡 2 + ⋯ + 𝑎𝑑 𝑡 2𝑗∑= 𝑑𝑗=1 𝑎𝑗 𝑡 2 ,

(1)

θs

∫ θi

λ(θ)dθ

(2)where θi [L3 L−3 ] denotes the initial water content of the soil
(assumed uniform), θs [L3 L−3 ] is the volumetric moisture
content at saturation, and λ(θ) [L T−1/2 ] is the so-called Boltzmann variable:λ(θ) = 𝑧(θ, 𝑡)𝑡−1∕2𝑎𝑗𝑆

where 𝑎1 [L T−1/2 ], 𝑎2 [L T−1 ], 𝑎3 [L T−3/2 ], and 𝑎𝑑 [L T−d/2 ]are soil-dependent constants, and t denotes time. Philip (1957)
showed that a1 is synonymous to the sorptivity, S [L T−1/2 ],a measure of the soil’s capacity to take up and release liquids
by capillarity, and a2 is equal to a multiple, c, of the saturated
hydraulic conductivity, 𝐾s [L T−1 ], which articulates a soil’s
ability to transmit water under the influence of gravity. Philip(1969a) expressed the soil sorptivity as follows
𝑆(θi ) =3 of 27(3)

The characteristic function, 𝑧(θ, 𝑡), expresses the relationshipbetween soil depth, z [L], and the volumetric water content, θ [L3 L−3 ], at time 𝑡 ≥ 0 during an infiltration event
under gravity-free conditions, where θi ≤ θ ≤ θs . For notational convenience, we omit the explicit functional dependence of the soil sorptivity, S, on the initial moisture content. Thus, the first two coefficients, 𝑎1 and 𝑎2 , of Philip’s
time series expansion exhibit a sound physical underpinning,
adjustable to the initial and boundary conditions of the infiltration event.
The finite series of expansion terms of Philip (1957) not
only benefits a solid theoretical underpinning but also enjoys amore practical advantage. As the cumulative infiltration, 𝐼(𝑡),is written as a multiple of the coefficients, 𝑎𝑗 (𝑗 = 1, … , 𝑑),of the d expansion terms, Equation 1 admits the application of
linear regression to determine the optimum parameter values
and their associated uncertainty from a dataset, {𝑡̃𝑖 , 𝐼̃𝑖 }𝑛𝑖=1 , of ncumulative infiltration measurements, 𝐼̃1 , 𝐼̃2 , … , 𝐼̃𝑛 , observed
at 𝑡̃1 , 𝑡̃2 , … , 𝑡̃𝑛 . Indeed, as Equation 1 belongs to the classof linear regression models, we can write Philip’s infiltration function as an inner product, 𝐼(𝑡) = 𝐝(𝑡)⊤ 𝐚, of a 𝑑 × 1
design vector, 𝐝(𝑡) = [𝑡1∕2 𝑡 𝑡3∕2 ⋯ 𝑡𝑑∕2 ]⊤ , and the d vector, 𝐚 = [𝑎1 𝑎2 𝑎3 ⋯ 𝑎𝑑 ]⊤ , of parameters. If we stack thedesign vectors, 𝐝(𝑡̃) of the n measurement times, 𝑡̃1 , 𝑡̃2 , … , 𝑡̃𝑛 ,in a n × d design matrix, 𝐃, then the least squares parameter values, ̂
𝐚 = [𝑎̂1 𝑎̂2 𝑎̂3 ⋯ 𝑎̂𝑑 ]⊤ , of the measured cumu𝐚=
lative infiltration data, ̃𝐈 = [𝐼̃1 𝐼̃2 … 𝐼̃𝑛 ]⊤ , will equal, ̂⊤−1⊤̃(𝐃 𝐃) 𝐃 𝐈.

>

( 𝑎 )𝑗−1
2(4)

𝑆

for all 𝑗 = (3, … , 𝑑). This constraint implies that all coefficients a3 , a4 , . . . , ad must exceed zero, and, as a result, Equation 1 will not converge to a constant infiltration rate at large t.
Instead, Equation 4 promotes significance of the higher orderterms at the potential risk of overfitting (Volpi et al., 2017).
This is also known as Runge’s phenomenon (polynomial wiggle) and cautions against the use of a large number of expansion terms. The most popular variant of Equation 1 uses only
two expansion terms:
𝐼(𝑡) = 𝑆𝑡1∕2 + 𝑐𝐾s 𝑡(5)where the multiplicative coefficient c is soil dependent.
This two-term formulation enjoys a rigorous mathematical–
physical underpinning. This includes recent work of Huntet al. (2017), who used percolation theory to interpret S and𝑐𝐾s . The use of 𝑑 = 2 expansion terms, however, has animportant side effect. As the higher order expansion termsof Equation 1 are discarded, Philip’s two-term equation cannot describe adequately the infiltration at later times when
the higher order terms have an increasingly larger impact on
the cumulative infiltration. In other words, Philip’s two-term
equation has a limited time validity, 𝑡valid [T] and, consequently, cumulative infiltration measurements beyond 𝑡valid
should not be used for parameter estimation purposes, as thesemeasurements will corrupt estimates of the soil sorptivity, S,
and saturated hydraulic conductivity, 𝐾s .
In his theoretic treatise, Philip (1957) derived a simple analytic expression for the so-called characteristic time, 𝑡char [T],of Equation 1:(𝑡char =

𝑆
𝐾s − 𝐾i)2(6)

at which gravity and capillary action contribute equally tothe instantaneous infiltration. In this equation, 𝐾i = 𝐾(θi ) [LT−1 ] is the soil hydraulic conductivity at the initial moisture
content. Philip (1957) used the characteristic time, 𝑡char , todemarcate the window for which Equation 1 is valid. Specifically, he postulated that the finite series of 𝑑 ≥ 3 expansion terms is valid only for infiltration times up to 2𝑡char to
4𝑡char . We should be careful, however, in using these estimates for Philip’s two-term infiltration (Equation 5). Without
doubt, the truncation of (1) to 𝑑 = 2 terms must deteriorate
the time validity of the two-term formulation. Indeed, practical guidelines are warranted on the time validity, 𝑡valid , ofEquation 5.

4 of 27

Vadose Zone Journal

In this paper, we introduce a new method for determining
the values of the dimensionless coefficient c, soil properties S
and 𝐾s , and time validity 𝑡valid of Philip’s two-term infiltration
equation. We cannot infer these four quantities simultaneously
from measured cumulative infiltration curves as this constitutes an ill-posed problem. Therefore, we divide the inferencein two separate steps using an approach coined parasite inversion. The first step takes advantage of the infiltration equationof Parlange et al. (1982), an quasi-exact solution of Richards’
equation to back out the values of S and 𝐾s from a collection
̃ data pairs using Bayesian inference withof measured (𝑡̃, 𝐼)
the DREAM(ZS) algorithm with sampling from past states(Vrugt, 2016). To this end, we present a robust, efficient, andexact numerical solution of Parlange’s infiltration equation,
which generalizes existing solutions of Latorre et al. (2015)
and Fernández-Gálveza et al. (2019) to warrant treatment andestimation of all its four parameters: S, 𝐾s , 𝐾i , and β. In a second step, the maximum a posteriori density values of S and𝐾s are admitted to Philip’s two-term infiltration equation toyield estimates of c and 𝑡valid . We benchmark, test, and evaluate our approach using synthetic infiltration data simulated by
HYDRUS-1D (Šimůnek et al., 2008). Our approach has someelements in common with Dohnal et al. (2010), who used theinfiltration equation of Haverkamp et al. (1994) with β = 0.6
to estimate S and 𝐾s from synthetic infiltration data of a threedimensional axisymmetric soil model. The optimized valuesof S and 𝐾s were entered in Philip’s two-term Equation 5 tosimulate one-dimensional infiltration.
The remainder of the paper is organized as follows. Section 2 introduces the cumulative infiltration data that is used totest, benchmark, and evaluate our two-step procedure. In Section 3, we detail the different building blocks of our methodology. This is followed in Section 4 with a presentation anddiscussion of our results. In Section 5, we briefly discuss
the main implications of this paper and present a simplified
implementation of our methodology. Finally, Section 6 concludes this paper with a summary of the main findings.

2

SYNTHETIC INFILTRATION DATA

To benchmark, test, and evaluate our method, we are in needof cumulative infiltration data and corresponding estimates ofthe soil sorptivity, S, and the saturated hydraulic conductivity, 𝐾s . These two soil properties should be known beforehand, otherwise we cannot verify that our estimates of S and𝐾s are unbiased and accurate. Practical experience suggeststhat it is not particularly easy to accurately determine the soil
sorptivity. This necessitates infiltration measurements undergravity-free conditions when water flow is largely controlled
by capillary adsorption or desorption. This requires considerable experimental effort, particularly if we wish to evaluate
our method for many different soil types. What is more, our

JAISWAL ET AL.

analysis demands a high temporal measurement resolution toaccurately pinpoint the time validity of Philip’s two-term infiltration equation. Therefore, we resort to synthetic infiltration
data in this paper, as this allows control on data quality andquantity. We did evaluate our method for measured data as
well using the more than 1,300 infiltration experiments documented in the Soil Water Infiltration Global (SWIG) database(Rahmati, Weihermüller, Vanderborght, et al., 2018; Rahmati,
Weihermüller, & Vereecken, 2018). Only a small number ofthese experiments have a data quality and quantity sufficientenough to warrant application of our methodology. These
results (not shown here) corroborate our findings for the synthetic infiltration data.
Vertical infiltration into a 200-cm-deep homogeneous
soil was simulated with HYDRUS-1D (Šimůnek et al.,2008, 2016) using the hydraulic functions of Mualem–van
Genuchten (MVG) (Mualem, 1976; van Genuchten, 1980).
The soil is initialized as dry by setting the hydraulic head of
the column equal to −15,000 cm (except for sand, we used
−1,000 cm). Cumulative infiltration is then simulated for aperiod of 240 h using a constant pressure head at the surface
and free drainage condition at the bottom of the soil profile(see also Rahmati et al., 2020). Table 1 lists the MVG parameter values of the HYDRUS-1D catalog that were used to simulate the infiltration data for each soil type.
As soil column discretization may affect the numericalresults, we used HYDRUS-1D to simulate horizontal infiltration for an ensemble of discretized profiles. The simulated infiltration was compared against an analytic solution for
infiltration without gravity, which is based on the Boltzmanntransform. For each discretized profile, the simulated watercontent distribution in the soil column was plotted against the
Boltzmann variable, λ, in Equation 2. A discretized profile
is deemed accurate if the (θ, λ) relationships at many different simulation times coalesce in one single curve. The highest simulation accuracy was achieved with a profile composedof 401 nodes using an entrance node of 10−6 cm and gradually increasing size of the elements in the profile up to 1 cm
for the last node. This profile was subsequently used to simulate vertical infiltration. To negate numerical errors and promote accuracy of the simulated vertical infiltration, the internal interpolation tables of (θ, ℎ) and (ℎ, 𝐾) were disabled,
and the default initial time step of HYDRUS-1D was adjusted
to satisfy convergence criteria for all different USDA soils.
What is more, a modified MVG model with air-entry valueof −2 cm was used for soils with 𝑛VG < 1.2 to avoid unrealistically large changes in the hydraulic conductivity near saturation (Schaap & van Genuchten, 2006; Vogel & Cislerova,
1988; Vogel et al., 2000).
Before we can analyze the vertical infiltration experimentsof HYDRUS-1D, we must first post-process the simulated
output (between 1,200 and 13,000 print times) so as to yield
a common data length and format for each soil. First, we0.380

0.430
0.380

0.045

0.170

0.111

0.066

0.090

0.104

0.266

0.197

Sand

Sandy clay

Sandy clay loam

Sandy loam

Silt

Silt loam

Silt clay

Silty clay loam

0.430

cm3 cm−3

0.089

0.070

0.067

0.034

0.065

0.100

0.100

0.045

0.057

0.078

0.095

0.068

𝛉𝐫

0.010

0.005

0.020

0.016

0.075

0.059

0.027

0.145

0.124

0.036

0.019

0.008

cm−1

𝛂

1.23

1.09

1.41

1.37

1.89

1.48

1.23

2.68

2.28

1.56

1.1

1.09

𝒏𝐕𝐆

0.52

0.35

1.65

1.34

3.3

1.60

0.78

9.21

6.20

2.19

1.45

1.02

cm h−1/2

𝑺̃

0.070

0.020

0.450

0.250

4.421

1.310

0.120

29.70

14.59

1.040

0.260

0.200

cm h−1𝑲̃ 𝐬

̃ which are of main interest in the present study. The tilde operator is used for both variables to signify “measured” quantities. The soil
Note. The last two columns list the saturated soil hydraulic conductivity, 𝐾̃ s , and sorptivity, 𝑆,
̃ is computed with Equation 2 using HYDRUS-1D simulated moisture profiles of a horizontal infiltration experiment with initial water content, θi . θs , saturated water content; θr , residual water content; α and 𝑛VG , vansorptivity, 𝑆,
Genuchten parameters.

0.430

0.360

0.450

0.460

0.410

0.390

0.410

0.088

0.057

0.410

Loamy sandClay loam

𝛉𝐬

Loam

0.271

0.150

Clay

𝛉𝐢

Values of the van Genuchten (1980) soil hydraulic parameters for the 12 USDA textural classes derived from Carsel and Parrish (1988)

Parameters

TA B L E 1

JAISWAL ET AL.
Vadose Zone Journal

5 of 27

6 of 27

JAISWAL ET AL.

Vadose Zone Journal

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

̃ data for a 200-cm homogeneous soil profile using the Mualem–van Genuchten
F I G U R E 1 HYDRUS-1D simulated cumulative infiltration, 𝐼,(MVG) hydraulic functions with parameter values of the 12 different textural classes documented in Table 1. A constant pressure head was assumed
at the soil surface, and a free drainage condition was imposed at the bottom boundary of the profile. The HYDRUS-1D output is post-processed tolimit the cumulative infiltration for each soil type to 5 cmsynchronize the print times of the datasets. To this end, we
discretize the 240-h simulation period into n = 100 equally
spaced print times on a logarithmic (base 10) scale. The cumulative infiltration at each print time (query point) is then computed from the raw HYDRUS-1D output using linear interpolation. Logarithmic scaling of the print times guarantees ahigh measurement resolution at the beginning of the infiltration experiment. This is deemed necessary to portray accurately the relatively high infiltration rates early on in the experiment. Figure 1 presents the post-processed infiltration curvesof each USDA soil type using 100 logarithmically scaled measurement times. The simulation period of 240 h is too large to
notice subtle differences in the functional shape of the infiltration curves at early times. Nevertheless, it is evident that thesoils differ substantially in the amount of cumulative infiltration. For example, sand (Figure 1E) has infiltrated more than
7,000 cm of water, whereas the clay soil (Figure 1A) barely
reaches 50 cm in the same time period of 240 h. In theory,
one would expect clay to infiltrate the least amount of water.
However, the sandy clay (Figure 1F), silty clay (Figure 1L),and silty clay loam (Figure 1L) soils infiltrate only 32.2, 7.2,
and 17.6 cm of water, respectively, over the 240-h period. To
understand this finding, we must look at the hydraulic properties of the different soils. According to Table 1, the sandy clay,
silty clay, and silty clay loam soils exhibit the lowest sorptivity, S, and saturated hydraulic conductivity, 𝐾s , of all 12 soils.

Thus, the MVG parameter values, as documented in Carsel
and Parrish (1988) and used in the HYDRUS-1D catalog, may
not portray accurately the hydraulic properties of each texturalclass. It should suffice to say that the post-processed infiltration curves are a nearly perfect match with the raw output of
HYDRUS-1D. Their distance is negligible small.
A last comment about the numerical accuracy of the simulated infiltration data is in order. Our experience suggests,
that it is not always easy to simulate accurately infiltration
experiments. For example, a careful examination of the simulated data of silty clay loam demonstrates a minor dip in the
cumulative infiltration at 𝑡 = 30 h. This artifact is small andmay not catch everyone’s attention; nevertheless, the subtle,
unexpected variations in the infiltration rate may explain inpart why silty clay loam is classified as one of our two dissonant soils. As we have devoted special care and attention
to numerical errors, we are confident that the current dataset
is robust and accurate for all soils, with possible exception of
sandy clay loam and silty clay loam.
The logarithmically scaled infiltration curves use synchronized measurement times yet correspond to an infiltration
event of 240 h. For most soils, this duration is excessive andfar exceeds the time at which a constant infiltration rate isreached. What is more, in field experiments, the cumulative
infiltration is typically limited, as the soil underneath may not
be homogeneous over a very large depth. Hence, we should

JAISWAL ET AL.

adjust the duration of the experiment so as to have a reasonable amount of infiltration for each soil type. For example,
if we assume a porosity ϕ = 0.4 and initial moisture content
θi = 0.1, then about 3 cm of cumulative infiltration would suffice to bring the top 10 cm of the soil to saturation. Continued
infiltration is not expected to yield more information about the
topsoil’s hydraulic properties.
For reasons detailed above, we do not proceed with thisfirst dataset but create a second dataset from the raw output
of HYDRUS-1D wherein we limit the cumulative infiltration
to 5 cm. In doing so, we divide the raw output of HYDRUS1D into equally spaced intervals of 0.05 cm. This results in
𝑛 = 100 measurements of the cumulative infiltration. The corresponding measurement times are derived via linear interpolation using the raw HYDRUS-1D output. Thus, the measurement times in this second dataset correspond to cumulative
infiltration values of 0.05, 0.10, 0.15, . . . , 5.00 cm. The soobtained measurement times differ substantially among thedifferent soil types. Sand, for instance, only needs a handful of minutes to infiltrate 5 cm of water, whereas the silty
clay soil requires about 150 h to do so. This second dataset isused throughout the remainder of this paper. In anticipation ofexperimental data, we refer to this second dataset as measured
data or measured infiltration data. This fixes the notation.
Now we have discussed the vertical infiltration data simulated with HYDRUS-1D, we are left with estimation of the“true” sorptivity, S, for each soil type. Unlike the saturated
hydraulic conductivity, 𝐾s , the sorptivity, S, is not invariant
but dependent on the soil’s initial moisture content. Hence,
S is not an innate property of the MVG model but a byproduct of the hydraulic parameters, θs , θr , α, 𝑛VG , and 𝐾slisted in Table 1 and the uniform initial moisture content,
θi , of the 200-cm-deep soil column. We can express 𝑆 =𝑓 (θs , θr , α, 𝑛VG , 𝐾s , θi ) in closed form (e.g., Moret-Fernández
et al., 2017), but experience suggests that this relationshipdoes not always yield accurate estimates of the soil sorptivity. Therefore, we resort to the definition of the sorptivity inEquation 2 and use numerical integration of the HYDRUS-1D
simulated λ(θ) curves for horizontal infiltration (no gravity).
The penultimate column of Table 1 documents the values ofthe sorptivity for all our twelve soil types. This concludes thedescription of the infiltration dataset.
To provide more insights into the hydraulic properties ofthe soils in our data set, please consider Figure 2, which displays the textural triangle with color coding of each texturalclass matching (a) the sorptivity, and (b) saturated hydraulicconductivity. Thus, the two triangles map the listed values of
S and 𝐾s in Table 1 to their respective textural class. It isevident that soil texture exerts a large control on the sorptivity and saturated hydraulic conductivity. In general, thesorptivity decreases with decreasing sand fraction. We noteagain that the sorptivities of the silty clay and silty clay loam
soils are smaller than their counterpart of the clay soil. This

Vadose Zone Journal

7 of 27

anomaly somewhat deteriorates the strength of the relationship between the sorptivity, S, and texture of the soil.
Qualitatively similar findings are observed for the saturated
hydraulic conductivity of the USDA collection of soils. The
color variations of the textural triangle are less pronounced.
The saturated hydraulic conductivities of the silty clay, sandy
clay, and silty clay loam soils appear rather low in comparison with the documented 𝐾s value of the clay sample. Consequently, these three soils cannot reach the infiltration rate of
the clay soil at late times, as is clearly evident from the shapeof the cumulative infiltration curves in our dataset. From the
textural triangles, we conclude that the silty clay and silty clay
loam soils exhibit the lowest values of the sorptivity, S, andsaturated hydraulic conductivity, 𝐾s . This explains why fromall USDA soils, these two soils infiltrate the least amount ofwater.
Before we proceed to the section below, we store the infiltration dataset, {𝑡̃𝑖 , 𝐼̃𝑖 }𝑛𝑖=1 , of each soil type in two separate vectors, namely, a 𝑛 × 1 vector, ̃𝐈 = [𝐼̃1 𝐼̃2 … 𝐼̃𝑛 ], of

cumulative infiltration measurements, and a 𝑛 × 1 vectorof associated print (measurement) times, ̃𝐭 = [𝑡̃1 𝑡̃2 … 𝑡̃𝑛 ].
The corresponding values of the sorptivity and saturated
hydraulic conductivity are referred to as 𝑆̃ and 𝐾̃ s , respectively. Although we use simulated data in the present analysis,
the tilde notation, ∼, emphasizes use of observed quantities asprecursor to measured infiltration data from laboratory and/or
field experiments.

3

THEORYIn this section, we present the different building blocks ofour methodology for assessing the time validity, 𝑡valid [L],of Philip’s two-term infiltration Equation 5 from measured
cumulative infiltration data. As part of this analysis, we also
yield estimates of the soil sorptivity, S [L T1/2 ], saturated
hydraulic conductivity, 𝐾s [L T−1 ], and dimensionless coefficient, c. As by-product, we also infer the marginal distributionof coefficient, β, in the infiltration equation of Parlange et al.
(1982).

3.1

The infiltration equation of Parlange

The infiltration equation of Parlange et al. (1982) serves as
principal foundation of our methodology. This quasi-exact
analytic solution of the water content form of Richards’ equation simulates cumulative vertical infiltration, I(t) [L], into a
homogeneous soil at uniform initial moisture content:(𝐾 −𝐾 )[𝐼(𝑡)−𝐾i 𝑡](𝐾s −𝐾i )2(1 − β)𝑡 = s i 𝑆 2
𝑆2
}
)
(
{
2β(𝐾s −𝐾i )[𝐼(𝑡)−𝐾i 𝑡]
β−1+− 12 log 1β exp
2
𝑆
β

(7)

8 of 27

JAISWAL ET AL.

Vadose Zone Journal

̃ in units of cm h−1/2 and (b)F I G U R E 2 Textural triangle with the 12 USDA soil types and associated measurements of the (a) sorptivity, 𝑆,−1saturated hydraulic conductivity, 𝐾s , in cm h , using a parula (blue–green–orange–yellow) color map. The colorbar at the right-hand side of both
triangles assigns values to the parula color palette

where t > 0 [T] is the elapsed time since the start of the infiltration event and β ∈ (0, 1) is a dimensionless integral shape
constant (Haverkamp et al., 1990). Equation 7 was generalized to nonponded conditions by Parlange et al. (1985) andHaverkamp et al. (1990) with an expression in explicit form
presented by Barry et al. (1995). The coefficient β can becomputed as follows (Fuentes et al., 1992; Haverkamp et al.,1994):
θ

β=2−2

∫θ s
i

[

𝐾(θ)−𝐾i
𝐾s −𝐾iθ](θs −θiθ−θi)𝐷w (θ)dθ∫θ s 𝐷w (θ)dθ(8)

iwhere θ [L3 L−3 ] denotes the volumetric soil moisture content at time 𝑡 > 0, 𝐾(θ) [L T−1 ] and 𝐷w (θ) [L2 T−1 ] signify the unsaturated soil hydraulic conductivity and soil water
diffusivity functions, respectively, and θi ≤ θ ≤ θs . Notwithstanding this theoretical definition for β, its value is difficult
to determine from measured infiltration data, as its impact
on the 𝐼(𝑡) relationship of Equation 7 is visible only at lateinfiltration times, well beyond the length of most infiltration experiments. It is common practice, therefore, to treat β
as a known constant. For example, Haverkamp et al. (1994)
obtained β = 0.563 from application of equation (8), whereas
a value of β = 0.6 is commonly used in the literature (Dohnal
et al., 2010; Latorre et al., 2015, 2018; Moret-Fernández et al.,2020; Rahmati et al., 2019). We treat β as an unknown quantity.
Equation 7 has several desirable properties for infiltration
modeling and soil hydraulic parameter estimation. First andforemost, Parlange’s equation is valid for all infiltration times.

Secondly, as Haverkamp et al. (1994) has shown for a threedimensional infiltration experiment, the equation is in close
agreement with measured data. Lastly, the parameters S, 𝐾s ,and β have a sound physical and/or mathematical underpinning, adjustable to the initial and boundary conditions of theinfiltration event (Haverkamp et al., 1994; Rahmati et al.,2020). Despite these desirable qualities, Parlange’s infiltration
Equation 7 has not entered into mainstream use for infiltration modeling and data analysis in lieu of Philip’s two-term
Equation 5. As Parlange’s equation is an implicit solution of
the infiltration process, some mathematical skill is required tonumerically solve Equation 7 in pursuit of the I(t) relationship.
To simplify practical application of Equation 7, Haverkamp
et al. (1994) derived three separate expressions for cumulative
infiltration at very short (VS), short (S), and long (L) times,
respectively:√
𝐼 VS (𝑡) = 𝑆 𝑡(9a)]√ [( 2 − β )𝐼 S (𝑡) = 𝑆 𝑡 +(𝐾s − 𝐾i ) + 𝐾i 𝑡
3

(9b)

1
𝑆2
log (β)2(1 − β) 𝐾s − 𝐾i(9c)

𝐼 L (𝑡) = 𝐾s 𝑡 −These analytic expressions are much easier to evaluate yet
trade difficulties with the numerical solution of Equation 7
with problems about how to establish exact time validities of

JAISWAL ET AL.

Equations 9a–9c and guarantee convergence of the simulatedI(t) relationship. Lassabatere et al. (2009) developed time limits for each of the equations and their dependency on soil texture. These limits are only an approximation of the true time
validity and require knowledge of the soil texture and valuesof S and 𝐾s . This undermines practical use. Note that the firstterm of Equations 9a and 9b is synonymous with the capillaryterm of Philip’s two-term infiltration Equation 5. The second
term at the right-hand side of Equation 9b must, therefore,
equal Philip’s gravity term. Thus, for very short to short infiltration times, we may write[(
𝑐𝐾s 𝑡 =

9 of 27

Vadose Zone Journal

2−β
3

3.2

]

)

(𝐾s − 𝐾i ) + 𝐾i 𝑡

(10)

For 𝐾i = 0, the expression above simplifies to the following
relationship (Moret-Fernández et al., 2017):
𝑐=

2 1
− β
3 3

regression. The anticipated difficulties with the inference ofthe coefficient β prompt use of a Bayesian method coupledwith the DREAM(ZS) algorithm (Vrugt et al., 2008, 2009;
Vrugt, 2016). Such methodology not only returns estimatesof the optimum parameter values but simultaneously also
assesses their associated uncertainty. Uncertainty quantification is key, as deterministic estimates are difficult to interpretby themselves without knowledge of their underlying confidence intervals.

(11)

between coefficients β and c in the infiltration equations of
Parlange et al. (1982) and Philip (1957) (two-term variant),
respectively.
In the Supplemental Information S1, we present theory
and implementation of an exact and robust numerical solution of Parlange’s infiltration equation. This solution considers time, 𝑡 ≥ 0 (h), as the independent (explanatory) variable
and cumulative infiltration, 𝐼 ≥ 0 (cm), as the dependent variable. Our solution builds on the residual form of Equation 7
and uses root finding with the secant method to iteratively findits zero-points. Unlike Newton’s method, this procedure does
not demand a closed-form expression for the derivative function but uses secant lines to converge the sequence of iterates
to the zero-point of Parlange’s residual function. Algorithm
S1.1 presents a detailed step-by-step recipe of our numericalsolution of Equation 7, and Appendix A presents its implementation in MATLAB. This subroutine called Parlange may
be executed from the MATLAB command prompt as I = Parlange(eta,t) and is synonymous to a vector-valued function,
𝐈 = 𝐟 (𝛈, ̃𝐭), of the Parlange parameters, 𝛈 = [𝑆 𝐾s β 𝐾i ]⊤ ,and n measurement times, ̃𝐭 = [𝑡̃1 𝑡̃2 … 𝑡̃𝑛 ]⊤ . The function
returns as output argument the 𝑛 × 1 vector of cumulative
infiltration values, 𝐈 = [𝐼1 𝐼2 … 𝐼𝑛 ]⊤ , of Equation 7. Benchmark experiments (not shown herein) demonstrate that the
MATLAB implementation is exact, robust, and computationally efficient requiring only a few milliseconds to complete a240-h infiltration experiment on a modern desktop computerwith print step of 0.1 h.
Parlange’s infiltration Equation 7 not only demands an iterative method to solve for its I(t) relationship but also requires
use of a search method to estimate its parameters from measured cumulative infiltration data. This amounts to nonlinear

Step 1: Bayesian inverse modeling

We estimate the parameters in Parlange’s infiltration equation
using Bayesian inference with the DREAM(ZS) algorithm.
Bayesian inference allows for an exact description of parameter uncertainty (and other sources of uncertainty) by treating
the parameters, η, of Equation 7 as probabilistic variables withjoint posterior probability density function, 𝑝(η|̃𝐈). This multivariate distribution, the so-called posterior parameter distribution, is the consequence of two antecedents, a prior distribution, 𝑝(η), which captures our initial degree of beliefs inthe values of the model parameters, and a likelihood function,𝐿(η|̃𝐈), which quantifies per the rules of probability theory
the level of confidence (= conditional belief) in the parameter values, η, in light of the cumulative infiltration data,̃𝐈 = [𝐼̃1 𝐼̃2 ⋯ 𝐼̃𝑛 ]⊤ , alone. Bayes’ theorem expresses mathematically the fundamental relationship between the prior, conditional, and posterior (= updated) beliefs of the parameters.
In the present case, this results in the following formulation:𝑝(η)𝐿(η|̃𝐈)∝ 𝑝(η)𝐿(η|̃𝐈)𝑝(η|̃𝐈) =
𝑝(̃𝐈)

(12)

where the denominator, 𝑝(̃𝐈), the so-called evidence or
marginal likelihood, acts as a normalizing constant𝑝(̃𝐈) =

∫𝐄𝑝(η)𝑝(̃𝐈|η)dη =∫𝐄

𝑝(η)𝐿(η|̃𝐈)dη =∫𝐄𝑝(η|̃𝐈)dη(13)
so that the posterior distribution, 𝑝(η|̃𝐈), integrates to unity
over the prior (feasible) parameter space, η ∈ 𝐄 ⊆ ℝ𝑑 , where
d signifies parameter dimensionality. Knowledge of 𝑝(̃𝐈) is
strictly necessary for hypothesis testing to select the mostplausible infiltration function from a suite of different candidate models deemed valid a priori. If, however, we rely on asingle hypothesis, then the denominator, 𝑝(̃𝐈), in Equation 12
is of no particular interest, as all statistical inferences about
the parameters of Parlange’s infiltration Equation 7 can bemade from the unnormalized posterior distribution, 𝑝(η|̃𝐈) ∝𝑝(η)𝐿(η|̃𝐈). This necessitates definition of a prior distribution,𝑝(η), for the Parlange model parameters, a mathematical form

10 of 27

JAISWAL ET AL.

Vadose Zone Journalof the likelihood function, 𝐿(η|̃𝐈), and a method that can infer
the posterior distribution, 𝑝(η|̃𝐈).
Before we move on to our choice of prior distribution,
we postulate that the hydraulic conductivity, 𝐾i , at the initial moisture content, θi , is equal to zero. This assumption
is certainly justifiable for our more coarse-textured soils, as
their initial moisture contents are relatively low and close toresidual saturation (see Table 1). For the more fine-textured
soils, the assumption that 𝐾i = 0 is questionable; nevertheless, numerical experiments have shown this has a negligible
effect on their posterior estimates of S, 𝐾s and β. With 𝐾i = 0,
Equation 7 simplifies to
(

𝐾s
𝑆

)2

𝐾 𝐼(𝑡) 1(1 − β)𝑡 = s 2 − log
2
𝑆

{

[
}
]
2β𝐾s 𝐼(𝑡)
β−1
1
exp
+
β
β
𝑆2

(14)
which leaves us with three unknown parameters, namely,η = [𝑆 𝐾s β]⊤ . Note that this simplification is not given bylimitations of our parameter estimation algorithm, as theDREAM(ZS) algorithm can handle a very large number ofparameters (Laloy et al., 2015; Linde & Vrugt, 2013; Qinet al., 2018; Zhang et al., 2020). Indeed, a comparison of themarginal prior and posterior distribution of each parameter
will convey how much information is available in the calibration data to constrain each parameter. Knowledge of the initialmoisture content is important when interpreting the inferredvalues of S.

3.2.1

Prior distribution

The prior distribution should encode all the “subjective” knowledge about the Parlange model parameters, η =[𝑆 𝐾s β]⊤ , before collection of the measured cumulative
infiltration curve, ̃𝐈. This distribution, often simply called theprior, expresses one’s beliefs about the parameters before thedata (also referred to as evidence) are taken into account. In
the absence of detailed prior information about the values ofS, 𝐾s , and β, we treat all parameter values as equally likely
a priori and assign a multivariate uniform prior distribution
to the parameter space demarcated by the upper and lower
bounds listed in Table 2. This 𝑑 = 3-dimensional hypercube
is synonymous to the feasible parameter space, 𝐄, and implies
a prior density, 𝑝(η) = 5 × 10−4 , so that ∫𝐄 𝑝(η) = 1.
The lower bounds of S, 𝐾s , and β are dictated by physicallimits and must be set to zero. Their upper limits, however, are
not unambiguous and warrant some discussion. The formal
definition of parameter β in Equation 8 implies that 0 < β < 2.
This exceeds the interval of β ∈ (0, 1) given by Parlange et al.
(1982) but has been used in other recent studies. Note thatβ ≠ 1, otherwise the left- and right-hand sides of Equation 7
equal zero. The upper bounds of the soil sorptivity and saturated hydraulic conductivity are not as clearly defined andare far more diffuse. We specify rather liberal upper values ofT A B L E 2 Description of the parameters of Parlange’s infiltration
Equation 14 including symbol, units, and lower and upper bounds

Parameter Description

Unit

Min.

Max.

S

Soil sorptivity

cm h−1/2

0

20

𝐾s

Saturatedhydraulicconductivity

cm h−1050βDimensionlesscoefficient

–

0

2

20 cm h−0.5 and 50 cm h−1 for S and 𝐾s , respectively, so as toguarantee that our uniform prior, 𝑝(η), accommodates a diversity of soils. If, for some reason, the upper limits of S and 𝐾s
truncate the posterior distribution, 𝑝(η|̃𝐈), then these bounds
can always be enlarged.
Note that one can exploit the SWIG database and construct
an informative prior for S, 𝐾s , and β for each soil type. Suchinformative (nonuniform) prior may simplify the inferenceof β, particularly for infiltration experiments with a limited
duration. This warrants further research. The consequence ofan uninformative prior is that the posterior parameter distribution, 𝑝(η|̃𝐈), will depend only on the cumulative infiltration data, ̃𝐈 = [ 𝐼̃1 𝐼̃2 ⋯ 𝐼̃𝑛 ]⊤ , through the likelihood function, 𝐿(η|̃𝐈). Indeed, with a uniform prior, the posterior density, 𝑝(η|̃𝐈), of the Parlange parameters, η = [𝑆 𝐾s β]⊤ , will
simply equal a constant multiple of 5 × 10−4 of the likelihood, 𝐿(η|̃𝐈), for all η ∈ 𝐄 ⊆ ℝ3 . Thus, with an uninformative prior distribution, Bayesian inference reduces to maximum likelihood estimation, but with the important distinctionthat our framework infers not only the peak of 𝑝(η|̃𝐈), but also
its underlying distribution.

3.2.2

Likelihood function

We are left with the definition of the likelihood function,𝐿(η|̃𝐈), in Equation 12. The likelihood was designated asmathematical quantity by Sir Ronald Fisher (1934) to measure our degree of belief (confidence) in simulated outcomes.
The likelihood principle states that all relevant information
about the parameters, η, is contained in the likelihood function
for the observed data given the assumed (statistical) model.
This offers a self-contained framework for statistical modeling and inference and has become the foundation of modern
statistics and probability theory with application to frequentist
and Bayesian inference.
In our application, the likelihood (function) quantifies inprobabilistic terms the distance between the observed cumulative infiltration data, ̃𝐈, and its counterpart, 𝐈 = 𝑓 (η, ̃𝐭), simulated by Parlange’s infiltration Equation 14 using the parameter values, η = [𝑆 𝐾s β]⊤ , at the n measurement times, ̃𝐭 =[ 𝑡̃1 𝑡̃2 ⋯ 𝑡̃𝑛 ]⊤ . The word likelihood is often used loosely inJAISWAL ET AL.

Vadose Zone Journal

common speech as synonym for probability, yet probability
is inadequate as mathematical quantity to express our preferences among different outcomes. In practice, the use of
the word probability is appropriate when describing possible future outcomes for fixed parameter values before data
become available. The term likelihood, on the contrary, should
be used to express the level of confidence for a given outcome(parameter vector) after the infiltration data have become
available.
If we wish to quantify the likelihood of a simulated outcome, 𝐼1 , in light of an observation, 𝐼̃1 , we need to define
a probability density function, 𝑓𝐼̃1 , for 𝐼̃1 . In theory, we
could hypothesize any probability distribution for 𝐼̃1 , and as
a result, many potential formulations of the likelihood function, 𝐿(η|̃𝐈). In practice, however, we often limit our choice tocommon probability distributions that are easy to manipulate
analytically and evaluate numerically. We refer the interested
reader to Vrugt and Massoud (2019) for the derivation of asuite of commonly used likelihood functions. In the present
application, we make the common assumption that the measurement errors, ϵ = [ϵ1 ϵ2 ⋯ ϵ𝑛 ]⊤ , of the 𝑛 × 1 vector of
cumulative infiltration data, ̃𝐈, are independent (uncorrelated),zero-mean normally distributed with a constant variance, σ2̃ .
𝐼Thus, we assume that ϵ𝑗 ̃  (0,σ2̃ ) for all 𝑗 = (1, 2, … , 𝑛),𝐼

where  (𝑎, 𝑏) signifies the normal distribution with mean,𝑎, and variance, 𝑏 and symbol ∼ means distributed according to. Our assumptions lead to the following definition of thelikelihood function:𝐿(η|̃𝐈) ∝

{ 𝑛
∑[
𝑗=1

2 }− 2 𝑛
1𝐼̃𝑗 − 𝐼(η, 𝑡̃𝑗 )

]

(15)

where 𝐼̃𝑗 denotes the observed cumulative infiltration at thejth measurement time and 𝐼(η, 𝑡̃𝑗 ) signifies its simulated counterpart at time 𝑡̃𝑗 using Equation 14 with parameter values,η = [𝑆 𝐾s β]⊤ . The summation term between the brackets is
the well-known sum of squared residuals (SSR) used in leastsquares regression methods. Note that 𝐿(η|̃𝐈) is inversely proportional to the SSR; hence, parameter values that minimize
the SSR will also maximize the likelihood. As a result, themaximum likelihood solution of the Parlange parameters issynonymous with the least squares solution derived from nonlinear regression methods. Therefore, in the remainder of thispaper, the maximum a posteriori density (MAP) solution of𝑝(η|̃𝐈) is also referred to as least squares solution.

3.2.3
The DiffeRential Evolution Adaptive
Metropolis (DREAM) algorithm
Now that the prior distribution and likelihood function havebeen defined, what is left is to summarize the posterior dis-

11 of 27

tribution, 𝑝(η|̃𝐈), of the Parlange model parameters in Equation 14. As this model is nonlinear in its parameters, theposterior distribution 𝑝(η|̃𝐈) cannot be obtained by analytical means nor by analytical approximation. We must therefore resort to iterative methods that approximate the posterior distribution by generating a large sample from this distribution. In this paper, we approximate the target distribution, 𝑝(η|̃𝐈), using Markov Chain Monte Carlo (MCMC) simulation with the DREAM(ZS) algorithm (Vrugt, 2016; Vrugt
et al., 2008, 2009). This multichain MCMC method has foundwidespread application and use in many different fields ofstudy, and it exhibits excellent sampling efficiencies on complex, high-dimensional, and multimodal target distributions(Laloy et al., 2015; Qin et al., 2018; Volpi et al., 2017;, Zhang
et al., 2020; to name a few). A detailed description of thismethod and its different extensions appears in Vrugt (2016),
and interested readers are referred to this paper for more
information.
For each soil, we approximate the posterior distribution,𝑝(η|̃𝐈), of the Parlange parameters, η = [𝑆 𝐾s β]⊤ , using 𝑁 =3 Markov chains and default settings of the algorithmic variables of the DREAM(ZS) algorithm. In doing so, we use afixed computational budget of 60,000 Parlange function evaluations and create 20,000 samples in each chain. This number of samples is much larger than required by convergencecriteria such as the multivariate scale-reduction factor, 𝑅̂ d ,of Brooks and Gelman (1998), which declared convergence,𝑅̂ d ≤ 1.2, after just a few thousand chain samples. Hence, we
can safely discard the first half of each chain (50% burn-in)
and use the last 10,000 samples in each chain to summarize
the posterior distribution of the Parlange parameters, S, 𝐾s ,and β for each different soil. This amounts to 30,000 posterior
parameter vectors, η, and requires about 2 min to complete
on a modern desktop computer. The MAP or least squares
solution is readily found among these posterior samples by
locating the parameter vector with highest posterior density,𝑝(η|̃𝐈).

3.3
Step 2: Estimation of c and 𝒕𝐯𝐚𝐥𝐢𝐝 in
Philip’s two-term equation
The trivariate posterior distribution of the soil sorptivity S,
saturated hydraulic conductivity 𝐾s , and coefficient β of Parlange’s infiltration equation serves as input to our second step.
Specifically, among the samples of the 𝑁 = 3 Markov chains
simulated by the DREAM(ZS) algorithm, we locate the solution that maximizes the posterior density, 𝑝(η|̃𝐈). Due to our
choice of prior distribution and likelihood function, the corresponding parameter values are maximum likelihood as wellas least squares estimates. We use this wording interchangê ⊤ , ofably. The least squares parameter values, ̂η = [𝑆̂ 𝐾̂ s β]Parlange’s infiltration Equation 14 are hard-wired in Philip’s12 of 27

JAISWAL ET AL.

Vadose Zone Journal

two-term infiltration equation to yield√
𝐼(𝑡) = 𝑆̂ 𝑡 + 𝑐 𝐾̂ s 𝑡(16)This leaves us with two unknowns for each soil, namely, thedimensionless coefficient, c, and the time validity, 𝑡valid , of theabove expression. The least squares value, 𝑐,̂ of Equation 16
can be derived by analytic means. For the first m observations,̃𝐈𝑚 = [𝐼̃1 𝐼̃2 ⋯ 𝐼̃𝑚 ]⊤ , of the n-record of cumulative infiltration
measurements, ̃𝐈, we yield𝑐̂ = (𝐝⊤ 𝐝)−1 𝐝⊤ (̃𝐈𝑚 − 𝑆̂ ̃𝐭𝑚 )
1∕2

(17)

where the 𝑚 × 1 design vector, 𝐝 = [ 𝐾̂ s 𝑡̃1 𝐾̂ s 𝑡̃2 … 𝐾̂ s 𝑡̃𝑚 ]⊤ , the𝑚 × 1 vector ̃𝐭𝑚 = [ 𝑡̃1 𝑡̃2 ⋯ 𝑡̃𝑚 ]⊤ stores the measurement timesof ̃𝐈𝑚 and 1 ≤ 𝑚 ≤ 𝑛. Estimation of the time validity, 𝑡valid ,of Philip’s two-term infiltration equation now boils down
to selection of an optimal value for the integer, m. For𝑡 ≤ 𝑡̃𝑚 , Equation 17 should provide an adequate descriptionof the measured cumulative infiltration curve. For infiltration times, 𝑡 > 𝑡̃𝑚 , Philip’s two-term infiltration equation isdeemed invalid, and thus its goodness-of-fit is expected todeteriorate substantially.
The Bayesian information criterion (BIC) will help find an
optimal trade-off between the quality-of-fit of Equation 17
and the length, m, of the measured infiltration curve. This metric is commonly used for model selection among a finite setof competing hypotheses. The BIC encodes a natural preference for simpler models. This parsimony principle, known as
Occam’s razor, is traceable to the works of philosophers suchas Aristotle (384–322 BC) and Ptolemy (circa AD 90 to circa
AD 168) and consistent with requirements of falsifiability inthe scientific method of Popper (1992). The BIC is defined asfollows in our application:][
BIC = −2 log 𝐿(η|̃𝐈𝑚 ) + log(𝑚)(18)̂ ⊤ signifies the least squares values ofwhere η = [𝑆̂ 𝐾̂ s 𝑐]
the Parlange parameters, 𝐿(η|̃𝐈𝑚 ) denotes the likelihood ofPhilip’s two-term equation for subset, ̃𝐈𝑚 , and, again, m is thelength of the cumulative infiltration dataset (𝑚 ≤ 𝑛). Lower
values of the BIC are preferred.
We formalize the above ideas in Algorithm 1, which
presents a step-by-step recipe on our use of the BIC for deter̂mining the time validity, 𝑡valid , and least squares estimate, 𝑐,of the dimensionless constant, c, of Philip’s two-term infiltration (Equation 5). In words, for each length, m, of the subset, ̃𝐈𝑚 , of cumulative infiltration measurements, we resort toEquation 17 to compute the least squares value of the multiplicative coefficient, c. This value, 𝑐,̂ along with 𝑆̂ and̂
𝐾s from Parlange’s infiltration equation, enter Philip’s twoterm infiltration equation and are used to simulate cumula-

tive infiltration, 𝐈𝑚 , at measurement times, ̃𝐭𝑚 . Next, we calculate the likelihood, 𝐿(η|𝐈𝑚 ), of the simulated infiltration
curve, 𝐈𝑚 , using Equation 15 and compute the BIC. We execute this recipe for 3 ≤ 𝑚 ≤ 𝑛 and look for the value of mthat minimizes the BIC. This equals the optimal trade-off
between data length and quality of fit and thus is synonymous to the time validity of Philip’s two-term infiltration
equation.
Algorithm 1. Determination of the time validity, tvalid
Input:
Measured
cumulative
infiltration
curve,
̃𝐈 = [𝐼̃1 𝐼̃2 ⋯ 𝐼̃𝑛 ]⊤ at times, ̃𝐭 = [𝑡̃1 𝑡̃2 ⋯ 𝑡̃𝑛 ]⊤Least squares values of sorptivity, 𝑆̂ > 0, and saturated
hydraulic conductivity, 𝐾s > 0
Output: Time validity, 𝑡valid , of Equation 14 and corresponding least squares value, 𝑐,̂ of c1: begin
2: Set BIC1 = 1010 , BIC2 = 1010 , 𝑐̂1 = 0, and 𝑐̂2 = 0
3: for 𝑚 ← 3 to 𝑛 do
4:
Define cumulative infiltration subset,
̃𝐈𝑚 = [𝐼̃1 𝐼̃2 ⋯ 𝐼̃𝑚 ]⊤ , at ̃𝐭𝑚 = [𝑡̃1 𝑡̃2 ⋯ 𝑡̃𝑚 ]⊤5:
Define the 𝑚 × 1 design vector, 𝐝, of coefficientc, that is, 𝐝 = [𝐾̂ 𝑡 𝑡̃1 𝐾̂ 𝑡 𝑡̃2 ⋯ 𝐾̂ 𝑡 𝑡̃𝑚 ]⊤6:
Compute the least squares value, 𝑐̂𝑚 , of
1∕2
coefficient c, 𝑐̂𝑚 = (𝐝⊤ 𝐝)−1 𝐝⊤ (̃𝐈𝑚 − 𝑆̂ ̃𝐭𝑚 )7:
Evaluate Philip’s two-term Equation 16 at ̃𝐭𝑚using η = [𝑆̂ 𝐾̂ s 𝑐̂𝑚 ]⊤8:Compute the likelihood, 𝐿(η|̃𝐈𝑚 ), of the simulatedinfiltration data, 𝐈𝑚 , using Equation 15
9:Compute the Bayesian information criterion,
BIC = −2 log[𝐿(η|̃𝐈𝑚 )] + log(𝑚)10. end for
11. Find the index, 𝑚,̂ of the minimum BIC, thus,𝑚̂ = arg min1≤𝑚≤𝑛 BIC(𝑚)12. Define the time validity, 𝑡valid = 𝑡̃𝑚̂
13: return 𝑡valid and 𝑐̂𝑚̂
14: end

4

RESULTS AND DISCUSSION

In this section, we present the results of our framework
for the 12 different soil types considered in our analysis.
We first present the results of Step 1 of our methodology,
after which we move on to our findings for the second
step.

4.1

Results of Step 1: Parlange’s equation

Figure 3 displays time series plots of the HYDRUS-1D simulated cumulative infiltration data (red circles) of each soil type
and their simulated counterparts (solid blue line) of the infil-

JAISWAL ET AL.

Vadose Zone Journal

13 of 27F I G U R E 3 Comparison of observed (red dots) and simulated cumulative infiltration, I, data using the semi-implicit equation of Parlange et al.
(1982) for the 12 different USDA textural classes. Each dataset is composed of 𝑛 = 100 measurements and differs in duration as the cumulative
infiltration is limited to 5 cm

tration equation of Parlange with least squares value of S, 𝐾s ,and β derived from Bayesian analysis using the DREAM(ZS)
algorithm. The simulated cumulative infiltration is in excellent agreement with the measured data. For all soil types, thesolid blue line goes nearly perfectly through the measured
̃ data pairs. This result may not come as a complete sur(𝑡̃, 𝐼)
prise, as both the measured and simulated data equal a solution of Richards’ equation, albeit in numerical and analytic
form and with and without specification of the initial moisture content, respectively.
To provide insights into the corresponding estimates of the
Parlange parameters, please consider Figure 4, which presents
marginal posterior distributions of the soil sorptivity, S, in cm
h−1/2 , saturated hydraulic conductivity, 𝐾s , in cm h−1 , andunitless coefficient, β, for a representative group of five soils
including (a) clay, (b) clay loam, (c) sand, (d) silt loam, and(e) silty clay loam. The relative frequencies on the 𝑦 axis are
scaled between 0 and 1 to yield a common empirical density
for each parameter and soil type. The least squares value ofeach parameter is separately indicated in each graph using
a red cross. The Parlange parameters, S, 𝐾s , and β, appear
well defined by calibration against the measured cumulative
infiltration data. The histograms occupy only a small region
interior to their uniform prior distributions. This is a testament to the information content of the cumulative infiltration
data. The posterior marginal distributions of S, 𝐾s , and β are
generally well described by a normal distribution with meanthat coincides with their respective least squares estimates andsmall dispersion. The only exception is parameter β for silty
clay loam (Figure 4E3). Its histogram is truncated at two by
the prior distribution and, consequently, differs markedly fromnormality. A similar result is observed for sandy clay loam
(not shown). Our numerical solution of Parlange’s infiltration
equation detailed in Supplemental Section S1 and Appendix
A accommodates β values larger than 2. Thus, we temporarily lift the restriction that 0 < β < 2 and use an upper bound,β = 5, in our next trials with the DREAM(ZS) algorithm. The
use of such enlarged prior for β only changes the posterior
distribution of the sandy clay loam and silty clay loam soils.
The marginal distributions of their β values (not shown) now
become Gaussian as well (as all other histograms) with posterior means of about 3.5 and 4.5, respectively, and small
variance. The desire of the posterior β values of sandy clay
loam and silty clay loam to exceed 2 defies the mathematical–
physical definition of β in Equation 8. Therefore, we revokethat 0 < β ≤ 2 and classify sandy clay loam and silty clay
loam as dissonant soils in the remainder of this paper.
Table 3 documents summary statistics of the posterior
parameter distribution, 𝑝(η|̃𝐈), sampled with the DREAM(ZS)
algorithm. Specifically, for each soil type, we list the posteriormean, median, and standard deviation of S, 𝐾s , and β. This
table reiterates the main findings of Figure 3. The Parlange
parameters are well defined with posterior standard deviations that appear rather small. The total cumulative infiltra-

0.3530.518

Silt clay

Silty clay loam

0.518

0.353

1.703

1.381

3.907

1.698

0.799

9.326

6.317

2.267

1.491

1.041

Median

0.002

0.001

0.006

0.005

0.005

0.004

0.003

0.035

0.023

0.008

0.005

0.004

SD

Note. Parameter S is sorptivity, 𝐾s is saturated hydraulic conductivity, and β is a coefficient.

1.381

1.703

Silt loam

Sandy loam

Silt

1.698

3.907

Sandy clay loam

9.326

0.799

Sandy clay

Loamy sand

Sand

2.267

6.316

Loam

1.041

1.491

Clay loam—cm h−1/2

Mean

Clay

Soils

S

0.058

0.021

0.377

0.204

4.443

1.261

0.136

31.86615.260

0.969

0.306

0.213

—cm h−1

Mean

𝑲𝐬

0.058

0.021

0.378

0.205

4.443

1.261

0.136

31.87915.262

0.969

0.307

0.213

Median

0.001

0.000

0.007

0.006

0.011

0.002

0.001

0.242

0.078

0.009

0.006

0.002

SD

1.975

1.698

1.686

1.759

1.053

1.988

1.366

0.703

0.837

1.508

1.500

1.636

Mean

𝛃

1.984

1.698

1.687

1.760

1.053

1.992

1.367

0.704

0.837

1.508

1.500

1.636

Median

0.026

0.048

0.064

0.070

0.014

0.012

0.050

0.035

0.031

0.049

0.063

0.047

SD

T A B L E 3 Summary statistics of the DREAM(ZS) –derived posterior distribution, 𝑝(η|̃𝐈), of the parameters, 𝑥 = [𝑆 𝐾s β]⊤ , of Parlange’s infiltration equation for the twelve different soils
considered in our analysis

14 of 27
Vadose Zone Journal

JAISWAL ET AL.

JAISWAL ET AL.

Vadose Zone Journal

15 of 27F I G U R E 4 Histograms of the marginal posterior distribution of the parameters sorptivity (S), saturated hydraulic conductivity (𝐾s ), and βcoefficient of Parlange’s semi-implicit infiltration Equation 14 for a representative set of five soils, including (a) clay, (b) clay loam, (c) sand, (d) silt
loam, and (e) silty clay loam. The least squares values of the parameters are separately indicated in each graph with the red cross symboltion of 5 cm is more than enough for each soil type to warranta precise determination of its sorptivity, saturated hydraulicconductivity, and dimensionless coefficient, β, in Parlange’s
infiltration Equation 14. Whether the posterior estimates of
S and 𝐾s are accurate or not will soon be investigated. Notethat sandy clay loam suffers a similar problem with its β value
as silty clay loam, our other dissonant soil. Fortunately, as
we will demonstrate soon, the truncation of β at a value of
2 hardly affects the accuracy of the inferred S and 𝐾s values.
Next, we investigate the induced correlation structure
between the posterior samples of the Parlange parameters.
Figure 5 presents scatter plots of the bivariate posterior samples of (1) (𝑆, 𝐾s ), (2) (𝑆, β) and (3) (𝐾s , β) for the (a) clay,(b) clay loam, (c) sand, (d) silt loam, and (e) silty clay loam
soils. The dashed black line presents the fit of a linear regression model to the scattered samples, and the red cross portrays the least squares solution of each parameter pair. The
bivariate scatter plots of the sampled posterior parameter values demonstrate the presence of considerable parameter correlation. Indeed, for all but our two dissonant soils (sandy clay
loam not shown), the Parlange parameter pairs (𝑆, 𝐾s ), (𝑆, β),and (𝐾s , β) exhibit an almost linear relationship with strengthof the correlation that depends on the spread of the samples
around the regression line and results in correlation coefficients that range between .80 and .95 (see Table 4). Thus, if
a small increment of the sorptivity is met by a proportionate
positive perturbation to parameters 𝐾s or β, then the cumula-

T A B L E 4 Pairwise linear correlation coefficients between theparameters of Parlange’s infiltration equation for each of the soil types
considered herein
Parameter pairs
Soil type(𝑺, 𝑲𝐬 )(𝑺, 𝛃)(𝑲𝐬 , 𝛃)
0.932

Clay

0.826

0.971

Clay loam

0.858

0.959

0.966

Loam

0.827

0.964

0.943

Loamy sand

0.817

0.954

0.949

Sand

0.849

0.950

0.970

Sandy clay

0.837

0.964

0.949

Sandy clay loam

−0.663

0.553

0.120

Sandy loam

0.815

0.955

0.947

Silt

0.858

0.962

0.963

Silt loam

0.851

0.964

0.957

Silt clay

0.842

0.969

0.945

Silty clay loam

−0.671

0.485

0.208

Note. The correlation coefficients are derived from the posterior samples of the
DiffeRential Evolution Adaptive Metropolis (DREAM) algorithm Parameter S issorptivity, 𝐾s is saturated hydraulic conductivity, and β is a coefficient.

tive infiltration simulated by the Parlange infiltration equation
may not change much, resulting in an approximately similarvalue of the likelihood, 𝐿(η|̃𝐈). The truncated marginal distribution of β for the silty clay loam and sandy clay loam soils is16 of 27

JAISWAL ET AL.

Vadose Zone Journal

(a1)

(b1)

(c1)

(d1)

(e1)

(a2)

(b2)

(c2)

(d2)

(e2)

(a3)

(b3)

(d3)

(c3)

(e3)

F I G U R E 5 Scatter plots of the bivariate posterior samples of the parameter pairs, (1) (𝑆, 𝐾s ), (2) (𝑆, β) and (3) (𝐾s , β) of the Parlange
infiltration Equation 14 for a representative group of five soil types, including (a) clay, (b) clay loam, (c) sand, (d) silt loam, and (e) silty clay loam.
Parameter S is sorptivity, 𝐾s is saturated hydraulic conductivity, and β is a coefficient. The dashed black line characterizes the fit of a linear
regression function to the sampled data pairs. The red cross signifies the least squares solution of the Parlange parameters

responsible for an enhanced dispersion of the posterior samples along the regression line. This increased spread not only
reduces considerably the strength of the relationship between
the Parlange parameters, S, 𝐾s and β, but also changes to anegative value the correlation between S and 𝐾s (see Figure 4E1).
The posterior correlation structure among S, 𝐾s , and β ofParlange’s infiltration equation signals a warning. We should
not expect a method that can simultaneously determine thetime validity, 𝑡valid , and the parameters, S, 𝐾s , and c, ofPhilip’s two-term infiltration equation. No method and/or
framework would allow us to do this with confidence. If nothing else, we may be able to determine the product of c and𝐾s , but their individual values are ill defined in calibration
against cumulative infiltration data. This testifies to the needof our two-step inference method using Parlange’s infiltration
function as a vehicle for the soil sorptivity and saturated
hydraulic conductivity.
Finally, we conclude the results of our first step with a pairwise comparison of the observed and least squares (maximum
likelihood) estimates of the (a) soil sorptivity, S, and (b) saturated hydraulic conductivity, 𝐾s (see Figure 6). For completeness, Table 5 lists the maximum likelihood values of S and 𝐾s
̂ The maximum
including corresponding estimates of β.
likelihood estimates of 𝑆̂ and 𝐾̂ s derived from theDREAM(ZS) algorithm are in excellent agreement withT A B L E 5 Maximum likelihood estimates of the parameters,
sorptivity (S), saturated hydraulic conductivity (𝐾s ), and β coefficient,of Parlange’s infiltration equation derived from the joint samples of theMarkov chains of the DREAM(ZS) algorithm
𝛃̂

𝑺̂

𝑲̂ 𝐬

cm h−1/2

cm h−1

Clay

1.042

0.214

1.644

Soil type

Clay loam

1.491

0.307

1.505

Loam

2.267

0.969

1.507

Loamy sand

6.315

15.261

0.836

Sand

9.327

31.8770.704

Sandy clay

0.799

0.136

1.363

Sandy clay loam

1.700

1.261

2.000

Sandy loam

3.907

4.443

1.053

Silt

1.381

0.205

1.766

Silt loam

1.704

0.379

1.699

Silt clay

0.353

0.021

1.702

Silty clay loam

0.519

0.058

1.999

Note. The “measured” values of the soil sorptivity, S, and saturated hydraulic conductivity, 𝐾s , are listed in Table 1.

their measured counterparts listed in the last two columnsof Table 1. Indeed, the data pairs of the 12 soil types lie
almost exactly on the 1:1 line that is going from left to rightJAISWAL ET AL.

Vadose Zone Journal

17 of 27F I G U R E 6 Scatter plots of the observed, ∼, and estimated, ˆ, values of the (a) soil sorptivity, S, in cm h−1/2 and (b) saturated soil hydraulicconductivity, 𝐾s , in units of cm h−1 . The solid black line characterizes the 1:1 relationship between the measured and estimated quantities

across the graphs. Thus, the maximum likelihood parametervalues of Parlange’s infiltration equation do not only describe
accurately the measured infiltration data but also match theirobserved values. Setting 𝐾i = 0 in Parlange’s infiltration
equation does not seem to deteriorate the accuracy of theinferred sorptivities and saturated soil hydraulic conductivities. The β values listed in Table 4 of Lassabatere et al. (2009)
for sand, loam, silt, and silt loam are in good agreement with
our estimates for these soil types.
Now, the MAP values of the sorptivity, S, and saturated
hydraulic conductivity, 𝐾s , are known for each soil type, we
can exploit this knowledge in Step 2 of our methodology anddetermine the multiplicative coefficient, c, and time validity,𝑡valid , of Philip’s two-term infiltration Equation 5.

4.2
Results of Step 2: Time validity ofPhilip’s two-term equation
Figure 7 visualizes the results of Algorithm 1 and displays
traces of the BIC (in blue) and least squares value of the multiplicative coefficient, c (in green), of Philip’s two-term infiltration equation for all 12 USDA soil types. All graphs have theexact same vertical axes (in blue and green) with tick labels
listed only for the three figures in the leftmost column. The
solid red dot at the intersection of the dashed red line andthe 𝑥 axis equals the time validity, 𝑡valid , of Philip’s two-term
infiltration equation. The different traces of the BIC and 𝑐̂
highlight several important findings. In the first place, noticethe strong similarities in the BIC and 𝑐̂ traces of the different soil types. The BIC decreases rapidly at early times of
infiltration, then assumes rather similar values at intermediate times before increasing again towards the end of the infiltration experiment. The characteristic banana-shaped form ofeach BIC trace is the net result of two competing model evaluation criteria, namely, the length of the data set, m, in Equa-tion 18, and the quality of fit (likelihood) of Philip’s two-term
infiltration equation. The BIC will go down if a unit increment of the length of the cumulative infiltration dataset ismet by a sufficiently large increase in the likelihood, 𝐿(η|̃𝐈𝑚 ).
With exception of the loamy sand and sandy loam soils (seeFigures 7D, H), the least squares values of the multiplicative coefficient, c, drop almost instantly to values between
0.05 and 0.30 at very early infiltration times and then increaseagain with each length increment of the cumulative infiltration
dataset. Secondly, the BIC traces of the different soil types
exhibit a single well-defined minimum. This characteristic of
the BIC metric is crucially important in our search for an optimum model complexity and allows for an accurate determination of the time validity of Philip’s two-parameter infiltration equation. Indeed, the minimum BIC presents an optimal
trade-off between the length of the cumulative infiltration data
set and the quality of fit of Philip’s two-term infiltration equation.
In the third place, notice that the time validity of each soil
is well within the time required to infiltrate 5 cm of water.
Whereas silty clay and silty clay loam warrant an infiltration
experiment of about 2–3 d, most other soils need only a few
hours. Furthermore, the inferred time validities correlate withsoil texture. This may help develop practical guidelines onwhen to use Philip’s two-term equation. Lastly, the values ofcoefficient c at 𝑡 = 𝑡valid are significantly smaller than one andwell within the ranges reported in the literature. This inspires
further confidence in our methodology and will be discussed
below.
Table 6 reports the least squares values of coefficient, c, andcorresponding estimates of 𝑡valid for Philip’s two-term infiltration equation. The last column lists the characteristic time,𝑡char , of Equation 6 of Philip (1957). The optimum value ofcoefficient c ranges between 0.3 and 0.5 for soils with a relatively large sand fraction to values of 0.20–0.25 for soils witha high loam and clay content and values of about 0.15–0.20

18 of 27

JAISWAL ET AL.

Vadose Zone Journal

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

F I G U R E 7 Relationship between the Bayesian information criterion, BIC (in blue) and least squares value, 𝑐̂ (in green), of the multiplicative
coefficient, c, of Philip’s two-term infiltration Equation 5 as function of the length of the infiltration experiment for all 12 soil types considered in our
analysis. The graphs in Columns 2, 3, and 4 share the same vertical axis (left and right) as Graphs a, e, and i in the leftmost column. The verticaldashed red line pinpoints the minimum of the BIC. The solid red dot at the intersection of the red line and the x axis signifies the time validity ofPhilip’s two-term infiltration equation
T A B L E 6 Least squares value of the multiplicative coefficient, c,with corresponding estimates of the time validity, 𝑡valid , of Philip’s
two-term infiltration equation
Soil type𝒄̂

𝒕𝐯𝐚𝐥𝐢𝐝

Clay

0.210

4.909

23.745𝒕𝐜𝐡𝐚𝐫hClay loam

0.235

3.799

23.616

Loam

0.245

1.169

5.471

Loamy sand

0.458

0.099

0.171

Sand

0.498

0.055

0.086

Sandy clay

0.282

7.800

34.464

Sandy clay loam

0.207

0.795

1.818

Sandy loam

0.395

0.369

0.773

Silt

0.157

5.541

45.451

Silt loam

0.181

2.949

20.197

Silt clay

0.192

54.454288.160

Silty clay loam

0.178

31.42579.799Note. The last column lists the characteristic time, 𝑡char , derived from Philip (1957)
using Equation 6.

for silty soils. These values are not uncommon for the different
soil types. We can confirm these estimates by application of
our method to the soils of the SWIG database. The time validity of each soil type is much smaller than the characteristic

time, 𝑡char , of Philip (1957) and the gravity time, 𝑡grav , of thecharacteristic time method of Rahmati et al. (2020). This reiterates the importance of an assumption-free and data-driven
methodology for assessing the time validities of infiltration
models. Nevertheless, 𝑡valid and 𝑡char (and for that matter 𝑡grav )exhibit a strong linear relationship with correlation coefficientthat exceeds .95.
We have focused our attention on an assessment of themultiplicative coefficient, c, and the time validity, 𝑡valid , ofPhilip’s two-term equation without recourse to an assessment of the quality of fit of this infiltration equation.
Figure 8 compares measured (red dots) and simulated infiltration data (solid blue line). Measured data beyond 𝑡valid
are separately displayed in light red. Philip’s two-term equation is in excellent agreement with the measured infiltration data for 𝑡 ≤ 𝑡valid . This is true for all soil types. For
infiltration times larger than 𝑡valid , Philip’s two-term equation is known to be deficient in describing the infiltration
process, and the simulated curves deviate more and morefrom the measured data. The cumulative infiltration at 𝑡 =𝑡valid ranges between 2.5 to 3.2 cm but is limited to about
1.8 cm for sandy clay loam. Thus, a cumulative infiltrationof about 2.8–3.0 cm may serve as a reasonable proxy forthe time validity of Philip’s two-term equation. In fact, at
𝑡 = 𝑡valid the soil has reached an approximately constant rate of
infiltration.

JAISWAL ET AL.

Vadose Zone Journal

19 of 27F I G U R E 8 Comparison of observed (red dots) and simulated cumulative infiltration, I, using Philip’s two-term infiltration Equation 5 withleast squares values of the sorptivity, S, and saturated hydraulic conductivity, 𝐾s , derived from Parlange’s equation (Step 1, see Table 5) andmultiplicative coefficient, c, and time validity, 𝑡valid , from Step 2 of our methodology (e.g., see Table 6). Each graph corresponds to a different soil
type. Measured cumulative infiltration data beyond the time validity of Philip’s two-term equation are displayed with a tint of red

F I G U R E 9 Textural triangle of the 12 USDA soil types withcolor coding of the least squares values of the time validity, 𝑡valid , ofPhilip’s two-term infiltration equation. The colorbar assigns values tothe palette of the time validity

Figure 9 links the time validity of Philip’s two-term infiltration equation to soil texture using color coding for eachsoil type. The rather large values of 𝑡valid for silty clay andsilty clay loam suppress the subtle color variations in the time

validity of the other soil types. Along with Table 6 it is not
too difficult to discern a trend in the inferred time validities
of Philip’s two-term infiltration equation. Indeed, the finer the
texture of the soil the larger the value of 𝑡valid will be. We may
turn this trend into a regression equation to anticipate the time
validity of Philip’s two-term expression for other soils, yet thisdemands textural information. If, however, the interest is onlyin the hydraulic properties of the soil, then Parlange’s infiltration equation will suffice to back out the soil’s sorptivity andsaturated hydraulic conductivity.
Next, we present the texture triangle with color of eachsoil type proportional to the maximum likelihood value ofthe multiplicative coefficient c (see Figure 10). The colored
triangle discloses a pattern in the optimum value of c inPhilip’s two-term infiltration equation. In general, the valueof c decreases along two of the principal axes of the triangle,
namely, the percentages of sand and clay. Not all soil types
honor this relationship—for example, sandy clay loam, oneof our two dissonant soils.
Now that we completed Step 2, we can investigate therelationship between the least squares value of coefficient cin Philip’s two-term infiltration equation, and the maximum
likelihood value of β in Parlange’s infiltration equation. Figure 11 visualizes this relationship.
The dashed black line plots the identity, 𝑐 = (2∕3) −(1∕3)β, of Equation 11, valid for very short to short infiltration

20 of 27

JAISWAL ET AL.

Vadose Zone Journal

F I G U R E 1 0 Textural triangle of the 12 USDA soil types withprojection of the least squares values of the multiplicative
dimensionless coefficient, c, of Philip’s two-term infiltration equation.
The colorbar assigns values to coefficient, ctimes. The data points of sandy clay loam and silty clay loam
are marked in red as their respective β̂ values have been trun̂ 𝑐̂ ) data pairs of these twocated at 2. Thus, we discard the (β,̂ 𝑐̂ ) scatterdissonant soils in our subsequent analysis. The (β,plot illustrates several interesting findings. First and foremost,the data pairs appear to lay on a line, confirming the presenceof a linear relationship between the dimensionless coefficients
β and c. Secondly, the linearity of the β(𝑐) relationship ofEquation 11 appears to extend beyond early infiltration times.
̂ 𝑐̂ ) data pairs do not match perfectly the dashed
Third, the (β,black line, valid for very short to short infiltration times. Lin̂ 𝑐̂ ) data without the two dissonant
ear least squares of the (β,soils produces the relationship 𝑐 = 0.723 − 0.317β between
coefficients β and c of Parlange and Philip’s infiltration equations, respectively. The slope and intercept of this regression function (dotted blue line) are only marginally different from their counterparts of Equation 11; nevertheless, theyextend up to the time validity of Philip’s two-term infiltration equation. Thus, the identity 𝑐 = 0.723 − 0.317β produces
an almost perfect match with the maximum likelihood values of Parlange’s β coefficient derived from the DREAM(ZS)
algorithm in Step 1 of our methodology. This new identity isapproximately similar to11𝑐= √ ψ− β
π
5

(19)

√
where ψ = (1∕2) + (1∕2) 5 ≈ 1.6180 denotes the so-called
golden ratio.

5

F I G U R E 1 1 Comparison of the DREAM(ZS) –derived maximum
likelihood value of parameter β of Parlange’s semi-implicit infiltration
equation (Step 1) and the least squares value of coefficient c of Philip’s
two-term infiltration equation derived in Step 2 of our methodology(values tabulated in Table 6). The solid black line portrays thetheoretical relationship of Equation 11 between β and c, valid at short
infiltration times. The two red colored data points correspond to our
dissonant soils, sandy clay loam, and silty clay loam, whose histogramsof β are truncated by the prior distribution (see Figures 4E3, 5E2,
and 5E3). The dotted blue line equals the least squares fit through thê 𝑐̂ ) pairs of the dissonant soils
data, but without the two (β,IMPLICATIONSThe scatter plot in Figure 11 confirms the presence of a linear𝑐(β) relationship within the time validity of Philip’s two-term
infiltration equation. We can take advantage of this identity tosimplify our methodology. After Step 1 is completed, we canturn the sampled posterior values of β into a marginal distribution of the dimensionless coefficient c, hereafter referred to as𝑐β , using the relationship of Equation 19. Next, we can evaluate Philip’s two-term infiltration equation for each augmented
vector, η+ = [𝑆 𝐾s β 𝑐β ]⊤ , of posterior parameter values andcompute the 2.5 and 97.5% simulation uncertainty ranges at
each measurement time.
Figure 12 presents the marginal distribution of 𝑐β and displays the associated 95% prediction ranges of Philip’s twoterm infiltration equation for (a) clay, (b) clay loam, (c) sand,(d) silt loam, and (e) silty clay loam. Measured infiltration
data beyond 𝑡valid use a different red tint. The most important results are as follows. The marginal distribution of 𝑐β centers on 𝑐̂ derived from Step 2 of our methodology and is well

JAISWAL ET AL.

Vadose Zone Journal

21 of 27F I G U R E 1 2 Top panel: Marginal distribution of the multiplicative coefficient, c, of Philip’s two-term infiltration equation derived from
Equation 19, for a representative group of five soil types, including (a) clay, (b) clay loam, (c) sand, (d) silt loam, and (e) silty clay loam. The leastsquares values of c from Step 2 of our methodology are separately indicated in each graph with the red cross symbol. Bottom panel: Comparison ofobserved (red dots) and 95% prediction ranges of the cumulative infiltration (I, solid blue lines) simulated with Philip’s two-term infiltration
Equation 5 using the posterior samples of the sorptivity, S, saturated hydraulic conductivity, 𝐾s , and coefficient, 𝑐β , derived from Parlange’s equationin Step 1 of our methodology. Thus, the two infiltration curves correspond to the 2.5 and 97.5% prediction quantiles derived from the posterior
samples. Measured data beyond the time validity of Philip’s two-term equation inferred in Step 2 of our methodology are displayed with a tint of red

described by a normal distribution for all but the two dissonantsoils. The histogram of 𝑐β for silty clay loam (Figure 12E1)
suffers from truncation at β = 2, resulting in a skewed distribution with right tail. Its least squares value appears outsidê 𝑐̂ ) data pairs of thethe marginal distribution of 𝑐β as the (β,two dissonant soils were omitted in the computation of the𝑐(β) relationship of Equation 19. This discrepancy between 𝑐̂and the histogram of 𝑐β values causes the 2.5 and 97.5% prediction quantiles of Philip’s two-term equation for the silty
clay loam soil to deviate from the measured infiltration data
before 𝑡valid is reached. Upon closer inspection of the residuals, one would anticipate a time validity for silty clay loamon the order of 7–8 h. For the other four soils, Philip’s 2.5 and97.5% prediction quantiles (Figure 12A2–D2) match perfectly
the measured infiltration data for 0 ≤ 𝑡 ≤ 𝑡valid . As expected,
at larger infiltration times, the two prediction quantiles deviate
̃ data pairs.
systematically from the measured (𝑡̃, 𝐼)
In Step 1, we made the assumption of normally distributed
infiltration data measurements errors. The variance of themeasurement errors was integrated out to yield a simple formulation of the likelihood function. If we specify the cumulative infiltration measurement error up front, then we can use
a χ2 test to turn each residual of Philip’s infiltration equation and the measured cumulative infiltration data into an
exceedance probability. This then can be used to pinpoint
the time validity of Philip’s two-term infiltration equation.
As we can repeat this exercise for each posterior simulation,
this gives rise to a probability distribution of 𝑡valid for eachsoil type. This refined method carries forward any assumptions we may make with respect to the measurement errors of

the cumulative infiltration data, from definition of the likelihood function, to the posterior samples of the DREAM(ZS)
algorithm, the marginal distribution of 𝑐β , the chi-squared
test statistic, and resulting time validities. Indeed, this refined
method address more explicitly the uncertainty inherent in ourestimates of 𝑡valid .
We would be remiss not to comment on the use of realworld infiltration data. In our preliminary investigations presented herein, we purposely used synthetic infiltration data sothat we can benchmark our sorptivity and saturated hydraulicconductivity estimates. Nothing prevents us from using eitherof the two methods to measured infiltration data. As Parlange’s infiltration equation does not suffer a limited time
validity, we can always use Step 1 of our methodology to infer
the posterior parameter distribution of S, 𝐾s , and β. Then, in
Step 2, we can iterate over the measurements to determine
the time validity of Philip’s infiltration model and associatedvalue of the multiplicative coefficient, c. Alternatively, one
can use the refined approach describe above as replacement ofStep 2, and turn the sampled β values into a marginal distribution of c values, followed by computation of the exceedance
probabilities of the measured infiltration data. This ultimately
provides a distribution of time validities for each measured
infiltration dataset. Indeed, experiments carried out on measured data from the SWIG database confirm our main conclusions. However, the inferred time validities are less subtle, as
the measured infiltration data generally exhibit a much poorer
temporal resolution.
We like to comment last on the practical usefulness of our
methodology. If the interest is only in the hydraulic properties22 of 27of the soil, then Step 1 of our methodology will suffice. The
MATLAB code of Appendix A coupled with the DREAM(ZS)
algorithm realizes the full potential of Parlange’s infiltration
equation for inverse estimation of soil hydraulic properties
and will help rectify the abundant use by practitioners of moreconvenient analytical infiltration functions that suffer physical underpinning and/or a limited time validity. The posterior
distribution of the Parlange parameters, S, 𝐾s , and β is a witness of data quality and quantity. With few or poor infiltration
measurements at early times, for example, the soil sorptivity
will suffer and exhibit a relatively large posterior uncertainty.
Short experiments, on the contrary, will complicate the inference of β as its impact on the simulated I(t) relationship is negligible at early times of infiltration. This is easy to verify fromthe sensitivity matrix of Parlange’s infiltration equation andpromotes a flat posterior marginal distribution of β. The posterior uncertainty of S, 𝐾s , and β can be propagated forward
to quantify prediction intervals. What is more, if we adopt the𝑐(β) relationship of Equation 19, then Step 1 yields a marginal
posterior distribution of coefficient c as well. The trivariate
distribution of S, 𝐾s , and c can enter Philip’s two-term infiltration equation. The second step of our methodology serves aperhaps more academic purpose. Nevertheless, knowledge ofthe coefficient c and time validity, 𝑡valid , of Philip’s two-term
expression is of utmost importance in the practical applicationof this infiltration equation.

6

JAISWAL ET AL.

Vadose Zone Journal

SUMMARY AND CONCLUSIONS√The two-term infiltration equation, 𝐼(𝑡) = 𝑆 𝑡 + 𝑐𝐾s 𝑡, ofPhilip (1957) is widely used to infer the soil sorptivity, S, andsaturated hydraulic conductivity, 𝐾s , from a n-record of measured cumulative infiltration data, {𝑡̃𝑖 , 𝐼̃𝑖 }𝑛𝑖=1 . This infiltration
equation is easy to use, supported by detailed mathematical–
physical analysis, and the optimum values of S and 𝐾s can beobtained from linear least squares via a closed-form expression. Despite this progress made, Philip’s two-term infiltration function has a limited time validity, 𝑡valid , and consẽ 𝑡̃), beyond 𝑡 = 𝑡valid ,quently, cumulative infiltration data, 𝐼(
will corrupt the least squares estimates of S and 𝐾s .
The contributions of this paper were twofold. First, we
introduced a method that can successfully infer the sorptivity,
S, saturated hydraulic conductivity, 𝐾s , dimensionless coefficient, c, and time validity, 𝑡valid , of Philip’s two-term infiltration equation.
As prerequisite of this methodology, the second contribution of this work is a robust, exact and efficient numericalsolution of Parlange’s implicit infiltration equation. The proposed methodology was tested on 12 USDA soil types using
HYDRUS-1D simulated infiltration data. Results demonstrated an excellent match between the inferred values of S
and 𝐾s and their “true” values with estimates of c and 𝑡validthat correlate well with soil texture and corroborate linearity of the 𝑐(β) relationship for infiltration times up to 𝑡valid .
Additional experiments carried out on measured infiltration
data from the SWIG database of Rahmati, Weihermüller,
Vanderborght, et al. (2018) and Rahmati, Weihermüller, andVereecken (2018) confirmed our main findings.
In the penultimate section of this paper, we described amodification to Step 2 of our method which takes advantageof the 𝑐(β) relationship to determine the posterior parameter
distribution of Philip’s two-term infiltration equation. Its time
validity can then be inferred using a χ2 test of the cumulative infiltration residuals of each posterior simulation. Thisapproach propagates forward assumptions made about themeasurement error of the infiltration data to a marginal distribution of time validities for each soil.
AC K N OW L E D G M E N T S
We appreciate the constructive comments of the three anonymous reviewers.
AU T H O R C O N T R I B U T I O N S
Parakh Jaiswal: Formal analysis; Investigation; Visualization. Yifu Gao: Data curation; Investigation; Software. Mehdi
Rahmati: Data curation; Investigation; Software. Jan Vanderborght: Conceptualization; Methodology. Jirka Šimůnek:Software; Validation. Harry Vereecken: Conceptualization.
Jasper A. Vrugt: Conceptualization; Formal analysis; Investigation; Methodology; Software; Supervision; Validation;
Visualization; Writing-original draft; Writing-review & editing.
CONFLICT OF INTEREST
The authors declare no conflict of interest.
D A T A A N D S O F T WA R E AVA I L A B I L I T Y
The data and software are available upon request fromthe corresponding author (jasper@uci.edu) and can be
downloaded at the following URL: https://www.dropbox.
com/sh/xpgpe8mu9pquyva/AAASIRsnJU7Q2FO6RAF_
MD4Ta?dl=0.
ORCID
Mehdi Rahmati https://orcid.org/0000-0001-5547-6442
Harry Vereecken https://orcid.org/0000-0002-8051-8517
Jasper A. Vrugt https://orcid.org/0000-0003-2599-1165
REFERENCES
Assouline, S. (2013). Infiltration into soils: Conceptual approaches andsolutions. Water Resources Research, 49, 1755–1772. https://doi.org/
10.1002/wrcr.20155
Barry, D. A., Parlange, J.-Y., Haverkamp, R., & Ross, P. J. (1995).
Infiltration under ponded conditions: 4. An explicit predictive infil-

JAISWAL ET AL.

tration formula. Soil Science, 160, 8–17. https://doi.org/10.1097/
00010694-199507000-00002
Basha, H. A. (1999). Multidimensional linearized nonsteady infiltration with prescribed boundary conditions at the soil surface.
Water Resources Research, 35, 75–93. https://doi.org/10.1029/
1998WR900015
Basha, H. A. (2000). Multidimensional quasi-linear steady infiltration
toward a shallow water table. Water Resources Research, 36(7), 1697–
1705. https://doi.org/10.1029/2000WR900064
Basha, H. A. (2002). Burger’s equation: A general nonlinear solutionof infiltration and redistribution. Water Resources Research, 38(11).
https://doi.org/10.1029/2001WR000954
Braester, C. (1973). Moisture variation at the soil surface and theadvance of the wetting front during infiltration at constant flux.
Water Resources Research, 9(11), 687–694. https://doi.org/10.1029/
WR009i003p00687
Brooks, S. P., & Gelman, A. (1998). General methods for monitoring
convergence of iterative simulations. Journal of Computational andGraphical Statistics, 7, 434–455.
Brutsaert, W. (1977). Vertical infiltration in dry soil. Water Resources
Research, 13, 363–368. https://doi.org/10.1029/WR013i002p00363
Carsel, R. F., & Parrish, R. S. (1988). Developing joint probability distributions of soil water retention characteristics. Water
Resources Research, 24(5), 755–769. https://doi.org/10.1029/
WR024i005p00755
Chen, C., & Gallipoli, D. (2004). Steady infiltration from buried point
source into heterogeneous cross-anisotropic unsaturated soil. International Journal for Numerical and Analytical Methods in Geomechanics, 28, 1033–1055. https://doi.org/10.1002/nag.370
Chen, X., Luo, Z., & Zhou, S. (2014). Influences of soil hydraulic andmechanical parameters on land subsidence and ground fissures causedby groundwater exploitation. Journal of Hydrodynamics, 26(1), 155–
164. https://doi.org/10.1016/S1001-6058(14)60018-4
Collis-George, N. (1977). Infiltration equations for simple soil systems. Water Resources Research, 13, 395–403. https://doi.org/10.
1029/WR013i002p00395
Dohnal, M., Dusek, J., & Vogel, T. (2010). Improving hydraulic conductivity estimates from minidisk infiltrometer measurements for soils
with wide pore-size distributions. Soil Science Society of America
Journal, 74, 804–811. https://doi.org/10.2136/sssaj2009.0099
Fernández-Gálveza, J., Pollacco, J. A. P., Lassabatere, L., AnguloJaramillo, R., & Carrick, S. (2019). A general Beerkan estimationof soil transfer parameters method predicting hydraulic parametersof any unimodal water retention and hydraulic conductivity curves:
Application to the Kosugi soil hydraulic model without using particle
size distribution data. Advances in Water Resources, 129, 118–130.
https://doi.org/10.1016/j.advwatres.2019.05.005
Fisher, R. A. (1934). Statistical methods for research workers. Oliver &Boyd.
Fuentes, C., Haverkamp, R., & Parlange, J.-Y. (1992). Parameter constraints on closed-form soil water relationships. Journal of Hydrology,
134, 117–142. https://doi.org/10.1016/0022-1694(92)90032-Q
Green, W. H., & Ampt, G. A. (1911). Studies on soil physics. The
Journal of Agricultural Science, 4(1), 1–24. https://doi.org/10.1017/
S0021859600001441
Grigorjev, V. Y., & Iritz, L. (1991). Dynamic simulation model of vertical infiltration of water in soil. Hydrological Sciences Journal, 36(2),171–179. https://doi.org/10.1080/02626669109492497

Vadose Zone Journal

23 of 27

Haverkamp, R., Parlange, J.-Y., Starr, J. L., Schmitz, G., & Fuentes, C.
(1990). Infiltration under ponded conditions: 3. A predictive equation
based on physical parameters. Soil Science, 149, 292–300. https://doi.
org/10.1097/00010694-199005000-00006
Haverkamp, R., Ross, P. J., Smettem, K. R. J., & Parlange, J. - Y. (1994).
Three-dimensional analysis of infiltration from the disc infiltrometer:
2. Physically based infiltration equation. Water Resources Research,30(11), 2931–2935. https://doi.org/10.1029/94WR01788
Hellin, J., & Schrader, K. (2003). The case against direct incentives andthe search for alternative approaches to better land management inCentral America. Agriculture, Ecosystems & Environment, 99(1), 61–81.
Holtan, H. N. (1961). A concept for infiltration estimates in watershedengineering. USDA-ARS.
Horton, R. E. (1933). The role of infiltration in the hydrologic cycle. Eos,Transactions American Geophysical Union, 14(1), 446–460. https://doi.org/10.1029/TR014i001p00446
Horton, R. E. (1941). An approach toward a physical interpretation of infiltration-capacity. Soil Science Society of America Journal, 5(100), 399–417. https://doi.org/10.2136/sssaj1941.
036159950005000C0075x
Huggins, L. F., & Monke, E. J. (1966). The mathematical simulation of
the hydrology of small watersheds (Indiana Water Resources ResearchCenter Technical Reports 8-1-1966). Purdue University.
Hunt, A. G., Holtzman, R., & Ghanbarian, B. (2017). A percolationbased approach to scaling infiltration and evapotranspiration. Water,
9(2). https://doi.org/10.3390/w9020104
Kostiakov, A. N. (1932). On the dynamics of the coefficient of waterpercolation in soils and the necessity of studying it from the dynamic
point of view for the purposes of amelioration. In O. Fauser (Ed.),
Transactions of Sixth Committee International Society of Soil Science(pp. 7–21). International Society of Soil Science.
Laloy, E., Linde, N., Jacques, D., & Vrugt, J. A. (2015). Probabilistic inference of multi-Gaussian fields from indirect hydrological data using circulant embedding and dimensionality reduction.
Water Resources Research, 51, 4224–4243. https://doi.org/10.1002/
2014WR016395
Lassabatere, L., Angulo-Jaramillo, R., Soria-Ugalde, J. M., Šimůnek, J.,& Haverkamp, R. (2009). Numerical evaluation of a set of analyticalinfiltration equations. Water Resources Research, 45(12). https://doi.
org/10.1029/2009WR007941
Latorre, B., Moret-Fernández, D., Lassabatere, L., Rahmati, M., López,M. V., Angulo-Jaramillo, R., Sorando, R., Comín, F., & Jiménez, J. J.
(2018). Influence of the beta parameter of the Haverkamp model on
the transient soil water infiltration curve. Journal of Hydrology, 564,
222–229. https://doi.org/10.1016/j.jhydrol.2018.07.006
Latorre, B., Peña, C., Lassabatere, L., Angulo-Jaramillo, R., & MoretFernández, D. (2015). Estimate of soil hydraulic properties from disc
infiltrometer three-dimensional infiltration curve. numerical analysis
and field application. Journal of Hydrology, 527, 1–12. https://doi.
org/10.1016/j.jhydrol.2015.04.015
Linde, N., & Vrugt, J. A. (2013). Distributed soil moisture from crosshole ground-penetrating radar travel times using stochastic inversion.
Vadose Zone Journal, 12(1), 1–16. https://doi.org/10.2136/vzj2012.
0101
Liu, H., Lei, T. W., Zhao, J., Yuan, C. P., Fan, Y. T., & Qu, L. Q.
(2011). Effects of rainfall intensity and antecedent soil water content
on soil infiltrability under rainfall conditions using the run off-on-

24 of 27

Vadose Zone Journal

out method. Journal of Hydrology, 396(1), 24–32. https://doi.org/10.
1016/j.jhydrol.2010.10.028
Mein, R. G., & Larson, C. L. (1971). Modeling the infiltration component of the rainfall-runoff process. Water Resources Research Center,
University of Minnesota.
Mein, R. G., & Larson, C. L. (1973). Modeling infiltration during a
steady rain. Water Resources Research, 9(2), 384–394. https://doi.org/
10.1029/WR009i002p00384
Moret-Fernández, D., Latorre, B., & Angulo-Martínez, M. (2017). Comparison of different methods to estimate the soil sorptivity froman upward infiltration curve. Catena, 155, 86–92. https://doi.org/10.
1016/j.catena.2017.03.009
Moret-Fernández, D., Latorre, B., López, M. V., Pueyo, Y., Lassabatere,
L., Angulo-Jaramilo, R., Rahmati, M., Tormo, J., & Nicolau, J.
M. (2020). Three- and four-term approximate expansions of the
Haverkamp formulation to estimate soil hydraulic properties fromdisc infiltrometer measurements. Hydrological Processes, 34(26),5543–5556. https://doi.org/10.1002/hyp.13966
Mualem, Y. (1976). A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resources Research, 12(3),513–522. https://doi.org/10.1029/WR012i003p00513
Parlange, J.-Y., Haverkamp, R., & Touma, J. (1985). Infiltration under
ponded conditions: 1. optimal analytical solution and comparison
with experimental observations. Soil Science, 139, 305–311. https://doi.org/10.1097/00010694-198504000-00003
Parlange, J.-Y., Lisle, I., Braddock, R. D., & Smith, R. E. (1982). The
three-parameter infiltration equation. Soil Science, 133(6), 337–341.
https://doi.org/10.1097/00010694-198206000-00001
Philip, J. R. (1957). The theory of infiltration: 1. The infiltration equation and its solution. Soil Science, 83(5), 345–358. https://doi.org/10.
1097/00010694-195705000-00002
Philip, J. R. (1969a). Moisture equilibrium in the vertical in swellingsoils. I. Basic theory. Soil Research, 7(2), 99–120. https://doi.org/10.
1071/SR9690099
Philip, J. R. (1969b). Theory of infiltration. Advances in Hydroscience,
5, 215–296. https://doi.org/10.1016/B978-1-4831-9936-8.50010-6
Popper, K. (1992). The logic of scientific discovery. Routledge.
Qin, H., Vrugt, J. A., Xie, X., & Zhou, Y. (2018). Improved characterization of underground structure defects from two-stage Bayesian inversion using crosshole GPR data. Automation in Construction, 95, 233–
244. https://doi.org/10.1016/j.autcon.2018.08.014
Rahmati, M., Latorre, B., Lassabatere, L., Angulo-Jaramillo, R., &Moret-Fernández, D. (2019). The relevance of Philip theory to
Haverkamp quasi-exact implicit analytical formulation and its uses
to predict soil hydraulic properties. Journal of Hydrology, 570, 816–
826. https://doi.org/10.1016/j.jhydrol.2019.01.038
Rahmati, M., Vanderborght, J., Šimůnek, J., Vrugt, J. A., MoretFernández, D., Latorre, B., Lassabatere, L., & Vereecken, H. (2020).
Soil hydraulic properties estimation from one-dimensional infiltration
experiments using characteristic time concept. Vadose Zone Journal,
19, 1–22. https://doi.org/10.1002/vzj2.20068
Rahmati, M., Weihermüller, L., Vanderborght, J., Pachepsky, Y. A.,
Mao, L., Sadeghi, S. H., Moosavi, N., Kheirfam, H., Montzka, C.,
Van Looy, K., Toth, B., Hazbavi, Z., Al Yamani, W., Albalasmeh,
A. A., Alghzawi, M. Z., Angulo-Jaramillo, R., Antonino, A. C. D.,
Arampatzis, G., Armindo, R. A., . . . , & Vereecken, H. (2018). Development and analysis of the Soil Water Infiltration Global database.
Earth System Science Data, 10(3), 1237–1263. https://doi.org/10.
5194/essd-10-1237-2018

JAISWAL ET AL.

Rahmati, M., Weihermüller, L., & Vereecken, H. (2018). Soil Water
Infiltration Global (SWIG) database [Supplement to: Rahmati, M
et al. (2018): Development and analysis of Soil Water Infiltration
Global database]. Earth System Science Data. https://doi.org/10.
1594/PANGAEA.885492
Salvucci, G., & Entekhabi, D. (1994). Explicit expressions for GreenAmpt (delta function diffusivity) infiltration rate and cumulative storage. Water Resources Research, 30, 2661–2663. https://doi.org/10.
1029/94WR01494
Schaap, M. G., & van Genuchten, M. Th. (2006). A modified Mualemvan Genuchten formulation for improved description of the hydraulicconductivity near saturation. Vadose Zone Journal, 5(1), 27–34. https://doi.org/10.2136/vzj2005.0005
Schoups, G., Hopmans, J. W., Young, C. A., Vrugt, J. A., Wallender,
W. W., Tanji, K. K., & Panday, S. (2005). Sustainability of irrigated agriculture in the San Joaquin Valley, California. Proceedingsof the National Academy of Sciences, 102(43), 15352–15356. https://doi.org/10.1073/pnas.0507723102
Šimůnek, J., van Genuchten, M. Th., & Sejna, M. (2008). Development
and applications of the HYDRUS and STANMOD software packages
and related codes. Vadose Zone Journal, 7(2), 587–600. https://doi.
org/10.2136/vzj2007.0077
Šimůnek, J., van Genuchten, M. Th., & Sejna, M. (2016). Recent developments and applications of the HYDRUS computer software packages. Vadose Zone Journal, 15(7). https://doi.org/10.2136/vzj2016.
04.0033
Singh, V., & Yu, F. X. (1990). Derivation of infiltration equation using
systems approach. Journal of Irrigation and Drainage Engineering,
116, 837–858. https://doi.org/10.1061/(ASCE)0733-9437(1990)116:
6(837)
Smith, R. E. (1972). The infiltration envelope: Results from a theoreticalinfiltrometer. Journal of Hydrology, 17(1), 1–22. https://doi.org/10.
1016/0022-1694(72)90063-7
Smith, R. E., & Parlange, J.-Y. (1978). A parameter-efficient hydrologic
infiltration model. Water Resources Research, 14, 533–538. https://doi.org/10.1029/WR014i003p00533
Stroosnijder, L. (1976). Cumulative infiltration and infiltration rate inhomogeneous soils. Agricultural Research, 847, 69–99.
Valiantzas, J. D. (2010). New linearized two-parameter infiltration equation for direct determination of conductivity and sorptivity. Journalof Hydrology, 384(1–2), 1–13. https://doi.org/10.1016/j.jhydrol.2009.
12.049van Genuchten, M. Th. (1980). A closed-form equation for predicting
the hydraulic conductivity of unsaturated soils. Soil Science Society ofAmerica Journal, 44(5), 892–898. https://doi.org/10.2136/sssaj1980.
03615995004400050002x
Vogel, T., & Cislerova, M. (1988). On the reliability of unsaturated
hydraulic conductivity calculated from the moisture retention curve.
Transport in Porous Media, 3(1), 1–15. https://doi.org/10.1007/
BF00222683
Vogel, T., van Genuchten, M. Th., & Cislerova, M. (2000). Effect of
the shape of the soil hydraulic functions near saturation on variablysaturated ow predictions. Advances in Water Resources, 24(2), 133–
144. https://doi.org/10.1016/S0309-1708(00)00037-3
Volpi, E., Schoups, G., & Fermani, G. (2017). Sworn testimony of themodel evidence: Gaussian mixture importance (GAME) sampling.
Water Resources Research, 53, 6133–6158. https://doi.org/10.1002/
2016WR020167

JAISWAL ET AL.

Vrugt, J. A. (2016). Markov chain Monte Carlo simulation using theDREAM software package: Theory, concepts, and MATLAB implementation. Environmental Modeling & Software, 75, 273–316.
Vrugt, J. A., & Massoud, E. C. (2019). Uncertainty quantification ofcomplex system models: Bayesian analysis. In Q. Duan, F. Pappenberger, A. Wood, H. L. Cloke, & J. C. Schaake (Eds.), Handbook ofhydrometeorological ensemble forecasting (pp. 563–636). Springer.
Vrugt, J. A., ter Braak, C. J. F., Clark, M. P., Hyman, J. M., & Robinson, B. A. (2008). Treatment of input uncertainty in hydrologic modeling: Doing hydrology backward with Markov chain Monte Carlosimulation. Water Resources Research, 44, W00B09. https://doi.org/
10.1029/2007WR006720
Vrugt, J. A., ter Braak, C. J. F., Diks, C., Robinson, B. A., Hyman, J.
M., & Higdon, D. (2009). Accelerating Markov chain Monte Carlosimulation by differential evolution with self-adaptive randomized
subspace sampling. International Journal of Nonlinear Sciences &Numerical Simulation, 10(3), 273–290.
Zhang, J., Vrugt, J. A., Shi, X., Lin, G., Wu, L., & Zeng, L. (2020).
Improving simulation efficiency of MCMC for inverse modelingof hydrologic systems with a Kalman-inspired proposal distribution. Water Resources Research, 56, 1–24. https://doi.org/10.1029/
2019WR025474
Zhu, J., & Mohanty, B. P. (2002). Analytical solutions for steady-statevertical infiltration. Water Resources Research, 38(8). https://doi.org/
10.1029/2001WR000398

S U P P O RT I N G I N F O R M AT I O N
Additional supporting information may be found in the online
version of the article at the publisher’s website.

How to cite this article: Jaiswal, P., Gao, Y.,
Rahmati, M., Vanderborght, J., Šimůnek, J.,
Vereecken, H., & Vrugt, J. A. (2021). Parasite
inversion for determining the coefficients andtime-validity of Philip’s two-term infiltration
equation. Vadose Zone J. e20166.
https://doi.org/10.1002/vzj2.20166

APPENDIX A
This appendix presents the implementation of our numerical solution of Parlange’s infiltration equation in MATLAB.
Based on the 4 × 1 vector, η, with values of the sorptivity,
S (cm h−1/2 ), saturated hydraulic conductivity, 𝐾s (cm h−1 ),unitless coefficient, β, and initial hydraulic conductivity, 𝐾i ,and the 𝑛 × 1 vector, t, with discretized values of time indimension of hour, the MATLAB function Parlange returns
the cumulative infiltration, I, and infiltration rate, i (in unitsof cm and cm h−1 , respectively), and an exit flag, flag, which

Vadose Zone Journal

25 of 27

reports the corresponding exit condition of the secant method(see Figure A1).
The user may specify two additional input arguments, tolfun and maxiter, which equal the desired tolerance on thefunction value at the final root and the maximum number of
iterations of the secant method, respectively. These last two
input arguments are optional. If omitted in the Parlange function call, tolfun and maxiter will be assigned default values of
10−12 and 20, respectively. Built-in functions are highlighted
with a underlining.
The Parlange function calls the function secant, whichsolves for the zero-point of Parlange’s residual function (S1.7)
at each time, t (see Figure A2).
The secant function has three required input arguments—
namely, fun, an anonymous function handle of Parlange’s
residual function, and Ii1 and Ii2, which define the initial
search bracket, [𝐼(0) , 𝐼(1) ], of the root of 𝑟(𝐼, 𝑡). The secant
method may not always be used in each iteration of the
Parlange function. For the secant method to work in practice, the residual function, 𝑟(𝐼, 𝑡), at the upper bound, 𝐼(1) ,of the initial search bracket, [𝐼(0) , 𝐼(1) ], cannot be infinity, nor be smaller than 10−10 , hence, 10−10 < 𝑟[𝐼(1) , 𝑡] <∞. If 𝑟[𝐼(1) , 𝑡] ≤ 10−10 , then our upper bound, 𝐼(1) , is very
close to the actual root. Hence, we can use, 𝐼(1) , as oursolution of 𝐼(𝑡). A value of 𝑟[𝐼(1) , 𝑡] equal to infinity isa numerical artifact. The term exp[2βξ(𝐼 − 𝐾i 𝑡)] in Parlange’s residual function (S1.7) may return infinity for certain (often implausible) combinations of S, 𝐾s , and β, as
the exponential function grows larger than its limit on a 64bit computer. This is called numerical overflow and is typically restricted to saturated conditions. At this time, the soil
has reached a constant infiltration rate and, thus, we canmake use of the identity, 𝐼(𝑘) = 𝐼(𝑘−1) + [𝑡(𝑘) − 𝑡(𝑘−1) ]𝐾s , tosolve for the remaining part of the cumulative infiltration
curve.
The Parlange function has been evaluated numerous times
for 𝐾i = 0 and 0 ≤ 𝑡 ≤ 240 using Monte Carlo simulation
with 𝑆 ∈ (0, 10] (cm h−1/2 ), 𝐾s ∈ (0, 50] (cm h−1 ), and β ∈(0, 2). The resulting 𝐼(𝑡) curves were compared against a
direction solution of 𝑡(𝐼) derived with 𝐾𝑖 = 0 using a swapof the independent and dependent variables. This comparisonhas shown that the simulated cumulative infiltration curvesof the Parlange function are almost exact with a negligible
numerical error. In 1, 000 separate trials, the code would terminate prematurely a handful of times in response to numerical underflow of the exponential function. This can be avoided
by setting a lower limit of, say, 10−4 on the values of β. Altogether, our experiments have shown that the Parlange function
is robust, efficient, and exact.

26 of 27

Vadose Zone Journal

FIGURE A1

MATLAB code of the Parlange function

JAISWAL ET AL.

JAISWAL ET AL.

FIGURE A2

Vadose Zone Journal

MATLAB code of the secant method

27 of 27