Experimental procedure
To observe the soil moisture flow in DI system under different conditions, experiments were carried out in the Water Resources Engineering Laboratory of Civil Engineering Department at NIT Patna, India, during January 2022. Experimental setup presented in Fig 1 consists of a rectangular box with length, width and height 1, 0.6 and 0.75 m, respectively. The front side of the box consists of a transparent glass. To prevent the preferential flow, the glass surface was roughened by applying glue and spraying a coarse sand over it by preserving the transparency of glass (
Kandelous and Simunek, 2010). Water was supplied from the storage tank and pumped with the help of a motor through the low-density polyethylene pipes. Constant pressure was maintained in the laterals throughout the irrigation supply. Valves were used to regulate the discharge to the experimental setup and additional water was diverted back to the tank. Emitter was placed at the center of the length and near to the transparent face of box.
The soil was collected from the field of NIT Patna and dried for 48 hrs to attain uniform moisture content. It was sieved through 2 mm size sieve to avoid lumps, stones and gravels. Soil texture was sandy loam and characteristics of soil are presented in Table 1. Fig 2 presents soil water retention curve, which was obtained using pressure plate apparatus in the above lab for the pressures between 0.1 to 15 bar. Parameters of van Genuchten characteristic curve (α, n) were computed using MATLAB and values are presented in Table 1. Soil was placed in the experimental setup with known bulk density. The initial moisture content was noted.
Three sets of experiment were conducted with different discharge rates 0.5, 0.85 and 2 lph for 20, 12, 5 hrs irrigation period, respectively to supply total 10 litres of irrigation in all sets. Wetting pattern of the soil was marked at regular intervals of time for 48 hrs.
Numerical procedure
The essential aspect in the design of DI is soil wetting pattern, which depends on the discharge rate, frequency, soil characteristics, plant characteristics, spacing of laterals and emitters
(Prajapat et al., 2020).
Governing equations
The soil moisture movement in the unsaturated zone is assumed as 1-D transient flow in vertical direction assuming incompressible fluid and isotropic incompressible porous medium,
i.e., the mixed form of 1-D RE,
Where,
θ = Moisture content (cm
3/cm
3).
Ψ = Pressure/suction head (cm).
K(Ψ) = Unsaturated hydraulic conductivity (cm/sec).
z = Vertical distance.
t = Time (sec).
Soil moisture characteristic equations of
Haverkamp et al. (1977) and
Van Genuchten, (1980) equations have been used.
Haverkamp equation:
Van Genuchten equation:
Where,
K
s = Saturated hydraulic conductivity (cm/sec).
Ψ = Suction head (cm).
θ(Ψ) = Moisture content depends on the pressure head (cm
3/cm
3).
θ
s = Saturated moisture content (cm
3/cm
3).
θ
r = Residual moisture content (cm
3/cm
3).
A, B, a, b = Haverkamp parameters.
α, m, n,
l = Van Genuchten parameters.
S
e = Effective water content.
Numerical solution
Governing equation is a non-linear partial differential equation and solved using fully implicit finite-difference scheme with appropriate initial and boundary conditions (BC). The subsurface flow domain is divided into a soil column of rectangular blocks with node size, Dz in vertical direction and specifying moisture content, q and pressure head, y at the centre of the block (node). In contrast, the flow velocity is defined at the interblock faces, as shown in Fig 3.
The finite-difference form of Eq. (1) is
Where,
superscripts ‘n’ = Known time level.
‘n+1’ = Unknown time level.
subscript ‘j’ = Block number in the z-direction.
V = The time-averaged value of the velocity.
Δt and Δz = Time spacing and nodal spacing in z-direction, respectively.
The time-averaged velocity is given as:
Where,
w = Time weighting factor.
w = 1.0 for a fully implicit scheme.
Velocity at any interblock was computed using the pressure heads at the adjoining nodes
Where,
K
II= Unsaturated hydraulic conductivity at interblock face II.
K
IV= Unsaturated hydraulic conductivity at interblock face IV.
The unsaturated hydraulic conductivity at the interblock face was computed using weighted arithmetic mean of unsaturated hydraulic conductivity at the adjoining nodes
i.e.,
By substituting Eqs. (9) to (13) in Eq. (8)
The residual equation [
i.e., Eq. (14)] is represented as
, which was written for all the nodes in the flow domain that results in a set of highly non-linear simultaneous algebraic equations in the unknown
. These equations were solved using an iteration method, Newton-Raphson technique, which was written for all the nodes by considering upper and lower boundary conditions. These equations form a tridiagonal matrix, which was solved using successive over Relaxation (SOR) method. At each iteration, new hydraulic parameters are simulated and solved to obtain the soil suction head and moisture content at all the nodes in flow domain. The iterative process will continue until all nodes reach a reasonable convergence (0.001). This process is repeated for all the time steps upto simulation time. Flow chart of the numerical solution is presented in Fig 4.
Initial condition
Initial moisture content, θ
ini was specified at all the nodes and corresponding values of Ψ and K(Ψ) were computed using soil moisture characteristic equations.
θ = θ
ini for 0<z<zmax
Boundary conditions
No flow was considered in lateral directions. The upper BC and lower BC have been applied at the top face of the first block and last block’s lower face, respectively. Upper BC was considered as the constant flux during irrigation period and no flux during no irrigation period. Lower BC was considered as free drainage.
Upper BC during irrigation
Where
q = Constant flux (cm/hr).
Q = Emitter discharge (cm
3/hr).
A = Surface area (cm
2).
Upper BC during no irrigation
Lower BC
Statistical evaluation
Statistical indices - root mean square error (RMSE), coefficient of determination (R
2) and Percentage Bias (PBIAS) between the depths of wetting fronts of experiments and present model at different times have been computed. Better model should have lower RMSE values, R
2 values near to 1 and smaller PBIAS values (
Ranjan and Singh, 2022).
Where,
X = Observed value of the present model.
Y = Observed value from experimental models.
n = Total number of observations.
i = 1, 2, 3,…n.