Published April 11, 2019

Special Section: Nonuniform

Flow across Vadose Zone Scales
BEST-2K Method for Characterizing
Core Ideas Dual-Permeability Unsaturated Soils
• We present a new method for the
hydraulic characterization of dual- with Ponded and Tension Iniltrometers
permeability soils. Laurent Lassabatere,* Simone Di Prima, Sofia Bouarafa,
• The BEST-2K method was validated
with synthetic soil data.
Massimo Iovino, Vincenzo Bagarello,
• The BEST-2K method was also vali- and Rafael Angulo-Jaramillo
dated with real experimental data
This study presents a new method (BEST-2K) that extends the existing BEST
for contrasting soils.
methods for use in characterizing the water retention and hydraulic conductiv-
ity functions of matrix and fast-flow regions in dual-permeability soils. BEST-2K
requires input information from two water infiltration experiments that are
performed under ponded (Beerkan) and unsaturated (tension infiltrometer) con-
ditions at the surface. Other required inputs include water content measurements
and the traditional BEST inputs (particle size distribution and bulk density). In this
study, first, a flowchart of the BEST-2K method was developed and illustrated
with analytically generated data for a synthetic dual-permeability soil. Next, a
sensitivity analysis was performed to assess the accuracy of BEST-2K and its sensi-
tivity to the quality of the inputs (water contents and cumulative infiltrations, and
the prior estimation of the volume ratio occupied by the fast-flow region). Lastly,
BEST-2K was applied to real experimental data to characterize three soils that are
prone to preferential flow. BEST-2K was found to be a particularly useful tool that
combines experimental and modeling approaches for characterizing dual-per-
meability soils and, more generally, soils prone to preferential flows.
Abbreviations: DP, dual-permeability; PSD, particle size distribution; PTF, pedotransfer function; SP, sin-
gle-permeability.

The hydraulic characterization of soils, particularly their unsaturated properties,

is a prerequisite for understanding water flow in the vadose zone. Several water infiltration
techniques have been developed to identify the soil hydraulic properties (Angulo-Jaramillo
L. Lassabatere, S. Di Prima, S. Bouarafa, and R. et al., 2016). Water infiltration techniques have several clear advantages: they are non-
Angulo-Jaramillo, Univ. Lyon, Univ. Claude Ber-
nard Lyon 1, CNRS, ENTPE, UMR5023 LEHNA, intrusive, do not involve sampling or disturbance of the connection to the surrounding
F-69518, Vaulx-en-Velin, France; M. Iovino and soil environment, and can be performed relatively quickly and at reasonable costs. Among
V. Bagarello, Dep. of Agricultural, Food and
Forest Sciences, Univ. of Palermo, Viale delle these methods, the BEST methods were developed to determine the whole set of hydraulic
Scienze, 90128 Palermo, Italy. *Correspon- parameters related to the hydraulic functions, namely the water retention function that
ding author ([email protected]).
describes the soil’s ability to retain water by capillarity and the unsaturated hydraulic
conductivity that describes the soil’s ability to allow water flow (Lassabatere et al., 2006;
Received 27 June 2018. Yilmaz et al., 2010; Bagarello et al., 2014). These methods use raw soil data (including water
Accepted 7 Dec. 2018.
Supplemental material online. contents, bulk density, and particle size distribution) and one Beerkan test (consisting of
a water infiltration experiment under ponded conditions through a ring). However, these
methods are suitable only for single-permeability soils, whereas many soils have been shown
Citation: Lassabatere, L., S. Di Prima, S. Bouarafa,
M. Iovino, V. Bagarello, and R. Angulo-Jaramillo. to have several pore structures and thus should be described using multimodal or at least
2019. BEST-2K method for characterizing dual- dual-permeability approaches (Durner, 1994).
permeability unsaturated soils with ponded
and tension infiltrometers. Vadose Zone J. In this study, we introduce the BEST-2K method for the hydraulic characterization
18:180124. doi:10.2136/vzj2018.06.0124 of dual-permeability (DP) soils. BEST-2K is based on the principle of the BEST method,
which was developed for the hydraulic characterization of single permeability (SP) soils
(Lassabatere et al., 2006). The objectives of this study were to introduce the BEST-2K
framework and test it with both numerical and real experimental data. The accuracy and
© 2019 The Author(s). This is an open access article
distributed under the CC BY-NC-ND license robustness of BEST-2K was evaluated using analytically generated (error-free) data for a
(http://creativecommons.org/licenses/by-nc- synthetic dual-permeability soil. BEST-2K was then applied to real data for several types
nd/4.0/).

Vadose Zone Journal | Advancing Critical Zone Science

of soils and land management conditions, including pasture, forest,
and orchards. where q 2K is the bulk volumetric water content of the DP soil, q f
and q m are the local volumetric water contents in the fast-flow and

6 BEST-2K Method matrix regions, respectively, K2K is the bulk hydraulic conductivity,
Kf and Km are the local hydraulic conductivities of the fast-flow
BEST-2K Water Retention and matrix regions, respectively, and wf is the volume fraction
and Hydraulic Conductivity Functions of the DP soil occupied by the fast-flow region. By analogy with
Dual-permeability soils are composed of two regions: a matrix BEST method (Lassabatere et al., 2006), both regions are assigned
flow region that hosts the smaller pores and a fast-flow region that the van Genuchten (1980) model with the Burdine condition for
hosts the larger pores. These regions are described as Darcean the description of their water retention functions and the Brooks
porous media with different local water retention and hydraulic and Corey (1964) model for the description of their unsaturated
conductivity functions (Fig. 1a and 1b). The bulk water content hydraulic conductivity functions:
and hydraulic conductivity correspond to the summation of the é æ nm ù -mm
contributions of the matrix and fast-flow regions, as suggested by ê ç h ö÷÷ ú
q m (h )= (q s,m -q r,m )ê1 + çç ÷ ú +q r,m [2a]
Gerke and van Genuchten (1993) (Fig. 1c): ê çè hg,m ø÷÷ ú
ëê ûú
q 2 K (h )= wf q f (h )+ (1- wf )q m (h ) [1a]
é æ ö÷nf ù -mf
ê ç h ÷÷ ú
K 2 K (q )= wf K f (q f )+ (1- wf )K m (q m ) [1b] q f (h )= (q s,f -q r,f )ê1 + çç ÷ ú +q r,f [2b]
ê çè hg,f ø÷ ú
êë úû

Fig. 1. Framework of (a) the dual-permeability (DP) approach (adapted from Gerke and van Genuchten, 1993), (b) regions and water contents in the
DP soil, (c) 2K hydraulic functions of the DP soil, and (d) BEST-2K water infiltration experiments.

VZJ | Advancing Critical Zone Science p. 2 of 20

æ q m - q r,m ÷öhm threshold corresponds to the deactivation of pores with a radius
K m (q m )= K s,m ççç ÷÷ [2c] >0.5 mm (Timlin et al., 1994). Practically, the value of the water
èç q s,m - q r,m ÷ø pressure head could be adapted to the expected size of macropores
hf constituting the fast-flow region.
æ q - q r,f ÷÷ö
K f (q f )= K s,f ççç f ÷÷ [2d] Figure 2 illustrates the framework of the BEST-2K method.
çè q s,f - q r,f ø The preprocessing functions (Fig. 2, A) use the aforementioned
experiments to determine the inputs, which would have been
2 obtained with the regular Beerkan method had each region been
mm = 1- [2e]
nm sampled separately. Next, the regular BEST method is applied to
the computed inputs to determine the hydraulic parameters for
2
mf = 1- [2f ] each region (Fig. 2, B). Note that the regular BEST methods used
nf in BEST-2K are referred to as BEST-1K, in which 1K stands for
single permeability, in analogy with BEST-2K in which 2K stands
2
hm = +2 + p [2g] for dual-permeability soils. Finally, the hydraulic parameters esti-
nm mm
mated for the matrix and fast-flow regions are used to build the
water retention and hydraulic conductivity functions of the DP
2
hf = +2 + p [2h] soils (Fig. 2, C).
nf mf
The following is a more detailed description of the BEST-2K
where m, n, and h are shape parameters, q r and q s are the residual method. The preprocessing functions divide the bulk PSD
and saturated water contents, respectively, hg is the scale parameter into two distributions, PSDm and PSD f, which describe the
for water pressure head, Ks is the saturated hydraulic conductiv- composition of the matrix and fast-flow regions, respectively.
ity, p is the tortuosity parameter fixed at 1.0 (Lassabatere et al., These functions use the measured water contents (q s,2K, q TI, q 0,B,
2006), and subscripts m and f refer to the matrix and the fast-flow and q 0,TI) to derive the saturated and initial water contents for
regions, respectively. Note that the subscript m refers to the matrix the matrix and fast-flow regions (q s,m, q 0,m, q s,f, and q 0,f ). The
component of the soil, whereas the variable m refers to the shape preprocessing functions also use ITI to compute the cumulative
parameter of the van Genuchten model, m = 1 − 2/n. BEST-2K infiltration (I B,m) that would have infiltrated into the matrix
estimates the whole set of the local hydraulic parameters for the alone under zero water pressure head (in agreement with the
matrix region (q r,m, q s,m, hg,m, Ks,m, nm, and h m) and the fast-flow Beerkan method). In this step, the inputs PSDm , q 0,m , q s,m ,
region (q r,f, q s,f, hg,f, Ks,f, nf, and h f ). and I B,m for the matrix and PSD f, q 0,f, and q s,f for the fast-
flow region are determined (Fig. 2, arrows d1). The set of data
BEST-2K Framework previously computed for the matrix (PSDm, q 0,m, q s,m, and IB,m)
The main concept of BEST-2K involves the combination of is processed by the regular BEST method, BEST-1K, to derive
water infiltration experiments performed at two water pressure the hydraulic parameters for the matrix (q r,m, q s,m, hg,m, Ks,m,
heads: one slightly positive and the other slightly negative. These nm , and h m) (Fig. 2, d 2). Then this set of parameters is used
experiments quantify water infiltration when the whole pore net- to compute the cumulative infiltration into the matrix during
work (both matrix and fast-flow regions) is activated and when the Beerkan experiment, I B,m(tB), and thus its contribution to
only the matrix is activated (Fig. 1d). The contribution of the fast- the total cumulative infiltration, I B , to derive the cumulative
flow region to water infiltration is determined by comparing the infiltration into the fast-flow region, I B,f (Fig. 2, d 3). The I B,f
resulting cumulative infiltrations obtained from both experiments. corresponds to the cumulative infiltration that would have
The main inputs for BEST-2K include two water infiltration infiltrated into the fast-flow region on its own, under zero water
experiments (Fig. 1d): a Beerkan experiment to quantify the bulk pressure head (Beerkan method, Braud et al., 2005). Next, the
cumulative water infiltration under ponded conditions, IB(tB), set of data computed for the fast-flow region (PSDf, q 0,f, q s,f, and
and a tension infiltrometer (TI) experiment performed at a suc- I B,f ) is processed by the BEST-1K method to derive the set of
tion high enough to quantify the cumulative infiltration into the hydraulic parameters for the fast-flow region (q r,f, q s,f, hg,f, Ks,f,
matrix alone, ITI(t TI). The other inputs include the initial water nf, and h f ) (Fig. 2, d4). Thus, the full set of hydraulic parameters
contents for the two tests, q 0,B and q 0,TI, the final water content for the DP soil are obtained, and the bulk water retention and
at the end of the tension infiltrometer experiment, q TI, the soil hydraulic conductivity functions can be determined.
particle size distribution (PSD), and the bulk density, r d, which
is used to derive the bulk saturated water content, q s,2K. Based on BEST-2K Preprocessing Functions
previous studies (Timlin et al., 1994), we determined that a suc- The novelty of BEST-2K is mostly in its use of the
tion of 30 mm (i.e., a water pressure head of −30 mm) is sufficient preprocessing functions, which provide the sets of inputs PSDm,
to deactivate the fast-flow region, thus restricting water flow to q 0,m, q s,m, and IB,m and PSDf, q 0,f, q s,f, and IB,f that are processed
the matrix region, as required for the TI experiment. In fact, this separately by BEST-1K. These functions are based on several

VZJ | Advancing Critical Zone Science p. 3 of 20

Fig. 2. Flowchart of BEST-2K method using the Beerkan and the tension infiltrometer (TI) methods. The letters A, B, C, and di refer to comments in
the text; the uppercase letters indicate the main steps, and di indicate the detailed steps.

assumptions that are described and discussed below. BEST-1K is the bimodality of the water retention curve, the pore size distribu-
also described below and summarized in Table 1. tion, and the PSD is assumed, and the following function is used
to represent the bimodal PSDs:
Particle Size Distributions é æ D öN f ù - M f
ê ÷÷ úú
g,f ÷
Based on the description of soils as fragmented fractal FF2 K (D )= t f ê1 +ççç
porous media, the water retention function, pore size, and PSD ê è D ø÷ ú
ç
ë û [3]
are expected to be consistent with one another (Rieu and Sposito,
é æ D öN m ù - M m
1991). The bimodal water retention functions represent the pro- ê g,m ÷ ú
+ (1-t f )ê1 +ççç ÷÷ ú
gressive activation of the two modes of the pore size distribution. ê èç D ø÷ ú
ë û
The matrix region, which is composed of small particles with small
pores, may activate readily due to higher capillary forces. When where FF2K (D) is the cumulative distribution of particle size, D
the water pressure head increases enough, the larger pores of the is the particle diameter, tf is the fractional contribution of the
fast-flow region may fill up and raise the water content closer to particles in the fast-flow region (unknown a priori), Dg is the
saturation. These larger pores are expected to surround the larger average diameter of either particle size mode, and N and M are tex-
particles. Relationships between the water retention functions and tural parameters related by the expression M = 1 − 2/N, which is
pore size (or particle size) distributions have already been suggested analogous to the Burdine condition. Calibration to the bulk PSD
(e.g., Arya and Paris, 1981). Consequently, a strong link between provides optimized values for N, M, and Dg , which determines

VZJ | Advancing Critical Zone Science p. 4 of 20

Table 1. Synthesis of BEST methods. the particle size distributions for the matrix and fast-flow regions,
PSDm and PSDf, respectively. The textural parameters (Nm, Mm)
BEST water retention and hydraulic conductivity functions
and (Nf, Mf ), which define PSDm and PSDf, respectively, are
q(h) = (q s − q r)[1 + (h/hg)n]−m + q r [t1a]
inputs to the BEST-1K pedotransfer functions (PTFs) that pro-
K(q) = Ks[(q − q r)/(q s − q r)] h [t1b] vide the shape parameters for the hydraulic functions of the matrix
qr = 0 [t1c] region (nm, h m) and the fast-flow region (nf, h f ).
q s = 1 − r d/r s [t1d]
m = 1 − 2/n [t1e] Initial and Saturated Water Contents
[t1f]
The initial water content is considered negligible in the fast-
h = 2/nm + 2 + p with p = 1
flow region, given the very little water retention by capillarity in
BEST pedotransfer functions for particle size distribution the large pores; thus, the water remains in the matrix, occupying
Fitting the particle size distribution: FF(D) = [1 + (Dg /D)N]−M 1 − wf of the total soil volume, leading to
[t2]
with M = 1 − 2/N
s root of (1 − e)s + e 2s = 1 with e = q s
q 0,TI
[t3] q 0,m = [4a]
1- wf
pm = [MN/(1 + M)](1 + k)−1 [t4]

with k = (2s − 1)/[2s(1 − s)] [t5] q 0,f = 0 [4b]

n = 2/{1 − (1/pm)[Ö(1 + pm 2) − 1]} [t6]
The bulk saturated water content, q s,2K, is equal to the soil
h = 2/mn + 3 [t7]
porosity, which was previously derived from the soil bulk and min-
cp = G(1 + 1/n){[G(mh − 1/n)/G(mh)] [t8] eral densities, r d and r s, respectively), as proposed by Lassabatere
+ [G(mh + m − 1/n)/G(mh + m)]}
et al. (2006) for the BEST methods:
BEST analytical models
rd
Transient state q s,2 K = 1- [5]
rs
IO(2)(t) = SÖt + (AS 2 + BKs)t [t9a]
qO(2)(t) = S/(2Öt) + (AS 2 + BKs) [t9b] The mineral density is often assumed to be 2.65 g cm−3 when
Steady state
no specific measurement can be obtained. To obtain the local satu-
rated water contents, the two volumetric water contents, q TI and
I+¥(t) = (AS 2 + Ks)t + C(S2/Ks) [t10a]
q s,2K, are expressed as functions of the volumetric water contents
q+¥ (t) = AS 2 + K s [t10b] in the matrix and the fast-flow regions, leading to
A = g/[rd(q s − q 0)] [t11]
q TI = wf q TI,f + (1- wf )q TI,m [6a]
B = (1/3)(2 − b)[1 − (q 0/q s) h] + (q 0/q s) h [t12]
C = 1/{2[1 − (q 0/q s) h](1 − b)}ln(1/b) [t13] q s,2 K = wf q s,f + (1- wf )q s,m [6b]
BEST fitting algorithm for estimating sorptivity, S
where q TI,m and q TI,f denote the local water content in the
Regression of the last data points for steady state:
[t14] matrix and the fast-flow regions, respectively, at the end of the
(texp, Iexp)i Î {ntot − 3, …, ntot} ® slope qsexp & intercept bsexp
TI experiments. Because the tension infiltrometer deactivates the
BEST-Steady
fast-flow region, we assume that no water remains in this region,
S best = Ö[qsexp/(A + C/bsexp)] [t15]
i.e., q TI,f = 0. We also assume that the water pressure head is suf-
BEST-Slope & BEST-Intercept fitting procedure:
ficient to activate all the pores in the matrix region for the Beerkan
BEST-Slope model
and tension infiltrometer experiments; thus, q TI,m = q s,m. These
IO(2)(t,S) = SÖt + [A(1 − B)S 2 + Bqsexp)]t [t16]
hypotheses define a set of four equations with four unknowns,
BEST-Intercept model
leading to the following equations for the saturated water contents
IO(2)(t,S) = SÖt + [AS 2 + BC(S 2/bsexp)]t [t17]
in the matrix and fast-flow regions:
k
Sopt(k) minimizes OF(S) = S i=1 {Iexp(i) − IO(2)[texp(i),S]}2 [t18a]
q TI
S best = Sopt{max[k|texp(k) ≤ tmax(k)]} q s,m = [7a]
with tmax(k) = 1/[4(1 − B)2][Sopt(k)/Ks,opt(k)]2
[t18b] 1 - wf

BEST estimation for saturated hydraulic conductivity, Ks q s,2 K -q TI

exp − AS 2
q s,f = [7b]
Ks,best = qs best BEST-Slope [t19a] wf
Ks,best = CS best2/bs exp BEST-Intercept [t19b]
Given these equations, the saturated water contents are
Ks,best = Cqsexp/(Absexp + C) BEST-Steady [t19c]
related to the volume fraction occupied by the fast-flow region,
BEST estimation for scale parameter for water pressure head, hg
w f ; this parameter is estimated a priori. However, w f must be
hg,best = S best2/{cpKs,best(q s − q 0)[1 − (q 0/q s) h]} [t20] chosen such that Eq. [7a–7b] predict values between zero and

VZJ | Advancing Critical Zone Science p. 5 of 20

one. Indeed, by definition, the water content is less than one For the TI experiment, it is assumed that water flows exclu-
because it is the ratio of the volume of water to the total volume sively in the matrix and that infiltration into the fast-flow region
of water, air, and solid. In addition, the saturated water content is negligible, i.e., ITI,f = 0. It is also assumed that water infiltra-
in the matrix region is often less than that in the fast-flow region tion does not change appreciably between the Beerkan and TI
(e.g., Gerke and van Genuchten, 1993), which leads to the fol- experiments for the matrix, i.e., ITI,m = IB,m. The combination
lowing inequality: of these hypotheses with Eq. [11a] allows a direct computation of
the cumulative infiltration into the matrix for the Beerkan experi-
q s,m £ q s,f <1 [8]
ments, IB,m (see Eq. [12a]). Then this infiltration is subtracted
This condition, along with Eq. [7a–7b], implies the following from the total cumulative infiltration to obtain the cumulative
conditions on wf: infiltration into the fast-flow region, IB,f, as
q s,2 K -q TI I TI
q s,2 K -q TI £ wf £ [9] I B,m = [12a]
q s,2 K 1- wf
The strategy outlined above defines the BEST-2K-A method. I B - I TI
However, it relies on a very accurate measure of the difference I B,f = [12b]
wf
q s,2K − q TI, which may be complicated given the uncertainty of
water content measurements in the field. Consequently, a second I B -(1 - wf )I B,m (t B )
method, BEST-2K-B, was defined based on the assumption that I B,f = [12c]
wf
the saturated water contents are the same in the two regions and
equal the total saturated water content, q s,2K: Very often, the cumulative infiltrations, ITI(t TI) and IB(tB),
are not obtained at the same times, t TI ¹ tB . Thus, ITI cannot
q s,f = q s,2K [10a]
be directly subtracted as shown in Eq. [12b]. To circumvent this
inconsistency, Eq. [12c] is used to compute I TI(tB) = (1 − wf )
q s,m = q s,2K [10b]
IB,m(tB) for the time dataset tB. The infiltration IB,m(tB) is directly
In this case, there is no constraint on the value of parameter computed for the times tB using the hydraulic parameters esti-
wf. Indeed, the saturated water contents are equal to the bulk satu- mated for the matrix (q r,m, q s,m, hg,m, Ks,m, mm, nm, h m) and the
rated water content, q s,2K, and do not depend on wf. Thus, there quasi-exact model proposed by Haverkamp et al. (1994):
are no risks of errors in the saturated water contents, unlike the g S m2
case of BEST-2K-A, where errors can occur when Eq. [7a–7b] are I B,m (t B )= I1D,m (t B )+ tB [13a]
rd Dq m
calculated using erroneous values of wf.
From a physical point of view, BEST-2K-B considers that the 2D K m2 1 æç 2D K m é
tB = ç I1D,m (t B )- K 0,m t B ùû
void ratio, i.e., the ratio of the volume of the voids to the volume S m2 1-b çèç S m2 ë
of the solids, does not evolve with the particle size, meaning that
the porosity remains the same in both the matrix (fine particles) - ln íï
ï { ( m m )ë 1D,m (t B )- K 0,m t B ùû +b -1 üïïïö÷÷÷
ìï exp b 2D K S 2 é I
} [13b]
ý÷÷
ïïï b ïï÷ ÷
and the fast-flow region (large particles). In contrast, BEST-2K-A î þï÷
ïø
considers that the void ratio (and thus, the porosity) may increase
with the size of the particles. The fast-flow regions could also be where tB is time for the Beerkan experiment, rd is the radius of the
viewed as an ensemble of capillary tubes of large dimensions with disk source, Dq m (= q s,m − q 0,m) is the difference between the
a porosity equal to one and the lack of any surrounding particles. final and initial water contents, DK [= Ks,m − K0,m, with K0,m =
However, such a concept is beyond this study that considers dual- K(q 0,m)] is the difference between the final and initial hydraulic
permeability systems with a combination of two porous media conductivities, and b and g are constants commonly set at 0.6 and
made of voids and particles. 0.75, respectively (Haverkamp et al., 1994; Smettem et al., 1994).
The sorptivity, Sm, can be computed from the water diffusivity, Dm,
Cumulative Iniltration into the Matrix and the initial and final water contents using Parlange’s approxi-
and Fast-Flow Regions mation (Parlange, 1975):
The cumulative infiltrations into the DP soil for both TI q s,m
and Beerkan experiments, I TI and I B , respectively, correspond S m2 (q 0,m , q s,m )= ò (q s,m + q - 2q 0,m ) Dm (q )d q [14a]
q 0,m
to the summation of the infiltrations into the matrix and fast-
flow regions: dh
Dm (q m )= K m (q m ) m [14b]
I TI = wf I TI,f + (1- wf ) I TI,m [11a] dqm

where q is a dummy variable. The diffusivity Dm(q m) can be

I B = wf I B,f + (1- wf ) I B,m [11b]
easily computed from the hydraulic conductivity function

VZJ | Advancing Critical Zone Science p. 6 of 20

K m(q m), and the derivative of the function h m(q m) with respect infiltration into each region multiplied by the volume fraction
to q m, written as dhm/dq m. The function h m(q m) corresponds occupied by each region, and each cumulative infiltration was
to the inverse function of the water retention function q m(h) computed using the quasi-exact implicit model proposed by
defined by Eq. [2a]. More details on the use and validation of Haverkamp et al. (1994):
the model proposed by Haverkamp et al. (1994) can be found in é g S2 ù
Lassabatere et al. (2009). I 3D,2 K (t )= wf êê I1D,f (t )+ f f t úú
ëê rd Dq f ûú

6 Materials and Methods

é ù [16a]
g S2
+ (1- wf )êê I1D,m (t )+ m m t úú
ëê rd D q m ûú
Synthetic Data for Analytical Validation of BEST-2K
The synthetic DP soil is a loamy matrix with macropores, 2D K f2 1 é 2D K
ê f é ù
which corresponds to the loam texture defined in the soil catalog S f2
t=
1 - bfê S 2 ë I1D,f (t )- K 0,f t û
êë f
of Carsel and Parrish (1988). The hydraulic functions of this soil æ ù
ç ì
ï 2D K f ü
ï ÷ö
were described using the same equations as in BEST, i.e., the van ï
ççç exp ï
íbf ë û ý +b f - 1 ÷÷÷úú
é I1D,f (t )- K 0,f t ù ï [16b]
ï
2
Sf ï
ï ÷÷ú
Genuchten (1980) relationship with the Burdine (1953) condition - ln ççç î þ ÷÷ú
çç bf ÷÷ú
for the water retention function and the Brooks and Corey (1964) çç ÷÷ú
çè ø÷úû
relationship for the hydraulic conductivity. The fast-flow region
corresponds to an ensemble of macropores of 1-mm average radius
2D K m2 1 é 2D K
and occupies 10% of the bulk soil (Vtot,f/Vtot = 0.1; Fig. 1). The ê m é ù
S m2
t=
1 - bm ê S 2 ë I1D,m (t )- K 0,m t û
scale parameter for water pressure head, hg,f, was derived from the êë m
æ ì 2D K
ï ü öù
pore radius, rg,f, using the Young–Laplace equation (Kutilek and ç ïb m é ïï +b - 1 ÷÷ú
I1D,m (t )- K 0,m t ùû ý
ççç exp í ë ÷ ú [16c]
ï ï ÷÷ú
m 2 m
Nielsen, 1994; Lenhard et al., 2005): - ln ççç ï
î Sm ï
þ ÷÷÷ú
ççç bm ÷÷÷úú
z çç ÷
hg,f =- [15a] è ø÷ûú
rg,f

2saw cos (b c ) where the subscripts m and f denote the parameters of the matrix
z= [15b] and the fast-flow regions, respectively, and the variables rd, Dq,
rw g
DK, b, and g are defined as in Eq. [13]. The sorptivity, S, of each
where s aw is the surface tension of the air–water interface, b c region was computed with Eq. [14] using the initial and final water
is the contact angle, r w is the water density, and g is the gravita- contents for each region.
tional acceleration constant, leading to a value of 14.9 mm 2 for We simulated a Beerkan experiment (zero water pressure
z for the case of pure water. The value of the hydraulic conduc- head at the surface) and a TI experiment with a water pressure
tivity, Ks,f, was computed from that of the loamy matrix, Ks,m, head fixed at −30 mm at the surface, which is twice the value
assuming a linear increase with the square of the pore radius, of the scale parameter for water pressure head for the fast-flow
as indicated by Poiseuille’s law (Sutera and Skalak, 1993) and region (hg,f = −14.9 mm) and should be enough to deactivate
suggested by several other studies (e.g., Watson and Luxmoore, the fast-flow region during the TI experiment. Besides, this
1986; Timlin et al., 1994). The residual water content, q r,f, was threshold is usually considered as a guide value for the hydraulic
set at zero and the saturated water content, q s,f, at a large value characterization of water infiltration with infiltrometers (Timlin
of 0.70. The shape parameter, nf, was set at 3.75 to induce a steep et al., 1994). The water content at the end of the Beerkan experi-
shape for the water retention functions, as commonly used for ment was equal to the bulk saturated water content, q s,2K. For the
coarse soils (Schaap et al., 2001). Note that the matrix region TI experiment, the initial and final water contents, q 0,TI and q TI,
and the bulk soil have ordinary porosities (approximately 43 and were computed from the local water retention functions using
45.7%, respectively), whereas the fast flow-region was assigned Eq. [1a], considering water pressure heads of −10 m and −30 mm,
a high porosity (70%), assuming an ensemble of macropores respectively. For the Beerkan experiment, the same initial water
surrounded by tiny walls made of few particles. In total, the content was used. The cumulative infiltrations were computed
fast-flow region occupies 10% of the bulk DP soil and its poros- for ideal conditions with a precise description of the transient
ity constitutes 15.3% of the bulk porosity (i.e., w f q s,f/q s,2K ). state and attainment of the steady state. For the TI experiment,
The studied synthetic soil was designed to exhibit a typical DP the cumulative infiltrations were computed in increments of
behavior (see below). 5 mm with a total cumulative infiltration of about 675 mm, and
Water infiltrations were analytically modeled using the ana- for the Beerkan experiment, the incremental and total cumula-
lytical model developed by Lassabatere et al. (2014), where the tive infiltrations were about 0.15 and 40 mm, respectively. The
three-dimensional, axisymmetric cumulative infiltration into total durations were 1500 and 10 min for the TI and Beerkan
DP soils was computed by the summation of the cumulative experiments, respectively.

VZJ | Advancing Critical Zone Science p. 7 of 20

To test the robustness of BEST-2K methods with respect to Sardinia site (agroforestry management; elevation of 333 m asl,
erroneous inputs, we simulated several scenarios. We inverted the near Villanova Monteleone in the province of Sassari) had a loam
analytical data using BEST-2K-A and BEST-2K-B and fixed the soil. This forest site had a low tree density (about 180 trees per
parameter wf at several values other than the nominal value of 10%, hectare); the dominant tree species was the evergreen Quercus
specifically, 2, 4, 6, 7, 8, 9, 20, 30, and 40%. For the water contents, suber L. (cork oak), and forage species such as Avena, Trifolium,
q TI was varied from 0.1 to 0.45 in steps of 0.01 for a nominal value and Lolium grew under the trees. The Palermo University site
of 0.405. The value of q s,2K was varied from 0.41 to 0.7 in steps (orchard; elevation of 38 m asl) had a sandy loam soil, and the
of 0.01, for a nominal value of 0.457. For all these scenarios, the trees were spaced in a 4- by 4-m grid.
estimated hydraulic functions were compared with the estimates For a given site, undisturbed soil cores (0.05 m in height
obtained with the nominal input values. and 0.05 m in diameter) were collected at 0- to 0.05- and 0.05-
In the estimation of the shape parameters, we considered to 0.10-m depths and at three randomly selected points. Three
that the PTFs applied to the bulk PSD lead to the target values disturbed soil samples (at 0–0.10-m depth) were also collected.
for shape parameters nm and nf. The bulk and local PSDs can be The undisturbed soil cores were used to determine the dry soil
computed from the shape parameters Nm and Nf that correspond bulk density, r d, and the initial soil water content, q 0, in the labo-
to the values of nm and nf with respect to the PTFs Eq. [t3–t6] in ratory; the bulk saturated water content, q s,2K, was determined
Table 1. In this analytical study, we did not question the accuracy from r d (Eq. [5]). The disturbed soil was used to determine the
of shape parameter estimates of BEST-2K and instead focused on PSD using conventional methods after H2O2 pretreatment to
the quality of the scale parameter estimates (q r,m, q s,m, hg,m, Ks,m eliminate organic matter and clay deflocculation using sodium
and q r,f, q s,f, hg,f, Ks,f ), assuming that the shape parameters (nm, h m metaphosphate and mechanical agitation (Gee and Bauder, 1986).
and nf, h f ) are perfectly estimated. This aspect will be the subject In particular, the fine size fractions were determined using the
of further studies. hydrometer method, and the coarse fractions were obtained
by mechanical dry sieving. According to the USDA standards,
Experimental Data for the Validation of BEST-2K three fractions, namely clay (0–2 mm), silt (2–50 mm), and sand
Three sites with different land uses were sampled for this (50–2000 mm), were also determined.
investigation (Fig. 3). Two of these sites, located in northwestern The cumulative infiltrations, IB(t) and ITI(t), as functions of
Sardinia (Italy), were sampled in November 2017. The third site, time, t, were determined at zero water pressure head and a suction
located in a citrus orchard at the Department of Agriculture, of 30 mm using the Beerkan method and an SW-080B infiltrom-
Food and Forest Sciences of Palermo University (Italy), was eter. At each site, two Beerkan infiltration experiments were
sampled in July 2016. The soils are described according to the performed using a ring with an inner diameter of 0.15 m inserted
USDA soil classification system. The first Sardinia site (pas- to a depth of about 0.01 m to avoid lateral loss of the ponded water,
ture; elevation of 20 m asl, near Alghero, about 6 km from the as recommended by Lassabatere et al. (2006). A known volume
Mediterranean Sea) had a sandy clay loam soil. The second of water (150 mL) was poured into the cylinder, and the elapsed

Fig. 3. Experimental sites: (a) locations of the study areas in Italy; and (b,c) pasture, (d,e) forest, and (f,g) orchard sites.

VZJ | Advancing Critical Zone Science p. 8 of 20

time during the infiltration was measured. When the water had the steady-state straight line described by the final points (Table
completely infiltrated, the same amount of water was poured into 1, Eq. [t14]). Finally, BEST Steady makes use of only the steady-
the cylinder again, and the time needed for the water to completely state straight line and its slope and intercept (Table 1, Eq. [t14]).
infiltrate was measured once more. This process was repeated until The saturated soil hydraulic conductivity, Ks, and sorptivity, S, are
the difference in infiltration times between three consecutive water then derived from these fitting procedures (Table 1, Eq. [t15–t19]).
supplies was negligible; this indicated that a practically steady infil- The scale parameter for the water retention function, hg , is derived
tration rate had been achieved. An SW-080B infiltrometer with a from these previous estimates (Table 1, Eq. [t20]). Thus, all unsat-
0.20-m-diameter porous plate was applied close to each Beerkan urated hydraulic parameters are estimated, and the complete water
experiment. The intermediate volumetric water content, q TI, was retention and hydraulic conductivity functions can be determined.
determined on a small amount of soil sampled under the porous In comparison to the application of BEST methods for SP
plate at the end of the water infiltration experiment. soils in the field, the use of BEST-1K methods within the frame-
Each site was characterized by single r d, q 0, q s,2K, q TI, and work of BEST-2K differs in the quantification of the inputs. The
PSD data obtained by averaging the individual measurements, bulk saturated water content, q s,2K, is still defined using the bulk
except for the orchard site where specific measurements were per- density. However, the local saturated water contents q s,m and q s,f
formed similar in manner to the TI and Beerkan experiments at are computed from q s,2K and qTI using the BEST-2K preprocessing
two different locations. The experimental infiltration data were functions. In the BEST-2K framework, BEST-1K no longer fits the
inverted using the BEST-2K-A method along with three meth- bulk PSD to the unimodal model of Eq. [t2] (Table 1). Instead, the
ods: Slope, Intercept, and Steady for BEST-1K. For the BEST-2K-A bulk PSD is fitted to a bimodal model (Eq. [3]) to derive the PSD
method, without additional information, wf was estimated as the shape parameters Nm, Mm, Nf, and Mf to be used in the BEST-1K
arithmetic mean of the constraints defined by Eq. [9]. PTFs (Eq. [t3–t7]) for the derivation of the shape parameters nm,
h m, nf, and h f. Last, the initial water contents (q 0,m, q 0,f ) and the
BEST-1K Methods Beerkan cumulative infiltration (IB,m, IB,f ), which are needed for
The BEST methods are synthesized in Table 1, using equa- the application of the BEST-1K fitting functions, are provided by
tions that begin with the letter t. As described above, these BEST-2K preprocessing functions from the bulk measurements
methods use the van Genuchten (1980) relation with Burdine con- (q 0, q TI, ITI, IB). More details are provided with the illustration of
dition for the water retention function, h(q), and the Brooks and the treatment of one dataset in the supplemental material.
Corey (1964) model for hydraulic conductivity, K(q) (Table 1, Eq.
[t1]). The assumptions of the BEST methods are described next.
The value of q r is assumed to be zero, and the saturated water con-
6 Results
tent q s equals the soil porosity, which is computed by Eq. [t1d]. The Validation of BEST-2K with Synthetic Data
shape parameters n and h are estimated from the PTFs that use the Synthetic Experimental Data
PSD of the size fraction <2 mm (Table 1). The PSD is fitted to Eq. Before presenting the application of BEST-2K, we present the
[t2] (Table 1), which is consistent with the water retention function analytical data (hydraulic curves, water contents, and cumulative
(Table 1, Eq. [t1a]), and the related shape parameters, M and N, are infiltrations) for the synthetic soil. The synthetic soil exhibits the
inserted in the PTFs (Table 1, Eq. [t3–t7]). Then the scale param- water retention and hydraulic conductivity functions that show
eters, Ks and hg , are derived from the analysis of the cumulative the effect of the fast-flow region (Fig. 4a and 4b). This region
infiltrations. Three BEST methods were developed, including the induces a steep increase in water content and hydraulic conduc-
original method (BEST Slope; Lassabatere et al., 2006), a method tivity close to saturation; in particular, the hydraulic conductivity
dedicated to coarse media (BEST Intercept; Yilmaz et al., 2010), increases by two orders of magnitude (Fig. 4b). This increase
and another relying only on steady state (Bagarello et al., 2014). draws an inflection point on the plots of the water retention and
These three methods use the same analytical model to fit the hydraulic conductivity functions (Fig. 4a and 4b, A) that marks
experimental data, which includes two approximate expansions the transition between the activation of the matrix alone and the
for the transient and steady states (BEST analytical models, Table concomitant activations of the matrix and fast-flow regions. This
1). These expansions were developed to approximate the quasi- graphical pattern is typical in DP systems (Durner, 1994) and must
exact implicit formulation proposed by Haverkamp et al. (1994) be considered as one of the crucial points for assessing the accuracy
and were validated numerically by Lassabatere et al. (2009). The of BEST-2K.
three methods differ in how they perform the fit (BEST fitting The occurrences of flow in both the matrix and fast-flow
algorithm for cumulative infiltration, Table 1). BEST Slope uses regions impact the soil response to hydraulic solicitations. The
the first part of the cumulative infiltration curve to fit the tran- water content measured at the end of the TI experiment is less
sient state expansion (Table 1, Eq. [t16]) and the slope of the final than the bulk saturated content due to the deactivation of the
part to describe steady state (Table 1, Eq. [t14]). BEST Intercept fast-flow region, which leads to q TI < q s,2K (Fig. 4c, DP). Instead,
uses the first part of the cumulative infiltration to fit to the same the water content would have remained constant for SP soils com-
transient expansion (Table 1, Eq. [t17]) but uses the intercept of posed of only the matrix, i.e., q TI » q s,m (Fig. 4c, SP). However,

VZJ | Advancing Critical Zone Science p. 9 of 20

Fig. 4. Hydraulic characteristic curves: (a) water retention and (b) hydraulic conductivity of the synthetic dual-permeability (DP) soil; (c) water con-
tents at time zero (q0) and at the end of the tension infiltrometer (TI) experiment (qTI), and bulk saturated water content (qs) for the single-permeability
(SP) soil (matrix alone, SP-1K) and the DP soil (DP-2K), (d) cumulative infiltration into the synthetic DP soil for the TI experiment (ITI) and for the
Beerkan experiment (IB) with the contribution of the matrix (matrix contr.). The letter A refers to comments in the text.

the difference in the water content remains within the range of 5% methods, BEST-2K-A and BEST-2K-B, to error-free analytical data
for DP soils. These results show that comparison of water contents provided sets of hydraulic parameters quite similar to each other,
measured after the TI experiments with the soil porosity could leading to close water retention and hydraulic conductivity func-
help in detecting the occurrence of DP behavior, provided that tions. Given the similarity of the results, we illustrate only the
the measurement uncertainty is low enough (on the order of 1%). results for the case of BEST-2K-A coupled with the BEST Slope
Regarding the cumulative infiltration cur ves, the method of BEST-1K.
contribution of the fast-flow region to the total cumulative The estimated functions are nearly identical to the tar-
infiltration is minimal for the TI experiment, with an alignment gets (Fig. 5a and 5b); in particular, the estimated hydraulic
of the matrix contribution to the bulk cumulative infiltration parameters are on the same order of magnitude as the target
(Fig. 4d, I TI). Conversely, for the Beerkan experiment, the fast- parameters listed in Table 2 (columns Target for the target and
flow region has a large contribution to the bulk cumulative (1) for BEST-2K-A). However, the accuracy of BEST-2K can be
infiltration (Fig. 4d, IB). Thus, the fast-flow region contributes improved. The shapes and, in particular, the bimodality of the
more to cumulative infiltration in the Beerkan experiments target water retention and hydraulic conductivity functions are
than in the TI experiments. Large differences in cumulative not properly depicted because the inflection point of the target
infiltration between the TI and Beerkan experiments could functions is more conspicuous (Fig. 5a and 5b). The analysis of
indicate the occurrence of DP behavior. the hydraulic parameters shows that the saturated water contents
and hydraulic conductivities are properly estimated, with relative
Application of BEST-2K Methods to Synthetic Data errors <10 to 20%, except for the scale parameter hg,f with rela-
A full description of the application of the BEST-2K method tive errors around 110% (Table 2). In fact, the scale parameter
to the analytical data is presented in the supplemental material. for water pressure head |hg | is slightly underestimated for the
For clarity and parsimony, we focus on the discussion of the matrix and largely overestimated for the fast-flow region. As a
BEST-2K results and their accuracy. One observation is that the result, the difference in |hg | between the matrix and fast-flow
three BEST-1K methods (BEST Slope, Intercept, or Steady) lead to regions decreases. This means that when the water pressure head
similar results and estimates. The application of the two BEST-2K increases from hg,m to hg,f, the smallest pores of the fast-flow

VZJ | Advancing Critical Zone Science p. 10 of 20

Fig. 5. Impacts of BEST-2K preprocessing functions on BEST-2K accuracy: (a–d) analysis of the sources of errors for BEST-2K-A by comparison of the
target curves (T) with those estimated by BEST-2K-A (1), of BEST-2K-A without errors related to the computation of initial water contents (2), saturated
water content (3), cumulative infiltrations (4), and without any error related to inputs (5); and influence of input accuracy on BEST-2K estimates for (e–h)
errors in estimation of parameter wf, (i–l) erroneous inputs for the measured water content after the tension infiltrometer (TI) experiment qTI, and (m–p)
erroneous inputs for the bulk saturated water content qs,2K. The large uppercase letters refer to comments in the text. The arrows show the direction of
change of the estimated hydraulic functions when the parameter is increased. The black lines refer to the hydraulic functions estimated by BEST-2K.

region begin to fill with water before the largest pores of the 1. All preprocessing function errors, i.e., the errors due to the
computation of the initial water content (Erq 0), saturated water
matrix. In other words, the more the activation of the matrix
content (Erq s), and cumulative infiltrations (ErI)
and fast-flow regions overlap, as well as their contributions to
the bulk hydraulic conductivity, the less evident is the bimodality. 2. Similar to (1) but without the error Erq 0
Thus, the BEST-2K methods underestimate the contrast between 3. Similar to (1) but without the error Erq s
the matrix and fast-flow regions, i.e., the magnitude of the DP 4. Similar to (1) but without the error ErI
behavior is underestimated. 5. The version with residual errors due to the application of
BEST-1K.
Analysis of Sources of Errors
The discrepancy between the target and estimated functions To perform this assessment, we replaced the inputs provided by
results from the errors that may be induced either by the BEST-2K the preprocessing functions, (q 0,m, q s,m, IB,m and q 0,f, q s,f, IB,f ),
preprocessing functions that provide input to BEST-1K or by the by the real values before processing them with BEST-1K (Fig.
application of BEST-1K itself. Here, we assess the errors due to 2). The real values were directly computed from the synthetic
the preprocessing functions by comparing the target function and hydraulic functions (Eq. [1]) for water contents and the analyti-
related parameters to those obtained with BEST-2K based on the cal models (Eq. [16b] and [16c]) for the cumulative infiltrations.
following scenarios: The removal of Erq 0 does not significantly improve the estimates

VZJ | Advancing Critical Zone Science p. 11 of 20

(Table 2, Scenario 2, Fig. 5c and 5d). However, the removal of Erq s the determination of q s,2K may also be subject to uncertainty.
or ErI (Table 2, Scenario 3 or 4, respectively) leads to a signifi- However, q TI and q s,2K, and the difference between them, are key
cant improvement. Very accurate estimates are obtained when all parameters of the BEST-2K preprocessing functions. Last, the
errors are removed, with residual relative errors mostly below 10% quality of cumulative infiltrations ITI and IB impacts the quality
(Table 2, Scenario 5) and a good agreement between the estimated of the cumulative infiltrations IB,f and IB,m, which in turn impacts
and target functions (Fig. 5c and 5d). In conclusion, the errors the degree of success of the treatment with BEST-1K. This is par-
produced by the preprocessing functions, and in particular Erq s ticularly true for the BEST Slope and BEST Intercept methods,
and ErI, explains most of the discrepancy between the target and which need a proper description of both the transient and steady
BEST-2K estimates. In addition, the residual error due to BEST-1K states. In this case, the problem is not with the BEST-2K approach
can be considered negligible in this particular case. Note that for but with the efficiency with which the BEST-1K methods handle
these results, as stated above, the shape parameters nm and nf are imprecise descriptions of the cumulative infiltrations. This issue
assumed to be perfectly estimated, although these may involve a is addressed below with a discussion of the field experimental data.
supplementary source of error. The investigation of the impact
of erroneous estimations of shape parameters on the accuracy of Sensitivity of BEST-2K Methods to Erroneous Estimates
BEST-2K will be the subject of further studies. of the Volume Ratio Occupied by the Fast-Flow Region
In the previous results, the experimental inputs (PSD, q 0, q TI, The parameter wf is used in most of the parts of the BEST-2K
q s,2K, ITI, and IB) were considered error free, and parameter wf (the methods, including all the preprocessing functions for computing
volume fraction occupied by the fast-flow region) was perfectly the inputs q 0,m, q s,m, IB,m and q 0,f, q s,f, IB,f, and the water retention
known. Here, we first assess the robustness of BEST-2K-A and and hydraulic conductivity functions used for computing q 2K (h)
BEST-2K-B with regard to erroneous estimations of wf, which is a and K2K (q) (Fig. 2). Despite its potentially significant influence,
key parameter for the BEST-2K preprocessing functions and one wf does not impact the water retention and hydraulic conductiv-
that is also complicated to estimate (Kodešová et al., 2010). In addi- ity functions predicted by BEST-2K-A (Fig. 5e and 5f, A). In fact,
tion, we investigate the robustness of these methods with regard the values of the hydraulic parameters vary with wf, yet they pro-
to erroneous inputs of q TI and q s,2K. Indeed, the measurement of duce similar hydraulic functions. However, for BEST-2K-B, wf has
water contents in the field is known to be subject to uncertain- a strong influence on the predicted water retention and hydraulic
ties in the soil moisture probe, which renders the estimation q TI conductivity functions (Fig. 5g and 5h). An increase in wf increases
difficult. Besides, the measure of bulk density that is needed for the bulk saturated hydraulic conductivity Ks,2K and the contrast

Table 2. Impact of preprocessing functions on the accuracy of BEST-2K-A: target hydraulic parameters corresponding to the synthetic soil (Target),
BEST-2K-A estimates (1), and BEST-2K-A estimates without errors in the computation of initial water content (2), final water contents (3), cumulative
infiltration (4), and without any errors on the computation of BEST-1K inputs (5).

(1) (2) (3) (4) (5)

Parameter† Target Erq 0 + Erq s + ErI Erq s + ErI Erq 0 + ErI Erq 0 + Erq s Residual errors
Matrix
q s,m 0.430 0.450 0.450 0.430 0.450 0.430
|hg,m|, mm 144.0 122.6 122.6 131.1 151.2 161.8
Ks,m, mm min−1 1.73 ´ 10−1 1.67 ´ 10−1 1.67 ´ 10−1 1.60 ´ 10−1 1.68 ´ 10−1 1.60 ´ 10−1
Matrix error, %
q s,m – 4.8 4.8 0.0 4.8 0.0
|hg,m| – 14.9 14.9 9.0 5.0 12.3
Ks,m – 3.9 3.9 8.0 3.2 7.6
Fast-flow region
q s,f 0.700 0.516 0.516 0.700 0.516 0.700
|hg,f |, mm 14.9 31.3 31.3 23.5 23.6 16.9
Ks,f, mm min−1 1.62 ´ 101 1.42 ´ 101 1.42 ´ 101 1.55 ´ 101 1.46 ´ 101 1.58 ´ 101
Fast-flow region error, %
q s,f – 26.3 26.3 0.0 26.3 0.0
|hg,f | – 109.8 109.8 57.8 58.2 13.7
Ks,f – 12.6 12.6 4.6 10.1 2.6
† q s, saturated water content; |hg|, scale parameter for water pressure head; Ks saturated hydraulic conductivity. The matrix and fast-flow regions are indicated by
subscripts m and f , respectively.

VZJ | Advancing Critical Zone Science p. 12 of 20

in parameter hg between the matrix and fast-flow regions, which erroneous input for q TI induces errors in the estimates of the
leads to a very conspicuous inflection point and bimodal shape of saturated water content for both the matrix and fast-flow regions,
the water retention and hydraulic conductivity functions (Fig. 5g which in turn impacts the estimates of the other hydraulic param-
and 5h, B). However, when the value of wf is too small, the esti- eters. When q TI increases, the parameter |hg| increases for the
mated hydraulic functions tend toward unimodal curves, which fast-flow region, thus approaching that of the matrix region. As a
correspond to SP soils (Fig. 5g and 5h, C). These results show that result, unimodal curves are obtained (Fig. 5i and 5j, I). For under-
the BEST-2K-B method is very sensitive to the prior estimation of estimated values of q TI, the fast-flow region is predicted as being
the volume ratio wf, in contrast to BEST-2K-A. Therefore, BEST- too dominant, which leads to extremely bimodal hydraulic func-
2K-A should be used when the volume ratio wf cannot be estimated tions and a complete separation of the fast-flow and matrix regions
properly. Another finding is that BEST-2K-A and BEST-2K-B pro- (Fig. 5i and 5j, J). In such a case, the estimated hydraulic functions
vided accurate fits for the cumulative infiltrations for all scenarios. differ significantly from the target functions. Clearly, particular
There was no degradation of fits when parameter wf was varied from attention must be paid to the accuracy of the input parameter q TI
its nominal value. Consequently, the quality of fits cannot be used when using BEST-2K-A.
to identify parameter wf. In general terms, parameter optimization The main conclusion of the sensitivity analysis is that, pro-
may be problematic when many parameters are considered (Angulo- vided the measurement uncertainty does not undermine the
Jaramillo et al., 2016). With reference to the results obtained here, quality of the BEST-2K inputs, either BEST-2K-A or BEST-
additional data are needed to properly estimate parameter wf prior 2K-B may be used. In addition, the estimated hydraulic functions
to running the BEST-2K methods because this parameter cannot be are not significantly different from the target functions. These
estimated by fitting the cumulative infiltration data. results validate the BEST-2K approach numerically. However,
because BEST-2K-B is very sensitive to the prior estimation of
Sensitivity of BEST-2K Methods to Erroneous Inputs w f, this method should not be used when the estimation of w f
for the Bulk Saturated Water Contents is uncertain. On the other hand, the main advantage of BEST-
For the BEST-2K-A method, erroneous values of q s,2K impact 2K-B over BEST-2K-A is that BEST-2K-B does not require the
only the computation of the saturated water content for the fast- measurement of water contents at the end of the TI experi-
flow region q s,f. Inputs q 0m, q s,m, and IB,m remain unchanged (see ments, whereas the quality of the water content measurement
Eq. [4a], [7b], and [12a]), leading to similar hydraulic parameters significantly influences the quality of estimates obtained using
for the matrix. However, the changes are drastic for the fast-flow BEST-2K-A. BEST-2K-B is preferable when there are large uncer-
region. When q s,2K is overestimated, the bulk hydraulic functions tainties in the water content measurements.
exhibit an increase in saturated water content q s,2K and saturated Finally, it should be noted that although BEST-2K-A and
hydraulic conductivity Ks,2K along with a more conspicuous BEST-2K-B use the initial water content q 0, the experiments per-
bimodality (Fig. 5m and 5n, D). On the other hand, when q s,2K is formed with this input are not presented because uncertainty in
underestimated, the hydraulic functions tend toward unimodal q 0 has no impact on the results. Furthermore, q TI and q s,2K are
curves (Fig. 5m and 5n, E). For the BEST-2K-B method, erroneous not subject to the same uncertainties in the field. Indeed, in con-
determinations of q s,2K directly impact the estimates of saturated trast to q TI, q s,2K is not directly measured but rather is inferred
water contents for both the matrix and fast-flow regions, q s,m and from the measurement of bulk density, which is probably subject
q s,f, which in turn affects the estimates of the other hydraulic to less uncertainty. Nevertheless, the value of q s,2K must be precise
parameters. When q s,2K is increased, the bulk saturated water because it affects the BEST-2K-A and BEST-2K-B methods.
content and saturated hydraulic conductivity increase, whereas
the scale parameter |hg| decreases. Consequently, water retention Application of BEST-2K Method to Real Data
functions move upward (Fig. 5o, F) and hydraulic conductivity Experimental Results
functions move toward the left (Fig. 5p, G), which indicates more The experimental data that were collected for the three sites
permeable soils with less capillarity. The bimodality of the curves are depicted in Fig. 6, and the physical parameters are listed in
remains similar across the range of q s,2K values (Fig. 5o and 5p, G). Table 3. The bulk saturated water content q s,2K has high values
The BEST-2K-B method was found to be more sensitive for the for the forest and orchard sites and lower values for the pasture
estimation of the water retention function than for the hydraulic site (Table 3). The volumetric water contents q TI are significantly
conductivity function (Fig. 5o vs. 5p). In terms of bimodality of lower than the saturated water contents, except for the pasture site.
the hydraulic conductivity function, BEST-2K-B is less sensitive The cumulative infiltrations vary with the site, with much larger
to the accuracy of q s,2K than BEST-2K-A. infiltrations for the orchard site and comparable infiltrations for
the forest and pasture sites (Fig. 6b and 6d). This difference may
Sensitivity of BEST-2K Methods to Erroneous Inputs be due to the soil texture and structure, with finer soils at the forest
for Tension Iniltrometer Water Contents and pasture sites and a coarse texture and more open structure at
BEST-2K-B does not use q TI and is therefore entirely insen- the orchard site. The cumulative infiltrations obtained with the
sitive to this parameter (Fig. 5k and 5l, H). For BEST-2K-A, an Beerkan method are much higher than those obtained with the

VZJ | Advancing Critical Zone Science p. 13 of 20

Fig. 6. Experimental data collected for three soils: (a) particle size distributions (PSD); (b) cumulative infiltration obtained with the tension infil-
trometer (ITI); (c) textural triangle and textural properties; (d) cumulative infiltration obtained for the Beerkan experiments (IB), with the cumulative
infiltration for the tension infiltrometer represented with dashed lines for the orchard soil (A).

TI (Fig. 6d, A), which indicates a potential for DP behavior for hydraulic parameters. As an example, Fig. 7a to 7h illustrate some
the three soils (Gonzalez-Sosa et al., 2010); the large increase cor- key findings for the specific case of the Orchard 2 trial. First, the
responds to the activation of the macropore network in these soils bimodal model for the PSD is more accurate than the unimodal
(Lassabatere et al., 2014). model (Fig. 7a vs. 7b), with a more precise modeling of intermedi-
ate points and the two modes around 10 and 103 mm (Fig. 7a, A).
Illustration of the Case of Orchard 2 (Orchard Site) Second, the fit of the model to the experimental infiltration curves
For the analysis of these experimental data, we considered was also accurate (Fig. 7c–7f). For the matrix region, the transient-
only BEST-2K-A, since the volume fraction occupied by the fast- state model was adjusted based on the complete dataset (Fig. 7c
flow region is not known a priori and erroneous estimations of this and 7d, B) because the transient state was longer than the total
parameter may undermine the quality of BEST-2K-B estimates. duration of the experiments. This result points to the difficulty in
However, as discussed above, the BEST-2K-A method requires reaching steady state, which indicates potential estimation errors
a very accurate measure of the water content q TI. In addition, if in the BEST-1K methods (i.e., Slope, Intercept, and Steady) that
BEST-2K-A is able to match the target hydraulic functions with require the attainment of a steady state. In contrast, for the fast-
any value of the parameter wf, the estimates for hydraulic param- flow region, the transient state was extremely short and the steady
eters may still be erroneous for erroneous values of wf (see above). state was attained in a very short time (Fig. 7e and 7f). The tran-
BEST-2K-A was successfully applied to all the trials, and sient state is represented by two points only, including the starting
it provided accurate fits and plausible estimates for most of the point (tB = 0, IB,f = 0) (Fig. 7f, C), whereas the steady-state model

VZJ | Advancing Critical Zone Science p. 14 of 20

matches most of the dataset (Fig. 7e, D). This pattern probably conductivity functions are unimodal without any inflection point
impacts the quality of BEST-1K and thus the accuracy of esti- (Fig. 7p, L). Similar but more pronounced trends are shown for
mations for the fast-flow region (e.g., Di Prima et al., 2016). In the pasture site (Fig. 7m and 7n, J and K). In these cases, BEST-2K
addition, this pattern was often encountered, which highlights predicts values of |hg| that are much larger for the fast-flow region
the necessity of having a precise description of the transient state (Table 4, Orchard 1, Pasture 1, and Pasture 2), which is not physi-
at very short times for the Beerkan test. After the completion of cally realistic. Indeed, larger values of |hg| are associated with soils
the infiltration experiments, the bulk hydraulic parameters and that have smaller pore sizes and thus can increase water retention
hydraulic functions were characterized (Fig. 7g and 7h). The by capillarity (Angulo-Jaramillo et al., 2016).
estimated hydraulic functions clearly show the activation of the Several hypotheses may be evoked for the difficulty in char-
fast-flow region close to saturation, with an increase in both the acterizing half of the trials. First, these trials may not involve any
water content and hydraulic conductivity. The graphs show that fast-flow region or preferential flow. Specific investigations are
the bimodality of the hydraulic conductivity curve is much more needed to examine what happens when BEST-2K is applied to SP
pronounced than that of the water retention curve (Fig. 7h, E). soils. Besides, the difficulty in estimating the hydraulic functions
and parameters may explain a part of the failure. As discussed
Water Retention and Hydraulic Conductivity above, proper estimation of hg for the matrix region requires a
Functions for the Three Sites steady-state water infiltration condition, which may take too long
The water retention and hydraulic conductivity functions to attain. On the other hand, the opposite conditions create diffi-
estimated for the three sites are depicted in Fig. 7i to 7p. The com- culties in estimating the parameter hg for the fast-flow region. Like
parison of the sites shows that the forest and orchard sites have the most permeable porous media with small sorptivity, the transient
highest values of saturated water content (Fig. 7i, F) and hydraulic states are very short, which gives too little time to properly sample
conductivity (Fig. 7j, G). The same hydraulic functions are shown the transient state and leads to errors in estimating the parameter
in Fig. 7k to 7p with the contribution of the matrix region for each hg. In addition, the estimation of the volume fraction occupied
site. For the two tests performed at the forest site, the bulk hydrau- by the fast-flow region, wf, remains problematic, with potential
lic conductivity functions align with the matrix contribution until impacts on the estimation of the hydraulic parameters. However,
an inflection point (Fig. 7l, H), after which the bulk hydraulic despite these shortcomings, the prospects of applying the BEST-2K
functions deviate due to the contribution of the fast-flow region. methods to real experimental data are quite promising. Indeed,
A similar pattern is obtained for Orchard 2, as depicted in Fig. 7o BEST-2K provided good fits with plausible results for several trials
and 7p (M). The inflection points and shifts in the curves reveal performed in this study. In half of the cases, the methods were able
the effect of the activation of the fast-flow region in boosting the to characterize the bimodality of the water retention and hydraulic
bulk water retention and hydraulic conductivity functions. For conductivity functions.
these cases, the BEST-2K method predicts a clear DP behavior
with a higher value of |hg| for the matrix region in comparison to Comparison of BEST-2K and BEST-1K
the fast-flow region (Table 4, Forest 1, Forest 2, and Orchard 2). BEST-1K methods were also used to retrieve water retention
In contrast, for Orchard 1, the bulk water retention and hydraulic and hydraulic conductivity functions corresponding to SP soils
considering the inputs PSD, q 0, q s,2K, and IB. The fit of the PSD
Table 3. Coordinates, clay, silt, and sand content (USDA classification had to be restricted to the first mode, which gave similar results. In
system) in the 0- to 10-cm depth range, soil textural classification, dry all cases, good fits were obtained for the cumulative infiltration IB.
soil bulk density (r d), initial volumetric soil water content (q 0), bulk
saturated volumetric water content (q s,2K), and volumetric water con- However, the derived hydraulic parameters and related water reten-
tent obtained at −30 mm (q TI) for the sampled soils at the pasture tion and hydraulic conductivity curves corresponded to different
(Alghero), forest (Villanova Monteleone), and orchard (Palermo) sites. hydraulic parameters and hydraulic functions. These estimated
Variable Pasture Forest Orchard water retention and hydraulic conductivity curves were used to
Coordinates 40°37¢33.7² N, 40°27¢5.0² N, 38°6¢25.7² N, compute the cumulative infiltration corresponding to an imposed
8°21¢0.4² E 8°30¢47.5² E 13°21¢7.7² E water head of −30 mm and to enable a comparison with the experi-
Clay, % 29.0 (0.39)† 22.2 (0.70) 21.0 (2.02) mental data ITI. The cumulative infiltrations corresponding to
Silt, % 21.6 (2.17) 37.6 (1.00) 23.8 (3.72) the BEST-1K hydraulic parameters were also compared with
Sand, % 49.4 (2.19) 40.2 (1.06) 54.3 (3.67) those obtained with BEST-2K. For the soils that are predicted by
Textural classification sandy clay loam loam sandy loam BEST-2K to have a clear DP behavior (bimodality of the water
1.640 (0.08) 0.842 (0.04) 1.085 (0.08)
retention and hydraulic conductivity curves), the TI cumulative
rd , g cm−3
infiltrations predicted by BEST-1K were less accurate than those
q 0, cm3 cm−3 0.244 (0.02) 0.300 (0.08) 0.110 (0.05)
predicted by BEST-2K. This is particularly the case for Orchard
q s,2K, cm3 cm−3 0.381 (0.03) 0.603 (0.02) 0.584 (0.02) 2. For this case, BEST-1K predicts a very permeable and draining
q TI, cm3 cm−3 0.347 (0.05) 0.312 (0.01) 0.323 (0.02) behavior, typical of coarse soil with very little water retention by
† Standard deviations are indicated in parentheses. capillarity. At −30 mm, many pores are already deactivated and do

VZJ | Advancing Critical Zone Science p. 15 of 20

VZJ | Advancing Critical Zone Science
p. 16 of 20

Fig. 7. (a–h) Application of BEST-2K-A to the Orchard 2 soil: PSD fitted with (a) the bimodal and (b) the unimodal models, (c,d) fit of Beerkan cumulative infiltration into the matrix IB,m and related
infiltration rate qB,m, (e,f ) fit of Beerkan cumulative infiltration into the fast-flow domain IB,f and related infiltration rate qB,f, and (g,h) estimated water retention (WR) and hydraulic conductivity (HC)
functions; (i–p) BEST-2K-A results for all cases, with details for the (k,l) forest, (m,n) pasture, and (o,p) orchard sites. The uppercase letters refer to comments in the text; the dashed lines represent the
contribution of the matrix region.
not conduct water, inducing a drastic drop in cumulative infiltra- precision. The test of the BEST-2K methods with experimental
tion. In contrast, the hydraulic functions predicted by BEST-2K data demonstrated that the cumulative infiltration obtained with
involve a matrix region that sustains enough water infiltration at the tension infiltrometer must be long enough to reach steady state
−30 mm. For the other cases, no logical trends were found. It was for the matrix. For the Beerkan experiment, the occurrence of the
found that BEST-1K fit only the Beerkan data and thus offers the fast-flow region may increase the infiltration rate and reduce the
best fits for Beerkan data, whereas BEST-2K offers the best accu- duration of the transient state at the same time, thus requiring an
racy considering both TI and Beerkan experiments. For real DP experimental setup to enable precise definition of the cumulative
soils, BEST-1K clearly provides a better consistency with regard to infiltration over a very short duration. Last, the analytical calcula-
the modeling of the two cumulative infiltrations. tions proved that any of the BEST methods (Slope, Intercept, or
Steady) should be capable of producing identical results, with the

6 Discussion and Limitations exception of a risk of failure for BEST Slope and BEST Intercept
when the transient state is not properly described.
We have validated the BEST-2K methods using analytically Despite the validation of BEST-2K with both numerical and
generated and real experimental data acquired in the field and experimental data, some limitations remain. First, half of the
highlighted the following crucial points. Under optimal condi- experimental data were not predicted to indicate DP behavior;
tions of use (error-free experimental data), both BEST-2K-A and this finding suggests that either the sampled soils were SP soils or
BEST-2K-B provide results close to the target curves. The analyti- that the BEST-2K method was unable to detect their DP behaviors.
cal sensitivity analysis proves that the robustness of estimates with Several opportunities for improvement emerged. First, we
respect to erroneous inputs varies with the selected method. BEST- need to improve our ability to detect DP behavior. In general, DP
2K-A seems more robust, in particular, with regard to the volume behaviors induce a sharp increase in the hydraulic conductivity
fraction occupied by the fast-flow region, wf. Even though the (Angulo-Jaramillo et al., 2016), with differences of various orders
estimates of the hydraulic parameters may differ from the target of magnitude between unsaturated and saturated hydraulic con-
values, the prediction of the hydraulic functions match closely ductivity (e.g., Watson and Luxmoore, 1986). However, it may be
with the target hydraulic functions. Consequently, BEST-2K-A is difficult to distinguish between SP systems with high saturated
more suitable when wf cannot be determined properly. BEST-2K-A hydraulic conductivity and very low water retention by capillar-
and BEST-2K-B are sensitive to the accuracy of the water measure- ity and soils that effectively exhibit DP behaviors. Indeed, both
ments, including the water content at the end of the TI experiment types of soils may experience a large increase in water content and
and the bulk saturated water content, which is derived from the hydraulic conductivity close to saturation, even if the increase is
bulk density. Water contents should be measured with a minimal expected to be even larger for SP soils (see comparison of BEST-2K
uncertainty. Regarding the cumulative infiltration, the experi- and BEST-1K above). For this objective, more detailed informa-
mental devices should be chosen so as to allow sufficient time to tion provided by multitension experiments could enhance our
reach steady state while describing the transient state with enough ability to distinguish between SP and DP systems, as suggested by

Table 4. Values of the hydraulic parameters† for the sampled soils at the forest, pasture, and orchard sites.
Site Domain w qr qs n h Ks |hg|
% 3 −3
———— cm cm ———— mm min−1 mm
Forest 1 matrix 61.3 0 0.509 2.198 13.1 3.59 ´ 10−2 98.5
fracture 38.7 0 0.752 2.230 11.7 2.14 ´ 10 0 12.5
Forest 2 matrix 61.3 0 0.509 2.198 13.1 3.56 ´ 10−2 76.2
fracture 38.7 0 0.752 2.230 11.7 2.89 ´ 10 0 10.1
Pasture 1 matrix 93.8 0 0.370 2.127 18.8 6.54 ´ 10−2 9.23
fracture 6.2 0 0.552 2.688 5.91 3.58 ´ 10 0 149.8
Pasture 2 matrix 93.8 0 0.370 2.127 18.8 5.86 ´ 10−2 4.92
fracture 6.2 0 0.552 2.688 5.91 3.26 ´ 10 0 260.7
Orchard 1 matrix 72.9 0 0.497 2.172 14.7 2.50 ´ 10 0 16.5
fracture 27.1 0 0.715 2.384 8.21 4.10 ´ 101 29.0
Orchard 2 matrix 62.6 0 0.496 2.172 14.6 1.27 ´ 10 0 35.0
fracture 37.4 0 0.740 2.384 8.36 1.39 ´ 101 2.68
† w, volume percentage occupied by the medium; q r, residual water content; q s, saturated water content; n, shape parameter for the water retention curve; h, exponent
of the relative hydraulic conductivity; Ks, saturated hydraulic conductivity; |hg|, scale parameter for water pressure head

VZJ | Advancing Critical Zone Science p. 17 of 20

Lassabatere et al. (2014). Further research, including numerical obtained for the matrix and fast-flow regions. BEST-1K meth-
or experimental investigations, is needed to properly define the ods may lead to unreliable results in the case of nonattainment
values of the water pressure heads to be imposed for the detection of steady-state conditions or inappropriate descriptions of the
of DP behaviors. transient state. This was probably the case with our experimental
Second, the accuracy of the BEST-2K method that estimates data. We obtained higher values of the scale parameter for water
the shape parameters from the bulk PSD should be investigated pressure head |hg| for the fast-flow region, which is not physically
more deeply. The impact of erroneous estimates of the shape realistic because it suggests that more water retention by capil-
parameter on the bulk hydraulic curves should be investigated larity occurs in the fast-flow region than in the matrix region.
in more detail. In addition to the errors resulting from the PTFs, The imprecise descriptions of cumulative infiltration, with a poor
decomposing the PSD into bimodal curves and assigning each description of the transient state for the fast-flow region, and the
mode to the matrix or fast-flow regions may be questionable. nonattainment of a steady state in the matrix region are probably
Indeed, if it is clear that soil structure may be induced by the soil the explanations for the errors found in estimation of the scale
texture (with the development of larger pores around the largest parameter |hg|. In future studies, the use of precise, simple, and
particles), other external factors may impact the soil structure. inexpensive automatic monitoring systems, such as that proposed
Bioturbation and shrinkage due to wetting–drying or freeze– by Di Prima (2015), could improve the quality of cumulative
thaw cycles may also create macropore systems, even in a soil that infiltration measurements over short durations. Similar devices
exhibit a unimodal PSD. In this situation, the derivation of shape could be designed for the application of tension infiltrometers
parameters from the PSD may be questionable, particularly for that would allow infiltration experiments to last long enough to
the fast-flow region. attain steady-state conditions.
Third, some parameters of the DP system have yet to be esti-
mated. The proposed methodology does not provide a specific way
to estimate the volume fraction occupied by the fast-flow region,
6 Conclusion
wf. The determination of wf remains very tricky and may require This study developed a new method (BEST-2K) for the char-
additional experiments. Micromorphological images could be acterization of DP soils on the basis of the BEST-1K methods
used to estimate wf from image analysis (Kodešová et al., 2009). previously developed for SP soils. BEST-2K needs only raw data
Tracer experiments may also be used to distinguish the porosity (PSD and initial and final water contents) and the cumulative
that is easily accessible by solutes (which correspond to the fast- infiltrations obtained at two different water pressure heads, i.e.,
flow region) from the inaccessible porosity (which correspond to −30 mm to activate only the matrix and 0 mm to activate the
stagnant water zones in the matrix) (Kodešová et al., 2010, 2012). entire pore network. With these data, BEST-2K provides a full
The use of dyes in the field may also help in the detection of pref- characterization of DP soils. Two methods were presented. BEST-
erential pathways and the characterization of the volume of soil 2K-A uses the two water contents measured at −30 and 0 mm to
affected by preferential flow (Cey and Rudolph, 2009; Cey et al., derive the saturated water contents of the matrix and fast-flow
2009). Besides wf, the BEST-2K methods also do not provide any regions. Accurate water measurements are then required. On the
estimate of the interfacial hydraulic conductivity Ks,a, which deter- other hand, BEST-2K-B is a simpler alternative that equates both
mines the water exchange between the matrix and fast-flow regions. local water contents to the soil bulk porosity. BEST-2K-B may
In fact, Lassabatere et al. (2014) suggested that water exchange be used when the uncertainties in the water measurements are
does not change water infiltration into the DP porous medium too high. The two methods were validated using both analyti-
for any water pressure head imposed at the soil surface. In other cal and real experimental data. In the absence of clear evidence
words, whatever the value of Ks,a, the TI and Beerkan experimental for the selection of one or the other method, the use of both
data are identical, which leads to similar characterization by the methods and subsequent merging of the obtained results may
BEST-2K methods. However, the numerical results of Lassabatere be recommended.
et al. (2014) require additional validation by experimental data. The tests of the two proposed methods highlighted their
Further numerical and experimental studies are needed (i) to inves- strengths and weaknesses. Notwithstanding the required research
tigate more deeply the effect of nonequilibrium between the matrix for additional improvements, the BEST-2K method marks a useful
and fast-flow regions and related effects on water infiltration, and first step toward the characterization of DP soils. The method
(ii) to design a proper strategy to estimate Ks,a. As for wf, additional requires only simple water infiltration tests and soil characteristics,
laboratory and field experiments, such as tracer or dye experiments, thus minimizing time and cost requirements and providing an
may be necessary to better describe the water exchange between exhaustive characterization of DP soils. Its association with the
the matrix and fast-flow regions to obtain more insight on Ks,a regular BEST-1K methods may offer a very interesting tool for
and the effect of nonequilibrium between these regions on the the hydraulic characterization of soils prone to preferential flow.
flow processes. Further research activities will be conducted to enhance this
Last, the quality of BEST-2K also strongly depends on approach for the detection and quantification of preferential flow
the efficiency of BEST-1K in treating cumulative infiltrations in soils.

VZJ | Advancing Critical Zone Science p. 18 of 20

List of Symbols h m, h f shape parameter of the local hydraulic
BEST-1K parameters conductivity function for the matrix and the
fast-flow regions
qr residual water content
rd, r s dry bulk density and soil particle density
qs saturated water content
PSD bulk particle size distribution
Ks saturated hydraulic conductivity PSDm, PSDf local particle size distribution for the matrix and
hg scale parameter for water pressure head the fast-flow regions
n shape parameter for water retention function FF2K bulk cumulative particle size distribution
h shape parameter for hydraulic conductivity tf fraction related to the contribution of the
function fast-flow region to the bulk particle size
S sorptivity distribution
e porosity Dg,m, Dg,f average diameter of the particles of the matrix
s fractal dimension and the fast-flow regions
Nm, Nf textural parameters for the matrix and the fast-
pm shape index related to water retention function flow regions
pM textural shape index related to particle size qTI bulk water content at the end of TI experiments
distribution
p tortuosity parameter qTI,m, qTI,f local water content in the matrix and the fast-
flow regions at the end of TI experiments
cp shape parameter for cumulative infiltration
q 0,TI, q 0,B bulk initial water content for TI and Beerkan
IO(2)(t) transient-state approximate expansion for experiments
cumulative infiltration q 0,m, q 0.f local initial water contents for the matrix and
I+¥(t) steady-state approximate expansion for the fast-flow regions
cumulative infiltration IB bulk cumulative infiltration related to Beerkan
qO(2)(t) transient-state approximate expansion for experiment
infiltration rate tB time dataset for the Beerkan experiment
q+¥(t) steady-state approximate expansion for
infiltration rate ITI bulk cumulative infiltration related to TI
experiment
A, B, C coefficients for the approximate expansions
t TI time dataset for the TI experiment
qsexp steady-state infiltration rate (slope of the steady
state asymptote) IB,m, IB,f cumulative infiltration for Beerkan experiments
sampling only the matrix and the fast-flow
bsexp intercept of the steady-state asymptote regions
Sopt optimized value for sorptivity IB,m(tB) cumulative infiltration for Beerkan experiments
Ks,opt optimized value for saturated hydraulic sampling only the matrix computed for time
conductivity dataset tB
tmax maximum time for the validity of the transient- I3D,2K (t) model for the computation of cumulative
state approximate expansion infiltrations into dual-permeability soils for TI
BEST-2K parameters and Beerkan experiments
wf volume fraction of the dual-permeability soil rd radius of the source for water infiltration
occupied by the fast-flow region experiments
q2K bulk water content for the dual-permeability soil Sm , Sf computed sorptivities for the matrix and the
fast-flow regions
q s,2K bulk saturated water content for the dual-
permeability soil Dm , D f water diffusivity for the matrix and the fast-flow
regions
q m, q f local water contents in the matrix and the fast-
flow regions b m, g m, bf, g f infiltration constants related to quasi-exact
implicit model and related approximate
q s,m, q s,f local saturated water contents for the matrix and expansions for the matrix and fast-flow regions
the fast-flow regions
rg,m, rg,f local averaged pore size for the matrix and the
q r,m, q r,f local residual water contents for the matrix and fast-flow regions
the fast-flow regions
saw surface tension of the air/water interface
K2K bulk hydraulic conductivity for the dual-
permeability soil bc contact angle
Ks,2K bulk saturated hydraulic conductivity for the rw water density
dual-permeability soil g gravitational acceleration constant
Km, Kf local hydraulic conductivities for the matrix and Erq 0, Erq s, errors of BEST-2K preprocessing function for
the fast-flow regions ErI the computation of initial water contents,
Ks,m, Ks,f local saturated hydraulic conductivities for the saturated water contents, and cumulative
matrix and the fast-flow regions infiltration, respectively
hg,m, hg,f scale parameter for water pressure head for the
matrix and the fast-flow regions Acknowledgments
nm, nf shape parameter of the local water retention We wish to thank the French National Research Agency (ANR) for its contribu-
function of the matrix and the fast-flow tion to the funding of the INFILTRON Project (ANR-17-CE04-0010) and, in
regions particular, this work and for providing scientific supervision of the research. We
also thank the Field Observatory in Urban Water Management (OTHU) for tech-
nical and scientific support. We would like to thank the associate editor, Majdi
Abou Najm, and the reviewers for their excellent suggestions and contributions.

VZJ | Advancing Critical Zone Science p. 19 of 20

These were crucial for the improvement of the manuscript, and we greatly appreci- Kodešová, R., K. Němeček, V. Kodeš, and A. Žigová. 2012. Using dye tracer
ate their advice. We also thank R. Marrosu for his assistance in the field activity. for visualization of preferential flow at macro- and microscales. Vadose
Zone J. 11(1). doi:10.2136/vzj2011.0088
Kodešová, R., J. Šimůnek, A. Nikodem, and V. Jirků. 2010. Estimation
References of the dual-permeability model parameters using tension disk
Arya, L.M., and J.F. Paris. 1981. A physicoempirical model to pre- infiltrometer and Guelph permeameter. Vadose Zone J. 9:213–225.
dict the soil moisture characteristic from particle-size distri- doi:10.2136/vzj2009.0069
bution and bulk density. Soil Sci. Soc. Am. J. 45:1023–1030. Kodešová, R., N. Vignozzi, M. Rohošková, T. Hájková, M. Kočárek, M. Pagliai,
doi:10.2136/sssaj1981.03615995004500060004x et al. 2009. Impact of varying soil structure on transport processes
Angulo-Jaramillo, R., V. Bagarello, M. Iovino, and L. Lassabatere. 2016. in different diagnostic horizons of three soil types. J. Contam. Hydrol.
Infiltration measurements for soil hydraulic characterization. Springer, 104:107–125. doi:10.1016/j.jconhyd.2008.10.008
Berlin. doi:10.1007/978-3-319-31788-5 Kutilek, M., and D.R. Nielsen. 1994. Soil hydrology. Catena Verlag, Cremlin-
Bagarello, V., S. Di Prima, and M. Iovino. 2014. Comparing alternative gen, Germany.
algorithms to analyze the Beerkan infiltration experiment. Soil Sci. Soc. Lassabatere, L., R. Angulo-Jaramillo, J.M. Soria Ugalde, R. Cuenca, I.
Am. J. 78:724–736. doi:10.2136/sssaj2013.06.0231 Braud, and R. Haverkamp. 2006. Beerkan estimation of soil transfer
Braud, I., D. De Condappa, J.M. Soria, R. Haverkamp, R. Angulo-Jaramillo, parameters through infiltration experiments: BEST. Soil Sci. Soc. Am. J.
S. Galle, and M. Vauclin. 2005. Use of scaled forms of the infiltra- 70:521–532. doi:10.2136/sssaj2005.0026
tion equation for the estimation of unsaturated soil hydraulic Lassabatere, L., R. Angulo-Jaramillo, J.M. Soria-Ugalde, J. Šimůnek,
properties (the Beerkan method). Eur. J. Soil Sci. 56:361–374. and R. Haverkamp. 2009. Numerical evaluation of a set of
doi:10.1111/j.1365-2389.2004.00660.x analytical infiltration equations. Water Resour. Res. 45:W12415.
Brooks, R., and A. Corey. 1964. Hydraulic properties of porous media. doi:10.1029/2009WR007941
Hydrol. Pap. 3. Colo. State Univ., Fort Collins. Lassabatere, L., D. Yilmaz, X. Peyrard, P.E. Peyneau, T. Lenoir, J. Šimůnek,
Burdine, N.T. 1953. Relative permeability calculations from pore size distri- and R. Angulo-Jaramillo. 2014. New analytical model for cumula-
bution data. J. Pet. Technol. 5(03):71–78. doi:10.2118/225-G tive infiltration into dual-permeability soils. Vadose Zone J. 13(12).
doi:10.2136/vzj2013.10.0181
Carsel, R.F., and R.S. Parrish. 1988. Developing joint probability distri-
butions of soil water retention characteristics. Water Resour. Res. Lenhard, R.J., J.H. Dane, and M. Oostrom. 2005. Immiscible fluids. In: D. Hil-
24:755–769. doi:10.1029/WR024i005p00755 lel, editor, Encyclopedia of soils in the environment. Elsevier, Amster-
dam. p. 239–247.
Cey, E.E., and D.L. Rudolph. 2009. Field study of macropore flow processes
using tension infiltration of a dye tracer in partially saturated soils. Parlange, J.-Y. 1975. On solving the flow equation in unsaturated soils by
Hydrol. Processes 23:1768–1779. doi:10.1002/hyp.7302 optimization: Horizontal infiltration. Soil Sci. Soc. Am. J. 39:415–418.
doi:10.2136/sssaj1975.03615995003900030019x
Cey, E.E., D. Rudolph, and J. Passmore. 2009. Influence of macroporosity
on preferential solute and colloid transport in unsaturated field soils. J. Rieu, M., and G. Sposito. 1991. Fractal fragmentation, soil porosity, and soil
Contam. Hydrol. 107:45–57. doi:10.1016/j.jconhyd.2009.03.004 water properties: II. Applications. Soil Sci. Soc. Am. J. 55:1239–1244.
doi:10.2136/sssaj1991.03615995005500050007x
Di Prima, S. 2015. Automated single ring infiltrometer with a low-cost
microcontroller circuit. Comput. Electron. Agric. 118:390–395. Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. ROSETTA:
doi:10.1016/j.compag.2015.09.022 A computer program for estimating soil hydraulic parameters
with hierarchical pedotransfer functions. J. Hydrol. 251:163–176.
Di Prima, S., L. Lassabatere, V. Bagarello, M. Iovino, and R. Angulo- doi:10.1016/S0022-1694(01)00466-8
Jaramillo. 2016. Testing a new automated single ring infiltrom-
eter for Beerkan infiltration experiments. Geoderma 262:20–34. Smettem, K.R.J., J.-Y. Parlange, P.J. Ross, and R. Haverkamp. 1994.
doi:10.1016/j.geoderma.2015.08.006 Three-dimensional analysis of infiltration from the disc infiltrom-
eter: 1. A capillary-based theory. Water Resour. Res. 30:2925–2929.
Durner, W. 1994. Hydraulic conductivity estimation for soils with doi:10.1029/94WR01787
heterogeneous pore structure. Water Resour. Res. 30:211–223.
doi:10.1029/93WR02676 Sutera, S.P., and R. Skalak. 1993. The history of Poiseuille’s law. Annu. Rev.
Fluid Mech. 25:1–20. doi:10.1146/annurev.fl.25.010193.000245
Gee, G.W., and J.W. Bauder. 1986. Particle-size analysis. In: A. Klute,
editor, Methods of soil analysis. Part 1. Physical and mineralogical Timlin, D.J., L.R. Ahuja, and M.D. Ankeny. 1994. Comparison
methods. SSSA Book Ser. 5. SSSA and ASA, Madison, WI. p. 383–411. of three field methods to characterize apparent mac-
doi:10.2136/sssabookser5.1.2ed.c15 ropore conductivity. Soil Sci. Soc. Am. J. 58:278–284.
doi:10.2136/sssaj1994.03615995005800020003x
Gerke, H.H., and M.Th. van Genuchten. 1993. A dual-porosity model
for simulating the preferential movement of water and solutes van Genuchten, M.Th. 1980. A closed-form equation for predicting
in structured porous-media. Water Resour. Res. 29:305–319. the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J.
doi:10.1029/92WR02339 44:892–898. doi:10.2136/sssaj1980.03615995004400050002x

Gonzalez-Sosa, E., I. Braud, J. Dehotin, L. Lassabatere, R. Angulo-Jaramillo, Watson, K.W., and R.J. Luxmoore. 1986. Estimating macroporosity in a
M. Lagouy, et al. 2010. Impact of land use on the hydraulic properties forest watershed by use of a tension infiltrometer. Soil Sci. Soc. Am. J.
of the top soil in a French catchment. Hydrol. Processes 24:2382–2399. 50:578–582. doi:10.2136/sssaj1986.03615995005000030007x
doi:10.1002/hyp.7640 Yilmaz, D., L. Lassabatere, R. Angulo-Jaramillo, D. Deneele, and M. Legret.
Haverkamp, R., P.J. Ross, K.R.J. Smettem, and J.-Y. Parlange. 1994. Three- 2010. Hydrodynamic characterization of basic oxygen furnace
dimensional analysis of infiltration from the disc infiltrometer: 2. slag through an adapted BEST method. Vadose Zone J. 9:107–116.
Physically-based infiltration equation. Water Resour. Res. 30:2931– doi:10.2136/vzj2009.0039.
2935. doi:10.1029/94WR01788

VZJ | Advancing Critical Zone Science p. 20 of 20