Irrig Sci (2014) 32:111–125
DOI 10.1007/s00271-013-0417-x

ORIGINAL PAPERHYDRUS simulations of the effects of dual-drip subsurface
irrigation and a physical barrier on water movement and solutetransport in soils
Mohammad N. El-Nesr • A. A. Alazba •Jiřı́ Šimůnek

Received: 9 March 2013 / Accepted: 16 September 2013 / Published online: 4 October 2013
 Springer-Verlag Berlin Heidelberg 2013

Abstract Subsurface drip irrigation systems, compared toother irrigation systems, enhance the delivery of water andnutrients directly into the root zone. However, in lighttextured soils, certain quantities of water may percolate
below the root zone due to the subsurface position of drip
lines and/or poor management of irrigation systems. The
main objective of this paper is to evaluate three technologies to enhance a spatial distribution of water and solutes in
the root zone and to limit downward leaching. The threetechnologies include (a) a physical barrier, (b) a dual-drip
system with concurrent irrigation, and (c) a dual-drip system with sequential irrigation. To achieve this objective,
we performed computer simulations using the HYDRUS(2D/3D) software for both bare and vegetated soils. The
results indicate that the physical barrier is more efficientthan dual-drip systems in enhancing the water distribution
in the root zone while preventing downward leaching. On
the other hand, the dual-drip system improves water distribution in sandy soils. Additionally, the dual-drip systemwith sequential irrigation, followed by the dual-drip systemwith concurrent irrigation, is the most efficient in limiting
downward leaching of solutes.

Communicated by N. Lazarovitch.
M. N. El-Nesr (&)  A. A. Alazba
Alamoudi Chair for Water Research, King Saud University,
Riyadh, Saudi Arabia
e-mail: drnesr@gmail.com; melnesr@ksu.edu.sa
J. Šimůnek
Department of Environmental Sciences, University of California
Riverside, Riverside, CA, USA

Introduction
Surface and subsurface drip irrigation (SDI) systems are
increasingly being used in arid regions with limited water
resources to irrigate agricultural crops. The SDI systems
have an especially promising future due to their many
advantages. The main advantages of SDI, compared tosurface drip irrigation (DI), include the following: (a) asignificant reduction in evaporation, (b) a direct applicationof water and fertilizers into the root zone (Phene et al.
1987; Reich et al. 2009), and (c) an easier and, for drip
lines, safer operation of machinery.
The proper design of SDI systems requires knowledge of
the water distribution patterns around the emitters thatmatch the root extraction patterns and minimize wetting of
the soil surface and deep percolation (Kandelous et al.
2011). The exact shape of the wetted volume and water
distribution depends on many factors, including soil
hydraulic characteristics, initial soil conditions, discharge
rate, application frequency, root characteristics, evaporation, and transpiration (Subbaiah 2013). Additionally, thewetting pattern depends on the location of the emitter withrespect to the soil surface. For uniform and homogeneous
soils, the theoretical wetting pattern is hemiellipsoidal or
ellipsoidal when the dripper is located on the soil surface orin the subsurface, respectively.
Generally, plants have higher root densities in the upper
part of the root zone. It was reported that in a surfaceirrigated soil profile with a uniform water content, mostplants extract about 40, 30, 20, and 10 % of their waterneeds from corresponding quarters of the root zone (Hansen et al. 1980; Ayers and Westcot 1985; Majumdar 2004).
However, the root extraction pattern from a subsurfaceirrigated soil profile may be different. Despite a similar
wetting pattern under DI and SDI systems, a considerable

123

112

amount of water escapes from the root zone by percolating
downward or dispersing laterally in the soil beyond thereach of roots. El-Berry (1989) reported that the main
avenue for water losses under SDI is deep percolation,
which is highest during the seedling stage and declines
with the growth of the root system.
Several investigators have tried to adjust the shape of
the wetted zone to better match the root extraction pattern(e.g., Phene et al. 1987; Barth 1995; Welsh et al. 1995;
Elawady et al. 2003; Ismail et al. 2006). Phene et al.
(1987) showed that the wetted pattern around the buried
emitter can be managed by adjusting the irrigation frequency and that more water moves toward the soil surface
when the irrigation frequency is increased. Other investigators (e.g., Barth 1995; Welsh et al. 1995; Brown et al.
1996) have suggested placing an impermeable barrier
below the lateral drip lines. This barrier can be made ofpolyethylene (Barth 1995; Brown et al. 1996) or metal(Welsh et al. 1995). It has been reported by several
investigators that the physical barrier leads to a significant
increase in crop yield in sandy soils. For example, anincrease in spinach yield by 18 % was reported by
Elawady et al. (2003) and in tomato yield by 141–190 %by Awady et al. (2008). ElNesr (2011) reported for a
Sinai sandy dunes soil an increase in tomato yield by
119 % and in Jerusalem artichoke yield by 131 %. Theselarge increases in crop yield due to the presence of aphysical barrier were attributed to the high infiltration rateof the sandy soil (e.g., 67 cm/h in the experiments of
ElNesr (2011)), with the physical barrier preventing
downward water percolation and keeping water in the rootzone. On the other hand, there are several problems
concerning this physical barrier. For example, there aretechnical and economic problems with digging a wide anddeep trench to install the physical barrier. There may also
be physiological problems for the root system, related toroot rot and/or shallow root disease. Additionally, there
are potential hazards of salt accumulation and other toxicity problems related to accumulation of fertilizers andother chemicals. Hence, the physical barrier should not be
used except in soils with very high infiltration rates.
Several investigators have attempted to achieve similargoals using capillary barriers, i.e., improved water content
distributions (Kampf et al. 1998 and Morris and Stormont
1998) and increased crop yields (Ityel et al. 2010, 2011). A
capillary barrier is an interface between two soil layers
having distinct differences in hydraulic characteristics,
such as a fine soil layer overlaying a coarse-textured soil
(Kampf et al. 1998). Ityel et al. (2010) found that thepresence of a capillary barrier increased volumetric soil
water contents by 20–70 %, depending on the soil texture
and a depth of a barrier, and corresponding yields of pepper
and lettuce by 25 and 36 %, respectively, while having a

123

Irrig Sci (2014) 32:111–125

negligible effect on the yields of tomatoes and melons(Ityel et al. 2011).
Ismail et al. (2006) attempted to modify the wetting
pattern by burying a secondary drip line beneath the primary one and by dividing the required water volume
between the two drip lines. They hypothesized that water
moves faster into the dry soil (due to the higher pressure
head gradient) than into the moist soil, and thus, when thesecondary drip line moistens the soil below the primarydrip line, it forces water from the upper drip line toredistribute upward and laterally, rather than moving
downward. Hence, they called this technique ‘‘a hydraulicbarrier’’. This technique requires no wider trenching than
the normal SDI trenching and prevents all the above-discussed risks to roots. Note that water applications may be
adjusted between the upper and lower emitters dependingon root depth and root density, and more water may be
applied through the upper emitter during the early growthstages when the plant roots are shallow. Ismail et al. (2006)
reported two complimentary irrigations during the first
2 weeks of the plant growth. Their results showed thatwhen applied in the field, the hydraulic barrier increased
the total and marketable yields of the Jerusalem artichokesby 12 and 48 %, respectively. These results clearly document the benefits of using such techniques to increase crop
yields under specific circumstances.
In order to efficiently design and manage SDI systems,
several models (analytical and empirical) have beendeveloped to describe water flow from an emitting source(a point or line source) in the soil (surface or subsurface)(e.g., Brandt et al. 1971; Warrick 1985; Revol et al. 1997;
Khalifa et al. 2004; Singh et al. 2006). One of the mostcomplete packages for simulating water, heat, and solute
movement in both two- and three-dimensional, variably
saturated, porous media is the HYDRUS software package(Šimůnek et al. 2008). Many investigators have used thismodel to evaluate either field or laboratory experiments, orother mathematical models (e.g., Skaggs et al. 2004; Lazarovitch et al. 2005; Provenzano 2007; Zhou et al. 2007;
Kandelous and Šimůnek 2010a, b; Ramos et al. 2012). The
HYDRUS model enables its users to trace the movement ofwater and solutes and the wetting patterns in both simple
and complex geometries for homogeneous or heterogeneous soils and for different combinations of initial andboundary conditions (BCs).
The main objectives of this study therefore are (a) tosimulate water flow and solute transport for an SDI system
while considering both the physical barrier and the dualdrip system using the HYDRUS package and (b) tonumerically evaluate how these techniques affect the
movement of water and solutes in the soil and whether or
not the dual-drip system can act as a hydraulic barrier for
movement of water and solutes.

Irrig Sci (2014) 32:111–125

113

Materials and methods

respectively, a(h) is a dimensionless water stress response
function for water uptake by plant roots (Feddes et al.
1978), b(r, z) is a function describing the spatial root distribution (Vrugt et al. 2001) [L-3], Tp is the potential
transpiration rate, and AT is the surface area associated withtranspiration [L2].
The advection–dispersion equation governing thetransport of a single non-reactive ion in a homogeneous,radially symmetric, porous medium can be described as:
ohc 1 oocoocoðqr cÞ oðqz cÞ
¼
rhD
hD

þ

ot
r or
or
oz
oz
or
oz

SDI modeling
Soil water and solute distributions around a buried drip
source were simulated using the numerical model HYDRUS (2D/3D) (Šimůnek et al. 2011). The model simultaneously solves numerically the Richards equation and theadvection–dispersion equation for variably saturated waterflow and solute transport in soils, respectively. The Richards equation governing water flow from a point source
through variably saturated, porous media can be written inaxisymmetric coordinates as follows:
oh 1 o
oh
o
oh
oKðhÞ
¼
rKðhÞ
KðhÞ
 SðhÞ
þ
ot r or
or
oz
oz
oz
ð1Þ
where h is the volumetric water content [L3L-3], t is thetime [T], h is the soil water pressure head [L], r is the radial(horizontal) coordinate [L], z is the vertical coordinate thatis positive upward [L], K(h) is the unsaturated hydraulicconductivity [LT-1], and S(h) is the sink term representing
root water uptake expressed as a volume of water removed
from a unit volume of soil per unit time [L3L-3T-1].
The soil water retention was modeled using the van
Genuchten equation (van Genuchten 1980),

ð5Þ
where c is the concentration in the liquid phase [ML-3],
D is the effective dispersion coefficient [L2T-1], and qr andqz are the volumetric flux densities in r and z directions,
respectively, [LT-1].
Modeled scenarios
To address our main objectives, we have chosen to evaluate the following alternative scenarios:
1.
2.
3.

Two soil textures (sand and loam).
Bare or vegetated soil profiles.
Different water applications:3.1h  hr1
Se ¼¼
hs  hr ð1 þ ðahÞn Þm

ð2Þ

where Se is the effective degree of saturation or the reduced
water content [–], hr and hs are the residual and saturated
water contents [L3L-3], respectively, a is an empiricalparameter [L-l] inversely related to the air entry value, andn and m are empirical constants affecting the shape of theretention curve [–]. The value of m is restricted by
m = 1 - 1/n.
The hydraulic conductivity as a function of Se wasassumed to be described using the closed form equation of
van Genuchten (1980), which combines the analyticalexpression (2) with the pore size distribution model of
Mualem (1976):
1 m 2m
K ð Se Þ ¼ K s S !
ð3Þ
e 1  1  Se
where  is the pore connectivity parameter estimated by
Mualem (1976) to be about 0.5 as an average for manysoils and Ks is the saturated hydraulic conductivity [LT-1].
The sink term S(h) was computed using the Feddes et al.
(1978) model adapted for a radially symmetric problem:SðhÞ ¼ aðhÞSp ¼ aðhÞbðr; zÞAT Tp

ð4Þ

where S and Sp are the actual and potential (during nostress periods) root water uptake rates [L3L-3T-1],

3.2
3.3

3.4
4.

A single emitter without a physical barrier(denoted below as case a)
Two emitters operating concurrently (case b).
Two emitters operating sequentially (the upper
emitter starts operating 30 min before the lowerone) (case c).
A single emitter with a physical barrier (case d).

Different solute applications:
4.1
4.2

4.3
4.4
4.5

Water and solute are applied from a single
emitter without a physical barrier (case a).
Water is applied from both emitters, and solute isapplied only from the upper emitter. The upper
emitter starts operating 30 min before the lower
one (case c).
Same as 4.2, but the upper emitter starts operating 120 min before the lower one (case cu).
Same as 4.2, but the upper emitter starts operating 120 min after the lower one (case cd).
Water and solute are applied from a single
emitter in the presence of a physical barrier (cased).

Note that the term ‘‘case’’ refers to the operation
sequence (i.e., one emitter with or without a physicalbarrier and two emitters operating either concurrently orsequentially). The same ‘‘case’’ indicates simulations with

123

114

Irrig Sci (2014) 32:111–125

Table 1 Physical properties of soils considered in HYDRUS simulations
Soil
textureBulkdensity(g cm-3)

Residual soil
water content(cm3 cm-3)

Saturated soil
water content(cm3 cm-3)

Parameter a in thesoil water retentionfunction (cm-1)

Parameter n in thesoil water retentionfunction (-)

Saturatedhydraulicconductivity, Ks(cm min-1)

Tortuosity parameter in the conductivity
function (-)

Sand

1.64

0.045

0.43

0.145

2.68

0.49500

0.5

Loam

1.50

0.078

0.43

0.036

1.56

0.01733

0.5

Table 2 Root water uptake parameters for analyzed crops
Crop

Values of the pressure head (cm) below which root water extraction …
…starts. (h1)…occur at the maximum
possible rate (h2)…starts to decline from the maximum
rate at the potential transpiration rateequals…
… RHigh(hHigh)2
3…stops (h4)

… RLow(hLow)
2
3

Limiting potential
transpiration rates (cm/min)Highest (RHigh)2

Lowest (RLow
2 )Tomatoes

-10

-25

-500

-600

-1,0000.008007

6.944e-5

Potatoes

-10

-25

-300

-500

-800

0.00766

6.944e-5

Table 3 Parameters describing a spatial root distribution for analyzed crops
Parameters
Maximum rooting depth (cm)*Potatoes

Tomatoes50

40Depth of maximum root uptake intensity (cm)1015

Maximum rooting radius (cm)80

60

Radius of maximum root uptake intensity (cm)1515Surface area associated with transpiration, AT
(cm2)

4,000

10,000* When no physical barrier is consideredor without solute. The timing of emitters operating
sequentially is indicated by subscripts u and d for case c.
Properties of soils, crops, and solutes consideredin HYDRUS simulations
Simulations were carried out for soils representing two
textural classes: sand and loam. Soil hydraulic parameters
for the two textural classes were taken from the soil catalog
provided by the HYDRUS software (Carsel and Parrish1988) and are given in Table 1.
Parameters representing two model crops, tomatoes(Solanum lycopersicum) and potatoes (Solanum tuberosum), were used in the simulations. The stress response
function of the root water uptake model for selected crops
was defined according to Feddes et al. (1978), and its
parameters are listed in Table 2. Parameters describing the
spatial root distribution according to the model of Vrugt
et al. (2001) were taken from ElNesr et al. (2012) and arelisted in Table 3. These parameters were considered

123

constant with time, since the current version of HYDRUS
does not support a time-variable spatial root distribution.
Also, the same root spatial distribution was considered forboth soil types. Although in reality, the root spatial distribution varies dynamically with the soil water content, theage of the plant, the soil type, and the emitter’s location(Coelho and Or 1999), since the duration of our simulations(either one irrigation cycle for water flow or 30 days forsolute transport scenarios) was relatively short, anassumption of a constant root spatial distribution seemed
adequate.
A non-reactive solute tracer was considered in solutetransport simulations. The soil profile was considered to be
initially solute free (zero concentration), and a unit concentration was applied during solute applications throughthe main (upper) emitter. The longitudinal and transverse
dispersivities were considered with values of 5 and 0.5 cm,
respectively, while the molecular diffusion was neglected(Radcliffe and Šimůnek 2010). Since our interest was
mainly in the transport pathways and leaching of thestudied solute, the root solute uptake was neglected as well.
Cauchy and Neuman BCs were considered at the emitterand outflow boundaries, respectively.
Transport domain, initial and boundary conditions
The HYDRUS (2D/3D) software package (version 2.02)(Šimůnek et al. 2008) is used in numerical simulations for
all modeled scenarios. Water flow, root water uptake, andsolute transport are the processes considered in the simulations. The transport domain for all scenarios was considered to be axisymmetrical around a vertical axis,

Irrig Sci (2014) 32:111–125

115

Fig. 1 Location of the emittersand the physical barrier in thetransport domain (discretizedusing unstructured finite
element mesh) considered in
HYDRUS simulations: adomain around a dripper ismagnified in excerpts.
Dimensions are given in cm

Fig. 2 The transport domain with applied boundary conditions

similarly as done by many other researchers (e.g., Gärdenäs
et al. 2005; Lazarovitch et al. 2005; Kandelous et al. 2011).
Figure 1 shows the detail of the upper left corner of thetransport domain, in which the emitters and a physicalbarrier are located. The transport domain was 100 cm wide
(radius) and 130 cm deep (depth), i.e., large enough so thatan overlap in water content profiles from neighboring
emitters (on other laterals) do not have to be considered.
The (upper) emitter was located 15 cm below the soil
surface (Fig. 1, left), while the secondary emitter (in scenarios b and c) was 25 cm below the soil surface (Fig. 1,
center). The physical barrier (when considered, scenario d)

was placed 27 cm below the soil surface and was considered to have a radius of 25 cm (Fig. 1, right). The soil
profile was considered to be homogeneous (with different
soil types considered in different simulations) and initially
uniformly dry, with the initial pressure head of -1,000 cm.
In the scenarios with solute transport, the soil profile wasconsidered to be initially solute free (zero concentration).
Soil evaporation from bare soil was neglected.
Figure 2 shows the BCs considered in different scenarios in this study. In all simulated scenarios, the upperboundary of the transport domain was subjected to atmospheric conditions, while the lower boundary of the domain
was free drainage. Boundaries at both vertical sides wereassigned a ‘‘no-flux’’ BC. Emitters were represented in allcases as half circles with a radius of 1 cm, located on theleft vertical boundary of the transport domain. The upper
emitter was assigned a ‘‘Variable Flux 1’’ BC. In scenariosb and c, in which two emitters were considered, the second
emitter was assigned a ‘‘Variable Flux 2’’ BC. In scenariod, the physical barrier was simulated as a 1-cm-thickimpermeable barrier 25 cm wide with a ‘‘no-flux’’ BC(Figs. 1d, 2). The transpiration rate was considered to beconstant with time and equal to 4.0 mm/day, i.e., daily
variations of the transpiration rate were not considered.
Since it is common to irrigate vegetables every 2–3 days
(e.g., Allen et al. 1998), simulations were carried out for
2,880 min (i.e., 2 days) when only water flow was considered. When solute transport was also simulated, thesame irrigation conditions (the main cycle) were repeated
15 times, resulting in the total time of 30 days.
Time-variable BCs were used to simulate DI. Since
vegetables are widely irrigated using 1-gph (gallon perhour; 3.7 L/h) pressure-compensating emitters (especiallyin Saudi Arabia where this research was carried out), thisdischarge was used for the dual-drip system(3.7 9 2 % 7.5 L/h). The total discharge (Q) was considered to be 7.5 L/h, which is equivalent to a steady-state flux

123

116

Irrig Sci (2014) 32:111–125

Fig. 3 Water mass balanceerror in a 30-day simulation
scenario (case d, bare sand andloam soils)

Fig. 4 Irrigation fluxes appliedin different scenariosof about 10 cm/min [q = Q/4pr2, where q is the boundary
flux (cm/min), Q is the emitter discharge (cm3/min), andr is the wetted radius of the emitter (cm)]. In scenarios witha single emitter (cases a and d), the irrigation flux was
10 cm/min. In scenarios with two emitters (cases b and c),the irrigation flux at each emitter was 5 cm/min, thus
maintaining the same total discharge in all scenarios. Since
this discharge flux is significantly larger than the saturated
hydraulic conductivity of the two selected soils, the positive pressure heads developed at boundaries representing
the emitters. Although HYDRUS has an option to reduce
the emitter discharge depending on this back pressure(Lazarovitch et al. 2005), in order to avoid the need for123

additional parameters and since we assume the emitters tobe pressure compensating, a constant discharge was considered also during this time period, allowing a buildup of
positive pressures at the emitters. The numerical massbalance was continuously checked in all simulations, with atarget not to exceed the mass balance error of 1.0 %. The
largest recorded mass balance error in all simulations was
0.68 %. An example of the development of the mass balance error in one 30-day simulation is presented in Fig. 3.
The operation sequence of different emitters in different
scenarios is shown in Fig. 4. In cases a and d with only one
emitter, this emitter operates for 60 min every 2 days. In
case b with two emitters, the two emitters have the same

Irrig Sci (2014) 32:111–125

117

720, 1,440, and 2,880 min after the beginning of infiltration. Displayed water content profiles thus cover both theinfiltration and redistribution parts of the numericalexperiment. We will discuss first the results for bare soils
and then for vegetated soils. Since we did not observe any
significant differences between results obtained for two
considered crops, only those obtained for tomatoes will bepresented and discussed below.
Bare, loamy soil

Fig. 5 Locations of analyzed cross sections (CS) in the transportdomainduration of operation, but with only 50 % flux applied toeach of them. Finally, for cases c, cu, and cd with twoemitters operating sequentially, Fig. 4 shows the pattern ofeach sequential operation. The secondary (deeper) emitter
starts operating either 30 or 120 min before the main
(upper) emitter in cases c and cd, respectively, or 120 min
after the upper emitter in case cu. Again, in these sequential
scenarios, the irrigation flux to each emitter is half of theflux in cases a and d.
In the ‘‘Results and discussion’’ section below, we will
provide graphical outputs for pressure heads along verticaland horizontal cross sections throughout the transportdomain. Five vertical cross sections are at radial distancesof 5, 10, 15, 25, and 35 cm away from the axis of symmetry, while seven horizontal cross sections were definedat depths of 5, 10, 20, 30, 40, 50, and 60 cm. The transportdomain and the analyzed cross sections are illustrated inFig. 5.

Results and discussion
Simulation results for water flow
In this section, we will discuss the results of numericalsimulations for different cases and for two soils. For eachsoil, we will discuss both water content profiles alonghorizontal and vertical cross sections for nine output times
with an increasing time interval, i.e., 1, 5, 20, 60, 90, 360,

Water distribution along vertical cross sections Figure 6
shows the vertical water content distributions at five vertical cross sections at distances of 5, 10, 15, 25, and 35 cm
from the emitter(s) for the bare, loamy soil for four analyzed scenarios (a–d). Multiple curves represent outputs at
different times.
The irrigation scheme can be clearly identified from thewater content profiles after 1 min in the 5 cm from theemitter cross section (Fig. 6, left). While water content
profiles for cases a and d clearly show infiltration from asingle emitter, for case b, they show infiltration from twoemitters, and for case c only from the bottom emitter, since
the upper emitter starts operating only after 30 min. In casea, the area around the emitter (depths of 8–25 cm) at a
distance of 5 cm reached saturation before 5 min. This area
(at a 5 cm distance from the emitter) continues to widen
until about 30 cm of the soil profile around the emitter issaturated at the end of irrigation. The water content profile
after 60 min already reflects the redistribution process,since the irrigation flux stops at 60 min in all cases except c(as shown in Fig. 4). The root zone is close to full saturation for the longest period of time for case c (in whichcase, irrigation from the upper dripper stops at 90 min),followed by case d.
Water content profiles are quite similar at a 10 cm distance from the emitter, except in case d, which shows somenoticeable differences in the upper 25 cm. Maximum watercontents are reached in case c at times of 60 and 90 min. At
a 15 cm distance from the emitter, the highest water contents are reached in case d in the upper 25 cm, while thethree other cases are quite similar. At a 25 cm distance,
only small increases in water contents can be observed overtime. While the first increase in the water content isobserved only after 60 min in cases b and c, this occurs
earlier (after 20 min) in cases a and d due to the higher
irrigation flux from a single emitter.
Interestingly, at a 25 cm distance, water content values
appear in reverse order versus time compared to watercontent values at shorter distances. While at this distancewater contents are increasing for larger times, at shorter
distances they are decreasing. This is due to the redistribution process, which drives water both vertically and

123

118

Irrig Sci (2014) 32:111–125

Fig. 6 Vertical water content
distributions at different
distances (5, 10, 15, 25, and 35)
away from the dripper fordifferent irrigation scenarios fora bare, loamy soil

laterally to distances further from emitters. Finally, noticethat water did not reach the 35 cm distance, except in cased between depths of 10 and 25 cm.
Water distribution along horizontal cross sections Water
content distributions for four analyzed cases are shown in
Fig. 7 at seven horizontal cross sections at depths of 5,
10, 20, 30, 40, 50, and 60 cm. The same output times as
in Fig. 6 are shown. In cases b and c, the highest watercontents are found at a depth of 20 cm, i.e., between themain and secondary emitters, which are placed at depthsof 15 and 25 cm, respectively. In this depth, the soil stays
saturated to a horizontal distance of about 20 cm for up to
about 90 min in case c. In case d, water contents arehighest compared to other scenarios at all times at depths
of 5, 10, and 20 cm, while only a small amount of new
water appears at a 30-cm depth at the far end of thephysical barrier at later times (at about 360 min). This isdue to the redistribution process, which is demonstrated inFig. 8. Figure 8 shows that the physical barrier is very
effective in preventing the leaching of water to deeper
depths, as only small quantity of water flows around thebarrier.
While water content profiles for the first three cases a, b,
and c are quite similar in the top 30 cm (Fig. 7), substantial
differences occur at a depth of 40 cm and below. While in123case a, at a depth of 40 cm, the maximum water content isonly about 0.2 cm3 cm-3, in both cases b and c, watercontents are as much as 0.38 at 90 min and 0.3 at 360 min.
Only cases b and c show increases in water contents (up to0.20) at later times due to the redistribution process (attimes of 720, 1,440, 2,880 min). None of the analyzed
cases delivered water down to a 60-cm depth.
Bare, sandy soil
Water distribution along vertical cross sections Sandy
soils have significantly higher hydraulic conductivities andinfiltration rates than loamy soils and significantly lower
macroscopic capillary length (e.g., Radcliffe and Šimůnek2010). This causes water to move deeper and much lesslaterally into the sandy soil profile than in the loamy soil.
Our simulations very clearly show this in Fig. 9. In cases b
and c, water infiltrated down to 60 and 90 cm at a 5 cm
distance at times 60 and 2,880 min, respectively. On theother hand, the rest of the profile is much less saturated.
Although water contents around emitters are higher than
0.40, water quickly redistributes when irrigation stops. The
highest water contents around the emitter at all times canbe observed in case d because of the physical barrier,which prevents redistribution of water downward, while
physical properties of the sandy soil prevent lateralIrrig Sci (2014) 32:111–125

119

Fig. 7 Horizontal water content
distributions at different depths(5, 10, 20, 30, 40, 50, and 60)
for different irrigation scenarios
for a bare, loamy soil

redistribution. In case c, the highest water contents arereached at 90 min when irrigation from the primary emitter
stops.
There are only very small differences in water content
profiles between the first three cases after 360 min. Thisshows that none of the analyzed cases, except case d,
improve the wetting pattern once redistribution in the
sandy soil starts. The same phenomena can be observed atcross sections at a distance of 10 and 15 cm. On the other
hand, in case d, one can observe the accumulation of water
above the physical barrier, as well as its bypassing during
the redistribution process (at the 25- and 35-cm crosssections). This shows that the physical barrier indeed
enhances water redistribution in the root zone.
Water distribution along horizontal cross sections One of
the main benefits of the SDI is that it keeps a dry soil
surface. This helps to reduce the loss of water due toevaporation and prohibits the growth of weeds at the soil
surface. Cases b and c have the lowest water contents in the
top 5 cm of the soil profile (Fig. 10). This is caused mainlyby the lower value of the discharge (5 cm/min) of the upper
emitter in the case of the dual-drip system compared to thesingle-drip system (10 cm/min). This may be considered
one of the benefits of the dual-drip system.
At a 10-cm depth, the highest water contents for cases aand d occur at 20 min, while for case c, the highest watercontents occur at times of 60 and 90 min. No noticeable
differences are visible between cases a–c from time360–2,880 min. In case d, water laterally spreads bothabove (at higher values) and below the physical barrier,moving more than 40 cm at times of 1,440 and 2,880 min,while in other cases, the lateral spreading is at most 30 cm.
At a 20-cm depth, the highest water contents were obtainedin case d at almost all times. This reflects the role of thephysical barrier. Case c has the second highest water

123

120

Fig. 8 Water content distributions at different times for case d for a bare, loamy soil

Fig. 9 Vertical water content
distributions at different
distances (5, 10, 15, 25, and 35)
away from the dripper fordifferent irrigation scenarios fora bare, sandy soil

123

Irrig Sci (2014) 32:111–125

Irrig Sci (2014) 32:111–125

121

Fig. 10 Horizontal watercontent distributions at different
depths (5, 10, 20, 30, 40, 50, and60 cm) for different irrigation
scenarios for a bare, sandy soil

content with a near-saturation state for up to 90 min. This
reflects the effect of the sequential operation of the dualdrip system. At depths of 30 and 40 cm, the highest watercontents were obtained in case c at 90 min, while thelowest water contents were obtained in case d due to theexistence of the physical barrier at a depth of 27 cm. In
case d, there was almost no water content increase at a
depth of 50 cm and no increase at all at a depth of 60 cm.
On the other hand, the dual-drip scenarios (cases b and c)show a better distribution of water at a depth of 50 cm thancase a up to 90 min.
Bare soil comparisons In the loamy soil, water neverreached, except in case d, the radial cross section of 35 cm(Fig. 6). Additionally, in none of the analyzed cases water

moved down to a depth of 60 cm (Fig. 7). On the other
hand, for the sandy soil, Fig. 9 clearly shows that waterreached (although in only small amounts) the radial distance of 35 cm at times 1,440 and 2,880 in all cases, andespecially in cases b and c. Moreover, water also moved
downward and reached depths of 60 and 90 cm (Fig. 10).
Therefore, the commonly observed fact that water moves
deeper and much less laterally in the sandy soil than in the
loamy soil is not exactly true for the conditions studied. In
the sandy soil, in addition to moving deeper in the soil
profile than in the loamy soil, water also moved morelaterally when the dual-drip subsurface system was used.
These findings show the advantage of the dual-drip subsurface systems, i.e., these systems improve the redistribution of water applied in sandy soils.

123

122

Irrig Sci (2014) 32:111–125

Fig. 11 A comparison of water content distributions for a vegetated soil (bottom) and a bare (top) loamy soil after 2 days. The upper left chartshows the root water uptake pattern for tomato crop

Vegetated soil
All simulations for bare soils discussed above were rerunfor soils vegetated with tomatoes and potatoes. Due to thespatial distribution of roots and a corresponding distribution of root water uptake (transpiration) (Fig. 11), crops
appear to affect water content distributions mainly in the
top layer of the soil profile. The root water uptake patternof tomatoes, obtained using parameters listed in Table 2, isshown in Fig. 11. Figure 11 shows that the root systemtakes up water with the highest intensity within a hemisphere with a radius of about 15 cm. This explains theeffects the tomato crop has on the water content patterns
for the 4 analyzed cases.
Figure 11 compares the water content distributions in
the bare (top) and vegetated (bottom) soils at a time of
2,880 min. At that particular time, optimal water contents
near field capacity (0.215 cm3 cm-3) were obtained inalmost an entire root zone for case d. This indicates thatthe frequency of irrigations could be extended to conserve
water. Cases b and c are almost identical, showing theleast amount of water in the top 30 cm of the soil.
Finally, in case a, there is an intermediate amount ofwater in the root zone, higher than in cases b and c andsmaller than in case d. Figure 11 also shows that root
water uptake significantly affects the overall distributionof water only in the upper 10–15 cm. Although Fig. 11
only compares results for a bare, loamy soil and a loamy

123

soil vegetated with tomatoes, similar trends were also
obtained for a bare, sandy soil and a sandy soil vegetatedwith potatoes.
Since all simulations only show small changes in watercontent distributions between the vegetated and bare soils,
this allows us to use bare soil simulations in the solutetransport simulations discussed next.
Simulation results for solute transport
While water content profiles for scenario with a dual-drip
system are different than water content profiles for a scenario with a single dripper, these differences do not seem tobe sufficient to explain large differences in reported yields
of Jerusalem artichokes (Ismail et al. 2006) and tomatoes
and potatoes (ElNesr et al. 2012). Therefore, we have also
explored whether the use of a dual-drip system and thepresence of a barrier affects solute transport and solute
distributions in the root zone. We have rerun many of thesimulations discussed above while also considering solutetransport. Below, we discuss simulations only carried out
for a loamy soil.
Figure 12 shows the distribution patterns of water andsolutes for different irrigation scenarios (defined in Fig. 4).
In addition to commonly used isolines and spectral maps topresent the results for solute concentrations, we also used a
new feature that is available in HYDRUS, i.e., the so-calledflowing particles. These particles show where a solute

Irrig Sci (2014) 32:111–125

123

Fig. 12 Water content and concentration distributions at 30 days calculated for different irrigation scenarios for a bare, loamy soil. Dashed lines
with red dots at the end represent flowing particles, indicating flow paths of a non-reactive solute during the infiltration process

particle released at time zero at the boundary of thetransport domain (here from the emitter) would be at aparticular time, as well as its trajectory. In general, theseparticles thus visualize the pure convective transport of anonreactive solute and thus the furthest possible extent of
the solute transport from the emitter. While they are notintended to represent any particular fertilizer or a potential
crop response, they do provide an indication whether solutes/nutrients will have a tendency to remain in the rootzone or, if highly mobile such as nitrates, to leach below
the root zone.
While the water content distributions (Fig. 12, bottom)
are relatively similar, the concentration patterns are substantially different (Fig. 12, top). Note that in all thesescenarios, solute is injected into the transport domain only
from the main emitter, while pure water is emitted fromthe secondary (lower) emitter. When the upper emitter
starts operating 120 min before the lower emitter (casecu), the concentration isolines are pushed upward forcing
the solute into the top soil layer and its lateral redistribution, i.e., into the soil layers with the highest density ofroots. Only a very small amount of solute leaches
downward. On the other hand, when the lower emitter
starts operating 120 min (case cd) before the main emitter,the concentration isolines show a different behavior. The
isolines tend to be less compressed than in the former
scenario, and a small amount of solute still leaches
downward. However, even in this scenario, most solute
can be found in the top soil layer and laterally redistributed. Scenario c, with almost concurrent operation of thetwo emitters, most closely resembles case a with one
emitter, although even in this case, the solute bulb isdeformed toward upper soil layers.
The effect of the physical barrier on the solute movement is shown in Fig. 13. Figure 13 shows that soluteconcentrations above the physical barrier are almost equal
to the application concentration. This can be beneficial if
applied solutes are useful for plants (such as fertilizers) or
harmful if applied solutes are damaging to plants (such as
emitter-cleaning chemicals or salts). The option to install aphysical barrier thus has to be considered from different
points of view.

123

124

Irrig Sci (2014) 32:111–125

Fig. 13 Water flow and solute transport in the presence of a physical barrier (at 30 days, a bare, loamy soil). Dashed lines with red dots
represent flowing particles (defined in Fig. 12)

Conclusions
Numerical simulations carried out in this study show that theapplication of a dual-drip system or the installation of aphysical barrier can significantly alter both the wetting
pattern and spatial distribution of applied solutes. It can beconcluded that both the physical barrier and the dual-drip
systems are more suitable for the sandy soils than for the
loamy soils, as they help solving some of the problems
caused by the high infiltration rate in the sandy soils.
Physical barriers simply prevent downward movement of all
substances (e.g., water, nutrients, and other chemicals), andthus, their usefulness depends on whether we want to retain
these substances in the root zone or we prefer to flush themout. A dual-drip irrigation system in addition to improving
the redistribution of water applied in sandy soils, it also
represents a powerful tool for manipulating the distributionof solutes in the root zone, especially if the two emitters canbe operated sequentially. Such a system allows growers tocontrol which solute to retain in the root zone and which one
to leach by simply altering the operation of the two drippers.
However, this technology requires much more researchthrough the evaluation of larger numbers of possible scenarios involving different solutes, soils, and operation scenarios, especially through evaluation under field conditions.
Acknowledgments The authors wish to express their gratitude to
the National Plan of Science and Technology at King Saud University
for funding this research by the research project, 10-WAT985-02.
Thanks are also due to Alamoudi Chair for Water Research, where
this research was being carried out. Many thanks and appreciation are123

to the associate editor and three anonymous reviewers for their very
helpful comments and suggestions.

References
Allen RG, Pereira LS, Raes D, Smith M (1998) Crop evapotranspiration: guidelines for computing crop requirements. Irrigationand drainage Paper No. 56. FAO, Rome
Awady MN, Wassif MA, Abd El-Salam MF, and El-Farrah MA(2008) Moisture distribution from subsurface dripping using
saline water in sandy soil. The 15th. Annual Conference of theMisr Society of Ag. Eng., 12–13 March, 2008, pp 477–496
Ayers RS, Westcot DW (1985) Water quality for agriculture.
Irrigation and drainage paper No. 29, FAO, RomeBarth HK (1995) Resource conservation and preservation through a
new subsurface irrigation system, In: Lamm FR (ed) Proceedingsof the 5th international microirrigation congress, April 2–6,
Orlando, FA, ASAE, 168–174
Brandt A, Bresler E, Diner N, Ben-Asher I, Heller J, Goldberg D(1971) Infiltration from a trickle source: I. mathematical model.
Soil Sci Soc Am J Proc 35:675–682
Brown KW, Thomas JC, Friedman S, Meiri A (1996) Wetting
patterns associated with directed subsurface irrigation. In: Camp
CR, Sadler EJ, Yoder RE (eds) Proceedings International
Conference on evapotranspiration and irrigation scheduling. St.
Joseph, Mich, ASAE, pp 806–811
Carsel RF, Parrish RS (1988) Developing joint probability distributions of soil water retention characteristics. Water Resour Res
24:755–769
Coelho EF, Or D (1999) Root distribution and water uptake patternsof corn under surface and subsurface drip irrigation. Plant Soil
206:123–136
Elawady MN, Abd El, Salam MF, Elnawawy MM, El-Farrah MA(2003) Surface and subsurface irrigation effects on Spinach andsorghum. The 4th Annual Conference of Misr Society ofAgricultural Engineers. pp 118–130, Oct 2003

Irrig Sci (2014) 32:111–125
El-Berry AM (1989) Design and utilization of subsurface drip
irrigation system for fodder production in arid lands. Misr J Agr
Eng 6(2):153–165
ElNesr MN (2011) Subsurface drip irrigation development andmodeling of wetting pattern. Lambert Academic Publishing.
ISBN-13: 978-3847339106
ElNesr MN, Alazba AA, Amin MT (2012) Assessing the effect of
three innovative techniques on water conservation and crop yieldthrough subsurface drip irrigation. Project report, King Saud
University
Feddes RA, Kowalik PJ, Zaradny H (1978) Simulation of field water
use and crop yield. Wiley, New York, NY
Gärdenäs AI, Hopmans JW, Hanson BR, Šimůnek J (2005) Twodimensional modeling of nitrate leaching for various fertigation
scenarios under micro-irrigation. Agric Water Manag 74(3):219–242
Hansen VE, Israelsen OW, Stringham GE (1980) Irrigation principles
and practices, 4th edn. Wiley, London, p 417
Ismail SM, Zien El-Abedin TK, Wassif MM, El-Nesr MN (2006)
Physical and hydraulic barriers under surface and subsurface drip
irrigation systems. In: The 14th annual conference of Misr J Agr
Eng MSAE, 23(4):1021–1034
Ityel E, Lazarovitch N, Silberbush M, Ben-Gal A (2010) An artificial
capillary barrier to improve root zone conditions for horticultural
crops: physical effects on water content. Irrig Sci 29(2):171–180.
doi:10.1007/s00271-010-0227-3
Ityel E, Lazarovitch N, Silberbush M, Ben-Gal A (2011) An artificial
capillary barrier to improve root-zone conditions for horticultural crops: response of pepper, lettuce, melon, and tomato. Irrig
Sci 30(4):293–301. doi:10.1007/s00271-011-0281-5
Kampf M, Holfelder T, Montenegro H (1998) Inspection andnumerical simulations of flow processes in capillary barrier
cover systems. In: Holz KP, Bechteler W, Wang SSY, Kawahara
M (eds) Advances in hydro-science and engineering, proceedings of the 3rd international conference on hydro-science and engineering. Brandenburg University, Cottbus
Kandelous MM, Šimůnek J (2010a) Comparison of numerical,
analytical, and empirical models to estimate wetting patterns
for surface and subsurface drip irrigation. Irrig Sci 28(5):435–444
Kandelous MM, Šimůnek J (2010b) Numerical simulations of water
movement in a subsurface drip irrigation system under field andlaboratory conditions using HYDRUS-2D. Agric Water Manag
97:1070–1076
Kandelous MM, Šimůnek J, van Genuchten MTh, Malek K (2011)
Soil water content distributions between two emitters of a
subsurface drip irrigation system. Soil Sci Soc Am J 75(2):488–497
Khalifa HE, El-Gindy AM, Sharaf GA, Allam KhA (2004) Simulatingwater movement in sandy soil under surface point source emitter,I. model development. Misr J Ag Eng 21(2):341–361
Lazarovitch N, Šimůnek J, Shani U (2005) System dependent
boundary condition for water flow from subsurface source. Soil
Sci Soc Am J 69(1):46–50
Majumdar DK (2004) Irrigation water management: principles andpractice. PHI Learning Pvt. Ltd, New Delhi
Morris CE, Stormont JC (1998) Evaluation of numerical simulationsof capillary barrier field tests. Geotech Geolog Eng 16:201–213

125
Mualem Y (1976) A new model for predicting the hydraulicconductivity of unsaturated porous media. Water Resour Res
12(3):513–522
Phene CJ, Davis KR, Hutmacher RB, McCormick RL (1987)
Advantages of subsurface drip irrigation for processing tomatoes. Acta Hortic 200:101–113
Provenzano G (2007) Using HYDRUS-2D simulation model toevaluate wetted soil volume in subsurface drip irrigation
systems. J Irrig Drain Eng 133:342–349
Radcliffe D, Šimůnek J (2010) Soil physics with HYDRUS: modelingand applications, CRC Press. Taylor and Francis Group, Boca
Raton, FL, p 373. ISBN 978-1-4200-7380-5
Ramos TB, Šimůnek J, Gonçalves MC, Martins JC, Prazeres A,
Pereira LS (2012) Two-dimensional modeling of water andnitrogen fate from sweet sorghum irrigated with fresh andblended saline waters. Agric Water Manag 111:87–104
Reich D, Godin R, Chavez JL, Broner I (2009) Subsurface drip
irrigation (SDI). Fort collins: Colorado State University.
http://www.ext.colostate.edu/pubs/crops/04716.pdf. Accessed:
25.08.2013
Revol PB, Clothier E, Mailhol JC, Vachaud G, Vauclin M (1997)
Infiltration from a surface point source and drip irrigation: 2. anapproximate time-dependent solution for wet front position. Wat
Resour Res 33(8):1869–1874
Šimůnek J, van Genuchten MTh, Šejna M (2008) Development andapplications of the HYDRUS and STANMOD software packages, and related codes. Vadose Zone J 7(2):587–600
Šimůnek J, van Genuchten MTh, Šejna M (2011) The HYDRUS
software package for simulating two- and three-dimensional
movement of water, heat, and multiple solutes in variablysaturated media, technical manual, version 2.0, PC progress,
Prague, Czech Republic, pp 258
Singh DK, Rajput TBS, Singh DK, Sikarwar HS, Sahoo RN, Ahmad
T (2006) Simulation of soil wetting pattern with subsurface drip
irrigation from line source. Agric Water Manag 83(1–2):
130–134
Skaggs TH, Trout TJ, Šimůnek J, Shouse PJ (2004) Comparison of
HYDRUS-2D simulations of drip irrigation with experimental
observations. J Irrig Drain Eng 130:304–310
Subbaiah R (2013) A review of models for predicting soil water
dynamics during trickle irrigation. Irrig Sci 31(3):225–258.
doi:10.1007/s00271-011-0309-x
van Genuchten MTh (1980) A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci Soc Am J
44:892–898
Vrugt JA, van Wijk MT, Hopmans JW, Šimůnek J (2001) One-, two-,and three-dimensional root water uptake functions for transientmodeling. Water Resour Res 37(10):2457–2470
Warrick AW (1985) Point and line infiltration calculation of thewetted soil surface. Soil Sci Soc Am J Proc 49:1581–1583
Welsh DF, Kreuter UP, Byles JD (1995) Enhancing subsurface drip
irrigation through vector flow. In: Lamm FR (ed) Proceedings of
the 5th International Microirrigation Congress, Orlando, FA,
ASAE, 2–6 April 1995, pp 688–693
Zhou Q, Kang S, Zhang L, Li F (2007) Comparison of APRI andHYDRUS-2D models to simulate soil water dynamics in avineyard under alternate partial root zone drip irrigation. Plant
Soil 291(1):211–223

123