Soil Physics

Soil Water Content Distributions between Two
Emitters of a Subsurface Drip Irrigation System
Maziar M. Kandelous*
Dep. of Environmental Sciences
Univ. of California–Riverside
Riverside, CA 92521
currently at
Dep. of Land, Air and Water Resources
Univ. of California–Davis
Davis, CA 95616

Jiří Šimůnek

Dep. of Environmental Sciences
Univ. of California–Riverside
Riverside, CA 92521

Subsurface drip irrigation (SDI) systems are increasingly being used in agriculture in attempts to use the available
water more efficiently. The proper design and management of SDI systems requires knowledge of precise
distribution of water around emitters. We conducted both field and numerical experiments to evaluate the soil
water content distributions between two neighboring emitters when their wetting patterns started to overlap. The
experiments involved SDI systems with emitters installed at different depths and with different spacings along thedrip lateral. The HYDRUS software package was used to analyze the field data, assuming modeling approaches
in which emitters were represented as (i) a point source in an axisymmetrical two-dimensional domain, (ii) a line
source in a planar two-dimensional domain, or (iii) a point source in a fully three-dimensional domain. Results
indicated that SDI systems can be accurately described using an axisymmetrical two-dimensional model only
before wetting patterns start to overlap, and a planar two-dimensional model only after full merging of the wetting
fronts from neighboring emitters. A fully three-dimensional model appears to be required for describing subsurface
drip irrigation processes in their entirety.
Abbreviations: SDI, subsurface drip irrigation.

M. Th. van Genuchten
Dep. of Mechanical Engineering
Federal Univ. of Rio de Janeiro
Rio de Janeiro, RJ 21945-970, Brazil

Keyvan Malek
Dep. of Irrigation
Bu-Ali Sina Univ.
Hamedan, IranS

ubsurface drip irrigation systems are becoming increasingly popular as a means
for supplying irrigation water, fertilizers, and pesticides more efficiently, andhave been used on a range of vegetables, field crops, and fruit trees (Camp, 1998;
Singh et al., 2006). Increasing the efficiency of water use is a major issue because
competition between agricultural, municipal, and industrial users is likely to increase in the near future. As noted by Skaggs et al. (2004), realizing the full potential of subsurface drip technologies requires optimizing the operational parametersthat are available to irrigators, such as the frequency and duration of irrigation, theemitter discharge rate and spacing, and the placement of the drip laterals. Thus, theproper design and management of SDI systems requires knowledge of the precise
distribution of water around the emitters to provide an optimal distribution ofwater in the crop root zone without undue wetting of the soil surface and drainage
to the groundwater.
During the past several decades, much research has been performed todescribe water flow from point and line sources into the soil in order to designand efficiently manage SDI systems. This research included the development ofboth empirical and analytical models to describe soil water content patterns orwetting fronts for both surface and subsurface drip irrigation systems (e.g., Philip,
1968; Warrick, 1974; Schwartzman and Zur, 1986; Chu, 1994; Moncef et al.,
2002; Cook et al., 2003). Due to increasing computer speed and the availabilityof more comprehensive numerical models for simulating flow in variably saturated
soils, numerical approaches are now increasingly being used for evaluating waterflow in SDI systems (e.g., Taghavi et al., 1984; Angelakis et al., 1993; Meshkat
et al., 1999; Schmitz et al., 2002; Ben-Asher and Phene, 1996; Cote et al., 2003;
Soil Sci. Soc. Am. J. 75:488–497
Posted online 3 Jan. 2011
doi:10.2136/sssaj2010.0181
Received 26 Apr. 2010.
*Corresponding author (mkandelous@ucdavis.edu; mkandelous@gmail.com).
© Soil Science Society of America, 5585 Guilford Rd., Madison WI 53711 USA
All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by
any means, electronic or mechanical, including photocopying, recording, or any information storage
and retrieval system, without permission in writing from the publisher. Permission for printing and forreprinting the material contained herein has been obtained by the publisher.

488

SSSAJ: Volume 75: Number 2 • March–April 2011

Skaggs et al., 2004; Lazarovitch et al., 2005, 2007; Provenzano,
2007; Kandelous and Šimůnek, 2010a,b). Several recent studies
used for this purpose the two-dimensional version of HYDRUS(Šimůnek et al., 1999, 2006), which is a Microsoft Windowsbased computer software package for simulating water, heat, andsolute movement in both two- and three-dimensional variably
saturated porous media. Some of these studies simulated the SDI
process as an equivalent line source (a lateral) (e.g., Skaggs et al.,
2004; Ben-Gal et al., 2004), while others approximated SDI bymeans of a point source (an individual emitter) (e.g., Provenzano,
2007; Lazarovitch et al., 2005, 2007; Kandelous and Šimůnek,2010a,b). All of these studies were performed using either planar
or axisymmetrical two-dimensional models, while in actuality
SDI is a fully three-dimensional flow problem.
In this study, we addressed the types of errors that are madewhen a fully three-dimensional SDI flow process is approximatedby an apparently two-dimensional process. Figure 1A shows
simulated water content distributions in a three-dimensional SDI
soil profile. For our simulations, we used a range of geometries andsoil hydraulic properties, which are discussed below. Figure 1B
illustrates that when water flow from a lateral is approximated asa two-dimensional process, it is assumed that water is discharged
uniformly along the entire length of the lateral. Water is, in fact,
discharged from individual emitters located along the lateral,
while no discharge occurrs from the lateral between neighboring
emitters (Fig. 1A). Therefore, the assumption of two-dimensional
flow is only an approximation of the real flow process. In fact, before
the wetting patterns of two adjacent emitters begin to overlap (Fig.
1C), the three-dimensional flow process could be approximated wellby an axisymmetrical two-dimensional process, and hence solved
using a two-dimensional code like the two-dimensional version of
HYDRUS (Šimůnek et al., 2006), the older HYDRUS-2D code(Šimůnek et al., 1999), or some other simulator. When two adjacent
wetting patterns begin to overlap, however, the problem becomes afully three-dimensional problem and hence needs to be solved usinga three-dimensional numerical simulator. Depending on emitter
spacing, irrigation duration, initial water content, and soil properties,
the wetting patterns of two adjacent emitters will eventually overlap
to a such extent that the water distribution between two emittersalong a lateral becomes relatively uniform (Skaggs et al., 2004), as
shown by plots IIA and IIB in Fig. 1. Only from this point onward
can the flow process be approximated as a geometrically twodimensional problem.
Several of these SDI geometry issues were explored indetail in this study, both experimentally and numerically, using
HYDRUS. The experimental studies produced soil water content
distributions of SDI systems with emitters installed at different
distances and depths. The observed water content distributions
were compared with numerical simulations performed assuming
approximate axisymmetrical two-dimensional and fully threedimensional geometries. We also compared the results ofhypothetical numerical simulations performed using planar
two-dimensional, axisymmetrical two-dimensional, and fully
three-dimensional geometries of the SDI system, with emittersinstalled at different depths and assuming different irrigation times,
once they were calibrated to the experimental data.
The two main objectives of this study were (i) to validate
HYDRUS-3D for modeling water movement in SDI systems by
evaluating the measured soil water content distributions fromseveral field experiments, and (ii) to compare simulated watercontent distributions for SDI systems with emitters installed at
different depths and with different irrigation times using planar
two-dimensional, axisymmetrical two-dimensional, and fully threedimensional geometries and to discuss any possible differences.

MATERIALS AND METHODS
Field Experiments
Eight field experiments (Table 1) were performed at the research
field of the College of Agricultural and Natural Resources of the
University of Tehran, Iran, on a clay loam soil (32.5% clay, 36.5% silt,
and 31.0% sand). Each experiment was conducted using two emittersinstalled at the same depth and having the same discharge rate. We
used different emitter depths (20, 25, and 30 cm), which are typicalof SDI systems for agronomic, turf, and tree crops (Camp, 1998).
Spacings between two neighboring emitters along a lateral were 30, 40,
and 45 cm, which are common for laterals produced by commercial
irrigation companies. Finally, discharge rates commonly used in SDI
systems were selected (e.g., Cote et al., 2003).
The field installation of the emitters used in this study wasdescribed in detail by Kandelous and Šimůnek (2010a). Water from
Fig. 1. Distribution of the soil water content in a subsurface drip irrigation
system simulated as (A) a three-dimensional system with two point sources, a reservoir (consisting of a 160-mm-diameter polyvinyl chloride pipe)(B) a two-dimensional system with a line source, and (C) an axisymmetrical was delivered to the emitters with millimeter precision using a small
two-dimensional system with a point source. Simulations with small (I, top) pump. The volume of water in the water reservoir was recorded every
and large (II, bottom) amounts of applied water are shown.

SSSAJ: Volume 75: Number 2 • March–April 2011489

Table 1. Physical conditions for Experiments 1 through 8 in the field.
Parameter
Emitter installation depth, cm
Distance between emitters, cm
Discharge variation, ×10−3 m3 h−1
Discharge duration, min
Irrigation amount, ×10−3 m3 (two drippers)

125
45
1.7–2.98
300
22.6

2
30
40
1.86–2.56359
25.6

30 min to measure variations in the water discharge rate. The average
discharge rate of individual experiments varied from 1.3 × 10−3 to 2.4 ×
10−3 m3 h−1. Discharge rates gradually decreased during most irrigation
experiments, probably due to the development of overpressure aroundemitter outlets (Lazarovitch et al., 2005). The emitter discharge
variations, durations of irrigation, and the amounts of applied water are
summarized in Table 1.
Soil surrounding the emitters was excavated immediately after
each irrigation experiment to expose a vertical soil profile through ahypothetical lateral between two emitters for gravimetric soil watercontent measurements. Soil samples for this purpose were taken using
a 3-cm-long steel soil sampler with a 3-cm inside diameter at locations
0, d/4, d/2, 3d/4, and d centimeters away from the emitters (see Fig. 2),
where d is the distance between two emitters, at the 10-, 20-, 30-, 40-,
and 50-cm depths. In each case, 25 samples were collected by pressing
the soil sampler horizontally into the soil profile at selected locations.
Because the wetting patterns were presumably horizontally symmetricin the vertical plane between two emitters, samples taken at mirror
positions were mixed after the individual sampling. Volumetric watercontents were determined by multiplying the gravimetric water contents
by the measured bulk density of 1.55 g cm−3.
Except for the saturated water content, θs, the soil hydraulicparameters (van Genuchten, 1980) were initially estimated using the
ROSETTA pedotransfer functions of Schaap et al. (2001) as implemented
in the HYDRUS codes. The saturated water content was estimated directlyin the laboratory using the collected field soil samples. To obtain better
simulations of the field data, we further optimized the van Genuchten
soil hydraulic parameter n, and the saturated hydraulic conductivity Ks,

Fig. 2. Schematic of the field experimental setup, showing emittersand sampling distances, where d is the distance between emitters.
490

3
25
40
1.86–3.24
473
36.6

4
20
40
2.12–2.88
450
35.8

5
30
30
0.54–1.92
363
15.8

6
30
30
1.8–2.5
276
22.2

7
25
30
2.13–2.58
450
32.9

8
25
30
2.04–2.96
335
26.3

against the measured soil water contents using data from a control drip
irrigation experiment by running the inverse parameter estimation optionin HYDRUS. The control experiment and parameter optimization were
described in detail by Kandelous and Šimůnek (2010a). The final estimatedsoil hydraulic parameters used in this study for the clay loam soil were:residual water content θr = 0.07 m3 m−3, θs = 0.38 m3 m−3, shape parameter
α = 1.0 m−1, n = 1.89, Ks = 5 × 10−6 m s−1, and shape parameter l = 0.5.
The measured initial volumetric water content was 0.13 m3 m−3.

Numerical Experiments
Numerical flow experiments were performed considering different
geometries of the SDI system. Soil water contents were simulated
assuming laterals installed at the 20-, 30-, or 40-cm depth and located100 cm apart, with emitter spacings of 50 cm (Fig. 3) and an emitter
discharge rate of 1.60 × 10−3 m3 h−1. Irrigation in all simulations wasapplied during a 24-h period. Soil hydraulic parameters and the initialconditions used in all numerical experiments were the same as thoseused for the field experiments. Observation nodes (27 nodes) werelocated at nine depths between the soil surface and a depth of 80 cm (atintervals of 10 cm) and at three horizontal distances from the lateral (0,12.5, and 25 cm) to obtain water contents at seven times (2, 4, 6, 8, 10,
12, and 24 h) for the statistical analysis comparing the numerical results
obtained for the different SDI scenarios.

Fig. 3. Schematic of the flow domain used in the numericalexperiments. The x and y axes are horizontal axes perpendicularand parallel to the direction of a lateral, respectively. The z axis is avertical axis positive upward.
SSSAJ: Volume 75: Number 2 • March–April 2011

Water flow in the variably saturated SDI soil systems was describedusing the following modified form of the Richards equation:öù¶θ ¶ é æ A ¶h=
+ K izA ÷ úê K ç K ij÷ú¶ t ¶ x i êë çè
¶x j
øû

[1]

where θ is the volumetric water content [L3 L−3], h is the pressure
head [L], xi (i = 1,2 for a two-dimensional model and i = 1,2,3 for athree-dimensional model) are spatial coordinates [L], t is time [T], KijA
and KizA are components of the dimensionless anisotropy hydraulicconductivity tensor KA (dimensionless), and K is the unsaturated
hydraulic conductivity function [L T−1] given by

K ( h, x i ) = K s ( x i ) K r ( h, x i )[2]where Kr is the relative hydraulic conductivity (dimensionless) and
Ks is the saturated hydraulic conductivity [L T−1]. The soil hydraulicproperties were described using the van Genuchten–Mualem
constitutive relationships (van Genuchten, 1980):

θs -θ r
ì
h <0
ïθr +
n m
1+ α hθ ( h )=í
ï
h ³0
îθs

(

)

m
K ( h ) = K s S el é1- ( 1- S e1/m ) ù
ë
û

2

[3]

[4]

where θs is the saturated water content [L3 L−3]; θr is the residual watercontent [L3 L−3]; α [L−1 ], n, and l are shape parameters with m = 1 −1/n; and Se is the effective saturation, given bySe =θ -θ r1, m =1θs -θ rn

[5]

For our studies, we assumed the soil to be isotropic, in which case thedimensionless anisotropy tensor, KA, becomes a unit matrix (i.e., having
values of one along the diagonal and zero values off-diagonal).
The governing flow Eq. [1] within HYDRUS was solved
numerically using Galerkin-type linear finite elements (Šimůnek et al.,1999, 2006, 2008) in both two and three dimensions. This numerical
scheme also allows a two-dimensional version of the Richards equation
to be solved for three-dimensional problems exhibiting radial symmetry
about the vertical axis (such as for flow from a point source). Becausesuch problems involve only two coordinates (i.e., a vertical axis z anda radial axis r), we refer to these problems as axisymmetrical twodimensional problems.

Possible Geometries of the Drip Irrigation System
As shown in Fig. 1, three different geometries could be applied
to a SDI system: a fully three-dimensional system with two dripper
sources (Fig. 1A), an approximate two-dimensional planar system
involving a line source (Fig. 1B), and an approximate two-dimensional
axisymmetrical system containing the dripper sources (Fig. 1C). Each isbriefly discussed below.

Planar Two-Dimensional Geometry
The simplest approach assumes that the SDI system can be
simulated using a planar two-dimensional model in which the lateralis represented as a point source in a two-dimensional domain (Fig. 1B).
The dripper lateral is considered to be an infinite line source in a direction
perpendicular to the simulated plane. This geometry has been used mostfrequently in numerical studies of subsurface drip systems (e.g., Skaggs
et al., 2004; Ben-Gal et al., 2004; Gärdenäs et al., 2005; Lazarovitch etal., 2005,2007; Hanson et al., 2008; Elmaloglou and Diamantopoulos,
2009; Roberts et al., 2009). The line source model, referred to below as
the L2D model (with L standing for line and 2D for two dimensions),
was used only for the numerical experiments and not for analysis of the
field data. The model is applicable to conditions in which the laterals are
farther apart than the emitters along a lateral. In practice, this geometry
can be used for conditions where the flow patterns from individual
emitters overlap and merge to create a combined and relatively uniform
flow pattern from an apparently single infinite line source (compare IIA
and IIB in Fig. 1). Once this happens, the originally radially symmetrical(or axisymmetrical) two-dimensional or three-dimensional flow process
can be simplified to an equivalent two-dimensional line source problem.
Because the flow field in a homogeneous subsurface drip irrigatedfield should be symmetric to the two vertical axes (planes) passing
through the lateral and halfway between the two laterals, only half thespace between two laterals needs to be simulated. The two-dimensionaltransport domain is thus rectangular and comprises half of the length
between the two laterals (in our case 50 cm) and is 80 cm deep. A small
semicircle on the right side of the transport domain represents the lateral(see Fig. 4B). The transport domain for our simulations was discretized
into 4311 triangular finite elements, with their sizes gradually increasing
with distance from the lateral. The smallest element around the lateral
was about 0.1 cm, and the largest on the opposite of the domain was
about 2 cm.

Axisymmetrical Two-Dimensional Geometry

Fig. 4. Flow domains used in the numerical experiments and theirdiscretizations assuming (A) a fully three-dimensional domain (P3D),(B) a planar two-dimensional domain (L2D), and (C) an axisymmetrical
two-dimensional domain (P2D).

SSSAJ: Volume 75: Number 2 • March–April 2011

A slightly more complicated geometry assumes that the SDI system
can be simulated using an axisymmetrical (or radially symmetrical) twodimensional model in which an individual emitter is represented as apoint source located on the axis of rotation (Fig. 1C). This geometry has
been used also in several SDI numerical studies (e.g., Provenzano, 2007;
Kandelous and Šimůnek, 2010a,b). This geometry, referred to as the P2Dmodel (with P standing for point and 2D for two dimensions), was used
for both the field and numerical experiments. The model is primarily491

applicable at early times when flow from any individual emitter can beconsidered to be radially symmetrical. Flow from each emitter remains
axisymmetrical until the wetting patterns from neighboring emitters begin
overlapping each other. After that, the axisymmetrical representation can
only be an approximation of the fully three-dimensional problem insofar
as the numerical results should increasingly depart from reality withincreasing time and irrigation volume.
The transport domain for simulations using this geometry wasagain rectangular (d/2 cm wide and 80 cm deep, where d is the distance
between emitters), except for a semicircle on the left side of the domain
representing the dripper (see Fig. 4C, right side, and Fig. 5, left side). The
transport domain for these simulations was discretized into between
2207 and 3268 triangular elements, depending on the width of thetransport domain. Observation nodes for the numerical experiments
were located at the same depths and radial distances from the source asin the planar two-dimensional model.

that flow was symmetrical about the y–z plane going through the emittersand about the x–z planes going through both the emitter and the middle
between two emitters (see Fig. 3). As a result, the flow domain was 50 cm
wide, 25(d/2) cm deep, and 80 cm high (Fig. 4A) and discretized into
177,000 triangular prismatic three-dimensional elements, again with theirsizes gradually increasing away from the emitter.
Observation nodes in the numerical experiments were, in eachcase, located at the same nine depths between the soil surface anda depth of 80 cm (with an interval of 10 cm) and at three horizontal
distances from the lateral (0, 12.5, and 25 cm in the x direction) as thoseutilized in the two-dimensional simulations. Data were collected at twoplanes in the y direction, i.e., for y = 0 and 25 cm (see Fig. 3). The latter
plane was used to evaluate the soil water content distributions in theplane between two emitters to show the effect of the two overlapping
wetting patterns (see Fig. 3).

Boundary Conditions
Full Three-Dimensional Geometry

A time-variable flux boundary condition was used along boundary
elements representing the emitter during water application in simulations
Subsurface drip irrigation was also simulated using a fully threeof the field experiment. A constant-flux boundary condition was useddimensional model in which each emitter was considered as a point
for the same purpose in the numerical experiments. The time-variable
source. This model is referred to as the P3D model (P for point andboundary fluxes represented the corresponding measured fluxes of the
3D for three dimensions). This model, presumably the most realistic forfield experiments. These fluxes were calculated by dividing the quantities
actual field conditions, was used for analysis of the field data and forof water applied in 30-min intervals by the length or surface area of theproviding comparisons with the approximate planar and axisymmetricalboundary that represents an emitter in the HYDRUS model. A freetwo-dimensional models. Results obtained with this model should
drainage boundary condition was applied along the bottom boundary.
be identical to those for the axisymmetrical two-dimensional modelAll the remaining boundaries were assigned a zero-flux boundary
during the early stages of the infiltration process. Once neighboring
condition. A zero-flux condition was also used along the soil surface
wetting patterns start to overlap or interact, however, flow becomes
because evaporation could be neglected due to the presence of a plastic
fully three dimensional, and only the three-dimensional model should
mulch during irrigation. Boundaries other than those along planes or
be able to describe the complete infiltration process involving multiple
axes of symmetry were assumed to be far enough away to not affect theoverlapping wetting patterns.
flow dynamics in the domain.
Slightly different transport domains were simulated with thisgeometry when analyzing the field experiments in comparison to thenumerical experiments. Field experiments involved
only two emitters, while an infinite line of emitters
was assumed for the numerical experiments. Hence,for the field experiments, we could not assume
that flow was symmetrical about the y–z planepassing through the emitter (because there is thesecond emitter only on one side) but only about
the x–z plane passing through the emitter andthe y–z plane passing through the center between
two emitters. The three-dimensional hexahedral
domain representing half of the domain displayedin Fig. 2 was thus 50 + d/2 cm wide, 50 cm deep,
and 80 cm high (Fig. 5, left), where d is the distance
between two emitters. The emitter was located in
the front plane of the flow domain. The domain
was discretized into 70,000 to 140,000 threedimensional triangular prismatic elements, with
their sizes gradually increasing with their distancefrom the emitter.
For the numerical experiments, for which we Fig. 5. Spatial discretization of the three-dimensional (left) and axisymmetrical twodimensional (right) flow domains used for analysis of the field subsurface drip irrigation
assumed an infinite line of emitters, we could assume
experiments, where d is the distance between emitters.

492

SSSAJ: Volume 75: Number 2 • March–April 2011

Statistical Analyses
Root mean square errors between the simulated and measured
volumetric water contents, or between the simulated water contents
obtained using different SDI geometries, were calculated to provide aquantitative comparison of the correspondence between the measured
and simulated data or between the results of different models. The
RMSE is given byå (y -f )RMSE =
n

i =1

i

2

i

[6]

n

where fi and yi are observed and simulated values, respectively, and n is
the number of observations.

RESULTS AND DISCUSSIONS
Field Experiments
Figure 6 shows the measured and simulated soil watercontent distributions in a vertical plane between two neighboring
emitters of a lateral for the eight field experiments with emittersinstalled at different depths and different spacings. Results
are given for the distributions at the end of the infiltration
experiments. Because of symmetry in the infiltration process, Fig.
6 shows only half of the area between the two emitters. Noticethat the location of the dripper is indicated in each figure andthat the size of each figure reflects the particular emitter spacingused in the various SDI experiments. The MATLAB 4 griddata
method (e.g., Higham and Higham, 2005) was used to grid
the measured water contents obtained at the sampling points.
Subsequently, the contouring scheme of MATLAB was used
to spatially interpolate the gridded water contents. Soil watercontent distributions were simulated for each field experimentusing both the fully three-dimension (P3D) and axisymmetrical
two-dimensional (P2D) versions of HYDRUS. The resultsin Fig. 6 show a qualitatively much better agreement between
the measured and calculated distributions for the full threedimensional simulations, especially for Exp. 3, 4, 7, and 8, which
involved relatively large quantities of infiltrated water and hence
showed more significant interactions between the individual
moisture fronts.
The results in Fig. 6 are confirmed by the plots in Fig. 7, whichshow comparisons between the measured soil water contents at
individual sampling points and those simulated using both thethree-dimensional and axisymmetrical two-dimensional models
at vertical transects 0, d/4, and d/2 away from the emitter,
where d is the distance between two emitters along the lateral.
The plots in Fig. 7 show that the three-dimensional model wasable to describe the collected field data much better than theaxisymmetrical two-dimensional model. This is especially true
for those experiments (notably Exp. 3, 4, 7, and 8) that produced
greater interaction between the moisture fronts coming fromtwo neighboring emitters.
Table 2 summarizes the statistical results of comparisons
between the measured and simulated water contents at the
sampling points at the end of the infiltration experiments.
Measured soil water contents are again compared with thosesimulated using both the three-dimensional (P3D) andaxisymmetrical two-dimensional (P2D) models. The RMSE
values characterizing differences between observed and simulated
water contents using P3D varied between 0.017 and 0.031, andFig. 6. Comparison of measured and simulated distributions of the volumetric soil water content (m3 m−3) between two emitters for Experiments
1 to 8 (from left to right): measured (I, top); simulated using the three dimensional model (II, middle), and simulated using the axisymmetrical two
dimensional model (III, bottom). The white semicircles indicate the location of a dripper.

SSSAJ: Volume 75: Number 2 • March–April 2011493

between 0.026 and 0.087 for P2D. These values are very similar
to the results reported by Skaggs et al. (2004) and Kandelous andŠimůnek (2010b) for their two-dimensional analyses.
The correspondence between the observed water contents
and those simulated using the three-dimensional model wasmuch better than with those simulated using the axisymmetrical
two-dimensional model. For this reason, we used the threedimensional simulations as our control and compared results
obtained with the P2D model against results obtained using theP3D model. Table 2 shows that RMSE values for the differences
between the P2D and P3D models varied from 0.018 to 0.087.
The very low RMSE value of 0.018 for Exp. 1 was expected
because the wetting patterns of the two neighboring emittersoverlapped only slightly at the end of the simulations, thus
reflecting the fact that the two models produced very similar
results. Much larger RMSE values of 0.027, 0.052, and 0.072
were obtained for Exp. 2, 3, and 4, respectively, where emitterswere installed 40 cm apart at decreasing depths. The larger
RMSE values for the shallower experiments are a direct resultof increasingly overlapping wetting patterns. Overlapping wascaused in part by the longer duration of irrigation for emittersinstalled at shallower depths and in part by the failure of thewetting front to follow an axisymmetrical development whenreaching the soil surface. The shallower installation depth of theemitters of Exp. 4, compared with Exp. 3, caused the wetting
front to reach the soil surface much earlier, causing it to deviate
increasingly from the predicted axisymmetrical distributions
typical of P2D simulations.

For Exp. 5, 6, 7, and 8, in which the emitters were all
installed 30 cm apart, the RMSE values were 0.052, 0.066,
0.080, and 0.087, respectively. These values also reflect theincreased overlapping of the wetting patterns caused by greater
quantities of applied irrigated water and consequently more
substantial differences between two- and three-dimensional
simulations. The larger quantity of applied water caused the
RMSE to increase from 0.052 to 0.066 between Exp. 5 and 6.
The larger RMSE value of Exp. 8 compared with Exp. 5 and 6
reflects the increasing effect of the soil surface because of theshallower emitter installation depth and the larger amount ofapplied water, both of which caused larger deviations between
the two models.
Finally, as shown in Table 1, the amount of applied water inExp. 7 was larger than that of Exp. 8, but the RMSE between theapproximate P2D and complete P3D simulations was smaller.
This is because the axisymmetrical P2D model predicted almost
fully saturated soil profiles at the end of Exp. 8, contrary to theresults gained with the three-dimensional model. The increased
amount of irrigated water for Exp. 7, and thus the higher watercontents predicted with the three-dimensional model, could
then only further reduce deviations between the two simulations.
Finally, a word of caution on interpreting the statisticalanalyses of the simulation results. The statistical analyses
were performed for water contents collected at the end of theinfiltration experiments. As such, they pertain only to watercontent distributions at these particular times, and not thefull transient flow process involved in getting to these final

Fig. 7. Comparison of measured (dots) and simulated (lines) soil water contents between two emitters evaluated at locations involving different
depths and distances from the emitter for Experiments 1 to 8 at the end of the irrigation experiment. The top (I), middle (II), and bottom (III) plots
represent distanced of 0, d/4, and d/2 away from the emitter, respectively, where d is the distance between two emitters. Simulated results wereobtained with the three-dimensional model (P3D, solid lines) and the two-dimensional axisymmetrical model (P2D, dashed lines).
494

SSSAJ: Volume 75: Number 2 • March–April 2011

Table 2. Statistical comparison of observed and simulated water contents at the endof the field infiltration experiments.

the drip line. Figure 8 presents the calculated
distributions for the three geometries for a
No. Emitter depth Emitter spacing Irrigation amount Comparison†
R2
RMSE‡
single emitter installation depth of 30 cm.
−33—————— cm ——————
×10 m
The results are for three vertical cross-sections
Obs. with P3D 0.82
0.021
1
25
45
22.6
0, 12.5, and 25 cm away from, but parallel to,
Obs. with P2D 0.84
0.026
the dripper line.
P3D with P2D 0.94
0.018
Figure 8 indicates no discernable
Obs. with P3D 0.85
0.017
2
30
40
25.6
differences between the results obtained
Obs. with P2D 0.80
0.027
P3D with P2D 0.95
0.027
using the P2D and P3D models at earlyObs. with P3D 0.68
0.022
3
25
40
36.6
simulation times before the wetting patterns
Obs. with P2D 0.57
0.053
of two adjacent emitters began overlapping
P3D with P2D 0.97
0.052
(after about 4 h). As the irrigation timeObs. with P3D 0.82
0.031
4
20
40
35.8
increased, the wetting patterns began toObs. with P2D 0.75
0.074
overlap, causing them to mutually interact.
P3D with P2D 0.97
0.072
From this point onward, water flow in
Obs. with P3D 0.85
0.025
5
30
30
15.8
the y direction along the lateral became
Obs. with P2D 0.82
0.048
different from water flow in the x direction
P3D with P2D 0.94
0.052
perpendicular to the lateral. Water flow was
Obs. with P3D 0.80
0.026
6
30
30
22.2
then no longer radially symmetric, causing
Obs. with P2D 0.62
0.069
a divergence between the P2D and P3DP3D with P2D 0.83
0.066
Obs. with P3D 0.64
0.023
results. This can also be demonstrated using
7
25
30
32.9
Obs. with P2D NA§ 0.080
the RMSE values in Table 3. These values
P3D with P2DNA
0.080
were initially very small for all simulated
Obs. with P3D 0.81
0.019
8
25
30
26.3
depths because the results simulated for theObs. with P2D 0.67
0.087
P2D and P3D geometries were essentially
P3D with P2D 0.90
0.087
identical. The RMSE values increased with† Obs, observed data; P3D, three-dimensional geometry; P2D, axisymmetric twotime, however, as the assumption of radiallydimensional geometry.
symmetric flow was increasingly violated. As‡ RMSE of calculated volumetric soil water contents.
§ NA, not applicable.
irrigation times increased from 2 to 24 h, theRMSE values increased from 0.002 to 0.118,
from 0.002 to 0.098, and from 0.002 to 0.081 for installation
states. If the entire wetting process had been considered, thendepths of 20, 30, and 40 cm, respectively.
the statistical measures used to characterize the comparisonsDifferent dynamics can be observed when comparing the
between the P2D and P3D simulations (notably R2 or RMSE)
results
obtained using the L2D and P3D models (Fig. 8; Table 3).
would have been much better (larger and smaller, respectively)
Notice that Fig. 8 presents two sets of lines for the P3D model.
because the simulated wetting patterns during the early stages of
One line (denoted P3D1) represents the soil water contents in a
infiltration would have been the same for the two models whileplane perpendicular to the lateral at the emitter, while a second
producing differences only when neighboring wetting patterns
line (denoted P3D2) represents the water contents in a planestarted to interact.
perpendicular to the lateral halfway between the two emitters.
Numerical Experiments
As expected, the water contents simulated using the L2D
model are between those in the P3D1 and P3D2 planes (Fig.
To provide greater insight into the effects of the assumed
8). While the L2D geometry led to only one line representing
flow geometry on the resulting soil water content distributions
water contents along the entire lateral, the P3D geometry showsin a subsurface drip irrigated field, a large number of simulations
that soil water contents varied from maximum values in thewere performed with different emitter installation depths andP3D1 plane (at the emitter) to minimum values in the P3D2
emitter spacings along the SDI lateral. Soil water contents wereplane (halfway between the emitters). At early times, water
simulated assuming three different geometric representations
contents in the P3D2 plane remained at their initial values. Anof the SDI system: a two-dimensional model (L2D), an
increase in water content in the P3D2 plane was always delayed
axisymmetric two-dimensional model (P2D), and a threein comparison with the P3D1 plane, at least during irrigation.
dimensional model (P3D). The RMSE values were used to
As the irrigation time increased and the wetting patterns fromquantify the degree to which the two-dimensional models (L2D
two emitters started overlapping, differences between the waterand P2D) approximated the results obtained with the threecontents in the P3D1 and P3D2 planes were reduced again.
dimensional model (P3D). The results for three different emitter
Contrary to the gradually increasing RMSE values between the
installation depths are summarized in Table 3 for distributions
P2D and P3D results, the RMSE for differences between thein a vertical cross-section through the emitter perpendicular to

SSSAJ: Volume 75: Number 2 • March–April 2011495

L2D and P3D models gradually decreased Table 3. Statistical comparison of simulated water contents at the end of the numericalinfiltration experiments using different flow geometries and irrigation durations of 2 to 24 h.
as the initially three-dimensional flow
RMSE‡
process became more two dimensional Emitter Scenarios used indepth
the comparison†
2
h
4
h
6
h
8h
10 h
12 h
24 h
due to the overlapping moisture fronts.
cm
This is shown by the RMSE values for theL2D vs. P3D1
0.032 0.030
0.026
0.022
0.017
0.014
0.009
20
differences between the L2D and P3DP2D vs. P3D1
0.002 0.006
0.016
0.027
0.036
0.044
0.118
models, which decreased as the irrigation
L2D vs. P3D1
0.032 0.029
0.025
0.023
0.019
0.016
0.010
30
times increased from 2 to 24 h from 0.032
P2D vs. P3D1
0.002 0.006
0.016
0.024
0.031
0.040
0.098
L2D vs. P3D1
0.030 0.025
0.024
0.022
0.018
0.016
0.009
to 0.009, from 0.032 to 0.0108, and from
40
P2D
vs.
P3D1
0.002
0.006
0.015
0.023
0.030
0.038
0.081
0.030 to 0.0091 for emitter installation† L2D, planar two-dimensional geometry representing the dripper as a line source; P2D,
depths of 20, 30, and 40 cm, respectively. axisymmetrical two-dimensional model representing the emitter as a point source; P3D1, threeFigure 8 and the results in Table 3 dimensional geometry representing each emitter as a point source along a lateral (values hold forindicate that water flow during the early a vertical cross-section through the drip lateral as shown in Fig. 3).
stages of infiltration can be described ‡ RMSE of calculated volumetric soil water contents.
well using the P2D model, with results
essentially identical to the P3D results. As infiltration proceeded
early stages of infiltration, while the L2D simulations were more
and the wetting fronts of the individual emitters started torealistic at later times. In general, however, the three-dimensional
overlap, the P2D results increasingly diverged from the P3D
simulation provided more realistic results.
results. At that stage, only the fully three-dimensional model
SUMMARY AND CONCLUSIONS
could describe the complex dynamics of the interacting moisture
Eight SDI experiments, performed under field conditions,fronts between the emitters. As irrigation proceeded further,were used to follow and study in detail the development of initiallyhowever, the P3D results started to resemble more and more
separate and later overlapping wetting fronts stemming fromthe L2D results for an infinite line source. The P3D and L2Dneighboring emitters. The observed distributions were simulated
simulations eventually became nearly identical, with the waterusing several geometry options available in the HYDRUS
contents showing nearly uniform distributions along the lateral.
software package, i.e., axisymmetrical two-dimensional flow
From the field and numerical experiments in this study,(P2D) and fully three-dimensional flow (P3D). Soil water
we conclude that the P2D simulations provided a good
approximation of the actual SDI wetting patterns during thecontents were obtained from soil samples collected at various

Fig. 8. Comparison of simulated soil water contents between two emitters obtained assuming three different subsurface drip irrigation domains.
Water contents are shown for emitters installed at a depth of 30 cm for three distances away from the emitter. The top (I), middle (II), and bottom(III) plots represent distances of 0, 12.5, and 25 cm away from the drip line, respectively, and seven different irrigation times (left to right).
Simulated results are for the planar two-dimensional geometry applicable to a line source (L2D, solid lines), axisymmetrical two-dimensional flow
involving a point source (P2D, dashed dotted lines), and a three-dimensional domain with point sources (P3D1 for a cross-section across the drip
line, dashed lines, and P3D2 for a cross-section halfway between two emitters perpendicular to the drip line, dotted lines).
496

SSSAJ: Volume 75: Number 2 • March–April 2011

locations from excavated soil profiles at the end of the irrigation
experiments. The hydraulic parameters that provided the best
fit were obtained using a combination of direct measurements,
pedotransfer function predictions, and inverse analyses of thecollected water content data.
Comparisons of the measured and simulated data showedthat the three-dimensional HYDRUS model described theexperimental data very well, providing levels of accuracy similar
to those reported by Skaggs et al. (2004) for planar geometries
and Kandelous and Šimůnek (2010b) for axisymmetrical twodimensional geometries. Because of overlapping wetting patterns
at the end of the infiltration experiments, the P2D model could
not describe the collected water content data between twoemitters as well as the three-dimensional model. A comparison of
results obtained with the P2D and P3D geometries showed thatthe more the wetting patterns overlapped, the less successful theaxisymmetrical two-dimensional model became in describing
full three-dimensional subsurface irrigation processes.
We also conducted a series of numerical experiments
to compare the results generated with the axisymmetrical(P2D) and planar (L2D) two-dimensional models with thoseobtained using the fully three-dimensional model (P3D). The
results were compared at multiple irrigation times and for threedifferent installation depths. The simulations showed that theaxisymmetrical P2D model could predict water contents wellduring the early stages of infiltration before individual water
patterns started overlapping, while the planar L2D model
performed well at later stages of infiltration when the wetting
patterns of individual emitters were fully merged. Only the fully
three-dimensional model was able to describe the complex soil
water SDI dynamics at all stages of the infiltration process.
We conclude that only the three-dimensional geometry
provided accurate simulations of the water content for theassumed SDI system in a study involving multiple drippers, while
two-dimensional axisymmetrical and line-source geometries
provided good descriptions only before the wetting fronts
started to overlap and after full merging of the wetting patterns,
respectively. Ultimately, the HYDRUS code was found to be avery useful tool for studying the various SDI flow geometries anddynamic wetting patterns.

ACKNOWLEDGMENTS
This paper was based on work supported in part by the BARD
(Binational Agricultural Research and Development Fund) ProjectIS-3823-06 and by the National Science Foundation Biocomplexity
programs no. 04-10055 and NSF DEB 04-21530.

REFERENCES
Angelakis, A.N., T.N. Kadir, and D.E. Rolston. 1993. Time-dependent soil-water
distribution under a circular trickler source. Water Resour. Manage. 7:225–235.
Ben-Asher, J., and C.J. Phene. 1996. Surface and subsurface drip irrigation: An
analysis by a numerical model. Res. Rep. Jacob Blaustein Inst. for Desert
Res., Ben Gurion Univ. of the Negev, Negev, Israel.
Ben-Gal, A., N. Lazarovitch, and U. Shani. 2004. Subsurface drip irrigation in
gravel-filled cavities. Vadose Zone J. 3:1407–1413.
Camp, C.R. 1998. Subsurface drip irrigation: A review. Trans. ASAE 41:1353–1367.
Chu, S.T. 1994. Green–Ampt analysis of wetting patterns for surface emitters. J.

SSSAJ: Volume 75: Number 2 • March–April 2011

Irrig. Drain. Eng. 120:414–421.
Cook, F.J., P.J. Thorburn, P. Fitch, and K.L. Bristow. 2003. WetUp: A software tool
to display approximate wetting pattern from drippers. Irrig. Sci. 22:129–134.
Cote, C.M., K.L. Bristow, P.B. Charlesworth, F.J. Cook, and P.J. Thorburn. 2003.
Analysis of soil wetting and solute transport in subsurface trickle irrigation.
Irrig. Sci. 22:143–156.
Elmaloglou, S., and E. Diamantopoulos. 2009. Simulation of water dynamics
under subsurface drip irrigation from line sources. Agric. Water Manage.
96:1587–1595.
Gärdenäs, A.I., J.W. Hopmans, B.R. Hanson, and J. Šimůnek. 2005. Twodimensional modeling of nitrate leaching for various fertigation scenarios
under micro-irrigation. Agric. Water Manage. 74:219–242.
Hanson, B.R., J. Šimůnek, and J.W. Hopmans. 2008. Leaching with subsurface
drip irrigation under saline, shallow ground water conditions. Vadose Zone
J. 7:810–818.
Higham, D.J., and N.J. Higham. 2005. MATLAB guide. 2nd ed. Soc. Ind. Appl.
Math., Philadelphia, PA.
Kandelous, M.M., and J. Šimůnek. 2010a. Comparison of numerical, analyticaland empirical models to estimate wetting pattern for surface and subsurface
drip irrigation. Irrig. Sci. 28:435–444.
Kandelous, M.M., and J. Šimůnek. 2010b. Numerical simulations of water
movement in a subsurface drip irrigation system under field and laboratory
conditions using HYDRUS-2D. Agric. Water Manage. 97:1070–1076.
Lazarovitch, N., J. Šimůnek, and U. Shani. 2005. System-dependent boundary
condition for water flow from subsurface source. Soil Sci. Soc. Am. J. 69:46–50.
Lazarovitch, N., A.W. Warrick, A. Furman, and J. Šimůnek. 2007. Subsurface
water distribution from drip irrigation described by moment analyses.
Vadose Zone J. 6:116–123.
Meshkat, M., R.C. Warner, and S.R. Workman. 1999. Modeling of evaporation
reduction in drip irrigation system. J. Irrig. Drain. Eng. 125:315–323.
Moncef, H., D. Hedi, B. Jelloul, and M. Mohamed. 2002. Approach forpredicting the wetting front depth beneath a surface point source: Theory
and numerical aspect. Irrig. Drain. 51:347–360.
Philip, J.R. 1968. Steady infiltration from buried point sources and sphericalcavities. Water Resour. Res. 4:1039–1047.
Provenzano, G. 2007. Using HYDRUS-2D simulation model to evaluate wetted
soil volume in subsurface drip irrigation systems. J. Irrig. Drain. Eng.
133:342–349.
Roberts, T., N. Lazarovitch, A.W. Warrick, and T.L. Thompson. 2009. Modelingsalt accumulation with subsurface drip irrigation using HYDRUS-2D. Soil
Sci. Soc. Am. J. 73:233–240.
Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. ROSETTA: A computer
program for estimating soil hydraulic properties with hierarchicalpedotransfer functions. J. Hydrol. 251:163–176.
Schmitz, G.H., N. Schutze, and U. Petersohn. 2002. New strategy for optimizing
water application under trickle irrigation. J. Irrig. Drain. Eng. 128:287–297.
Schwartzman, M., and B. Zur. 1986. Emitter spacing and geometry of wetted soilvolume. J. Irrig. Drain. Eng. 112:242–253.
Šimůnek, J., M. Šejna, and M.Th. van Genuchten. 1999. The HYDRUS-2D
software package for simulating two-dimensional movement of water, heatand multiple solutes in variably saturated media. Version 2.0. IGCWMCTPS 53. Int. Ground Water Model. Ctr., Colorado School of Mines, Golden.
Šimůnek, J., M. Šejna, and M.Th. van Genuchten. 2006. The HYDRUS software
package for simulating two- and three-dimensional movement of water,
heat, and multiple solutes in variably-saturated media. Technical Manual,
Version 1.0. PC Progress, Prague, Czech Republic.
Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2008. Development andapplications of the HYDRUS and STANMOD software packages andrelated codes. Vadose Zone J. 7:587–600.
Singh, D.K., T.B.S. Rajput, D.K. Singh, H.S. Sikarwar, R.N. Sahoo, and T.
Ahmad. 2006. Simulation of soil wetting pattern with subsurface drip
irrigation from line source. Agric. Water Manage. 83:130–134.
Skaggs, T.H., T.J. Trout, J. Šimůnek, and P.J. Shouse. 2004. Comparisonof HYDRUS-2D simulations of drip irrigation with experimental
observations. J. Irrig. Drain. Eng. 130:304–310.
Taghavi, S.A., M.A. Marino, and D.E. Rolston. 1984. Infiltration from a trickle
irrigation source. J. Irrig. Drain. Eng. 110:331–341.
van Genuchten, M.Th. 1980. A closed-form equation for predicting hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.
Warrick, A.W. 1974. Time-dependent linearized infiltration: I. Point sources.
Soil Sci. Soc. Am. J. 38:383–386.

497