Published online 16 August 2005
Published in Vadose Zone J 4:838-847 (2005)
DOI: 10.2136/vzj2004.0160
© 2005 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
ORIGINAL RESEARCH
Numerical Analysis of Flow and Transport from Trickle Sources on a Spatially Heterogeneous Hillslope
David Russo*,
Jacob Zaidel and
Asher Laufer
AMEC Earth & Environmental Ltd., 160 Traders Blvd. East, Suite 110, Mississauga, ON, Canada, L4Z 3K7
* Corresponding author (vwrosd{at}volcani.agri.gov.il)
Received 9 November 2004.
 |
ABSTRACT
|
|---|
The purpose of the present study was to analyze movement and spreading of water and a passive solute (chloride) on a hillslope under surface drip irrigation, taking into account the texture of the soil, the slope of the terrain, the spatial heterogeneity of the soil hydraulic properties, and water extraction by plant roots. Results of the present investigation suggest that under surface drip irrigation, the movement and spread of water and solutes are affected mostly by the soil texture and less by the terrain slope. Increasing terrain slope is shown to increase the deflection of the trajectories of the centroid of the solute mass from the vertical axis, particularly in fine-textured soils associated with low saturated conductivity and considerable capillary forces, in which the interaction between adjacent drip line laterals may be appreciable. Transient flows originating from periodic water application and water uptake by plant roots are shown to enhance the effect of the terrain slope on water flow and solute movement, particularly in fine-textured soils. Implications regarding the problem of sensor placement with respect to drip irrigation management are briefly discussed.
Abbreviations: AW, asymptotic weighting • CDE, convection–dispersion equation • PPCG, polynomial preconditioned conjugate gradient • RSF, random space functions
|
INTRODUCTION
|
|---|
DRIP IRRIGATION as a method of irrigating large fields has become quite common practice in the last three decades. During this period, drip irrigation has been adapted to almost all types of crop production (strawberries [Fragaria X ananassa Duch.] to large-scale sugarcane [Saccharum officinarum L.]), land (level land planting of row and tree crops to steep hill-side planting of fruit crops) and irrigation water (nonsaline high quality water to brackish water with more than 2000 mg L–1 total salts), even those which previously could not be used for successful agriculture (Bresler, 1977). The possibility of applying water at a relatively low rate makes drip irrigation suitable for maintaining high-frequency water application, thus minimizing fluctuations in the soil water content in the effective root zone during the irrigation cycle (Rawlins, 1973). Furthermore, under drip irrigation, because the surface area across which water infiltration takes place is relatively small compared with the total soil surface of the irrigated field, the relatively high water velocity in the vicinity of the drippers during infiltration displaces the solutes beyond the main effective root zone.
Quantitative descriptions of water and solute movement and spreading, originating from a trickle source located at the soil surface, are essential for the design, management and scheduling of surface drip irrigation systems, with the aim of maximizing crop yields, and minimizing the potential of drainage water to pollute underlying groundwater supplies. To the best of our knowledge, there was no attempt to quantify water and solute movement and spreading under surface drip irrigation on a spatially-heterogeneous hillslope. The general objective of the present study, therefore, is to investigate movement and spreading of water and passive solute (chloride) originating from trickle sources, taking into account the texture of the soil (i.e., the magnitude of saturated conductivity and the capillary forces), the slope of the terrain, the spatial heterogeneity of the soil hydraulic properties, and water uptake by the plant roots.
To pursue this objective, we analyze the flow and the transport problem by means of physically based flow and transport models and a stochastic presentation of the soil properties that affect water flow and solute transport. The approach undertaken here is a "single realization" approach (e.g., Ababou, 1988; Russo, 1991; Polmann et al., 1991; Russo et al., 1994, 1998; Tseng and Jury, 1994), that may be understood as quantifying the macroscopic spreading of a single solute mass for particular site-specific applications.
The present study is a numerical experiment that may provide detailed information on the consequences of characteristics of the soil and the irrigation water for the movement and spreading of water and solute under realistic conditions. At the price of reduced generality, this approach circumvents most of the stringent assumptions of theoretical studies and facilitates analysis of simplified yet realistic situations, at a fraction of the cost of physical experiments.
|
THEORETICAL CONSIDERATIONS
|
|---|
The Physical Domain and the Simulated Scenario
Using a Cartesian coordinate system (x1,x2,x3), where x1 is directed vertically downward, a spatially-heterogeneous soil on a hillslope is considered. It is assumed that the slope of the soil surface is aligned parallel to the x2 axis such that the soil surface elevation, z, is given by the power law (Carson and Kirkby, 1972) as:
 | [1] |
where a and L are the height and length of the vertical and horizontal projections of the hillslope, respectively, and l and m are "shape" parameters.
The irrigated field on the hillslope consists of Nr rows of a crop, which, in turn, are aligned parallel to the x3 coordinate. Each row is irrigated with a single drip line lateral whose emitters are spaced very close to each other such that their ponding areas overlap very rapidly; consequently one can assume that the total ponding area of the emitters along the drip line lateral has the form of an infinite strip of width 
2(t), which, in turn, is oriented parallel to the x3 axis.
Governing Partial Differential Equations
We considered water flow and solute movement in a flow domain, 
, originating from a set of N surface drip line laterals that are aligned parallel to the x3 axis and are spaced at equal distance interval, d. We assume that the water flow is described locally by the Richards' equation, the physical parameters of which are visualized as realizations of stationary random space functions (RSFs). It is further assumed that the transport of the passive solute is described locally by the classical, one-region, convection–dispersion equation (CDE).
Viewing the drip line laterals as a set of infinite ponded strips of width 2(t), which, in turn, are oriented parallel to the x3 axis, the flow and the transport are independent of the x3 coordinate. If we consider water uptake by the plant roots, the "mixed" form of the Richards' equation that governs saturated-unsaturated flow in the vertical x1x2 plane is:
 | [2] |
where t is time; 
= (x,t) is the pressure head; 
= (x,t) is the volumetric water content; Kii, i = 1,2, are the principal components of the hydraulic conductivity tensor, taken as a symmetrical tensor of rank two with zero off-diagonal components, and Sw = Sw(x,t) is a sink term representing water uptake by plant roots.
Similarly, the CDE that governs transport of a passive solute in the vertical x1x2 plane is:
 | [3] |
where c = c(x,t) is the resident solute concentration, expressed as mass per unit volume of soil solution; ui (i = 1,2) is the Eulerian velocity vector, and Dij (i,j = 1,2) is the pore-scale dispersion tensor, given (Bear, 1972) as:
 | [4a] |
where 
L and T are the longitudinal and the transverse pore-scale dispersivities; 
ij is the Kronecker delta; |u| = 
1/2; and Dm is the effective molecular diffusion coefficient, given (Jury et al., 1983) as
 | [4b] |
where Dlw is the liquid diffusion coefficient in water, and s is the saturated volumetric water content.
Water is infiltrated into the soil from each of the drip line laterals through the ponding strips that develop along the drip line laterals, by an imposed time-dependent discharge, [F
= Q, t'j < t < t''j, F(t) = 0, elsewhere), where Q is the emitter discharge per unit length of the ponding strip, 
tI = t''j – t'j is the duration of an irrigation event, ti = t'j – t''j–1 is the time interval between successive events, j = 1 to NI, NI is the number of events, and tI + ti is the duration of an irrigation cycle.
Water may leave the flow system by uptake by plant roots, and by evaporation and drainage through the upper and lower horizontal faces of the flow domain, respectively. There is no water flow across the vertical faces of the flow domain, for which the normal derivatives of the pressure head vanish. In each irrigation, chloride with concentration, C0, invades the soil via the ponding strips along the drip line laterals. There is no solute transport across the upper horizontal face of the flow domain, nor across the vertical faces of the flow domain for which the normal derivatives of the chloride concentration vanish. In addition, zero-gradient-boundary is specified for the chloride concentration at the lower horizontal face of the flow domain.
Boundary and Initial Conditions
For the flow, the boundary and initial conditions for the N-drip line lateral system are:
 | [5a] |
 | [5b] |
 | [5c] |
 | [5d] |
 | [5e] |
where E is the evaporation rate, Gk, k = 1 to N, is given by Gk = k2, and k2 = k22–k21 is the time-dependent width of the ponded strip that develops along the kth drip line lateral, which, in turn, is oriented along the x3 axis.
For the transport, the boundary and initial conditions for the resident concentration, c, for the N-drip line lateral system are:
 | [6a] |
 | [6b] |
 | [6c] |
 | [6d] |
 | [6e] |
The Sink Term in Equation [2]
It is assumed here that locally, the rate of water uptake by plant roots is proportional to the unsaturated conductivity and to the difference between the total pressure head at the root–soil interface, 
r and the reduced water pressure head of the soil, + 
, where is the osmotic pressure head of the soil solution. According to this approach (e.g., Nimah and Hanks, 1973; Feddes et al., 1974, Bresler, 1987), the sink term Sw in Eq. [2] is
 | [7a] |
where Re(x,t), is the root effectiveness function, given by Re(x,t) = Rek(x) for x2 
bk, and Re(x,t) = 0 for x2 
bk, Rek is the root effectiveness function associated with the kth plant row, given by the time-invariant, Gaussian bivariate semilognormal spatial distribution (Coelho and Or, 1996) as
 | [7b] |
where rk = [(xk2 – x2)2]1/2, xk2 is the coordinate location of the kth plant row at the soil surface, x1 = z(x2), µh, µv, and 
2h and 2v are the mean values and the variances, respectively, of Rek(x) in the horizontal direction, and in the log-transformed vertical direction, respectively, ß is a scaling parameter, and bk = 
– 
is the so-called range of Rek (along the x2 axis), which shows the lateral extent of the influence of the roots of the kth plant row.
Characterization of the Flow and the Transport Parameters
It is assumed here that the local K() and () relationships are described by the van Genuchten (1980) parametric expressions. Ignoring hysteresis and local anisotropy, considering the pressure head, as the dependent variable, they read
 | [8a] |
 | [8b] |
where 
= ( – r)/(s – r) is the effective water saturation; s and r are the saturated and residual volumetric water contents, respectively; Ks is the saturated conductivity; 
VG and m are parameters which are related to the soil pore size distribution, and n = 1/(1 – m).
It is assumed further that each parameter in Eq. [8], denoted by p(x), is a second-order stationary, statistically anisotropic RSF, characterized completely by a constant mean, <p(x)> , independent of the spatial position, and a covariance function, Cpp(x,x'), that, in turn, depends on the separation vector, 
= x – x', and not on x and x' individually. A three-dimensional, exponential covariance is adopted here for p(x); that is,
 | [9] |
where ' = (x – x')/Ip is the scaled separation vector, ' = |'|; 2p and Ip = (Ip1,Ip2,Ip3) are the variance and correlation scales of p(x), respectively. Axisymmetric anisotropy is adopted for Cpp(), that is, Ipv = Ip1 and Iph = Ip2 = Ip3 are the characteristic length scales of p(x) in the vertical (longitudinal) direction, and in the horizontal plane, respectively. The values of 2p for the various soil parameters of Eq. [8] were adopted from Russo and Bouton (1992), and are summarized in Table 1; the correlation length scales of the various parameters of Eq. [8] were assumed as equal to those of logKs, adopted from Russo et al. (1997); that is, Iph = Ih = 0.8 m and Ipv = Iv = 0.2 m. Mean values of the soil parameters of Eq. [8] representing coarse- and fine-textured soils, (sandy and clay soils, respectively), were adopted from Mishra et al. (1989), and are summarized in Table 1. In addition, the values of L = 2 x 10–3 m and T = 1 x 10–4 m (Perkins and Johnston, 1963) were adopted for the two soils.
Generation of the Flow and Transport Parameters Field
In the present study we considered a vertical slice of the three-dimensional, actual heterogeneous flow domain, and thus we generated a hypothetical two-dimensional flow domain in the vertical x1x2 plane. The flow domain which extends over 3 and 5 m along the x1 and the x2 axes, respectively, was discretized into equalized numerical cells, 
, measuring 
1 = 0.010 m and 2 = 0.0167 m in the vertical and the horizontal directions, respectively. Using the covariance Eq. [9] and the values of the variance 2p (Table 1) and the aforementioned values of the correlation length scales, Ipv and Iph, the turning bands method (Tompson et al., 1989), coupled with a bilinear interpolation scheme, was used to generate realizations of logKs, logVG, n, s, and r for each of the 300 by 300 cells of the flow domain, taking into account the cross-correlation coefficients (Table 2) between the various soil properties.
For a given discretization of the flow domain, using Eq. [8], two-dimensional tables of relative conductivity, Kr
= K/Ks and relative water content, r
= 
/
were constructed as functions of 'VG = 
VG and n = 1.1 + ( – 1)n, where 'VG = VG||, = 1 to nVG, = 1 to nn, nVG = 50000, nn = 80, VG = 0.01 and n = 0.05. Employing a bilinear interpolation scheme, these tables were used to calculate values of Kr and = r(s – r) + r for each cell of the flow domain by means of the generated realizations of VG, n, s and r and the simulated pressure head, . Assuming local isotropy, the principal components of the hydraulic conductivity tensor in Eq. [2] are given by Kii = Ks(x)Kr(,x), i = 1,2,3. Hydraulic conductivity between cells, the so-called interblock conductivity, was estimated from the generated realization of Ks(x) and the calculated Kr(,x), using the asymptotic weighting (AW) scheme proposed by Zaidel and Russo (1992).
Implementation
A modified version of the numerical code of Russo et al. (1998) was used to solve Eq. [2] and [3], subject to Eq. [5] and [6], respectively, considering numerical cells above the elevation, z(x2), to be inactive for both flow and transport. Boundary conditions [5a] and [5b] and [6a] and [6b] were applied at the uppermost active cells in each vertical model column. In this code, the "mixed" form of Richards' equation governing two-dimensional water flow (Eq. [2]) was approximated by a fully implicit Euler scheme with a truncation error of O
. This scheme is convergent and unconditionally stable for the linear diffusion equation; it has been found appropriate for flow problems in highly heterogeneous porous media, in which both steep gradients and saturated regions can be developed. The resulting system of nonlinear algebraic equations with respect to the pressure head, , is solved iteratively, by applying the so-called modified Picard method (Celia et al., 1990). Picard (external) iterations are applied for both capillary and gravity terms and the resulting system of linear algebraic equations was solved by the polynomial preconditioned conjugate gradient (PPCG) method (e.g., Hageman and Young, 1981).
Inasmuch as the width of the ponded strips along the drip line laterals (Eq. [5a]) is not known a priori, the numerical code was combined with an iterative procedure. In this procedure, for each time-step during an irrigation event, [i.e., when, F(t) = Q], the width of the ponded strips which allows the fulfillment of the equality in Eq. [5a] was determined by appropriate adjustments. The iterative procedure continues until a prescribed error tolerance is met, or until the number of iterations has exceeded a predetermined limit.
For each time step, t, the CDE (Eq. [3]) governing solute transport in the two-dimensional flow domain was approximated by an operator-splitting approach (e.g., Konikow and Bredehoeft, 1978; Wheeler and Dawson, 1987). Under this approach, during the first stage of computation, only the advective part of the equation is solved over the time interval t. During the second stage, the solution obtained in the first stage becomes the initial condition for solving a pure dispersion equation over the same t. A numerical solution of the advection equation was obtained by means of a second-order accurate, explicit finite difference scheme proposed by Zaidel and Levi (1980). This scheme is capable of accurately describing both smooth and steep concentration profiles without introducing artificial numerical dispersion. Dispersive fluxes were approximated by a standard central difference scheme. For more details see the appendix in Russo et al. (1994).
Parameters of the time-invariant root effectiveness function Eq. [7] pertinent to corn (Zea mays L.) were adopted from Coelho and Or (1996); that is, b = 1, µh = 0, µv = –1.5, 2h = 0.09 m2 and 2v = 0.49. Water uptake by plant roots was implemented by a maximization iterative approach (Neuman et al., 1975). In this approach, the rate of transpiration per unit area of the soil surface at time t, given by Tr(t) = 
Sw(x,t)dx, is maximized subject to two requirements: (i) the actual rate of transpiration is not allowed to exceed the potential rate of transpiration, Tp and (ii) the total pressure head at the root–soil interface, r, is not allowed to fall below a critical value, c, equivalent to the so-called wilting point of the soil–plant system (e.g., c = –150 m). Meteorological data from the weather station at Bet Dagan were used to estimate potential evapotranspiration rates, ETp(t). The estimated cumulative potential evapotranspiration was 1005 mm yr–1. To simplify the problem, potential soil evaporation, Ep(t), was set to Ep(t) = ETp(t)/10 (ETp = Tp + Ep).
For each of the two soils, hillslope profiles, z(x2), with a linear slope (l = 1, m = 1), an increasing slope (l = 0, m = 1) and a decreasing slope (l = 2, m = 2) were considered in the simulations. The initial pressure head, i, for the two different soils was selected as the pressure head corresponding to <K(,x)>/<Ks(x)> = 0.0001, that is, i = –0.4 and –6.5 m, for the sandy and the clay soils, respectively. For the sandy and the clay soils, the flow and transport simulations proceeded for 60 and 240 d, respectively, taking into account four crop rows spaced 1-m apart with one drip line lateral per row, considering daily irrigations (with amounts equal to 1.1 times the daily amount of evapotranspitarion) and with chloride concentration, C0 = 10 molc m–3. In line with irrigation practice, the emitter discharges were selected as soil-texture dependent. The emitter discharges per unit length of the ponding strip were taken as Q = 312.5 and 1250 cm3 cm–1 h–1, for the clay and the sandy soils, respectively. For each soil, both the duration of the irrigation event, t''j – t'j, and the time interval between successive events, t'j – t''j–1, were adjusted to obtain the same daily amount of applied water.
|
RESULTS AND DISCUSSION
|
|---|
Analysis of the Flow
Water flow on a hillslope is controlled by the mean values and the spatial heterogeneity of the inherent properties of the soil (e.g., saturated hydraulic conductivity, and the parameters which relate the unsaturated conductivity, K, and water saturation, , to the pressure head, , (Eq. [8]), by the slope of the terrain, by the boundary and initial conditions imposed on the flow domain and by water uptake by plant roots. The movement and spread of a passive solute (chloride) on a hillslope are determined by pore-scale dispersion and by the spatial variability of the velocity vector, which, itself, is controlled by the spatial heterogeneity of inherent properties of the soil, by the slope of the terrain, by the boundary and initial conditions imposed on the flow and the transport domain, by water uptake by plant roots and by the spatial distribution of the water content, (x,t), which, itself, is also controlled by the aforementioned entities.
Flow regimes associated with the different soils and a given terrain slope in the presence of water uptake by plant roots and evaporation, are illustrated in Fig. 1, 2, and 3
. Figure 1 displays contour lines of the simulated pressure head distribution in the upper 2 m of the soil profile in the vicinity of the two inner drip line laterals, just after the cessation of the 14th irrigation event.
View larger version (69K):
[in this window]
[in a new window]
|
Fig. 1. Contours of the simulated pressure head, , in the vertical x1x2 plane at the end of an irrigation event. Results are depicted for both the sandy (top) and the clay (bottom) soils for Q = 312.5 and 1250 cm3 cm–1 h–1, respectively. Terrain slope is 20%
|
|
View larger version (40K):
[in this window]
[in a new window]
|
Fig. 2. Velocity fields in the vertical x1x2 plane just after an irrigation event for both the sandy (top) and the clay (bottom) soils. Note that the direction of the arrow is equal to the direction of the vector field at its base point; the longest vector is scaled to a length of 1, and all other vectors are scaled proportionately. Profiles of the components of the mean velocity vector, averaged over the two inner drip laterals, <u1> and <u2> , are also included in this figure. Terrain slope is 20%.
|
|
View larger version (46K):
[in this window]
[in a new window]
|
Fig. 3. Velocity fields in the vertical x1x2 plane just before an irrigation event for both the sandy (top) and the clay (bottom) soils. Note that the direction of the arrow is equal to the direction of the vector field at its base point; the longest vector is scaled to a length of 1, and all other vectors are scaled proportionately. Profiles of the components of the mean velocity vector, averaged over the two inner drip laterals, <u1> and <u2> , are also included in this figure. Terrain slope is 20%.
|
|
Note that the wriggles in Fig. 1 are a consequence of local variations in the hydraulic parameters and not of the grid discretization. Note also that because the emitter discharges associated with the different soils were adjusted to the soils' texture, the gap between the responses of the different soils is smaller than the gap that one would expect under the same emitter discharge.
In addition, water uptake by plant roots in response to the same transpiration demand imposed on the flow system, also narrows the gap between the responses of the different soils. Nevertheless, Fig. 1 clearly demonstrates the effect of soil texture on the pattern of the pressure head distribution. In the sandy, coarse-textured soil, with relatively high saturated conductivity and small capillary forces, gravity has a stronger effect. Consequently, the lateral extent of the ponded area in the vicinity of the drip line laterals is relatively small, the interaction between adjacent drip line laterals is small, and the flow is more concentrated along the vertical axis.
Velocity fields after an irrigation event and before a subsequent irrigation event, are depicted in Fig. 2 and 3, respectively, as patterns of arrows representing the velocity vectors in the vertical x1x2 plane. Profiles of the components of the mean velocity vector, averaged over the two inner drip laterals, <u1> and <u2>, are also included in these figures. The effect of soil texture on the velocity vector field after an irrigation event (Fig. 2) is clearly demonstrated. Detailed examination of the velocity vector field in this figure suggests that at the end of an irrigation event, because of the inclination of the terrain, at the ponded strips along the drip line laterals, the flow is directed uphill (i.e., the flow is more in a direction normal to the inclined soil surface). Below the wetted zones, however, the magnitude of the velocity vector decreases rapidly with increasing vertical distance, particularly when the soil is coarse textured. Note that the uphill flow (i.e., negative <u2> ) in the clay soil persists to a larger soil depth than in the sandy soil. At sufficiently large vertical distance from the drip line laterals, in both soils, the flow is essentially unidirectional vertical, with small local deviations due to the spatial heterogeneity in the soil properties.
Following a redistribution period between successive irrigations (Fig. 3), the velocity vector field is relaxed from the restraint imposed by the influx along the drip line laterals. The resultant flow pattern is rather complicated; it is clearly two dimensional and depends on the soil texture. Detailed examination of the velocity vector field close to the soil surface suggests that following a redistribution period, in the case of the sandy soil, the flow, particularly between the drip line laterals, is directed downhill (i.e., the flow is more in a direction parallel to the inclined soil surface). It may change its direction slightly with increasing soil depth, approaching a unidirectional vertical flow at sufficiently large distance from the soil surface. On the other hand, in the case of the clay soil, although in the vicinity of the soil surface the flow is directed downhill, it changes its direction with increasing soil depth considerably, particularly between the drip line laterals, leading to an uphill mean flow. At sufficiently large vertical distance from the soil surface, the flow approaches unidirectional vertical flow.
Additional simulations using different shapes for the slope of the terrain (not shown here), suggest that under transient flow conditions, the pattern of the pressure head distribution within the wetted zones in the vicinity of the drip line laterals is essentially robust to the shape of the slope of the terrain as long as the mean slope is kept the same.
Analysis of the Transport
Solute spreading patterns in the vicinity of the two inner drip line laterals, associated with the different soils and a given terrain slope are illustrated in Fig. 4
. This figure displays contour lines of the simulated solute resident concentrations, c, in the soil profile for the sandy and the clay soils at elapsed times, t = 10 and 56 d, respectively, which, in turn, correspond to the times at which the center of mass of the solute body is displaced to a vertical distance of 0.45 m.
View larger version (63K):
[in this window]
[in a new window]
|
Fig. 4. Contours of the chloride concentration, c, in the vertical x1x2 plane. Results are depicted for the sandy soil (top) and the clay soil (bottom), for elapsed times of t = 10 d and t = 56 d, respectively. Terrain slope is 20%.
|
|
Figure 4 clearly demonstrates the effect of soil texture on the pattern of the solute concentrations. In a coarse-textured sandy soil with relatively high saturated conductivity and small capillary forces, both solute movement and solute spreading in the longitudinal (vertical) direction is larger, while solute spreading in the transverse (horizontal) direction is smaller than in a fine-textured clay soil with low saturated conductivity and significant capillary forces. Consequently, for a given cumulative amount of applied water, solute leaching is more significant in the sandy soil than in the clay soil. In turn, the enhanced solute leaching associated with the sandy soil compensates in part for the effect of both water uptake by plant roots and evaporation to concentrate the soil solution.
This is further demonstrated in Fig. 5
in which chloride concentration profiles below the two inner drip line laterals and midway between the two drip line laterals, are depicted for the different soils. Note that for each soil, the concentration peak associated with the downhill drip line lateral is slightly smaller than its counterpart. Below the drip line laterals the concentration profiles exhibit symmetrical Fickian distribution, with concentration peak that increases as the soil texture becomes finer. On the other hand, midway between the two adjacent drip line laterals, the concentration profiles become more irregular and less symmetric, with an apparent bimodality, particularly when the soil texture is finer. The latter stems from the evaporation from the soil surface and the subsequent salt accumulation, and, concurrently, the buildup of a second concentration peak close to the soil surface.
View larger version (23K):
[in this window]
[in a new window]
|
Fig. 5. Profiles of the chloride concentration below the two inner drip line laterals and midway between the two drip line laterals, corresponding to the chloride concentrations in Fig. 4.
|
|
To quantify the effect of the terrain slope on the transport, integrated measures of the solute transport were calculated. After Aris (1956), the time-dependent, spatial moments of the distribution of the point values of the resident concentration, c(x,t), for the two inner drip line laterals, are given by
 | [10a] |
 | [10b] |
 | [10c] |
where i = 1,2, x'22 = x22 – d/2, x'32 = x32 + d/2, x22 and x32 are the positions of the second and the third (inner) drip line laterals along the x2 axis, respectively, d is the lateral spacing, Mc(t) is the dissolved mass of the solute, Rci (i = 1,2) is the coordinate vector of the centroid of the solute body, and S'ij
(i,j = 1,2) are second spatial moments, proportional to the moments of inertia of the solute body. In other words, Eq. [10a], [10b], and [10c], respectively, provide measures of the mass, location, and spread of the solute body.
Results of the analysis of the time-dependent solute concentration spatial distribution using Eq. [10], suggest that the solute spread about its centroid, expressed in terms of the principal components of the spatial covariance tensor, Sij = S'ij – S'ij
, depends on the soil texture, and, essentially, is independent of the terrain slope. On the other hand, the trajectories of the centroid of the solute bodies depicted in Fig. 6
suggest that the coordinate location of the centroid of the solute body depends on the terrain slope, particularly in the fine-textured soil associated with considerable capillary forces in which the interaction between adjacent drip line laterals may be appreciable. In this case, increasing terrain slope is shown to increase the uphill deflection of the trajectories of the centroid of the solute body and to increase the vertical distance at which they start to deflect downhill. On the other hand, when the soil is relatively coarse-textured, the coordinate location of the centroid the solute body only slightly depends on the terrain slope.
View larger version (32K):
[in this window]
[in a new window]
|
Fig. 6. Trajectories of the centroid of the solute body associated with the two inner drip line laterals under transient flow conditions and in the presence of water uptake by plant roots. Negative values of the horizontal coordinate indicate uphill deflection of the centroid.
|
|
The behavior of the trajectories of the centroid of the solute body depicted in Fig. 6 is in agreement with Fig. 2 and 3, which show that for a given terrain slope, the persistence of the flow that is directed uphill, which, in turn, is originating from the inclined drip line laterals at the soil surface, and the persistence of flow that is directed downhill, at soil depths below the main root zone, are larger in the clay soil than in the sandy soil. Additional simulations using different shapes for the slope of the terrain (not shown here), suggest that for the horizontal length scale considered here (few meters), under transient flow conditions, the pattern of the solute concentration distribution in the vicinity of the drip line laterals is essentially robust to the shape of the slope of the terrain as long as the mean slope is kept the same.
Comparison with Transport under Steady-State Flow
Trajectories of the centroid of the solute bodies associated with steady-state flows in the absence of water uptake by plant roots, are depicted in Fig. 7
for the two soils and different terrain slopes. Figure 7 suggests that the coordinate location of the centroid of the solute body depends on both terrain slope and soil texture. In the sandy soil, increasing terrain slope is shown to increase the uphill deflection of the centroid trajectories at relatively shallow soil depth, and to postpone the tendency of the centroid to coincide with the vertical axis at deeper soil depth. On the other hand, in the clay soil, for a given terrain slope, close to the soil surface, the uphill deflection of the centroid trajectories is smaller than that in the sandy soil. Furthermore, in the clay soil, a secondary uphill deflection of the centroid trajectories occurs at deeper soil depth.
View larger version (27K):
[in this window]
[in a new window]
|
Fig. 7. Trajectories of the centroid of the solute body associated with the two inner drip line laterals under steady-state flow conditions and in the absence of water uptake by plant roots. Negative values of the horizontal coordinate indicate uphill deflection of the centroid.
|
|
The displacement of the solute body depicted in Fig. 7 is in agreement with the steady-state velocity vector fields associated with the two soils (not shown here). In both soils, because of the inclination of the terrain, at the ponded strips along the drip line laterals the mean flow is directed uphill (i.e., clockwise from the negative vertical axis); however, it rotates counter clockwise with increasing depth. On the other hand, in the clay soil in which the interaction between adjacent drip line laterals may be appreciable, at a soil depth which decreases with increasing terrain slope, the mean flow is affected by a flow that originates from the second (downhill) drip line lateral which is directed uphill. At sufficiently large vertical distance from the ponded strips, in both soils, the mean flow coincides with the vertical coordinate.
Comparison of Fig. 6 and 7 suggests that as compared with steady-state flow, transient flow originating from periodic water application and water uptake by plant roots may enhance the effect of the terrain slope on water flow and solute movement, particularly in fine-textured soils associated with low saturated conductivity and significant capillary forces.
Applications
Monitoring of the temporal changes in soil water and solute concentration in the root zone and their proper interpretation, is a prerequisite to an accurate drip irrigation scheduling and management. With respect to the problem of sensor placement, one possible location for the sensor probe is the position of the centroid of the volume of water extracted by plant roots per unit length, given by the coordinate vector, Rui(t) (i = 1,2); that is,
 | [11a] |
where w(x,t) is the daily-average rate of water uptake by plant roots, and Mu(t) is the volume of water extracted by the plant roots per unit length, given by
 | [11b] |
The rate of water uptake by plant roots, Sw(x,t), and, concurrently, w(x,t) and Ru(t), are affected by the spatial and temporal distribution of both the pressure head and solute concentration, which, in turn, are affected by the mean soil hydraulic properties (i.e., soil type) and their spatial distribution, the terrain slope, the drip line lateral discharge, the root effectiveness function, and the time-dependent potential evapotranspiration.
Components of the coordinate vector of the centroid of the volume of water extracted by the plant roots, Ru, associated with the two inner drip line laterals, are depicted in Fig. 8
as functions of time, for the different soils and different slopes of the terrain. In the case of the sandy soil, for a given terrain slope, as time increases, the centroid of the water extracted by plant roots moves to a shallower soil depth and is deflected uphill. Note that in the sandy soil, the vertical position of the centroid of the water extracted by plant roots is independent of the terrain slope (Fig. 8a), while its uphill deflection slightly increases with increasing terrain slope (Fig. 8c).
View larger version (21K):
[in this window]
[in a new window]
|
Fig. 8. Longitudinal (a,b) and transverse (c,d) components of the coordinate vector of the centroid of the volume of water extracted by the plant roots, Ru, associated with the two inner drip line laterals, as functions of time, for the different soils and different slopes of the terrain. Negative values of the horizontal coordinate indicate uphill deflection of the centroid.
|
|
In the case of clay soil, for a given terrain slope, initially, the centroid of the water extracted by plant roots moves to a deeper soil depth and is deflected downhill. As time further increases, the vertical position of the centroid of the water extracted by plant roots approaches a constant depth, while the extent of its downhill deflection decreases with increasing time. Note that in the case of the clay soil, both the vertical position of the centroid of the water extracted by plant roots and the extent of its downhill deflection slightly increase with increasing terrain slope (Fig. 8b, 8d).
Considering a set of drip line laterals that are aligned perpendicular to the terrain slope, the results in Fig. 8 suggest that in the case of a sandy soil, the sensor probe should be placed at a soil depth which decreases with increasing time, uphill from the vertical mid-plane between two adjacent drip line laterals, at a transverse distance which slightly increases with increasing terrain slope. On the other hand, in the case of the clay soil, the sensor can be placed at a time-invariant fixed soil depth, downhill from the vertical mid-plane between two adjacent drip line laterals, at a transverse distance which slightly increases with increasing terrain slope.
|
SUMMARY AND CONCLUSIONS
|
|---|
The main purpose of the present study was to investigate movement and spreading of water and a passive solute (chloride) on a hillslope under surface drip irrigation, taking into account the texture of the soil, the slope of the terrain, the spatial heterogeneity of the soil hydraulic properties, and water extraction by the plant roots. This is a rather complicated flow system in which the flow, and, concurrently, the transport, are affected by the terrain slope, by capillary and gravity forces, by the mean values and the spatial heterogeneity of the soil hydraulic properties, and by the extent of the (partial) wetting of the soil surface (i.e., the magnitude of the interaction between adjacent drip line laterals) on the one hand, and by the periodic infiltration and water extraction by plant roots, on the other hand.
Results of the present investigation suggest that under surface drip irrigation, the movement and spread of water and solutes are affected mostly by the soil texture and less by the terrain slope. Increasing terrain slope is shown to increase the deflection of the trajectories of the centroid of the solute mass from the vertical axis, particularly in fine-textured soils associated with relatively low saturated conductivity and considerable capillary forces, in which the interaction between adjacent drip line laterals may be appreciable. Transient flows originating from periodic water application and water uptake by plant roots are shown to enhance the effect of the terrain slope on water flow and solute movement, particularly in fine-textured soils.
We would like to stress that the numerical experiments conducted in the present study provide detailed information on the consequences of soil and irrigation water characteristics for the water flow and solute transport under quite realistic conditions; information that, in general, cannot be obtained in practice from field investigations. We would like to emphasize, however, that the conclusions drawn from the present study should be considered with caution, inasmuch as the numerical results presented here are based on analyses of single realizations of the medium properties. Nevertheless, we believe that the simulated results are sufficiently accurate to indicate appropriate trends.
|
ACKNOWLEDGMENTS
|
|---|
This is contribution 625/04 from the Institute of Soils, Water and Environmental Sciences, ARO, the Volcani Center, Bet Dagan, Israel.
|
REFERENCES
|
|---|
- Ababou, R. 1988. Three-dimensional flow in random porous media. Ph. D. thesis. Massachusetts Inst. of Technol., Cambridge.
- Aris, R. 1956. On the dispersion of a solute in a fluid through a tube. Proc. R. Soc. London A 235:67–78.
- Bear, J. 1972. Dynamics of fluids in porous media. Elsevier, New York.
- Bresler, E. 1977. Trickle-drip irrigation: Principles and application to soil water management. Adv. Agron. 29:343–393.
- Bresler, E. 1987. Application of conceptual model to irrigation requirements and salt tolerance of crops. Soil Sci. Soc. Am. J. 51:788–793.[Abstract/Free Full Text]
- Carson, M.A., and M.J. Kirkby. 1972. Hillslope form and process. Cambridge University Press, Cambridge, UK.
- Celia, M., E. Bouloutas, and R. Zarba. 1990. A general mass-conservative numerical solution of the unsaturated flow equation. Water Resour. Res. 26:1483–1496.[CrossRef]
- Coelho, F.E., and D. Or. 1996. A parametric model for two-dimensional water uptake intensity by corn roots under drip irrigation. Soil Sci. Soc. Am. J. 60:1039–1049.[Abstract/Free Full Text]
- Feddes, R.E., E. Bresler, and S.P. Neuman. 1974. Field test of a modified numerical model for water uptake by root systems. Water Resour. Res. 10:1199–1206.[CrossRef]
- Hageman, L.A., and D.M. Young. 1981. Applied iterative methods. Academic Press, New York.
- Jury, W.A., W.F. Spencer, and W.J. Farmer. 1983. Behavior assessment model for trace organic in soil. I. Model description. J. Environ. Qual. 12:558–564.[Abstract/Free Full Text]
- Konikow, L.F., and J.D. Bredehoeft. 1978. Computer models of two-dimensional solute transport in groundwater, in Techniques of Water-Resources Investigations of the U.S Geological Survey. Chap. C2. U.S. Gov. Print. Office, Washington, DC.
- Mishra, S., J.C. Parker, and N. Singhal. 1989. Estimation of soil hydraulic properties and their uncertainties from particle size distribution data. J. Hydrol. (Amsterdam) 108:1–18.
- Neuman, S.P., R.A. Feddes, and E. Bresler. 1975. Finite element analysis of two-dimensional flow in soils considering water uptake by roots. 1. Theory. Soil Sci. Soc. Am. Proc. 39:224–230.
- Nimah, M.N., and R.J. Hanks. 1973. Model for estimating soil water, plant and atmospheric relationships. 1. Description and sensitivity. Soil Sci. Soc. Am. Proc. 37:522–527.
- Perkins, T.K., and O.C. Johnston. 1963. A review of diffusion and dispersion in porous media. Soc. Pet. Eng. J. 3:70–84.
- Polmann, D.J., D. Mclaughlin, S. Luis, L.W. Gelhar, and R. Ababou. 1991. Stochastic modeling of large-scale flow in heterogeneous unsaturated soils. Water Resour. Res. 27:1447–1458.[CrossRef]
- Rawlins, S.L. 1973. Principles of managing high frequency irrigation. Soil Sci. Soc. Am. Proc. 37:626–629.
- Russo, D. 1991. Stochastic analysis of vadose-zone solute transport in a vertical cross section of heterogeneous soil during nonsteady water flow. Water Resour. Res. 27:267–283.[CrossRef]
- Russo, D., and M. Bouton. 1992. Statistical analysis of spatial variability in unsaturated flow parameters. Water Resour. Res. 28:1911–1925.[CrossRef]
- Russo, D., J. Zaidel, and A. Laufer. 1994. Stochastic analysis of solute transport in partially saturated heterogeneous soil: I. Numerical experiments. Water Resour. Res. 30:769–779.[CrossRef]
- Russo, D., I. Russo, and A. Laufer. 1997. On the spatial variability of parameters of the unsaturated hydraulic conductivity. Water Resour. Res. 33:946–956.
- Russo, D., J. Zaidel, and A. Laufer. 1998. Numerical analysis of flow and transport in a three-dimensional partially saturated heterogeneous soil. Water Resour. Res. 34: 1451–1468.[CrossRef]
- Sposito, G., R.E. White, P.R. Darrah, and W.A. Jury. 1986. A transfer function model of solute transport through soil. 3. The convection-dispersion equation. Water Resour. Res. 22:255–262.
- Tompson, A.F.B., R. Ababou, and L.W. Gelhar. 1989. Implementation of the three-dimensional turning bands random field generator. Water Resour. Res. 25:2227–2243.[CrossRef]
- Tseng, P.H., and W.A. Jury. 1994. Comparison of transfer function and deterministic modeling of area-average solute transport in heterogeneous field. Water Resour. Res. 30:2051–2064.[CrossRef]
- 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]
- Wheeler, M.F., and C.N. Dawson. 1987. An operator splitting method for advection-diffusion-reaction problems. Tech. Rep. 87–9. Dep. of Math. Sci., Rice Univ., Houston, TX.
- Zaidel, J.M., and B.I. Levi. 1980. Numerical schemes for the solution of two-phase, three-component flow and transport in porous medium equations. Numer. Methods Continuum Mechanics 11(5):75–89 [in Russian].
- Zaidel, J., and D. Russo. 1992. Estimation of finite difference interblock conductivities for infiltration into initially dry soils. Water Resour. Res. 28:2285–2295.[CrossRef]
This article has been cited by other articles:
|
|
|
|
 
D. Russo, J. Zaidel, and A. Laufer
Numerical Analysis of Solute Transport from Trickle Sources in a Combined Desert Soil-Imported Soil Flow System
Vadose Zone J.,
January 23, 2008;
7(1):
53 - 66.
[Abstract]
[Full Text]
[PDF]
|
|
|
TOP
ABSTRACT
INTRODUCTION
of the hydraulic parameters of the van Genuchten model Ks,