water
ArticleNumerical Analysis of Soil Water Dynamics during Spinach
Cultivation in a Soil Column with an Artificial Capillary Barrier
under Different Irrigation Managements
Davy Sao 1 , Hirotaka Saito 1, * , Tasuku Kato 1
12*

Citation: Sao, D.; Saito, H.; Kato, T.;Šimůnek, J. Numerical Analysis ofSoil Water Dynamics during Spinach
Cultivation in a Soil Column with an
Artificial Capillary Barrier underDifferent Irrigation Managements.
Water 2021, 13, 2176. https://doi.org/10.3390/w13162176
Academic Editor: Guido D’Urso

and Jirka Šimůnek 2

Department of Agricultural and Environmental Engineering, United Graduate School of Agricultural Science,
Tokyo University of Agriculture and Technology, 3-8-5 Saiwai-cho, Fuchu-shi, Tokyo 183-8538, Japan;
davy_sao@hotmail.com (D.S.); taskkato@cc.tuat.ac.jp (T.K.)
Department of Environmental Science, University of California, Riverside, CA 92508, USA; jsimunek@ucr.edu
Correspondence: hiros@cc.tuat.ac.jp

Abstract: Artificial capillary barriers (CBs) are used to improve root zone conditions as they cankeep water and nutrients in the root zone by limiting downward percolation. Numerical analysis is
one of the promising tools for evaluating CB systems’ performance during the cultivation of leafy
vegetables. This study aims to investigate the effects of the CB system on soil water dynamics during
spinach cultivation in a soil column under different irrigation scenarios using HYDRUS (2D/3D)
by comparing uniform (UNI), line-source (LSI), and plant-targeted (PTI) irrigations combined withalternative irrigation schedules. Simulation results of volumetric soil water contents were generally
corresponding to measured data. Simulation results with various hypothetical irrigation scenarios
exhibited that the CB was an effective system to diminish percolation losses and improve the rootzone’s soil water storage capacity. On the other hand, evaporation loss can be increased as morewater is maintained near the surface. While this loss can be significantly minimized by reducing thewater application area, the irrigation amount must be carefully defined because applying water in asmaller area may accelerate downward water movement so that the water content at the CB interface
can reach close to saturation. In addition to the malfunction of the CB layer, it may also cause areduction of plant root water uptake (RWU) because the root zone is too wet. Among evaluatedirrigation scenarios, irrigating every two days with PTI was the best scenario for the spinach as water
use efficiency was greatly improved.
Keywords: HYDRUS (2D/3D); capillary barrier; irrigation scenarios; soil water dynamics; root water
uptake; water use efficiency

Received: 23 May 2021
Accepted: 5 August 2021
Published: 9 August 2021

1. Introduction
Publisher’s Note: MDPI stays neutral
with regard to jurisdictional claims inpublished maps and institutional affiliations.

Copyright: © 2021 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms andconditions of the Creative Commons
Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Global demand for agricultural production, which is projected to rise as populationand incomes are increasing [1], will directly affect agricultural water usage. Since water
resources have been increasingly under pressure due to competition with other sectors,
climate variability, socio-economic changes, overpopulation, together with unsustainableagricultural practices will put greater pressure on these limited resources [1]. To maintain or
increase agricultural production, especially in water-scarce regions where water shortage isa major concern, it is important to develop reliable, simple, and cost-effective technologies
to improve crop water use efficiency. An artificial capillary barrier (CB), a natural barrierfor soil water flow under unsaturated conditions, is one such option as it can improve cropwater use efficiency when properly designed and installed [2,3].
A CB happens at the joint boundary between two soil layers with different hydraulicproperties, such as an interface of a fine-textured soil overlaid a coarse-textured soil.
Because of differences in pore size distributions, capillary forces in the finer layer are greaterthan those in the coarser layer. When a finer particle soil overlays a coarser particle soil, thisWater 2021, 13, 2176. https://doi.org/10.3390/w13162176

https://www.mdpi.com/journal/water

Water 2021, 13, 2176

2 of 18

phenomenon causes downward movement of infiltrating water to stop temporally at theinterface until the matric potential at the interface reaches a value at which the coarser layer
becomes more conductive. The effects of capillary barriers on water flow in soil profiles
have been extensively studied. Stormont and Morris [4] and Morris and Stormont [5] used
CB systems in landfill covers to redirect and prevent water from infiltrating into a sensitive
sub-layer. Ityel et al. [2] installed a CB system to regulate infiltration rate and increase theroot zone water retention capacity. The system was also used by Rooney et al. [3] to avoid
the capillary rise of saline water from contaminated lower soil layer into the topsoil byplacing a CB layer in between.
Numerical simulations have been recognized as a cost-effective tool and are increasingly used to examine soil water dynamics for optimizing irrigation management [6,7].
Among numerical models, the HYDRUS models are one of the most reliable and complete
software packages that have been widely used and thus validated over the years [8]. The
HYDRUS models have been successfully employed to study the effects of CBs on unsaturated soil water flow processes. For instance, Ityel et al. [2] confirmed HYDRUS (2D/3D)
accuracy in assessing soil water contents in the CB soil system and used the model toidentify the impacts of the CB on improving root zone conditions for the cultivation ofhorticultural crops. Kämpf et al. [9] fruitfully used HYDRUS-2D to conceptualize andparameterize flow processes in the CB system. Fala et al. [10] employed HYDRUS-2D to numerically simulate flow and the influence of CB in unsaturated waste rock piles. Henry [11]applied the HYDRUS-2D model to evaluate the degradation of CB performance due tothe contaminant-induced surface tension reduction. Wongkaew et al. [12] confirmed theperformance of HYDRUS-1D in predicting soil water contents and evaluated the impacts ofalternative irrigation schedules under uniform irrigation on root zone moisture conditions
in a CB soil column. Most recently, Al-Mayahi et al. [13] studied the influences of soil
substrates made of Smart Capillary Barrier Wick, containing silt loam blocks enclosed bysand-sheaths and irrigated with a sand wick cylinder.
HYDRUS (2D/3D) has also been a popular model for assessing various irrigation
systems [8]. The 3D model has been extensively used in modeling line-source andpoint-source irrigations in both homogeneous and heterogeneous soils. For instance,
Elmaloglou et al. [14] successfully used the 2D model to estimate soil water fluxes undersurface drip irrigation from equidistant line sources. Shakhali et al. [15] used HYDRUS-2D
to simulate changes in volumetric soil water contents (VWCs) underline sources withdifferent irrigation water salinities and reported an excellent performance of the modelcompared with field observations. Arbat et al. [16] also confirmed the reliability of themodel to assess soil water dynamics, such as deep drainage and plant water uptake forvarious dripline depths, soil types, and irrigation schedules. Honari et al. [17] also reported
high accuracy of the 3D model in modeling soil water dynamics under subsurface drip
irrigation. Dabach et al. [18] successfully used HYDRUS (2D/3D) to evaluate the optimal
placement of tensiometer for high-frequency subsurface drip irrigation in heterogeneous
soils, while Mubarak et al. [19] demonstrated the accuracy of the 2D model in the study on
the influences of high-frequency pulsing of drip irrigation in heterogeneous soils. However, the CB system’s effects under different irrigation management (e.g., line-source andplant-targeted irrigations) have not yet been investigated.
This study investigates the effects of different irrigation practices, including uniform,line-source, and plant-targeted irrigations, on soil water dynamics during Japanese mustard
spinach (Brassica rapa var. perviridis) cultivation in a soil profile with a CB layer installed.
Alternative irrigation schedules were also evaluated to determine whether water use
efficiency could be further improved while maintaining root water uptake by the plant.
In a previous study, HYDRUS-1D successfully predicted soil water dynamics in the CB
soil profile under spinach cultivation [12]. However, the one-dimensional model cannot
simulate spatially non-uniform irrigations, such as line-source irrigation and plant-targeted
irrigation. To achieve these objectives, numerical simulations for different irrigation types
and alternative irrigation schedules were performed in a 3D domain. The applicability

Water 2021, 13, 2176

CB soil profile under spinach cultivation [12]. However, the one-dimensional model cannot simulate spatially non-uniform irrigations, such as line-source irrigation and planttargeted irrigation. To achieve these objectives, numerical simulations for different irrigation types and alternative irrigation schedules were performed in a 3D domain. The3 of
ap18
plicability of the HYDRUS (2D/3D) model was first evaluated by simulating changes involumetric soil water contents (VWCs) and water fluxes in the experimental CB systemof the the
HYDRUSmodel
wasvarious
first evaluatedby hypotheticalsimulating changes
in volumetric
underspinach(2D/3D)
cultivation
beforeadditional
alternative
irrigationsoil waterwere
contents(VWCs) and water fluxes in the experimental CB system under the
scenarios
simulated.
spinach cultivation before various additional hypothetical alternative irrigation scenarios
weresimulated.
2.
Materials
and Methods
2.1. Cultivation Experiment
2. Materials and Methods
The cultivation
experiment was carried out from 17 December 2013, to 14 January
2.1. Cultivation
Experiment
2014, in a controlled environment with a temperature of 25 °C during daytime and 15 °C
The cultivation experiment was carried out from 17 December 2013, to 14 January 2014,
during nighttime [20]. Japanese spinach, one of the◦ common leafy vegetables
widelyin a controlled environment with a temperature of 25 C during daytime and 15 ◦ C during
grown in Japan, was cultivated in four containers with a dimension of 24 cm wide × 36 cm
nighttime [20]. Japanese spinach, one of the common leafy vegetables widely grown in Japan,
long × 23 cm deep. One container, used as a reference, was filled with Tottori Dune sand,was cultivated in four containers with a dimension of 24 cm wide × 36 cm long × 23 cm deep.
having an average particle size of 0.35 mm with 85% of particles between 0.106 to 0.425
One container, used as a reference, was filled with Tottori Dune sand, having an averagemm, to a total depth of 16 cm (Figure 1b). Three other containers contained a CB system,
particle size of 0.35 mm with 85% of particles between 0.106 to 0.425 mm, to a total depth ofas
in Figure
1a. A other
CB layer,consistingof a 2-cm
of coarse
sandinwith
the1a.
par16shown
cm (Figure1b). Threecontainers
containeda CBlayer
system,as shown
FigureA
ticle
size between
0.85
2.00
mmofanda 4-cmlayer
gravelwithsizethebetween
particle 0.85
size and
beCB layer,
consistingof anda 2-cmlayer
coarsesand
withofthe
particletween
4.75and
anda9.00
was
top andTottoriDune
sandlayers2.00 mm
4-cmmm,layerof placedgravel between
with the theparticlesizebottom
between
4.75 and
9.00
mm,
was
makingathickness
of
6
cm
of
the
Dune
sand
at
the
top
for
plant
root
development
4
placed between the top and bottom Tottori Dune sand layers making a thickness of 6 and
cm of
cm
of
the
Dune
sand
remainingatthebottom.
The
top
6-cm
Tottori
Dune
sand
layer
in
the Dune sand at the top for plant root development and 4 cm of the Dune sand remaining
both
were
mixed
with
fertilizers
seeding.
A weekafter
at thereferencebottom. and
TheCB
topcontainers
6-cm Tottori
Dune
sand
layer
in both before
referenceand CB
containers
seeding,liquidfertilizers
wereadditionallyweek together
with irrigation
were mixed
withfertilizers
beforeseeding.applied
A weektwice
after aseeding,liquid fertilizers
werewater.
additionally applied twice a week together with irrigation water.

Figure 1. Soil profiles in the cultivation experiment of Miyake et al. [20] with (a) and without (b) a
Figure 1. Soil profiles in the cultivation experiment of Miyake et al. [20] with (a) and without (b) acapillary barrier layer.
capillary barrier layer.

Irrigation at different rates was applied daily to each container uniformly at the surface
Irrigation at different rates was applied daily to each container uniformly at the surusing a sprinkling pot. The amount of water applied to the three containers with CB wasface
a sprinkling
pot.day,
Thehereinafter
amount of called
water CB15,
appliedto thethreecontainers
with TheCB
0.15,using
0.30, and
0.45 cm per
CB30,
and
CB45,
respectively.
was
0.15,
0.30,
and
0.45
cm
per
day,
hereinafter
called
CB15,
CB30,
and
CB45,
respectively.
irrigation rate at the reference container (referred to as Non-CB30) was 0.30 cm per day.
The
rate
thecm
referencecontainer(referred
to asday
Non-CB30)
wasafter
0.30(thecm first
per
The irrigationwater depth
of at
0.50
was appliedon boththe seeding
and the day
day.
The water
depth
of 0.50 cm
was appliedon boththe seeding
and attheallday
after
day after
seeding)to promote
germination,while
no irrigation
wasdayapplied
plots
on
(the
first
day
after
seeding)topromotegermination,
while
noirrigation
wasappliedatallthe second day after seeding. The CB45 data were not analyzed further in this study since
the irrigation amount in the CB45 was too high and caused the break in a CB, as reportedby Miyake et al. [20] and also confirmed by the model simulation.
2.2. Numerical Simulations
The HYDRUS (2D/3D) model [8,21,22], a computer software package for simulatingwater, solute, and heat movement in variably saturated media, was used in this study

Water 2021, 13, 2176

4 of 18

to evaluate spatial and temporal changes in volumetric soil water contents, soil water
fluxes, and root water uptakes (RWU) under different irrigation designs and schedules
in a three-dimensional (3D) domain. HYDRUS (2D/3D) numerically solves the Richards
equation for saturated and unsaturated water flow in the 3D domain based on a finite
element scheme. A sink term is included in the flow equation to consider plant root water
uptake, which can be calculated using the Feddes RWU model [23]. Details about model
simulation steps, as well as others required input data/parameters, are discussed below.
2.2.1. Soil Hydraulic Properties
In this study, similarly as in Wongkaew et al. [12], the single-porosity (SP) model of van
Genuchten [24] was used to describe soil hydraulic properties of Totorri Dune sand, while
the dual-porosity (DP) model of Durner [25] based on the van Genuchten [24] model wasused for coarse sand and gravel. The SP model may not be appropriate to describe waterflow in the dry range of relatively coarse materials, where water flows as non-capillary typeflow after capillary water is drained. Therefore, the DP model, which can be used to mimic
non-capillary type flow by suppressing the sharp decrease in the hydraulic conductivity inthe dry range, was employed for coarse sand and gravel in this study.
To determine soil water retention curves of these materials, a modified suction table method [26] was used to measure water contents at predetermined matric potential heads ranging from 0 to −100 cm (Figure 2a). A pore size distribution model of
Mualem [27] was used in both models to estimate unsaturated hydraulic conductivities (Figure 2b). The soil hydraulic properties, i.e., the soil water retention curve, Se(h),and the unsaturated hydraulic conductivity function, K(h), for the SP model of van
Genuchten [24] and the DP model of Durner [25] were described using Equations (1) and (2)and Equations (3) and (4), respectively.
Se ( h ) =θ − θr=θ s − θr1
m ; h < 0
1 + |αh|nh m i2K (h) = Ks Sel 1 − 1 − Se1/m
Se ( h ) =K ( h ) = Ks

θ − θr1
1= w1
 m1 + w2
mn1
θ s − θr
1 + | α1 h |
1 + | α 2 h | n2 2 m1 m2 2
1
2+ w2 α2 1 − 1 − Se1/m(w1 Se1 + w2 Se2 )l w1 α1 1 − 1 − Se1/m
2
1( w1 α 1 + w2 α 2 ) 2

(1)

(2)
(3)

(4)

where Se is the effective saturation (-), θ is the volumetric soil water content (VWC)(cm3/cm3), θr and θs are the residual and saturated volumetric soil water contents (cm3/cm3),
respectively, α is the shape coefficient (1/cm > 0), h is the soil water pressure head (cm),n is a pore-size distribution index (-), m is equal to 1−1/n (0 < m < 1) (-) (the parameters
α, n, and m are fitting parameters affecting the shape of the retention curve), l is a poreconnectivity parameter assumed to be 0.5 (-) [27], Ks is the saturated hydraulic conductivity(cm/d), w1 and w2 are weighting factors for the two overlapping regions (-), and α1 , n1 ,m1 (= 1 − 1/n1 ), α2 , n2 and m2 (= 1 − 1/n2 ) are fitting parameters of the sub-curves.
The constant head method was used to measure saturated hydraulic conductivities(Ks ) of dune sand and coarse sand. For gravel, Ks could not be measured accurately using
the same method because of the too high seepage rate. Therefore, it was calculated using
the available empirical model of Hazen [28], following Wongkaew et al. [12]. Based on the
Hazen equation, Ks (cm/s) can be calculated using Equation (5).
2
Ks = CH D10(5)

Water 2021, 13, 2176

5 of 18

where CH is an empirical coefficient (1/(cm/s)), which ranges from 1 to 1000 [29], and D10
is the effective particle size (cm). Following Wongkaew et al. [12], CH was set to 1, and D10
was a geometric mean of ten different D10 values, ranging from 0.475 to 0.485 cm. The basic
5 of 20aresoil hydraulic parameters of the materials used in the HYDRUS (2D/3D) simulations
listed in Table 1.

Water 2021, 13, x

FigureMeasuredandandfitted
fittedretention
retentioncurves
curves andand (b)(b) correspondingofof
Tottori
Dune
Figure2. 2.
(a)(a)
Measuredcorrespondinghydraulichydraulicconductivity
conductivityfunctions
functions
TottoriDunesand, coarse sand, and gravel. The single (SP) and dual (DP) porosity models with n₂ = 3 were used [12].
sand, coarse sand, and gravel. The single (SP) and dual (DP) porosity models with n2 = 3 were used [12].

The constant head method was used to measure saturated hydraulic conductivities

Table 1. Soil hydraulic parameters for the single-porosity (SP) and dual-porosity (DP) models for Tottori Dune sand, coarse(Ks) of dune sand and coarse sand. For gravel, Ks could not be measured accurately using
sand, and gravel [12].
Materials

θr

Tottori Dune sandCoarse sandGravel

0.016
0.004
0.005

the same method because of the too high seepage rate. Therefore, it was calculated using
the available[28], following
Based onntheθsαempirical(1/cm) modeln of Hazen
Ks (cm/s)
w1 Wongkaew
w2 et al.
α2[12].
(1/cm)
2
Hazen equation, Ks (cm/s) can be calculated using Equation (5).
0.384
0.371
0.274

0.030
0.133
0.204

6.66
20.8
15.5

0.020
0.081
𝐾 = 𝐶 𝐷0.99
0.230
0.99

0.01
0.01

0.25
0.25

3(5)
3

where CH is an empirical coefficient (1/(cm/s)), which ranges from 1 to 1000 [29], and D10
is the effective particle size (cm). Following Wongkaew et al. [12], CH was set to 1, and D10
2.2.2.
Watermean
Uptake(RWU)
Model
was aRoot
geometricof tendifferent
D10 values, ranging from 0.475 to 0.485 cm. The basic
soil Inhydraulicparameters
theformmaterials
in theequation
HYDRUSwas(2D/3D)
simulationsare
HYDRUS(2D/3D),oftheof theused
Richards
modifiedand adapted
listedTable 1. flow equation with the consideration of 2D/3D isothermal Darcian flowforthe ingoverningof water in a variably-saturated rigid porous medium by assuming that the air phase does
Table 1. Soil hydraulic parameters
the single-porosity(SP)liquidand dual-porosity(DP)modelsfor Tottori
Dune sand,not play aforsignificant
role in theflow process.
Theequationis describedas follows:coarse sand, and gravel [12].

Materials
Tottori Dune sandCoarse sandGravel"!#∂ 𝑲𝒔 (cm/s) A ∂h𝒘𝟏𝜽𝒓𝒏𝟐
𝜽𝒔
𝜶 (1/cm) ∂θ 𝒏=
K (h) K
−𝟐 S 𝜶𝟐 (1/cm)(6)+ KizA 𝒘∂t 6.66 ∂xi 0.020 ij ∂x j0.016
0.384
0.030
0.004
0.371
0.133
20.8
0.081
0.99
0.01
0.25
3
0.005where 0.274
15.5spatial 0.230
0.01
0.25
3 (cm),xi (i = 1, 2, 0.204
and 3) are thecoordinates0.99
of the Cartesian
coordinate
systemt is time (day), Kij A are components of a dimensionless anisotropy tensor KA , and S is a
2.2.2.
Root(1/d).
WaterAnisotropy
Uptake (RWU)
Model
sink
term
in soilhydraulic properties was not considered in this study.
The
sink
term,
S,
represents
the
water
volume removed
the soil inanda unitvolumeIn HYDRUS (2D/3D), the form of
the Richards
equationfrom
was modified
adapted
per
unit
time
because
of
plant
water
uptakeandisdescribed
byFeddes
et
al.
[23]
as:
for the governing flow equation with the consideration of 2D/3D isothermal Darcian flowof water in a variably-saturated rigid porous medium by assuming that the air phase doesh) =process.
α(h)S pThe equation is described as follows:(7)not play a significant role in the liquidS(flow𝜕𝜕𝜃 function𝜕ℎwhere the water stress response(-)=
𝐾(ℎ)α(h)𝐾 is a prescribed+𝐾− 𝑆 dimensionless function(6)
𝜕𝑡 h 𝜕𝑥
𝜕𝑥 Sp is the potential water uptake rateof the soil water pressure head(0 ≤ α ≤ 1), and(1/d).
water
uptakeis assumed
to beofzero
to saturation(i.e.,
wetterthan
whereNoticexi (i = 1,that
2, and3) arethe spatial
coordinates
the close
Cartesiancoordinate
system(cm),
t is time (day), KijA are components of a dimensionless anisotropy tensor KA, and S is a sinkterm (1/d). Anisotropy in soil hydraulic properties was not considered in this study.
The sink term, S, represents the water volume removed from the soil in a unit volume
per unit time because of plant water uptake and is described by Feddes et al. [23] as:Water 2021, 13, 2176

6 of 18“anaerobiosis point”, h1 ). For h < h4 (the wilting point pressure head), the water uptake isalso assumed to be zero. Water uptake is considered optimal between pressure head h2 andh3 , whereas for pressure head between h3 and h4 (or h1 and h2 ), water uptake decreases (orincreases) linearly with h. The variable Sp is equal to the water uptake rate during periodsof no water stress when α(h) = 1.
In this study, because of the similarity in morphology with Japanese spinach, theFeddes parameters for lettuce were selected from the HYDRUS (2D/3D) database andadapted for the simulations. The initial threshold pressure heads of the lettuce’s stressresponse function and their modified values used in this study are listed in Table 2. Root
distribution parameters for the vertical and horizontal root distributions based on Vrugt’s
root distribution model [30] for Non-CB30, CB15, and CB30 scenarios shown in Table 3
were defined based on root data from Miyake et al. [20] and manually specified in themodel. The HYDRUS (2D/3D) version used in this study did not allow to prescribe timevariable root growth, but only maximum values (e.g., the maximum rooting depth, themaximum rooting radius, etc.). The following Equation (8) describes the three-dimensionalroot distribution function of Vrugt et al. [30] implemented into HYDRUS.




Py
Pz
Px
∗
∗
∗
x
y
z
b( x, y, z) = 1 −
(8)
1−
1−
e−( Xm | x − x|+ Ym |y −y|+ Zm |z −z|)
Xm
Ym
Zm
where Xm , Ym , and Zm are the maximum rooting lengths in the x-, y-, and z- directions (cm),respectively; x, y, and z are distances from the origin of the tree in the x-, y-, and z- directions(cm), respectively; Px (-), Py (-), Pz (-), x ∗ (cm), y∗ (cm), and z∗ (cm) are empirical parameters(x ∗ , y∗ , and z∗ are the coordinates of maximum root uptake intensity); and b( x, y, z) denotes
the three-dimensional spatial distribution of the potential root water uptake (-).
Table 2. Threshold parameters of the Feddes model used in the HYDRUS simulations for spinach.
Threshold Parametersof the Feddes Model

Initial Parametersfor Lettuce (cm)Modified Parameters (cm)h1 (P0+)h2 (P0pt)
h3 H (P2H)
h3 L (P2L)
h4 (P3)

−10
−25
−400
−600
−1200

−10
−25
−50
−50
−100

As displayed in the HYDRUS GUI.

2.2.3. Domain for Simulation
The numerical simulations were conducted in a 3D simulation domain 16 cm deep,
6 cm wide, and 6 cm long (Figure 3). The size of the simulation domain, which was muchsmaller than the container used in the experiments, was determined to match the size ofthe soil profile corresponding to growing one spinach plant. The Non-CB domain consistedof a homogeneous uniform 16-cm layer of dune sand. The CB domain contained a 6-cm
thick CB layer (a 2-cm layer of coarse sand overlying a 4-cm layer of gravel) located at a
depth of 6 cm following the setup of the cultivation experiment (Figures 1 and 3). The
spinach roots were assigned around the center of the simulation domain. The root density
distributions at the cross-section of the center of the 3D model are depicted in Figure 3a.
To compare with the observed VWCs, simulated observation points at the corresponding
depths of 3 cm and 14 cm were assigned at the center of the transport domain.

Table3.
Parameters
for
the
Vrugt
root
distribution
model
[30]
for
different
experimental
scenarios(Non-CB30—A
homoTable
Table3.
3.Parameters
Parametersfor
forthetheVrugt
Vrugtroot
rootdistribution
distributionmodel
model[30][30]forfordifferent
differentexperimental
experimentalscenarios
scenarios(Non-CB30—A(Non-CB30—Ahomohomogeneousreferencecolumnanddaily
irrigation
of
0.30
cm,
CB15—A
column
with
CB
and
daily
irrigation
of
0.15
cm,
and
geneous
geneousreference
referencecolumn
columnandanddailydailyirrigation
irrigationof
of0.300.30cm,cm,CB15—A
CB15—Acolumn
columnwith
withaaaCB
CBand
anddailydailyirrigation
irrigationof
of0.150.15cm,cm,andand
CB30—A
column
with
CB
and
daily
irrigation
of
0.30
cm)
based
on
root
data
from
Miyakeetal.
[20].
CB30—A
CB30—Acolumn
columnwith
withaaaCB
CBand
anddailydailyirrigation
irrigationof
of0.300.30cm)cm)based
basedononrootrootdata
datafrom
fromMiyake
Miyakeetetal.
al.[20].
[20].
Water 2021, 13, 2176

7 of 18

Root
Distribution
Parameters
Root
RootDistribution
DistributionParameters
Parameters

Non-CB30
Non-CB30
Non-CB30

CB15
CB15
CB15

CB30
CB30
CB30

Horizontal
Horizontal
DistribuDistribution—Y
tion—Y

Horizontal
Horizontal
Distribution—X
Distribution—X

VerticalVerticalDistribuDistribution
tion

Maximumrooting
depth
7.2
cm
2.8
cm
4.3
cm
Maximum
depth
7.2
2.8
cm
Maximumrooting
depth
7.2cm
cm experimental scenarios
2.8cm
cm (Non-CB30—A4.3
4.3cm
Table 3. Parameters for
the Vrugtrooting
root distribution
model [30] for different
homogeneous reference column
and
daily
irrigation
of
0.30
cm,
CB15—A
column
with
a
CB
and
daily
irrigation
of
0.15
cm,
Depth
of
maximum
intensity
cm
cm
cm
Depth
222cm
222cm
222cm
Depthofofmaximum
maximumintensity
intensity
cm
cm
cmand
CB30—A column with a CB and daily irrigation of 0.30 cm) based on root data from Miyake et al. [20].
Parameterpz
Parameter
111
111
111
Parameterpz
pz
Root Distribution Parameters
Non-CB30
CB15
CB30
Maximumrootingradiuscm
cm
cm
Maximum
222cm
222cm
222cm
Maximumrooting
rootingradius
radiuscm
cm
cm
Maximum rooting depth
7.2 cm
2.8 cm
4.3 cmDepthof maximum
intensitycm
Vertical Distribution Radius
Radiusofmaximum
intensity
cm
cm
cm
of
intensity
0200cm0200cm0020cmRadiusofmaximum
maximum
intensitycm
cm
Parameter pz
1
1
1
Parameterpx
Parameter
111
111
111
Parameterpx
px
Maximum rooting radius
2 cm
2 cm
2 cm
Radiusmaximum intensity
0 cm
0 0cm
0 00cmCentercoordinate
00
Centercoordinate
000
Centerof
coordinate
0
Horizontal
Parameter px
1
1
1
Distribution—X
Maximumrootingradiuscm
cm
Maximumrootingradius
222cm
222cm
222cm
Maximumrootingradiuscm
cm
Centercoordinate
0cm
0cm
0
Horizontal
Distribution—Y

MaximumrootingradiusRadiusofmaximum
intensity
Radiusintensity
Radiusof
ofmaximum
maximum
intensity
Radius of maximum intensity
ParameterpyParameterpyParameterpy pyParameter

0200cm
cm
0 cm
0000

0200 cm
cm
0 cm
0000

cm
cm
0020cm
cm
0 cm
0000

Graphical
distribution
of
the
root
systemGraphical
Graphicaldistribution
distributionof
ofthe
therootrootsystemsystemWater 2021, 13, x

8 of 202.2.3.
DomainforSimulation
2.2.3.
2.2.3.Domain
Domainfor
forSimulation
SimulationThe
numerical
simulations
wereconductedin
3D
simulation
domain
16
cm
deep,
The
Thenumericalnumericalsimulations
simulationswere
wereconductedconductedin
inaaa3D
3Dsimulation
simulationdomain
domain1616cm
cmdeep,deep,666
cm
wide,
and
cm
long
(Figure
3).
The
size
of
the
simulation
domain,
which
was
much
cm
cm wide,wide, andand 666 cm
cm longlong (Figure(Figure 3).
3). TheThe size
size ofof thethe simulation
simulation domain,
domain, which
which was
was muchmuchsmallerthanthecontainer
used
in
the
experiments,
wasdetermined
to
match
the
size
of
smallersmaller than
than thethe container
container used
used inin thethe experiments,
experiments, waswas determined
determined to
to matchmatch thethe size
size ofof
thesoil
profilecorresponding
togrowingone
spinachplant.
The
Non-CB
domainconthethe soil
soil profile
profile corresponding
corresponding toto growing
growing one
one spinach
spinach plant.
plant. The
The Non-CB
Non-CB domain
domain conconsistedofhomogeneous
uniform
16-cm
layer
of
dune
sand.
The
CB
domaincontained
sistedsisted ofof aaa homogeneous
homogeneous uniform
uniform 16-cm
16-cm layer
layer ofof dunedune sand.
sand. TheThe CBCB domain
domain contained
contained aaa
6-cm
thick
CB
layer
(a
2-cm
layer
of
coarsesandoverlying
4-cm
layer
of
gravel)located
6-cm
6-cmthickthickCBCB layerlayer(a(a2-cm2-cmlayer
layerofof coarsecoarsesandsandoverlying
overlyingaaa4-cm
4-cmlayer
layerofofgravel)gravel)locatedlocated
at
depth
of
cm
following
the
setup
of
the
cultivation
experiment(Figures
and
3).
The
at
ataaadepth
depthofof666cmcmfollowing
followingthe
thesetup
setupofofthe
thecultivation
cultivationexperimentexperiment(Figures(Figures111andand3).
3).TheThe
spinachroots
wereassignedaroundthecenterofthesimulation
domain.
Therootdensityspinachspinachroots
rootswere
wereassigned
assignedaroundaroundthe
thecenter
centerof
ofthe
thesimulation
simulationdomain.
domain.The
Theroot
rootdensity
density
distributionsatthecross-sectionofthecenter
of
the
3D
model
are
depictedinFigure3a.
distributions
distributions at
at thethe cross-section
cross-section ofof thethe center
center ofof thethe 3D
3D modelmodel areare depicted
depicted inin Figure
Figure 3a.
3a.
To
comparewiththeobserved
VWCs,
simulated
observation
pointsatthecorrespondingTo
Tocomparecomparewith
withthetheobserved
observedVWCs,
VWCs,simulated
simulatedobservation
observationpoints
pointsat
atthe
thecorresponding
corresponding
depths
of
cm
and
14
cm
were
assignedatthecenterofthetransport
domain.
depthsdepthsofof333cm
cmand
and14
14cm
cmwerewereassigned
assignedat
atthe
thecenter
centerof
ofthe
thetransporttransportdomain.
domain.

Figure 3.
3. Schematic
Schematicof
ofrootroot density
density distribution
distributionat
at thethe center
center cross-section
cross-section ofof thethe simulation
simulation domain
domain (a),(a), andand irrigation
irrigation
Figureapplicationtypes
types(uniform(uniform (b(b11),),line-sourceline-source(b(b22),),and
andplant-targeted
plant-targeted(b(b3)).
3 )).
application

2.2.4.Numerical
NumericalAnalysis
Analysisof
ofCultivation
CultivationExperiments
Experiments
2.2.4.
In the numericalthethecultivation
experiments,
the observed
irrigation
fluxesInnumericalanalysis
analysisofof
cultivation
experiments,
the observed
irrigation
were
used
as
the
surfaceboundary
condition.
Since
the
experiments
were
conducted
influxes were used as the surface boundary condition. Since the experiments were cona
controlledenvironment(agreenhouse),
there
was
no
rainfall.
The
actualevaporation
ducted in a controlled environment (a greenhouse), there was no rainfall. The actual evaporation (E) and transpiration (T) rates (Figure 4) were adopted from Wongkaew et al. [12]by dividing the actual evapotranspiration (ET) data acquired from Miyake et al. [20] using
the Surface Cover Fraction (SCF), which was measured once a week. Alternative irrigations were evaluated by assessing if these actual E and T rates were achieved when usedas potential rates in the model simulations. The model also requires the minimum allowedWater 2021, 13, 2176

Water 2021, 13, x

8 of 18(E) and transpiration (T) rates (Figure 4) were adopted from Wongkaew et al. [12] bydividing the actual evapotranspiration (ET) data acquired from Miyake et al. [20] using theSurface Cover Fraction (SCF), which was measured once a week. Alternative irrigations
were evaluated by assessing if these actual E and T rates were achieved when used as
potential rates in the model simulations. The model also requires the minimum allowed
pressure head value at the soil surface (hCritA), which defines the end of the first stage of
evaporation. For numerical reasons, according to HYDRUS help, this value needs to be
defined to make the corresponding water content at least 0.005 bigger than the residual
water content, which is necessary particularly for coarse-textured soils, because of thesteep slopes observed in their retention curves. A big value of hCritA may therefore be
required for the stability of the numerical solution for the soils with coarse grain size. Also,
when root water uptake is considered, a value of hCritA smaller (when negative) than P3(i.e., the wilting point pressure head, h4 ) should be taken. In this study, hCritA was set
to −250 cm. For the boundary condition of the sidewalls and the bottom of the domain,a no-flux was assigned. No flow boundary at the sidewalls of the simulation domain
was adapted, assuming that no lateral water flow interaction between plants occurred
9 of 20 as
uniform irrigation was applied in the experiments.

Figure4.4.Evaporation
Evaporationand
andtranspirationtranspirationrates
ratesforforNon-CB30
Non-CB30(top),(top),CB15
CB15(middle),(middle),and
andCB30
CB30(bottom)(botFiguretom) adopted from Wongkaew et al. [12] by dividing the actual evapotranspiration rate obtained
adopted from Wongkaew et al. [12] by dividing the actual evapotranspiration rate obtained fromfrom Miyake et al. [20] for each simulation.
Miyake et al. [20] for each simulation.

As an initial condition for the pressure head, a constant value of -200 cm was alloAs an initial condition for the pressure head, a constant value of −200 cm was allocatedcated to the entire domain. The model simulations were conducted from seeding until theto the entire domain. The model simulations were conducted from seeding until the
final cultivation day, i.e., a total of 28 days. The initial time step was adjusted to a tiny
final cultivation
day, i.e., a total of 28 days. The initial time step was adjusted to a tiny
value of 10−9 days to mitigate the effects of highly nonlinear soil hydraulic properties ofvalue of 10−9 days to mitigate the effects of highly nonlinear soil hydraulic propertiesofthe coarse-textured soils. The minimum and maximum allowed time steps were set to 10−11
and 10−3 days, respectively.
The simulated time series of VWCs at two observation nodes (3 cm and 14 cm deep)
were compared with the observed VWC values at corresponding locations to assess themodel’s accuracy. The squared Pearson correlation coefficient (R2) and the root meansquared error (RMSE) were computed [31]. The R2 value, varying between 0 (no correla-Water 2021, 13, 2176

9 of 18

the coarse-textured soils. The minimum and maximum allowed time steps were set to10−11 and 10−3 days, respectively.
The simulated time series of VWCs at two observation nodes (3 cm and 14 cm deep)
were compared with the observed VWC values at corresponding locations to assess themodel’s accuracy. The squared Pearson correlation coefficient (R2 ) and the root meansquared error (RMSE) were computed [31]. The R2 value, varying between 0 (no correlation)and 1 (perfect correlation), defines how the observed and predicted values share thevariance. The RMSE determines how well the model performs by assessing the difference
between observed and predicted values. The RMSE is greater than or equal to 0, though
the value of 0 (a perfect fit) is rarely attained in real life. Generally, the smaller the
RMSE, the better the performance. The R2 and RMSE values were calculated according toEquations (9) and (10), respectively.
2∑i Oi − O Pi − P
R =
2
2
∑i Oi − O ∑i Pi − P
2"RMSE =

1 N( Pi − Oi )2
N i∑
=1

(9)

#0.5
(10)

where Pi is the ith predicted VWC value, Oi is the ith observed VWC value, P is the mean
predicted VWC value, O is the mean observed VWC value, and N is the total number ofobservations. The R2 is unitless, while the RMSE units are the same as the units of P and O.
2.3. Analysis of Alternative Hypothetical Irrigation Scenarios
The effects of the CB system on soil water dynamics and root water uptake during
spinach cultivation were evaluated numerically under three different types of irrigation;
uniform, line-source, and plant-targeted (referred to as UNI, LSI, and PTI, respectively).
For each type of irrigation, three different irrigation schedules (referred to as Irri0, Irri1, andIrri2) were simulated to evaluate the influence of further reducing the amount of irrigationwater on soil water dynamics. All simulations were conducted using the same simulation
domain in 3D as done in the analysis of the cultivation experiments. In all simulations, afree bottom drainage boundary condition was used to allow water to drain freely at thebottom boundary and to quantify water loss to deep drainage.
Description of the three irrigation schedules mentioned above are as below:•••

Irri0: Daily irrigation of 0.3 cm of water from Day 4 to the last day of the cultivation period;
Irri1: Bi-daily irrigation of 0.3 cm of water during the first half of plant development stage from Day 4 to Day 14, and then daily irrigation until the last day of thecultivation period;
Irri2: Bi-daily irrigation of 0.3 cm of water during the entire cultivation period fromDay 4 to Day 28.

In all hypothetical irrigation scenarios, similarly as in the experiments, 0.5 cm of water
was applied on Day 1 and Day 2, and no irrigation was applied on day 3.
In UNI, irrigation was applied uniformly on the whole surface (Figure 3b1 ). In LSI,
the same amount of water as in the UNI was applied around the plant on a surface areaof 3 cm × 6 cm (Figure 3b2 ). In PTI, the same amount of water was applied around theplant on a surface area of 3 cm × 3 cm (Figure 3b3 ). The results of all these scenarios
are discussed based on water balances and VWCs in the following sections. Note that
UNI with Irri0 corresponds to either CB30 or Non-CB30 in the experiment, dependingon whether a CB layer is installed or not. The only differences are the bottom boundary
conditions; a free-drainage BC for hypothetical scenario simulations and a no-flux BC forexperimental simulations. A summary of all simulation setups with their descriptions inboth cultivation experiments and hypothetical irrigation scenarios is given in Table 4.

Water 2021, 13, 2176

10 of 18

Table 4. Summary of all simulation setups.
Experiments/Scenarios

Cultivation
Experiments

Bottom Boundary
Condition

Irrigation

Simulation
Non-CB30

No-flux

Uniform

CB15
CB30
Non-CB_Irri0

Hypothetical Irrigation
Scenarios

Free drainage

Uniform,line-source, andplant-targeted

CB_Irri0
CB_Irri1

CB_Irri2

Description
A homogeneous column (without CB)
and daily irrigation of 0.30 cm
With CB and daily irrigation of 0.15 cm
With CB and daily irrigation of 0.30 cm
A homogeneous column (without CB)
and daily irrigation of 0.30 cm
With CB and daily irrigation of 0.30 cm
With CB and bi-daily irrigation of
0.30 cm from Day 4 to Day 14, andthen daily irrigation of 0.30 cm until
the last day
With CB and bi-daily irrigation of
0.30 cm from Day 4 to Day 28

Since the Non-CB30 and CB30 cases were used as the benchmarks, the root distribution
parameters and evapotranspiration rates obtained from the Non-CB30 experiment wereused in the hypothetical scenario Non-CB_Irri0, while those from the CB30 experimentwere used in hypothetical scenarios CB_Irri0, CB_Irri1, and CB_Irri2.
3. Results and Discussions
3.1. Simulation of Cultivation Experiments
3.1.1. Volumetric Soil Water Contents
Figure 5 shows simulated VWC profiles at the center of the 3D simulation domain
for Non-CB30, CB15, and CB30, every 4 days between the start of the seeding period andthe end of the experiment (Day 28). Observed VWCs at two depths of 3 cm and 14 cm arealso plotted for comparison. The results indicated that, up to Day 4, the simulated VWCs
in all experimental setups matched the observed data well at both observation depths.
For the Non-CB30 setup, while the model predicted VWCs well at the lower observation
point throughout the cultivation period, the model over-estimated VWCs at the upper
observation point after Day 8. VWCs at the lower observation point were much higher than
VWCs at the upper observation point, indicating that water accumulated at the bottom of
the container in the Non-CB30 setup. In the CB15 and CB30 setups, the simulated VWCs at
the lower observation points remained nearly constant at a relatively small value during
the entire simulation period, while the observed VWCs slightly increased until Day 24 for
CB15 and notably increased for CB30 until the last day. In the upper observation points,
VWCs were generally overestimated for both CB15 and CB30.
Figure 6 shows scatter plots of simulated VWCs against observed VWCs at bothobservation points. For the Non-CB30 setup, simulations resulted in a small correlation(R2 = 0.044) and a relatively high RMSE (=0.036) at the upper observation point. In contrast,
simulations achieved a high correlation (R2 = 0.908) and a small error (RMSE = 0.017) forthe lower observation point for Non-CB30. For the CB15 system, the results showed a goodcorrelation (R2 = 0.885) with a reasonable error (RMSE = 0.027) at a depth of 3 cm and avery small error (RMSE = 0.013) at a depth of 14 cm with almost no correlation. The resultsfor CB30 showed a large error (RMSE = 0.063) and a reasonable correlation (R2 = 0.607) at a
depth of 3 cm and relatively large errors (RMSE = 0.048) with almost no correlation at a
depth of 14 cm, the depth which was beneath the CB layer.

Water 2021, 13, x

12 of 20

Water 2021, 13, 2176

11 of 18

Figure 5. Simulated and observed volumetric soil water contents at different times (0, 4, 8, 12, 16, 20, 24, and 28 d; from
Figure 5. Simulated and observed volumetric soil water contents at different times (0, 4, 8, 12, 16, 20, 24, and 28 d; fromleft toto right)
right) andand forfor different
different experimental
homogeneous reference
reference columnofleftexperimental setups
setups (Non-CB30—A(Non-CB30—A homogeneouscolumn andand daily
daily irrigation
irrigation of
0.30
cm(top),
CB15—A
column
with
a
CB
and
daily
irrigation
of
0.15
cm
(middle),and
CB30—A
column
with
a
CB
and
0.30 cm (top), CB15—A column with a CB and daily irrigation of 0.15 cm (middle), and CB30—A column with a CB anddaily irrigation
irrigation ofdaily
of 0.30
0.30 cmcm (bottom)).
(bottom)).

Althoughsimulated
results
that infiltratingwater
stopped VWCs
at the CB
Figure6 shows
scatterplotsindicatedof simulated
VWCs against
observed
at interface
both ob(the
VWCs
were
high
in
the
top
layer
and
almostzero
inthebottomlayer),
theobserved
servation points. For the Non-CB30 setup, simulations resulted in a small correlation
(R2
VWCs
at
the
bottomofthecontainers
increasedslowlyduringtheexperiment,
especially= 0.044) and a relatively high RMSE (=0.036) at the upper observation point. In contrast,in CB30. Thismay have
beencorrelation
due to wateralong
the containers’
sidewalls
simulations
achieved
a high(R2 leaching= 0.908) anda smallerror (RMSE= 0.017)into
forthe bottom layer, which usually occurs in lysimeter studies [32–34]. Moreover, vapor flow
the lower observation point for Non-CB30. For the CB15 system, the results showed a
towards the bottom containers may have happened because of the downward temperature
good correlation (R2 = 0.885) with a reasonable error (RMSE = 0.027) at a depth of 3 cm andgradient [20,35]. The soil surface temperature increased during the daytime while thea very small error (RMSE = 0.013) at a depth of 14 cm with almost no correlation. The
sidewalls and the bottom boundary of all containers were insulated with styrene foams. In
results for CB30 showed a large error (RMSE = 0.063) and a reasonable correlation (R2 =addition, although visual inspection could not detect obvious leakage of water through the
0.607) at a depth of 3 cm and relatively large errors (RMSE = 0.048) with almost no correCB layer, there could be liquid water flow through the CB layer in the experiments.
lation at a depth of 14 cm, the depth which was beneath the CB layer.
3.1.2. Cumulative Fluxes
Simulated cumulative fluxes (i.e., infiltration, evaporation, transpiration) for NonCB30, CB15, and CB30 are summarized in Table 5. The Non-CB30 and CB30 setups
produced similar water balance components: more than 50% of applied water was taken
up by plant roots, almost 35% was lost via evaporation, and less than 15% was stored in
the soil. For CB15, most of the irrigation water was lost due to evaporation (about 95%),
while plant roots absorbed only less than 3%, and about 2% remained in the soil. These
results imply that the CB system may enhance water losses through evaporation. When asmall amount of irrigation water is applied, most of it can be returned to the atmosphereby evaporation.

Water 2021, 13, x

13 of 20

Water 2021, 13, 2176

12 of 18

Figure6.Scatter
Scatterplots
betweensimulatedand observed
volumetricsoilcontents
water contentsof 314(top)andFigure 6.
plots
betweensimulatedand observed
volumetricsoil waterat depthsatof depths
3 (top) and(bottom)14
cm fordifferent experimental
setups (Non-CB30—A
referencecolumnand ofdaily
cm(bottom)
for different
experimental
setups (Non-CB30—A
homogeneous homogeneousreference column
and daily
irrigation
0.30irrigation
cm (left),of0.30 cmcolumn(left), CB15—A
columnwithirrigationa CB andof
daily
of 0.15cmCB30—A(middle),and CB30—A
columnwith irrigationa CB and
CB15—A
with a CB
and daily
0.15irrigation
cm (middle),
andcolumnwith a CB
and daily
daily
of 0.30 cm (right)).
of 0.30irrigation
cm (right)).

simulatedindicated
thatevaporation,
infiltrating root
water
stoppedthe CBand
intervalues of results
cumulative
infiltration,water
uptakeat(RWU),
soil
TableAlthough
5. Simulated
face
(the
VWCs
were
high
in
the
top
layer
and
almostzero
inthebottomlayer),
theobwater storage for different experimental scenarios (Non-CB30—A homogeneous reference column
servedVWCs
at thebottomof the containers
increasedslowlytheand daily
irrigationof 0.30
cm, CB15—A
column witha CB anddaily during
irrigationof experiment,
0.15 cm, andespeciallyin CB30.
been dueofto0.30
waterCB30—A column
withThisa CBmay
and have
daily irrigationcm).leaching along the containers’ sidewalls into the bottom layer, which usually occurs in lysimeter studies [32–34]. Moreover,
Cumulative
Fluxes
Non-CB30may(cm) have happened
CB15 (cm)CB30(cm)
vapor flow
towards
the bottom containers
because of thedownward
Infiltration(Irrigation)
8.50daytime
temperature
gradient[20,35]. The soil 8.50
surface temperature4.75
increased during the
Evaporation
2.78
4.53
2.97
while the sidewalls and the bottom boundary of all containers were insulated with styrene
RootInwater
uptakealthough(RWU) visual inspection
4.75
0.13 obvious leakage
4.47of water
foams.
addition,
could not detect
Soil water storage
0.97
0.10
1.06
through the CB layer, there could be liquid water flow through the CB layer in the experiments.
To further illustrate the relationship between plant growth and simulated RWU,
spinachheights (Figure
3.1.2.
Cumulative
Fluxes7a) and spinach’s leaf area index (LAI) (Figure 7b) for all experimental
setups(Non-CB30,
CB15, (i.e.,
and CB30)
at theevaporation,
2nd, 3rd, andtranspiration)
4th week of the
Simulated cumulative fluxes
infiltration,
for plant
Nongrowth
obtainedfrom
Miyake
et
al.
[20]
are
plottedagainstcorresponding
simulated
cumuCB30, CB15, and CB30 are summarized in Table 5. The Non-CB30 and CB30 setups prolative RWU
at correspondingdays. The graphs
show
goodofagreement
betweenduced
similarwater balance components:
morethana50%
applied water
wassimulated
taken uproot water uptake and the corresponding spinach growth, having high positive correlations
by plant roots, almost 35% was lost via evaporation,2 and less than 15% was stored in the
between these two pairs of the variables, with an R value of 0.93 for simulated RWU vs.
soil. For CB15, most of theirrigation water was lost due to evaporation (about 95%), whilespinach height and an R2 value of 0.99 for simulated RWU vs. LAI. For instance, simulated
plant roots absorbed only less than 3%, and about 2% remained in the soil. These results
RWU in Non-CB30 was similar to that in CB30, while the actual spinach growth conditions
imply that the CB system may enhance water losses through evaporation. When a small
in the two experiments were nearly the same. A more obvious example is that simulated
amount of irrigation water is applied, most of it can be returned to the atmosphere by
RWU in CB15 was only about 3% of that in CB30 or Non-CB30 (Table 5), while leaf area
evaporation.
index (LAI) of the spinach in CB15 was no more than 3% of those in CB30 or Non-CB30(Figure 7b). These results indicated that the evaluation of spinach growth based on thesimulated RWU is dependable. Moreover, it was also observed that when the spinachWater 2021, 13, 2176

ach growth conditions in the two experiments were nearly the same. A more obvious example is that simulated RWU in CB15 was only about 3% of that in CB30 or Non-CB30(Table 5), while leaf area index (LAI) of the spinach in CB15 was no more than 3% of thosein CB30 or Non-CB30 (Figure 7b). These results indicated that the evaluation of spinach
growth based on the simulated RWU is dependable. Moreover, it was also observed that
13 of 18
when the spinach developed more leaves in CB30 and Non-CB30 and the surface cover
fraction correspondingly increased, simulated evaporation rates in those two experimental simulations were low, meaning that water loss is mainly due to transpiration. On
developed more leaves in CB30 and Non-CB30 and the surface cover fraction correspondthe other hand, when the spinach grew extremely small in CB15, simulated evaporation
ingly increased, simulated evaporation rates in those two experimental simulations werein that treatment was quite large.
low, meaning that water loss is mainly due to transpiration. On the other hand, when
As simulated cumulative RWU reflects the spinach growth well, simulated RWU is
the spinach grew extremely small in CB15, simulated evaporation in that treatment wasused in the remainder of the manuscript to discuss the performance of various irrigation
quite large.
scenarios.

FigureScatterplot
cumulative
RWU
simulatedmodel
versusspinachheightsindex
obFigure
7. 7.Scatter
plot
ofof
cumulative
RWU
simulatedbyby
thethemodel
versus(a) (a)spinachheights
andand(b) (b)leafleaf
areaareaindex
obtainedtainedfrom
Miyake
et
al.
[20]
at
days
14,
21,
and
28
for
all
experimental
setups(Non-CB30,
CB15,
and
CB30).
from Miyake et al. [20] at days 14, 21, and 28 for all experimental setups (Non-CB30, CB15, and CB30).

As simulated cumulative RWU reflects the spinach growth well, simulated RWU
is used in the remainder of the manuscript to discuss the performance of various irrigation scenarios.
3.2. Simulation Results for Hypothetical Scenarios
3.2.1. Water Balances for Different Irrigation Scenarios
Figure 8 shows cumulative water fluxes for different irrigation scenarios simulatedby the HYDRUS (2D/3D) model in a 3D domain. Simulated water balance components
were found to vary depending on irrigation types (whether uniform, line-source, or planttargeted), especially in the soil profiles with a CB. For Non-CB with Irri0 (Figure 8a),the three different irrigation types, UNI, LSI, and PTI, did not influence much simulated
cumulative fluxes and produced almost no change in soil water storage (SWS). Whilecumulative evaporation was reduced from 1.3 cm for UNI to 0.8 cm for LSI, in which
irrigation water was applied on a smaller part of the surface, both cumulative root water
uptake (RWU) and cumulative bottom drainage (BD) slightly increased. When the sameamount of water was directly applied around the plant in PTI, cumulative evaporation was
further reduced to 0.5 cm with slight increases in RWU and BD.
On the other hand, simulated cumulative fluxes in soil profiles with a CB were strongly
affected not only by irrigation types but also by irrigation schedules (Figure 8b–d). For
instance, for CB_Irri0 with UNI, in which daily irrigation of 0.3 cm was applied uniformly,
cumulative evaporation, cumulative RWU, and cumulative BD were 3.0 cm, 4.5 cm, and0 cm, respectively, and SWS increased by 1.1 cm (Figure 8b). For LSI, while cumulative
evaporation decreased to 1.5 cm and cumulative RWU dropped to 1.2 cm, cumulative BD
and SWS increased considerably to 3.5 cm and 2.3 cm, respectively. The water balancecomponents varied further for PTI, except for SWS. Cumulative evaporation for PTI was
reduced by half compared to that for LSI. Additionally, while cumulative RWU decreased

Water 2021, 13, x

Water 2021, 13, 2176

15 of 20

3.2. Simulation Results for Hypothetical Scenarios

14 of 18

3.2.1. Water Balances for Different Irrigation Scenarios
Figure 8 shows cumulative water fluxes for different irrigation scenarios simulatedtocumulative
BD increased
4.4 domain.
cm for PTI.
Althoughirrigating
cm of waterby1.0
thecm,HYDRUS(2D/3D)model in to
a 3D
Simulatedwaterbalance0.3
components
everyday (CB_Irri0)
was good
UNI astypes
the waterbalanceratioline-source,
was improved
withwere foundto vary depending
onfor
irrigation(whetheruniform,or planthigher
cumulative
RWU
andsmallerwater
lossescomparedtoNon-CB,
it
might
not
targeted), especially in the soil profiles with a CB. For Non-CB with Irri0 (Figure 8a), thebe
appropriate
for irrigation
LSI and PTI,
as aUNI,
hugeLSI,amountof waterwasinfluence
lost throughdrainage
three different
types,and PTI,
did not
muchbottom
simulated
cuwithlower
cumulative
RWUs.
The
local
applicationsof
thisamount
of
water
could
lead tomulative fluxes and produced almost no change in soil water storage (SWS). While cumuan
exceedanceof thethresholdthe wetof LSI,the Feddes
transpiration
lativeevaporation
wasreducedpressure
from 1.3heads
cm forinUNI
to 0.8range
cm forin which
irrigation
function.
This
indicates
that
even
with
a
CB
systeminstalled,applying
irrigationa smallwater was applied on a smaller part of the surface, both cumulative root wateronuptake
area
with
an
excessivewater
amountcould
reduceroot
water
uptakebecause
the
root
zone
(RWU) and cumulative bottom drainage (BD) slightly increased. When the same amountisoftoowetwasanddirectlyroots areunder
stress.theRisingat theinterface may
producewaterappliedaroundplantwaterin PTI,content
cumulative
evaporation
was further
breakageof 0.5
thecm
CBwithbarrierandincreases
water losses
throughbottom drainage.
reduced toslightin RWUand BD.

Figure8.8. Cumulative
Cumulative water
water fluxes
(UNI,
LSI,
Figurefluxes forfor uniform,
uniform,line-source,line-source,and
andplant-targeted
plant-targetedirrigation
irrigation
(UNI,
LSI,
and PTI, respectively) for different simulation setups: (a) Non-CB_Irri0, (b) CB_Irri0, (c) CB_Irri1,
and PTI, respectively) for different simulation setups: (a) Non-CB_Irri0, (b) CB_Irri0, (c) CB_Irri1,
and (d) CB_Irri2. Free bottom drainage was used fot the bottom boundary condition. Note: Irri—and(d) CB_Irri2. Free bottom drainage was used fot the bottom boundary condition. Note: Irri—irrigation; Eva—evaporation; RWU—root water uptake; BD—bottom drainage; and SWS—soil wairrigation;
ter storage.Eva—evaporation; RWU—root water uptake; BD—bottom drainage; and SWS—soilwater storage.

On the other hand, simulated cumulative fluxes in soil profiles with a CB were
For CB_Irri1
irrigation
water
was
applying(Figure
water every
strongly
affected (Figure
not only8c),
by whenirrigation
types butalso
by reduced
irrigationby
schedules
8b–
two
days
during
the
first
half
of
the
plant
growthstage,
cumulative
bottomdrainage
didd). For instance, for CB_Irri0 with UNI, in which daily irrigation of 0.3 cm was applied
not
occur
for
UNI.
Cumulative
evaporation
andcumulative
RWU
slightly
decreased
uniformly, cumulative evaporation, cumulative RWU, and cumulative BD were 3.0 cm,to2.7
cm, respectively,
respectively, and
while
SWS
was remarkably
to 0.1cmLSI,
compared
4.5cm
cm,andand4.10 cm,
SWS
increasedby 1.1 cmreduced(Figure 8b).
For
while
to
CB_Irri0 with
UNI. Notable
resultswere
when it RWU
came dropped
to LSI withcumulative
evaporation
decreased
to 1.5
cmobtained
and cumulativeto CB_Irri1.
1.2 cm,Cumulative
furtherconsiderably
reduced to as
as 1.5
cumulative RWU
cumulative evaporation
BD and SWSwas
increased
to low
3.5 cmandcm,
2.3while
cm, respectively.
The
and SWS increased to 4.5 cm and 1.0 cm, respectively, similar to the values obtainedfor UNI with CB_Irri0. However, for PTI, although cumulative evaporation was further
reduced by half compared to cumulative evaporation for LSI and the SWS increased to
2.4 cm, cumulative RWU was considerably reduced to 2.6 cm, and cumulative BD increased

Water 2021, 13, 2176

15 of 18

to 1.3 cm. Like CB_Irri0 with PTI, applying irrigation over a smaller area resulted inmalfunctioning of a CB, as excess infiltrating water accumulated at the CB interface.
For CB_Irri2 with UNI (Figure 8d), in which irrigation water was further reduced byapplying water every two days throughout the cultivation period, cumulative evaporationand cumulative RWU decreased to 2.5 cm and 2.3 cm, respectively. Cumulative BD and
SWS were the same as for CB_Irri1 with the same irrigation type. For LSI with CB_Irri2,
although there was no effect on cumulative BD and SWS, cumulative evaporation was
reduced to 1.5 cm while cumulative RWU increased to as much as 3.3 cm. For PTI, on
the other hand, even though the total amount of water applied was remarkably reduced
compared to Irri0 and Irri1, cumulative RWU was attained at 4.0 cm, while cumulative
evaporation decreased to 0.7 cm, and cumulative BD and SWS remained unchanged.
This result indicates that because the total amount of water applied was much smaller,even though the water was applied locally around the plant, there was no excess water
accumulated at the CB interface, and thus the CB functioned as designed.
Interestingly, although in CB_Irri2 with UNI, the same total amount and the same
type of irrigation were applied as in CB15, cumulative RWU was much larger than that in
CB15. This may be explained by the fact that applying 0.30 cm of water once in two days in
CB_Irri2 resulted in the reduction of water losses by evaporation as water could infiltratemuch deeper. When applying 0.15 cm of water every day, water loss by evaporation
was enhanced.
This study demonstrated that in addition to drastically reducing water loss throughevaporation, applying irrigation water around the plant, as done in LSI and PTI, could improve soil water storage capacity compared to UNI. However, when the water application
area was reduced, the soil at the interface of the CB layer can be quickly saturated withexcess infiltration water to cause malfunctioning of the CB system. Increased saturation in
the soil may further inhibit RWU by reducing oxygen in the root zone. The plant growth
can be thus depressed. The numerical analysis of hypothetical irrigation scenarios suggeststhat the total amount of water needs to be reduced if water is applied more locally when
the CB system is installed to avoid its malfunctioning and plant growth depression.
3.2.2. Vertical VWCs by Different Scenarios
To further understand how irrigation scenarios affect soil water dynamics, temporalchanges in simulated VWC profiles in the center of the simulation domain for all irrigation
scenarios were plotted in Figure 9. VWCs were almost identical between different irrigation
types for Non-CB_Irri0 as differences in irrigation application areas (either 6 cm × 6 cm
for UNI, 3 cm × 6 cm for LSI, or 3 cm × 3 cm for PTI) produced only minimal differencesin temporal changes in the VWC profiles (they were mostly overlapping as shown in thisfigure). The same trend was observed for cumulative fluxes (Figure 8a). On the contrary,
for the soil profiles having a CB, the VWC profiles for the three irrigation types, UNI, LSI,
and PTI, were visibly different from Day 4, regardless of the irrigation schedule (i.e., Irri0,
Irri1, or Irri2). Figure 9 shows that the VWC at the top of the CB layer for PTI with CB_Irri0
reaches a level that would cause water to percolate through the CB system (to break thefunction of the CB system and to increase downward water percolation) after Day 8. The
CB function was broken for LSI with CB_Irri0 from Day 12, while it was never broken for
UNI with this irrigation schedule. Therefore, no bottom drainage for UNI was exhibited in
Figure 8b.
The VWC profiles were all different for the three irrigation types with CB_Irri1. Only
PTI with CB_Irri1 had water percolation through the CB layer after Day 20. Consequently,
bottom drainage was only observed for PTI in this irrigation schedule, as presented in
Figure 8c. The CB layer remained fully functioning regardless of the irrigation type with
CB_Irri2. PTI provided the highest VWC in the root zone (i.e., the soil layer above the
CB layer) until Day 24. At the end of the cultivation period (i.e., Day 28), the entire
domain became completely dry for all irrigation types with Irri2. As a result, there was no

Water 2021, 13, 2176

figure). The same trend was observed for cumulative fluxes (Figure 8a). On the contrary,
for the soil profiles having a CB, the VWC profiles for the three irrigation types, UNI, LSI,
and PTI, were visibly different from Day 4, regardless of the irrigation schedule (i.e., Irri0,
Irri1, or Irri2). Figure 9 shows that the VWC at the top of the CB layer for PTI with CB_Irri0
reaches a level that would cause water to percolate through the CB system (to break the
18
function of the CB system and to increase downward water percolation) after Day 16
8. of
The
CB function was broken for LSI with CB_Irri0 from Day 12, while it was never broken for
UNI with this irrigation schedule. Therefore, no bottom drainage for UNI was exhibitedbottom drainage throughout the simulation period and a relatively small soil water storage
in Figure 8b.
regardless of the irrigation type (Figure 8d).

Figure9. Profiles
Profilesofofvolumetricvolumetricsoil
soilwater
water
contents
centerof thedomains
obtained
the model
for uniform,
Figure 9.
contents
at at
thethecenterof thedomains
obtained
fromfrom
the model
for uniform,
lineline-source,
plant-targeted
irrigation
(UNI,
LSI,
and
PTI,
respectively)atdifferent
times
(0,
4,
8,
12,
16,
20,
24
and
28 day)source, plant-targeted irrigation (UNI, LSI, and PTI, respectively) at different times (0, 4, 8, 12, 16, 20, 24 and 28
day)
and
and
for different
simulation
setupsbottomdrainage(Non-CB_Irri0,
CB_Irri0,
CB_Irri1,
CB_Irri2).
for different
simulation
setups
withwith
free freebottomdrainage(Non-CB_Irri0,
CB_Irri0,
CB_Irri1,
andandCB_Irri2).

4. Conclusions
The HYDRUS (2D/3D) model was employed to evaluate soil water dynamics insoil profiles with a capillary barrier system installed during a cultivation experiment ofJapanese spinach under different hypothetical irrigation scenarios. The performance of the
HYDRUS (2D/3D) model was first evaluated by simulating the soil water dynamics in a 3D
domain based on the cultivation experiment of Miyake et al. [20]. The results showed goodcorrespondence between observed and simulated VWCs and confirmed the effectiveness of
the CB layer to improve root zone conditions by suppressing downward water percolation.
Based on the numerical analysis of hypothetical irrigation scenarios, the following
conclusions are drawn: (1) while a CB generally reduced water losses due to deep percolation and improved root zone water storage capacity during the cultivation, water losses
due to evaporation could be enhanced as the water was retained in the topsoil layer by
the CB system, (2) evaporation loss could be significantly reduced in the CB system byshrinking the irrigation application area, and (3) reduction of RWU and malfunctioning of
the CB could occur when too much irrigation water was applied around the plant as theWater 2021, 13, 2176

17 of 18

excess water accumulated at the CB interface, thus causing breakage of the CB and stressto plant roots due to wetness.
The study suggests that irrigation types, amounts, and schedules need to be carefully
determined to maximize water use efficiency during crop cultivation in a layered soil
consisting of a CB layer. In this context, numerical simulations are helpful and can be used
as a decision support system.
Author Contributions: Conceptualization, H.S. and D.S.; methodology, H.S. and D.S.; model simulations, D.S.; validation, D.S. and H.S; formal analysis, D.S. and H.S.; investigation, D.S. and H.S.;
resources, H.S.; data curation, H.S.; writing—original draft preparation, D.S.; writing—review andediting, H.S., J.Š., T.K. and D.S.; visualization, D.S.; supervision, H.S., T.K. and J.Š.; project administration, H.S.; funding acquisition, H.S. All authors have read and agreed to the published version of
the manuscript.
Funding: This study was partially supported by the Joint Research Program of Arid Land Research Center, Tottori University, and the JSPS KAKENHI Grant Numbers JP17H03885, JP18H02295,
and JP20H03097.
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: The data presented in this study are available on request from thecorresponding author.
Acknowledgments: The authors would like to thank the anonymous reviewers for their comments
to improve the initial manuscript.
Conflicts of Interest: The authors declare no conflict of interest. The funders had no role in the
study’s design, in the collection, analyses, or interpretation of data, in the writing of the manuscript,
or in the decision to publish the results.

References
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.

FAO. The State of the World’s Land and Water Resources for Food and Agriculture (SOLAW)—Managing Systems at Risk; Food andAgriculture Organization of the United Nations: Rome, Italy; Earthscan: London, UK, 2011; ISBN 978-1-136-49888-6.
Ityel, E.; Lazarovitch, N.; Silberbush, M.; Ben-Gal, A. An Artificial Capillary Barrier to Improve Root Zone Conditions for
Horticultural Crops: Physical Effects on Water Content. Irrig. Sci. 2011, 29, 171–180. [CrossRef]
Rooney, D.J.; Brown, K.W.; Thomas, J.C. The Effectiveness of Capillary Barriers to Hydraulically Isolate Salt Contaminated Soils.
Water. Air. Soil Pollut. 1998, 104, 403–411. [CrossRef]
Stormont, J.C.; Morris, C.E. Unsaturated Drainage Layers for Diversion of Infiltrating Water. J. Irrig. Drain. Eng. 1997, 123, 364–366.
[CrossRef]
Morris, C.E.; Stormont, J.C. Evaluation of Numerical Simulations of Capillary Barrier Field Tests. Geotech. Geol. Eng. 1998,
16, 201–213. [CrossRef]
Meshkat, M.; Warner, R.C.; Workman, S.R. Modeling of Evaporation Reduction in Drip Irrigation System. J. Irrig. Drain. Eng.
1999, 125, 315–323. [CrossRef]
Cote, C.M.; Bristow, K.L.; Charlesworth, P.B.; Cook, F.J.; Thorburn, P.J. Analysis of Soil Wetting and Solute Transport in Subsurface
Trickle Irrigation. Irrig. Sci. 2003, 22, 143–156. [CrossRef]Šimůnek, J.; van Genuchten, M.T.; Šejna, M. Recent Developments and Applications of the HYDRUS Computer Software Packages.
Vadose Zone J. 2016, 15, 1539–1663. [CrossRef]
Kämpf, M.; Holfelder, T.; Montenegro, H. Identification and Parameterization of Flow Processes in Artificial Capillary Barriers.
Water Resour. Res. 2003, 39, 1276. [CrossRef]
Fala, O.; Molson, J.; Aubertin, M.; Bussière, B. Numerical Modelling of Flow and Capillary Barrier Effects in Unsaturated Waste
Rock Piles. Mine Water Environ. 2005, 24, 172–185. [CrossRef]
Henry, E.J. Analytical and Numerical Modeling Assessment of Capillary Barrier Performance Degradation Due to ContaminantInduced Surface Tension Reduction. J. Geotech. Geoenviron. Eng. 2007, 133, 231–236. [CrossRef]
Wongkaew, A.; Saito, H.; Fujimaki, H.; Šimůnek, J. Numerical Analysis of Soil Water Dynamics in a Soil Column with an Artificial
Capillary Barrier Growing Leaf Vegetables. Soil Use Manag. 2018, 34, 206–215. [CrossRef]
Al-Mayahi, A.; Al-Ismaily, S.; Al-Maktoumi, A.; Al-Busaidi, H.; Kacimov, A.; Janke, R.; Bouma, J.; Šimůnek, J. A Smart Capillary
Barrier-Wick Irrigation System for Home Gardens in Arid Zones. Irrig. Sci. 2020, 38, 235–250. [CrossRef]
Elmaloglou, S.; Soulis, K.X.; Dercas, N. Simulation of Soil Water Dynamics under Surface Drip Irrigation from Equidistant Line
Sources. Water Resour. Manag. 2013, 27, 4131–4148. [CrossRef]

Water 2021, 13, 2176

15.

16.

17.

18.
19.
20.
21.
22.

23.
24.
25.
26.

27.
28.
29.
30.
31.
32.
33.
34.
35.

18 of 18

Javadzadeh Shakhali, F.; Khaledian, M.; Navabian, M.; Shahinrokhsar, P. Simulation of Soil Volumetric Water Content by
HYDRUS-2D under Line Source with Different Salinities of Irrigation Water and Comparison with Field Observations. J. Water
Soil Sci. 2016, 20, 129–140. [CrossRef]
Arbat, G.; Cufí, S.; Duran-Ros, M.; Pinsach, J.; Puig-Bargués, J.; Pujol, J.; Ramírez de Cartagena, F. Modeling Approaches for
Determining Dripline Depth and Irrigation Frequency of Subsurface Drip Irrigated Rice on Different Soil Textures. Water 2020,12, 1724. [CrossRef]
Honari, M.; Ashrafzadeh, A.; Khaledian, M.; Vazifedoust, M.; Mailhol, J.C. Comparison of HYDRUS-3D Soil Moisture Simulationsof Subsurface Drip Irrigation with Experimental Observations in the South of France. J. Irrig. Drain. Eng. 2017, 143, 04017014.
[CrossRef]
Dabach, S.; Shani, U.; Lazarovitch, N. Optimal Tensiometer Placement for High-Frequency Subsurface Drip Irrigation Management in Heterogeneous Soils. Agric. Water Manag. 2015, 152, 91–98. [CrossRef]
Mubarak, I.; Mailhol, J.C.; Angulo-Jaramillo, R.; Bouarfa, S.; Ruelle, P. Effect of Temporal Variability in Soil Hydraulic Properties
on Simulated Water Transfer under High-Frequency Drip Irrigation. Agric. Water Manag. 2009, 96, 1547–1559. [CrossRef]
Miyake, M.; Otsuka, M.; Fujimaki, H.; Inoue, M.; Saito, H. Evaluation of an Artificial Capillary Barrier to Control Infiltration and
Capillary Rise at Root Zone Areas. J. Arid Land Stud. 2015, 25, 117–120. [CrossRef]Šimůnek, J.; van Genuchten, M.T.; Šejna, M. Development and Applications of the HYDRUS and STANMOD Software Packages
and Related Codes. Vadose Zone J. 2008, 7, 587–600. [CrossRef]Šimůnek, J.; van Genuchten, M.T.; Šejna, M. The HYDRUS Software Package for Simulating the Two- and Three-Dimensional Movementof Water, Heat, and Multiple Solutes in Variably-Saturated Porous Media; Technical Manual, Version 2.0; PC Progress: Prague,
Czech Republic, 2011.
Feddes, R.A.; Kowalik, P.J.; Zaradny, H. Simulation of Field Water Use and Crop Yield; Simulation Monographs; Wiley: New York,NY, USA, 1978; ISBN 978-0-470-26463-8.
van Genuchten, M.T. A Closed-Form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils. Soil Sci. Soc. Am. J.
1980, 44, 892–898. [CrossRef]
Durner, W. Hydraulic Conductivity Estimation for Soils with Heterogeneous Pore Structure. Water Resour. Res. 1994, 30, 211–223.
[CrossRef]
Dane, J.H.; Hopmans, J.W. Methods of Soil Analysis: Part 4 Physical Methods: Laboratory Determination of Water Retention; SSSA
Book Series; Dane, J.H., Clarke Topp, G., Eds.; Soil Science Society of America: Madison, WI, USA, 2002; pp. 675–720,
ISBN 978-0-89118-893-3.
Mualem, Y. A New Model for Predicting the Hydraulic Conductivity of Unsaturated Porous Media. Water Resour. Res. 1976,
12, 513–522. [CrossRef]
Hazen, A. Some Physical Properties of Sands and Gravels, with Special Reference to Their Use in Filtration; 24th Annual Report;
Massachusetts State Board of Health: Boston, MA, USA, 1892; pp. 539–556.
Carrier, W.D. Goodbye, Hazen; Hello, Kozeny-Carman. J. Geotech. Geoenviron. Eng. 2003, 129, 1054–1056. [CrossRef]
Vrugt, J.A.; van Wijk, M.T.; Hopmans, J.W.; Šimunek, J. One-, Two-, and Three-Dimensional Root Water Uptake Functions forTransient Modeling. Water Resour. Res. 2001, 37, 2457–2470. [CrossRef]
Skaggs, T.H.; Trout, T.J.; Šimůnek, J.; Shouse, P.J. Comparison of HYDRUS-2D Simulations of Drip Irrigation with Experimental
Observations. J. Irrig. Drain. Eng. 2004, 130, 304–310. [CrossRef]
Corwin, D.L. Evaluation of a Simple Lysimeter-Design Modification to Minimize Sidewall Flow. J. Contam. Hydrol. 2000, 42, 35–49.
[CrossRef]
Till, A.R.; McCabe, T.P. Sulfur Leaching and Lysimeter Characterization. Soil Sci. 1976, 122, 44–47. [CrossRef]
Saffigna, P.G.; Keeney, D.R.; Tanner, C.B. Lysimeter and Field Measurements of Chloride and Bromide Leaching in an UncultivatedLoamy Sand. Soil Sci. Soc. Am. J. 1977, 41, 478–482. [CrossRef]
Sakai, M.; Toride, N.; Šimůnek, J. Water and Vapor Movement with Condensation and Evaporation in a Sandy Column. Soil Sci.
Soc. Am. J. 2009, 73, 707–717. [CrossRef]