Published in Soil Sci. Soc. Am. J. 69:46-50 (2005).
© 2005 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
Division S-1—Soil Physics
System-Dependent Boundary Condition for Water Flow from Subsurface Source
N. Lazarovitcha,
J. 
im
nekb and
U. Shania,*
a Dep. of Soil and Water Sciences, Faculty of Agricultural, Food and Environmental Sciences, The Hebrew University of Jerusalem, P.O. Box 12, Rehovot 76100, Israel
b Dep. of Environmental Sciences, University of California, Riverside, CA USA
* Corresponding author (shuri{at}agri.huji.ac.il).
 |
ABSTRACT
|
|---|
Infiltration rate of water from a subsurface cavity is affected by many 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) is larger 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 this work was to develop a boundary condition that describes this process while considering the source characteristics, and to implement this boundary condition into the transient numerical model Hydrus-2D. 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 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.
|
INTRODUCTION
|
|---|
WHEN A PREDETERMINED DISCHARGE of a subsurface source (e.g., a subsurface emitter) is larger than the soil infiltration capacity, the pressure head in the source 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 that the degree of the pressure difference reduction is larger for soils having lower hydraulic conductivity and that the positive pressure in the source vicinity increases rapidly at the beginning of the infiltration event and approaches a final value after only several minutes (Shani et al., 1996; Warrick and Shani, 1996). The source discharge is also affected by the source characteristics and the cavity size of the source outlet (Shani et al., 1996).
Analytical solutions for water flow from a subsurface cavity are available for calculating spatial distribution of the pressure head in the soil, 
(L). However, since these solutions were derived only for steady-state water flow 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 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 et al., 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 from subsurface water sources, and (v) to compare obtained experimental results with a numerical model.
|
THEORY
|
|---|
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]):
 | [1] |
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, and c (–) 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 physically required 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:
 | [2] |
where r0 (L) is the cavity (source) radius, Q
(L3 T–1) is the steady state discharge, KS (L T–1) is the saturated hydraulic conductivity 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 the solution presented in Eq. [2]:
 | [3] |
where is the pressure head (L). Calculation of Q is then done by simultaneous solution of Eq. [1] and [2] (Shani et al., 1996; Warrick and Shani, 1996).
In our analysis we used a numerical solution of the Richards' equation as implemented in the HYDRUS-2D code (imnek 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 the following mixed form of the Richards' equation:
 | [4] |
where 
is the volumetric water content (L3 L–3), K is the hydraulic 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 and the outer boundaries do not affect the flow field near the source, this boundary condition allows the calculation of the source discharge while considering properties of the source, and the inlet and back pressures:
 | [5] |
where 
iz is the Kronecker delta (–), 
is the boundary of the source, 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:
 | [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 (van Genuchten, 1980) was selected for the numerical simulations:
 | [7] |
 | [8] |
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 empirical shape parameters.
A flow domain was selected such that the outer boundaries do not affect the flow field near the cavity. The flow domain (50 x 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 (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 and that uses an algorithm based on the Delaunay's triangulation (imnek et al., 1999). All other boundaries, except for the free drainage boundary condition at the bottom of the flow domain and the source boundary condition (Eq. [5]) at the source-soil interface, had a no flow boundary condition.
View larger version (97K):
[in this window]
[in a new window]
|
Fig. 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 the other boundaries.
|
|
Experimental
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 was to 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 for comparison 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):
 | [9] |
of the two hydraulic models and substituting G = KS x 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 in Table 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 pressure transducer was inserted into the source to measure the pressure head in the source cavity. The initial water content was measured at several points and depths using time domain reflectometry (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, pressure transducers, and TDRs were attached to the data logger (CR10x, Campbell Scientific, Logan, UT). Measurements were taken every 0.1 s and averaged over a second. Each experiment (infiltration event) was replicated three times.
|
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 to be constant throughout the flow domain and equal to 0.07. For both sources s increased rapidly during the first 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.
View larger version (21K):
[in this window]
[in a new window]
|
Fig. 2. Measured (triangles and circles) and calculated (lines) soil pressure heads as a function of time for two subsurface sources (4 and 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.
|
|
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 s increase 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 x 10–6 to 1.73 x 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 the subsurface source discharge and the significant decrease of discharge due to development of the back pressure.
View larger version (15K):
[in this window]
[in a new window]
|
Fig. 3. Measured (triangles) and calculated (lines) discharge from a subsurface source as a function of time for an 8 L h–1 source and 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.
|
|
Philip's (1992) solution for steady-state water flow from a subsurface cavity implies a linear relationship between s and Q (see Eq. [2]). The back pressure s as a function of Q for two soils (Magal clay loam and Arava sandy loam) and three sources (2, 4, and 8 L h–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 by assuming the same criteria for approaching steady-state flow 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 back pressures were 3.55 and 0.27 m for the Magal and Arava soils, respectively, with Q of 1.67 x 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 the Gardner'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.
View larger version (18K):
[in this window]
[in a new window]
|
Fig. 4. Back pressures calculated using analytical [1,2] and numerical [4,5] solutions as a function of the final discharge rate Q under steady state conditions for two soils (Magal clay loam and Arava sandy loam), the source radius r0 of 0.01 m, the inlet pressure Pin of 10 m, and the nominal discharge rates Q0 of 2, 4, and 8 L h–1.
|
|
The effect of the cavity size is presented in Fig. 5
, where the linear Q–s relationships are plotted for two source radii r0 of 0.0025 and 0.01 m. Calculations were performed using the inlet pressure Pin of 10 m and the Arava 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.
View larger version (17K):
[in this window]
[in a new window]
|
Fig. 5. Back pressures calculated using analytical [1,2] and numerical [4,5] solutions as a function of the final discharge rate Q under steady state conditions for two radii r0 of 0.0025 and 0.01 m, the inlet pressure Pin of 10 m, and the Arava soil.
|
|
|
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 found between transient back pressures measured in the laboratory experiment and those calculated using the updated numerical model. The modified model allows the use 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 enables considering complex processes in the root zone involving short-time irrigation and root water uptake.
Previous infiltration models that describe transient water flow from subsurface sources neglect the phenomena where the source discharge can be limited by the discharge at the infiltration boundary due to the effects of the soil hydraulic properties that can limit the soil infiltration capacity. This phenomenon can lead to an increase in the water pressure head at the source, and a subsequent decrease in the source discharge rate. If neglected, the dependence of the subsurface source discharge on the immediate soil hydraulic properties can lead to an application of smaller than designed irrigation volumes when the irrigation amount is preset by time and 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).
Received for publication January 5, 2004.
|
REFERENCES
|
|---|
- Gardner, W.R. 1958. Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation from a water table. Soil Sci. 85:228–232.
- 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 and spherical 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.
- imnek, J., M. Sejna, and M. T. van Genuchten, 1999. The HYDRUS-2D software package for simulating two-dimensional movement of 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. imnek, 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. Soil limiting 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.[Abstract/Free Full Text]
- 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 infiltration from spheroids. In J.R. Philip Tribute, Geophysical Monogr. Ser. American Geophysical Union, Washington, DC.
This article has been cited by other articles:
|
|
|
|
 
J. Simunek, M. Th. van Genuchten, and M. Sejna
Development and Applications of the HYDRUS and STANMOD Software Packages and Related Codes
Vadose Zone J.,
May 27, 2008;
7(2):
587 - 600.
[Abstract]
[Full Text]
[PDF]
|
|
|
|
|
|
|
N. Lazarovitch, A. W. Warrick, A. Furman, and J. Simunek
Subsurface Water Distribution from Drip Irrigation Described by Moment Analyses
Vadose Zone J.,
January 24, 2007;
6(1):
116 - 123.
[Abstract]
[Full Text]
[PDF]
|
|
|
TOP
ABSTRACT
INTRODUCTION
