Spherical Interpretation of Infiltration from Trickle Irrigation

Abstract

The hypothesis of this paper is that infiltration into drip irrigated soils can be described by simple spherical considerations as well as two dimensional (2D) numerical modeling. The major goal was to test a very simple model based on geometry of a sphere formulas, and compare it with elaborated numerical solutions and field experiments. Detailed analysis of soil–water infiltration under trickle regimes is shown to be pre-requisite in the search for the optimal design of system layout. Optimality and simplicity are sought by modeling a sphere for subsurface trickle/drip (SDI) and hemisphere (DI) pattern of moisture distributions during infiltration. Numerical simulations by MATLAB software were used to describe the distribution of soil water. The data produced by this simulation were successfully compared with analytical models and numerical results of Panoche clay loam. To simulate the four discharge rates (0.5, 1, 2, 3 $\frac{L}{h}$) under DI and SDI we used the input of Panoche soil properties, i.e., hydraulic conductivity function (K_θ) and soil water retention curve (ψ_θ). The resulting regression equation of numerical analysis (N) vs. spherical interpretation (S) was N = 0.97 × S − 19.1; r² = 0.98. This result exposes the novelty of the approach by showing that infiltration from a drip/trickle source can be described by simple spherical radial symmetry in addition to analytical or numerical simulations. An example of a design parameter for 3000 $\frac{{\mathrm{cm}}^3}{\mathrm{h}}$ suggested more emitters per meter laterals for SDI than for DI (100 vs. 77 unites, respectively) due to the shorter distance between SDI emitters that are required in order to maintain wetting continuity. At a discharge of around 500 $\frac{{\mathrm{cm}}^3}{\mathrm{h}}$ of three different soils’ SDI, positive pressure was detected near the orifice and it caused discharge reduction. This is a self-compensating property of SDI that regulates individual emitters according to the soil hydraulic properties. In conclusion SDI is associated with larger capital investment compared to DI, but it can be compensated by improving the water use efficiency due to increased productivity while reducing losses of water through evaporation, but this option should be investigated as part of specific research.

1. Introduction

Drip irrigation has been successfully used in the last few decades all over the world. It is considered a revolutionary method, and its contribution to agricultural production is one the most important in the last 200 years [1,2,3]. Today, the two dripping approaches are surface drip irrigation (DI) and subsurface drip irrigation (SDI).

There are many aspects used to make the comparison between DI and SDI. They include engineering, mathematical, physical, agronomic, economic, industrial, environmental, etc.

This paper is focused on the mathematical differences and engineering. Mathematically, on the DI surface, a shallow pond is created that is freely increased with wetting time. In mathematical terms, it is necessary to build a model with moving boundary conditions. On the SDI, water does not enter freely into the soil, but conversely a pressure buildup by the soil outside the orifice may reduce or even prevent water outlet away from the emitter.

From the engineering point of view, the depth of the SDI is an important factor. In shallow subsurface systems (<20 cm depth), water reaches the surface, and some is lost by evaporation, then replaced annually via thin walls laterals of DI, similar to surface drip pipes. Multiyear SDIs are installed below the tillage depth. SDI requires concentrated and consistent management of both water and nutrients to assure adequate crop performance, and they are also less forgiving than DI.

Economically, SDI irrigation systems may have a higher initial investment cost than DI.

Most users work with DI due to its simplicity and cost effectiveness. However, the use of SDI is increasing due to its agronomic advantages and, in water-limited regions, its safe application of wastewater to crops. In both DI and SDI, the particular problem of transient infiltration is still under intensive research and analysis. To simplify it, we present here a new model based on geometrical considerations to analyze water infiltration patterns under DI and SDI.

Here, we intend to re-examine the physical properties of SDI, to extend our understanding and compare it with DI. Pioneering research on SDI was published by Fok [4], who focused on spherical geometry to formulate design criteria, and Ben-Asher and Phene [5], who used the SUTRA numerical solution that was developed by Voss [6] to examine DI and SDI. The highlights of their study were that SDI maximizes transpiration at the expense of evaporation, and hence increases yield. Compared to DI, the disadvantage of SDI is its small wetted radius that requires more emitters. In 1996, Shani et al. [7] and Warrick and Shani [8] strengthened another finding of [5] by showing experimentally and theoretically that the larger the ratio between the designed emitters’ discharge and soil hydraulic conductivity, the stronger the reduction in the designed discharge due to pressure buildup in the surrounding soil. This is a unique feature of SDI compared to DI and other irrigation systems. The potential of SDI to facilitate emitter–soil interface interaction enables self-compensation in the system discharge, and obviates the industrial pressure compensation within the dripper itself. This is of high practical importance and deserves the re-evaluation carried out in this study.

Advanced computer technologies have been developed in the last two decades and used to deepen trickle irrigation research. Many current modelers are using HYDRUS 2 and 3D [9] to describe water and solute transport from a source of DI and SDI in the soil [9,10,11]. The results indicate that 2D for SDI has acceptable accuracy, and they have thus become a quasi-reference for other mathematical models, such as the moment analysis [12] of SDI and artificial neural networks [13] that compare it to DI, wherein HYDRUS is used as a reference.

The 2D infiltration studies based on this innovative mathematical method compared their results to HYDRUS, whose reliability has been proven in many studies, but they lack a comparison to the infiltration process in natural soil.

Therefore, in order to strengthen this study, we chose to combine the spherical radial model proposed here with a numerical model that also tested infiltration on Panoche clay loam.

Our hypothesis is that modeling infiltration into drip-irrigated soils can be simplified by spherical interpretation as well as two-dimensional (2D) numerical modeling. The major goal was to test the simple model based on the geometry of a sphere formula combined with a numerical solution and field experiments. It was assumed that the soil is a stable isotropic and homogenous porous media, DI is an ideal hemisphere, and SDI is an ideal sphere. This goal was obtained by following certain specific objectives.

Specific objectives:

(A)

Simulate DI and SDI using a MATLAB numerical model;

(B)

Validate the model by comparing it to the data of Panoche clay loam [14];

(C)

Formulate a simplified sphere (SDI) and hemisphere (DI) geometrical model;

(D)

Use experiments to validate and test the applicability of the sphere formulas using a 2D numerical model as the database;

(E)

Re-evaluate the uniqueness of SDI and compare it to DI.

2. Materials and Methods

This simulation was performed by solving the two-dimensional advection–diffusion equation, which is essentially an approximation of the Richards [15] equation. In this approach, two exponential function forms are used to represent soil hydraulic conductivity and volumetric water content. Numerical solutions of the advection–diffusion equation are obtained by means of the finite-difference scheme. Validation of the mathematical model was performed by means of its comparison with the finite–analytic–numerical solution. To mimic the 3D process we used a 2D mathematical model, and in order to validate its suitability for 3D reality, we conducted a field experiment with four actual water sources that were inserted into the soil. A compensating 2 Lhr⁻¹ “mother” dripper was evenly split into four spaghetti-like pipes; each one of them provided 500 $\frac{{\mathrm{cm}}^3}{\mathrm{h}}$ and was inserted into a clay loam soil (loess) at a depth of 50 cm. Water was applied for six hours to provide 3000 cm³, producing a real 3D sphere. Careful exposure of the entire spherical volume was done by excavation in all directions around it at the end of the experiment. It will be shown later that the results justified the use of radial symmetry for spherical geometry, and the wetting shape of the 2D model can properly mimic the 3D reality.

The simulations were performed for a square soil profile with $Z=X=85 \mathrm{cm}$. As in

Table 1. Emitters discharge, time table for each simulation and total water application of DI and SDI.

Emitters Discharge for DI and SDI $\left(\frac{{\mathbf{cm}}^3}{\mathbf{h}}\right)$	Application Durations (h)	Total Application (cm3)
500	6	3000
1000	3	3000
2000	1.5	3000
3000	1	3000

, the program in MATLAB was run four times under the same conditions for both DI and SDI sources.

Table 1 contains suitable details for the DI and SDI simulation scenarios. The total volume of applied water was 3000 cm³ for all scenarios. In the numerical experiments we used four discharge rates—500, 1000, 2000 and 3000 $\frac{{\mathrm{cm}}^3}{\mathrm{h}}$. Times for the scenarios (DI and SDI) ranged from 1 to 6 h.

3. Theory

3.1. Formulation of 2D Trickle Irrigation Model

Water flow in unsaturated soil is formulated with the functions given by Warrick [14]: (1) the soil is assumed to be a stable, isotropic homogeneous and non-swelling porous medium; (2) both the volumetric water content and the hydraulic conductivity of the soil are single-valued and continuous functions of the matrix flux potential [9]; (3) Darcy’s law, applicable in its modified form, given by Richards (1931) [15].

The 2D governing equation is

$$C\left(\psi \right)\frac{\partial \psi }{\partial t}=\frac{\partial }{\partial x}(k\left(\theta \right)\frac{\partial \psi \left(x,z,t\right)}{\partial x})+\frac{\partial }{\partial z}(k\left(\theta \right)\frac{\partial \psi \left(x,z,t\right)}{\partial z})-\frac{\partial k\left(\theta \right)}{\partial z}+R\left(x_0,z_0,t\right)$$

(1)

where t, x and z represent time, horizontal and the vеrtiсаl сооrdinаtеs, respectively. The vertical axis Z is directed downwards from the surface toward the lowest water content (zero). The horizontal axis $X$ shows the length of the problem. Symbols $\theta ,k,\psi ,R\left(x_0,y_0\right)$ represent the volumetric water content, hydraulic conductivity, рrеssurе head and water sink/source, respectively. In this paper we simulated only infiltration, such that R = 0.

The set of constitutive relations describing the dependence of the hydraulic conductivity and moisture content on the pressure head are employed in this study by the matrix flux potential of Gardner (1958) ([16] Equations (2)–(5))

$$k\left(\psi \right)=K_{s}exp(\alpha \psi ) \mathrm{for} \psi <0$$

(2)

$$k\left(\psi \right)=K_s \mathrm{for} \psi \ge 0$$

(3)

$$\theta \left(\psi \right)={\theta }_r+\left({\theta }_s-{\theta }_r\right)exp(\alpha \psi )$$

(4)

$$C\left(\psi \right)=\frac{\partial \theta }{\partial \psi }=\alpha \left({\theta }_s-{\theta }_r\right)exp(\alpha \psi )$$

(5)

Here, $k_s,\alpha ,{\theta }_s,{\theta }_r$ represent saturated hydraulic conductivity; α is a parameter representing soil pore–size distribution and θ_s and θ_r are the saturated and residual volumetric water contents.

Notice that the effective saturation or reduced water content $Sc=\frac{\left(Q-Q_r\right)}{\left(Q_s-Q_r\right)}$ varies between zero and one. It can be derived from Equation (4). On θ_r it can be said that it is somewhat arbitrarily defined, as the water content at which the corresponding hydraulic conductivity is appreciable is zero. However, very often, θ_r is used as an empirical constant when fitting measured hydraulic conductivity to Equation (4). Equation (5) is defined as the specific storage. It expresses the volume of water that is released from the soil profile per unit of decline in hydraulic head.

With the above relationships, the governing Equation (1) can be linearized and expressed as given below by Equation (6):

$$\frac{\partial \theta \left(x,z,t\right)}{\partial t}=D(\frac{{\partial }^2\theta \left(x,z,t\right)}{\partial x^2}+\frac{{\partial }^2\theta \left(x,z,t\right)}{\partial z^2})-V\frac{\partial \theta \left(x,z,t\right)}{\partial z}+R\left(x_0,z_0,t\right)$$

(6)

Equation (6) is also called the “advection–diffusion equation”. It is linear in $\theta$ with a given constant advection velocity $V=\frac{k_s}{\left({\theta }_s-{\theta }_r\right), }$, constant diffusivity $D=\frac{k_s}{\left(\alpha \left({\theta }_s-{\theta }_r\right)\right)}$, and a time-dependent inner source term $R\left(x_0,z_0,t\right)$.

3.2. Initial and Boundary Conditions

We consider the above governing equation applied to the soil profile defined as $\Omega \subset \left[0,X\left]x\right[0,Z\right]$, subjected to the predetermined initial and boundary conditions.

The general initial condition is:

$$\theta ={\theta }_i\left(x,z\right) \mathrm{if} 0\le x\le X; 0\le z\le Z \mathrm{and} t=0$$

(7)

Under the initial conditions, the function ${\theta }_i\left(x,z\right)$ is assumed to be constant in such a way that ${\theta }_r\le \theta \le {\theta }_s$, where ${\theta }_r$ and ${\theta }_s$ are the residual and saturated soil water contents.

We also have the no-flow conditions at the lateral boundaries of the soil profile:

$$\frac{\partial \theta \left(x,z,t\right)}{\partial x}=0, \mathrm{if}=0, \mathrm{and} x=X, 0\le z\le Z, t\ge 0$$

(8)

At the soil surface within the DI solution, the water is applied by the emitter with a constant flow rate during the operational time $t_p$, followed by the zero flow rate during the rest time $t_r$, i.e.,

$$-D\frac{\partial \vartheta \left(x,z,t\right)}{\partial z}+V\theta =Q_p-E_p, \mathrm{if} x=x_e, \mathrm{for} \mathrm{the} 0<t<t_p$$

(9)

and

$$-D\frac{\partial \vartheta \left(x,z,t\right)}{\partial z}+V\theta =-E_p \mathrm{if} 0\le x\le x_e \mathrm{and} x<x_e<X$$

(10)

where E_p is the rate of evaporation and V is the flow velocity of the advection in the $z$ direction.

At the bottom of the profile, the zero flux condition is predetermined:

$$-D\frac{\partial \theta \left(x,z,t\right)}{\partial z}+V\theta =0, \mathrm{if} 0\le x\le X,z=0,t\ge 0$$

(11)

The above mathematical model must be expanded to take into consideration the ponding phenomenon, i.e., the development of the radial area of the ponded water with the soil water content $\theta ={\theta }_s$ in the vicinity of the trickle irrigation source.

3.3. Numerical Simulation

The simulations were performed for a square soil domain of depth and width $\mathit{Z}=\mathit{X}=85 cm$. For the numerical solution of Equation (6), the simulated region is divided into a network of equally spaced grid-points, where each point is designated by the subscripts $\left(\mathit{i},\mathit{j}\right)$, referring to the horizontal and vertical directions, respectively. Additionally, the subscript $\mathit{n}$ is used to specify the $\mathit{n}-\mathbf{th}$ time level, and $\Delta \mathit{t}$ is the time-step. The solution for the time level $\left(\mathit{n}+1\right)$ is obtained from the solution for the time-level $\mathit{n}$ by means of the finite-difference scheme following the approach given by (Wang and Anderson, 1983 [17]):

$${\theta }_{ij}{}^{n+1}=(1-\frac{4D\Delta t}{\Delta x^2}){\theta }_{ij}{}^n+\frac{4D\Delta t}{\Delta x^2}(\frac{{\theta }^n{}_{i+1,j}+{\theta }^n{}_{i-1,j}+{\theta }^{}{}_{i,j+1}+{\theta }^{}{}_{i,j-1}}{4})+\frac{R_{ij}\left(x_0,z_0,t^n\right)}{\Delta x^2}-V\frac{\Delta t}{2\Delta x}\left({\theta }^n{}_{i,j+1}-{\theta }^n{}_{i,j-1}\right)$$

(12)

A computer program written in MATLAB was developed to simulate the water distribution patterns due to the trickle source placed at the soil surface (DI) or at the center of the soil profile $\Omega$(SDI).

To treat the ponding phenomenon within the present numerical scheme, it is assumed that if ${\theta }_{ij}>{\theta }_s$ then ${\theta }_{ij}={\theta }_s$. To check the relevance of this treatment, mass–balance control was performed through the simulation.

3.4. Capillary Length (λc)

We used λc to estimate the radius near the emitter (r₀) a short time after the irrigation began. As regards the capillary length, λc (L), it can be said that it substantiates the relative importance of capillarity over gravity forces during water movement in unsaturated soil. Low λc values (e.g., 0 < λc ≤ 1 cm) indicate the dominance of gravity over capillarity, and are typically found in coarse-textured soils. Alternately, high λc (i.e., 12.5 < λc < 100 cm) indicates capillarity dominance, as for example in Panoche clay loam, where pore radius affects flow.

Equation (13) approximates the final radius that maintains infiltration area into the soil at the same rate of dripper discharge.

This was thoroughly analyzed by Warrick [18], the bottom line of whose analysis is Equations (6)–(84), as given in Equation (13):

$$\lambda c=\frac{\theta s\ast D}{K_0}~\frac{208.2\ast 0.5}{4.2}=25 \mathrm{cm}$$

(13)

This equation is acceptable for a linearized case and the approximated capillary length λc (neglecting the very small dry hydraulic conductivity and dry θ). In Equation (13), D is the constant diffusivity (208.2 cm²/h), K₀ is the saturate hydraulic conductivity (4.2 cm/h) and θ_s is the saturated water content of the ponded region (0.5), assuming K_dry~0 and θ_dry~0.

With a known λc (25 cm) and q (for example: 500 cm³h⁻¹), the value of the short time r₀ was calculated from Equation (14) (i.e., Equations (6)–(82) in [18])

r₀ = q/(4π ∗ λc ∗ K₀)~500/(4 × 3.14 × 25 × 4.2) = 0.38

(14)

A summary of the data for short-term DI analysis is given in

Table 2. Short-times radius of saturated pond using surface drip irrigation.

Soil Parameters	K0$\left(\frac{\mathbf{cm}}{\mathbf{hr}}\right)$ = 4.2	λ (cm) = 25		Θs = 0.5
Dripper discharge q =	$\frac{{\mathrm{cm}}^3}{\mathrm{hr}}$	500	1000	2000	3000
Short time radius region r0 =	cm	0.39	0.75	1.5	2.27

.

4. Model Validations

4.1. Surface Drip Validation with Experimental Data

Whey-Fone et al. [19] provided an analytical solution for a strip source of constant flux, and solved the transient 2D water infiltration pattern.

The initial water content for the solution of the strip source infiltration into Panoche clay loam (Warrick [14]) is uniformly distributed with very dry conditions ${\theta }_i=2.3\times {10}^{-5}$ or $\psi =-250$ cm. The exponential forms of $K\left(\psi \right)$ and $\theta \left(\psi \right)$ are given by $K\left(\psi \right)=K_{0}exp(\alpha \psi )$ and $\theta \left(\psi \right)=\frac{K_0}{A_o}\left(exp(\alpha \psi \right))$.

Here, $K_o=0.0696 \frac{\mathrm{cm}}{\min }$ and $A_0=0.1388 \frac{\mathrm{cm}}{\min }$ are the set of the empirical constants obtained for the Panoche clay loam [14]. These parameters were also used in the work of Celia et al. [20]. Calculated contour plots are given for the hydraulic head when t = 36 min from the start of the infiltration (Figure 1A).

To compare the results of Figure 1A with the experimental data, we plotted a vertical cross-section, calculated by our model, and compared it successfully with the data of Celia et al. [20] in Figure 1B. The fairly good agreement supports the validity of the numerical solution.

4.2. Subsurface Drip with Experimental Data

Water infiltration from the SDI emitter is displayed in Figure 2A–C. This is an example of a comparison between an analytical solution [19] and our numerical solution that considers an inner source in the square domain, $\Omega =\left[0,X\left]x\right[0,Z\right]$:

$$\frac{\partial \theta \left(x,z,t\right)}{\partial t}=D(\frac{{\partial }^2\theta \left(x,z,t\right)}{\partial x^2}+\frac{{\partial }^2\theta \left(x,z,t\right)}{\partial z^2})-V\frac{\partial \theta \left(x,z,t\right)}{\partial z}$$

(15a)

$$\theta \left(x,z,0\right)=exp(-\frac{{(x-L/2)}^2}{D}-\frac{{(z-L/2)}^2}{D})$$

(15b)

Here, $\theta \left(x,z,t\right)$ is the transported scalar variable, $V$ is the speed of advection in the $z$ direction and $D$ is the constant diffusivity in $\Omega$. We solved the above equation with the initial condition at

Figure 2. (A) Distribution of soil water content after 36 min of infiltration. Applying 1800 cm³ (0.6 × 3000 $\frac{{\mathrm{cm}}^3}{\mathrm{hr}}$). Numerical values of θ are shown in (B). (B) displays water content after fifteen redistribution minutes (0.25 h). The solid horizontal line connects (A,B) to emphasize the 10 cm downward movement of the center of mass. (C) compares changes in the analytical and numerical solution of wetness $\frac{ \left(Q-Q_r\right)}{\left(Q_s-Q_r\right)}$ taken at different time steps at a single point located about 25 cm from the center of mass.

Figure 2. (A) Distribution of soil water content after 36 min of infiltration. Applying 1800 cm³ (0.6 × 3000 $\frac{{\mathrm{cm}}^3}{\mathrm{hr}}$). Numerical values of θ are shown in (B). (B) displays water content after fifteen redistribution minutes (0.25 h). The solid horizontal line connects (A,B) to emphasize the 10 cm downward movement of the center of mass. (C) compares changes in the analytical and numerical solution of wetness $\frac{ \left(Q-Q_r\right)}{\left(Q_s-Q_r\right)}$ taken at different time steps at a single point located about 25 cm from the center of mass.

4.3. Radial Interpretation of the DI and SDI Infiltration

Radial interpretation is suggested here as an approach to modeling water distribution from trickle emitters. It is a simplified infiltration model of DI and SDI, interpreting them as transient sphere and hemisphere.

In the analysis of 2D and 3D water flow in variably saturated soils, there are several approaches to mimicking the flow processes. For example, surface drip infiltration is characterized by a wetted bulb of the ponded water in the vicinity of its emitter. A short time after infiltration starts, the pond appears as a point that increases with water application time. When the dripper discharge “q” on a saturated dripping point is greater than the saturated soil hydraulic conductivity (K₀), a large part of the water does not penetrate the soil, but instead it forms a disk-like ponding surface. Then, all infiltrating water comes from this disk, and it is distributed around and below the saturated pond. The flow rate of the wetting pond into the soil is quoted from Warrick [18] in Equation (16) (See Equations (6)–(103) in Warrick [18])

$$q = \pi r_0^2 K_0 \left(1+\frac{4}{\pi \alpha r_0}\right)$$

(16)

where α is a constant for the analytical determination of hydraulic conductivity (L⁻¹). For a large disk, q is approximately π r₀² K₀. The second term corresponds to an edge effect and r₀ is the saturated radius of the ponded area around the dripper.

The radial modeling of the pond simplifies the numerical solution via Equations (17a)–(17c).

The maximum value of “r₀” is obtained by Equation (17a),

$$r_0 ~\sqrt{\frac{q}{K_0\pi }}r(\mathrm{t})=\sqrt{\frac{Qa}{4K_0\pi }}r=\sqrt{\frac{q}{K_0\pi }}$$

(17a)

The maximum surface area S of the pond around the emitter is

S = π r₀²

(17b)

The volume of saturated hemisphere below the ponded surface is θ_s~0.5; V_w = 0.5 × Vs, where Vw is the water volume of the hemisphere below the surface and within Vs

Vw = θs × π ∗ r₀³ × 2/3

(17c)

A quantitative summary of Equation (17) and its comparison with the numerical solution for Panoche clay loam is given in

Table 3. Large-time radius of ponded surface under DI, its area and volume of water as a function of drippers discharge. In all cases the volume of water applied was 3000 cm.

Effect of Dripper Discharge	Units					Applied cm3
Dripper discharge q =	$\frac{{\mathrm{cm}}^3}{\mathrm{hr}}$	500	1000	2000	3000	3000
Geometrically calculated ponded r0. Equation (17a)	cm	6.18	8.75	12.4	15.1	3000
Numerically calculated ponded r0 =	cm	5.5	9	12	15.5	3000
Geometrically calculated ponded area Equation (17b)	cm2	129	240	483	716	3000
Numerically calculated ponded area	cm2	95	459	452	706	3000
H2O amount in the ponded volume (geometrical Equation (17c))	cm3	247.5	701.5	1996.5	3605.5	3000
Numerical H2O amount in the ponded volume	cm3	174	763	1809.5	3534	3000

.

In Table 3, the radial calculations of “r₀” are based on dripper discharge rates, “q”, and the two soil parameters K₀ (in Equation (17a)) and θs (in Equation (17c)), thus, the total apparent soil volume was multiplied by θs. For the numerical simulation of infiltration processes, we used the set of equations specified in Section 3 (Theory). For all four discharges the average difference between numerical and radial calculations ((r₀ numerical − r₀ geometrical)/(r₀ numerical + r₀ geometrical)/2) was only about 1% of the calculated values. Concurrently, other hemispherical properties that are functions of “r₀” also differ from each other by low percentages.

4.4. Numerical Simulations of DI

The simulations of infiltration from the DI source were performed for the square soil profile with $Z=X=85 \mathrm{cm}$. The program in MATLAB was run four times (Figure 3A–D) under the same conditions without transpiration. We applied the same total volume of water (3000 cm³); therefore, the run time table was changed according to Table 1.

In Figure 3 the saturated pond is surrounded by a closed isoline. It can be viewed as a hemisphere with a radius “r₀” and water content $\theta$ = θs = 0.5, ($\psi =0$).

4.5. Subsurface Drip Irrigation (SDI)

The objectives of this paragraph were pursued to study the following unique properties of SDI:

(A) The spherical 3D reality, exposed experimentally. (B) The similarity between 2D numerical simulation and 3D reality. (C) The radial model of the sphere as compared with the 2D numerical simulation. (D) A re-examination of the pressure build-up caused by SDI.

Generally, the agreement between radial and numerical simulation is acceptable (

Table 4. Large-time radius of SDI sphere, its volume and its surface area as a function of drippers discharge. Total water application was 3000 cm³.

Dripper Discharge	$\mathbf{Units}: \frac{{\mathbf{cm}}^3}{\mathbf{hr}}$	500	1000	2000	3000
Geometrically calculated ponded (r0)	cm	5	6.2	7.8	8.9
Numerically calculated ponded (r0)	cm	3.5	6.0	7.5	8.5
Geometrically calculated transect area (s)	cm2	78.5	120.7	191.1	248.8
Numerically calculated transect area (s)	cm2	32.1	113.1	176.7	226.9
Applied H2O for ponded volume (v)	cm3	500	1000	2000	3000
Simulated ponded H2O volume(v)	cm3	179.5	998	1767	2572

), with one exception at a low discharge rate (500 $\frac{{\mathrm{cm}}^3}{\mathrm{hr}}$). In spite of this fact, the other simulations reinforce the approach presented here, and show that a model of drip infiltration can be approximated by radial spherical considerations.

In Table 4, the radius of the ponded zone is calculated by measuring the distance from the center of the sphere (i.e., 42.5 cm) to its saturated isoline. Taking into account all variables in Table 4 (r, s, and v), the average difference between the two modeling approaches of SDI was only 9% (without lowest discharge). The excavated sphere (Figure 5A,B) resulted from the SDI irrigation of clay loam displayed at the saturated zone an “ideal” sphere, but beyond it, in Figure 4 and Figure 5B, a set of unsaturated isolines displays reduced water content. It is assumed that the belt around the central saturated sphere of Figure 5B resulted from the lower water content of the isolines around it, and the gravity effect created isolines with an egg shape below it.

Figure 4A and Figure 5 describe the water distribution at the same discharge rates, and the wetting front of both creates an “ideal sphere” of the saturated zone (θ = 50%). The belt around it shows a gradual decrease in soil water content from 50% to 30%, to 15%, etc. (see Figure 4B–D). Its distribution and the wetting front depend on the dripper discharge. The larger the discharge, the greater the volume of the sphere. According to Figure 4A–D, high water contents advanced uniformly, close to the point source, but the unsaturated belt around it inclined downwards due to gravity, not only in Figure 4 but also in Figure 5B.

4.6. Soil-Limiting Flow from SDI Due to Point Pressure Build Up

One of the unique properties of SDI close to the dripper orifice is the pressure increase that opposes the outlet pressure from the dripper and reduces the flow of water movement from the dripper to the soil (5,7,8. This phenomenon causes the self-compensation of the SDI drippers’ discharge, and obviates the industrial compensation inside itself. At any point in the field, the water flow rate from the dripper can adjust itself to the ability of the soil to absorb the irrigated water. It should also be mentioned that in SDI, the soil particles around the dripper interfere with the water outlet, and can even clog it. In order to allow water infiltration into the soil, the pressure inside the emitter orifice must be greater than the opposing pressure of the soil water. This opposing pressure of the soil is shown in Figure 6.

The macroscopic observation of soil water potential near the orifice shows that it increases fast after the start of infiltration, before asymptotically approaching zero potential around a discharge of 200 $\frac{{\mathrm{cm}}^3}{\mathrm{hr}}$. However, microscopic observation shows that from 500 to 600 $\frac{{\mathrm{cm}}^3}{\mathrm{hr}}$, the soil water potential continues to increase above the zero potential line.

The internal figure within Figure 6 shows that for specific soil properties (Panoche clay loam with Ks = 4.2 $\frac{{\mathrm{cm}}^{}}{\mathrm{hr}}$), the soil water potential is zero when the emitter discharge is 560 $\frac{{\mathrm{cm}}^3}{\mathrm{hr}}$, and it is positive (around 2 cm) at a discharge of 600 $\frac{{\mathrm{cm}}^3}{\mathrm{hr}}$. Thus, in this situation, the build-up of pressure in the soil may reduce or even prevent water outlet away from the emitter. Shani et al. [7] measured this back-pressure and showed that it result in a substantial decrease in emitter discharge, especially if the soil has low conductivity, and for non-compensating emitters. Warrick and Shani [8] showed the significant deviation of the actual discharge from designed discharge. In their study, it was reduced to 0.905, 0.825, and 0.704 for designed discharges of 1, 2, and 4 L/h, respectively. Ben-Asher [21] provided complementary conclusions regarding the pressure build-up phenomena by emphasizing its practical implications. First, there is a potential for elimination of pressure-compensating emitters, because the variability of soil hydraulic properties creates a self-compensating environment. Second, each emitter may have individual flow rates according to its soil properties, and thus allow for longer laterals. Third, the decrease in emitter discharge due to the hydrostatic pressure opposing the the water pressure inside the emitter theoretically enables the use of emitters with a greater orifice, consequently enabling a greater flow rate and reduced orifice-clogging.

5. Discussion

5.1. General: Comparisons between Spherical/Radial and Numerical Modeling of DI and SDI Are Displayed in Figure 6

As anticipated, the slope of the regression in Figure 7 is close to a full unit (0.97), implying that the difference between the two methods of simulation is not large, and its correlation coefficient is high (r² = 0.98). This supports the hypothesis that infiltration into drip-irrigated soils can be described by simple spherical considerations when the dominant factor is the radius of the wetted pond. This is not only the case for the DI hemisphere, but also for the SDI sphere.

The radii in Table 3 for DI are consistently larger than the radii of SDI in Table 4. Both become increasingly larger with increased discharge.

Notice that the dashed line in Figure 6 may render ordinary least squares inadequate to account for the relationship between the spherical and numerical solutions. It was therefore important for this discussion to include a second linear equation (dashed line) that was based on the three parameters that belong to r, s and v, which were predicted by the lowest discharge in our study (500 $\frac{{\mathrm{cm}}^3}{\mathrm{hr}}$). The deviation from the main solid line was small (0.97), even when combining these points with the majority of the points in Figure 6 (solid line).

The importance of this result mainly relates to the design of drip irrigation systems. The radial model may provide a first approximation for planning a drip-irrigated field, assuming that an appropriate radius ensures the overlapping of the wet area between two adjacent emitters. We propose to maintain the overlapping of the wetting volume along the laterals network as a criterion for drip layout. This can be obtained when two adjacent emitters are placed at the same diameter interval. Accordingly, for a known saturated hydraulic conductivity (K₀), and with a wetting front approximated via Equation (17), a discharge of 3000 $\frac{{\mathrm{cm}}^3}{\mathrm{hr}}$ would result in a spacing interval of about 30 cm on DI (see Table 3), and only about 20 cm on SDI (see Table 4), provided that the layout criterion is based on the saturated wetting front.

5.2. Pressure Buildup and Capillary Length

Pressure build-up and capillary length are two of the unique properties of SDI that are calculated by the spherical interpretation of drip irrigation infiltration, and therefore deserve a special paragraph in the discussion. Capillary length λc was used here to estimate the radius of the pond created near the emitter (r₀) a short time after the irrigation began. When flow volume increases, the r₀ of the pond increases (Table 2), and the pressure head of the soil next to the emitter increases (Figure 6). It is most important to test the pressure build-up in soils other than Panoche clay loams. The results of this test are displayed in Figure 8.

Several calculation steps were carried out in order to obtain the required data set. This includes the calculation of K_(θ) as in Saxtion et al. [23], in which soil texture was provided by Noppol et al. [22]. The K₀ was 16, 7 and 4 cm hr⁻¹ for loam, sandy and clay loam, respectively.

It can be seen from Figure 7 that the dripper discharge of SDI is the most important factor affecting the pressure buildup within the dripper surrounding waterThe larger the discharge, the stronger the back-pressure.

The second factor is the type of soil, which affects saturated hydraulic conductivity. The greater this is, the stronger the counter pressure at a given flow rate.

The relation of r₀ to inflow rate, q (cm³/h), was observed via a typical plot (Warrick and Shani 8; not shown here) of r₀ $\mathrm{vs} q/\pi r^2$. When this relationship was plotted for all radii of the three soils on the same graph, the correlation of the equation was

$$r^0 = 13.4 \times {(q/\pi r^2)}^{-0.5};r^2 = 0.99$$

(18a)

and vice versa, $q/\pi r^2$ vs. 1/r₀ yields

$$q/\pi r^2=-4395+2035\times 1/r;r^2 = 0.98$$

(18b)

The correlations show that r₀ the radius of the water cavity around the dripper in SDI depends on the above ratio between the flow rate q and the projection area $\pi r^2$ inside the cavity.

6. Conclusions

More than 50 years have passed since the introduction of drip irrigation into agricultural farming and research. In spite of this, SDI, which is a version of drip irrigation, is still in its infancy, and therefore deserves additional research. The pressure build-up close to the orifice of the SDI dripper is a unique phenomenon that may have a positive impact on precision due to its ability to adjust itself to soil spatial variability via self-compensation. It may also have a negative impact due to its low irrigation uniformity. Another aspect that deserves more attention is the impact of capillary length variability, especially under SDI. The radial approximation of infiltration from a drip irrigation source agrees with the numerical simulation of DI and SDI. In addition, the 2D numerical simulation of the wetting bulb model properly mirrors the 3D reality.

Under the same conditions, there was an insignificant difference between the saturated water volumes in DI and SDI, simulated numerically and radially, excluding the very low (500 $\frac{{\mathrm{cm}}^3}{\mathrm{hr}}$) discharge rate.

The practical conclusions of this study are listed below:

(A)

When 3000 cm³ of water was applied to DI and SDI, the radius of the wetted surface of DI was larger than that of SDI because the infiltrating surface area of DI is not limited by soil hydraulic properties, whereas the transect of the SDI sphere is limited by the soil hydraulic properties and also by the pressure build-up, which reduces dripper discharge and is unique to SDI;

(B)

The simplicity of the radial model can be utilized as a first approximation for designing the network layout of DI and SDI. Moreover, to avoid evaporation from the soil surface of SDI, the calculated radius of the sphere may inform the choice of the installation depth.

During this analysis, several important aspects were not considered, in spite of their being especially relevant to the unique characteristics of drip irrigation. Subjects that were not studied here included redistribution, improving nutrient use efficiency by fertigation, and improving water use efficiency by reduced evapotranspiration. General water resources can be extended via the use of saline water, and this can also lead to a significant increase in agronomic productivity. Even though these missing aspects are important to decisions regarding whether drip irrigation can contribute to agronomic or economic productivity, this paper deals only with infiltration from a point source, and it can thus contribute to the design of a network layout.