Agricultural Science Digest

  • Chief EditorArvind kumar

  • Print ISSN 0253-150X

  • Online ISSN 0976-0547

  • NAAS Rating 5.52

  • SJR 0.156

Frequency :
Bi-monthly (February, April, June, August, October and December)
Indexing Services :
BIOSIS Preview, Biological Abstracts, Elsevier (Scopus and Embase), AGRICOLA, Google Scholar, CrossRef, CAB Abstracting Journals, Chemical Abstracts, Indian Science Abstracts, EBSCO Indexing Services, Index Copernicus
Agricultural Science Digest, volume 43 issue 3 (june 2023) : 311-318

Numerical Model to Simulate Soil Wetting Pattern under Drip Irrigation System

Kowkuntla Rama Krishna Reddy1,*, Vivekanand Singh1
1Department of Civil Engineering, National Institute of Technology Patna-800 005, Bihar, India.
Cite article:- Reddy Krishna Rama Kowkuntla, Singh Vivekanand (2023). Numerical Model to Simulate Soil Wetting Pattern under Drip Irrigation System . Agricultural Science Digest. 43(3): 311-318. doi: 10.18805/ag.D-5619.
Background: A properly designed drip irrigation system increases the efficiency of irrigation water and its design is based on the soil wetting pattern.

Methods: In this study, a numerical model has been developed using one-dimensional Richard’s Equation to simulate soil wetting pattern under drip irrigation. A fully implicit finite-difference scheme has been used to solve the governing equations with suitable initial and boundary conditions. An experiment was also conducted to analyse the wetting pattern in sandy loam soil. Wetting front computed by numerical model has been compared with the experimental wetting front at different time interval of drip irrigation system. The present model is also validated with the results available in the literature.

Result: Statistical indices - root mean square error, percentage bias and coefficient of determination have been used to evaluate the model’s performance. The computed results are closely matched with the and experimental results with RMSE is less than 3, PBIAS less than 6% and R2 greater than 0.96 and also with the earlier models with RMSE is less than 0.01, PBIAS less than 4% and R2 greater than 0.98, It shows that the present model can simulate soil wetting pattern in the drip irrigation system.
Increasing water scarcity affects the irrigation practices in arid and semi-arid regions. Drip irrigation (DI) system is best alternative to other irrigation systems (Kadasiddappa and Rao, 2018). In DI, water is applied adjacent to the plant root to wet only the active root zone to maintain sufficient moisture for the plant. Since water is conveyed through a pipeline and directly applied to the root, it reduces the water loss during irrigation (Kumavat et al., 2016). The wetting pattern of soil is the most important parameter to be considered in the design of spacing between laterals and emitters under DI.

Several empirical, analytical and numerical models are developed to determine the wetting pattern under DI (Subbiah, 2013). In Empirical models, wetted depth and radius were found using the dimensional analysis and regression techniques by using experimental and field data (Al-Ogaidi et al., 2015). Empirical models are unreliable due to high variation in soil hydraulic properties. The analytical models (Kilic, 2020) are quick and simple but it assumes homogenous soil and specific soil and water properties. Numerical models (Simunek et al., 2008) consider actual conditions of soil and flow properties, that makes model more realistic and reliable (Warrick, 2003).

Most of the numerical and analytical models are based on the solutions of Richards Equation (RE) (Richard, 1931). Due to limitations of h-based and θ-based equations, the mixed form of RE is more reliable (Celia et al., 1990). The mixed form of RE provides mass balance and has no limitations to solve field problems (Singh and Bhallamudi, 1998). Unsaturated soil characteristics are required to solve the mixed form of the RE. Haverkamp et al. (1977) validated the experimental data of Vanchaud and Thony, (1971) for soil moisture movement in unsaturated zone. Gottardi and Venutelli, (1992) has used RE for infiltration and movement of moisture in the soil. Though the numerical models are very beneficial but, they are not easy to use. Numerical models are sensitive to the initial and boundary conditions, have difficulties in convergence and require good knowledge of subject.

In this paper, a One-Dimensional (1-D) numerical model has been developed to simulate soil wetting pattern using mixed form of RE. This model is simple to program, easy to implement any boundary conditions and can be applied for homogenous and heterogeneous soil types. Further, this model can be used with simple knowledge of the subject. Thus, the present model can help in designing a DI system for different top boundary conditions. The performance of the model is evaluated by using statistical indices.
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.

Fig 1: Layout of experimental setup.

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.

Table 1: Physical characteristics and Van Genuchten hydraulic parameters of soil.

Fig 2: Soil water retention curve of the sandy loam soil.

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,

θ = Moisture content (cm3/cm3).
Ψ = 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:



Ks = Saturated hydraulic conductivity (cm/sec).
Ψ = Suction head (cm).
θ(Ψ) = Moisture content depends on the pressure head (cm3/cm3).
θs = Saturated moisture content (cm3/cm3).
θr = Residual moisture content (cm3/cm3).
A, B, a, b = Haverkamp parameters.
α, m, n, l = Van Genuchten parameters.
Se = 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.

Fig 3: Finite-difference grid for one-dimensional subsurface soil column.

The finite-difference form of Eq. (1) is
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:

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

KII= Unsaturated hydraulic conductivity at interblock face II.
KIV= 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
q = Constant flux (cm/hr).
Q = Emitter discharge (cm3/hr).
A = Surface area (cm2).
Upper BC during no irrigation


Lower BC
Statistical evaluation
Statistical indices - root mean square error (RMSE), coefficient of determination (R2) 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, R2 values near to 1 and smaller PBIAS values (Ranjan and Singh, 2022).
X = Observed value of the present model.
Y = Observed value from experimental models.
n = Total number of observations.
i = 1, 2, 3,…n.
Verification of the model

Present model is verified using experimental results by comparing the wetted depth of both the results. Further, developed model was tested with the experimental results of Haverkamp et al., (1977) and Gottardi and Venutelli, (1992).

Verification with the experimental data

In laboratory experiment, discharge rates of 0.5, 0.85 and 2 lph were applied at the top and resulting wetting depths were measured at different time and presented in Fig 5. During an irrigation period of 20, 12 and 5 hrs for discharge rates 0.5, 0.85 and 2.0 lph, wetted depth were observed at 24.6, 23.5 and 19.5 cms, respectively. The wetted depth observed at 48 hrs was 27.6, 32.5 and 25 cms for discharge rate 0.5, 0.85 and 2.0 lph respectively after the irrigation was stopped. 

A soil column of 60 cm depth was taken with a grid size of Δz = 2 cm. Uniform initial moisture content was considered throughout the domain. The wetting front at different time intervals were simulated for 48 hrs with three discharge rates 0.5, 0.85 and 2 lph for 20, 12, 5 hrs, respectively. Time step,  Δt = 10 sec, tolerance of convergence = 0.001 and weighting factor, w = 1 was considered for simulation.

The numerical and experimental results of wetting front depth were compared at different times as in Fig 5 and Fig 6. The observed and simulated results are complementing to each other with small deviations. Statistical evaluation has been done for all the 3 discharge rates and results were presented in Table 2. Comparison shows that numerical model of 0.5 lph is slightly underestimating during irrigation and overestimating after irrigation was stopped. For 0.85 lph discharge rate, simulated wetting depth matches well with experiment, whereas, results of 2 lph shows good agreement during irrigation period and slight deviation after irrigation. These slight deviations may be due to the flow in horizontal direction in the experiment.

Fig 5: Observed and simulated soil moisture profiles at different times and discharge rates(a) 0.5 lph, (b) 0.85 lph and (c) 2 lph.

Fig 6: comparison of wetting front depth between observed experimentally and simulated numerically for discharge rate (a) 0.5 lph (b) 0.85 lph and (c) 2 lph.

Table 2: Statistical parameters for comparison of wetting front depth of the present numerical model and experimental result.

Verification with the results of past study
Haverkamp et al., (1977) simulated numerically the experimental results conducted by Vachaud and Thony, (1971). Parameters used for soil characteristics are presented in Table 3. Soil column of 80 cm with uniform initial moisture content 0.1 cm3/ cm3 was considered. Constant flux of 13.69 cm/hr was applied as upper BC. Present model has been simulated with the same parameters of Haverkamp et al., (1977). Grid size, Δz = 4 cm,  time step, Δt = 10 sec were considered in the simulation. The results of present model are presented in Fig 7, which match well with Haverkamp et al., (1977). The wetting depth reaches 18, 50 and 70 cm depths at time 0.1, 0.5 and 0.8 hr, respectively. Statistical indices- RMSE, R2 and PBIAS, were computed and presented in Table 4. RMSE values indicates better performance of the model, R2 values were close to 1, indicates no degree of variance. The PBIAS was less than 5%, which shows small variation in the results.

Fig 7: Verifications of the model with Haverkamp et al., (1977) for wetting fronts at different times along depth.

Table 3: Soil hydraulic parameters for Haverkamp model (Haverkamp et al., 1977) and Van Genuchten model (Gottardi and Venutelli 1992).

Present model has also been validated using the results of Gottardi and Venutelli, (1992) to check the applicability and versatility with different types of soil. Soil characteristic parameters are presented in Table 3. A soil column of 100 cm with the initial moisture content of 0.242 cm3/ cm3 from top to 6 cm depth and 0.143 cm3/ cm3 from depth 6 cm onwards was considered. The constant flux rate of 10.7 cm/hr was applied and grid size, Δz = 4 cm, time step, Δt = 10 sec were considered. Fig 8 shows the wetting fronts computed using the present model and Gottardi and Venutelli, (1992) and observed that both reached 30, 56 and 90 cm depth at time 0.28, 0.67 and 1.33 hrs. RMSE, R2 and PBIAS were computed and presented in Table 4. Statistical evaluation shows satisfactory performance of the model.

Sensitivity analysis has also been carried out by varying time step (Δt) and nodal spacing (Δz). The results are satisfactory upto Δt = 60 sec and Dz = 5 cm.

Fig 8: Verifications of the model with Gottardi and Venutelli, (1992) for wetting fronts at different times along depth.

Table 4: Statistical parameters for performance evaluation of the present model with past literature.

The 1-D numerical model was developed to simulate the wetting front in DI system by considering the soil hydraulic properties, discharge rate and appropriate initial and BCs. The model is validated with the experiments conducted in. The model’s performance has also been evaluated using statistical indices-RMSE, R2 and PBAIS. RMSE value was less than 5 in all the tests, which indicated better model performance. The R2 was always near to 1, which shows a strong relationship between the results of the present numerical model with experimental results and earlier models. The PBIAS was less than 10%, which indicates small deviation. Thus, it can be concluded that the present model gives satisfactory results and can be used to analyse the distribution of soil wetting front in DI system, which reduces the cost and labor for the field experiments.

  1. Al-Ogaidi, A.A., Wayayok, A., Kamal, M.R., Abdullah, A.F. (2015). A modified empirical model for estimating the wetted zone dimensions under drip irrigation. Jurnal Teknologi. 76(15): 69-73. 

  2. Celia, M.A., Bououtas, E.T., Zarba, R.L. (1990). A general mass- conservative numerical solution for the unsaturated flow equation. Water Resources Research. 26(7): 1483-1496. 

  3. Gottardi, G. and Venutelli, M. (1992). Moving finite element model for one-dimensional infiltration in unsaturated soil. Water Resources Research. 28: 3259-3267.

  4. Haverkamp, R., Vauclin, M., Touma, J., Wierenga, P.J., Vachaud, G. (1977). A comparison of numerical simulation models for one-dimensional infiltration. Soil Science Society of American Journal. 41: 286-294.

  5. Kadasiddappa, M.M. and Rao, V.P. (2018). Irrigation scheduling through drip and surface methods- A critical review on growth, yield, nutrient uptake and water use studies of rabi maize. Agricultural Reviews. 39: 300-306.

  6. Kandelous, M.M. and Šimùnek, J. (2010). Comparison of numerical, analytical and empirical models to estimate wetting patterns for surface and subsurface drip irrigation. Irrigation Science. 28(5): 435-444.

  7. Kilic, M. (2020). A new analytical method for estimating the 3D volumetric wetting pattern under drip irrigation system. Agricultural Water Management. 228: 105-898. 

  8. Kumawat, P.D., Kacha, D.J., Dahima, N.U. (2016). Effect of crop geometry and drip irrigation levels on sugarcane in south Saurashtra region of India. Indian Journal of Agricultural Research. 50(4): 366-369.

  9. Prajapat, A.L., Saxena, R., Choudhary, R.R. and Kumhar, M. (2020). Consumptive use, water use efficiency, moisture depletion pattern of wheat under semi-arid eastern plain zone of Rajasthan. Agricultural Science Digest. 40(4): 387-391.

  10. Ranjan, S. and Singh, V. (2022). HEC-HMS based rainfall-runoff model for Punpun river basin. Water Practice and Technology.  17(5): 986-1001.

  11. Richards, L.A. (1931). Capillary conduction of liquids through porous mediums. Physics 1: 318-333.

  12. Simunek, J., Van Genuchten, M.T., Sejnac, M. (2008). Development and applications of the hydrus and stanmod software packages and related codes. Vadose Zone Journal. 7: 587-600.

  13. Singh, V. and Bhallamudi, S.M. (1998). Conjunctive surface-subsurface modelling of overland flow. Advances in Water Resources. 21: 567-579.

  14. Subbaiah, R. (2013). A review of models for predicting soil water dynamics during trickle irrigation. Irrigation Science. 31: 225-258. 

  15. Vachaud, G. and Thony, J. (1971). Hysteresis during infiltration and redistribution in a soil column at different initial water contents. Water Resources Research. 7(1): 111-127.

  16. Van Genuchten, M.T. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal. 44: 892-889.

  17. Warrick, A.W. (2003). Soil Water Dynamics. Oxford University Press, New York.

Editorial Board

View all (0)