Analysis of Water Flow under Trickle Irrigation:
I. Theory and Numerical Solution
F. Lafolie, R. Guennelon, and M. Th. van Genuchten*
ABSTRACT
Estimating the size of the ponded surface area is a major sourceof error in modeling water movement under trickle irrigation. After
reviewing available approaches and models for simulating water flow
under trickle irrigation, we present a new numerical solution of thegoverning partial differential equations that improves the predictionof the size of the ponded area at the soil surface. A common formulation could be used for both Cartesian and axicylindrical flowgeometries. The spatial operator was approximated using finite differences, and time integration was performed with an implicit
method. Nonlinearity was treated by means of a fixed point iterativealgorithm. The proposed model is applicable to saturated-unsaturated flow in layered soils having anisotropic hydraulic properties.
Accuracy and stability of the solution are checked and comparisons
made with previous numerical solutions for several trickle irrigationand other unsaturated flow problems.

NCREASED IRRIGATION EFFICIENCY, reduced water
Itilizers
stress, ease of automation, the ability to apply ferand other chemicals with irrigation water, andlow economic cost are some of the factors that havecontributed to the increased popularity of trickle irrigation over the last few decades (Bucks et al., 1982).
Trickle irrigation is now applied to a wide variety ofcrops grown on soils with a broad range in soil texture,soil structure, and hydraulic properties. Notwithstanding its popularity, surprisingly little attention has
been paid to estimating soil water distributions during
trickle irrigation under realistic field conditions. Lackof understanding of how soil water content distributions are affected by the unsaturated soil hydraulicproperties has sometimes resulted in suboptimal management and low water-use efficiency (Hillel, 1985).
The shape and total volume of the wetted soil region
below a trickle emitter varies widely with irrigationand soil hydraulic parameters. Using the number ofdrippers, the discharge rate, and the irrigation frequency as adjustable parameters, a trickle irrigation
F. Lafolie and R. Guennelon, Station de Science du Sol, INRA,
Domaine St-Paul, B.P. 91, 84140 Montfavet, France; M.Th. van
Genuchten, USDA-ARS, U.S. Salinity Lab., 4500 Glenwood Drive,
Riverside, CA 92501. Received 2 May 1988. Corresponding author.
Published in Soil Sci. Soc. Am. J. 53:1310-1318 (1989).

system should be designed such that the wetted soilvolume matches as closely as possible the crop rooting
pattern.
Soil hydraulic properties usually exhibit considerable spatial variability at the field scale. An importanttask in the design of a drip irrigation network is toevaluate how spatial variability affects water content
distributions in the field. Particularly useful would bea sensitivity analysis which shows how deterministic
soil water contents are affected by variations in relevant model parameters. One application of modelingis to find an optimal combination of irrigation andhydraulic parameters that minimizes variability problems and results in field-scale water content distributions which are as homogeneous as possible.
While much has been learned about saturated-unsaturated flow in soils under drip irrigation, a numberof areas in need of further investigation remain. As
pointed out by Hillel (1985), these include developing
realistic methods for predicting temporal and spatial
variations of soil water under drip irrigation for different crop, weather, and soil conditions (including
vertically and horizontally heterogeneous soil); determining the minimal and optimal fraction of soil volume needed for various crops; assessing and controlling downward seepage and leaching rates from theroot zone under drip irrigation; and adjusting water
discharge rates to the soil's infiltration capacity so as
to aid the infiltration of water under the drippers infine-textured or crusted soils, particularly on sloping
fields. While field experiments are important for understanding the physics of flow during trickle irrigation, the above research problems cannot be solved
by experimentation only. Simulation models can help
greatly in the design and analysis of field experiments,
and in determining the most important processes andproperties affecting the performance of a trickle irrigation system.
Mathematically, water flow under drip irrigation isessentially a three-dimensional saturated-unsaturated
flow problem with a moving boundary separating theponded and unsaturated areas of the soil surface. Determination of this transient boundary condition is amajor difficulty in accurately modeling drip irrigation.
Hence, the main objective of this paper is to present
a numerical solution of the flow equations with a rig-

1311

LAFOLIE ET AL.: ANALYSIS OF WATER FLOW UNDER TRICKLE IRRIGATION. I.

orous formulation and solution of the soil surface
boundary problem.
REVIEW OF PREVIOUS WORK
As compared with a relatively large number of published
one-dimensional unsaturated flow models, few studies thus
far have focused on trickle irrigation. Most models havebeen analytical or quasi-analytical, with a lesser number ofnumerical models. Analytical models generally involve
steady-state or transient line-source or point-source infiltration problems (Philip, 1971; Raats, 1971, 1972; Warrick,
1974; Lomen and Warrick, 1974; Warrick and Lomen, 1976;
Ben-Asheret al., 1978; Warrick, 1985). Wooding (1968) also
gave a solution for steady-state flow from a pond, while
Lockington et al. (1984) proposed a solution for point-source
infiltration with ponding, but without gravity. Most analytical solutions have only limited applicability because of various assumptions, including uniform initial water contents,
simplified boundary conditions, nonponding surface conditions, soil homogeneity, steady-state flow, and/or exponential relationships between the pressure head and the hydraulic conductivity to enable linearization of the Kirchhoff
equation.
The first and perhaps most quoted numerical solution to
the three-dimensional trickle irrigation problem was given
by Brandt et al. (1971). Their finite difference model has
been widely used to simulate various laboratory and field
experiments. For example, Bresler et al. (1971) found good
agreement between results obtained with this model and experimental data for two soils, although some discrepancies
were observed at relatively high trickle application rates.
Levin et al. (1979) and Mostaghimi and Mitchell (1983) used
the same model to study the effects of pulsed trickle irrigations on water content distributions in two different soils.
Finite difference (CSMP) solutions were used by van der
Ploeg and Benecke (1974), Ragab et al. (1984) and Fletcher
Armstrong and Wilson (1983). These solutions did not significantly improve the initial work of Brandt et al. (1971);
treatment of the surface boundary condition was in particular less accurate. Ababou (1981) later published a finite
difference model similar to the one by Brandt et al., butallowing also for root water uptake. He obtained relatively
good agreement with experimental data for a sand, but lessaccurate results for a loamy soil.
Brandt et al. (1971), Ababou (1981) and Ragab et al.
(1984) applied Kirchhoff transformations to linearize thenonlinear Richards' equation. While Brandt et al. used tabulated hydraulic functions after smoothing with a spline
function, Ababou and Ragab et al. assumed exponential relationships between the water content and the hydraulic conductivity so that an explicit relationship between the Kirchhoff variable and the pressure head could be obtained. The
models of van der Ploeg and Benecke (1974) and Fletcher
Armstrong and Wilson (1983) used Richards' pressure-headbased flow equation. We note that the Kirchhoff transformation is not easily applied to layered soil profiles becauseof the discontinuities at the interfaces.
One major problem with most numerical solutions published thus far is an inaccurate treatment of the boundary
condition imposed at the soil surface. This problem may
explain at least some of the discrepancies which have been
observed between predictions based on previous models andexperimental data for soil water contents and saturated areasat the soil surface. In this paper, we present a more rigorous
solution of the soil surface boundary wetting problem.
THEORETICAL
Governing Flow Equations
In our analysis, we assume that Darcy's law applies tosaturated-unsaturated flow in layered soils, that the soil is

nonswelling, and that hysteresis in the soil hydraulic properties can be neglected. We allow the conductivity tensor to
be anisotropic in the vertical (z) and horizontal (x,y) coordinates. However, the conductivity tensor is assumed to be
isotropic horizontally. Combining Darcy's law and the continuity equation leads, then, to[1]

in which \p is the soil-water pressure head (L); H = 1$ — z
is the total head (L); C is the soil-water capacity (L- ), approximated by the slope of the soil-water retention curve,0(^); 6 is the volumetric water content; Kh and Ky are thehorizontal
and vertical components of the conductivity tensor (LT~}), respectively; t is time (T); z is distance (L), considered to be positive downwards; and x and y are coordinates in the horizontal plane (L). Rather than using thepressure head, \l/, we found it more convenient to work withthe total head, H. Using H = \j/ - z and noting that d\p/dt= dH/dt, Eq. [1] becomes[2]

Trickle irrigation may lead to two possible flow geometries: plane symmetry for a line-source system, and axicylindrical symmetry for an isolated point source. Axicylindrical symmetry holds as long as there is no interaction
between neighboring point sources, while plane symmetry
may always be invoked between two neighboring line
sources. The flow equation for plane symmetry reduces towith (x,z) f fi
to[0,X] X [0,Z], and for cylindrical symmetry

,

h,

r

dr

[4]

with (r,z) < [0,R] X [0,Z] in which r is the radial coordinate(Fig. 1). Notice that Eq. [3] and [4] are defined on only halfof the real domain. Equation [4] can be simplified by thesubstitutionsc(t,z,r) = rC(t,z)

M</<,z,r) = rK^,z) . [5a,b]

Which leads to the following equation for cylindrical symmetry

Equation [3] and [6] are now both of the typeand differ only because of the different forms of the "effective" soil water capacity and hydraulic conductivity coeffi-

1312

SOIL SCI. SOC. AM. J., VOL. 53, SEPTEMBER-OCTOBER 1989

cients C and K. Hence, the same numerical algorithm anddiscretization can be used for both geometries.
Initial and Boundary Conditions
The general initial condition is- z.
H(r,z,t)\t=0 = ifo[7]
Because of symmetry (Fig. 1) we have a no-flow condition
on boundary r,, i.e.,

A no-flow condition also applies to boundary F2mWhile Eq. [9] is valid without restriction for plane symmetry, the condition may not hold for cylindrical geometry
if changes occur along that boundary during the calculations.
The problem becomes then fully three-dimensional and theassumption of having axicylindrical symmetry can no longer
be invoked to reduce the three-dimensional problem to theset of equations given above.
Three different conditions may be imposed along the bottom boundary (F3) of the soil profile. If a water table ispresent at a depth Z, then a Dirichlet type boundary condition applies[10]
H(r,Z,t) = -ZA free-draining profile results when the condition= 0 r—I-1[11]°r r3z'2 = :is imposed. Finally the presence of an impervious layer at
Z is simulated with a no-flow boundary conditionz[12]

At the soil surface, we must account for interactions between the applied trickle rate, evaporation, and infiltration
HORIZONTAL DISTANCE, r
R s ot(t)

T4,2

r=R

Grid point

Internodal Conductivity

-r,

t
UJ

Q = -Kv(t,r,z]

dzz

[13]

where Q — Q0/S — E, Q0 is the trickle discharge rate, E isthe evaporation rate from the ponded area, and S is a small
area over which Q is evenly distributed. However, when theapplied flow rate exceeds the soil's infiltration capacity, asaturated area develops on the soil surface below the emitter.
This ponded area is assumed to be circular for axicylindricalsymmetry and an infinite strip for plane symmetry. We assume that there is no storage of water on the soil surface
during ponding. The following Dirichlet condition then applies to the saturated zone
[14]
H(r,0,t) = 0
re[0,Rs(t)]where /?,(?) is the radius of the ponded circular area (point
source) or the half-width of the strip (line source). Evaporation from the remaining part of the soil surface is described with a Neuman-type boundary condition of the form
air

~EP
[15]
0
where ^Um is some arbitrary lower limit of \f>, and Ep is thepotential evaporation rate. The main, and significant, problem remaining is to quantify Rs(t). This is discussed next.
Determination of the Ponded Area
Let q(r,t) denote the distribution of vertical fluxes acrossthe saturated zone at the soil surface. We define for both
symmetries Qi«(0 as tne total amount of water crossing thesaturated zone per unit time, i.e.,{= 2ir rq(r,t) dr (cylindrical) [16a]Jo

dT

Trickle Source

into the soil. The following condition holds when the soil
surface remains unsaturated (no ponding)Q/r»(Lt(0 = 2(•R,

q(r,t) dr (Cartesian) .
[16b]
Jo
Evaporation from the saturated area is assumed to occur at
its potential rate, Ep, so that we can write
QJt)+ TrRl(f)EP= Qo (cylindrical)[17a]QJit) + 2Rs(t)Ep = Q0 (Cartesian) .

[17b]

For notational convenience below, we drop the explicit timedependence of Q^ and Rs. Combining Eq. [16] and [17]gives
2ir P rq(r,f) dr = Q0 - irR2sEp (cylindrical)
Jo

[18a]

2 p q(r,t) dr=Q0- 2RSEP (Cartesian)Jo
where q(r,f) is given by Darcy's law as
affq(r,t) = -Kv,

[18b]

The time derivatives of Eq. [18a,b] are.dRs

r;dR*.

=

dt
z=ZFig. 1. Schematic cross section of the flow domain and the imposed

finite difference grid system for simulating trickle irrigation.

[19]

dt
R

+ 2ir

C * r^(r,t)drda(cylindrical)

Jo

dt

[20a]

1313

LAFOLIE ET AL.: ANALYSIS OF WATER FLOW UNDER TRICKLE IRRIGATION. I.

dR,~

more accurate procedure for estimating Rs(t) is described
below.

rJJQo

dr(rt)dr (Cartesian)[20b]

which leads to the differential equations

An Alternative Method for Estimating Rs(t)For any part of o> of the wetted domain (Fig. 2), we canwrite the mass conservation equationf ^ da = - f div(f)fi?o>

Judt

r

r—(r,t)dr (cylindrical)
ot
Jo

dt

[21a]

}„

where f = [K]Grad(W), and [K] is the conductivity tensor.
The divergence theorem states thatfdiv(f) du = I f -n dj[Ep + q(Rs,t)]

J <a^(r,t)dr (Cartesian). [21 b]
Jo-/'
Equation [21a,b] shows that a decrease in the flux throughthe saturated zone (dq/dt < 0) results in an increase of thesaturated area at the surface (dRJdt > 0). The initial condition for Eq. [21a,b] is the radius R,# of the saturated area
corresponding to the impact of a single drop. This area isnumerically approximated by the area S corresponding to
the first node of the finite difference grid. If we denote tp as
the time at which ponding occurs, the initial condition becomesp it \ — i>mi
•KsUp.) ~~ •'Ho •l"JThus, we obtain a set of equations describing saturatedunsaturated flow in the soil, and a differential equation predicting the evolution of Rs as a function of time. The resulting Stefan-type problem can be solved using an iterativealgorithm (Mitchell and Griffith, 1980) involving alternate
solutions of the flow equation, Eq. [6], and Eq. [21a,b] for
R, until convergence is reached. Two problems remain.
First, accurate imposition of a Dirichlet boundary conditionon [0,.RS(0] requires modification of the finite difference grid
system each time a new estimate for Rs(t) is obtained. Thisimplies a need to interpolate the potential field each time a
new grid is defined, with the likelihood of introducing additional errors. Also, a good approximation of the functionq(r,t) over the interval [0,Rs(t)\ is crucial for accurately solving Eq. [21]. The standard method for estimating water
fluxes through the saturated zone is to use Eq. [19] in which
the gradient (dH/dz) is approximated with a two- or threepoint uncentered scheme (Bresler et al., 1971; van der Ploeg
and Benecke, 1974; Ababou, 1981; Ragab et al., 1984) according to

123]
Az" "
Az
A7/ = -3//t + 47/2 [24]Az
2Az
where Hlt H2, and H3 are the total head values at the first
three finite difference nodes along a vertical starting fromthe soil surface.
Numerical experiments using two types of soils were carried out to investigate the accuracy of Eq. [24a,b] in approximating q(r,t). Unacceptable errors appeared, especiallywith relatively large step sizes, Az, and/or low initial watercontents. These results may explain why, for coarse-textured
soils where small ponded areas and very small gradientsnear the surface prevail, most models give relatively goodresults, and why, for fine-textured soils, they fail to accurately predict the ponded area on the soil surface. Hence,previous models often led to inaccurate water-content distributions in the profile, and also exhibited relatively large
mass balance errors when used on fine-textured soils. A

[25]

[26]

J d<a

where n is the outward direction normal to 3<o. Equation[26] can be expanded to give the following equation relatingwater content variations in &> to the various boundary fluxes

cde ,

r 7-

\-rdw = - f -ndy
JadtJdu
- - f f-nrfy + f f-krfy + f f-Trfy
JdF
JH
Jv

[27]

where k and i are unit vectors for the z and r coordinates,
respectively; H = [Q,RS], V=[Q,Z{] and dF is an arbitrary
curve inside n joining the points (Rs,0) and 0,Zf). The last
integral of Eq. [27] is zero because of a no-flow boundary
condition on r\. Rearranging Eq. [27] gives thenT'f'k dx=-

I 1 -n dy .

[28]

JHF
Jo
Hence, the discharge rate through the saturated zone (lefthand side of Eq. [28]) can be calculated from the rate ofchange in total water content in w, and the flux through someinternal boundary 8F. For cylindrical symmetry, this boundary is generated by rotating 8F around the axis (0,ik), andfor plane symmetry by translating 8F along a direction orthogonal to the plane (0,i,k). Note that when 6 in o> does
not change with time (steady flow in a), the discharge rates
through the saturated zone and dF must be equal.
We required dF to be such that the integration of 6 overto, as well as the estimation of the fluxes through dF, would
be as accurate as possible. For convenience, dF was denned
TrickleSourceGrid point

Fig._2. Approximation of the smooth curve, dF, by a broken line,
df, along the finite difference grid system.

1314

SOIL SCI. SOC. AM. J., VOL. 53, SEPTEMBER-OCTOBER 1989

such that the pressure head in w was greater than some preset
limit cltise to saturation (e.g., a few centimeters pressurehead), thus guaranteeing a smoothly changing potential field.
Numerical implementation of this algorithm was facilitatedby approximating the original smooth boundary, dF, with abroken line, dF, that follows the grid system as shown in
Fig. 2. This approximation does not lead to additional errors
since the divergence theorem holds for any piecewise continuous bdundary. Our numerical algorithm automaticallyconstructs the curve dF, calculates the required fluxes, andevaluates the total moisture change in w.
The computer algorithm to implement the above scheme
is described schematically in Table 1. In this table, /„, indicates the set of nodes belonging to w, while the expression"calculate the horizontal flux through (ip jj)" means to calculate the flux from node (if,.//) to node (/} + 1,y}). Similarly,"calculate the vertical flux through (if,j$ means to calculate
the flux from node (if,jj) to node O/, y/ + 1). Two-point centered finite difference equations were used to calculate all
fluxes through dF.
The accuracy of the method for both symmetries was verified by simulating two-dimensional infiltration through afixed wetted zone using several application rates. The accuracy during each time step (A/) was expressed by the relative mass balance variable

[29]

v,

where V, is the amount of water applied to the soil surface
during A?, and
80 = I dd(x,y) dxdy
Jv>

[31]

JiF

Table 1. Algorithm for calculating the curve dfand the flux throughthe unsaturated zone at the soil surface._____________________________________________

a. Initialize (if,jf) = (4,1) where (4) is the first unsaturated node at the topof the soil (see Fig. 2).
b. Calculate the horizontal flux between the nodes (if — \Jf) and (ifjf).
c. Let (ifjf) = (irjf + I),
1. If (ifjf) t /,„ thena) iel(if,jf) = ( i f - IJf)p) l f ( i f < 0)stop.
7) If(i>7»«/ M then—calculate the vertical flux through (if + ljf — I),—calculate the horizontal flux through (ifjf),—go to c.
5) If (ifjf) tf /„, or if (if = 0) then—calculate the vertical flux through (if + \Jf — I),—let (ifjf) = (ifjf — 1),—go to c.

where Rs'(t) is the current estimate ofRs(t), and Qs(t) is theflux through the saturated area corresponding to Rs(t). If weapply boundary condition Eq. [14] to the saturated area
[0,Rs'(t)] rather than to [0,R,(t)], a difference, D(t), between
the amount of water applied and the amount entering thesoil profile will result. This difference is given by D(i) = Q0
— Qs(0- To account for D(t) in the model, a flux boundary
condition is imposed on the first unsaturated node, /'„ of thegrid system (Fig. 2). This gives

Vr = I Q-n dy

in which d0 is the local water content change during A/, andQ is the flux through dF as given by Darcy's law. The variables SO and VT represent the change in water storage in toand the flux through the boundary OF, respectively. For several soils and application rates, we always found the ratio
Rm to be between 0.995 and 1.005. Cumulative mass balances errors during the simulations generally also remained
less than 0.5%.
Discharge rates through the saturated zone calculated inthis fashion take advantage of the accuracy of space-centeredapproximations in a region where the pressure head always
varies smoothly. Also, the saturated hydraulic conductivity(Ks) is not needed to calculate hydraulic fluxes through theponded area on the soil surface. Actually, fluxes are calculated using conductivities, which often are one to three orders of magnitude smaller than Ks, thus further reducing
errors when gradients are approximated by discrete equations.
Some authors (e.g., Bresler et al., 1971; Ragab et al., 1984)
have pointed out that poor estimates of Ks could have

Algorithm

caused their deviations between calculated and experimental
data. The As is known to be a highly variable parameter inthe field, with variations of two orders of magnitude or more
not being uncommon in the same field (e.g., Russo and Bresler, 1980). Hence, if one is interested in evaluating the effectsof spatial variability in Ks on soil water distributions, thena model should be used in which numerical errors are notvery sensitive to this parameter. Previous models may wellhave induced large errors in this respect. Because the hydraulic conductivity at small negative values of ty is often
an order of magnitude less than at saturation, and usuallyalso exhibits less variability at the field scale, our approachshould be more reliable than previous models.
Finally, we note that it is not necessary to obtain Rs(t) as
a continuous function of t. It is much simpler and more
economical to use a step function Rs'(t) which approaches
Rs(t) such that Rs'(t) always coincides with a finite difference
nodal point at the soil surface. For cylindrical geometry,
Rs(t) is then given bywhich ensures a correct mass balance. In Eq. [31], QT(t) is
the residual flux applied at surface node is, and given by
QT(t) = D(t)/Si?, where S^ is the soil surface area between
nodes *s_, and is of the finite difference grid. As the infiltration rate through the saturated zone decreases, the residual
flux, Qr(t), increases until saturation is reached, at whichtime Rs'(f) jumps from node i, — 1 to node /,.
NUMERICAL APPROXIMATION OF THE
FLOW EQUATION
The flow domain (Fig. 1) is divided into an unequally spaced, rectangular finite difference grid fih =(r,, ... ,rnx) X (zi, ... ,znz) with arbitrary node (ij)located at (r,,z,). Another set of points in Q is denned

by
[32a]

X

such that

2,nx
2,nz.

[32b]
[32c]

To more easily include no-flow boundary conditions, we choose r{ and rnx such that r{ = 0 andrnx+\ = R- F°r the horizontal boundaries, T3 and T4,
we have z, = 0 and znz = Z so that z{ = — z2 andznz+i = 2znz — znz. The continuous function H (r,z,t)
is approximated on the finite difference grid by thetime-dependent vector [Ht>J].
Approximation of Spatial Derivatives

At any point (ij), spatial derivatives are approximated with three-point centered finite difference equations as follows

1315

LAFOLIE ET AL.: ANALYSIS OF WATER FLOW UNDER TRICKLE IRRIGATION. I.

- Hu

dzvdz

— Zi

[33]

z

i ~

which yield a spatial approximation of quadratic order(Godounov, 1973). Internodal conductivities, kh andfcv, are calculated by geometric means (Vauclin et al.,1979), to give

Neuman-type boundary conditions were readily included by two-point centered approximations, while
Dirichlet-type nodes were directly eliminated from thefinal matrix equations. Flux conditions on T3 and T4
required virtual nodes outside the domain; thesenodes give contributions that cancel with approximations from inside the domain. Combining the different spatial approximations ultimately leads to thefollowing set of nonlinear differential equations(C(H)]= ~[K(H)] [H] + {S(H)}[34]

in which [C] is the diagonal matrix of nodal soil watercapacities, [K(H)] is a multidiagonal matrix associatedwith the spatial derivatives, {S(H)} is a constraint vector accounting for Neuman-type boundary conditions,and [dH/dt] is the vector of temporal pressure derivatives at the grid nodes.
Integration in Time
Equation [34] was integrated by means of a onestep A(a)-stable method (Richtmeyer and Morton,
1967). For a differential equation of the form y'(t) =f(y,t) this results in the following scheme
yM = yj + &(af(yj+l,tj+l) + (1 - ctffytj)] [35]where the indices j and 7+1 refer to the old and newtime levels, respectively. For a equal to 1/2 and 1, thewell-known Crank-Nicholson and fully implicit
schemes result. Applying Eq. [35] to Eq. [34] leads to+ [[7] - (1 - a)A/[C-i]/[AK) [H}>[36]
where [/] is the identity matrix. The nonlinearity ofthis set of equations is treated by means of a fixedpoint algorithm. The iterative process continues until
the following criterion is satisfied
IIWi - WIU < «[37]where k is the fixed-point iteration index and t. anerror limit. Numerical tests showed that a value of 0.1
cm for « was adequate; results generally changed little
for e-values ranging from 0.1 to 0.001. The above al-

gorithm was also used to control the time step A? during the calculations. The A? was increased if the number of iterations during the preceding step was lessthan three, remained the same if the number of iterations was three, four or five, and was reduced formore than five iterations. The final set of linear equations is of relatively large size, and direct inversion
methods can be quite time consuming. Because of theblock-diagonal properties of the matrix [K\, the system
is well suited for iterative equation solvers such as theGauss-Seidel or SOR algorithms. We used the GaussSeidel
5 method with a relative stopping criterion ofio- .

NUMERICAL VERIFICATION

Convergence and stability of the numerical solution
were first studied by using increasingly smaller time
(A?) and spatial (Ar,Az) increments and assuming afully implicit integration in time (a = 1). By letting
Az and Ar go to zero, simulations showed that convergence to the exact solution was quadratic everywhere in the domain, provided that Az and Ar near
the wetting front were relatively small. As expected,
the fully implicit scheme showed a linear rate of convergence with A;. We successfully used a Richardsontype extrapolation to the limit (Fairweather, 1978)
with a geometric series of time increments: Af/2', / =0, 1 ... Thus, we concluded that the numerical scheme
was stable and convergent in space and time.
Partial verification of the numerical scheme could
have also been carried out against available analyticalsolutions for point sources (Warrick, 1974; Lomen
and Warrick, 1 974), or strip and disc sources (Warrickand Lomen, 1976). Unfortunately, these solutions arefor the quasi-linear Kirchhoff equationdt

dz

in which U is the matric flux potential, and k and 0
are constants such that K=KS exp (/%<) and dK/dB =k. Because linearized solutions do not provide a convincing numerical check of the complete nonlinear
model, we decided to compare our results with previously published solutions on problems involving
both Dirichlet (Example 1) and flux-type (Example 2)
boundary conditions. Several trickle irrigation problems are discussed under Example 3.
Example 1. To verify the model for Dirichlet-typeconditions, we simulated the problem of infiltrationfrom an equally spaced furrow irrigation system. The
flow domain and boundary conditions are shown in
Fig. 3. The problem was solved previously by Zyvoloski et al. (1976) using a finite element method, by
Caussade et al. (1973) using an iterative linear/nonlinear finite difference approximation after applicationof Kirchhoffs transformation, and by Selim and Kirkham (1973) and Firdaouss and Ta-Phuc-Loc (1983)
using alternating-direction implicit (ADI) methods
with and without iteration, respectively. We refer tothe studies of Zyvoloski et al. (1976) or Caussade etal. (1979) for detailed descriptions of the hydraulicproperties. In our calculations, we used a regular finite
difference grid with vertical and horizontal nodal
spacings of 2 cm. Figure 4 compares our calculated
water content distribution after 50 min of infiltration

1316

SOIL SCI. SOC. AM. J., VOL. 53, SEPTEMBER-OCTOBER 1989
10

20

30

• rfcm)

dr =0

LJ
g

I.U4

Sand DEK
Line Source< 1.0203-pQp-dSX-fln O D O O O tP

15
5

25-

0.98 Q

UU

Q

(

~

o Q0=|. 5 l-h" rri'

p 0.96 •o Q0=30 Itf'ni1

LJ

8

35-

cP nQ-m _ - _'———_l 0.94

10;. . , , i .
20

15

TIME (h)Fig. 5. Calculated mass balance errors during line-source trickle
irrigation on DEK sand for two irrigation rates.

45-

Z
(cm)

1.20

Fig. 3. Schematic cross section and boundary conditions for the
furrow irrigation example.

Loess MC4Line Source

LJ

"1.16
ti

20

30

»rfcm)-

2oco 1.08

Grid Size:

21.04

IS

25-

D: Az = I.O cm
O: &z = 2.0cm

UJ

>1.00

OD OD o OD

J0.96

o:
0.92

6

12

18

24

30

36

42

48

TIME (h)Fig. 6. Calculated mass balance errors during line-source trickle
irrigation on a loess soil for two finite difference nodal spacings.

— —-- Zyvoloski——— Selim———— This Study
Fig. 4. Comparison of calculated water content distributions for the
furrow irrigation problem.

with the results of Zyvoloski et al. (1976) and Selimand Kirkham (1973). Clearly, our results agree closelywith those of Zyvoloski et al., but differ from the ADI
solution of Selim and Kirkham.
Example 2. The time dependent behavior of the
global mass balance often gives a good estimate of theoverall accuracy of a numerical solution. Figure 5 and6 show calculated mass balances for trickle irrigationfrom a line source using the same rectangular geometry as in Fig. 1. Results in Fig. 5 are for two irrigation
rates involving a sandy soil, and assuming a regulargrid system with constant nodal spacings of 2 cm. Figure 6 gives similar results for a loess soil, assuming
grid spacings of 1 and 2 cm. These examples werepreviously considered by Ababou (1981), who also
gave a detailed description of the unsaturated hydraulic properties of the two soils. No surface ponding
occurred for both soil types. Figure 5 shows that massbalance errors for the DEK sand of Ababou becameless than about 0.5% after only a few minutes of infiltration. Much larger errors occurred with the MC4loess soil (Fig. 6), which had a very low saturated hydraulic conductivity; errors in this case remained

about 2% during the simulation. Improved accuracy
in this case was only possible by decreasing the nodal
spacing. The errors shown in Fig. 5 and 6 were only
about one-fourth of those obtained by Ababou (1981).
Example 3. We briefly present two other cases fortrickle irrigation. The first case was carried out for thesame loess soil as in Example 2. We again compared
our results with those obtained by Ababou (1981) using an ADI-based model. The following experiments
were simulated: a line1 sourcewith discharge rates of
0.5, 1.5, and 3.0 L h" nr1, and a point1 source withirrigation rates of 0.5, 1.5, and 3.0 L h" . An irregular
finite difference grid was now used for the calculations;the vertical nodal spacing increased from 2 cm at thesoil surface to 4 cm at the lower boundary, while aconstant spacing of 2 cm was used in the horizontal
direction. The calculated mass balance errors and saturated soil surface areas are summarized in Table 2.
Note the large differences in mass balance between ourscheme and the ADI scheme of Ababou (1981). Figure
7 shows a typical example of the type of water content
distributions obtained for a line source. As for theponded infiltration case (Example 1), the moisture
front calculated with the ADI model is located at ashallower depth than those obtained with our finite
difference solution.
Model results were also compared with those obtained by Brandt et al. (1971) for the more fine-textured Gilat loam, which was expected to show muchlarger discrepancies compared with their Nahal Sinai
sand. Simulations were carried out
for two trickle discharge rates (0.495 and 0.983 cm3 crrr1 min-') as used

1317

LAFOLIE ET AL.: ANALYSIS OF WATER FLOW UNDER TRICKLE IRRIGATION. I.
Table 2. Calculated mass balances and saturated soil surface areas
for three trickle irrigation rates involving a line source and a point
source. Results of this study are compared with those obtained
by Ababou (1981) using an ADI scheme.
Line source
discharge rate

Mass balance %this study
Mass balance %Ababou (1981)

1.5

3.0

<1.0 <1.0 <1.0

10

20

30

40

Point source
discharge rate

Go (L h-' m-)

0.5

HORIZONTAL DISTANCE (cm)

Co (L h-)1.5

3.0

<1.0 <1.0

<1.0

0.5

-3.5

-5.6

-9.8

-3.0

-6.2

-8.6

K,

2.0

23.0

47.5

18.0

39.0

58.0

R,

1.5

20.0

46.0

16.0

36.0

53.0

this study, cm
Ababou (1981), cm

HORIZONTAL DISTANCE, r(cm)
O

10

20 _

30_

40

Fig. 8. Calculated water content distributions obtained in this study(solid lines) and by Brandt et al., 1971 (dashed lines) for linesource trickle irrigation on Gilat loam.

Fig. 7. Calculated water content distributions resulting from linesource trickle irrigation on a loess soil.

by Brandt et al., and assuming plane symmetry. The
same irregular finite difference grid as before was applied. We obtained much smaller saturated zones atthe soil surface, especially for the higher trickle irrigation rate (Fig. 8). After applying 3 L of water, themodel of Brandt et al. yielded saturated areas corresponding to Rs values of 11 and 24 cm for the twodischarge rates, while we obtained R^ values of 3 and9.5 cm, respectively. Water content distributions in
the profile were also different, with our calculations
showing deeper infiltration of water. Similar but
smaller differences between the ADI and finite difference/finite element solutions were already noted in
Example 1, as well as in the preceding case involving
trickle irrigation.
SUMMARY AND CONCLUSIONSAn improved numerical model for simulating saturated-unsaturated water flow in general and trickle
irrigation in particular has been described. The timedependent nature of the interface between the saturated and unsaturated zones at the soil surface wassimulated by means of a moving boundary-type dif-ferential equation. This equation was coupled with theunsaturated flow equation through the water flux
across the saturated zone at the soil surface. A newand computationally efficient algorithm was used to
avoid the calculation of fluxes across the saturated
surface area by uncentered finite difference equations.
The numerical solution was successfully tested against
two other numerical finite difference/finite element solutions for furrow irrigation subject to Dirichletboundary conditions, and for trickle irrigation involving Neuman-type conditions at the soil surface.
We also compared our results with those obtained
by Ababou (1981) and Brandt et al. (1971). Comparisons for trickle irrigation involving surface ponding
were carried out using various soil hydraulic data sets.
Our algorithm for calculating the saturated area at thesoil surface produced excellent mass balances as compared with these previous solutions. The algorithm
also resulted in improved estimates of the saturated
surface area, and in improved predictions of the soil
water content in the profile. The numerical solution
will be used in Part 2 of this study (Lafolie et al., 1989)
to analyze a set of field experiments involving trickle
irrigation.

1318

SOIL SCI. SOC. AM. J., VOL. 53, SEPTEMBER-OCTOBER 1989