Received: 21 March 2020

Accepted: 30 August 2020

DOI: 10.1002/vzj2.20068
Vadose Zone Journal

ORIGINAL RESEARCH ARTICLE

Soil hydraulic properties estimation from one-dimensional
infiltration experiments using characteristic time concept
Mehdi Rahmati1,2

Jan Vanderborght2

David Moret-Fernández6
Harry Vereecken2

Borja Latorre6

Jirka Šimůnek3

Jasper A. Vrugt4,5Laurent Lassabatere7

1 Dep. of Soil Science and Engineering, Faculty of Agriculture, Univ. of Maragheh, Maragheh, Iran
2 Institute of Bio- and Geosciences: Agrosphere (IBG-3), Forschungszentrum Jülich, Jülich, Germany
3 Dep. of Environmental Sciences, Univ. of California, Riverside, CA 92507, USA
4 Dep. of Civil and Environmental Engineering, Univ. of California, Irvine, CA 92697, USA
5 Dep. of Earth System Science, Univ. of California, Irvine, CA 92697, USA
6 Dep. de Suelo y Agua, Estación Experimental de Aula Dei, Consejo Superior de Investigaciones Científicas (CSIC), PO Box 13034, Zaragoza 50080,

Spain
7 Univ. Lyon, Univ. Claude Bernard Lyon 1, CNRS, ENTPE, UMR5023 LEHNA, Vaulx-en-Velin F-69518, France

Correspondence
Mehdi Rahmati, Dep. of Soil Science andEngineering, Faculty of Agriculture, Univ.
of Maragheh, Maragheh, Iran; Institute ofBio- and Geosciences: Agrosphere (IBG3), Forschungszentrum Jülich, Jülich,
Germany.
Email: mehdirmti@gmail.com.
Harry Vereecken, Institute of Bio- andGeosciences: Agrosphere (IBG-3),
Forschungszentrum Jülich, Jülich, Germany.
Email: h.vereecken@fz-juelich.de

Abstract
Many different equations ranging from simple empirical to semi-analytical solutions of the Richards equation have been proposed for quantitative description ofwater infiltration into variably saturated soils. The sorptivity, S, and the saturated
hydraulic conductivity, Ks , in these equations are typically unknown and have tobe estimated from measured data. In this paper, we use so-called characteristic
time (tchar ) to design a new method, referred to as the characteristic time method(CTM) that estimates S, and Ks , from one-dimensional (1D) cumulative infiltration data. We demonstrate the usefulness and power of the CTM by comparing
it with a suite of existing methods using synthetic cumulative infiltration data
simulated by HYDRUS-1D for 12 synthetic soils reflecting different USDA textural classes, as well as experimental data selected from the Soil Water Infiltration
Global (SWIG) database. Results demonstrate that the inferred values of S and
Ks are in excellent agreement with their theoretical values used in the synthetically simulated infiltration experiments with Nash–Sutcliffe criterion close tounity and RMSE values of 0.04 cm h−1/2 and 0.05 cm h−1 , respectively. The CTM
also showed very high accuracy when applied on synthetic data with added measurement noise, as well as robustness when applied to experimental data. Unlike
previously published methods, the CTM does not require knowledge of the time

Abbreviations: CF2, two-term curve-fitting model; CF3, three-term curve-fitting model; CTM, characteristic time method; 1D, one-dimensional;
QEI, quasi-exact implicit; SH, Sharma; SWIG, Soil Water Infiltration Global; VGM, van Genuchten–Mualem.
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided theoriginal work is properly cited.
© 2020 The Authors. Vadose Zone Journal published by Wiley Periodicals LLC on behalf of Soil Science Society of America
Vadose Zone J. 2020;19:e20068.
https://doi.org/10.1002/vzj2.20068

wileyonlinelibrary.com/journal/vzj2

1 of 22

2 of 22

RAHMATI et al.

Vadose Zone Journal

validity of the applied semi-analytical solution for infiltration and, therefore, isapplicable to infiltrations with durations from 5 min to several days. A script written in Python of the CTM method is provided in the supplemental material.

1

INTRODUCTION

Accurate and reliable estimation of hydraulic propertiesof variably saturated soils from infiltration experiments
is still a major challenge in many research fields suchas hydrology, irrigation, civil engineering, and environmental science (Moret-Fernández, Latorre, & AnguloMartínez, 2017; Valiantzas, 2010). In the past decades, alarge number of theoretical, semiempirical, and empirical models have appeared in the literature that quantitatively describe water infiltration into the soil (Green& Ampt, 1911; Haverkamp, Ross, Smettem, & Parlange,
1994; Parlange, Lisle, Braddock, & Smith, 1982; Philip, 1957;
Swartzendruber, 1987). These models were used for thecharacterization of soil hydraulic properties, including soil
sorptivity (S) and saturated hydraulic conductivity (Ks ).
According to their theoretical definitions, S is a measurable physical quantity, which expresses the capacity of aporous medium to take up and release liquids by capillarity (Philip, 1957), whereas Ks measures the soil’s ability totransmit water under the influence of gravity.
In general, there are two major approaches to estimate S
and Ks from one-dimensional (1D) infiltration experiments
considered in the literature: (a) linearization approachesand (b) inverse estimation of S and Ks using curve-fitting
methods. Table 1 summarizes the most often used linearization methods to predict S and Ks from 1D infiltration
curves using Philip’s (1957) two-term equation (I = St1/2 +cKs t, where t represents time and c is a constant).
The linearity of the linearized model proposed by Smiles
and Knight (1976) for long-time intervals indicates anappropriate description of the infiltration process by thetwo-term equation and that values of S and Ks can be determined from the y intercept and the slope of the line, respectively. However, there is still substantial arbitrariness indeciding when the linear part stops and nonlinearity starts,
and in detection of a plausibly linear relationship, leading

to uncertainty in estimated values depending on the choicemade.
The linearization method proposed by Sharma, Gander, and Hunt (1980) is motivated by the fact that capillary forces, and thus S, dominate the infiltration process atearly stages, whereas at later stages, the impact of capillary
forces decreases and the gravitational force becomes dominant. Therefore, they proposed to determine S from infiltration measurements at early time steps where
one can√assume a linear relationship between I and 𝑡 in Philip’s
two-term equation by setting Ks to zero. They also indicated that the I(t) vs. t plot usually shows a linear relationship for long times when infiltration reaches steadystate conditions. However, if the duration of infiltration
measurements is not long enough to reach such a steadystate condition, this method will, in most cases, overestimate Ks . They also indicated that the slope of the I(t) vs. t
plot is equal to Ks at the end. In fact, for short times (transient state), the slope is approximately equal to (2 − β)/3
× Ks (where β = 0.6; Haverkamp et al., 1994). When we
get closer to the steady state, the value of the slope tends
towards Ks . On the other hand, the steady-state infiltration corresponds to an intercept + Ks × t that may not be
negligible. This is especially the case for fine-textured soils
where the intercept, which is related to the ratio between
S2 and Ks (Haverkamp et al., 1994), takes large values.
We do not discuss the differentiated linearization
method of Vandervaere, Vauclin, and Elrick (2000), and
it is excluded from further comparisons. Indeed, AnguloJaramillo et al. (2000) found that this method is only
applicable for short to medium times, and its validity
may be questionable when a steady state is reached too
quickly. Latorre, Peña, Lassabatere, Angulo-Jaramillo, andMoret-Fernández (2015) also used this method to determine S and Ks from infiltration measurements using a disk
infiltrometer and demonstrated that it is inaccurate toderive soil hydraulic parameters. On the other hand, theT A B L E 1 Most often used linearization methods to predict soil sorptivity (S) and saturated hydraulic conductivity (Ks ) fromone-dimensional infiltration curves using Philip’s (1957) two-term equation (I = St1/2 + cKs t)
Procedure

Linearized model

Proposed by

Cumulative linearization

𝐼𝑡−1∕2 = 𝑆 + 𝑐𝐾s 𝑡1∕2

Smiles and Knight (1976)

Two-part linearization
Differentiated linearization

𝐼 = 𝑆𝑡1∕2

for early times; 𝐼 = 𝐾s 𝑡 for later times
1√√
Δ𝐼(𝑡)∕Δ 𝑡 = 𝑆 + 2𝐴2 ( 𝑡i 𝑡i+1 ) 2

Sharma et al. (1980)
Vandervaere et al. (2000)

Note. I, infiltration; t, time; c and A2 , infiltration constants; i, subscript ranging from 1 to the number of data points in the infiltration curve.

RAHMATI et al.

Vandervaere et al. (2000) method was mainly intended
for disc infiltrometer measurements, where a contact sandlayer is placed between the soil surface and the base disc.
According to Moret-Fernández, Blanco, Martínez-Chueca,
and Bielsa (2013), without a contact sand layer, the Smiles
and Knight (1976) method performs much better thanthe Vandervaere et al. (2000) one. Therefore, since weare working with 1D infiltration curves generated without
using a contact sand layer, it was rational to exclude thismethod from further analysis.
The second major approach uses curve fitting (Bonell &Williams, 1986; Bristow & Savage, 1987; Marquardt, 1963;
Vandervaere et al., 2000) of nonlinear infiltration equations in time to estimate soil hydraulic properties. Typically, an objective function that involves squared differences between measured and predicted infiltration rates
or cumulative infiltration volumes is minimized by changing S and Ks . Philip’s (1957) two-term equation, due to its
simplicity, is usually used with the curve fitting method todetermine S and Ks . However, Vandervaere et al. (2000)
pointed out that although the determination of hydraulicproperties (S and Ks ) using Philip’s two-term equation
is a well-posed problem, since it has only two degrees
of freedom√ for scale and shape, the obvious correlation
between 𝑡 and t makes it a relatively ill-conditioned problem where an increase in one parameter can be compensated by a decrease in the other parameter. Additionally,
the curve fitting of Philip’s two-term equation only provides physically consistent parameters if the time validity
domain of the equation is considered.
Although the curve-fitting procedure usually results in avery good agreement (R2 > .9) between measured and predicted infiltration curves, the predicted parameters may be
far from their true values. Equifinality, the principle thatthe minimum of the objective function can be obtainedby a broad set of parameter values (Beven & Freer, 2001),
is the main reason why no one can guarantee the results
obtained using the curve-fitting method. Vandervaere et al.
(2000) criticized this method by simply stating that “thebest fit does not necessarily mean good fit.”
In order to overcome the arbitrariness of the linearization method and equifinality of the curve-fitting method,
we present a new procedure based on the use of the characteristic time to uniquely predict S and Ks from infiltration experiments while implicitly accounting for the time
domain validity. Therefore, we use in our proposed method(hereafter called the characteristic time method [CTM])
the gravity time, tgrav , defined as the time when gravity
and capillarity have exactly the same impact on infiltration, to predict S and Ks . In a similar manner, tgrav was also
used to give an idea about the transition time that separates
the steady from the transient regimes of water infiltration(Lassabatere et al., 2006; Philip, 1969b). We test the proce-

3 of 22

Vadose Zone JournalCore Ideas∙ A new method named CTM is proposed to estimate the S and Ks from 1D infiltration curve.
∙ The CTM takes the advantages of the characteristic time in prediction of S and Ks .
∙ The CTM was more accurate than the nonlinear
curve fitting method in simultaneous predictionof S and Ks .
∙ The CTM was applicable to any infiltration
duration of 15 min to several days.

dure against the curve-fitting method using two- (CF2) or
three-term (CF3) models and the method of Sharma et al.
(1980) (SH). The CTM and SH methods that we compare
are not limited in time because they both treat the dataover the entire course of the infiltration (including transient and steady states).

2

THEORETICAL BACKGROUND

2.1

New procedure development

In order to characterize the 1D infiltration process, Philip(1957) derived a semianalytical solution of the Richards(1931) equation for ponded infiltration in unsaturated soils
for the transient state. Philip’s time series expansion for
cumulative infiltration can be written as
3

5𝐼(𝑡) = 𝐴1 𝑡0.5 + 𝐴2 𝑡 + 𝐴3 𝑡 2 + 𝐴4 𝑡2 + 𝐴5 𝑡 2
+

𝑛
∑

𝑖

𝐴𝑖 𝑡 2(1)𝑖=6where I [L] denotes the cumulative infiltration, t [T] signifies time, A1 [L T−1/2 ] to A5 [L T−1/2 ] are coefficients, and𝑛∑𝑖𝐴𝑖 𝑡 2𝑖=6represent higher order terms that are much smaller thanthe lower orders and can therefore be neglected. Philip(1957) showed that A1 is equal to the sorptivity S [L T−1/2 ],A2 [L T−1 ] is proportional to the saturated hydraulic conductivity Ks [L T−1 ], and the coefficients with n > 2 mustsatisfy
𝐴𝑛
>
𝑆

(

𝐴2
𝑆

)𝑛−1
(2)

4 of 22

RAHMATI et al.

Vadose Zone Journal

It is important to note that Equation 1 remains valid over
the time interval [0, tmax ], where tmax increases with n.
Accounting for more orders allows the extension of thetime validity interval. However, in all cases, the expansions
must be restricted to the given validity time intervals andthus address the case of transient state. Philip (1957) introduced the gravity time, tgrav [T], for the identification of thevalidity time for the case of the two-term expansion.
To remedy problems with the time validity of infiltration models, Haverkamp et al. (1994) derived a quasiexact implicit (QEI) solution of the Richards equation
to describe 1D cumulative infiltration, which was firstly
proposed by Parlange et al. (1982) and then redefined by
Haverkamp, Parlange, Starr, Schmitz, and Fuentes (1990):2(𝐾s −𝐾i )𝑆222 (𝐾 −𝐾 )[𝐼(𝑡)−𝐾 𝑡]

𝑆

θ

β=2−2

i

(

)(

𝐾s −𝐾iθs −θi(8)

𝐼(𝑡) = 𝑐(1)𝑡 1∕2 + 𝑐(2)𝑡 + 𝑐(3)𝑡 3∕2+ 𝑐(4)𝑡2 + 𝑐(5)𝑡5∕2(9)where c(1) to c(3) are defined as before, and c(4) and c(5)are defined as𝑐(4) =

𝐷(θ)dθ
(4)

where K(θ) and D(θ) are hydraulic conductivity and diffiusivty, respectively, of soil water at a given soil water content
of θ.
Haverkamp et al. (1994) proposed an approximate time
expansion of Equation 3 for the description of the transientstate considering a two-term equation in line with Philip(1957) two-term equation:
𝐼(𝑡) = 𝑐(1)𝑡0.5 + 𝑐(2)𝑡

(5)

where
𝑐(1) = 𝑆2−β𝑐(2) =𝐾s(6)

3

As stated above, the Equation 5 is valid for short tointermediate times only. Recently, Rahmati, Latorre, Lassabatere, Angulo-Jaramillo, and Moret-Fernández (2019)provided the three-term expansion of the QEI solution ofHaverkamp et al. (1994), which is valid for longer times
compared with Equation 5 and is identical to the Philip(1957) three-term equation:
𝐼(𝑡) = 𝑐(1)𝑡0.5 + 𝑐(2)𝑡 + 𝑐(3)𝑡 3∕2

(7)

𝐾s3
𝑆2

(10)

𝑐(𝑖)𝑡𝑖∕2, 𝑖 = 1, 2, … , 5
𝐼(𝑡)

(11)

𝑊1 + 𝑊2 + 𝑊3 + 𝑊4 + 𝑊5 = 1

(12)

𝑊𝑖 =

θ
∫θ s 𝐷(θ)dθi(β − 2) (β + 1) (1 − 2β)

If we divide both the right and left sides of Equation 9 byI(t), we get the following equation depicting the contribution (W) of each term to the infiltration process:

)

θ−θi

2

135
)2 𝐾 41 ( 2β − β + 1 3s
𝑐(5) =270𝑆

β𝐾(θ)−𝐾i) 𝐾21( 2β −β+1 s
9𝑆

Similarly to Rahmati et al. (2019), one can apply the Taylor series expansion (in powers of 0.5) to the QEI formulation to introduce the fifth (or any other) term of theapproximate expansion. The following equation presents
the approximate expansion up to the fifth term:(3)

where Ki and Ks are soil hydraulic conductivities at initialand saturated water contents of θi and θs , respectively, andβ (dimensionless) is an integral infiltration constant that iseither fixed at 0.6 or that can be computed soil specificallyas follows (Haverkamp et al., 1994):∫θ s𝑐(3) =

1

s
i
i
𝑡=
−
𝑆2
( 1−β{}1−β)
12β(𝐾s −𝐾i )[𝐼(𝑡)−𝐾i 𝑡]
β−1
ln exp+2βwhere c(1) and c(2) are as previously defined, and c(3) iswhere

In Equation 12, we assume that the contribution of thehigher order terms (orders > 5) are negligible and that we
remain in the validity time interval related to the five-term
approximate expansion. By plotting the contributions ofdifferent components to a simulated infiltration curve for
known values of S and Ks (Figure 1), we can outline thefollowing facts:
I. At the start of the infiltration process (t = 0), thecontribution of the first term to the infiltration process is maximal and equal to one, whereas the integrated contribution of remaining terms is equal to zero
(Figures 1b and 1c).
II. At a time equal to tgrav , the contribution of the firstterm is equal to the integrated contribution of theremaining terms, both being equal to 0.5 (Figure 1c).
III. Similar to second term, the third and subsequentterms act as gravity components since they have zero
contributions on infiltration at initial time steps andnonzero contribution at later time steps (Figure 1b)
where the contributions are increased by advancing intime.

RAHMATI et al.

Vadose Zone Journal

5 of 22F I G U R E 1 The (a) simulated infiltration curve for known values of soil sorptivity (S) and saturated hydraulic conductivity (Ks ), (b) temporal variations of the contributions (W) of different terms of the five-term equation to the infiltration process, and (c) temporal variations of
the first term contribution vs. the contribution of remaining terms. The β is an integral infiltration constant, t is time, I is cumulative infiltration
data, W is the contribution of the sorptivity and gravity components to the infiltration process, and c(1) to c(5) are constants used in approximate
expansions of the Haverkamp et al. (1994) formulation

IV. The contribution of the fourth and higher terms
is negligible compared with the contribution of thefirst three components (Figure 1b), and therefore one
can exclude the fourth and fifth components fromanalysis.
V. The variations of the contributions of the sorptivity and gravity components vs. time are symmetric(Figure 1c). This means that by advancing in time, any
decrease in the contribution of sorptivity component
results in a corresponding increase in the contributions of the gravity components.
VI. When the contribution of the sorptivity (or gravity)
component vs. time is plotted, an exponential form of
below is obtained (Figures 1b and 1c):
𝑊 = 𝑒α𝑡(13)where α (h−1 ) is a soil-dependent shape factor being negative for the sorptivity component and positive for the gravity component.
As a consequence of above facts, once tgrav has been
identified, the S and Ks can be estimated from thecumulative infiltration data using the iterative procedure
below.

2.2
CTM-I: Iterative procedure toestimate S and Ks
CTM-I is the first part of the CTM method. It considers agiven characteristic time (tchar ) falling between zero andtgrav with its related weight, ω:
at 0 < √𝑡char ≤ 𝑡grav ⇒ 𝑊1 = 1 − ω⇒

𝑆 𝑡char
=1−ω
𝐼char

(14)

where 1 − ω is the contribution of the first component to
the infiltration process, and Ichar is the cumulative infiltration at time equal to tchar . In the above equation, if tchar
is equal to tgrav (tchar = tgrav ), ω will be equal to 0.5. For
tchar < tgrav , smaller values of ω will be applied, whereas
for tchar ≈ 0, ω = 0 will be applied. The tchar can be related
to tgrav in the following manner:𝑡char = κ𝑡grav ,where 0 < κ ≤ 1(15)

Rearranging Equation 14 leads to the final solution for S,
enabling us to predict S from the knowledge of tchar :𝐼
𝑆 = (1 − ω) √char
𝑡char

(16)

6 of 22

RAHMATI et al.

Vadose Zone Journal

Considering that W4 and W5 in Equation 12 are negligible in comparison with W2 + W3 (Figure 1c), ω correspondswith W2 + W3 , leading to
at 0 < 𝑡char ≤ 𝑡grav ⇒ 𝑊2 + 𝑊3 = ω
3∕2⇒𝑐(2)𝑡char +𝑐(3)𝑡char

(17)

=ω

𝐼char

where c(2) and c(3) are defined as below:
2−β𝑐(2) =
𝐾s
( 2 3
) 𝐾2
𝑐(3) = β − β + 1 s

(18)

9𝑆

The rearrangement of Equation 17 based on Ks leads to thefollowing equation:
𝑎𝐾s2 + 𝑏𝐾s + 𝑐 = 0

(19)

where
𝑎=

1
9(1−ω)( 2) 𝑡2
β − β + 1 char

𝑏=

2−β
3

𝐼char

𝑡char

(20)

𝑐 = −ω𝐼char
where β can be fixed at 0.6 (Haverkamp et al., 1994).
Finally, solving for Ks in Equation 19 gives
𝐾s =

−𝑏 +

√

𝑏2 − 4𝑎𝑐
2𝑎(21)

where a, b, and c are defined as above.
In order to use Equations 16 and 21, one only needs toestimate from infiltration data the tchar − Ichar set as wellas ω. The following iterative procedure is applied to determine the tchar − Ichar set and ω from field- or laboratorymeasured 1D infiltration data in order to compute S and
Ks . The procedure is based on the fact that when estimated tchar (tchar,est ) is misestimated, the plot of the firstweight W1 leads to inconsistent values in the vicinity of
t = 0. In the case of underestimation of tchar , W1 tends
towards a value <1, whereas in the case of overestimationof tchar , W1 tends towards a value >1 (see illustrative examples in Figure 2). This implies that a correct estimate oftchar can be found by equaling tchar to the different valuesof the observed times (from cumulative infiltration measurements) and selecting tchar,est , for which West (t = 0) converges to 1.
The proposed procedure (Figure 3) aims at determining precisely the right tchar with its cumulated infiltration,
Ichar , and weight 1 − ω.
Step 1: Set k = 1 and ω equal to 0.5, meaning thattchar = tgrav as a first choice.

F I G U R E 2 The temporal variations of the first component’s
contribution estimated using different characteristic time (tchar ) values. This plot is obtained by setting tchar to values lower than gravity
time (tgrav ), equal to tgrav , and larger than tgrav and then by calculating the contribution of the sorptivity component at each data point oftime series of analytically generated infiltration data and plotting it
vs. time. The t is time, I is cumulative infiltration data, W1 is the contribution of the sorptivity component to the infiltration process, andc(1) is the constant used in approximate expansions of the Haverkamp
et al. (1994) formulation

Step 2: Set tchar (k) and Ichar (k) equal to t(k) and I(k) of theinfiltration measurement and estimates preliminary valueof S(k) using the following equations:𝐼(𝑘)
𝑆 (𝑘) = (1 − ω) √char
𝑡char (𝑘)(22)

Step 3: Compute the contribution of the sorptivity (W1 )component for all data points in the time series as follows:√
𝑆 (𝑘) 𝑡𝑊1 (𝑡) =
𝐼(𝑡)

(23)

Step 4: Determine the maximum value of W1 . Theoretically, if tchar (k) and Ichar (k) are set correctly, the maximum
W1 should be equal or very close to 1. Otherwise, maximum
W1 will be lower or higher than 1. To increase the flexibility of the procedure, we define the following criteria as aselection index:
1 − 𝑑 ≤ max(𝑊1 ) ≤ 1 + 𝑑(24)

where d is a positive value between 0 and 1. The lower thevalue of d, the more precise the estimations. In our case, ad value of 0.001 is used to obtain results with the required
precision (i.e., 10−3 ). In the case of noisy experimental infil-

RAHMATI et al.

Vadose Zone Journal

7 of 22F I G U R E 3 A flowchart illustrating the characteristic time method coupled with the iterative procedure (CTM-I) to find the characteristic
time (tchar ) and its corresponding cumulative infiltration value (Ichar ) and to predict soil sorptivity (S) and saturated hydraulic conductivity (Ks ).
The recommended value for both flexibility constants of ε and d is 0.001. The β is an integral infiltration constant, t is time, I is cumulative
infiltration data, and W1 and ω, respectively, are the contributions of the sorptivity and gravity components to the infiltration process.

tration data, this criterion may be too low, and one mayneed to widen the range defined in Equation 24 by using
higher values of d (e.g., a d value of 0.05 or 0.1 is proposed).
If the criterion of Equation 24 is achieved, end the process by finalizing the tchar − Ichar set and calculate the final
solutions for S and Ks using Equations 16 and 21. Otherwise, proceed to Step 5.
Step 5: Set k = k + 1 and tchar (k)= t(k) and Ichar (k)= I(k) and estimate updated value of S(k) using
Equation 22.
Step 6: Repeat Steps 3–5 until the criterion defined in
Equation 24 is achieved, allowing to select the right valueof the index k ≤ N (N is the number of observed datapoints). If a value of k can be found for ω = 0.5, it meansthat the observed infiltration data encompasses tgrav anddescribes both the transient and the steady states. Otherwise, if Equation 24 is never fulfilled with West (t = 0) =1, it means that the tgrav value is larger than the maximum
time, t(N), and that the observations describe the transientstate only. In this case, ω needs to be reduced and Steps
2–6 need to be repeated until tchar is found within the measurement time window [tchar ≤ t(N)]. It is recommended todecrease ω by removing small increment (e.g., 0.001) from
0.5 to the appropriate value.

2.3
CTM-Ks : A procedure for accurate
estimation of Ks
It is demonstrated and discussed in Section 4 that CTM-I
always provides accurate predictions for S but fails to provide accurate prediction for Ks when the final ω < 0.5. Inthat case, the following CTM-Ks procedure must be applied
for correct predictions of Ks . To this end, the prediction of
S from CTM-I (SCTM-I ) is used to compute the true contribution of the first component as follows:√
𝑆CTM−I 𝑡
𝑊1 (𝑡) =
𝐼(𝑡)

(25)

Then, with known time (t) vector from experimental/simulated infiltration data and contribution (W1 ) vector from above, one can fit the Equation 13 [which is W1= exp(αt)] on t–W set to predict α (<0). Then, the predicted α value can be used to determine tgrav for which
W1 = 1/2:

𝑡grav =

ln (0.5)α(26)

8 of 22

RAHMATI et al.

Vadose Zone Journal

F I G U R E 4 A flowchart illustrating the usage of the characteristic time method (CTM) for accurate predictions of both soil sorptivity (S)and saturated hydraulic conductivity (Ks ) parameters. ω is the contribution of the gravity component to the infiltration process

Then, one could calculate Igrav by rearranging the Equation 13 and putting W = 0.5 as below:
𝐼grav =√
𝑆CTM−I 𝑡grav
0.5

(27)

Finally, when tgrav and Igrav values were predicted as above,
Equation 21 can be used to determine Ks by setting ω = 0.5,
tchar = tgrav , and Ichar = Igrav .
Figure 4 briefly shows the flowchart of how CTM-I and
CTM-Ks should be used to ensure accurate prediction ofboth S and Ks .
A script written in Python applying the CTM is provided
as supplemental material. The Appendix also provides the
script in text format. A script written in MATLAB or scilabalso could be provided upon to request.

3
3.1

MATERIALS AND METHODS
Numerically simulated test data

3.1.1
Numerical simulations with
HYDRUS-1D
To test the proposed procedure, synthetic infiltration
curves were simulated by using HYDRUS-1D (Šimůnek,van Genuchten, & Šejna, 2008, 2016) for 12 USDA soil
textural classes for a period of 240 h to reach the steadystate for most soils. The average soil hydraulic parametersof the van Genuchten–Mualem (VGM) (van Genuchten,1980) model for each textural class (Table 2) were obtained
from the HYDRUS-1D soil catalog (Carsel & Parrish, 1988).
A pressure head (equal to zero) boundary condition was
imposed at the soil surface, and a free drainage boundary
condition was specified at the bottom of the soil profile.

The 200-cm-deep soil profile was assumed homogeneous
and was discretized into 401 nodes. Since the profile discretization may affect the numerical results, the accuracyof the time and space discretization was tested against theanalytical solution for infiltration without gravity, which
is based on the Boltzmann transform. To do this, we simulated water content profiles for 20 different times between
0.01 and 240 h using different spatial resolutions for the
profile discretization. Then, water contents were plottedagainst the Boltzmann variable λ as defined in the following equation (Philip, 1969a):λ(θ) = 𝑍(θ, 𝑡)𝑡 −1∕2(28)

where Z(θ,t) is a characteristic function, which quantifies
the depth at which the volumetric water content equatesto θ at time t during an infiltration event under gravityfree conditions (or during horizontal water absorption).
In the case of accurate numerical simulations, all θ vs. λcurves at different times should coalesce on one curve. We
tested several spatial discretization scenarios and achieved
the highest accuracy by applying a spatial resolution of
10−6 cm for the first node and gradually increasing theresolution with depth up to a spatial resolution of 1 cm
for the last node. We also adjusted the initial time steps
in order to achieve the convergence during numericalsimulations.
The internal interpolation tables of HYDRUS-1D weredisabled to increase the precision of simulated results. A
modified VGM model with an air-entry value of −2 cm
was used for soils with the n parameter smaller than 1.2.
This has been recommended in the literature (Schaap &van Genuchten, 2006; Vogel & Cislerova, 1988; Vogel, vanGenuchten, & Cislerova, 2000) in order to avoid an unrealistically sharp decrease of the VGM hydraulic conductivity
function close to saturation.

RAHMATI et al.

Vadose Zone Journal

9 of 22T A B L E 2 Average soil hydraulic parameters for the van Genuchten (1980) model for 12 USDA textural classes (Carsel & Parrish, 1988)
and the sorptivity (S) value obtained from the horizontal infiltration simulation
Textural class

θr

θs

θi

cm3 cm−3

α

n

m

cm−1

Ks

S

cm h−1cm h−1/2

β

Clay

0.068

0.380

0.008

1.09

0.083

0.271

0.20

1.02

1.92

Clay loam

0.095

0.410

0.019

1.31

0.237

0.150

0.26

1.46

1.58

Loam

0.078

0.430

0.036

1.56

0.359

0.088

1.04

2.20

1.27

Loamy sand

0.057

0.410

0.124

2.28

0.561

0.057

14.6

6.22

0.80

Sand

0.045

0.430

0.145

2.68

0.627

0.045

29.7

9.23

0.60

Sandy clay

0.100

0.380

0.027

1.23

0.187

0.170

0.12

0.79

1.70

Sandy clay loam

0.100

0.390

0.059

1.48

0.324

0.111

1.31

1.61

1.36

Sandy loam

0.065

0.410

0.075

1.89

0.471

0.066

4.42

3.84

0.99

Silt

0.034

0.460

0.016

1.37

0.270

0.090

0.25

1.35

1.50

Silt loam

0.067

0.450

0.020

1.41

0.291

0.104

0.45

1.66

1.44

Silt clay

0.070

0.360

0.005

1.09

0.083

0.266

0.02

0.35

1.92

Silty clay loam

0.089

0.430

0.010

1.23

0.187

0.197

0.07

0.53

1.70

Note. θi , θs , and θr are initial, saturated, and residual water contents, respectively; α, n, and m are parameters of van Genuchten’s (1980) model; Ks is saturated
hydraulic conductivity, S is soil sorptivity, and β is an infiltration constant defined by Haverkamp et al. (1994)

3.1.2

Soil sorptivity inference

In order to test the estimated S values obtained using
different methods analyzed in this paper, we needed toknow their true values for each examined soil. Rather thanusing simplified equations to estimate S, we integrated λ(θ)
curves numerically simulated by HYDRUS-1D for infiltration without gravity and calculated the S values using thefollowing equation (Philip, 1969a):
θs

𝑆=

∫θ

λ(θ)dθ

(29)iwhere θi [L3 L−3 ] and θs [L3 L−3 ] are the initial and saturated volumetric water contents, respectively, and λ(θ)[L T−1/2 ] is the so-called Boltzmann variable defined in
Equation 28. The obtained values of S are reported inTable 2.

3.2

Model comparisons

We compared four different methods, including the characteristic time method (CTM) as proposed above, theSharma et al. (1980) method (SH), and a nonlinear curvefitting method (Bonell & Williams, 1986; Bristow & Savage,
1987; Marquardt, 1963; Vandervaere et al., 2000) with two(CF2, Equation 5) and three-term (CF3, Equation 7) equations, regarding their ability to estimate Ks and S.
In the case of SH method, we used simulated and measured infiltration data lasting for <30 min to estimate S.

The 30 min is selected as a critical√ time to ensure the lininitial time steps.
ear relationship between I and 𝑡 at √Therefore, the I was linearly plotted vs. 𝑡 for data points
obtained between start and 30 min. Then, the slope provides S when intercept is set to zero. For Ks predictions
by SH method, we set it to be Ks = (Iend – Iend – 1 )/(tend –tend – 1 ). In the case of CF2 and CF3, the “lsqnonlin” function is used to predict S and Ks .
We truncated simulated infiltration curves to obtain
infiltration curves of different durations from 5 min to
10 h, corresponding to a typical infiltration measurement
window derived from the Soil Water Infiltration Global(SWIG) database (Rahmati, Weihermüller, Vanderborght,
et al., 2018; Rahmati, Weihermüller, & Vereecken, 2018).
We used a variable infiltration duration in order to assess
the capability of CTM in finding the characteristic time andinfiltration independently of the duration of the infiltration experiments.
In addition to the analysis mentioned above, we also
used synthetic data with measurement noise as well as
experimental data selected from SWIG database (Rahmati,
Weihermüller, Vanderborght, et al., 2018; Rahmati, Weihermüller, & Vereecken, 2018) to verify the results obtained
from perfect and error-free infiltration data obtained fromnumerical simulation by HYDRUS-1D.
In the case of data with measurement noise, we first
assumed that the water level in the system can be measured with an accuracy of 1 mm, in accordance with mostexperimental devices and protocols. We account for thisby rounding up the simulated cumulative infiltration data
expressed in centimeters to one digit after the comma. We

10 of 22

RAHMATI et al.

Vadose Zone Journal

then selected the infiltration times corresponding to eachmillimeter increment until the whole cumulative infiltration. The cumulative infiltration (I) curve was then differentiated with respect to time to obtain the infiltration rate(i). Then we added random noise on the infiltration rates
using following equation:𝐢noised = 𝐢 + 𝑟𝛔

3.3
Evaluation of the estimatedhydraulic parameters
The accuracy of the CTM and other selected methods wasevaluated using the RMSE and Nash and Sutcliffe (1970)
efficiency, E, criteria between measured and predicted Ks
and S values:√(30)
RMSE =

where inoised is the noised infiltration rate vectors, r is a random number between −1 and 1 picked from normal distribution, and σ is the vector of standard deviations. We calculated the standard deviation as a function of the nominal
values of the infiltration rate i, leading to
𝛔 = 𝑒𝐢(31)

where e is the magnitude of error considered in analysis.
We used three different values of 0.05, 0.1, and 0.2 indicating an error of 5, 10, and 20%, respectively. Since the rangeof r in Equation 30 is −1 and 1, a value of e = 1 will impose
100% error in each point. Therefore, the parameter e is used
to control the maximum allowed error. Finally, after calculating the infiltration rates containing noise (inoised ), we
calculated the noised cumulative infiltration (Inoised ) curveby integration:
𝐼noised (𝑗) = 𝐼noised (𝑗 − 1) + 𝐢noised (𝑗)[𝑡(𝑗) − 𝑡(𝑗 − 1)],
𝑗 = 1, 2, … , 𝑛(32)
By applying the above procedure, we are able topropagate the error and generate more “realistic” timeseries that correspond better with the measurement process. The error was first defined for the infiltration rate
and then propagated to the cumulative infiltration. Thisapproach only considers the measurement error and not
errors caused by violation of the underlying theoretical assumptions such as the presence of soil layering,
nonuniform infiltration, and infiltration in (dead-end)
macropores.
In the case of experimental data, we selected 1D infiltration data available in the SWIG database obtained
from infiltration experiments using a zero-water potential
imposed at surface. In overall, 614 infiltration data wereselected for final analysis. As no independent measurements of Ks and S are reported in SWIG for the selected
data, we estimated these values using the QEI formulation(Latorre et al., 2015), which we consider as a benchmark
for comparison with the estimates from other approaches.
It should be noted the QEI is an implicit formulation thatslows down the computations and then the predictions of
Ks and S.

∑(

𝑋m − 𝑋p

)2

𝑛

)2
∑(
𝑋m − 𝑋p
𝐸 =1− ∑
2
(𝑋m − 𝑋̄ m )

(33)

(34)

where Xm and Xp are the logarithmic values of known
and predicted parameters (S and/or Ks ), respectively. In
the case of the SWIG database, the predictions provided
by the QEI model were considered as the known values.
The values of RMSE near zero denote a great accuracy,
whereas E takes values in the interval [−∞, 1]. A value ofE equal to unity shows a perfect match and an E value ofzero means that the average of the observed values is as
good as predicted values. Negative E values mean that themodel is worse for prediction than taking only the mean of
the observed values. E values higher than 0.9 show that theslope between observed and predicted values is near unity.

4

RESULTS AND DISCUSSION

4.1
CTM for infiltration data in thetransient and steady-state infiltration
regimes
In order to evaluate the accuracy of the CTM for infiltration data in transient and steady-state infiltration regimes,
we applied it on all 12 synthetic soils with infiltration timesof 240 h. We evaluated the procedure by comparing theestimated Ks and S values obtained from CTM with theirknown values and calculated the RMSE and E metrics(Figure 5). The proposed procedure showed very highaccuracy in predicting both S and Ks , respectively, with
RMSE values of 0.04 cm h−1/2 and 0.05 cm h−1 , and E values of 0.992 and 0.999 (Figure 5).
As seen from Table 3, all simulated infiltration curves
except for silty clay soil predicted tchar values with ω valueof 0.5. The silty clay soil had ω value of 0.246, indicatingthat tgrav for this soil is higher than the duration of the simulated infiltration curves (532 vs. 240 h, Table 3). Therefore,
excluding silty clay soil, the CTM-I is used to predict Ks forRAHMATI et al.

Vadose Zone Journal

11 of 22F I G U R E 5 Comparison between known and predicted values of soil sorptivity (S) and saturated hydraulic conductivity (Ks ) of examinedsoils using the characteristic time method (CTM) on simulated infiltration curves with 240-h duration. E is the Nash–Sutcliffe efficiency

T A B L E 3 Applied values of ω (the contribution of the gravity
component to the infiltration process), as well as obtained
characteristic (tchar ) and gravity (tgrav ) times and the soil-dependent
shape factor of α obtained from application of characteristic time
method (CTM) over simulated infiltration curves with 240-h
duration
Soil

ω

h
Clay

0.500

tchar

α
−1

−0.009

tgrav
h

81

81

Clay loam

0.500−0.007

93

93

Loam

0.500

−0.03216

16

Loamy sand

0.500

−0.600

0.72

0.72

Sand

0.500

−1.0320.48

0.48

Sandy clay

0.500

−0.006124

124

Sandy clay loam

0.500

−0.068

5.67

5.67

Sandy loam

0.500

−0.174

2.60

2.60

Silt

0.500

−0.006104

104

Silt loam

0.500−0.01150

50Silt clay

0.242−0.001

239

532

Silty clay loam

0.500

−0.003

212

212

Avg.

0.500

−0.16277

102

all soils. In the case of silty clay soil, CTM-Ks is used to provide accurate prediction for Ks . The differences between
CTM-I and CTM-Ks is investigated in Section 4.2 applying different infiltration times where the attainment of thesteady-state regime in infiltration data is not guaranteed
and the use of CTM-Ks is necessary.
Note that the reported tgrav values (Table 3) show thatcoarse-textured soils have low values of tgrav , whereas finetextured soils have higher values of tgrav , which is in linewith the physics of water infiltration into soils (with morecapillarity-driven flow in fine-textured soils).

The shape factor α values were easily obtained withexcellent fits to Equation 13 and considering the full infiltration curves with both the transient and the steady-state
infiltration regimes, except for silty clay soil. The corresponding values are reported in Table 3. It is obvious from
Table 3 that α values are soil dependent, with smaller values for coarser textured soils. There are very high correlations between α values and known values of S (with Pearson correlation of −.975) and Ks (with Pearson correlationof −.997). Extra research is needed to investigate the natureof the α parameter which is beyond the scope of the current
paper.

4.2
Evaluating the CTM as a function of
the expriment duration
In order to evaluate the accuracy of CTM in estimating Ks and S from infiltration measurements that do
not necessarily reach the steady-state regime, we applied
the method on truncated infiltration curves with different times. We limited the simulated infiltration time tobe <24 h because infiltration measurements are usually
performed in shorter times, typically up to a few hours.
We therefore expect that the steady-state regime is not
obtained in most real-world experiments requiring the useof the CTM-Ks procedure.
The outcomes demonstrate very high accuracy of the
CTM in the case of S predictions showing an average E
value of about unity (0.992 ± 0.0001) and average RMSE
value of 0.041 ± 0.0002 cm h−1/2 , and that the obtained
accuracy does not depend on the duration of infiltration,
which is already understandable from very small value ofstandard deviations (with SD values of 0.0001 and 0.0002
for E and RMSE, respectively). The reason why CTM
results in very high and time-independent accuracy in the12 of 22

Vadose Zone Journal

RAHMATI et al.

F I G U R E 6 Variations of soil-averaged RMSE and Nash–Sutcliffe efficiency (E) criteria of predicted saturated hydraulic conductivity (Ks )over different infiltrations durations using the characteristic time method coupled with the iterative procedure (CTM-I) and the CTM coupledwith the Ks procedure (CTM-Ks )case of S predictions is that CTM relies on the contribution of the first component for S predictions, which isalready dominant at initial time steps. Therefore, as soon
as enough information to retrieve S from infiltration data
is provided at initial time steps, the further times steps will
not be necessary for calculation.
Contrary to S, CTM-I fails to predict Ks accurately whenusing a ω value smaller than 0.5 (tchar < tgrav )(Figure 6).
Figure 6 shows that for ω values between 0 and 0.5,
Ks prediction was completely dependent on infiltration
time, leading to inaccurate estimates for shorter infiltration times and moderately accurate estimates for longer
infiltration times. On the other hand, using an iterative
procedure (CTM-I), accurate predictions of Ks in finetextured soils required infiltration times >24 h, whichmakes the method impractical. This happens because theiterative procedure is developed based on the calculationof the contributions of the sorptivity and gravity terms to
the infiltration process and their usage for S and Ks predictions. In fact, the procedure uses the contribution of
the sorptivity term to predict S and then the contribution of the gravity terms to predict Ks . Thus, with infiltration data lasting for short times (couple of minutes),
the variation of sorptivity contribution occurs around one
while at the same time the variation of the contributionof the gravity terms occurs around zero (see Figure 1c).
Particularly in the case of fine-textured soils, the contribution of gravity terms stays around zero for a longer
period, preventing the procedure from accurate predictionof the Ks . In other words, the iterative procedure requires

obtaining infiltration data for a steady state, which may
take a very long to attain (see Lassabatere, AnguloJaramillo, Soria-Ugalde, Šimůnek, and Haverkamp (2009)
for examples of times required to reach steady state).
To overcome this drawback, CTM-Ks was developed to
improve Ks predictions for the case of short-time infiltration data. We demonstrate its efficiency for the case of infiltration durations varying between a couple of minutes and
24 h, covering the typical infiltration measurement window (5 min to 10 h for most experiments; Angulo-Jaramillo
et al., 2000; Rahmati, Weihermüller, Vanderborght, et al.,2018; Rahmati, Weihermüller, & Vereecken, 2018).
Similar to CTM-I, the CTM-Ks also showed a timedependent, but high, accuracy for Ks predictions(Figure 6), showing an average E value of 0.921 ±
0.101 and average RMSE value of 0.288 ± 0.205 cm h−1 .
Exclusion of the infiltration durations shorter than 15 min(tmax < 15 min) resulted in even higher accuracy and lessdependency on infiltration duration in Ks predictions,
showing an average E value of 0.975 ± 0.026 and average
RMSE value of 0.179 ± 0.084 cm h−1 . Thus, we recommend using the CTM-Ks to predict Ks over infiltration data
lasting more than 15 min. The recommended duration
is within the typical time window of the infiltration
measurement (5 min to 10 h). Therefore, when CTM is
mentioned hereafter, it refers to use of the CTM-Ks in thecase of Ks predictions to ensure accurate prediction. Notethat to be applied, CTM-Ks needs a precise estimation of S
provided by the previous run of CTM-I, which is always
the case.

RAHMATI et al.

Vadose Zone Journal

13 of 22F I G U R E 7 RMSE and Nash–Sutcliffe efficiency (E) criteria obtained between known and predicted saturated hydraulic conductivity (Ks )values using the characteristic time method (CTM) with a constant (β = 0.6) and soil-dependent β parameter over simulated infiltration curvesof different durations

4.3
The impact of constant andsoil-dependent β on the results of the CTM
Since the CTM is developed based on the approximate
expansion of Haverkamp et al. (1994) model involving
three input parameters (Ks , S, and β), in this section, we
aimed to evaluated the remaining β parameter’s effects
on predictions of S and Ks . However, we excluded S fromthis analysis because as seen from Equation 16, β has noeffect on S predictions. Therefore, we only illustrates theresults for Ks predictions in this section. Lassabatere et al.
(2009) showed an excellent agreement between the numerically generated infiltration data and the model developed
by Haverkamp et al. (1994), provided that the values of
β are adapted to the soil type, ranging between values of
0.3 and 1.7 for sand and silt soils, respectively. However,more recent studies (Latorre et al., 2015; Rahmati et al.,2019) proved that the choice of the value of β had only aslight effect on the estimates of S and Ks when using Equation 3 and that the default value proposed by Haverkamp
et al. (1994) (i.e., 0.6) could be considered as a good option.
Therefore, in order to check the effect of a constant or soildependent β parameter on the efficiency of CTM in estimating Ks and S, we used a constant β value of 0.6 anda soil-dependent β parameter for time-variant infiltration
curves (0–24 h). In the case of soil-dependent β, we used
Equation 4 to estimate β before deriving Ks and S with the
CTM method.

In the case of a soil-dependent β parameter, the RMSE
and E criteria were slightly higher and lower, respectively,
compared with the RMSE and E values obtained with aconstant β (Figure 7), showing average RMSE value of
0.316 vs. 0.288 cm h−1 and E value of 0.909 vs. 0.921. One
might expect better accuracy for a soil-dependent β parameter compared with a constant β. However, this does not
happen because, as Lassabatere et al. (2009) has alreadydiscussed, Equation 4 does not necessarily provide thebest values for soil-dependent β parameter, particularly inthe case of fine-textured soils. We conducted the Mann–
Whitney–Wilcoxon test to examine whether the use of aconstant or soil-dependent β parameter led to differencesin Ks estimates. The result (with a p value of .75) showedno significant difference between these two approaches,
despite significant differences between the computed anddefault values for β. We therefore recommend the use of
CTM with a constant β value of 0.6. Previously Latorre et al.
(2018) also showed that the quasi-exact implicit solutionof Haverkamp et al. (1994) and its approximate expansions
are less sensitive to β.

4.4Models comparison using error-freeand noised synthetic infiltration curves
To compare CTM with other methods, we examined thegain in accuracy of CTM in estimating S and Ks compared14 of 22

Vadose Zone Journal

RAHMATI et al.

F I G U R E 8 Soil-averaged RMSE between logarithms of known and predicted saturated hydraulic conductivity (Ks ) (right) and soil sorptivity (S) (left) values using the characteristic time method (CTM), Sharma method (SH), and nonlinear curve fitting methods of two- (CF2)and three-terms (CF3) equations over original and noised synthetic datawith SH, CF2, and CF3. To do this, we used the synthetic
infiltration curves with a duration between 15 min and 10 h,corresponding to typical infiltration experiments, as discussed above. In addition to the original error-free data, we
also added some random noises over synthetic infiltration
data to make them more representative of field conditions
and to see the consequences on CTM estimations. Readers
are invited to refer Section 3 for the details of noise addition
on infiltration data. Figure 8 illustrates the RMSE values
obtained between known and predicted values of Ks and S
over different times, and Table 4 reports the average values
for each method.
Considering the error-free data, the CTM method
showed the best performance in estimating both Ks andS showing average RMSE values of 0.204 ± 0.086 cm h−1
and 0.040 ± 0.0001 cm h−1/2 , respectively (Table 4). The
SH method showed the lowest performance in estimating
the S (with average RMSE values of 0.118 ± 0.011 cm h−1/2 ),whereas it ranked second in estimating the Ks showing an
RMSE value of 0.379 ± 0.178 cm h−1 (Table 4). The CF2 andCF3 methods showed the lowest performance in the caseof Ks estimates, showing RMSE values>1.2 cm h−1 . In thecase of S estimates, the CF3 was competitive with CTM,
showing RMSE values of 0.034 ± 0.020 cm h−1/2 vs. 0.040 ±0.001 cm h−1/2 for CTM (Table 4). The CF2 ranked third in
the case of S estimations, showing an RMSE value of 0.106± 0.062 cm h−1/2 (Table 4). Therefore, the results revealthat although SH and CF3 methods show good accuracy

in estimation of Ks and S, respectively, they fail in accurate prediction of the other parameter. This happens while
the CTM shows very high accuracy in simultaneous estimation of both S and Ks parameters. The accuracy of bothpredicted parameters using CTM is nearly time independent, which makes the CTM more practical to apply overinfiltration data lasting for a couple of minutes to dozens
of hours.
In addition to error-free data, overall, 5, 10, and 20%
maximum errors were added on infiltration data, and thenoise impact on S and Ks predictions was analyzed. The
results (Figure 8, Table 4) revealed that in the case of Ks
predictions, the accuracy of the CTM even with noisy data
was considerably higher than that of CF2 and CF3, as well
as SH, showing average RMSE values of 0.341 cm h−1 vs.
0.407 cm h−1 in the case of SH and >2.9 cm h−1 in the caseof CF2 and CF3 methods (Table 4). For S predictions, CTM
and CF3 showed similar accuracies, with average RMSE
values of 0.07 cm h−1/2 (Table 4). The average RMSE values of CF2 and SH are nearly double those of CTM and
CF3 (Table 4).
A simple comparison between the error-free and noised
data also revealed that CTM works nicely, even with 20%
noise on data, showing a slight increase in RMSE values
due to noise addition. In contrast with CTM, the accuracy of CF2 and CF3 considerably degrades with increasein noise magnitude (particularly in the case of Ks predictions), whereas the average RMSE value increases from

RAHMATI et al.

T A B L E 4 Average RMSE between logarithms of known andpredicted saturated hydraulic conductivity (Ks ) and soil sorptivity(S) values using the characteristic time method (CTM), Sharma
method (SH), and nonlinear curve fitting method of two- (CF2) andthree-term (CF3) equations over original and noised synthetic data
for all soils
Data
Error-free
data

Noised data
(5%)

Noised data
(10%)

Noised data
(20%)

15 of 22

Vadose Zone Journal

MethodKs

Scm h−1cm h−1/2

T A B L E 5 Average RMSE and Nash and Sutcliffe (1970)
efficiency (E) criteria between logarithms of known and predictedsaturated hydraulic conductivity (Ks ) and soil sorptivity (S) values
using the characteristic time method (CTM), Sharma method (SH),and nonlinear curve fitting method of two- (CF2) and three-term(CF3) equations over experimental data of the Soil Water
Infiltration Global (SWIG) database
Criterion

MethodKs

Scm h−1cm h−1/2

CTM

0.204 ± 0.086

0.040 ± 0.001
CTM

0.298

0.307

SH

0.379 ± 0.178

0.118 ± 0.011

SH

0.202

0.184

CF2

1.296 ± 0.962

0.106 ± 0.062

CF2

1.685

1.354

CF3

1.391 ± 1.079

0.034 ± 0.020

CF3

1.769

0.279

CTM

0.335 ± 0.177

0.068 ± 0.007

CTM

0.891

0.736

SH

0.938

0.885

CF2

−2.490−4.139CF3−2.844

0.782

RMSE

SH

0.411 ± 0.232

0.144 ± 0.009

CF2

2.925 ± 0.694

0.136 ± 0.057

CF3

3.191 ± 0.851

0.067 ± 0.012

CTM

0.340 ± 0.182

0.067 ± 0.008

SH

0.408 ± 0.229

0.144 ± 0.008

CF2

2.897 ± 0.933

0.135 ± 0.056

CF3

3.299 ± 1.041

0.066 ± 0.015

CTM

0.349 ± 0.173

0.066 ± 0.006

SH

0.401 ± 0.232

0.143 ± 0.009

CF2

2.910 ± 0.758

0.128 ± 0.058

CF3

3.252 ± 0.671

0.066 ± 0.013

around 1.3 cm h−1 up to 3.2 cm h−1 . Although SH is lessaccurate than CTM, the results were shown to be stable byincrease in magnitude of the noise. In the case of S prediction, the RMSE value over noised data is nearly doubled compared with that of error-free data in both CTM
and CF3, regardless of the magnitude of the noise.

4.5
Evaluating the CTM overexperimental data
Up to this point in the study, the CTM method was tested
and verified over synthetic data obtained from HYDRUS1D. Experimental data obtained in field campaigns areusually noisy and not perfect. Although we added somerandom noises over those synthetic data to evaluate themethod performance over the noisy data as well, such anoise addition will not represent other important error
sources such as soil layering, nonuniform infiltration, andinfiltration in (dead-end) macropores. Therefore, a performance evaluation over experimental data is yet required.
However, experimental data usually lack the true values ofE

S and Ks , making it difficult to evaluate the accuracy of themethods (i.e., the comparison of estimated values against
targeted real values). On the other hand, the measured values of Ks mostly occur in disturbed or undisturbed soil
samples, which are mostly not representative of field conditions under infiltration measurements. Therefore, the Ks
values predicted from infiltration data rarely ever correspond to those values obtained from small soil samples.
As a conclusion of the above, we used field experimental data taken from the SWIG database (Rahmati,
Weihermüller, Vanderborght, et al., 2018; Rahmati, Weihermüller, & Vereecken, 2018) to verify the feasibility of
the application of the proposed method to real data andthe potential impact of experimental noise. In order tocover the lack of the true and benchmark values of the
S and Ks , we used QEI formulation of the Haverkamp
et al. (1994), which is valid for entire time to be solved
numerically to retrieve the benchmark values of S and Ks .
Then, we compared the CTM predicted values of S and Ks ,as well as those predicted by SH, CF2, and CF3, with thebenchmark values. Although the QEI-retrieved values of
S and Ks are considered as the most accurate predictionsof parameters, it is impractical for routine use because thenumerical calculations of the parameters last for coupleof hours in most cases.
As shown in Figure 9 and summarized in Table 5,both S and Ks values predicted by CTM are considerablycomparable with benchmark values. The CTM showed a
RMSE value of 0.202 cm h−1 in the case of Ks predictions
and 0.307 cm h−1/2 in the case of S prediction. The SH
method was even more robust compared with CTM, showing RMSE values of 0.202 cm h−1 and 0.184 cm h−1/2 inthe case of Ks and S predictions, respectively. Although,the SH method shows slightly better robustness than CTM,

16 of 22

RAHMATI et al.

Vadose Zone Journal

F I G U R E 9 Logarithms of predicted saturated hydraulic conductivity (Ks ) and sorptivity (S) values using the characteristic time method(CTM) and Sharma et al. (1980) (SH) method, as well as curve fitting methods of two- (CF2) and three-term (CF3) equations vs. logarithms of
their most accurate estimates using quasi-exact implicit (QEI) solution of Haverkamp et al. (1994)

it still requires arbitrary selection of the data for accurate
prediction of S. When using the SH method, one first needs
to decide up to which point of measured data the relationship between I and t1/2 is linear, and then the selected partshould be used for S predictions. This adds some arbitrariness to S predictions, leading to uncertainty in estimatedvalues depending on the choice made. In our calculations,after several examinations, we finally decided that whena critical time of 30 min is used to separate the first partof data for S predictions, it results in the highest accuracy of the S predictions overall. However, the use of CTM
omits the necessity of such arbitrariness in selection, showing very high robustness compared with QEI formulation.
Contrary to CTM, SH, CF2, and CF3 show unacceptableaccuracy for Ks predictions. Although CF3 shows an accuracy comparable with CTM and SH in the case of S predictions, CF2 fails in this case too. It must be noted that if one
excludes the outliers from CF2 and CF3 predictions, theaccuracy might be close to that obtained from CTM andSH (Figure 9). However, the outliers reveal that CF2 andCF3 might have physically unconstrained fitting over infiltration data in several cases, which challenges the reliable
use of them for S and Ks predictions. This happens while
CTM provides physically constrained fitting of data. This
could be the case for SH method as well, only in a condition that the separation of the linear part of the infiltration
data for S predictions is done correctly.

4.6

Practical interest of the CTM

The CTM method has several practical advantages since it
overcomes several shortcomings encountered when S and
Ks are estimated using non-linear curve-fitting methods. It
is built on the framework of the approximate expansions of
the quasi-exact implicit model developed by Haverkamp
et al. (1994), providing the most general and accurate
description of 1D water infiltration compared with classical choices based on, for example, Philip (1957) or empiricalequations (Holtan, 1961; Horton, 1941; Kostiakov, 1932). In
addition, CTM addresses the right physics (i.e., the physics
of water infiltration under the fixed water pressure head at
the surface, as is the case for regular devices and protocols).
Nonlinear curve fitting methods require an appropriate
selection of data, which is complicated since the characteristic time that divides the transient and steady statesis unknown and depends on unknown values of S and
Ks . To alleviate this problem, iterative procedures wereused to select a priori data, estimate values of S and Ks ,and then compute the characteristic time to validate theselection posteriori, as done in the Beerkan estimationof soil transfer (BEST) methods (Angulo-Jaramillo et al.,2019). Haverkamp et al. (1994) developed an implicit formulation, which is valid for all times and can be used
to invert experimental data (Fernández-Gálvez, Pollacco,
Lassabatere, Angulo-Jaramillo, & Carrick, 2019; Latorre

RAHMATI et al.

et al., 2015). However, its use is quite complicated due toits implicit feature, and thus direct models (dedicated toeither transient or steady-state regimes) are still preferred.
The fit itself requires paying attention to several additional crucial points. The fit involves the choice of an
objective function that is barely developed. The objectivefunctions can be of two main types: distance-based or
weak form-based functions, with the distance-based functions being the most popular in soil science and hydrology(Guinot, Cappelaere, Delenne, & Ruelland, 2011). Several
choices are available among the distance-based objectivefunctions, such as the Nash–Sutcliffe efficiency, the sum ofsquared errors (SSE), the RMSE, the mean absolute error(MAE), the volumetric efficiency (VE), or the generalized Nash–Sutcliffe efficiency (GNSE) (Guinot et al., 2011).
However, little attention is paid to the choice of the rightobjective function among these options. In addition, thesetypes of functions are known to introduce local minima inthe model response, resulting in a strong dependence on
the choice of the initial values. The model response andobjective functions should be tested to avoid equifinality
and nonuniqueness (Pollacco et al., 2013).
The CTM method avoids all these shortcomings outlined
above. Its sole complexity involves the use of an iterative
process that is necessary since the knowledge of the transition time between transient and steady regime depends
on unknown values of S and Ks . Despite this difficulty, nofits are required, and the computation involves only simple
algebraic operations. Its adaptation to all types of datasets(encompassing a part of the transient state or including
both transient and steady states) makes it a useful tool. In
addition to its advantages, our study demonstrates its efficiency in reaching its goal. All these aspects make the CTM
method a promising method for the estimation of S and Ks .

Vadose Zone Journal

17 of 22(1980) (SH). We used HYDRUS-1D to simulate the synthetic infiltration curves for soils with different textures
ranging from sandy to clayey soils. We found that CTM provided accurate predictions of Ks and S compared with theother methods when tested against numerically generated
infiltration curves with a duration of couple of minutes todozens of hours. The accuracy of the CTM was better thanother examined methods. The most important advantage
of CTM is that it needs no prior knowledge about the time
domain validity of the applied equation(s). In contrast, theother methods failed when applied to time-variant infiltration curves. The CTM was found to be applicable to any
infiltration duration within a typical time window used
for infiltration experiments (5 min to 10 h), and even up
to 240 h. The accuracy of Ks and S predictions using CTM
showed no dependency on the duration of infiltration. We
also evaluated the performance of the CTM in a comparison with other methods over noised synthetic data, as
well as experimental data selected from the SWIG database(Rahmati, Weihermüller, Vanderborght, et al., 2018; Rahmati, Weihermüller, & Vereecken, 2018), to verify the usageof the CTM with real noised data. The results showedthat even with noised and experimental data, the CTM
still was better than CF2 and CF3, showing considerably
higher accuracy, particularly in the case of Ks predictions.
The CTM method was also comparable with the classicalmethod of SH owing to the advantage of no arbitrariness inselection of data for accurate predictions of S and Ks . Further research is needed to corroborate the estimated valuesof S and Ks using CTM with independently determined values of S and Ks . Also, improved estimates of S and Ks may
prove to be valuable in constraining estimates of the van
Genuchten (1980) model’s parameters, such as α and n.
CONFLICT OF INTEREST
The authors declare no conflict of interest.

5

CONCLUSIONSIn this paper, we proposed a new method, the characteristic time method (CTM), based on the concept of characteristic time (tchar ) to estimate the soil sorptivity S andthe saturated hydraulic conductivity Ks in a physically constrained manner using the two- and three-term approximate expansion developed by Haverkamp et al. (1994) fromcumulative infiltration measurements. The CTM method
requires no knowledge of the infiltration duration, forwhich Philip’s (1957) infiltration equation and the two- or
three-term equations of Haverkamp et al. (1994) are valid.
We analyzed the accuracy of CTM using simulated andexperimental infiltration curves and by comparing CTM
with methods typically used in the literature, including
nonlinear curve fitting of two- (CF2) and three-term (CF3)equations, and the methods proposed by Sharma et al.

ORCID
Mehdi Rahmati https://orcid.org/0000-0001-5547-6442
Laurent Lassabatere https://orcid.org/0000-0002-86255455
Harry Vereecken https://orcid.org/0000-0002-8051-8517
REFERENCES
Angulo-Jaramillo, R., Bagarello, V., Di Prima, S., Gosset, A., Iovino,M., & Lassabatere, L. (2019). Beerkan estimation of soil transfer
parameters (BEST) across soils and scales. Journal of Hydrology,
576, 239–261. https://doi.org/10.1016/j.jhydrol.2019.06.007
Angulo-Jaramillo, R., Vandervaere, J. P., Roulier, S., Thony, J. L.,
Gaudet, J. P., & Vauclin, M. (2000). Field measurement of soil surface hydraulic properties by disc and ring infiltrometers: A review
and recent developments. Soil and Tillage Research, 55, 1–29.
https://doi.org/10.1016/S0167-1987(00)00098-2

18 of 22

Vadose Zone Journal

Beven, K., & Freer, J. (2001). Equifinality, data assimilation, anduncertainty estimation in mechanistic modelling of complex
environmental systems using the GLUE methodology. Journalof Hydrology, 249, 11–29. https://doi.org/10.1016/S0022-1694(01)00421-8
Bonell, M., & Williams, J. (1986). The two parameters of the Philip
infiltration equation: Their properties and spatial and temporal heterogeneity in a red earth of tropical semi-arid Queensland. Journal of Hydrology, 87, 9–31. https://doi.org/10.1016/00221694(86)90112-5
Bristow, K. L., & Savage, M. J. (1987). Estimation of parameters
for the Philip two-term infiltration equation applied to field soil
experiments. Soil Research, 25, 369–375. https://doi.org/10.1071/
SR9870369
Carsel, R. F., & Parrish, R. S. (1988). Developing joint probability distributions of soil water retention characteristics. Water Resources
Research, 24, 755–769. https://doi.org/10.1029/WR024i005p00755
Fernández-Gálvez, J., Pollacco, J., Lassabatere, L., Angulo-Jaramillo,
R., & Carrick, S. (2019). A general Beerkan estimation of soiltransfer parameters method predicting hydraulic parameters of
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
Green, W. H., & Ampt, G. A. (1911). Studies on soil physics.
Journal of Agricultural Science, 4, 1–24. https://doi.org/10.1017/
S0021859600001441
Guinot, V., Cappelaere, B., Delenne, C., & Ruelland, D. (2011).
Towards improved criteria for hydrological model calibration:
Theoretical analysis of distance- and weak form-based functions.
Journal of Hydrology, 401, 1–13. https://doi.org/10.1016/j.jhydrol.
2011.02.004
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.
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, 2931–2935. https://doi.org/10.1029/94WR01788
Holtan, H. N. (1961). Concept for infiltration estimates in watershedengineering. Washington, DC: USDA-ARS.
Horton, R. E. (1941). An approach toward a physical interpretation of infiltration-capacity. Soil Science Society of America Journal, 5, 399–417. https://doi.org/10.2136/sssaj1941.036159950005
000C0075x
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 Transactions of 6th Committee International Society of Soil Science (pp.
7–21). International Society of Soil Science.
Lassabatere, L., Angulo-Jaramillo, R., Soria Ugalde, J., Cuenca, R.,
Braud, I., & Haverkamp, R. (2006). Beerkan estimation of soiltransfer parameters through infiltration experiments: BEST. Soil
Science Society of America Journal, 70, 521–532. https://doi.org/10.
2136/sssaj2005.0026
Lassabatere, L., Angulo-Jaramillo, R., Soria-Ugalde, J., Šimůnek, J.,& Haverkamp, R. (2009). Numerical evaluation of a set of analytical infiltration equations. Water Resources Research, 45(12).
https://doi.org/10.1029/2009WR007941

RAHMATI et al.

Latorre, B., Moret-Fernández, D., Lassabatere, L., Rahmati, M.,
López, M. V., Angulo-Jaramillo, R., . . . Jiménez, J. J. (2018). Influence of the β parameter of the Haverkamp model on the transientsoil 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 fromdisc 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
Marquardt, D. W. (1963). An algorithm for least-squares estimationof nonlinear parameters. Journal of the society for Industrial andApplied Mathematics, 11, 431–441. https://doi.org/10.1137/0111030
Moret-Fernández, D., Blanco, N., Martínez-Chueca, V., & Bielsa, A.
(2013). Malleable disc base for direct infiltration measurements
using the tension infiltrometry technique. Hydrological Processes,
27, 275–283. https://doi.org/10.1002/hyp.9206
Moret-Fernández, D., Latorre, B., & Angulo-Martínez, M. (2017).
Comparison of different methods to estimate the soil sorptivity
from an upward infiltration curve. Catena, 155, 86–92. https://doi.
org/10.1016/j.catena.2017.03.009
Nash, J. E., & Sutcliffe, J. V. (1970). River flow forecasting throughconceptual models part I: A discussion of principles. Journalof Hydrology, 10, 282–290. https://doi.org/10.1016/0022-1694(70)90255-6
Parlange, J. Y., Lisle, I., Braddock, R. D., & Smith, R. E. (1982). The
three-parameter infiltration equation. Soil Science, 133, 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, 345–358.
Philip, J. R. (1969a). Moisture equilibrium in the vertical in swellingsoils. I. Basic theory. Soil Research, 7, 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
Pollacco, J. A. P., Nasta, P., Soria-Ugalde, J. M., Angulo-Jaramillo, R.,
Lassabatere, L., Mohanty, B. P., & Romano, N. (2013). Reduction offeasible parameter space of the inverted soil hydraulic parameter
sets for Kosugi model. Soil Science, 178, 267–280.
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., Weihermüller, L., Vanderborght, J., Pachepsky, Y. A.,
Mao, L., Sadeghi, S. H., . . . Vereecken, H. (2018). Development andanalysis of the Soil Water Infiltration Global database. Earth System Science Data, 10, 1237–1263. https://doi.org/10.5194/essd-101237-2018
Rahmati, M., Weihermüller, L., & Vereecken, H. (2018) Soil Water
Infiltration Global (SWIG) database. PANGAEA. https://doi.org/
10.1594/PANGAEA.885492
Richards, L. A. (1931). Capillary conduction of liquids throughporous mediums. Physics, 1, 318–333. https://doi.org/10.1063/1.
1745010
Schaap, M. G., & van Genuchten, M. Th. (2006). A modified Mualem–van Genuchten formulation for improved description of thehydraulic conductivity near saturation. Vadose Zone Journal, 5, 27–
34. https://doi.org/10.2136/vzj2005.0005

RAHMATI et al.

Sharma, M. L., Gander, G. A., & Hunt, C. G. (1980). Spatial variabilityof infiltration in a watershed. Journal of Hydrology, 45, 101–122.
https://doi.org/10.1016/0022-1694(80)90008-6
Šimůnek, J., van Genuchten, M. Th., & Šejna, M. (2008). Development and applications of the HYDRUS and STANMOD software packages and related codes. Vadose Zone Journal, 7, 587–600.
https://doi.org/10.2136/vzj2007.0077
Šimůnek, J., van Genuchten, M. Th., & Šejna, 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
Smiles, D. E., & Knight, J. H. (1976). A note on the use of the Philip
infiltration equation. Soil Research, 14, 103–108. https://doi.org/10.
1071/SR9760103
Swartzendruber, D. (1987). A quasi-solution of Richards’ equation
for the downward infiltration of water into soil. Water Resources
Research, 23, 809–817. https://doi.org/10.1029/WR023i005p00809
Valiantzas, J. D. (2010). New linearized two-parameter infiltration
equation for direct determination of conductivity and sorptivity.
Journal of Hydrology, 384, 1–13. https://doi.org/10.1016/j.jhydrol.
2009.12.049
van Genuchten, M. Th. (1980). A closed-form equation for predicting
the hydraulic conductivity of unsaturated soils. Soil Science Societyof America Journal, 44, 892–898. https://doi.org/10.2136/sssaj1980.
03615995004400050002x
Vandervaere, J. P., Vauclin, M., & Elrick, D. E. (2000). Transient flowfrom tension infiltrometers I. The two-parameter equation. Soil

19 of 22

Vadose Zone Journal

Science Society of America Journal, 64, 1263–1272. https://doi.org/
10.2136/sssaj2000.6441263x
Vogel, T., & Cislerova, M. (1988). On the reliability of unsaturated
hydraulic conductivity calculated from the moisture retention
curve. Transport in Porous Media, 3, 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 flow predictions. Advances in Water Resources, 24, 133–
144. https://doi.org/10.1016/S0309-1708(00)00037-3

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 online
in the Supporting Information section at the end of thearticle.

How to cite this article: Rahmati M,
Vanderborght J, Šimůnek J, et al. Soil hydraulicproperties estimation from one-dimensional
infiltration experiments using characteristic time
concept. Vadose Zone J. 2020;19:e20068.
https://doi.org/10.1002/vzj2.20068

Appendix
The Python code applying the CTM method on infiltration data
# Description# The function finds the S and Ks from infiltration curve supplied in data using CTM procedure# The DATA is a matrix with two column and n rows providing [t, I]# Classes
import pandas as pandas # Check if pandas is installed on your system. If not, install it by: pip install pandas –userimport numpy as numpy # Check if numpy is installed on your system. If not, install it by: pip install numpy –userfrom scipy.stats import linregress # Check if scipy is installed on your system. If not, install it by: pip install scipy –userimport scipy as scipy # Check if scipy is installed on your system. If not, install it by: pip install scipy –user# Reading data from EXCEL
File_Location = ‘SimulatedData.xlsx’ # Change the direction based on your file location
Sheet_Name = ‘Clay’ # Does the calculation for Clay soil. Change it accordingly
DataLabels = [‘T’, ‘I’] # Change the labels based on your input file.
DATA = pandas.read_excel(File_Location, Sheet_Name)
HEADERs = DATA.columns
NumberOFHeaders = len(HEADERs)
for i in range(NumberOFHeaders):
DataLabels[i] = DATA[HEADERs[i]]T = DataLabels[0]I = DataLabels[1]T = T[T<=240](Continues)

20 of 22

Vadose Zone Journal

RAHMATI et al.

I = I[0:len(T)]# Setting initial values
beta = 0.6 # beta value of Haverkamp model. Change this if you want to use soil-specific value
epsilon = 0.001 # The arbitrary value between zero and 0.5 used to adjust omega. Change it to speed up the procedured = 0.001 # An arbitrary value between 0 & 1. Change it to increase the flexibility of the procedure# CTM Codedef CTM(T, I, beta):t_char, I_char, Sp, Kp1, omega = CTM_I(T, I, beta)
Kp2, t_grav, I_grav, alpha = CTM_Ks(Sp, t_char, I_char, omega, beta, T, I)
return Sp, Kp1, Kp2, t_char, I_char, t_grav, I_grav, omega, alpha# CTM Core Codedef CTM_I(T, I, beta):# Default values
t_char = numpy.nan
I_char = numpy.nan
Sp = numpy.nan
Kp = numpy.nan
omega = 0.5 + epsilon
t_char = numpy.nan# Calculation
while numpy.isnan(t_char):omega = omega - epsilon
t_char, I_char, Sp, Kp = Iterative_procedure(T, I, beta, omega)
if omega <= 0.001:
break
return t_char, I_char, Sp, Kp, omega
def Iterative_procedure(T, I, beta, omega):# Default values
t_char = numpy.nan
I_char = numpy.nan
Sp = numpy.nan
Kp = numpy.nan# Calculation
data_lenght = len(T)for k in range(1, data_lenght):
tc = T[k]Ic = I[k]SS = (1-omega)*(Ic/numpy.sqrt(tc)) # Computes a preliminary stimate for S
W = []
W = (SS*numpy.sqrt(T))/I # Computes contribution of the sorptivity component
if numpy.max(W) > = 1 - d and numpy.max(W) <= 1 + d: # Checks the selection index of the proceduret_char = tc
I_char = Ic
aa = (1/(9*(1-omega)))*(beta**2 - beta + 1) * (t_char**2)/(I_char)(Continues)

RAHMATI et al.

Vadose Zone Journal

21 of 22bb = ((2 - beta) / 3) * t_char
cc = -1 * omega * I_char
Kp = (numpy.sqrt(bb**2 - 4 * aa * cc) - bb) / (2 * aa) # Computes the Ks
Sp = (1-omega)*(I_char/numpy.sqrt(t_char))break
return t_char, I_char, Sp, Kp
def CTM_Ks(Sp, t_char, I_char, omega, beta, T, I):
W = Sp*numpy.sqrt(T)/Idel W[0]del T[0]## Use linregress if you do not want to set the intercept to zero: ln(W) = a + bT
Stat = scipy.stats.linregress(T, numpy.log(W))
alpha = Stat.slope## Use lin_fit if you want to set the intercept to zero: ln(W) = bT
#alpha = lin_fit(T, numpy.log(W))#alpha = float(alpha)if omega = = 0.5:t_grav = t_char
I_grav = I_char
else:
t_grav = numpy.log(0.5)/alphaI_grav = (Sp*numpy.sqrt(t_grav))/0.5# Calculate Ks
a1 = (2*(beta**2 - beta +1)* t_grav**2)/(9*I_grav)
b1 = ((2-beta)/3)*t_grav
c1 = -1*0.5*I_grav
Kp = (numpy.sqrt(b1**2 - 4*a1*c1)-b1)/(2*a1)
return Kp, t_grav, I_grav, alphadef lin_fit(x, y):#‘”Fits a linear fit of the form mx+b to the data”’
fitfunc = lambda params, x: params[0] * x #create fitting function of form mx
errfunc = lambda p, x, y: fitfunc(p, x) - y
init_a = -1

#create error function for least squares fit#find initial value for a (gradient)

init_p = numpy.array((init_a))#bundle initial values in initial parameters#calculate best fitting parameters (i.e. m and b) using the error function
p1, success = scipy.optimize.leastsq(errfunc, init_p.copy(), args = (x, y))#f = fitfunc(p1, x)#create a fit with those parameters

return p1# Predicting S and Ks
Sp, Kp1, Kp2, t_char, I_char, t_grav, I_grav, omega, alpha = CTM(T, I, beta)# Printing the results(Continues)

22 of 22

Vadose Zone Journal

print(‘——————————–’)
print (‘———–> Results <———-’)
print(‘——————————–’)
print(‘Predicted S =’, round(Sp, 2), ‘[L/sqrt(T)]’)
print(‘Predicted Ks by CTM-Ks’, round(Kp2,2), ‘[L/T]’)
print(‘Predicted Ks by CTM-I’, round(Kp1,2), ‘[L/T]’)
print(’t_char =‘, round(t_char,2), ’[T]’)print(’I_char =‘, round(I_char,2), ’[L],’)
print(’t_grav =‘, round(t_grav,2), ’[T]’)print(’I_grav =‘, round(I_grav,2), ’[L],’)
print(’omega =‘, round(omega,2), ’[-]’)
print(’alpha =‘, round(alpha,4), ’[1/T]’)

RAHMATI et al.