Published online January 24, 2007

Subsurface Water Distribution from Drip Irrigation Described by Moment Analyses
N. Lazarovitch,* A. W. Warrick, A. Furman, and J. Šimůnek
1986; Schwartzman and Zur, 1986; Angelakis et al.,
1993; Chu, 1994; Zur, 1996; Dasberg and Or, 1999;
Hammami et al., 2002; Thorburn et al., 2003; Cooket al., 2003). The volume of the wetted soil represents
the amount of soil water stored in the root zone. The
domain of interest should be consistent with the anticipated depth of the root system, while its width is associated with the spacing between emitters and lines(Zur, 1996).
A comprehensive method of characterizing spatial–
temporal distributions is through moment analyses. Thisapproach is widely used to describe solute transport in thevadose zone (e.g., Barry and Sposito, 1990; Toride andLeij, 1996; Srivastava et al., 2002). With respect to water, Yeh et al. (2005) and Ye et al. (2005) calculated thezeroth, first, and second moments of a three-dimensional
water content plume and defined an ellipsoid that described the average shape and orientation of the plumefor each observation period. This led to snapshots of theobserved water content plume under transient flow conditions, which was used to derive a three-dimensional effective hydraulic conductivity tensor.
The objective of this study was to implement moment
analyses to describe subsurface water distribution resulting from drip irrigation and to test this approach withnumerically generated data. This includes infiltration
for both line and point sources in contrasting soils. The
considerable advantages gained by the moment analyses over alternatives is detailed.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

ABSTRACT
Moment analysis techniques are used to describe spatial and temporal subsurface wetting patterns resulting from drip emitters. The
water added is considered a ‘‘plume’’ with the zeroth moment representing the total volume of water applied. The first moments lead tothe location of the center of the plume, and the second moments relate
to the amount of spreading about the mean position. We tested thisapproach with numerically generated data for infiltration from surface
and buried line and point sources in three contrasting soils. Ellipses(in two dimensions) or ellipsoids (in three dimensions) can be depicted
about the center of the plume. Any fraction of water added can be
related to a ‘‘probability’’ curve relating the size of the ellipse (or ellipsoid) that contains that amount of water. Remarkably, the probability curves are identical for all times and all of the contrasting soils.
The consistency of the probability relationships can be exploited to
pinpoint the extent of subsurface water for any fraction of the volume
added. The new method can be immediately applied to the vitalquestion of how many sensors are needed and where to install themto capture the overall water distribution under drip irrigation. For
example, better agreement with the ‘‘exact’’ solution occurs with increasing the number of observation points from 6 to 9 and no significant improvement when increasing from 9 to 16. The method can also
be applied to parameter estimation of soil hydraulic properties, which
we uniquely reproduced for generated data.

Ddrip irrigation systems involves selection of an appropriate combination of emitter discharge rate and spacing between emitters for any given
set of soil, crop, and climatic conditions, as well as understanding the wetted zone pattern around the emitter(Bresler, 1978; Lubana and Narda, 2001). Water distribution is affected by many factors, including soil hydrauliccharacteristics, initial conditions, emitter discharge rate,application frequency, root characteristics, evaporation,and transpiration. A traditional way to visualize spatial
and temporal soil water distributions includes determination of the water content at points around the emitter and drawing contours between these points. Practicalinformation includes estimation of the position andshape of the wetted volume (Dasberg and Or, 1999). Previous investigators have used descriptions of the extent of wetting, including the surface wetted diameter,
wetted depth, and wetted volume (Ben-Asher et al.,ESIGNING

THEORYThe three-dimensional spatial moments for a waterplume, Mijk, are defined as (Yeh et al., 2005)
V

Mijk (t) 5

VV# # # u (x,y,z,t)x y z dxdydz
i j k

diff

[1]

2V 2V 2V

with udiff (x,y,z,t) 5 u(x,y,z,t) 2 uinit (x,y,z,t), whereu(x,y,z,t) is the water content at a given time t at alocation x,y,z and uinit(x,y,z,t) is the initial water content,
typically small. (Later we will include an example where
the initial water content is high enough that redistribution of the antecedent water should be taken into
account). The zeroth, first, and second spatial moments
correspond to i 1 j 1 k 5 0, 1, or 2, respectively. The
zeroth moment, M000 is the total volume of water
applied to the domain. The first moments, M100, M010,
and M001 lead to the location of the center of the plume,and the second moments, M200, M020, and M002 relate to
the amount of spreading about its mean position. The
center of mass is at xC,yC,zC, given by

N. Lazarovitch, Wyler Dep. of Dryland Agriculture, Jacob BlausteinInst. for Desert Research, Ben-Gurion Univ. of the Negev, BeerSheva 84990, Israel; A.W. Warrick, Soil Water and Environmental
Science, Univ. of Arizona, Tucson, AZ 85721; A. Furman, Soil, Water
and Environmental Sciences, Agricultural Research Organization,
Volcani Center, Bet Dagan, 50250 Israel; J. Šimůnek, Dep. of Environmental Science, Univ. of California, Riverside, CA 92521.
Received 31 Mar. 2006. *Corresponding author (lazarovi@bgu.ac.il).
Published in Vadose Zone Journal 6:116–123 (2007).
Original Researchdoi:10.2136/vzj2006.0052ª Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

xC 5
116

M100
M000

yC 5

M010
M000

zC 5

M001
M000

[2]

117

www.vadosezonejournal.org

The spread of the plume about its center is describedby the spatial variance in the x, y, and z directions (sx2,sy2, and sz2):

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

M200
s2x 5
2 x2C
M000

M020
s2y 5
2 y2C
M000

M002
s2z 5
2 z2C
M000

[3]

For axial symmetry about r 5 0, where r is the radial coordinate (r 2 5 x 2 1 y2), Eq. [1] can be rewritten by integrating over r and z for the relevant moments. Also by
symmetry, M200 is equal to M020. Thus, M200 and M002 can be
V V

##

2p
M200 (t) 5
udiff (r,z,t)(x2 1 y2 )rdrdz
2 2V 0# # u (r,z,t)r drdz
3

diff

[4]

2V 0
V V

M002 5 2p# # u (r,z,t)rz drdz
2

diff

]h
50
]z

[5]

2V 0

[6]

where u is the volumetric water content [L3 L23], h is the soil
matric head [L], K is the hydraulic conductivity [L T21], x andz are spatial coordinates [L] in transversal and longitudinal
directions, respectively (z positive downward), and t is time [T].
For the surface sources, a specialized version of the code
allowed the spatial distribution of the discharge rate at the soil
surface to change with time (Gardenas et al., 2005). This wasdone by switching from a Neumann (flux) to a Dirichlet (head)boundary condition if the surface pressure head required toaccommodate the specified flux for a surface node is .0. A
sufficient number of surface nodes are switched in an iterative
way until the entire irrigation flux is accounted for, and thepond width is obtained. Since the infiltration flux into the drysoil is larger for early times, the pond width continuously
increases as irrigation proceeds. Other initial and boundary
conditions were used to solve numerically Eq. [1]:2]h(x,z,t)
51
]z

t50

X 1 x0X 2 x0
,x,
,
2
2

z5Z

[10]

u 2 ur5 [1 1 (2ah)n ]2m ,us 2 ur

m 5 121n

K 5 KS Se0:5 [1 2 (1 2 Se 1/m )m ]2

[12]
[13]

where Se is the effective fluid saturation [unitless], h is the soil
matric head [L], us the water contents at saturation [L3 L23],
ur the residual water content [L3 L23], K the hydraulic conductivity [L T21], KS the hydraulic conductivity at saturation[L T21], and a [L21], n [unitless], and m [unitless] are empirical parameters.
The computational domains were selected such that the wetting fronts did not reach the side and the bottom boundaries.
The flow domains (X 5 1 m and Z 5 1.5 m in the twodimensional plan and R 5 1 m and Z 5 1.5 m in the threedimensional axial-symmetric plan) were discretized into
1236 nodes with significantly greater detail around the source.
Three homogenous soils with contrasting hydraulic properties were considered. The hydraulic properties of these soils
were taken from Carsel and Parrish (1988) and summarizedin Table 1. The discharge rate (per unit length) of the surface
line source, Q [L2 T21], was 0.003 m2 h21 and the discharge rate
for the point source, Q [L3 T21], was 0.001 m3 h21. The homogeneous initial effective fluid saturation, Se, was set to 0.01
for all the simulations (except for the cases of where initialconditions were investigated). Simulated application duration was 20 h for the line sources and 60 h for the point sources
with data saved at 20 evenly spaced times. These application
durations lead to an applied volume of 0.06 m3 per 1-m bed inboth line and point sources (assuming one emitter per meter),which is the approximate amount of water applied whenirrigating for the first time for germination purposes. For the
Table 1. Saturated hydraulic conductivity (KS), the van Genuchten–
Mualem fitting parameters a and n, and saturated and residual volumetric water contents (uS and ur, respectively) for thethree representative soils (after Carsel and Parrish, 1988).
KS[7]
z 5 0 [8][11]

where r is a radial coordinate [L]. The initial and boundary
conditions were analogous to the two-dimensional (x,z) cases.
The van Genuchten–Mualem soil hydraulic properties
model (Mualem, 1976; van Genuchten, 1980) was selected forthe numerical simulations:

MATERIALS AND METHODSIn the analyses, we used numerical solutions of the Richards
equation as implemented in the HYDRUS-2D code (Šimůnek
et al., 1999). HYDRUS-2D has been previously used to successfully simulate water flow for drip irrigation systems (Skaggs
et al., 2004; Lazarovitch et al., 2005). We demonstrate themoment analyses using four scenarios, including two-dimensional vertical planes with surface line and buried cylindricalcavity sources and three-dimensional axial-symmetrical geometries with surface point and buried spherical cavity sources.
For the surface line and buried cylindrical cavity sources,
we used the governing flow equation for two-dimensional
isothermal Darcian flow in a variably saturated isotropic
rigid porous medium presented by the mixed form of the
Richards equation:

h(x,z,t) 5 hi

[9]where x0 is the pond width [L], hi is the initial soil pressure
head [L], and Z and X are the length and width of the domain,
respectively [L].
For the subsurface cavity sources, we used the code from
Lazarovitch et al. (2005). This version allows calculation of thesource discharge while considering source properties, inlet
pressure, and effects of the soil hydraulic properties on theactual discharge.
In the axisymmetrical three-dimensional cases we used

Se 5





]u
]
]h
]
]h
5
K
1
K
11
]t
]x
]x
]z
]x

x2 1 z 2 Y V





]u
1]
]h
]
]h
]K
5
rK
1
K
2
]t
r ]r
]r
]z
]z
]z

V V

5p

h(x,z,t) 5 hi

21

Sand
Sandy loam
Loam

ms
25
8.25 3 1025
1.23 3 1026
2.89 3 10

a

n

uS

2.68
1.89
1.56

m m
0.43
0.045
0.41
0.065
0.43
0.078

21

m
14.5
7.5
3.6

ur
3

23

118

VADOSE ZONE J., VOL. 6, FEBRUARY 2007

where B(a,b) is the complete beta function, a and b are parameters, and u 5 k/kmax (an appropriate value of kmax 5 3 isused below). Results for the probability function and its interpretation will be given below.
For calculating Mijk, we assumed that the initial water content, uinit, was sufficiently small that redistribution of theantecedent water was negligible over the time frames considered. For wetter initial conditions, however, significant redistribution of antecedent water can occur due to gravity
drainage. This can lead to negative values for the integrals in
Eq. [1]. For this condition, we propose to use the backgroundwater content (ubg) rather than uinit for calculating udiff in
Eq. [1]. The background water content is time and depth
dependent and is defined as the equivalent water content after
redistribution starting from the same initial conditions. The
calculations are performed for the same domain and initialand boundary conditions but without adding any water.

Data Processing
After completing the simulations, we calculated Mijk, zC, sx,
and sz for the two-dimensional plane using Eq. [1] to [3] (xC
was equal to zero). For the axially symmetric case, Mijk, zC, sx,
and sz were found from Eq. [1] to [5]. A separate program wasused to compute the moments directly from the HYDRUS-2D
output. For expediency, an equally spaced grid was definedand values of water content in the grid points were assumed tobe that of the closest finite element node. The moments werethen computed using the gridded values. To check for massbalance, the values of M000 in all cases were verified to agree
within 0.1% of the applied water.
After calculating the moments for each case and time, we
defined ellipses around the center of mass (0, zC) in the twodimensional plane. The semi-axes of the ellipses are analogous
to the semi-axes of a binormal probability distribution (e.g.,Morrison, 1976, Ch. 3). For the two-dimensional case, we used
ksx and ksz, where k is the number of standard deviations.
Each ellipse was set asx2
z2
1 2 2 51
2
2
k sx
k sz

RESULTS AND DISCUSSION
For the two-dimensional domain, the center of mass, zC,
as a function of time for the three representative soils ispresented in Fig. 1A. As time progresses, water advancesfrom the surface source deeper into the profile, well
depicted by following values of zC. In the beginning of theirrigation event, zC advances fast and after several hours
the rate of the advance becomes nearly constant. The
contrast in texture leads to different retention capacities
for the soils. The coarse-textured sand has the lowest
retention capacity and the deepest zC; the fine-textured
loam has the highest retention capacity with shallowest zC.
The semi-axes in x and z directions, sx and sz (as afunction of time) are depicted in Fig. 1B and 1C, respectively. While sz exhibits differences between the soils,
sx reveals almost the same values for the sandy and the
sandy loam soils but slightly higher values for the loamy
soil at any given time.
Wetting patterns and the ellipses for one and two
standard deviations about zC are illustrated in Fig. 2 at
20 h. Of course, these ellipses are consistent with Fig. 1.
For example, zC 5 20.39 m, sx 5 0.14 m, and sz 50.22 m for the sandy soil with t 5 20 h. These ellipses canbe shown also for early times using the same procedure.
The shape of the ellipses varies with soil texture. Whilethe ellipses in the sandy soil are highly elongated, theellipses in the sandy loam soil are close to circles.
The center of mass, zC, for surface point sources forthe three representative soils is presented in Fig. 3A.

[14]

In the axially symmetric cases, we defined spheroids (symmetrical about the z axis with semi-axes sx, sy 5 sx and sz):

x2
y2
z2
r2
z2
1 2 2 1 2 2 5 2 2 1 2 2 51
2
2
k sx
k sx
k sz
k sx
k sz

[15]

The amount of water within an ellipse was also computed using
the gridded-water content values with Eq. [1]. This was doneusing the same program used to generate the moment values.
Dividing the mass of added water within any ellipse by M000
gives a corresponding fraction of the total water applied. By
repeating the calculations for increasing values of k, the result
is a cumulative probability function, P, going from P 5 0 for
k 5 0 (none of the added water is included within a point) toP 5 1 as k becomes large (all of the added water is included
as an ellipse or spheroid becomes large). Finally, cumulative
beta distribution functions were fitted with

(1 2 u)b21 ua21
B(a,b)

[16]

t (h)
0

5

10

t (h)
15

20

0
0.24

-0.1

0.20

σx (m)

0

-0.2
-0.3
-0.4

Loam
Sandy loam
Sand

0.16

5

10

t (h)
15

20
0.25

Loam
Sandy loam
Sand

0.20

0.12
0.08

A

0.04

0

σz (m)

P(u) 5zc (m)

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

cavity sources, the nominal discharge, Q0, was 0.003 m2 h21 forthe two-dimensional cases and 0.001 m3 h21 for the threedimensional cases. The radius and depth of the sources were
set to 0.01 and 0.15 m, respectively. For expediency, simulations with the subsurface cavity sources were only for the
sandy loam soil.

5

10

15

20

Loam
Loamy sand
Sand

0.15
0.10
0.05

B

C

-0.5
0.00
0.00
Fig. 1. (A) The vertical center of gravity, zC, for a surface point source in a two-dimensional plane as a function of time, t, for three representativesoils; (B) spatial variance in the x direction, sx, as a function of t ; and (C) spatial variance in the z direction, sz, as a function of t.

119Fig. 2. Wetting patterns and ellipses for one (solid line) and two (dashed line) standard deviations about the vertical center of gravity (small filled
circle), for three representative soils under a surface line source.

poses, the same. Increasing the size of the ellipses (orspheroids) includes a higher fraction of the water in theplume. Essentially, all of the added water resides within
the ellipses or spheroids defined by the number of standard deviations k 5 3. The solid curves are fitted using
the beta distribution Eq. [21]. The fitted parameters
for the two-dimensional case are a 5 3.15, b 5 3.98,
and for the three-dimensional case are a 5 4.33 andb 5 3.46.
The fact that the curve for the three-dimensional caseis below the two-dimensional results in Fig. 6 could havebeen anticipated based on an analogy to general probability distributions. To demonstrate this for a simple
case, compare normal, binormal, and trinormal distributions with independent variables. For the normal distribution the probability that a random variable |X| (ofmean 0) is less than the standard deviation is well known
to be 0.683; for two standard deviations, the probability
is 0.954 (see Table 2). Analogous results for the binormal distribution (assuming independence of variables) are found and represent the probability spacewithin a circle of radius ks. The circular space is conveniently approximated by taking a square that has thesame area as the circle or taking |X| and |Y| , 0.5p0.5k.
If the two variables are assumed independent, theindividual probabilities are multiplied together, giving
the expression in the first column of Table 2 and numerical values of 0.39 and 0.853 for one and two stan-The semi-axes in z and x directions, sx and sz, arepresented in Fig. 3B and 3C, respectively. The zC in thethree-dimensional cases shows the same properties as in
the two-dimensional cases in that the coarser the texture, the deeper the zC. Also, the horizontal variances
are nearly independent of texture, whereas the verticalvariance shows differences similar to those for the twodimensional case.
Wetting patterns and the resulting spheroids for thesurface point applications are illustrated in Fig. 4 at
the end of irrigation (60 h). Results are shown in thex–z plane, with the resulting ellipses (these are crosssections of the corresponding spheroids) for one and two
standard deviations (alternatively, the x axis may beconsidered as an r axis). As in the previous examples,the finer texture (loam) results in a distribution patternnearer to the surface, whereas the coarser textured sandhas pattern that is deeper in the profile.
For the subsurface cavity sources, distribution patterns and the resulting ellipses were calculated only forthe sandy loam. Both the two- and three-dimensional
results are shown in Fig. 5 for 40 h.
Cumulative probability as a function of the number ofstandard deviations for various soils, times, and dimensions of the irrigation event is presented in Fig. 6. Although zC, sz, and sx differ for each soil and time, thefraction of the amount of water within ellipses or spheroids from the total applied water is, for practical purt (h)
20

t (h)
40

60

0

0

t (h)
40

60

0
0.40

0.24
0.20

σx (m)

-0.2
-0.4
-0.6

20

Loam
Sandy loam
Sand

0.16

Loam
Sandy loam
Sand

0.30

σz (m)0zc (m)

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

0.12

20

40

60

Loam
Loamy sand
Sand

0.20

0.08
0.10

A

0.04

B

C

0.00
0.00
Fig. 3. (A) The vertical center of gravity, zC, for a surface point source in a three-dimensional plane as a function of time, t, for three representativesoils; (B) spatial variance in the x direction, sx, as a function of t ; and (C) spatial variance in the z direction, sz, as a function of t.
-0.8

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

120

VADOSE ZONE J., VOL. 6, FEBRUARY 2007

Fig. 4. Cross-sections of the wetting patterns and the spheroids for one (solid line) and two (dashed line) standard deviations about the verticalcenter of gravity (small filled circle), for three representative soils under a surface point source.

dard deviations. For the trinormal distribution, a similar
approximation is made that a sphere of radius ks can be
replaced by a cube of the same volume leading to theprobability given in the first column of Table 2 with theresulting values of 0.195 and 0.712. Values from Fig. 6
are shown in parentheses and are slightly smaller for
k 5 1 (e.g., 0.303 compared with 0.390 for the twodimensional case) and slightly larger for k 5 2 (0.892
compared with 0.853). The uniqueness of the probability
curves defined by the ellipses was not anticipated. This
allows very detailed definitions of the subsurface waterusing the time distributions of only three parameters: zC,
sz, and sx. In fact, a reasonable approximation may be
to use a single curve (with time) for the horizontal
variances as evidenced by Fig. 1B and 3B.
We used ubg for calculating zC for a surface point
source with discharge rate of 0.001 m3 h21 in a sandy soil
having initial saturations, Se, of 0.01, 0.05, and 0.1. We
chose the sandy soil since it yields the highest changes
in water contents when the initial water is redistributing in the soil during the irrigation event. The results
are shown in Fig. 7. For Se 5 0.01 the results are consistent with Fig. 3 but for Se 5 0.05 and 0.1, zC is deeper
as expected.

A possible use of the moment analysis method is toaddress the vital question of how many sensors areneeded and where to install them to capture the overall plume behavior in drip irrigation. The traditional
method is to place one to three sensors (tensiometers,
time domain reflectometry probes, etc.), hoping to capture the entire root zone water dynamics (e.g., Dasberg
and Or, 1999). Rarely more than one sensor is placed,
and in such cases the water content is simply averaged.
We suggest here that by calculating spatial moments
and generating ellipses or spheroids as described above,
the true distribution of the water as it resides in thesubsurface can be captured.
Figure 8 shows spheroids about zC for 6, 9, and 16
observation points that were distributed evenly in thedomain of interest. These spheroids’ curves (dashedlines) are comparable to the ‘‘exact’’ spheroids (solidlines) for 50 and 90% of the water in the domain. The
results are for the surface point source into the loamy
soil with a discharge rate of 0.0005 m3 h21 with Se 5 0.01
after 40 h. The ‘‘exact’’ curve was calculated using all thenodal information. In the other cases, the moments were
calculated by assigning an area (calculated with rectangles and presented in Fig. 8 with dashed lines) to eachFig. 5. Wetting patterns and ellipses for one (solid line) and two (dashed line) standard deviations about the vertical center of gravity, zC (small fullcircle), for a subsurface cavity source in two dimensions (2D) and the analogous results for a three-dimensional (3D) cavity for a loamy sand soil.

121

www.vadosezonejournal.org

t (h)

1

10

20

30

40

0

0.8

Se = 0.01

2D

Se = 0.05

-0.2
0.6

zc (m)

Cumulative probability

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

0

3D

0.4

Se = 0.10
-0.4
-0.6

0.2

-0.8

0
0

1

2

3

k

and the unknown soil. The sum of squared errors (F)was calculated with

2D
3D
Loam
20 h 60 h
Sand""
Sandy loam ""Loam
10 h 30 h
Sand""
Sandy loam ""Sandy loam subsurface

i51

2

Ci

Ci

2

zi

zi

ri

2

ri

[17]

observation point (presented as a triangle) and assuming that in this area the water content is constant andequal to the water content of the observation point.
Better agreement with the ‘‘exact’’ solution occurs withan increased number of observation points from six tonine, but no significant improvement when increasingfrom 9 to 16.
The moment analyses can also be applied to parameter estimation of the soil hydraulic properties. Parameter estimation methods have been extensively applied
to determine soil hydraulic properties (Russo et al.,
1991; Hopmans et al., 2002). Usually the input data forthe parameter estimation consist of water content, water
potential, or infiltrations. The knowledge from the moment analyses can serve as an alternative input for theparameter estimation. In the next analysis, we tried toestimate a and n of the loamy soil. We assumed that weknow the rest of the hydraulic parameters of this soil
(i.e., Ks, us, and ur) and we also know the information
about zC(t), sr(t) and sz(t). We used 64 (8 3 8) simulations of surface point sources with Q 5 0.001 m3 h21 tocreate a response surface between the simulated data
Table 2. Probability of being within one (k 5 1) and two (k 5 2)
standard deviations of the mean for normal, binormal, andtrinormal normal distributions. The values in parentheses arefrom Fig. 6.
Normal P(| X | , ks)
Binormal {P[| X | , (p0.5 ks)/2]}2
3
Trinormal {P[| X | , (p/6)1/3 ks]}

O (z 2 z *) 1 (s 2 s *) 1 (s 2 s *)b

F5Fig. 6. Cumulative probability as a function of the number of standard
deviations, k, for various soils, times, source locations, and two (2D)or three dimensions (3D).

Distribution†

Fig. 7. The vertical center of gravity, zC, for the surface point source
as a function of time, t, for the sandy soil with three initial saturations (Se).

k51

k52

0.683
0.390 (0.303)
0.195 (0.118)

0.954
0.853 (0.892)
0.712 (0.724)† P 5 probability, X 5 random variable, s 5 spatial variance.

where b is the number of print times, the asterisks denote the ‘‘measured point,’’ and the subscript i denotes
and index for a particular calculated point. Figure 9shows F as a function of a and n. The minima of theerror is exactly with a 5 3.6 m and n 5 1.56. This uniquereproduction of a and n for the loamy soil suggests theuse of the moment analyses with inverse estimation ofsoil hydraulic properties.

SUMMARY AND CONCLUSIONS
Spatial moments describing subsurface water distribution from drip sources are presented for numericallygenerated data. These analyses allow a straightforward,
physically meaningful description of the general patternof water distribution around surface and subsurface
sources. As opposed to traditional methods, which require detailed information of the water content distribution, moment analyses can accurately describe thewater content distribution in a statistical manner withjust three numbers (zC, sz, and sx). Remarkably, theprobability function relating fraction of water applied toan appropriate encapsulating ellipse (or spheroid) is
not soil or time dependent. Thus, once zC, sz, and sx are
known, the overall boundaries and shape of the wetted
volume can be approximated with a high accuracy. For
any specified fraction of the volume of water added, an
ellipse (or spheroids) is easily defined within which acorresponding proportion of the added water resides.
In implementing moment analyses to subsurface
water distribution, the emphasis has been primarily on
definitions and evaluation of the moments. Unlike otherexamples of moment analysis in vadose zone hydrology,
the results are unique in that here the total water in thesubsurface plume is increasing due to additional water
being added to the system. Numerous questions arise as
to how to address technical problems and how to best
exploit the techniques. Although considerable informa-

122

VADOSE ZONE J., VOL. 6, FEBRUARY 2007

4.8

r (cm)
0

10

20

30

40

0

4.4

4

α

z (cm)

3.6
-20

3.2
-30

2.8

A

-40

2.4

r (cm)
0

10

20

30

40

0

2
1.2

1.3

1.4

1.5

1.6

1.7

1.8

1.9

n
Fig. 9. The sum of squared errors, %, as a function of the fitting parameters n and a (m21) for the hypothetical case of the loamy soil.

z (cm)-10

ACKNOWLEDGMENTS
-20

-30

B

-40

r (cm)
10

20

30

40

0

-10

-20

-30

-40

This work was supported by The United States–Israel
Binational Agricultural Research and Development fund
(BARD), Project Grant Agreement no. US-3662-05R andWestern Research Project W-1188. Grateful acknowledgments
are made to Don Myers and Jim Yeh for helpful suggestions
with respect to the moment definitions.

REFERENCES
0z (cm)

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

-10

C

Fig. 8. Cross-sections of the spheroids (dashed lines) about the verticalcenter of gravity for (A) 6, (B) 9, and (C) 16 observation points(triangles) compared with the ‘‘exact’’ results (solid lines) under asurface point source. The associated Thiessen polygons are showed
as dashed lines.

tion follows from knowing zC, sz, and sx, their values
change with time; hence, more direct methods forevaluation would be extremely useful. The two applications (sensor evaluation and inverse parameter estimation) deserve further attention in the future, which isbeyond the scope of this paper. Other geometries andprocesses, such as flow from channels and boreholes, aswell as for water redistribution and water uptake byplants, also appear to be fruitful topics to receive attention in the future.

Angelakis, A.N., D.E. Rolston, T.N. Kadir, and V.N. Scott. 1993. Soil-water
distribution under trickle source. J. Irrig. Drain. Eng. 119:484–500.
Barry, D.A., and G. Sposito. 1990. Three-dimensional statistical moment analysis of the Stanford/Waterloo Borden tracer data. Water
Resour. Res. 26:1735–1747.
Ben-Asher, J., C. Charach, and A. Zemel. 1986. Infiltration and water
extraction from trickle irrigation source. The effective hemispheremodel. Soil Sci. Soc. Am. J. 50:882–887.
Bresler, E. 1978. Analysis of trickle irrigation with application to
design problems. Irrig. Sci. 1:3–17.
Carsel, R.F., and R.S. Parrish. 1988. Developing joint probability
distributions of soil water retention characteristics. Water Resour.
Res. 24:755–769.
Chu, S.T. 1994. Green–Ampt analysis of wetting patterns for surface
emitters. J. 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 patterns from drippers. Irrig. Sci. 22:129–134.
Dasberg, S., and D. Or. 1999. Drip irrigation. Springer Verlag, Berlin.
Gardenas, A.I., J.W. Hopmans, B.R. Hanson, and J. Šimůnek. 2005.
Two-dimensional modeling of nitrate leaching for various fertigation scenarios under micro-irrigation. Agric. Water Manage. 74:
219–242.
Hammami, M., H. Daghari, J. Balti, and M. Maalej. 2002. Approachfor predicting the wetting front depth beneath a surface point
source: Theory and numerical aspect. Irrig. Drain. 51:347–360.
Hopmans, J.W., J. Šimůnek, N. Romano, and W. Durner. 2002. Simultaneous determination of water transmission and retention properties. Inverse methods. p. 963–1008. In J.H. Dane and G.C. Topp(ed.) Methods of soil analysis. Part 4. SSSA Book Ser. 5. SSSA,
Madison, WI.
Lazarovitch, N., J. Šimůnek, and U. Shani. 2005. System-dependent
boundary condition for water flow from subsurface source. Soil Sci.
Soc. Am. J. 69:46–50.
Lubana, P.P.S., and N.K. Narda. 2001. Soil and water modelling. Soilwater dynamics under trickle emitters—A review. J. Agric. Eng.
Res. 78:217–232.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

Morrison, D.F. 1976. Multivariate statistical methods. McGraw-Hill,
New York.
Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:
513–522.
Russo, D., E. Bresler, U. Shani, and J.C. Parker. 1991. Analyses of infiltration events in relation to determining soil hydraulic properties by inverse problem methodology. Water Resour. Res. 27:
1361–1373.
Schwartzman, M., and B. Zur. 1986. Emitter spacing and geometry ofwetted soil volume. J. Irrig. Drain. Eng. 112:242–253.
Šimůnek, J., M. Sejna, and M.Th. van Genuchten. 1999. The
HYDRUS-2D software package for simulating two-dimensional
movement of water, heat, and multiple solutes in variably saturated
media. Version 2.0. Rep. IGWMC-TPS-53. Int. Ground Water
Model. Ctr., Golden, CO.
Skaggs, T.H., T.J. Trout, J. Šimůnek, and P.J. Shouse. 2004. Comparisonof HYDRUS-2D simulations of drip irrigation with experimental
observations. J. Irrig. Drain. Eng. 130:304–310.
Srivastava, R., P.K. Sharma, and M.L. Brusseau. 2002. Spatial mo-

123

ments for reactive transport in heterogeneous porous media. J.
Hydrol. Eng. 7:336–341.
Thorburn, P.J., F.J. Cook, and K.L. Bristow. 2003. Soil-dependent wetting from trickle emitters: Implications for system design and management. Irrig. Sci. 22:121–127.
Toride, N., and F.J. Leij. 1996. Convective–dispersive stream tube
model for field-scale solute transport: I. Moment analysis. Soil Sci.
Soc. Am. J. 60:342–352.
van Genuchten, M.Th. 1980. A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J.
44:892–898.
Ye, M., R. Khaleel, and T.C.J. Yeh. 2005. Stochastic analysis of moisture plume dynamics of a field injection experiment. Water Resour.
Res. 41(3):W03013, doi:10.1029/2004WR003735.
Yeh, T.C.J., M. Ye, and R. Khaleel. 2005. Estimation of effective
unsaturated hydraulic conductivity tensor using spatial momentsof observed moisture plume. Water Resour. Res. 41(3):W03014,
doi:10.1029/2004WR003736.
Zur, B. 1996. Wetted soil volume as a design objective in trickle irrigation. Irrig. Sci. 16:101–105.