System-Dependent Boundary Condition for Water Flow from Subsurface Source
Reproduced from Soil Science Society of America Journal. Published by Soil Science Society of America. All copyrights reserved.

N. Lazarovitch, J. Šimůnek, and U. Shani*
ABSTRACT

also restricted to the Gardner’s (1958) hydraulic conductivity function, which represents a significant simplification of the dependence of the hydraulic conductivity on
the pressure head. Gardner’s model reasonably predicts
soil hydraulic properties at high water contents but fails
to capture the behavior at low water contents (Or etal., 2000). Since numerical models, on the other hand,
can implement any model of the soil hydraulic properties and can provide solutions for transient conditions,they can overcome many restrictions imposed by analytical models.
The objectives of this work were (i) to implement
into a transient numerical model a system-dependent
boundary condition that describes infiltration from subsurface sources and considers the source characteristic
curves, (ii) to evaluate the dependence of the transient
source discharge on the soil and source hydraulic properties, (iii) to compare numerical results with existing
steady-state analytical solutions, (iv) to carry out laboratory experiments involving transient water flow fromsubsurface water sources, and (v) to compare obtained
experimental results with a numerical model.

Infiltration rate of water from a subsurface cavity is affected bymany factors, including the pressure in the cavity, its size and geometry,
and the hydraulic properties of the surrounding soil. When a predetermined discharge of a subsurface source (e.g., a subsurface emitter) islarger than the soil infiltration capacity, the pressure head in the source
outlet increases and becomes positive. The built up pressure may significantly reduce the source discharge rate. The main objective of thiswork was to develop a boundary condition that describes this process
while considering the source characteristics, and to implement thisboundary condition into the transient numerical model Hydrus-2D.
This new, system-dependent boundary condition allows calculationof the source discharge while considering source properties, inlet pressure, and effects of the soil hydraulic properties. The updated numerical model was verified against existing analytical solutions for simplifiedsteady-state conditions and validated against transient experimental
data. Good agreement was found between transient cavity pressures
measured in laboratory experiments and those calculated using the
updated numerical model. The modified program allows using any
analytical model that describes the soil hydraulic properties and source
characteristics, simulating both short and long duration infiltration
events, as well as considering various geometrical shapes of subsurface sources.

THEORYW

hen a predetermined discharge of a subsurface
source (e.g., a subsurface emitter) is larger than
the soil infiltration capacity, the pressure head in thesource outlet increases and becomes positive (Philip, 1992;
Shani et al., 1996). This pressure buildup in the soil reduces the pressure difference across the source orifice
and, subsequently, decreases the source discharge rate
(Warrick and Shani, 1996). It was previously shown thatthe degree of the pressure difference reduction is larger
for soils having lower hydraulic conductivity and thatthe positive pressure in the source vicinity increases
rapidly at the beginning of the infiltration event and approaches a final value after only several minutes (Shaniet al., 1996; Warrick and Shani, 1996). The source discharge is also affected by the source characteristics andthe cavity size of the source outlet (Shani et al., 1996).
Analytical solutions for water flow from a subsurface
cavity are available for calculating spatial distributionof the pressure head in the soil, ␺ (L). However, sincethese solutions were derived only for steady-state waterflow conditions (Philip, 1992; Warrick, 1993a, 1993b;
Warrick and Or, 2001), the use of these models is limited
to relatively long infiltration events. These models are

Warrick and Shani (1996) suggested the following relationship between the source discharge, Q (L3 T ⫺1), and the inlet
pressure of the source (Pin [L]):冢P P⫺ ␾ 冣

Q ⫽ Qo

c

in

s

[1]

in

where φs (L) is the hydraulic head at the source-soil interface,often called the back pressure, Q0 (L3 T ⫺1) is the nominal
discharge of the source for the reference inlet pressure Pin(usually being 10 m) and the back pressure equal to zero, andc (⫺) is an empirical constant that reflects the flow characteristics within the emitter. Normally, c ⫽ 0.5 corresponds to a
turbulent flow emitter and c ⫽ 1 to a laminar one. This equation states that when ␺s increases, the pressure difference
between the soil and the source inlet decreases and the discharge rate decreases correspondingly.
Philip (1992) approximated solution of pressure distribution
near a small spherical subsurface water source by matching
fluxes from the saturated region surrounding the physicallyrequired cavity, with fluxes predicted by his (Philip, 1968)
quasi-linear solution for unsaturated flow from a point source.
Shani and Or (1995) used Philip’s (1992) solution to relate ␺s(L) (the cavity pressure head) to the soil hydraulic properties
and the source discharge:

N. Lazarovitch and U. Shani, Dep. of Soil and Water Sciences, Facultyof Agricultural, Food and Environmental Sciences, The Hebrew University of Jerusalem, P.O. Box 12, Rehovot 76100, Israel; J. Šimůnek,
Dep. of Environmental Sciences, University of California, Riverside,
CA USA. Received 5 Jan. 2004. *Corresponding author (shuri@agri.
huji.ac.il).

␺s ⫽

⫺␣ r
1
Q ⫺
冢28␲K
r 冣
␣
G 0

S 0

∞

[2]

G

where r0 (L) is the cavity (source) radius, Q∞ (L3 T ⫺1) is thesteady state discharge, KS (L T ⫺1) is the saturated hydraulicconductivity and ␣G(L⫺1) is a soil parameter.
Gardner’s (1958) hydraulic model, slightly modified to account also for positive pressure heads was assumed for thesolution presented in Eq. [2]:Published in Soil Sci. Soc. Am. J. 69:46–50 (2005).
© Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

46

LAZAROVITCH ET AL.: SYSTEM-DEPENDENT BOUNDARY CONDITION

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

K(␺) ⫽ min [KS , KS exp(␣G␺)]

47

[3]

where ␺ is the pressure head (L). Calculation of Q∞ is thendone by simultaneous solution of Eq. [1] and [2] (Shani etal., 1996; Warrick and Shani, 1996).
In our analysis we used a numerical solution 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 subsurface irrigation systems(Skaggs et al., 2004). The governing flow equation for axisymmetrical three-dimensional isothermal Darcian flow in a variably saturated isotropic rigid porous medium is given by thefollowing mixed form of the Richards’ equation:

冢

冣

冢 冣

⳵␪ 1 ⳵
⳵␺
⳵ ⳵␺
⳵K
⫽
rK
⫹
K
⫺
⳵t
r ⳵r
⳵r
⳵z ⳵z
⳵z

[4]

where ␪ is the volumetric water content (L3 L⫺3), K is thehydraulic conductivity (L T ⫺1), r is a radial coordinate (L), z
is a vertical coordinate (L) positive downward, and t is time(T). Equation [4] is solved numerically using a finite element
method for given initial and boundary conditions applicable
to a subsurface infiltration experiment.
In addition to existing boundary conditions that are currently available in the Hydrus-2D code, Eq. [1] was implemented as a new system-dependent boundary condition. For
conditions where the source surrounding is homogenous andthe outer boundaries do not affect the flow field near thesource, this boundary condition allows the calculation of thesource discharge while considering properties of the source,
and the inlet and back pressures:冢P P⫺ ␾ 冣 ⫽ ⫺冮 K冢⳵␺⳵r ⫹ ⳵z⳵␺ ⫺ ␦ 冣n d⌫ [5]Q ⫽ QoFig. 1. Finite element grid and flow domain. Water flows from the half
sphere on the left (zoomed). Boundary condition [1] is specified on
the sphere surface, while no water flow condition is used on theother boundaries.

c

in

s

iz

in

i

⌫

where ␦iz is the Kronecker delta (⫺), ⌫ is the boundary of thesource, and ni are the components of the outward unit vector
normal to the boundary. At each time step during calculations,the discharge rate is adjusted while φs is calculated by averaging
pressure head (from the previous time step) over the nodes
that describe the source boundary as follows:␾s ⫽1· 冮␾sd⌫⌫(Fig. 1). The half circle of the source was represented with 12
nodes. Unstructured finite element mesh was generated using
automatic triangulation that is implemented in Hydrus-2D andthat uses an algorithm based on the Delaunay’s triangulation(Šimůnek et al., 1999). All other boundaries, except for thefree drainage boundary condition at the bottom of the flow
domain and the source boundary condition (Eq. [5]) at thesource-soil interface, had a no flow boundary condition.

Experimental[6]⌫

Note that during initial times φs can be negative. Equations[5] and [6] are valid also in this situation.

MATERIALS AND METHODS
Numerical Calculations
The Mualem–van Genuchten soil hydraulic model (vanGenuchten, 1980) was selected for the numerical simulations:␪ ⫺ ␪rn⫺1[7]⫽ [1 ⫹ (⫺␣␺)n]⫺m; m ⫽␪s ⫺ ␪rnK(␺) ⫽ KS S e0.5[1 ⫺ (1 ⫺ S e1/m)m]2
[8]
Se ⫽

where Se is the effective fluid saturation (⫺), ␪r and ␪s denote
the residual and saturated water contents (L3 L⫺3), respectively; and ␣(L⫺1), n (⫺), and m (⫽ 1 ⫺ 1/n ) (⫺) are empiricalshape parameters.
A flow domain was selected such that the outer boundaries
do not affect the flow field near the cavity. The flow domain(50 ⫻ 30 cm) was discretized into 1183 triangular finite elements with triangles significantly smaller around the source
and then smoothly increasing with distance from the source

Non-compensated drippers were alternatively packed inside perforated plastic spheres (the source) with radii of 0.0025
and 0.01 m. The purpose of the perforated plastic sphere wasto maintain a rigid cavity. The opening diameter was about
0.005 m, and included ⬎80% of the sphere surface area. The
cavity radii were small enough so we assumed φs ⫽ ␺s forcomparison with the analytical Eq. [2] predictions. The source
was connected to a flexible tube and buried 0.25 m deep in
two soils (Arava sandy loam, coarse-loamy calcareous Typic
Torrifluvent; and Magal clay loam, from the Netafim Experimental Farm, Magal, Israel), with the soil around it subsequently packed. Sources with nominal discharges Q0 of 2, 4,
and 8 L h⫺1 were used, while the value of the c exponent in
Eq. [1] was the same for all sources and equal to 0.5. Hereafter
we denote the source by its Q0 value.
The saturated hydraulic conductivity was estimated using
the horizontal soil column constant head method (Hillel, 1971,
p. 82–83, 106–107, 228–230). The retention curve ␪(␺) was
measured using a ceramic plate suction cell apparatus (Klute,1986). Parameters of the hydraulic conductivity function K(␺)
for the Mualem–van Genuchten model were estimated using
the least-square optimization technique of the RETC program(van Genuchten et al., 1991). The ␣G parameter for the Gardner model was estimated by equating the Kirchhoff potential,
F, (Gardner, 1958):

48

SOIL SCI. SOC. AM. J., VOL. 69, JANUARY–FEBRUARY 2005

Table 1. Textural properties and hydraulic parameters of studied soils.

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

Site
Soil class
Sand, %
Silt, %
Clay, %␣, m⫺1n (⫺)
K s m s⫺ 1
␣ G , m⫺ 1
␪r
␪s

Arava

Magal

sandy loam
83
8
9
6.1
1.495
3.7 ⫻ 10⫺5
19.2
0.04
0.36

clay loam
35
20
45
3.01
1.57
3.47 ⫻ 10⫺6
7.9
0.05
0.51

0

F ⫽ 冮 K(␺)d(␺)⫺∞[9]of the two hydraulic models and substituting ␣G ⫽ KS ⫻ F ⫺1.
The integration of Eq. [9] for the Mualem–van Genuchten
model Eq. [8] was executed numerically. Hydraulic and textural properties of the two soils used are summarized inTable 1. In our calculations below, we assumed that we could
neglect hysteresis and use retention curve determined using
the drying process for the infiltration process.
The subsurface source was connected to a Mariotte apparatus that maintained a constant input pressure, Pin. Discharge Q
was measured by weighing the Mariotte apparatus. A pressuretransducer was inserted into the source to measure the pressure head in the source cavity. The initial water content wasmeasured at several points and depths using time domainreflectometry (TDR) (TDR-100, Campbell Scientific, Logan,
UT). The average initial water content was slightly higher than
the residual water content. All measurements were performed
until discharge or the back pressure ␺s changes were smaller
than 0.5% during an interval of 5 min, assuming steady-state
conditions were reached. All sensors, such as scale, pressuretransducers, and TDRs were attached to the data logger(CR10x, Campbell Scientific, Logan, UT). Measurements weretaken every 0.1 s and averaged over a second. Each experiment(infiltration event) was replicated three times.

Fig. 2. Measured (triangles and circles) and calculated (lines) soilpressure heads as a function of time for two subsurface sources (4and 8 L h⫺1) in the Magal clay loam with the source radius r0 of
0.01 m and the inlet pressure Pin of 10 m. Error bars represent ⫾standard deviation.

RESULTS AND DISCUSSION
Measured and calculated back pressures ␺s as a function of time for two sources with the nominal discharge
Q0 of 4 and 8 L h⫺1 are presented in Fig. 2 for the Magal
clay loam with the radius r0 of 0.01 m and the inlet pressure
Pin of 10 m. The initial water content was considered tobe constant throughout the flow domain and equal to0.07. For both sources ␺s increased rapidly during thefirst 10 to 15 min and then gradually approached a final
value as time proceeded. The final value of ␺s was 2
and 3.7 m for the 4 and 8 L h⫺1 sources, respectively.
There was an excellent agreement between measurements and calculations.
The transient decrease of the source discharge is depicted for the Magal clay loam and for the source with
the nominal discharge Q0 of 8 L h⫺1, the inlet pressure
Pin of 10 m, and radius r0 of 0.01 m in Fig. 3. Following ␺sincrease from 0 to 3.7 m (Fig. 2), the pressure difference
between the back pressure and the source inlet decreased, and the discharge rate dropped from 2.22 ⫻10⫺6 to 1.73 ⫻ 10⫺6 m3 s⫺1. The numerical code again
predicted the measurements very well. The results presented in Fig. 3 demonstrate the transient nature of thesubsurface source discharge and the significant decreaseof discharge due to development of the back pressure.
Philip’s (1992) solution for steady-state water flowfrom a subsurface cavity implies a linear relationshipbetween ␺s and Q∞ (see Eq. [2]). The back pressure ␺s
as a function of Q∞ for two soils (Magal clay loam andArava sandy loam) and three sources (2, 4, and 8 Lh⫺1) are presented in Fig. 4. Values of ␺s and Q∞ were
calculated for known Q0, Pin, and r0 analytically using
Eq. [1] and [2] and numerically with Hydrus-2D byassuming the same criteria for approaching steady-stateflow as mentioned above. In both soils the back pressure␺s was higher for the larger nominal discharge Q0. The
back pressure ␺s for the Magal clay loam having lower
KS was larger than ␺s calculated for the Arava sandy
loam soil with the higher KS. For example, the backFig. 3. Measured (triangles) and calculated (lines) discharge from a
subsurface source as a function of time for an 8 L h⫺1 source andthe Magal clay loam with the source radius r0 of 0.01 m and the inlet
pressure Pin of 10 m. Error bars represent ⫾ standard deviation.

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

LAZAROVITCH ET AL.: SYSTEM-DEPENDENT BOUNDARY CONDITION

Fig. 4. Back pressures calculated using analytical [1,2] and numerical[4,5] solutions as a function of the final discharge rate Q∞ understeady state conditions for two soils (Magal clay loam and Arava
sandy loam), the source radius r0 of 0.01 m, the inlet pressure Pinof 10 m, and the nominal discharge rates Q0 of 2, 4, and 8 L h⫺1.

pressures were 3.55 and 0.27 m for the Magal and Arava
soils, respectively, with Q∞ of 1.67 ⫻ 10⫺6 m3 s⫺1. The
comparison between results obtained using the analytical and numerical models is again very good. That was
expected as both analytical and numerical solutions
solve the same governing Eq. [4], with similar boundary
conditions. The only difference, except for the solution
method, is the model of the soil hydraulic properties
considered. While the analytical solution considers theGardner’s model Eq. [3], the numerical model uses van
Genuchten–Mualem model Eq. [6] and [7]. Both models
use the same KS and φ, which explains the very good
agreement between them. Analysis of Eq. [2] reveals
that for given discharge Q0 and radius r0, the back pressure ␺s is much more sensitive to KS than to ␣G.
The effect of the cavity size is presented in Fig. 5,
where the linear Q∞–␺s relationships are plotted for twosource radii r0 of 0.0025 and 0.01 m. Calculations were
performed using the inlet pressure Pin of 10 m and theArava soil. Decreasing r0 from 0.01 to 0.0025 m resulted
in increase of the slope, suggesting faster increase of
the back pressure ␺s with the nominal discharge Q0,
resulting into a lower Q∞/Q0 ratio.

CONCLUSIONS
The main objective of this work was to implement
source characteristics to a transient numerical model.
This new system-dependent boundary condition allows
calculation of the source discharge while considering
source properties, inlet pressure, and effects of the soil
hydraulic properties. The updated numerical model was
verified against existing analytical solutions for simplified steady-state conditions and validated against collected experimental data. Good agreement was foundbetween transient back pressures measured in the labo-

49

Fig. 5. Back pressures calculated using analytical [1,2] and numerical[4,5] solutions as a function of the final discharge rate Q∞ understeady state conditions for two radii r0 of 0.0025 and 0.01 m, theinlet pressure Pin of 10 m, and the Arava soil.

ratory experiment and those calculated using the updated numerical model. The modified model allows theuse of any hydraulic model for soil and source properties, simulation of both short and long duration infiltration events, and consideration of various source geometrical shapes. The modified numerical model enablesconsidering complex processes in the root zone involving short-time irrigation and root water uptake.
Previous infiltration models that describe transientwater flow from subsurface sources neglect the phenomena where the source discharge can be limited by thedischarge at the infiltration boundary due to the effectsof the soil hydraulic properties that can limit the soil
infiltration capacity. This phenomenon can lead to anincrease in the water pressure head at the source, anda subsequent decrease in the source discharge rate. If
neglected, the dependence of the subsurface source discharge on the immediate soil hydraulic properties canlead to an application of smaller than designed irrigation
volumes when the irrigation amount is preset by timeand can lead to a non-uniformity in water application(Warrick and Shani, 1996).
The newly implemented system-dependent boundary
condition is solved for the case of homogenous and isotropic conditions near the source and outer boundaries
that are far enough so they do not affect the flow fields
near the source. These conditions are met with subsurface drip irrigation systems. These drippers have small
diameters and spatially variable conditions near the source
can therefore be neglected. A future work can include
non-homogenous conditions near the source (i.e., different soils or closer outer boundaries).
REFERENCES
Gardner, W.R. 1958. Some steady-state solutions of the unsaturated
moisture flow equation with application to evaporation from awater table. Soil Sci. 85:228–232.

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

50

SOIL SCI. SOC. AM. J., VOL. 69, JANUARY–FEBRUARY 2005

Hillel, D. 1971. Soil and water physical principles and processes. Academic Press, New York.
Klute, A. 1986. Water retention: Laboratory methods. p. 635–662. In
A Klute (ed.) Methods of soil analysis. Part 1. 2nd ed. ASA and
SSSA, Madison, WI.
Or, D., U. Shani, and A.W. Warrick. 2000. Subsurface tension permeametry. Water Resour. Res. 36:2043–2053.
Philip, J.R. 1968. Steady infiltration from buried point sources andspherical cavities. Water Resour. Res. 4:1039–1047.
Philip, J.R. 1992. What happens near a quasi-linear point source?
Water Resour. Res. 28:47–52.
Šimůnek, J., M. Sejna, and M. T. van Genuchten, 1999. The HYDRUS2D software package for simulating two-dimensional movementof water, heat, and multiple solutes in variably saturated media,
Ver. 2.0. Rep. IGWMC-TPS-53, Int. Ground Water Model. Cent.,
Colo. School of Mines, Golden.
Skaggs, T.H., T.J. Trout, J. Šimůnek, and P.J. Shouse. 2004. Comparison of HYDRUS-2D simulations of drip irrigation with experimental observations. J. Irrig. Drain. Eng. 130:304–310.
Shani, U., and D. Or. 1995. In situ method for estimating subsurface

unsaturated hydraulic conductivity. Water Resour. Res. 21:1863–
1870.
Shani, U., Xue, S., Gordin-Katz, R. and Warrick, A. W. 1996. Soillimiting flow from subsurface emitters. 1: Pressure measurements.
J. Irrig. Drain Eng. 122: 291–295.
van Genuchten, M.Th. 1980. A closed-form equation for predicting
the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am.
J. 44:892–898.
van Genuchten, M.Th., F.J. Leij, and S.R. Yates. 1991. The RETC
code for quantifying the hydraulic functions of unsaturated soils.
USEPA Rep. 600/2–91–065 (IAGDW1293–3-934). U.S. Salinity
Laboratory, USDA-ARS, Riverside, CA.
Warrick, A.W. 1993a. Comment on “What happens near a quasi-linear
point source?” by J. R. Philip. Water Resour. Res. 29:3299–3300.
Warrick, A.W. 1993b. Unsaturated-saturated flow near a quasi-linear
line source. Water Resour. Res. 29:3759–3762.
Warrick, A.W., and Shani, U. 1996. Soil limiting flow from subsurface
emitters. 2: Effect on uniformity. J. Irrig. Drain Eng. 122:296–300.
Warrick, A., and D. Or. 2001. Effect of gravity on steady infiltrationfrom spheroids. In J.R. Philip Tribute, Geophysical Monogr. Ser.
American Geophysical Union, Washington, DC.